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
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
|
/* mpfr.h -- Include file for mpfr.
Copyright (C) 1999 Free Software Foundation.
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 Library General Public License as published by
the Free Software Foundation; either version 2 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 Library General Public
License for more details.
You should have received a copy of the GNU Library 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. */
#include <stdio.h>
/* Cygnus does not know about *rand48 functions */
#ifdef __CYGWIN32__
#define mrand48 rand
#define drand48 rand
#define lrand48 rand
#define srand48 srand
#endif
/* Definition of rounding modes */
#define GMP_RNDN 0
#define GMP_RNDZ 1
#define GMP_RNDU 2
#define GMP_RNDD 3
/* Definition of constants */
#define LOG2 0.69314718055994528622 /* log(2) rounded to zero on 53 bits */
/* Definitions of types and their semantics */
typedef unsigned long int mp_prec_t; /* easy to change if necessary */
typedef int mp_rnd_t; /* preferred to char */
typedef struct {
mp_prec_t _mp_prec; /* WARNING : for the mpfr type, the precision */
/* should be understood as the number of BITS,*/
/* not the number of mp_limb_t's. This means */
/* that the corresponding number of allocated
limbs is >= ceil(_mp_prec/BITS_PER_MP_LIMB) */
mp_size_t _mp_size; /* MPFR_ABSSIZE(.) is the number of allocated
limbs the field _mp_d points to.
The sign is that of _mp_size.
The number 0 is such that _mp_d[k-1]=0
where k = ceil(_mp_prec/BITS_PER_MP_LIMB) */
mp_exp_t _mp_exp;
mp_limb_t *_mp_d;
}
__mpfr_struct;
/*
The number represented is
sign(_mp_size)*(_mp_d[k-1]/B+_mp_d[k-2]/B^2+...+_mp_d[0]/B^k)*2^_mp_exp
where k=ceil(_mp_prec/BITS_PER_MP_LIMB) and B=2^BITS_PER_MP_LIMB.
For the msb (most significant bit) normalized representation, we must have
_mp_d[k-1]>=B/2, unless the number is zero (in that case its sign is still
given by sign(_mp_size)).
We must also have the last k*BITS_PER_MP_LIMB-_mp_prec bits set to zero.
*/
typedef __mpfr_struct mpfr_t[1];
typedef __mpfr_struct *mpfr_ptr;
typedef __gmp_const __mpfr_struct *mpfr_srcptr;
/* 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)))
/* Prototypes */
#ifndef _PROTO
#if defined (__STDC__) || defined (__cplusplus)
#define _PROTO(x) x
#else
#define _PROTO(x) ()
#endif
#endif
/* bit 31 of _mp_size is used for sign,
bit 30 of _mp_size is used for Nan flag,
bit 29 of _mp_size is used for Inf flag,
remaining bits are used to store the number of allocated limbs */
#define MPFR_CLEAR_FLAGS(x) (((x) -> _mp_size &= ~(3 << 29)))
#define MPFR_IS_NAN(x) (((x)->_mp_size >> 30)&1)
#define MPFR_SET_NAN(x) ((x)->_mp_size |= (1<<30))
#define MPFR_IS_INF(x) (((x)->_mp_size >> 29)&1)
#define MPFR_SET_INF(x) ((x)->_mp_size |= (1<<29))
#define MPFR_RESET_INF(x) ((x)->_mp_size &= (~0 - (1<<29)))
#define MPFR_IS_FP(x) ((((x) -> _mp_size >> 29) & 3) == 0)
#define MPFR_SET_IS_FP(x) ((x) -> _mp_size &= 2684354559UL); /* 1001111...1 */
#define MPFR_ABSSIZE(x) ((x)->_mp_size & ((1<<29)-1))
#define MPFR_SIZE(x) ((x)->_mp_size)
#define MPFR_EXP(x) ((x)->_mp_exp)
#define MPFR_MANT(x) ((x)->_mp_d)
#define MPFR_SIGN(x) (((x)->_mp_size >> 31) ? -1 : 1)
#define MPFR_ISNONNEG(x) (MPFR_NOTZERO((x)) && MPFR_SIGN(x)>=0)
#define MPFR_ISNEG(x) (MPFR_NOTZERO((x)) && MPFR_SIGN(x)==-1)
#define MPFR_CHANGE_SIGN(x) (MPFR_SIZE(x) ^= (((mp_size_t)1)<<31))
#define MPFR_SET_SAME_SIGN(x, y) if (MPFR_SIGN((x)) != MPFR_SIGN((y))) { MPFR_CHANGE_SIGN((x)); }
#define MPFR_PREC(x) ((x)->_mp_prec)
#define MPFR_NOTZERO(x) (MPFR_MANT(x)[(MPFR_PREC(x)-1)/BITS_PER_MP_LIMB])
#define MPFR_IS_ZERO(x) ((MPFR_NOTZERO(x))==0)
#define MPFR_SET_ZERO(x) (MPFR_MANT(x)[(MPFR_PREC(x)-1)/BITS_PER_MP_LIMB] = 0)
#define mpfr_sgn(x) ((MPFR_NOTZERO(x)) ? MPFR_SIGN(x) : 0)
/* reallocates the mantissa of x to q bits and sets the precision to q */
#define _mpfr_realloc(x, q) { \
(x)->_mp_d = (mp_ptr) (*_mp_reallocate_func) \
((x)->_mp_d, (MPFR_PREC(x)-1)>>3, (q+7)>>3); \
MPFR_PREC(x) = q; }
#if defined (__cplusplus)
extern "C" {
#endif
void mpfr_init2 _PROTO ((mpfr_ptr, mp_prec_t));
void mpfr_init _PROTO ((mpfr_ptr));
int mpfr_round_raw _PROTO ((mp_limb_t *, mp_limb_t *, mp_prec_t, int,
mp_prec_t, mp_rnd_t));
int mpfr_round_raw2 _PROTO((mp_limb_t *, mp_prec_t, int, mp_rnd_t, mp_prec_t));
void mpfr_round _PROTO ((mpfr_ptr, mp_rnd_t, mp_prec_t));
int mpfr_can_round _PROTO ((mpfr_ptr, mp_prec_t, mp_rnd_t, mp_rnd_t,
mp_prec_t));
int mpfr_can_round_raw _PROTO ((mp_limb_t *, mp_prec_t, int, mp_prec_t,
mp_rnd_t, mp_rnd_t, mp_prec_t));
void mpfr_set_d _PROTO ((mpfr_ptr, double, mp_rnd_t));
int mpfr_set_z _PROTO ((mpfr_ptr, mpz_srcptr, mp_rnd_t));
mp_exp_t mpz_set_fr _PROTO ((mpz_ptr, mpfr_srcptr));
void mpfr_set_q _PROTO ((mpfr_ptr, mpq_srcptr, mp_rnd_t));
double mpfr_get_d _PROTO ((mpfr_srcptr));
double mpfr_get_d2 _PROTO ((mpfr_srcptr, long));
void mpfr_set_f _PROTO ((mpfr_ptr, mpf_srcptr, mp_rnd_t));
void mpfr_set_si _PROTO ((mpfr_ptr, long, mp_rnd_t));
void mpfr_set_ui _PROTO ((mpfr_ptr, unsigned long, mp_rnd_t));
void mpfr_print_raw _PROTO ((mpfr_srcptr));
void mpfr_random _PROTO ((mpfr_ptr));
void mpfr_random2 _PROTO ((mpfr_ptr, mp_size_t, mp_exp_t));
void mpfr_urandomb _PROTO ((mpfr_ptr, gmp_randstate_t));
void mpfr_clear _PROTO ((mpfr_ptr));
void mpfr_set_str_raw _PROTO ((mpfr_ptr, char *));
int mpfr_set_str _PROTO ((mpfr_ptr, char *, int, mp_rnd_t));
int mpfr_init_set_str _PROTO ((mpfr_ptr, char *, int, mp_rnd_t));
size_t mpfr_inp_str _PROTO ((mpfr_ptr, FILE *, int, mp_rnd_t));
void mpfr_get_str_raw _PROTO ((char *, mpfr_srcptr));
char* mpfr_get_str _PROTO ((char *, mp_exp_t *, int, size_t, mpfr_srcptr, mp_rnd_t));
size_t mpfr_out_str _PROTO ((FILE *, int, size_t, mpfr_srcptr, mp_rnd_t));
void mpfr_mul _PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mp_rnd_t));
int mpfr_pow_ui _PROTO ((mpfr_ptr, mpfr_srcptr, unsigned long int, mp_rnd_t));
int mpfr_ui_pow_ui _PROTO ((mpfr_ptr, unsigned long int, unsigned long int,
mp_rnd_t));
mp_limb_t mpn_divrem_n _PROTO ((mp_limb_t *, mp_limb_t *, mp_limb_t *, mp_size_t));
mp_size_t mpn_sqrtrem_new _PROTO ((mp_limb_t *, mp_limb_t *, mp_limb_t *, mp_size_t));
void mpfr_div _PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mp_rnd_t));
void mpfr_agm _PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mp_rnd_t));
int mpfr_sqrt _PROTO ((mpfr_ptr, mpfr_srcptr, mp_rnd_t));
int mpfr_sqrt_ui _PROTO ((mpfr_ptr, unsigned long, mp_rnd_t));
void mpfr_add _PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mp_rnd_t));
void mpfr_add_ui _PROTO ((mpfr_ptr, mpfr_srcptr, unsigned long, mp_rnd_t));
void mpfr_sub_ui _PROTO ((mpfr_ptr, mpfr_srcptr, unsigned long, mp_rnd_t));
void mpfr_add_one_ulp _PROTO ((mpfr_ptr));
void mpfr_sub _PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mp_rnd_t));
void mpfr_ui_sub _PROTO ((mpfr_ptr, unsigned long, mpfr_srcptr, mp_rnd_t));
void mpfr_reldiff _PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mp_rnd_t));
void mpfr_set4 _PROTO ((mpfr_ptr, mpfr_srcptr, mp_rnd_t, int));
void mpfr_const_pi _PROTO ((mpfr_ptr, mp_rnd_t));
void mpfr_const_log2 _PROTO ((mpfr_ptr, mp_rnd_t));
int mpfr_log _PROTO ((mpfr_ptr, mpfr_srcptr, mp_rnd_t));
int mpfr_exp _PROTO ((mpfr_ptr, mpfr_srcptr, mp_rnd_t));
int mpfr_zeta _PROTO ((mpfr_ptr, mpfr_srcptr, mp_rnd_t));
int mpfr_sin_cos _PROTO ((mpfr_ptr, mpfr_ptr, mpfr_srcptr, mp_rnd_t));
void mpfr_mul_ui _PROTO((mpfr_ptr, mpfr_srcptr, unsigned long int, mp_rnd_t));
void mpfr_set_machine_rnd_mode _PROTO ((mp_rnd_t));
int mpfr_cmp3 _PROTO ((mpfr_srcptr, mpfr_srcptr, long int));
int mpfr_cmp_ui_2exp _PROTO ((mpfr_srcptr, unsigned long int, int));
int mpfr_cmp_si_2exp _PROTO ((mpfr_srcptr, long int, int));
int mpfr_cmp2 _PROTO ((mpfr_srcptr, mpfr_srcptr));
void mpfr_mul_2exp _PROTO((mpfr_ptr, mpfr_srcptr, unsigned long int,mp_rnd_t));
void mpfr_div_2exp _PROTO((mpfr_ptr, mpfr_srcptr, unsigned long int,mp_rnd_t));
void mpfr_set_prec _PROTO((mpfr_ptr, mp_prec_t));
void mpfr_set_prec_raw _PROTO((mpfr_ptr, mp_prec_t));
void mpfr_set_default_prec _PROTO((mp_prec_t));
extern mp_prec_t __gmp_default_fp_bit_precision;
extern mp_rnd_t __gmp_default_rounding_mode;
char * mpfr_print_rnd_mode _PROTO((mp_rnd_t));
void mpfr_neg _PROTO((mpfr_ptr, mpfr_srcptr, mp_rnd_t));
void mpfr_sub_one_ulp _PROTO((mpfr_ptr));
int mpfr_div_ui _PROTO((mpfr_ptr, mpfr_srcptr, unsigned long int, mp_rnd_t));
void mpfr_ui_div _PROTO((mpfr_ptr, unsigned long int, mpfr_srcptr, mp_rnd_t));
mp_prec_t mpfr_get_prec _PROTO((mpfr_srcptr));
void mpfr_set_default_rounding_mode _PROTO((mp_rnd_t));
int mpfr_eq(mpfr_srcptr, mpfr_srcptr, unsigned long);
void mpfr_floor _PROTO((mpfr_ptr, mpfr_srcptr));
void mpfr_trunc _PROTO((mpfr_ptr, mpfr_srcptr));
void mpfr_ceil _PROTO((mpfr_ptr, mpfr_srcptr));
void mpfr_extract _PROTO((mpz_ptr, mpfr_srcptr, unsigned int));
#if defined (__cplusplus)
}
#endif
#define mpfr_cmp_ui(b,i) mpfr_cmp_ui_2exp(b,i,0)
#define mpfr_cmp_si(b,i) mpfr_cmp_si_2exp(b,i,0)
#define mpfr_set(a,b,r) mpfr_set4(a,b,r,MPFR_SIGN(b))
#define mpfr_abs(a,b,r) mpfr_set4(a,b,r,1)
#define mpfr_cmp(b,c) mpfr_cmp3(b,c,1)
#if (BITS_PER_MP_LIMB==32)
#define MPFR_LIMBS_PER_DOUBLE 2
#elif (BITS_PER_MP_LIMB==64)
#define MPFR_LIMBS_PER_DOUBLE 1
#endif
#define mpfr_init_set_si(x, i, rnd) \
do { mpfr_init(x); mpfr_set_si((x), (i), (rnd)); } while (0)
#define mpfr_init_set_ui(x, i, rnd) \
do { mpfr_init(x); mpfr_set_ui((x), (i), (rnd)); } while (0)
#define mpfr_init_set_d(x, d, rnd) \
do { mpfr_init(x); mpfr_set_d((x), (d), (rnd)); } while (0)
#define mpfr_init_set(x, y, rnd) \
do { mpfr_init(x); mpfr_set((x), (y), (rnd)); } while (0)
#define mpfr_init_set_f(x, y, rnd) \
do { mpfr_init(x); mpfr_set_f((x), (y), (rnd)); } while (0)
|
