summaryrefslogtreecommitdiff
path: root/mpfr-impl.h
blob: b8d35e38e972a270bff7f9146a669f5fa9f27529 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
/* Utilities for MPFR developers, not exported.

Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc.

This file is part of the MPFR Library.

The MPFR Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or (at your
option) any later version.

The MPFR Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
License for more details.

You should have received a copy of the GNU Lesser General Public License
along with the MPFR Library; see the file COPYING.LIB.  If not, write to
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */

/* Definition of MPFR_LIMB_HIGHBIT */

#ifdef GMP_LIMB_HIGHBIT
#define MPFR_LIMB_HIGHBIT GMP_LIMB_HIGHBIT
#else
#define MPFR_LIMB_HIGHBIT MP_LIMB_T_HIGHBIT
#endif

#if GMP_NAIL_BITS != 0
#error "MPFR doesn't support nonzero values of GMP_NAIL_BITS"
#endif

/* Test if X (positive) is a power of 2 */

#define IS_POW2(X) (((X) & ((X) - 1)) == 0)
#define NOT_POW2(X) (((X) & ((X) - 1)) != 0)

/* This unsigned type must correspond to the signed one defined in gmp.h */
#if defined (_CRAY) && ! defined (_CRAYMPP)
typedef unsigned int            mp_exp_unsigned_t;
typedef unsigned int            mp_size_unsigned_t;
#else
typedef unsigned long int       mp_exp_unsigned_t;
typedef unsigned long int       mp_size_unsigned_t;
#endif

#define MP_EXP_T_MAX ((mp_exp_t) ((~ (mp_exp_unsigned_t) 0) >> 1))
#define MP_EXP_T_MIN (-MP_EXP_T_MAX-1)

#define MP_LIMB_T_ONE ((mp_limb_t) 1)

#if (BITS_PER_MP_LIMB & (BITS_PER_MP_LIMB - 1))
#error "BITS_PER_MP_LIMB must be a power of 2"
#endif

#define MPFR_INTPREC_MAX (ULONG_MAX & ~(unsigned long) (BITS_PER_MP_LIMB - 1))

/* Assertions */

/* Compile with -DWANT_ASSERT to check all assert statements */

/* MPFR_ASSERTN(expr): assertions that should always be checked */
/* #define MPFR_ASSERTN(expr) ASSERT_ALWAYS(expr) */
#define MPFR_ASSERTN(expr)  ((expr) ? (void) 0 : (void) ASSERT_FAIL (expr))

/* MPFR_ASSERTD(expr): assertions that should be checked when testing */
/* #define MPFR_ASSERTD(expr) ASSERT(expr) */
#if WANT_ASSERT
#define MPFR_ASSERTD(expr)  ASSERT_ALWAYS (expr)
#else
#define MPFR_ASSERTD(expr)  ((void) 0)
#endif

/* Definition of constants */

#define LOG2 0.69314718055994528622 /* log(2) rounded to zero on 53 bits */
#define ALPHA 4.3191365662914471407 /* a+2 = a*log(a), rounded to +infinity */

/* Safe absolute value (to avoid possible integer overflow) */
/* type is the target (unsigned) type */

#define SAFE_ABS(type,x) ((x) >= 0 ? (type)(x) : -(type)(x))

/* macros for doubles, based on gmp union ieee_double_extract */

typedef union ieee_double_extract Ieee_double_extract;

#define DOUBLE_ISNANorINF(x) (((Ieee_double_extract *)&(x))->s.exp == 0x7ff)
#define DOUBLE_ISINF(x) (DOUBLE_ISNANorINF(x) && \
			 (((Ieee_double_extract *)&(x))->s.manl == 0) && \
                         (((Ieee_double_extract *)&(x))->s.manh == 0))
#define DOUBLE_ISNAN(x) (DOUBLE_ISNANorINF(x) && \
			 ((((Ieee_double_extract *)&(x))->s.manl != 0) || \
                         (((Ieee_double_extract *)&(x))->s.manh != 0)))

#define DBL_POS_INF (1.0/0.0)
#define DBL_NEG_INF (-1.0/0.0)
#define DBL_NAN (0.0/0.0)

/* bit 31 of _mpfr_size is used for sign,
   bit 30 of _mpfr_size is used for Nan flag,
   bit 29 of _mpfr_size is used for Inf flag,
   remaining bits are used to store the number of allocated limbs */
#define MPFR_CLEAR_FLAGS(x) \
  (((x) -> _mpfr_size &= ~((mp_size_unsigned_t) 3 << 29)))
#define MPFR_IS_NAN(x) (((x)->_mpfr_size >> 30) & 1)
#define MPFR_SET_NAN(x) ((x)->_mpfr_size |= ((mp_size_unsigned_t) 1 << 30))
#define MPFR_CLEAR_NAN(x) \
  (((x) -> _mpfr_size &= ~((mp_size_unsigned_t) 1 << 30)))
#define MPFR_IS_INF(x) (((x)->_mpfr_size >> 29) & 1)
#define MPFR_SET_INF(x) ((x)->_mpfr_size |= ((mp_size_unsigned_t) 1 << 29))
#define MPFR_CLEAR_INF(x) ((x)->_mpfr_size &= ~((mp_size_unsigned_t) 1 << 29))
#define MPFR_IS_FP(x) ((((x) -> _mpfr_size >> 29) & 3) == 0)
#define MPFR_ABSSIZE(x) \
  ((x)->_mpfr_size & (((mp_size_unsigned_t) 1 << 29) - 1))
#define MPFR_SET_ABSSIZE(x, n) \
  ((x)->_mpfr_size = ((x)->_mpfr_size & ((mp_size_unsigned_t) 7 << 29)) \
                     | (mp_size_unsigned_t) (n))
#define MPFR_SIZE(x) ((x)->_mpfr_size)
#define MPFR_EXP(x) ((x)->_mpfr_exp)
#define MPFR_MANT(x) ((x)->_mpfr_d)
#define MPFR_ISNONNEG(x) (MPFR_NOTZERO((x)) && MPFR_SIGN(x) >= 0)
#define MPFR_ISNEG(x) (MPFR_NOTZERO((x)) && MPFR_SIGN(x) == -1)
#define MPFR_SET_POS(x) (MPFR_SIZE(x) &= ~(((mp_size_unsigned_t) 1) << 31))
#define MPFR_SET_NEG(x) (MPFR_SIZE(x) |= (((mp_size_unsigned_t) 1) << 31))
#define MPFR_CHANGE_SIGN(x) (MPFR_SIZE(x) ^= (((mp_size_unsigned_t) 1) << 31))
#define MPFR_SET_SAME_SIGN(x, y) \
  (MPFR_SIGN((x)) != MPFR_SIGN((y)) && (MPFR_CHANGE_SIGN((x)), 0))
#define MPFR_PREC(x) ((x)->_mpfr_prec)
#define MPFR_NOTZERO(x) \
  (MPFR_MANT(x)[(MPFR_PREC(x)-1)/BITS_PER_MP_LIMB] != (mp_limb_t) 0)
#define MPFR_IS_ZERO(x) \
  (MPFR_MANT(x)[(MPFR_PREC(x)-1)/BITS_PER_MP_LIMB] == (mp_limb_t) 0)
#define MPFR_SET_ZERO(x) \
  (MPFR_MANT(x)[(MPFR_PREC(x)-1)/BITS_PER_MP_LIMB] = (mp_limb_t) 0)
#define MPFR_ESIZE(x) \
  ((MPFR_PREC((x)) - 1) / BITS_PER_MP_LIMB + 1); 
#define MPFR_EVEN_INEX 2

/* When returning the ternary inexact value, ALWAYS use one of the
   following two macros, unless the flag comes from another function
   returning the ternary inexact value */
#define MPFR_RET(I) return \
  (I) ? ((__mpfr_flags |= MPFR_FLAGS_INEXACT), (I)) : 0
#define MPFR_RET_NAN return (__mpfr_flags |= MPFR_FLAGS_NAN), 0

/* The following macro restores the exponent range and the flags,
   checks that the result is in the exponent range and returns the
   ternary inexact value. */
#define MPFR_RESTORE_RET(inex, x, rnd_mode) \
  do \
    { \
      int inex_cr; \
      mpfr_restore_emin_emax(); \
      inex_cr = mpfr_check_range(x, rnd_mode); \
      if (inex_cr) \
        return inex_cr; /* underflow or overflow */ \
      MPFR_RET(inex); \
    } \
  while(0)

/* Memory gestion */

/* temporary allocate s limbs at xp, and initialize mpfr variable x */
#define MPFR_INIT(xp, x, p, s) \
  (xp = (mp_ptr) TMP_ALLOC((size_t) (s) * BYTES_PER_MP_LIMB), \
   MPFR_PREC(x) = (p), \
   MPFR_MANT(x) = (xp), \
   MPFR_SIZE(x) = (s), \
   MPFR_EXP(x) = 0)
/* same when xp is already allocated */
#define MPFR_INIT1(xp, x, p, s) \
  (MPFR_PREC(x) = (p), MPFR_MANT(x) = (xp), MPFR_SIZE(x) = (s))

#ifndef _PROTO
#if defined (__STDC__) || defined (__cplusplus)
#define _PROTO(x) x
#else
#define _PROTO(x) ()
#endif
#endif

#if defined (__cplusplus)
extern "C" {
#endif  

extern mpfr_t __mpfr_const_log2;
extern mp_prec_t __mpfr_const_log2_prec;

extern mpfr_t __mpfr_const_pi;
extern mp_prec_t __mpfr_const_pi_prec;

int mpfr_set_underflow _PROTO ((mpfr_ptr, mp_rnd_t, int));
int mpfr_set_overflow _PROTO ((mpfr_ptr, mp_rnd_t, int));
void mpfr_save_emin_emax _PROTO ((void));
void mpfr_restore_emin_emax _PROTO ((void));
int mpfr_add1 _PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr,
                       mp_rnd_t, mp_exp_unsigned_t));
int mpfr_sub1 _PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr,
                       mp_rnd_t, int));
int mpfr_round_raw_generic _PROTO ((mp_limb_t *, mp_limb_t *, mp_prec_t, int,
				    mp_prec_t, mp_rnd_t, int *, int));
int mpfr_can_round_raw _PROTO ((mp_limb_t *, mp_size_t, int, mp_exp_t,
				mp_rnd_t, mp_rnd_t, mp_prec_t));
double mpfr_get_d3 _PROTO ((mpfr_srcptr, mp_exp_t, mp_rnd_t));
int mpfr_cmp2 _PROTO ((mpfr_srcptr, mpfr_srcptr, mp_prec_t *));
long _mpfr_ceil_log2 _PROTO ((double));
long _mpfr_floor_log2 _PROTO ((double));
double _mpfr_ceil_exp2 _PROTO ((double));
unsigned long _mpfr_isqrt _PROTO ((unsigned long));
unsigned long _mpfr_cuberoot _PROTO ((unsigned long));
int mpfr_exp_2 _PROTO ((mpfr_ptr, mpfr_srcptr, mp_rnd_t));
int mpfr_exp3 _PROTO ((mpfr_ptr, mpfr_srcptr, mp_rnd_t));
int mpfr_powerof2_raw _PROTO ((mpfr_srcptr));

#define mpfr_round_raw(yp, xp, xprec, neg, yprec, r, inexp) \
  mpfr_round_raw_generic((yp), (xp), (xprec), (neg), (yprec), (r), (inexp), 0)

#define mpfr_round_raw2(xp, xn, neg, r, prec) \
  mpfr_round_raw_generic(0, (xp), (xn) * BITS_PER_MP_LIMB, (neg), \
			 (prec), (r), 0, 1); 

#if defined (__cplusplus)
}
#endif