diff options
| -rw-r--r-- | add.c | 29 | ||||
| -rw-r--r-- | agm.c | 17 | ||||
| -rw-r--r-- | div_ui.c | 14 | ||||
| -rw-r--r-- | exp.c | 18 | ||||
| -rw-r--r-- | get_str.c | 53 | ||||
| -rw-r--r-- | log.c | 13 | ||||
| -rw-r--r-- | mul.c | 14 | ||||
| -rw-r--r-- | mul_ui.c | 18 | ||||
| -rw-r--r-- | neg.c | 10 |
9 files changed, 109 insertions, 77 deletions
@@ -25,7 +25,7 @@ MA 02111-1307, USA. */ #include "mpfr.h" extern void mpfr_sub1 _PROTO((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, - unsigned char, int)); + mp_rnd_t, int)); #define ONE ((mp_limb_t) 1) @@ -36,13 +36,13 @@ extern void mpfr_sub1 _PROTO((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, void #if __STDC__ mpfr_add1(mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, - unsigned char rnd_mode, int diff_exp) + mp_rnd_t rnd_mode, int diff_exp) #else mpfr_add1(a, b, c, rnd_mode, diff_exp) - mpfr_ptr a; - mpfr_srcptr b; - mpfr_srcptr c; - unsigned char rnd_mode; + mpfr_ptr a; + mpfr_srcptr b; + mpfr_srcptr c; + mp_rnd_t rnd_mode; int diff_exp; #endif { @@ -70,7 +70,7 @@ mpfr_add1(a, b, c, rnd_mode, diff_exp) bn = (PREC(b)-1)/mp_bits_per_limb+1; /* number of significant limbs of b */ EXP(a) = EXP(b); - if (SIGN(a)!=SIGN(b)) CHANGE_SIGN(a); + if (MPFR_SIGN(a) * MPFR_SIGN(b) < 0) CHANGE_SIGN(a); /* case 1: diff_exp>=prec(a), i.e. c only affects the last bit through rounding */ @@ -364,14 +364,13 @@ mpfr_add1(a, b, c, rnd_mode, diff_exp) void #if __STDC__ -mpfr_add(mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, - unsigned char rnd_mode) +mpfr_add(mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) #else mpfr_add(a, b, c, rnd_mode) - mpfr_ptr a; - mpfr_srcptr b; - mpfr_srcptr c; - unsigned char rnd_mode; + mpfr_ptr a; + mpfr_srcptr b; + mpfr_srcptr c; + mp_rnd_t rnd_mode; #endif { int diff_exp; @@ -384,7 +383,7 @@ mpfr_add(a, b, c, rnd_mode) if (!NOTZERO(c)) { mpfr_set(a, b, rnd_mode); return; } diff_exp = EXP(b)-EXP(c); - if (SIGN(b) != SIGN(c)) { /* signs differ, it's a subtraction */ + if (MPFR_SIGN(b) * MPFR_SIGN(c) < 0) { /* signs differ, it's a subtraction */ if (diff_exp<0) { mpfr_sub1(a, c, b, rnd_mode, -diff_exp); } @@ -393,7 +392,7 @@ mpfr_add(a, b, c, rnd_mode) diff_exp = mpfr_cmp3(b,c,-1); /* if b>0 and diff_exp>0 or b<0 and diff_exp<0: abs(b) > abs(c) */ if (diff_exp==0) SET_ZERO(a); - else if (diff_exp*SIGN(b)>0) mpfr_sub1(a, b, c, rnd_mode, 0); + else if (diff_exp * MPFR_SIGN(b)>0) mpfr_sub1(a, b, c, rnd_mode, 0); else mpfr_sub1(a, c, b, rnd_mode, 0); } } @@ -34,16 +34,17 @@ MA 02111-1307, USA. */ void #ifdef __STDC__ -mpfr_agm(mpfr_ptr r, mpfr_srcptr op2, mpfr_srcptr op1, unsigned char rnd_mode) +mpfr_agm(mpfr_ptr r, mpfr_srcptr op2, mpfr_srcptr op1, mp_rnd_t rnd_mode) #else mpfr_agm(r, a, b, rnd_mode) - mpfr_ptr r; - mpfr_srcptr a; - mpfr_srcptr b; - unsigned char rnd_mode; + mpfr_ptr r; + mpfr_srcptr a; + mpfr_srcptr b; + mp_rnd_t rnd_mode; #endif { - int s, p, q, go_on; + int s, go_on; + mp_prec_t p, q; double uo, vo; mp_limb_t *up, *vp, *tmpp, *tmpup, *tmpvp, *ap, *bp; mpfr_t u, v, tmp, tmpu, tmpv, a, b; @@ -57,13 +58,13 @@ mpfr_agm(r, a, b, rnd_mode) /* If a or b is negative, the result is NaN */ - if ((SIGN(op1)<0)||(SIGN(op2)<0)) + if ((MPFR_SIGN(op1) < 0) || (MPFR_SIGN(op2) < 0)) { SET_NAN(r); return; } /* If a or b is 0, the result is 0 */ - if ((SIGN(op1)==0)||(SIGN(op2)==0)) + if ((NOTZERO(op1) && NOTZERO(op2)) == 0) { SET_ZERO(r); return; } @@ -32,13 +32,13 @@ MA 02111-1307, USA. */ /* returns 0 if result exact, non-zero otherwise */ int #ifdef __STDC__ -mpfr_div_ui(mpfr_ptr y, mpfr_srcptr x, unsigned long u, unsigned char rnd_mode) +mpfr_div_ui(mpfr_ptr y, mpfr_srcptr x, unsigned long int u, mp_rnd_t rnd_mode) #else mpfr_div_ui(y, x, u, rnd_mode) - mpfr_ptr y; + mpfr_ptr y; mpfr_srcptr x; - unsigned long u; - unsigned char rnd_mode; + unsigned long int u; + mp_rnd_t rnd_mode; #endif { int xn, yn, dif, sh, i; mp_limb_t *xp, *yp, *tmp, c, d; @@ -54,7 +54,7 @@ mpfr_div_ui(y, x, u, rnd_mode) xp = MANT(x); yp = MANT(y); EXP(y) = EXP(x); - if (SIGN(x)!=SIGN(y)) CHANGE_SIGN(y); + if (MPFR_SIGN(x) * MPFR_SIGN(y) < 0) CHANGE_SIGN(y); dif = yn+1-xn; #ifdef DEBUG @@ -108,10 +108,10 @@ printf("y="); mpfr_print_raw(y); putchar('\n'); case GMP_RNDZ: return 1; /* result is inexact */ case GMP_RNDU: - if (SIGN(y)>0) mpfr_add_one_ulp(y); + if (MPFR_SIGN(y)>0) mpfr_add_one_ulp(y); return 1; /* result is inexact */ case GMP_RNDD: - if (SIGN(y)<0) mpfr_add_one_ulp(y); + if (MPFR_SIGN(y)<0) mpfr_add_one_ulp(y); return 1; /* result is inexact */ case GMP_RNDN: if (d < ((mp_limb_t)1 << (sh-1))) return 1; @@ -35,12 +35,12 @@ MA 02111-1307, USA. */ */ int #if __STDC__ -mpfr_exp(mpfr_ptr y, mpfr_srcptr x, unsigned char rnd_mode) +mpfr_exp(mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) #else mpfr_exp(y, x, rnd_mode) - mpfr_ptr y; - mpfr_srcptr x; - unsigned char rnd_mode; + mpfr_ptr y; + mpfr_srcptr x; + mp_rnd_t rnd_mode; #endif { int n, expx, K, precy, q, k, l, expr, err; @@ -57,12 +57,12 @@ mpfr_exp(y, x, rnd_mode) /* if x > (2^31-1)*ln(2), then exp(x) > 2^(2^31-1) i.e. gives +infinity */ if (expx > 30) { - if (SIGN(x)>0) { printf("+infinity"); return 1; } + if (MPFR_SIGN(x)>0) { printf("+infinity"); return 1; } else { SET_ZERO(y); return 1; } } /* if x < 2^(-precy), then exp(x) i.e. gives 1 +/- 1 ulp(1) */ - if (expx < -precy) { int signx = SIGN(x); + if (expx < -precy) { int signx = MPFR_SIGN(x); mpfr_set_ui(y, 1, rnd_mode); if (signx>0 && rnd_mode==GMP_RNDU) mpfr_add_one_ulp(y); else if (signx<0 && (rnd_mode==GMP_RNDD || rnd_mode==GMP_RNDZ)) @@ -87,7 +87,7 @@ mpfr_exp(y, x, rnd_mode) /* if n<0, we have to get an upper bound of log(2) in order to get an upper bound of r = x-n*log(2) */ - mpfr_log2(s, (n>=0) ? GMP_RNDZ : GMP_RNDU); + mpfr_const_log2(s, (n>=0) ? GMP_RNDZ : GMP_RNDU); #ifdef DEBUG printf("n=%d log(2)=",n); mpfr_print_raw(s); putchar('\n'); #endif @@ -102,7 +102,7 @@ mpfr_exp(y, x, rnd_mode) printf(" ="); mpfr_print_raw(r); putchar('\n'); #endif mpfr_sub(r, x, r, GMP_RNDU); - if (SIGN(r)<0) { /* initial approximation n was too large */ + if (MPFR_SIGN(r)<0) { /* initial approximation n was too large */ n--; mpfr_mul_ui(r, s, (n<0) ? -n : n, GMP_RNDZ); if (n<0) mpfr_neg(r, r, GMP_RNDD); @@ -111,7 +111,7 @@ mpfr_exp(y, x, rnd_mode) #ifdef DEBUG printf("x-r=%1.20e\n",mpfr_get_d(r)); printf(" ="); mpfr_print_raw(r); putchar('\n'); - if (SIGN(r)<0) { fprintf(stderr,"Error in mpfr_exp: r<0\n"); exit(1); } + if (MPFR_SIGN(r)<0) { fprintf(stderr,"Error in mpfr_exp: r<0\n"); exit(1); } #endif mpfr_div_2exp(r, r, K, GMP_RNDU); /* r = (x-n*log(2))/2^K */ mpfr_set_ui(s, 1, GMP_RNDU); @@ -37,19 +37,19 @@ MA 02111-1307, USA. */ */ #if __STDC__ char *mpfr_get_str(char *str, mp_exp_t *expptr, int base, size_t n, - mpfr_srcptr op, unsigned char rnd_mode) + mpfr_srcptr op, mp_rnd_t rnd_mode) #else char *mpfr_get_str(str, expptr, base, n, op, rnd_mode) - char *str; + char *str; mp_exp_t *expptr; int base; size_t n; - mpfr_srcptr op; - unsigned char rnd_mode; + mpfr_srcptr op; + mp_rnd_t rnd_mode; #endif { double d; long e, q, div, p, err, prec, sh; mpfr_t a, b; mpz_t bz; - char *str0; unsigned char rnd1; int f, pow2, ok=0, neg; + char *str0; mp_rnd_t rnd1; int f, pow2, ok=0, neg, n0=n; if (base<2 || 36<base) { fprintf(stderr, "Error: too small or too large base in mpfr_get_str: %d\n", @@ -57,11 +57,11 @@ char *mpfr_get_str(str, expptr, base, n, op, rnd_mode) exit(1); } - neg = (SIGN(op)<0) ? 1 : 0; + neg = (MPFR_SIGN(op)<0) ? 1 : 0; if (!NOTZERO(op)) { if (str==NULL) str0=str=(*_mp_allocate_func)(neg + n + 2); - if (SIGN(op)<0) *str++ = '-'; + if (MPFR_SIGN(op)<0) *str++ = '-'; for (f=0;f<n;f++) *str++ = '0'; *expptr = 1; return str0; @@ -86,11 +86,17 @@ char *mpfr_get_str(str, expptr, base, n, op, rnd_mode) y*base^(f-n) <= x*2^(e-p) < (x+1)*2^(e-p) <= (y+1)*base^(f-n) which implies 2^(EXP(op)-PREC(op)) <= base^(f-n) */ - n = f + (int) ceil(((double)PREC(op)-e)*log(2.0)/log((double)base)); + n = f + (int) floor(((double)PREC(op)-e)*log(2.0)/log((double)base)); } +#ifdef DEBUG + printf("f=%d n=%d EXP(op)=%d PREC(op)=%d\n", f, n, e, PREC(op)); +#endif /* now the first n digits of the mantissa are obtained from rnd(op*base^(n-f)) */ prec = (long) ceil((double)n*log((double)base)/log(2.0)); +#ifdef DEBUG + printf("prec=%d\n", prec); +#endif err = 5; q = prec+err; /* one has to use at least q bits */ @@ -126,16 +132,26 @@ char *mpfr_get_str(str, expptr, base, n, op, rnd_mode) mpfr_div(a, b, a, rnd_mode); } /* now a is an approximation by default of 1/base^(f-n) */ +#ifdef DEBUG + printf("a=%1.20e\n", mpfr_get_d(a)); +#endif mpfr_mul(b, op, a, rnd_mode); } } if (neg) CHANGE_SIGN(b); /* put b positive */ - +#ifdef DEBUG + printf("p=%d b=%1.20e\n", p, mpfr_get_d(b)); + printf("q=%d 2*prec+BITS_PER_MP_LIMB=%d\n", q, 2*prec+BITS_PER_MP_LIMB); +#endif if (q>2*prec+BITS_PER_MP_LIMB) { /* happens when just in the middle between two digits */ n--; q-=BITS_PER_MP_LIMB; if (n==0) { - fprintf(stderr, "cannot determine leading digit\n"); exit(1); + fprintf(stderr, "Error in mpfr_get_str: cannot determine leading digit\n"); + printf("base=%d, digits=%u, mode=%s, op=", base, n0, + mpfr_print_rnd_mode(rnd_mode)); + mpfr_print_raw(op); putchar('\n'); + exit(1); } } ok = pow2 || mpfr_can_round(b, q-err, rnd_mode, rnd_mode, prec); @@ -159,19 +175,34 @@ char *mpfr_get_str(str, expptr, base, n, op, rnd_mode) else if ((rnd_mode==GMP_RNDU && neg==0) || (rnd_mode==GMP_RNDD && neg)) mpfr_sub(a, b, a, rnd_mode); else mpfr_add(a, b, a, rnd_mode); +#ifdef DEBUG + printf("a=%1.20e = ", mpfr_get_d(a)); mpfr_print_raw(a); putchar('\n'); + printf("b=%1.20e = ", mpfr_get_d(b)); mpfr_print_raw(b); putchar('\n'); +#endif /* check that a and b are rounded similarly */ prec=EXP(b); +#ifdef DEBUG + printf("%d %d\n", EXP(a), prec); +#endif if (EXP(a) != prec) ok=0; else { mpfr_round(b, rnd_mode, prec); mpfr_round(a, rnd_mode, prec); +#ifdef DEBUG + printf("a=%1.20e = ", mpfr_get_d(a)); mpfr_print_raw(a); putchar('\n'); + printf("b=%1.20e = ", mpfr_get_d(b)); mpfr_print_raw(b); putchar('\n'); +#endif if (mpfr_cmp(a, b)) ok=0; } } if (ok==0) { /* n is too large */ n--; if (n==0) { - fprintf(stderr, "cannot determine leading digit\n"); exit(1); + fprintf(stderr, "Error in mpfr_get_str: cannot determine leading digit\n"); + printf("base=%d, digits=%u, mode=%s, op=", base, n0, + mpfr_print_rnd_mode(rnd_mode)); + mpfr_print_raw(op); putchar('\n'); + exit(1); } q -= BITS_PER_MP_LIMB; } @@ -47,22 +47,23 @@ MA 02111-1307, USA. */ int #if __STDC__ -mpfr_log(mpfr_ptr r, mpfr_srcptr a, unsigned char rnd_mode) +mpfr_log(mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode) #else mpfr_log() mpfr_ptr r; mpfr_srcptr a; - unsigned char rnd_mode; + mp_rnd_t rnd_mode; #endif { - int p, m, q, bool, size, cancel; + int m, bool, size, cancel; + mp_prec_t p, q; mpfr_t cst, rapport, agm, tmp1, tmp2, s, mm; mp_limb_t *cstp, *rapportp, *agmp, *tmp1p, *tmp2p, *sp, *mmp; double ref; TMP_DECL(marker); /* If a is NaN or a is negative or null, the result is NaN */ - if (FLAG_NAN(a) || (SIGN(a)<=0)) + if (FLAG_NAN(a) || (NOTZERO(a)==0) || (MPFR_SIGN(a)<0)) { SET_NAN(r); return 1; } /* If a is 1, the result is 0 */ @@ -109,9 +110,9 @@ mpfr_log() mpfr_div(rapport,tmp2,s,GMP_RNDN); /* I compute 4/s, err <= 2 ulps */ mpfr_agm(agm,tmp1,rapport,GMP_RNDN); /* AG(1,4/s), err<=3 ulps */ mpfr_mul_2exp(tmp1,agm,1,GMP_RNDN); /* 2*AG(1,4/s), still err<=3 ulps */ - mpfr_pi(cst, GMP_RNDN); /* I compute pi, err<=1ulp */ + mpfr_const_pi(cst, GMP_RNDN); /* compute pi, err<=1ulp */ mpfr_div(tmp2,cst,tmp1,GMP_RNDN); /* pi/2*AG(1,4/s), err<=5ulps */ - mpfr_log2(cst,GMP_RNDN); /* I compute log(2), err<=1ulp */ + mpfr_const_log2(cst,GMP_RNDN); /* compute log(2), err<=1ulp */ mpfr_mul(tmp1,cst,mm,GMP_RNDN); /* I compute m*log(2), err<=2ulps */ cancel = EXP(tmp2); mpfr_sub(cst,tmp2,tmp1,GMP_RNDN); /* log(a), err<=7ulps+cancel */ @@ -31,13 +31,13 @@ MA 02111-1307, USA. */ void #if __STDC__ -mpfr_mul(mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, unsigned char rnd_mode) +mpfr_mul(mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) #else mpfr_mul(a, b, c, rnd_mode) - mpfr_ptr a; - mpfr_srcptr b; - mpfr_srcptr c; - unsigned char rnd_mode; + mpfr_ptr a; + mpfr_srcptr b; + mpfr_srcptr c; + mp_rnd_t rnd_mode; #endif { unsigned int bn, cn, an, tn, k; int cc; @@ -49,7 +49,7 @@ mpfr_mul(a, b, c, rnd_mode) if (FLAG_NAN(b) || FLAG_NAN(c)) { SET_NAN(a); return; } if (!NOTZERO(b) || !NOTZERO(c)) { SET_ZERO(a); return; } - sign_product = SIGN(b) * SIGN(c); + sign_product = MPFR_SIGN(b) * MPFR_SIGN(c); bn = (PREC(b)-1)/mp_bits_per_limb+1; /* number of significant limbs of b */ cn = (PREC(c)-1)/mp_bits_per_limb+1; /* number of significant limbs of c */ tn = (PREC(c)+PREC(b)-1)/mp_bits_per_limb+1; @@ -72,7 +72,7 @@ mpfr_mul(a, b, c, rnd_mode) ap[an-1] = (mp_limb_t) 1 << (BITS_PER_MP_LIMB-1); } EXP(a) = EXP(b) + EXP(c) + b1 - 1 + cc; - if (sign_product * SIGN(a)<0) CHANGE_SIGN(a); + if (sign_product * MPFR_SIGN(a)<0) CHANGE_SIGN(a); TMP_FREE(marker); return; } @@ -27,13 +27,13 @@ MA 02111-1307, USA. */ void #if __STDC__ -mpfr_mul_ui(mpfr_ptr y, mpfr_srcptr x, unsigned long u, unsigned char RND_MODE) +mpfr_mul_ui(mpfr_ptr y, mpfr_srcptr x, unsigned long int u, mp_rnd_t rnd_mode) #else -mpfr_mul_ui() - mpfr_ptr y; - mpfr_srcptr x; - unsigned long u; - unsigned char RND_MODE; +mpfr_mul_ui(y, x, u, rnd_mode) + mpfr_ptr y; + mpfr_srcptr x; + unsigned long int u; + mp_rnd_t rnd_mode; #endif { mp_limb_t carry, *my, *old_my; unsigned long c; @@ -59,8 +59,8 @@ mpfr_mul_ui() if (carry) count_leading_zeros(cnt, carry); else cnt=BITS_PER_MP_LIMB; - c = mpfr_round_raw(my, my, PREC(x), (SIGN(x)<0), - PREC(y)-BITS_PER_MP_LIMB+cnt, RND_MODE); + c = mpfr_round_raw(my, my, PREC(x), (MPFR_SIGN(x)<0), + PREC(y)-BITS_PER_MP_LIMB+cnt, rnd_mode); /* If cnt = 1111111111111 and c = 1 we shall get depressed */ if (c && (carry == (((mp_limb_t)1) << (BITS_PER_MP_LIMB - cnt)) - 1)) @@ -79,6 +79,6 @@ mpfr_mul_ui() EXP(y) = EXP(x) + BITS_PER_MP_LIMB - cnt; if (ysize < xsize) MPN_COPY(old_my, my, ysize); /* set sign */ - if (SIGN(y) != SIGN(x)) CHANGE_SIGN(y); + if (MPFR_SIGN(y) * MPFR_SIGN(x) < 0) CHANGE_SIGN(y); TMP_FREE(marker); } @@ -26,15 +26,15 @@ MA 02111-1307, USA. */ void #if __STDC__ -mpfr_neg(mpfr_ptr a, mpfr_srcptr b, unsigned char rnd_mode) +mpfr_neg(mpfr_ptr a, mpfr_srcptr b, mp_rnd_t rnd_mode) #else mpfr_neg(a, b, rnd_mode) - mpfr_ptr a; - mpfr_srcptr b; - unsigned char rnd_mode; + mpfr_ptr a; + mpfr_srcptr b; + mp_rnd_t rnd_mode; #endif { - if (a != b) mpfr_set4(a, b, rnd_mode, -SIGN(b)); + if (a != b) mpfr_set4(a, b, rnd_mode, -MPFR_SIGN(b)); else CHANGE_SIGN(a); return; } |