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
|
/* Utilities for MPFR developers, not exported.
Copyright (C) 1999-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. */
/* 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) */
#define SAFE_ABS(type,x) ((x) >= 0 ? (unsigned type)(x) : -(unsigned 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)))
/* 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
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_d2 _PROTO ((mpfr_srcptr, mp_exp_t));
mp_size_t mpn_sqrtrem_new _PROTO ((mp_limb_t *, mp_limb_t *, mp_limb_t *, mp_size_t));
int mpfr_cmp_abs _PROTO ((mpfr_srcptr, mpfr_srcptr));
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));
#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
|
