diff options
| -rw-r--r-- | Makefile.in | 101 | ||||
| -rw-r--r-- | add.c | 102 | ||||
| -rw-r--r-- | add_ui.c | 4 | ||||
| -rw-r--r-- | add_ulp.c | 18 | ||||
| -rw-r--r-- | agm.c | 12 | ||||
| -rw-r--r-- | clear.c | 2 | ||||
| -rw-r--r-- | cmp.c | 30 | ||||
| -rw-r--r-- | cmp_ui.c | 18 | ||||
| -rw-r--r-- | div.c | 40 | ||||
| -rw-r--r-- | div_2exp.c | 2 | ||||
| -rw-r--r-- | div_ui.c | 22 | ||||
| -rw-r--r-- | eq.c | 12 | ||||
| -rw-r--r-- | exp.c | 10 | ||||
| -rw-r--r-- | exp2.c | 26 | ||||
| -rw-r--r-- | exp3.c | 24 | ||||
| -rw-r--r-- | extract.c | 12 | ||||
| -rw-r--r-- | generic.c | 6 | ||||
| -rw-r--r-- | get_str.c | 20 | ||||
| -rw-r--r-- | init.c | 8 | ||||
| -rw-r--r-- | log.c | 14 | ||||
| -rw-r--r-- | log2.c | 10 | ||||
| -rw-r--r-- | mpfi.c | 14 | ||||
| -rw-r--r-- | mpfr.h | 32 | ||||
| -rw-r--r-- | mpfr.texi | 16 | ||||
| -rw-r--r-- | mpz_set_fr.c | 10 | ||||
| -rw-r--r-- | mul.c | 24 | ||||
| -rw-r--r-- | mul_2exp.c | 2 | ||||
| -rw-r--r-- | mul_ui.c | 22 | ||||
| -rw-r--r-- | neg.c | 2 | ||||
| -rw-r--r-- | out_str.c | 4 | ||||
| -rw-r--r-- | pi.c | 12 | ||||
| -rw-r--r-- | print_raw.c | 18 | ||||
| -rw-r--r-- | random.c | 8 | ||||
| -rw-r--r-- | random2.c | 6 | ||||
| -rw-r--r-- | reldiff.c | 4 | ||||
| -rw-r--r-- | round.c | 40 | ||||
| -rw-r--r-- | set.c | 12 | ||||
| -rw-r--r-- | set_d.c | 32 | ||||
| -rw-r--r-- | set_f.c | 12 | ||||
| -rw-r--r-- | set_prc_raw.c | 4 | ||||
| -rw-r--r-- | set_prec.c | 10 | ||||
| -rw-r--r-- | set_q.c | 6 | ||||
| -rw-r--r-- | set_si.c | 36 | ||||
| -rw-r--r-- | set_str.c | 6 | ||||
| -rw-r--r-- | set_str_raw.c | 4 | ||||
| -rw-r--r-- | set_z.c | 12 | ||||
| -rw-r--r-- | sin_cos.c | 24 | ||||
| -rw-r--r-- | sqrt.c | 44 | ||||
| -rw-r--r-- | sqrt_ui.c | 4 | ||||
| -rw-r--r-- | sub.c | 100 | ||||
| -rw-r--r-- | sub_ui.c | 2 | ||||
| -rw-r--r-- | trunc.c | 26 | ||||
| -rw-r--r-- | ui_div.c | 4 | ||||
| -rw-r--r-- | ui_sub.c | 2 | ||||
| -rw-r--r-- | urandomb.c | 6 | ||||
| -rw-r--r-- | zeta.c | 2 |
56 files changed, 554 insertions, 501 deletions
diff --git a/Makefile.in b/Makefile.in index eb202fec1..ce2a64db7 100644 --- a/Makefile.in +++ b/Makefile.in @@ -76,7 +76,7 @@ include_HEADERS = mpfr.h mpf2mpfr.h lib_LIBRARIES = libmpfr.a -libmpfr_a_SOURCES = cputime.h mpfr.h mpfr-impl.h sin_cos.c extract.c add.c add_ui.c add_ulp.c agm.c clear.c cmp.c cmp_ui.c div_2exp.c div.c div_ui.c dump.c eq.c exp2.c exp3.c exp.c get_str.c init.c inp_str.c karasqrt.c log2.c log.c mul_2exp.c mul.c mul_ui.c neg.c out_str.c pi.c pow.c print_raw.c print_rnd_mode.c random2.c random.c reldiff.c rnd_mode.c round.c set.c set_d.c set_dfl_prec.c set_rnd.c set_f.c set_prc_raw.c set_prec.c set_q.c set_si.c set_str.c set_str_raw.c set_z.c sqrt.c sqrt_ui.c sub.c sub_ui.c trunc.c ui_div.c ui_sub.c zeta.c urandomb.c +libmpfr_a_SOURCES = cputime.h mpfr.h mpfr-impl.h sin_cos.c extract.c add.c add_ui.c add_ulp.c agm.c clear.c cmp.c cmp_ui.c div_2exp.c div.c div_ui.c dump.c eq.c exp2.c exp3.c exp.c get_str.c init.c inp_str.c log2.c log.c mul_2exp.c mul.c mul_ui.c neg.c out_str.c pi.c pow.c print_raw.c print_rnd_mode.c random2.c random.c reldiff.c rnd_mode.c round.c set.c set_d.c set_dfl_prec.c set_rnd.c set_f.c set_prc_raw.c set_prec.c set_q.c set_si.c set_str.c set_str_raw.c set_z.c sqrt.c sqrt_ui.c sub.c sub_ui.c trunc.c ui_div.c ui_sub.c zeta.c urandomb.c sqrtrem.c mpz_set_fr.c libmpfr_a_LIBADD = ceil.o floor.o info_TEXINFOS = mpfr.texi @@ -89,15 +89,16 @@ LIBRARIES = $(lib_LIBRARIES) DEFS = @DEFS@ -I. -I$(srcdir) CPPFLAGS = @CPPFLAGS@ LIBS = @LIBS@ -libmpfr_a_DEPENDENCIES = ceil.o floor.o $(MPFR_OBJECTS) +libmpfr_a_DEPENDENCIES = ceil.o floor.o libmpfr_a_OBJECTS = sin_cos.o extract.o add.o add_ui.o add_ulp.o agm.o \ clear.o cmp.o cmp_ui.o div_2exp.o div.o div_ui.o dump.o eq.o exp2.o \ -exp3.o exp.o get_str.o init.o inp_str.o karasqrt.o log2.o log.o \ -mul_2exp.o mul.o mul_ui.o neg.o out_str.o pi.o pow.o print_raw.o \ -print_rnd_mode.o random2.o random.o reldiff.o rnd_mode.o round.o set.o \ -set_d.o set_dfl_prec.o set_rnd.o set_f.o set_prc_raw.o set_prec.o \ -set_q.o set_si.o set_str.o set_str_raw.o set_z.o sqrt.o sqrt_ui.o sub.o \ -sub_ui.o trunc.o ui_div.o ui_sub.o zeta.o urandomb.o +exp3.o exp.o get_str.o init.o inp_str.o log2.o log.o mul_2exp.o mul.o \ +mul_ui.o neg.o out_str.o pi.o pow.o print_raw.o print_rnd_mode.o \ +random2.o random.o reldiff.o rnd_mode.o round.o set.o set_d.o \ +set_dfl_prec.o set_rnd.o set_f.o set_prc_raw.o set_prec.o set_q.o \ +set_si.o set_str.o set_str_raw.o set_z.o sqrt.o sqrt_ui.o sub.o \ +sub_ui.o trunc.o ui_div.o ui_sub.o zeta.o urandomb.o sqrtrem.o \ +mpz_set_fr.o AR = ar COMPILE = $(CC) $(DEFS) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) CCLD = $(CC) @@ -117,6 +118,20 @@ DISTFILES = $(DIST_COMMON) $(SOURCES) $(HEADERS) $(TEXINFOS) $(EXTRA_DIST) TAR = gtar GZIP_ENV = --best +DEP_FILES = .deps/add.P .deps/add_ui.P .deps/add_ulp.P .deps/agm.P \ +.deps/clear.P .deps/cmp.P .deps/cmp_ui.P .deps/div.P .deps/div_2exp.P \ +.deps/div_ui.P .deps/dump.P .deps/eq.P .deps/exp.P .deps/exp2.P \ +.deps/exp3.P .deps/extract.P .deps/get_str.P .deps/init.P \ +.deps/inp_str.P .deps/log.P .deps/log2.P .deps/mpz_set_fr.P .deps/mul.P \ +.deps/mul_2exp.P .deps/mul_ui.P .deps/neg.P .deps/out_str.P .deps/pi.P \ +.deps/pow.P .deps/print_raw.P .deps/print_rnd_mode.P .deps/random.P \ +.deps/random2.P .deps/reldiff.P .deps/rnd_mode.P .deps/round.P \ +.deps/set.P .deps/set_d.P .deps/set_dfl_prec.P .deps/set_f.P \ +.deps/set_prc_raw.P .deps/set_prec.P .deps/set_q.P .deps/set_rnd.P \ +.deps/set_si.P .deps/set_str.P .deps/set_str_raw.P .deps/set_z.P \ +.deps/sin_cos.P .deps/sqrt.P .deps/sqrt_ui.P .deps/sqrtrem.P \ +.deps/sub.P .deps/sub_ui.P .deps/trunc.P .deps/ui_div.P .deps/ui_sub.P \ +.deps/urandomb.P .deps/zeta.P SOURCES = $(libmpfr_a_SOURCES) OBJECTS = $(libmpfr_a_OBJECTS) @@ -124,9 +139,9 @@ all: all-redirect .SUFFIXES: .SUFFIXES: .S .c .dvi .info .o .ps .s .texi .texinfo .txi $(srcdir)/Makefile.in: Makefile.am $(top_srcdir)/configure.in $(ACLOCAL_M4) - cd $(top_srcdir) && $(AUTOMAKE) --gnu --include-deps Makefile + cd $(top_srcdir) && $(AUTOMAKE) --gnu Makefile -Makefile: $(srcdir)/Makefile.in $(top_builddir)/config.status +Makefile: $(srcdir)/Makefile.in $(top_builddir)/config.status $(BUILT_SOURCES) cd $(top_builddir) \ && CONFIG_FILES=$@ CONFIG_HEADERS= $(SHELL) ./config.status @@ -170,9 +185,6 @@ uninstall-libLIBRARIES: rm -f $(DESTDIR)$(libdir)/$$p; \ done -.c.o: - $(COMPILE) -c $< - .s.o: $(COMPILE) -c $< @@ -452,6 +464,11 @@ distdir: $(DISTFILES) -rm -rf $(distdir) mkdir $(distdir) -chmod 777 $(distdir) + here=`cd $(top_builddir) && pwd`; \ + top_distdir=`cd $(distdir) && pwd`; \ + distdir=`cd $(distdir) && pwd`; \ + cd $(top_srcdir) \ + && $(AUTOMAKE) --include-deps --build-dir=$$here --srcdir-name=$(top_srcdir) --output-dir=$$top_distdir --gnu Makefile @for file in $(DISTFILES); do \ d=$(srcdir); \ if test -d $$d/$$file; then \ @@ -473,6 +490,38 @@ distdir: $(DISTFILES) fi; \ done $(MAKE) $(AM_MAKEFLAGS) top_distdir="$(top_distdir)" distdir="$(distdir)" dist-info + +DEPS_MAGIC := $(shell mkdir .deps > /dev/null 2>&1 || :) + +-include $(DEP_FILES) + +mostlyclean-depend: + +clean-depend: + +distclean-depend: + -rm -rf .deps + +maintainer-clean-depend: + +%.o: %.c + @echo '$(COMPILE) -c $<'; \ + $(COMPILE) -Wp,-MD,.deps/$(*F).pp -c $< + @-cp .deps/$(*F).pp .deps/$(*F).P; \ + tr ' ' '\012' < .deps/$(*F).pp \ + | sed -e 's/^\\$$//' -e '/^$$/ d' -e '/:$$/ d' -e 's/$$/ :/' \ + >> .deps/$(*F).P; \ + rm .deps/$(*F).pp + +%.lo: %.c + @echo '$(LTCOMPILE) -c $<'; \ + $(LTCOMPILE) -Wp,-MD,.deps/$(*F).pp -c $< + @-sed -e 's/^\([^:]*\)\.o[ ]*:/\1.lo \1.o :/' \ + < .deps/$(*F).pp > .deps/$(*F).P; \ + tr ' ' '\012' < .deps/$(*F).pp \ + | sed -e 's/^\\$$//' -e '/^$$/ d' -e '/:$$/ d' -e 's/$$/ :/' \ + >> .deps/$(*F).P; \ + rm -f .deps/$(*F).pp info-am: $(INFO_DEPS) info: info-recursive dvi-am: $(DVIS) @@ -513,25 +562,27 @@ distclean-generic: maintainer-clean-generic: mostlyclean-am: mostlyclean-libLIBRARIES mostlyclean-compile \ - mostlyclean-aminfo mostlyclean-tags mostlyclean-generic + mostlyclean-aminfo mostlyclean-tags mostlyclean-depend \ + mostlyclean-generic mostlyclean: mostlyclean-recursive clean-am: clean-libLIBRARIES clean-compile clean-aminfo clean-tags \ - clean-generic mostlyclean-am + clean-depend clean-generic mostlyclean-am clean: clean-recursive distclean-am: distclean-libLIBRARIES distclean-compile distclean-aminfo \ - distclean-tags distclean-generic clean-am + distclean-tags distclean-depend distclean-generic \ + clean-am distclean: distclean-recursive -rm -f config.status maintainer-clean-am: maintainer-clean-libLIBRARIES \ maintainer-clean-compile maintainer-clean-aminfo \ - maintainer-clean-tags maintainer-clean-generic \ - distclean-am + maintainer-clean-tags maintainer-clean-depend \ + maintainer-clean-generic distclean-am @echo "This command is intended for maintainers to use;" @echo "it deletes files that may require special tools to rebuild." @@ -549,12 +600,14 @@ uninstall-exec-recursive installdirs-recursive uninstalldirs-recursive \ all-recursive check-recursive installcheck-recursive info-recursive \ dvi-recursive mostlyclean-recursive distclean-recursive clean-recursive \ maintainer-clean-recursive tags tags-recursive mostlyclean-tags \ -distclean-tags clean-tags maintainer-clean-tags distdir info-am info \ -dvi-am dvi check check-am installcheck-am installcheck install-exec-am \ -install-exec install-data-am install-data install-am install \ -uninstall-am uninstall all-redirect all-am all installdirs-am \ -installdirs mostlyclean-generic distclean-generic clean-generic \ -maintainer-clean-generic clean mostlyclean distclean maintainer-clean +distclean-tags clean-tags maintainer-clean-tags distdir \ +mostlyclean-depend distclean-depend clean-depend \ +maintainer-clean-depend info-am info dvi-am dvi check check-am \ +installcheck-am installcheck install-exec-am install-exec \ +install-data-am install-data install-am install uninstall-am uninstall \ +all-redirect all-am all installdirs-am installdirs mostlyclean-generic \ +distclean-generic clean-generic maintainer-clean-generic clean \ +mostlyclean distclean maintainer-clean ceil.o: trunc.c mpfr.h @@ -52,46 +52,46 @@ mpfr_add1(a, b, c, rnd_mode, diff_exp) TMP_DECL(marker); TMP_MARK(marker); - ap = MANT(a); - bp = MANT(b); - cp = MANT(c); + ap = MPFR_MANT(a); + bp = MPFR_MANT(b); + cp = MPFR_MANT(c); if (ap == bp) { - bp = (mp_ptr) TMP_ALLOC(ABSSIZE(b) * BYTES_PER_MP_LIMB); - MPN_COPY (bp, ap, ABSSIZE(b)); + bp = (mp_ptr) TMP_ALLOC(MPFR_ABSSIZE(b) * BYTES_PER_MP_LIMB); + MPN_COPY (bp, ap, MPFR_ABSSIZE(b)); if (ap == cp) { cp = bp; } } else if (ap == cp) { - cp = (mp_ptr) TMP_ALLOC (ABSSIZE(c) * BYTES_PER_MP_LIMB); - MPN_COPY(cp, ap, ABSSIZE(c)); + cp = (mp_ptr) TMP_ALLOC (MPFR_ABSSIZE(c) * BYTES_PER_MP_LIMB); + MPN_COPY(cp, ap, MPFR_ABSSIZE(c)); } - an = (PREC(a)-1)/BITS_PER_MP_LIMB+1; /* number of significant limbs of a */ + an = (MPFR_PREC(a)-1)/BITS_PER_MP_LIMB+1; /* number of significant limbs of a */ - sh = an*BITS_PER_MP_LIMB-PREC(a); /* non-significant bits in low limb */ - bn = (PREC(b)-1)/BITS_PER_MP_LIMB + 1; /* number of significant limbs of b */ - cn = (PREC(c)-1)/BITS_PER_MP_LIMB + 1; - EXP(a) = EXP(b); + sh = an*BITS_PER_MP_LIMB-MPFR_PREC(a); /* non-significant bits in low limb */ + bn = (MPFR_PREC(b)-1)/BITS_PER_MP_LIMB + 1; /* number of significant limbs of b */ + cn = (MPFR_PREC(c)-1)/BITS_PER_MP_LIMB + 1; + MPFR_EXP(a) = MPFR_EXP(b); - if (MPFR_SIGN(a) * MPFR_SIGN(b) < 0) CHANGE_SIGN(a); + if (MPFR_SIGN(a) * MPFR_SIGN(b) < 0) MPFR_CHANGE_SIGN(a); /* case 1: diff_exp>=prec(a), i.e. c only affects the last bit through rounding */ - dif = PREC(a)-diff_exp; + dif = MPFR_PREC(a)-diff_exp; #ifdef DEBUG - printf("diff_exp=%u dif=PREC(a)-diff_exp=%d\n", diff_exp, dif); + printf("diff_exp=%u dif=MPFR_PREC(a)-diff_exp=%d\n", diff_exp, dif); printf("b= "); mpfr_print_raw(b); putchar('\n'); printf("c="); for (k=0;k<diff_exp;k++) putchar(' '); if (MPFR_SIGN(c)>0) putchar(' '); mpfr_print_raw(c); putchar('\n'); #endif if (dif<=0) { - /* diff_exp>=PREC(a): c does not overlap with a */ - /* either PREC(b)<=PREC(a), and we can copy the mantissa of b directly - into that of a, or PREC(b)>PREC(a) and we have to round b+c */ + /* diff_exp>=MPFR_PREC(a): c does not overlap with a */ + /* either MPFR_PREC(b)<=MPFR_PREC(a), and we can copy the mantissa of b directly + into that of a, or MPFR_PREC(b)>MPFR_PREC(a) and we have to round b+c */ - if (PREC(b)<=PREC(a)) { + if (MPFR_PREC(b)<=MPFR_PREC(a)) { MPN_COPY(ap+(an-bn), bp, bn); /* fill low significant limbs with zero */ @@ -102,22 +102,22 @@ mpfr_add1(a, b, c, rnd_mode, diff_exp) if (rnd_mode==GMP_RNDN) { /* to nearest */ - /* if diff_exp > PREC(a), no change */ + /* if diff_exp > MPFR_PREC(a), no change */ - if (diff_exp==PREC(a)) { + if (diff_exp==MPFR_PREC(a)) { /* if c is not zero, then as it is normalized, we have to add one to the lsb of a if c>1/2, or c=1/2 and lsb(a)=1 (round to even) */ - if (NOTZERO(c)) { + if (MPFR_NOTZERO(c)) { /* c is not zero */ /* check whether mant(c)=1/2 or not */ cc = *cp - (ONE<<(BITS_PER_MP_LIMB-1)); if (cc==0) { - bp = cp+(PREC(c)-1)/BITS_PER_MP_LIMB; + bp = cp+(MPFR_PREC(c)-1)/BITS_PER_MP_LIMB; while (cp<bp && cc==0) cc = *++cp; } @@ -126,18 +126,18 @@ mpfr_add1(a, b, c, rnd_mode, diff_exp) } } } - else if ((ISNONNEG(b) && rnd_mode==GMP_RNDU) || - (ISNEG(b) && rnd_mode==GMP_RNDD)) { + else if ((MPFR_ISNONNEG(b) && rnd_mode==GMP_RNDU) || + (MPFR_ISNEG(b) && rnd_mode==GMP_RNDD)) { /* round up */ - if (NOTZERO(c)) goto add_one_ulp; + if (MPFR_NOTZERO(c)) goto add_one_ulp; } /* in the other cases (round to zero, or up/down with sign -/+), nothing to do */ } else { - /* PREC(b)>PREC(a) : we have to round b+c */ + /* MPFR_PREC(b)>MPFR_PREC(a) : we have to round b+c */ k=bn-an; /* first copy the 'an' most significant limbs of b to a */ @@ -181,8 +181,8 @@ mpfr_add1(a, b, c, rnd_mode, diff_exp) } } if (cout==0) cout=(cc!=0); - sign = (ISNONNEG(b) && rnd_mode==GMP_RNDU) - || (ISNEG(b) && rnd_mode==GMP_RNDD) || (rnd_mode==GMP_RNDN); + sign = (MPFR_ISNONNEG(b) && rnd_mode==GMP_RNDU) + || (MPFR_ISNEG(b) && rnd_mode==GMP_RNDD) || (rnd_mode==GMP_RNDN); /* round towards infinity if dif=1, towards zero otherwise */ if ((sign==1) && (cout>0)) goto add_one_ulp; else if (rnd_mode==GMP_RNDN && cout==0 && (*ap & (ONE<<sh))) @@ -191,20 +191,20 @@ mpfr_add1(a, b, c, rnd_mode, diff_exp) } } else { - /* diff_exp < PREC(a) : c overlaps with a by dif bits */ + /* diff_exp < MPFR_PREC(a) : c overlaps with a by dif bits */ /* first copy upper part of c into a (after shift) */ unsigned char overlap; k = (dif-1)/BITS_PER_MP_LIMB + 1; /* only the highest k limbs from c have to be considered */ - cn = (PREC(c)-1)/BITS_PER_MP_LIMB + 1; + cn = (MPFR_PREC(c)-1)/BITS_PER_MP_LIMB + 1; MPN_ZERO(ap+k, an-k); /* do it now otherwise ap[k] may be destroyed in case dif<0 */ #ifdef DEBUG - printf("PREC(c)=%d\n", PREC(c)); + printf("MPFR_PREC(c)=%d\n", MPFR_PREC(c)); #endif - if (dif<=PREC(c)) { + if (dif<=MPFR_PREC(c)) { /* c has to be truncated */ dif = dif % BITS_PER_MP_LIMB; dif = (dif) ? BITS_PER_MP_LIMB-dif-sh : -sh; @@ -244,22 +244,22 @@ mpfr_add1(a, b, c, rnd_mode, diff_exp) /* then put high limbs to zero */ /* now add 'an' upper limbs of b in place */ - if (PREC(b)<=PREC(a)) { + if (MPFR_PREC(b)<=MPFR_PREC(a)) { overlap += 2; cc = mpn_add_n(ap+(an-bn), ap+(an-bn), bp, bn); } else - /* PREC(b) > PREC(a): we have to truncate b */ + /* MPFR_PREC(b) > MPFR_PREC(a): we have to truncate b */ cc = mpn_add_n(ap, ap, bp+(bn-an), an); if (cc) { /* shift one bit to the right */ - c3 = (ap[0]&1) && (PREC(a)%BITS_PER_MP_LIMB==0); + c3 = (ap[0]&1) && (MPFR_PREC(a)%BITS_PER_MP_LIMB==0); mpn_rshift(ap, ap, an, 1); ap[an-1] += ONE<<(BITS_PER_MP_LIMB-1); - EXP(a)++; + MPFR_EXP(a)++; } /* remains to do the rounding */ @@ -267,16 +267,16 @@ mpfr_add1(a, b, c, rnd_mode, diff_exp) #ifdef DEBUG printf("overlap=%d\n", overlap); #endif - if (rnd_mode==GMP_RNDN || (ISNONNEG(b) && rnd_mode==GMP_RNDU) - || (ISNEG(b) && rnd_mode==GMP_RNDD)) { + if (rnd_mode==GMP_RNDN || (MPFR_ISNONNEG(b) && rnd_mode==GMP_RNDU) + || (MPFR_ISNEG(b) && rnd_mode==GMP_RNDD)) { int kc; /* four cases: overlap = - (0) PREC(b) > PREC(a) and diff_exp+PREC(c) <= PREC(a) - (1) PREC(b) > PREC(a) and diff_exp+PREC(c) > PREC(a) - (2) PREC(b) <= PREC(a) and diff_exp+PREC(c) <= PREC(a) - (3) PREC(b) <= PREC(a) and diff_exp+PREC(c) > PREC(a) */ + (0) MPFR_PREC(b) > MPFR_PREC(a) and diff_exp+MPFR_PREC(c) <= MPFR_PREC(a) + (1) MPFR_PREC(b) > MPFR_PREC(a) and diff_exp+MPFR_PREC(c) > MPFR_PREC(a) + (2) MPFR_PREC(b) <= MPFR_PREC(a) and diff_exp+MPFR_PREC(c) <= MPFR_PREC(a) + (3) MPFR_PREC(b) <= MPFR_PREC(a) and diff_exp+MPFR_PREC(c) > MPFR_PREC(a) */ switch (overlap) { mp_limb_t cout; @@ -341,7 +341,7 @@ mpfr_add1(a, b, c, rnd_mode, diff_exp) */ /* c3=1 whenever b+c gave a carry out in most significant limb and the least significant bit (shifted right) was 1. - This can occur only when BITS_PER_MP_LIMB divides PREC(a), + This can occur only when BITS_PER_MP_LIMB divides MPFR_PREC(a), i.e. sh=0. */ if (sh) { @@ -371,7 +371,7 @@ mpfr_add1(a, b, c, rnd_mode, diff_exp) cc = mpn_add_1(ap, ap, an, ONE<<sh); if (cc) { ap[an-1] = (mp_limb_t)1 << (BITS_PER_MP_LIMB-1); - EXP(a)++; + MPFR_EXP(a)++; } end_of_add: @@ -392,14 +392,14 @@ mpfr_add(a, b, c, rnd_mode) { int diff_exp; - if (FLAG_NAN(b) || FLAG_NAN(c)) { - SET_NAN(a); return; + if (MPFR_IS_NAN(b) || MPFR_IS_NAN(c)) { + MPFR_SET_NAN(a); return; } - if (!NOTZERO(b)) { mpfr_set(a, c, rnd_mode); return; } - if (!NOTZERO(c)) { mpfr_set(a, b, rnd_mode); return; } + if (!MPFR_NOTZERO(b)) { mpfr_set(a, c, rnd_mode); return; } + if (!MPFR_NOTZERO(c)) { mpfr_set(a, b, rnd_mode); return; } - diff_exp = EXP(b)-EXP(c); + diff_exp = MPFR_EXP(b)-MPFR_EXP(c); 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); @@ -408,7 +408,7 @@ mpfr_add(a, b, c, rnd_mode) else { /* diff_exp=0 */ 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); + if (diff_exp==0) MPFR_SET_ZERO(a); 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); } @@ -25,7 +25,7 @@ MA 02111-1307, USA. */ #include "longlong.h" #include "mpfr.h" -#define MON_INIT(xp, x, p, s) PREC(x)=p; MANT(x)=xp; SIZE(x)=s; +#define MON_INIT(xp, x, p, s) MPFR_PREC(x)=p; MPFR_MANT(x)=xp; MPFR_SIZE(x)=s; void #if __STDC__ @@ -46,7 +46,7 @@ mpfr_add_ui(y, x, u, rnd_mode) MON_INIT(up, uu, BITS_PER_MP_LIMB, 1); count_leading_zeros(cnt, (mp_limb_t) u); *up = (mp_limb_t) u << cnt; - EXP(uu) = BITS_PER_MP_LIMB-cnt; + MPFR_EXP(uu) = BITS_PER_MP_LIMB-cnt; mpfr_add(y, x, uu, rnd_mode); } @@ -23,7 +23,7 @@ MA 02111-1307, USA. */ #include "gmp-impl.h" #include "mpfr.h" -/* sets x to x+sign(x)*2^(EXP(x)-PREC(x)) */ +/* sets x to x+sign(x)*2^(MPFR_EXP(x)-MPFR_PREC(x)) */ int #if __STDC__ mpfr_add_one_ulp(mpfr_ptr x) @@ -34,11 +34,11 @@ mpfr_add_one_ulp(x) { int xn, sh; mp_limb_t *xp; - xn = 1 + (PREC(x)-1)/BITS_PER_MP_LIMB; - sh = xn*BITS_PER_MP_LIMB - PREC(x); - xp = MANT(x); + xn = 1 + (MPFR_PREC(x)-1)/BITS_PER_MP_LIMB; + sh = xn*BITS_PER_MP_LIMB - MPFR_PREC(x); + xp = MPFR_MANT(x); if (mpn_add_1(xp, xp, xn, (mp_limb_t)1<<sh)) { - EXP(x)++; + MPFR_EXP(x)++; mpn_rshift(xp, xp, xn, 1); xp[xn-1] += (mp_limb_t)1<<(BITS_PER_MP_LIMB-1); } @@ -50,13 +50,13 @@ int mpfr_sub_one_ulp(mpfr_ptr x) { int xn, sh; mp_limb_t *xp; - xn = 1 + (PREC(x)-1)/BITS_PER_MP_LIMB; - sh = xn*BITS_PER_MP_LIMB-PREC(x); - xp = MANT(x); + xn = 1 + (MPFR_PREC(x)-1)/BITS_PER_MP_LIMB; + sh = xn*BITS_PER_MP_LIMB-MPFR_PREC(x); + xp = MPFR_MANT(x); mpn_sub_1(xp, xp, xn, (mp_limb_t)1<<sh); if (xp[xn-1] >> (BITS_PER_MP_LIMB-1) == 0) { /* was an exact power of two: not normalized any more */ - EXP(x)--; + MPFR_EXP(x)--; mpn_lshift(xp, xp, xn, 1); *xp |= ((mp_limb_t)1 << sh); } @@ -53,24 +53,24 @@ mpfr_agm(r, a, b, rnd_mode) /* If a or b is NaN, the result is NaN */ - if (FLAG_NAN(op1) || FLAG_NAN(op2)) - { SET_NAN(r); return; } + if (MPFR_IS_NAN(op1) || MPFR_IS_NAN(op2)) + { MPFR_SET_NAN(r); return; } /* If a or b is negative, the result is NaN */ if ((MPFR_SIGN(op1) < 0) || (MPFR_SIGN(op2) < 0)) - { SET_NAN(r); return; } + { MPFR_SET_NAN(r); return; } /* If a or b is 0, the result is 0 */ - if ((NOTZERO(op1) && NOTZERO(op2)) == 0) - { SET_ZERO(r); + if ((MPFR_NOTZERO(op1) && MPFR_NOTZERO(op2)) == 0) + { MPFR_SET_ZERO(r); return; } /* precision of the following calculus */ - q = PREC(r); + q = MPFR_PREC(r); p = q + 15; @@ -33,5 +33,5 @@ mpfr_clear (m) #endif { /* be careful to always free an entire number of limbs */ - (*_mp_free_func) (MANT(m), ABSSIZE(m) * BYTES_PER_MP_LIMB); + (*_mp_free_func) (MPFR_MANT(m), MPFR_ABSSIZE(m) * BYTES_PER_MP_LIMB); } @@ -53,26 +53,26 @@ mpfr_cmp3(b, c, s) unsigned long bn, cn; mp_limb_t *bp, *cp; - if (!NOTZERO(b)) { - if (!NOTZERO(c)) return 0; else return -(MPFR_SIGN(c)); + if (!MPFR_NOTZERO(b)) { + if (!MPFR_NOTZERO(c)) return 0; else return -(MPFR_SIGN(c)); } - else if (!NOTZERO(c)) return MPFR_SIGN(b); + else if (!MPFR_NOTZERO(c)) return MPFR_SIGN(b); s = s * MPFR_SIGN(b) * MPFR_SIGN(c); if (s<0) return(MPFR_SIGN(b)); /* now signs are equal */ - diff_exp = EXP(b)-EXP(c); + diff_exp = MPFR_EXP(b)-MPFR_EXP(c); s = (MPFR_SIGN(b) > 0) ? 1 : -1; if (diff_exp>0) return(s*(1+diff_exp)); else if (diff_exp<0) return(s*(-1+diff_exp)); /* both signs and exponents are equal */ - bn = (PREC(b)-1)/BITS_PER_MP_LIMB+1; - cn = (PREC(c)-1)/BITS_PER_MP_LIMB+1; - bp = MANT(b); cp = MANT(c); + bn = (MPFR_PREC(b)-1)/BITS_PER_MP_LIMB+1; + cn = (MPFR_PREC(c)-1)/BITS_PER_MP_LIMB+1; + bp = MPFR_MANT(b); cp = MPFR_MANT(c); while (bn && cn) { if (bp[--bn] != cp[--cn]) @@ -86,7 +86,7 @@ mpfr_cmp3(b, c, s) } /* returns the number of cancelled bits when one subtracts abs(c) from abs(b). - Assumes b>=c, which implies EXP(b)>=EXP(c). + Assumes b>=c, which implies MPFR_EXP(b)>=MPFR_EXP(c). if b=c, returns prec(b). */ int @@ -105,17 +105,17 @@ mpfr_cmp2(b, c) printf("b="); mpfr_print_raw(b); putchar('\n'); printf("c="); mpfr_print_raw(c); putchar('\n'); #endif - if (NOTZERO(c)==0) return (NOTZERO(b)) ? 0 : PREC(b); - d = EXP(b)-EXP(c); + if (MPFR_NOTZERO(c)==0) return (MPFR_NOTZERO(b)) ? 0 : MPFR_PREC(b); + d = MPFR_EXP(b)-MPFR_EXP(c); k = 0; /* result can be d or d+1 if d>1, or >= d otherwise */ /* k is the number of identical bits in the high part, then z is the number of possibly cancelled bits */ #ifdef DEBUG - if (d<0) { printf("assumption EXP(b)<EXP(c) violated\n"); exit(1); } + if (d<0) { printf("assumption MPFR_EXP(b)<MPFR_EXP(c) violated\n"); exit(1); } #endif - bn = (PREC(b)-1)/BITS_PER_MP_LIMB; - cn = (PREC(c)-1)/BITS_PER_MP_LIMB; - bp = MANT(b); cp = MANT(c); + bn = (MPFR_PREC(b)-1)/BITS_PER_MP_LIMB; + cn = (MPFR_PREC(c)-1)/BITS_PER_MP_LIMB; + bp = MPFR_MANT(b); cp = MPFR_MANT(c); /* subtract c from b from most significant to less significant limbs, and first determines first non zero limb difference */ if (d) @@ -133,7 +133,7 @@ mpfr_cmp2(b, c) while (bn>=0 && (cc=bp[bn--])==0) k+=BITS_PER_MP_LIMB; } /* now bn<0 or cc<>0 */ - if (cc==0 && bn<0) return(PREC(b)); + if (cc==0 && bn<0) return(MPFR_PREC(b)); } /* the first non-zero limb difference is cc, and the number @@ -44,11 +44,11 @@ mpfr_cmp_ui_2exp (b, i, f) if (MPFR_SIGN(b) < 0) return -1; /* now b>=0 */ - else if (!NOTZERO(b)) return((i) ? -1 : 0); + else if (!MPFR_NOTZERO(b)) return((i) ? -1 : 0); /* now b>0 */ else if (i==0) return 1; else { /* b>0, i>0 */ - e = EXP(b); /* 2^(e-1) <= b < 2^e */ + e = MPFR_EXP(b); /* 2^(e-1) <= b < 2^e */ if (e>f+BITS_PER_MP_LIMB) return 1; c = (mp_limb_t) i; @@ -58,8 +58,8 @@ mpfr_cmp_ui_2exp (b, i, f) /* now k=e */ c <<= (f+BITS_PER_MP_LIMB-k); - bn = (PREC(b)-1)/BITS_PER_MP_LIMB; - bp = MANT(b) + bn; + bn = (MPFR_PREC(b)-1)/BITS_PER_MP_LIMB; + bp = MPFR_MANT(b) + bn; if (*bp>c) return 1; else if (*bp<c) return -1; @@ -89,13 +89,13 @@ mpfr_cmp_si_2exp(b, i, f) si = (i<0) ? -1 : 1; /* sign of i */ if (MPFR_SIGN(b) * i < 0) return MPFR_SIGN(b); /* both signs differ */ - else if (!NOTZERO(b) || (i==0)) { /* one is zero */ - if (i==0) return ((NOTZERO(b)) ? MPFR_SIGN(b) : 0); + else if (!MPFR_NOTZERO(b) || (i==0)) { /* one is zero */ + if (i==0) return ((MPFR_NOTZERO(b)) ? MPFR_SIGN(b) : 0); else return si; /* b is zero */ } else { /* b and i are of same sign */ - e = EXP(b); /* 2^(e-1) <= b < 2^e */ + e = MPFR_EXP(b); /* 2^(e-1) <= b < 2^e */ if (e>f+BITS_PER_MP_LIMB) return si; c = (mp_limb_t) ((i<0) ? -i : i); @@ -105,8 +105,8 @@ mpfr_cmp_si_2exp(b, i, f) /* now k=e */ c <<= (f+BITS_PER_MP_LIMB-k); - bn = (PREC(b)-1)/BITS_PER_MP_LIMB; - bp = MANT(b) + bn; + bn = (MPFR_PREC(b)-1)/BITS_PER_MP_LIMB; + bp = MPFR_MANT(b) + bn; if (*bp>c) return si; else if (*bp<c) return -si; @@ -53,19 +53,19 @@ mpfr_div (r, u, v, rnd_mode) char can_round = 0; TMP_DECL (marker); - if (FLAG_NAN(u) || FLAG_NAN(v)) { SET_NAN(r); return; } + if (MPFR_IS_NAN(u) || MPFR_IS_NAN(v)) { MPFR_SET_NAN(r); return; } - usize = (PREC(u) - 1)/BITS_PER_MP_LIMB + 1; - vsize = (PREC(v) - 1)/BITS_PER_MP_LIMB + 1; + usize = (MPFR_PREC(u) - 1)/BITS_PER_MP_LIMB + 1; + vsize = (MPFR_PREC(v) - 1)/BITS_PER_MP_LIMB + 1; sign_quotient = ((MPFR_SIGN(u) * MPFR_SIGN(v) > 0) ? 1 : -1); - prec = PREC(r); + prec = MPFR_PREC(r); - if (!NOTZERO(u)) { SET_ZERO(r); return; } + if (!MPFR_NOTZERO(u)) { MPFR_SET_ZERO(r); return; } - if (!NOTZERO(v)) + if (!MPFR_NOTZERO(v)) vsize = 1 / v->_mp_d[vsize - 1]; /* Gestion des infinis ? */ - if (!NOTZERO(v)) + if (!MPFR_NOTZERO(v)) { r->_mp_exp = 0; MPN_ZERO(r->_mp_d, r->_mp_size); @@ -93,13 +93,13 @@ mpfr_div (r, u, v, rnd_mode) } else { mult = (mult < 0 ? 1 : 0); } - rsize = (PREC(r) + 3)/BITS_PER_MP_LIMB + 1; - rrsize = PREC(r)/BITS_PER_MP_LIMB + 1; + rsize = (MPFR_PREC(r) + 3)/BITS_PER_MP_LIMB + 1; + rrsize = MPFR_PREC(r)/BITS_PER_MP_LIMB + 1; /* Three extra bits are needed in order to get the quotient with enough precision ; take one extra bit for rrsize in order to solve more easily the problem of rounding to nearest. */ - /* ATTENTION, USIZE DOIT RESTER > A VSIZE !!!!!!!! */ + /* ATTENTION, UMPFR_SIZE DOIT RESTER > A VMPFR_SIZE !!!!!!!! */ do { @@ -181,10 +181,10 @@ mpfr_div (r, u, v, rnd_mode) } can_round = (mpfr_can_round_raw(rp, rrsize, sign_quotient, err, - GMP_RNDN, rnd_mode, PREC(r)) + GMP_RNDN, rnd_mode, MPFR_PREC(r)) || (usize == rsize && vsize == rsize && mpfr_can_round_raw(rp, rrsize, sign_quotient, err, - GMP_RNDZ, rnd_mode, PREC(r)))); + GMP_RNDZ, rnd_mode, MPFR_PREC(r)))); /* If we used all the limbs of both the dividend and the divisor, then we have the correct RNDZ rounding */ @@ -204,12 +204,12 @@ mpfr_div (r, u, v, rnd_mode) /* MAIS IL FAUT AJOUTER LE BOUT QUI MANQUE DE usize A rsize */ oldrrsize = rrsize; - rrsize = (PREC(r) - 1)/BITS_PER_MP_LIMB + 1; + rrsize = (MPFR_PREC(r) - 1)/BITS_PER_MP_LIMB + 1; if (can_round) { cc = mpfr_round_raw(rp, rp, err, (sign_quotient == -1 ? 1 : 0), - PREC(r), rnd_mode); + MPFR_PREC(r), rnd_mode); } else { /* Use the remainder to find out the correct rounding */ @@ -225,7 +225,7 @@ mpfr_div (r, u, v, rnd_mode) while (k >= 0) { if (tp[k]) break; k--; } if (k >= 0) { - t = PREC(r) & (BITS_PER_MP_LIMB - 1); + t = MPFR_PREC(r) & (BITS_PER_MP_LIMB - 1); if (t) { cc = mpn_add_1(rp, rp, rrsize, @@ -239,7 +239,7 @@ mpfr_div (r, u, v, rnd_mode) else if (rnd_mode == GMP_RNDN) /* even rounding */ { - rw = (PREC(r) + 1) & (BITS_PER_MP_LIMB - 1); + rw = (MPFR_PREC(r) + 1) & (BITS_PER_MP_LIMB - 1); if (rw) { rw = BITS_PER_MP_LIMB - rw; nw = 0; } else nw = 1; if ((rw ? (rp[nw] >> (rw + 1)) & 1 : (rp[nw] >> (BITS_PER_MP_LIMB - 1)) & 1)) @@ -258,7 +258,7 @@ mpfr_div (r, u, v, rnd_mode) } - if (sign_quotient * MPFR_SIGN(r) < 0) { CHANGE_SIGN(r); } + if (sign_quotient * MPFR_SIGN(r) < 0) { MPFR_CHANGE_SIGN(r); } r->_mp_exp = rexp; if (cc) { @@ -267,9 +267,9 @@ mpfr_div (r, u, v, rnd_mode) r->_mp_exp++; } - rw = rrsize * BITS_PER_MP_LIMB - PREC(r); - MPN_COPY(MANT(r), rp, rrsize); - MANT(r)[0] &= ~(((mp_limb_t)1 << rw) - 1); + rw = rrsize * BITS_PER_MP_LIMB - MPFR_PREC(r); + MPN_COPY(MPFR_MANT(r), rp, rrsize); + MPFR_MANT(r)[0] &= ~(((mp_limb_t)1 << rw) - 1); TMP_FREE (marker); } diff --git a/div_2exp.c b/div_2exp.c index 0af678f7c..3264a5adf 100644 --- a/div_2exp.c +++ b/div_2exp.c @@ -37,7 +37,7 @@ mpfr_div_2exp(y, x, n, rnd_mode) { /* Important particular case */ if (y != x) mpfr_set(y, x, rnd_mode); - EXP(y) -= n; + MPFR_EXP(y) -= n; return; } @@ -49,17 +49,17 @@ mpfr_div_ui(y, x, u, rnd_mode) int xn, yn, dif, sh, i; mp_limb_t *xp, *yp, *tmp, c, d; TMP_DECL(marker); - if (FLAG_NAN(x)) { SET_NAN(y); return 1; } + if (MPFR_IS_NAN(x)) { MPFR_SET_NAN(y); return 1; } if (u==0) { fprintf(stderr, "division by zero\n"); exit(1); } TMP_MARK(marker); - xn = (PREC(x)-1)/BITS_PER_MP_LIMB + 1; - yn = (PREC(y)-1)/BITS_PER_MP_LIMB + 1; + xn = (MPFR_PREC(x)-1)/BITS_PER_MP_LIMB + 1; + yn = (MPFR_PREC(y)-1)/BITS_PER_MP_LIMB + 1; - xp = MANT(x); - yp = MANT(y); - EXP(y) = EXP(x); - if (MPFR_SIGN(x) * MPFR_SIGN(y) < 0) CHANGE_SIGN(y); + xp = MPFR_MANT(x); + yp = MPFR_MANT(y); + MPFR_EXP(y) = MPFR_EXP(x); + if (MPFR_SIGN(x) * MPFR_SIGN(y) < 0) MPFR_CHANGE_SIGN(y); dif = yn+1-xn; #ifdef DEBUG @@ -77,27 +77,27 @@ mpfr_div_ui(y, x, u, rnd_mode) /* patch for bug in mpn_divrem_1 for GMP 2.xxx */ count_leading_zeros(sh, c); c <<= sh; - EXP(y) += sh; + MPFR_EXP(y) += sh; #endif c = mpn_divrem_1(tmp, dif, xp, xn, c); } else /* dif < 0 i.e. xn > yn */ c = mpn_divrem_1(tmp, 0, xp-dif, yn, c); - if (tmp[yn]==0) { tmp--; sh=0; EXP(y) -= BITS_PER_MP_LIMB; } + if (tmp[yn]==0) { tmp--; sh=0; MPFR_EXP(y) -= BITS_PER_MP_LIMB; } /* shift left to normalize */ count_leading_zeros(sh, tmp[yn]); if (sh) { mpn_lshift(yp, tmp+1, yn, sh); yp[0] += tmp[0] >> (BITS_PER_MP_LIMB-sh); - EXP(y) -= sh; + MPFR_EXP(y) -= sh; } else MPN_COPY(yp, tmp+1, yn); #ifdef DEBUG printf("y="); mpfr_print_raw(y); putchar('\n'); #endif - sh = yn*BITS_PER_MP_LIMB - PREC(y); + sh = yn*BITS_PER_MP_LIMB - MPFR_PREC(y); /* it remains sh bits in less significant limb of y */ d = *yp & (((mp_limb_t)1 << sh) - 1); @@ -44,17 +44,17 @@ mpfr_eq (u, v, n_bits) uexp = u->_mp_exp; vexp = v->_mp_exp; - usize = (PREC(u)-1)/BITS_PER_MP_LIMB + 1; - vsize = (PREC(v)-1)/BITS_PER_MP_LIMB + 1; + usize = (MPFR_PREC(u)-1)/BITS_PER_MP_LIMB + 1; + vsize = (MPFR_PREC(v)-1)/BITS_PER_MP_LIMB + 1; /* 1. Are the signs different? */ if (MPFR_SIGN(u) == MPFR_SIGN(v)) { /* U and V are both non-negative or both negative. */ - if (!NOTZERO(u)) - return !NOTZERO(v); - if (!NOTZERO(v)) - return !NOTZERO(u); + if (!MPFR_NOTZERO(u)) + return !MPFR_NOTZERO(v); + if (!MPFR_NOTZERO(v)) + return !MPFR_NOTZERO(u); /* Fall out. */ } @@ -46,16 +46,16 @@ mpfr_exp(y, x, rnd_mode) int n, expx, K, precy, q, k, l, expr, err; mpfr_t r, s, t; - if (FLAG_NAN(x)) { SET_NAN(y); return 1; } - if (!NOTZERO(x)) { mpfr_set_ui(y, 1, GMP_RNDN); return 0; } + if (MPFR_IS_NAN(x)) { MPFR_SET_NAN(y); return 1; } + if (!MPFR_NOTZERO(x)) { mpfr_set_ui(y, 1, GMP_RNDN); return 0; } - expx = EXP(x); - precy = PREC(y); + expx = MPFR_EXP(x); + precy = MPFR_PREC(y); /* if x > (2^31-1)*ln(2), then exp(x) > 2^(2^31-1) i.e. gives +infinity */ if (expx > 30) { if (MPFR_SIGN(x)>0) { printf("+infinity"); return 1; } - else { SET_ZERO(y); return 1; } + else { MPFR_SET_ZERO(y); return 1; } } /* if x < 2^(-precy), then exp(x) i.e. gives 1 +/- 1 ulp(1) */ @@ -105,19 +105,19 @@ mpfr_exp2(y, x, rnd_mode) mpfr_t r, s, t; mpz_t ss; TMP_DECL(marker); - if (FLAG_NAN(x)) { SET_NAN(y); return 1; } - if (!NOTZERO(x)) { mpfr_set_ui(y, 1, GMP_RNDN); return 0; } + if (MPFR_IS_NAN(x)) { MPFR_SET_NAN(y); return 1; } + if (!MPFR_NOTZERO(x)) { mpfr_set_ui(y, 1, GMP_RNDN); return 0; } - expx = EXP(x); - precy = PREC(y); + expx = MPFR_EXP(x); + precy = MPFR_PREC(y); #ifdef DEBUG - printf("EXP(x)=%d\n",expx); + printf("MPFR_EXP(x)=%d\n",expx); #endif /* if x > (2^31-1)*ln(2), then exp(x) > 2^(2^31-1) i.e. gives +infinity */ if (expx > 30) { if (MPFR_SIGN(x)>0) { printf("+infinity"); return 1; } - else { SET_ZERO(y); return 1; } + else { MPFR_SET_ZERO(y); return 1; } } /* if x < 2^(-precy), then exp(x) i.e. gives 1 +/- 1 ulp(1) */ @@ -141,7 +141,7 @@ mpfr_exp2(y, x, rnd_mode) q = precy + err + K + 3; mpfr_init2(r, q); mpfr_init2(s, q); mpfr_init2(t, q); /* the algorithm consists in computing an upper bound of exp(x) using - a precision of q bits, and see if we can round to PREC(y) taking + a precision of q bits, and see if we can round to MPFR_PREC(y) taking into account the maximal error. Otherwise we increase q. */ do { #ifdef DEBUG @@ -193,7 +193,7 @@ mpfr_exp2(y, x, rnd_mode) mpz_mul(ss, ss, ss); exps <<= 1; exps += mpz_normalize(ss, ss, q); } - mpfr_set_z(s, ss, GMP_RNDN); EXP(s) += exps; + mpfr_set_z(s, ss, GMP_RNDN); MPFR_EXP(s) += exps; if (n>0) mpfr_mul_2exp(s, s, n, GMP_RNDU); else mpfr_div_2exp(s, s, -n, GMP_RNDU); @@ -206,7 +206,7 @@ mpfr_exp2(y, x, rnd_mode) K += k; #ifdef DEBUG printf("after mult. by 2^n:\n"); - if (EXP(s)>-1024) printf("s=%1.20e\n",mpfr_get_d(s)); + if (MPFR_EXP(s)>-1024) printf("s=%1.20e\n",mpfr_get_d(s)); printf(" ="); mpfr_print_raw(s); putchar('\n'); printf("err=%d bits\n", K); #endif @@ -229,7 +229,7 @@ mpfr_exp2(y, x, rnd_mode) return 1; } -/* s <- 1 + r/1! + r^2/2! + ... + r^l/l! while EXP(r^l/l!)+EXPR(r)>-q +/* s <- 1 + r/1! + r^2/2! + ... + r^l/l! while MPFR_EXP(r^l/l!)+MPFR_EXPR(r)>-q using naive method with O(l) multiplications. Return the number of iterations l. The absolute error on s is less than 3*l*(l+1)*2^(-q). @@ -286,7 +286,7 @@ mpfr_exp2_aux(s, r, q, exps) return l; } -/* s <- 1 + r/1! + r^2/2! + ... + r^l/l! while EXP(r^l/l!)+EXPR(r)>-q +/* s <- 1 + r/1! + r^2/2! + ... + r^l/l! while MPFR_EXP(r^l/l!)+MPFR_EXPR(r)>-q using Brent/Kung method with O(sqrt(l)) multiplications. Return l. Uses m multiplications of full size and 2l/m of decreasing size, @@ -311,12 +311,12 @@ mpfr_exp2_aux2(s, r, q, exps) TMP_DECL(marker); /* estimate value of l */ - l = q / (-EXP(r)); + l = q / (-MPFR_EXP(r)); m = (int) sqrt((double) l); TMP_MARK(marker); R = (mpz_t*) TMP_ALLOC((m+1)*sizeof(mpz_t)); /* R[i] stands for r^i */ expR = (int*) TMP_ALLOC((m+1)*sizeof(int)); /* exponent for R[i] */ - sizer = 1 + (PREC(r)-1)/BITS_PER_MP_LIMB; + sizer = 1 + (MPFR_PREC(r)-1)/BITS_PER_MP_LIMB; mpz_init(tmp); MY_INIT_MPZ(rr, sizer+2); MY_INIT_MPZ(t, 2*sizer); /* double size for products */ @@ -58,7 +58,7 @@ int m; mpz_t* P,*S; mpz_t* ptoj; int diff,expo; - int precy = PREC(y); + int precy = MPFR_PREC(y); int * mult; int prec_i_have; int *nb_terms; @@ -131,7 +131,7 @@ int m; mpz_tdiv_q(S[0], S[0], P[0]); mpfr_set_z(y,S[0], GMP_RNDD); - EXP(y) += expo; + MPFR_EXP(y) += expo; mpfr_div_2exp(y, y, r*(i-1),GMP_RNDN); for (i=0;i<=m;i++) { mpz_clear(P[i]); mpz_clear(S[i]); mpz_clear(ptoj[i]); } @@ -173,8 +173,8 @@ mp_rnd_t rnd_mode; int iter; int logn; /* commencons par 0 */ - if (FLAG_NAN(x)) { SET_NAN(y); return 1; } - if (!NOTZERO(x)) { + if (MPFR_IS_NAN(x)) { MPFR_SET_NAN(y); return 1; } + if (!MPFR_NOTZERO(x)) { mpfr_set_ui(y, 1, GMP_RNDN); return 0; } @@ -182,23 +182,23 @@ mp_rnd_t rnd_mode; /* on commence par ecrire x = 1.xxxxxxxxxxxxx ----- k bits -- */ prec_x = (int) ceil(log - ((double) (PREC(x)) / (double) BITS_PER_MP_LIMB) + ((double) (MPFR_PREC(x)) / (double) BITS_PER_MP_LIMB) /log(2.0)); logn = (int) ceil(log - ((double) prec_x+PREC(y)) + ((double) prec_x+MPFR_PREC(y)) /log(2.0)); if (logn < 2) logn = 2; - ttt = EXP(x); - mpfr_init2(x_copy,PREC(x)); + ttt = MPFR_EXP(x); + mpfr_init2(x_copy,MPFR_PREC(x)); mpfr_set(x_copy,x,GMP_RNDD); /* on fait le shift pour que le nombre soit inferieur a 1 */ if (ttt > 0) { shift_x = ttt; mpfr_mul_2exp(x_copy,x,-ttt, GMP_RNDN); - ttt = EXP(x_copy); + ttt = MPFR_EXP(x_copy); } - realprec = PREC(y)+logn; + realprec = MPFR_PREC(y)+logn; while (!good){ Prec = realprec+shift+2+shift_x; k = (int) ceil(log @@ -243,7 +243,7 @@ mp_rnd_t rnd_mode; #endif #ifdef DEBUG fprintf(stderr, "fin\n"); - mpfr_out_str(stderr, 2, PREC(y), t, GMP_RNDD); + mpfr_out_str(stderr, 2, MPFR_PREC(y), t, GMP_RNDD); fprintf(stderr, "\n ii --- ii \n"); #endif twopoweri <<= 1; @@ -252,7 +252,7 @@ mp_rnd_t rnd_mode; for (loop= 0 ; loop < shift_x; loop++) mpfr_mul(tmp,tmp,tmp,GMP_RNDD); mpfr_clear(t); - if (mpfr_can_round(tmp, realprec, GMP_RNDD, rnd_mode, PREC(y))){ + if (mpfr_can_round(tmp, realprec, GMP_RNDD, rnd_mode, MPFR_PREC(y))){ mpfr_set(y,tmp,rnd_mode); mpfr_clear(tmp); good = 1; @@ -23,30 +23,30 @@ int i; int size; int j; - if (ABSSIZE(p) < two_i){ + if (MPFR_ABSSIZE(p) < two_i){ int j; y->_mp_alloc=two_i_2 ; y->_mp_size=two_i_2 ; PTR(y) = (*_mp_allocate_func)(two_i_2 * sizeof(mp_limb_t)); MPN_ZERO (PTR(y), two_i_2); assert(y->_mp_d!=NULL); - for(j= 0; j < ABSSIZE(p) - two_i_2 ; j++){ - y->_mp_d[j + two_i - ABSSIZE(p)] = p->_mp_d[j]; + for(j= 0; j < MPFR_ABSSIZE(p) - two_i_2 ; j++){ + y->_mp_d[j + two_i - MPFR_ABSSIZE(p)] = p->_mp_d[j]; } } else { PTR(y) = (*_mp_allocate_func)(two_i_2 * sizeof(mp_limb_t)); - memcpy(y -> _mp_d, p->_mp_d+ABSSIZE(p) - two_i, two_i_2 * sizeof(mp_limb_t)); + memcpy(y -> _mp_d, p->_mp_d+MPFR_ABSSIZE(p) - two_i, two_i_2 * sizeof(mp_limb_t)); y->_mp_alloc=two_i_2 ; y->_mp_size=two_i_2 ; } - size = ABSSIZE(y); + size = MPFR_ABSSIZE(y); for (j = 0; j < size; j++) { if (y->_mp_d[j]) { - if ISNEG(p) + if MPFR_ISNEG(p) mpz_neg(y,y); return 0; } @@ -54,7 +54,7 @@ int m; mpfr_t tmp; #endif int diff,expo; - int precy = PREC(y); + int precy = MPFR_PREC(y); n = 1 << m; P = (mpz_t*) (*_mp_allocate_func) ((m+1) * sizeof(mpz_t)); S = (mpz_t*) (*_mp_allocate_func) ((m+1) * sizeof(mpz_t)); @@ -178,12 +178,12 @@ int m; mpz_tdiv_q(S[0], S[0], P[0]); mpfr_set_z(y,S[0], GMP_RNDD); - EXP(y) += expo; + MPFR_EXP(y) += expo; #ifdef R_IS_RATIONAL /* division exacte */ mpz_div_ui(qtoj[m], qtoj[m], r); - i = (PREC(y)); + i = (MPFR_PREC(y)); mpfr_init2(tmp,i); mpfr_set_z(tmp, qtoj[m] , GMP_RNDD); mpfr_div(y, y, tmp,GMP_RNDD); @@ -61,7 +61,7 @@ char *mpfr_get_str(str, expptr, base, n, op, rnd_mode) neg = (MPFR_SIGN(op)<0) ? 1 : 0; - if (!NOTZERO(op)) { + if (!MPFR_NOTZERO(op)) { if (str==NULL) str = (*_mp_allocate_func)(neg + n + 1); str0 = str; if (MPFR_SIGN(op)<0) *str++ = '-'; @@ -77,7 +77,7 @@ char *mpfr_get_str(str, expptr, base, n, op, rnd_mode) /* if pow2 <> 0, then base = 2^pow2 */ /* first determines the exponent */ - e = EXP(op); + e = MPFR_EXP(op); d = fabs(mpfr_get_d2(op, 0)); /* the absolute value of op is between 1/2*2^e and 2^e */ /* the output exponent f is such that base^(f-1) <= |op| < base^f @@ -88,11 +88,11 @@ char *mpfr_get_str(str, expptr, base, n, op, rnd_mode) /* performs exact rounding, i.e. returns y such that for GMP_RNDU for example, we have: x*2^(e-p) <= y*base^(f-n) */ - n = (int) ((double)PREC(op)*log(2.0)/log((double)base)); + n = (int) ((double)MPFR_PREC(op)*log(2.0)/log((double)base)); if (n==0) n=1; } #ifdef DEBUG - printf("f=%d n=%d EXP(op)=%d PREC(op)=%d\n", f, n, e, PREC(op)); + printf("f=%d n=%d MPFR_EXP(op)=%d MPFR_PREC(op)=%d\n", f, n, e, MPFR_PREC(op)); #endif /* now the first n digits of the mantissa are obtained from rnd(op*base^(n-f)) */ @@ -142,7 +142,7 @@ char *mpfr_get_str(str, expptr, base, n, op, rnd_mode) mpfr_mul(b, op, a, rnd_mode); } } - if (neg) CHANGE_SIGN(b); /* put b positive */ + if (neg) MPFR_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); @@ -163,17 +163,17 @@ char *mpfr_get_str(str, expptr, base, n, op, rnd_mode) case GMP_RNDD: rnd_mode=GMP_RNDU; break; } - if (ok) mpfr_round(b, rnd_mode, EXP(b)); + if (ok) mpfr_round(b, rnd_mode, MPFR_EXP(b)); - prec=EXP(b); /* may have changed due to rounding */ + prec=MPFR_EXP(b); /* may have changed due to rounding */ /* now the mantissa is the integer part of b */ mpz_init(bz); q=1+(prec-1)/BITS_PER_MP_LIMB; _mpz_realloc(bz, q); sh = prec%BITS_PER_MP_LIMB; - e = 1 + (PREC(b)-1)/BITS_PER_MP_LIMB-q; - if (sh) mpn_rshift(PTR(bz), MANT(b)+e, q, BITS_PER_MP_LIMB-sh); - else MPN_COPY(PTR(bz), MANT(b)+e, q); + e = 1 + (MPFR_PREC(b)-1)/BITS_PER_MP_LIMB-q; + if (sh) mpn_rshift(PTR(bz), MPFR_MANT(b)+e, q, BITS_PER_MP_LIMB-sh); + else MPN_COPY(PTR(bz), MPFR_MANT(b)+e, q); bz->_mp_size=q; /* computes the number of characters needed */ @@ -41,10 +41,10 @@ mpfr_init2 (x, p) xsize = (p - 1)/BITS_PER_MP_LIMB + 1; - PREC(x) = p; - MANT(x) = (mp_ptr) (*_mp_allocate_func) (xsize * BYTES_PER_MP_LIMB); - SIZE(x) = xsize; - EXP(x) = 0; /* avoids uninitialized memory reads for zero */ + MPFR_PREC(x) = p; + MPFR_MANT(x) = (mp_ptr) (*_mp_allocate_func) (xsize * BYTES_PER_MP_LIMB); + MPFR_SIZE(x) = xsize; + MPFR_EXP(x) = 0; /* avoids uninitialized memory reads for zero */ } void @@ -63,16 +63,16 @@ mpfr_log(r, a, rnd_mode) TMP_DECL(marker); /* If a is NaN or a is negative or null, the result is NaN */ - if (FLAG_NAN(a) || (NOTZERO(a)==0) || (MPFR_SIGN(a)<0)) - { SET_NAN(r); return 1; } + if (MPFR_IS_NAN(a) || (MPFR_NOTZERO(a)==0) || (MPFR_SIGN(a)<0)) + { MPFR_SET_NAN(r); return 1; } /* If a is 1, the result is 0 */ if (mpfr_cmp_ui_2exp(a,1,0)==0){ - SET_ZERO(r); + MPFR_SET_ZERO(r); return 0; /* only case where the result is exact */ } - q=PREC(r); + q=MPFR_PREC(r); ref=mpfr_get_d(a)-1.0; if (ref<0) @@ -92,7 +92,7 @@ mpfr_log(r, a, rnd_mode) printf("p=%d\n", p); #endif /* Calculus of m (depends on p) */ - m=(int) ceil(((double) p)/2.0) -EXP(a)+1; + m=(int) ceil(((double) p)/2.0) -MPFR_EXP(a)+1; /* All the mpfr_t needed have a precision of p */ TMP_MARK(marker); @@ -116,9 +116,9 @@ mpfr_log(r, a, rnd_mode) mpfr_div(tmp2,cst,tmp1,GMP_RNDN); /* pi/2*AG(1,4/s), err<=5ulps */ 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); + cancel = MPFR_EXP(tmp2); mpfr_sub(cst,tmp2,tmp1,GMP_RNDN); /* log(a), err<=7ulps+cancel */ - cancel -= EXP(cst); + cancel -= MPFR_EXP(cst); #ifdef DEBUG printf("cancelled bits=%d\n", cancel); printf("approx="); mpfr_print_raw(cst); putchar('\n'); @@ -63,12 +63,12 @@ mpfr_const_aux_log2(mylog, rnd_mode) mpfr_ptr mylog; mp_rnd_t rnd_mode; mpz_t cst; int good = 0; int logn; - int prec_i_want = PREC(mylog); + int prec_i_want = MPFR_PREC(mylog); int prec_x; mpz_init(cst); logn = (int) ceil(log - ((double) PREC(mylog)) + ((double) MPFR_PREC(mylog)) /log(2.0)); prec_x = prec_i_want + logn; while (!good){ @@ -112,7 +112,7 @@ mpfr_const_aux_log2(mylog, rnd_mode) mpfr_ptr mylog; mp_rnd_t rnd_mode; return 0; } -/* set x to log(2) rounded to precision PREC(x) with direction rnd_mode +/* set x to log(2) rounded to precision MPFR_PREC(x) with direction rnd_mode use formula log(2) = sum(1/k/2^k, k=1..infinity) @@ -134,7 +134,7 @@ mpfr_const_log2(x, rnd_mode) mpfr_ptr x; mp_rnd_t rnd_mode; { int N, oldN, k, precx; mpz_t s, t, u; - precx = PREC(x); + precx = MPFR_PREC(x); /* has stored value enough precision ? */ if (precx <= __mpfr_const_log2_prec) { @@ -174,7 +174,7 @@ mpfr_const_log2(x, rnd_mode) mpfr_ptr x; mp_rnd_t rnd_mode; } #endif mpfr_set_z(x, s, rnd_mode); - EXP(x) -= N; + MPFR_EXP(x) -= N; mpz_clear(s); mpz_clear(t); mpz_clear(u); } else { @@ -28,7 +28,7 @@ MA 02111-1307, USA. */ #define MPFR_SIGN_PART(x) mpfr_cmp_ui_2exp(x,0,0) #define MPFI_ISPOS(x) ((MPFR_SIGN_PART((&(x->left)))>=0) && (MPFR_SIGN_PART((&(x->right)))>0)) -#define MPFI_ISNEG(x) ((MPFR_SIGN_PART((&(x->left)))<0) && (MPFR_SIGN_PART((&(x->right)))<=0)) +#define MPFI_MPFR_ISNEG(x) ((MPFR_SIGN_PART((&(x->left)))<0) && (MPFR_SIGN_PART((&(x->right)))<=0)) #define MPFI_ISNULL(x) ((MPFR_SIGN_PART((&(x->left)))==0) && (MPFR_SIGN_PART((&(x->right)))==0)) #define MPFI_HASZERO(x) ((MPFR_SIGN_PART((&(x->left)))<0) && (MPFR_SIGN_PART((&(x->right)))>0)) @@ -200,11 +200,11 @@ void mpfi_mul (mpfi_ptr a, mpfi_srcptr b, mpfi_srcptr c) mpfr_mul(&(a->right), &(b->right), &(c->right), MPFI_RNDU); } else { - if (MPFI_ISNEG(c)) { + if (MPFI_MPFR_ISNEG(c)) { in_place = (b->right)._mp_d == (a->right)._mp_d; if (!in_place) u[0] = b->right; else { - mpfr_init2 (u, PREC(&(b->right))); + mpfr_init2 (u, MPFR_PREC(&(b->right))); mpfr_set (u, &(b->right), GMP_RNDD); } mpfr_mul (&(a->right), &(b->left), &(c->right), MPFI_RNDU); @@ -240,7 +240,7 @@ void mpfi_div (mpfi_ptr a, mpfi_srcptr u, mpfi_srcptr c) mpfr_div(&(a->right),&(b->right),&(c->left),MPFI_RNDU); } else { - if (MPFI_ISNEG(c)) { + if (MPFI_MPFR_ISNEG(c)) { mpfr_div(&(a->right),&(b->left),&(c->left),MPFI_RNDU); mpfr_div(&(a->left),&(b->right),&(c->right),MPFI_RNDD); } @@ -289,11 +289,11 @@ void mpfi_sub_ui (mpfi_ptr a, mpfi_srcptr b, unsigned int c) void mpfi_ui_div (mpfi_ptr a, unsigned int b, mpfi_srcptr c) { - if (MPFI_ISPOS(c) || MPFI_ISNEG(c)) { + if (MPFI_ISPOS(c) || MPFI_MPFR_ISNEG(c)) { mpfr_t tmp; int in_place = (a->left)._mp_d == (c->left)._mp_d; if (in_place) { - mpfr_init2 (tmp, PREC(&(a->left))); + mpfr_init2 (tmp, MPFR_PREC(&(a->left))); mpfr_set (tmp, &(a->left), GMP_RNDN); } else tmp[0] = a->left; @@ -350,7 +350,7 @@ void mpfi_neg(mpfi_ptr a, mpfi_srcptr b) { mpfr_t tmp; - mpfr_init2 (tmp, PREC(&(b->left))); + mpfr_init2 (tmp, MPFR_PREC(&(b->left))); mpfr_set (tmp, &(b->left), MPFI_RNDD); mpfr_neg (&(a->left), &(b->right), MPFI_RNDD); mpfr_neg (&(a->right), tmp, MPFI_RNDU); @@ -47,7 +47,7 @@ typedef struct { /* 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; /* ABSSIZE(.) is the number of allocated + 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 @@ -89,26 +89,26 @@ typedef __gmp_const __mpfr_struct *mpfr_srcptr; /* bit 31 of _mp_size is used for sign, bit 30 of _mp_size is used for Nan flag, remaining bits are used to store the number of allocated limbs */ -#define FLAG_NAN(x) (((x)->_mp_size >> 30)&1) -#define SET_NAN(x) ((x)->_mp_size |= (1<<30)) -#define ABSSIZE(x) ((x)->_mp_size & ((1<<30)-1)) -#define SIZE(x) ((x)->_mp_size) -#define EXP(x) ((x)->_mp_exp) -#define MANT(x) ((x)->_mp_d) +#define MPFR_IS_NAN(x) (((x)->_mp_size >> 30)&1) +#define MPFR_SET_NAN(x) ((x)->_mp_size |= (1<<30)) +#define MPFR_ABSSIZE(x) ((x)->_mp_size & ((1<<30)-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 ISNONNEG(x) (MPFR_SIGN(x)>=0) -#define ISNEG(x) (MPFR_SIGN(x)==-1) -#define CHANGE_SIGN(x) (SIZE(x) = SIZE(x) ^ (((mp_size_t)1)<<31)) -#define PREC(x) ((x)->_mp_prec) -#define NOTZERO(x) (MANT(x)[(PREC(x)-1)/BITS_PER_MP_LIMB]) -#define SET_ZERO(x) (MANT(x)[(PREC(x)-1)/BITS_PER_MP_LIMB] = 0) -#define mpfr_sgn(x) ((NOTZERO(x)) ? MPFR_SIGN(x) : 0) +#define MPFR_ISNONNEG(x) (MPFR_SIGN(x)>=0) +#define MPFR_ISNEG(x) (MPFR_SIGN(x)==-1) +#define MPFR_CHANGE_SIGN(x) (MPFR_SIZE(x) = MPFR_SIZE(x) ^ (((mp_size_t)1)<<31)) +#define MPFR_PREC(x) ((x)->_mp_prec) +#define MPFR_NOTZERO(x) (MPFR_MANT(x)[(MPFR_PREC(x)-1)/BITS_PER_MP_LIMB]) +#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, (PREC(x)-1)>>3, (q+7)>>3); \ - PREC(x) = q; } + ((x)->_mp_d, (MPFR_PREC(x)-1)>>3, (q+7)>>3); \ + MPFR_PREC(x) = q; } void mpfr_init2 _PROTO ((mpfr_ptr, mp_prec_t)); void mpfr_init _PROTO ((mpfr_ptr)); @@ -905,7 +905,7 @@ therefore the @var{SIGN} macro does not distinguish zero from non-zero numbers. @end deftypefn -@deftypefn Macro int NOTZERO (mpfr_t @var{op}) +@deftypefn Macro int MPFR_NOTZERO (mpfr_t @var{op}) Returns zero when its argument is zero, and a non-zero value otherwise. @end deftypefn @@ -1072,31 +1072,31 @@ Negative random numbers are generated when @var{max_size} is negative. These functions were mainly designed for the implementation of @code{mpfr}, but may be useful for users too. -@deftypefun Macro int FLAG_NAN (mpfr_t x) +@deftypefun Macro int MPFR_IS_NAN (mpfr_t x) Returns a non-zero value iff its argument is @samp{Not a Number}. @end deftypefun -@deftypefun Macro int SET_NAN (mpfr_t x) +@deftypefun Macro int MPFR_SET_NAN (mpfr_t x) Sets its argument to @samp{Not a Number}. @end deftypefun -@deftypefun Macro int ABSSIZE (mpfr_t x) +@deftypefun Macro int MPFR_ABSSIZE (mpfr_t x) Returns the number of limb the mantissa of @var{x} points to. @end deftypefun -@deftypefun Macro int EXP (mpfr_t x) +@deftypefun Macro int MPFR_EXP (mpfr_t x) Returns the exponent of @var{x}. @end deftypefun -@deftypefun Macro mp_limb_t* MANT (mpfr_t x) +@deftypefun Macro mp_limb_t* MPFR_MANT (mpfr_t x) Returns a pointer to the mantissa of @var{x}. @end deftypefun -@deftypefun Macro mp_prec_t PREC (mpfr_t x) +@deftypefun Macro mp_prec_t MPFR_PREC (mpfr_t x) Returns the precision of @var{x}. @end deftypefun -@deftypefun Macro int CHANGE_SIGN (mpfr_t x) +@deftypefun Macro int MPFR_CHANGE_SIGN (mpfr_t x) Changes the sign of @var{x}. @end deftypefun diff --git a/mpz_set_fr.c b/mpz_set_fr.c index 7f4d79c4d..7c01399af 100644 --- a/mpz_set_fr.c +++ b/mpz_set_fr.c @@ -37,7 +37,7 @@ mpz_set_fr (z, f) { int fn, sh; - fn = 1 + (PREC(f)-1)/BITS_PER_MP_LIMB; + fn = 1 + (MPFR_PREC(f)-1)/BITS_PER_MP_LIMB; /* check whether allocated space for z is enough */ if (ALLOC(z) < fn) { @@ -46,11 +46,11 @@ mpz_set_fr (z, f) ALLOC(z) = fn; } - sh = fn*BITS_PER_MP_LIMB - PREC(f); - if (sh) mpn_rshift(PTR(z), MANT(f), fn, sh); - else MPN_COPY(PTR(z), MANT(f), fn); + sh = fn*BITS_PER_MP_LIMB - MPFR_PREC(f); + if (sh) mpn_rshift(PTR(z), MPFR_MANT(f), fn, sh); + else MPN_COPY(PTR(z), MPFR_MANT(f), fn); SIZ(z) = fn; - return EXP(f)-PREC(f); + return MPFR_EXP(f)-MPFR_PREC(f); } @@ -25,8 +25,8 @@ MA 02111-1307, USA. */ #include "mpfr.h" /* Remains to do: -- do not use all bits of b and c when PREC(b)>PREC(a) or PREC(c)>PREC(a) - [current complexity is O(PREC(b)*PREC(c))] +- do not use all bits of b and c when MPFR_PREC(b)>MPFR_PREC(a) or MPFR_PREC(c)>MPFR_PREC(a) + [current complexity is O(MPFR_PREC(b)*MPFR_PREC(c))] */ void @@ -41,18 +41,18 @@ mpfr_mul(a, b, c, rnd_mode) #endif { unsigned int bn, cn, an, tn, k; int cc; - mp_limb_t *ap=MANT(a), *bp=MANT(b), *cp=MANT(c), *tmp, b1; + mp_limb_t *ap=MPFR_MANT(a), *bp=MPFR_MANT(b), *cp=MPFR_MANT(c), *tmp, b1; long int sign_product; TMP_DECL(marker); /* deal with NaN and zero */ - if (FLAG_NAN(b) || FLAG_NAN(c)) { SET_NAN(a); return; } - if (!NOTZERO(b) || !NOTZERO(c)) { SET_ZERO(a); return; } + if (MPFR_IS_NAN(b) || MPFR_IS_NAN(c)) { MPFR_SET_NAN(a); return; } + if (!MPFR_NOTZERO(b) || !MPFR_NOTZERO(c)) { MPFR_SET_ZERO(a); return; } sign_product = MPFR_SIGN(b) * MPFR_SIGN(c); - bn = (PREC(b)-1)/BITS_PER_MP_LIMB+1; /* number of significant limbs of b */ - cn = (PREC(c)-1)/BITS_PER_MP_LIMB+1; /* number of significant limbs of c */ - tn = (PREC(c)+PREC(b)-1)/BITS_PER_MP_LIMB+1; + bn = (MPFR_PREC(b)-1)/BITS_PER_MP_LIMB+1; /* number of significant limbs of b */ + cn = (MPFR_PREC(c)-1)/BITS_PER_MP_LIMB+1; /* number of significant limbs of c */ + tn = (MPFR_PREC(c)+MPFR_PREC(b)-1)/BITS_PER_MP_LIMB+1; k = bn+cn; /* effective nb of limbs used by b*c */ TMP_MARK(marker); tmp = (mp_limb_t*) TMP_ALLOC(k*BYTES_PER_MP_LIMB); @@ -62,17 +62,17 @@ mpfr_mul(a, b, c, rnd_mode) /* now tmp[0]..tmp[k-1] contains the product of both mantissa, with tmp[k-1]>=2^(BITS_PER_MP_LIMB-2) */ - an = (PREC(a)-1)/BITS_PER_MP_LIMB+1; /* number of significant limbs of a */ + an = (MPFR_PREC(a)-1)/BITS_PER_MP_LIMB+1; /* number of significant limbs of a */ b1 >>= BITS_PER_MP_LIMB-1; /* msb from the product */ if (b1==0) mpn_lshift(tmp, tmp, k, 1); cc = mpfr_round_raw(ap, tmp+bn+cn-tn, - PREC(b)+PREC(c), (sign_product<0), PREC(a), rnd_mode); + MPFR_PREC(b)+MPFR_PREC(c), (sign_product<0), MPFR_PREC(a), rnd_mode); if (cc) { /* cc = 1 ==> result is a power of two */ ap[an-1] = (mp_limb_t) 1 << (BITS_PER_MP_LIMB-1); } - EXP(a) = EXP(b) + EXP(c) + b1 - 1 + cc; - if (sign_product * MPFR_SIGN(a)<0) CHANGE_SIGN(a); + MPFR_EXP(a) = MPFR_EXP(b) + MPFR_EXP(c) + b1 - 1 + cc; + if (sign_product * MPFR_SIGN(a)<0) MPFR_CHANGE_SIGN(a); TMP_FREE(marker); return; } diff --git a/mul_2exp.c b/mul_2exp.c index 3420c96e0..4e0d4df28 100644 --- a/mul_2exp.c +++ b/mul_2exp.c @@ -37,7 +37,7 @@ mpfr_mul_2exp(y, x, n, rnd_mode) { /* Important particular case */ if (y != x) mpfr_set(y, x, rnd_mode); - EXP(y) += n; + MPFR_EXP(y) += n; return; } @@ -41,9 +41,9 @@ mpfr_mul_ui(y, x, u, rnd_mode) TMP_DECL(marker); TMP_MARK(marker); - my = MANT(y); ex = EXP(x); - ysize = (PREC(y)-1)/BITS_PER_MP_LIMB + 1; - xsize = (PREC(x)-1)/BITS_PER_MP_LIMB + 1; + my = MPFR_MANT(y); ex = MPFR_EXP(x); + ysize = (MPFR_PREC(y)-1)/BITS_PER_MP_LIMB + 1; + xsize = (MPFR_PREC(x)-1)/BITS_PER_MP_LIMB + 1; old_my = my; @@ -53,7 +53,7 @@ mpfr_mul_ui(y, x, u, rnd_mode) } else dif=ysize-xsize; - carry = mpn_mul_1(my+dif, MANT(x), xsize, u); + carry = mpn_mul_1(my+dif, MPFR_MANT(x), xsize, u); MPN_ZERO(my, dif); /* WARNING: count_leading_zeros is undefined for carry=0 */ @@ -62,12 +62,12 @@ mpfr_mul_ui(y, x, u, rnd_mode) /* Warning: the number of limbs used by x and the lower part of y may differ */ - sh = (PREC(x)+BITS_PER_MP_LIMB-1)/BITS_PER_MP_LIMB - - (PREC(y)+cnt-1)/BITS_PER_MP_LIMB; + sh = (MPFR_PREC(x)+BITS_PER_MP_LIMB-1)/BITS_PER_MP_LIMB + - (MPFR_PREC(y)+cnt-1)/BITS_PER_MP_LIMB; /* Warning: if all significant bits are in the carry, one has to be careful */ - if (cnt + PREC(y) < BITS_PER_MP_LIMB) + if (cnt + MPFR_PREC(y) < BITS_PER_MP_LIMB) { /* Quick 'n dirty */ @@ -86,8 +86,8 @@ mpfr_mul_ui(y, x, u, rnd_mode) carry = 0; cnt = BITS_PER_MP_LIMB; } - c = mpfr_round_raw(my+sh, my, PREC(x), (MPFR_SIGN(x)<0), - PREC(y)-BITS_PER_MP_LIMB+cnt, rnd_mode); + c = mpfr_round_raw(my+sh, my, MPFR_PREC(x), (MPFR_SIGN(x)<0), + MPFR_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) << (cnt ? BITS_PER_MP_LIMB - cnt : 0)) @@ -104,9 +104,9 @@ mpfr_mul_ui(y, x, u, rnd_mode) mpn_rshift(my, my, ysize, BITS_PER_MP_LIMB - cnt); my[ysize - 1] |= (carry << cnt); } - EXP(y) = ex + BITS_PER_MP_LIMB - cnt; + MPFR_EXP(y) = ex + BITS_PER_MP_LIMB - cnt; if (ysize < xsize) MPN_COPY(old_my, my, ysize); /* set sign */ - if (MPFR_SIGN(y) * MPFR_SIGN(x) < 0) CHANGE_SIGN(y); + if (MPFR_SIGN(y) * MPFR_SIGN(x) < 0) MPFR_CHANGE_SIGN(y); TMP_FREE(marker); } @@ -35,6 +35,6 @@ mpfr_neg(a, b, rnd_mode) #endif { if (a != b) mpfr_set4(a, b, rnd_mode, -MPFR_SIGN(b)); - else CHANGE_SIGN(a); + else MPFR_CHANGE_SIGN(a); return; } @@ -40,8 +40,8 @@ mpfr_out_str (stream, base, n_digits, op, rnd_mode) { char *s,*s0; size_t l; mp_exp_t e; - if (FLAG_NAN(op)) { fprintf(stream, "NaN"); return 3; } - if (!NOTZERO(op)) { fprintf(stream, "0"); return 1; } + if (MPFR_IS_NAN(op)) { fprintf(stream, "NaN"); return 3; } + if (!MPFR_NOTZERO(op)) { fprintf(stream, "0"); return 1; } s = mpfr_get_str(NULL, &e, base, n_digits, op, rnd_mode); s0 = s; @@ -49,12 +49,12 @@ mpfr_pi_machin3(mylog, rnd_mode) #endif { int prec, logn, prec_x; - int prec_i_want=PREC(mylog); + int prec_i_want=MPFR_PREC(mylog); int good = 0; mpfr_t tmp1, tmp2, result,tmp3,tmp4,tmp5,tmp6; mpz_t cst; logn = (int) ceil(log - ((double) PREC(mylog)) + ((double) MPFR_PREC(mylog)) /log(2.0)); prec_x = prec_i_want + logn + 5; mpz_init(cst); @@ -123,7 +123,7 @@ mpfr_pi_machin3(mylog, rnd_mode) } /* -Set x to the value of Pi to precision PREC(x) rounded to direction rnd_mode. +Set x to the value of Pi to precision MPFR_PREC(x) rounded to direction rnd_mode. Use the formula giving the binary representation of Pi found by Simon Plouffe and the Borwein's brothers: @@ -167,7 +167,7 @@ mpfr_const_pi(x, rnd_mode) { int N, oldN, n, prec; mpz_t pi, num, den, d3, d2, tmp; mpfr_t y; - prec=PREC(x); + prec=MPFR_PREC(x); /* has stored value enough precision ? */ if ((prec==__mpfr_const_pi_prec && rnd_mode==__mpfr_const_pi_rnd) || @@ -197,7 +197,7 @@ mpfr_const_pi(x, rnd_mode) */ for (n=0; n<N; n++) { /* num(n)-num(n-1) = 240*n+31 */ - mpz_add_ui(num, num, 240*n+31); /* no overflow up to PREC=71M */ + mpz_add_ui(num, num, 240*n+31); /* no overflow up to MPFR_PREC=71M */ /* d2(n) - d2(n-1) = 12288*(n-1) */ if (n>0) mpz_add_ui(d2, d2, 12288*(n-1)); else mpz_sub_ui(d2, d2, 12288); @@ -216,7 +216,7 @@ mpfr_const_pi(x, rnd_mode) if (mpfr_cmp(x, y) != 0) { fprintf(stderr, "does not converge\n"); exit(1); } - EXP(x) -= 4*N; + MPFR_EXP(x) -= 4*N; mpz_clear(pi); mpz_clear(num); mpz_clear(den); mpz_clear(d3); mpz_clear(d2); mpz_clear(tmp); mpfr_clear(y); } else diff --git a/print_raw.c b/print_raw.c index b6e680a8e..77821a61c 100644 --- a/print_raw.c +++ b/print_raw.c @@ -36,9 +36,9 @@ mpfr_get_str_raw(digit_ptr, x) { mp_limb_t *mx, wd, t; long ex, sx, k, l, p; - mx = MANT(x); - ex = EXP(x); - p = PREC(x); + mx = MPFR_MANT(x); + ex = MPFR_EXP(x); + p = MPFR_PREC(x); if (MPFR_SIGN(x) < 0) { *digit_ptr = '-'; digit_ptr++; } sprintf(digit_ptr, "0."); digit_ptr += 2; @@ -71,21 +71,21 @@ mpfr_print_raw(x) { char *str; - if (FLAG_NAN(x)) printf("NaN"); - else if (!NOTZERO(x)) printf("0"); + if (MPFR_IS_NAN(x)) printf("NaN"); + else if (!MPFR_NOTZERO(x)) printf("0"); else { /* 3 char for sign + 0 + binary point - + ABSSIZE(x) * BITS_PER_MP_LIMB for mantissa + + MPFR_ABSSIZE(x) * BITS_PER_MP_LIMB for mantissa + 2 for brackets in mantissa + 1 for 'E' + 11 for exponent (including sign) - = 17 + ABSSIZE(x) * BITS_PER_MP_LIMB + = 17 + MPFR_ABSSIZE(x) * BITS_PER_MP_LIMB */ - str = (char *) (*_mp_allocate_func) ((17 + ABSSIZE(x) * BITS_PER_MP_LIMB)*sizeof(char)); + str = (char *) (*_mp_allocate_func) ((17 + MPFR_ABSSIZE(x) * BITS_PER_MP_LIMB)*sizeof(char)); mpfr_get_str_raw(str, x); printf("%s", str); - (*_mp_free_func) (str, (17 + ABSSIZE(x) * BITS_PER_MP_LIMB)*sizeof(char)); + (*_mp_free_func) (str, (17 + MPFR_ABSSIZE(x) * BITS_PER_MP_LIMB)*sizeof(char)); } } @@ -25,7 +25,7 @@ MA 02111-1307, USA. */ #include "longlong.h" #include "mpfr.h" -/* Computes a random mpfr in [0, 1[ with precision PREC */ +/* Computes a random mpfr in [0, 1[ with precision MPFR_PREC */ extern long random _PROTO((void)); extern int srandom _PROTO((unsigned int)); @@ -45,9 +45,9 @@ mpfr_random(x) mpfr_ptr x; #endif { - mp_limb_t *xp; unsigned long xn, i, cnt, prec=PREC(x); + mp_limb_t *xp; unsigned long xn, i, cnt, prec=MPFR_PREC(x); - xp = MANT(x); + xp = MPFR_MANT(x); xn = (prec-1)/BITS_PER_MP_LIMB + 1; for (i = 0; i < xn; i++) @@ -60,7 +60,7 @@ mpfr_random(x) count_leading_zeros(cnt, xp[xn - 1]); if (cnt) mpn_lshift(xp, xp, xn, cnt); - EXP(x) = -cnt; + MPFR_EXP(x) = -cnt; cnt = xn*BITS_PER_MP_LIMB - prec; /* cnt is the number of non significant bits in the low limb */ @@ -49,8 +49,8 @@ mpfr_random2 (x, size, exp) mp_exp_t exp; #endif { - mp_size_t xn; unsigned long cnt; mp_ptr xp = MANT(x); - mp_size_t prec = (PREC(x) - 1)/BITS_PER_MP_LIMB; + mp_size_t xn; unsigned long cnt; mp_ptr xp = MPFR_MANT(x); + mp_size_t prec = (MPFR_PREC(x) - 1)/BITS_PER_MP_LIMB; xn = ABS (size); if (xn != 0) @@ -66,7 +66,7 @@ mpfr_random2 (x, size, exp) count_leading_zeros(cnt, xp[xn - 1]); if (cnt) mpn_lshift(xp, xp, xn, cnt); - EXP(x) = exp-cnt; + MPFR_EXP(x) = exp-cnt; cnt = xn*BITS_PER_MP_LIMB - prec; /* cnt is the number of non significant bits in the low limb */ xp[0] &= ~((((mp_limb_t)1)<<cnt) - 1); @@ -36,9 +36,9 @@ mpfr_reldiff(a, b, c, rnd_mode) mp_rnd_t rnd_mode; #endif { - if (FLAG_NAN(b) || FLAG_NAN(c)) { SET_NAN(a); return; } + if (MPFR_IS_NAN(b) || MPFR_IS_NAN(c)) { MPFR_SET_NAN(a); return; } - if (!NOTZERO(b)) /* reldiff = abs(c)/c = sign(c) */ + if (!MPFR_NOTZERO(b)) /* reldiff = abs(c)/c = sign(c) */ mpfr_set_ui(a, MPFR_SIGN(c), rnd_mode); else { @@ -189,28 +189,28 @@ mpfr_round(x, rnd_mode, prec) signx = MPFR_SIGN(x); /* check if x has enough allocated space for the mantissa */ - if (nw > ABSSIZE(x)) { - MANT(x) = (mp_ptr) (*_mp_reallocate_func) - (MANT(x), ABSSIZE(x)*BYTES_PER_MP_LIMB, nw * BYTES_PER_MP_LIMB); - SIZE(x) = nw; /* new number of allocated limbs */ + if (nw > MPFR_ABSSIZE(x)) { + MPFR_MANT(x) = (mp_ptr) (*_mp_reallocate_func) + (MPFR_MANT(x), MPFR_ABSSIZE(x)*BYTES_PER_MP_LIMB, nw * BYTES_PER_MP_LIMB); + MPFR_SIZE(x) = nw; /* new number of allocated limbs */ } TMP_MARK(marker); tmp = TMP_ALLOC (nw * BYTES_PER_MP_LIMB); - carry = mpfr_round_raw(tmp, MANT(x), PREC(x), (MPFR_SIGN(x)<0), prec, + carry = mpfr_round_raw(tmp, MPFR_MANT(x), MPFR_PREC(x), (MPFR_SIGN(x)<0), prec, rnd_mode); if (carry) { mpn_rshift(tmp, tmp, nw, 1); tmp [nw-1] |= (((mp_limb_t)1) << (BITS_PER_MP_LIMB - 1)); - EXP(x)++; + MPFR_EXP(x)++; } - SIZE(x)=nw; - if (signx * MPFR_SIGN(x)<0) CHANGE_SIGN(x); - PREC(x) = prec; - MPN_COPY(MANT(x), tmp, nw); + MPFR_SIZE(x)=nw; + if (signx * MPFR_SIGN(x)<0) MPFR_CHANGE_SIGN(x); + MPFR_PREC(x) = prec; + MPN_COPY(MPFR_MANT(x), tmp, nw); TMP_FREE(marker); } @@ -219,7 +219,7 @@ mpfr_round(x, rnd_mode, prec) /* assuming b is an approximation of x in direction rnd1 - with error at most 2^(EXP(b)-err), returns 1 if one is + with error at most 2^(MPFR_EXP(b)-err), returns 1 if one is able to round exactly x to precision prec with direction rnd2, and 0 otherwise. @@ -239,7 +239,7 @@ mpfr_can_round(b, err, rnd1, rnd2, prec) mp_prec_t prec; #endif { - return mpfr_can_round_raw(MANT(b), (PREC(b) - 1)/BITS_PER_MP_LIMB + 1, + return mpfr_can_round_raw(MPFR_MANT(b), (MPFR_PREC(b) - 1)/BITS_PER_MP_LIMB + 1, MPFR_SIGN(b), err, rnd1, rnd2, prec); } @@ -287,12 +287,12 @@ mpfr_can_round_raw(bp, bn, neg, err, rnd1, rnd2, prec) switch (rnd1) { - case GMP_RNDZ: /* b <= x <= b+2^(EXP(b)-err) */ + case GMP_RNDZ: /* b <= x <= b+2^(MPFR_EXP(b)-err) */ tmp = TMP_ALLOC(tn*BYTES_PER_MP_LIMB); cc = (bp[bn-1]>>l1) & 1; cc ^= mpfr_round_raw2(bp, bn, neg, rnd2, prec); - /* now round b+2^(EXP(b)-err) */ + /* now round b+2^(MPFR_EXP(b)-err) */ cc2 = mpn_add_1(tmp+bn-k, bp+bn-k, k, (mp_limb_t)1<<l); /* if carry, then all bits up to err were to 1, and we can round only if cc==0 and mpfr_round_raw2 returns 0 below */ @@ -303,13 +303,13 @@ mpfr_can_round_raw(bp, bn, neg, err, rnd1, rnd2, prec) TMP_FREE(marker); return (cc == cc2); - case GMP_RNDU: /* b-2^(EXP(b)-err) <= x <= b */ + case GMP_RNDU: /* b-2^(MPFR_EXP(b)-err) <= x <= b */ tmp = TMP_ALLOC(tn*BYTES_PER_MP_LIMB); /* first round b */ cc = (bp[bn-1]>>l1) & 1; cc ^= mpfr_round_raw2(bp, bn, neg, rnd2, prec); - /* now round b-2^(EXP(b)-err) */ + /* now round b-2^(MPFR_EXP(b)-err) */ cc2 = mpn_sub_1(tmp+bn-k, bp+bn-k, k, (mp_limb_t)1<<l); /* if borrow, then all bits up to err were to 0, and we can round only if cc==0 and mpfr_round_raw2 returns 1 below */ @@ -320,12 +320,12 @@ mpfr_can_round_raw(bp, bn, neg, err, rnd1, rnd2, prec) TMP_FREE(marker); return (cc == cc2); - case GMP_RNDN: /* b-2^(EXP(b)-err-1) <= x <= b+2^(EXP(b)-err-1) */ + case GMP_RNDN: /* b-2^(MPFR_EXP(b)-err-1) <= x <= b+2^(MPFR_EXP(b)-err-1) */ if (l==0) tn++; tmp = TMP_ALLOC(tn*BYTES_PER_MP_LIMB); /* this case is the same than GMP_RNDZ, except we first have to - subtract 2^(EXP(b)-err-1) from b */ + subtract 2^(MPFR_EXP(b)-err-1) from b */ if (l) { l--; /* tn=bn */ @@ -337,14 +337,14 @@ mpfr_can_round_raw(bp, bn, neg, err, rnd1, rnd2, prec) mpn_sub_1(tmp+tn-k, tmp+tn-k, k, (mp_limb_t)1<<l); } - /* round b-2^(EXP(b)-err-1) */ + /* round b-2^(MPFR_EXP(b)-err-1) */ /* we can disregard borrow, since we start from tmp in 2nd case too */ cc = (tmp[tn-1]>>l1) & 1; cc ^= mpfr_round_raw2(tmp, tn, neg, rnd2, prec); if (l==BITS_PER_MP_LIMB-1) { l=0; k--; } else l++; - /* round b+2^(EXP(b)-err-1) = b-2^(EXP(b)-err-1) + 2^(EXP(b)-err) */ + /* round b+2^(MPFR_EXP(b)-err-1) = b-2^(MPFR_EXP(b)-err-1) + 2^(MPFR_EXP(b)-err) */ cc2 = mpn_add_1(tmp+tn-k, tmp+tn-k, k, (mp_limb_t)1<<l); /* if carry, then all bits up to err were to 1, and we can round only if cc==0 and mpfr_round_raw2 returns 0 below */ @@ -35,12 +35,12 @@ mpfr_set4(a, b, rnd_mode, signb) int signb; #endif { - int carry, an, preca = PREC(a), sh; mp_limb_t *ap = MANT(a); + int carry, an, preca = MPFR_PREC(a), sh; mp_limb_t *ap = MPFR_MANT(a); - if (FLAG_NAN(b)) { SET_NAN(a); return; } + if (MPFR_IS_NAN(b)) { MPFR_SET_NAN(a); return; } - carry = mpfr_round_raw(ap, MANT(b), PREC(b), (signb<0), preca, rnd_mode); - EXP(a) = EXP(b); + carry = mpfr_round_raw(ap, MPFR_MANT(b), MPFR_PREC(b), (signb<0), preca, rnd_mode); + MPFR_EXP(a) = MPFR_EXP(b); if (carry) { an = (preca-1)/BITS_PER_MP_LIMB + 1; sh = an * BITS_PER_MP_LIMB - preca; @@ -49,7 +49,7 @@ mpfr_set4(a, b, rnd_mode, signb) } mpn_rshift(ap, ap, an, 1); ap[an-1] |= (mp_limb_t) 1 << (BITS_PER_MP_LIMB-1); - EXP(a)++; + MPFR_EXP(a)++; } - if (MPFR_SIGN(a) * signb < 0) CHANGE_SIGN(a); + if (MPFR_SIGN(a) * signb < 0) MPFR_CHANGE_SIGN(a); } @@ -245,29 +245,29 @@ mpfr_set_d(r, d, rnd_mode) { int signd, sizer; unsigned int cnt; - if (d == 0) { SET_ZERO(r); return; } - else if (isnan(d)) { SET_NAN(r); return; } + if (d == 0) { MPFR_SET_ZERO(r); return; } + else if (isnan(d)) { MPFR_SET_NAN(r); return; } signd = (d < 0) ? -1 : 1; d = ABS (d); - sizer = (PREC(r)-1)/BITS_PER_MP_LIMB + 1; + sizer = (MPFR_PREC(r)-1)/BITS_PER_MP_LIMB + 1; /* warning: __mpfr_extract_double requires at least two limbs */ if (sizer < MPFR_LIMBS_PER_DOUBLE) - EXP(r) = __mpfr_extract_double (MANT(r), d, 0); + MPFR_EXP(r) = __mpfr_extract_double (MPFR_MANT(r), d, 0); else - EXP(r) = __mpfr_extract_double (MANT(r) + sizer - MPFR_LIMBS_PER_DOUBLE, d, 1); + MPFR_EXP(r) = __mpfr_extract_double (MPFR_MANT(r) + sizer - MPFR_LIMBS_PER_DOUBLE, d, 1); if (sizer > MPFR_LIMBS_PER_DOUBLE) - MPN_ZERO(MANT(r), sizer - MPFR_LIMBS_PER_DOUBLE); + MPN_ZERO(MPFR_MANT(r), sizer - MPFR_LIMBS_PER_DOUBLE); - count_leading_zeros(cnt, MANT(r)[sizer-1]); - if (cnt) mpn_lshift(MANT(r), MANT(r), sizer, cnt); + count_leading_zeros(cnt, MPFR_MANT(r)[sizer-1]); + if (cnt) mpn_lshift(MPFR_MANT(r), MPFR_MANT(r), sizer, cnt); - EXP(r) -= cnt; - if (MPFR_SIGN(r)*signd<0) CHANGE_SIGN(r); + MPFR_EXP(r) -= cnt; + if (MPFR_SIGN(r)*signd<0) MPFR_CHANGE_SIGN(r); - mpfr_round(r, rnd_mode, PREC(r)); + mpfr_round(r, rnd_mode, MPFR_PREC(r)); return; } @@ -285,14 +285,14 @@ mpfr_get_d2(src, e) mp_ptr qp; int negative; - if (FLAG_NAN(src)) { + if (MPFR_IS_NAN(src)) { #ifdef DEBUG printf("recognized NaN\n"); #endif return NaN; } - if (NOTZERO(src)==0) return 0.0; - size = 1+(PREC(src)-1)/BITS_PER_MP_LIMB; - qp = MANT(src); + if (MPFR_NOTZERO(src)==0) return 0.0; + size = 1+(MPFR_PREC(src)-1)/BITS_PER_MP_LIMB; + qp = MPFR_MANT(src); negative = (MPFR_SIGN(src) < 0); /* Warning: don't compute the abs(res) and set the sign afterwards, @@ -324,6 +324,6 @@ mpfr_get_d(src) mpfr_srcptr src; #endif { - return mpfr_get_d2(src, EXP(src)); + return mpfr_get_d2(src, MPFR_EXP(src)); } @@ -38,16 +38,16 @@ mpfr_set_f(y, x, rnd_mode) mp_limb_t *my, *mx, *tmp; unsigned long cnt, sx, sy; TMP_DECL(marker); - sx = ABS(SIZ(x)); sy = ABSSIZE(y); - my = MANT(y); mx = MANT(x); + sx = ABS(SIZ(x)); sy = MPFR_ABSSIZE(y); + my = MPFR_MANT(y); mx = MPFR_MANT(x); if (sx==0) { /* x is zero */ - SET_ZERO(y); return; + MPFR_SET_ZERO(y); return; } count_leading_zeros(cnt, mx[sx - 1]); - if (SIZ(x)*MPFR_SIGN(y)<0) CHANGE_SIGN(y); + if (SIZ(x)*MPFR_SIGN(y)<0) MPFR_CHANGE_SIGN(y); if (sy < sx) { @@ -56,7 +56,7 @@ mpfr_set_f(y, x, rnd_mode) tmp = (mp_limb_t*) TMP_ALLOC(xprec); if (cnt) mpn_lshift(tmp, mx, sx, cnt); else MPN_COPY(tmp, mx, sx); - mpfr_round_raw(my, tmp, xprec, (SIZ(x)<0), PREC(y), rnd_mode); + mpfr_round_raw(my, tmp, xprec, (SIZ(x)<0), MPFR_PREC(y), rnd_mode); } else { @@ -66,7 +66,7 @@ mpfr_set_f(y, x, rnd_mode) /* no rounding necessary, since y has a larger mantissa */ } - EXP(y) = EXP(x) * BITS_PER_MP_LIMB - cnt; + MPFR_EXP(y) = MPFR_EXP(x) * BITS_PER_MP_LIMB - cnt; TMP_FREE(marker); } diff --git a/set_prc_raw.c b/set_prc_raw.c index 8ae4fd636..a0be55303 100644 --- a/set_prc_raw.c +++ b/set_prc_raw.c @@ -37,12 +37,12 @@ mpfr_set_prec_raw (x, p) fprintf(stderr, "*** cannot set precision to 0 bits\n"); exit(1); } - if (p > ABSSIZE(x) * BITS_PER_MP_LIMB) { + if (p > MPFR_ABSSIZE(x) * BITS_PER_MP_LIMB) { fprintf(stderr, "*** precision too large for allocated space\n"); exit(1); } - PREC(x) = p; + MPFR_PREC(x) = p; } diff --git a/set_prec.c b/set_prec.c index 98ef02b61..bc1fac979 100644 --- a/set_prec.c +++ b/set_prec.c @@ -41,13 +41,13 @@ mpfr_set_prec (x, p) xsize = (p - 1)/BITS_PER_MP_LIMB + 1; /* new limb size */ - if (xsize > ABSSIZE(x)) { + if (xsize > MPFR_ABSSIZE(x)) { x -> _mp_d = (mp_ptr) (*_mp_reallocate_func) - (x -> _mp_d, ABSSIZE(x)*BYTES_PER_MP_LIMB, xsize * BYTES_PER_MP_LIMB); - SIZE(x) = xsize; /* new number of allocated limbs */ + (x -> _mp_d, MPFR_ABSSIZE(x)*BYTES_PER_MP_LIMB, xsize * BYTES_PER_MP_LIMB); + MPFR_SIZE(x) = xsize; /* new number of allocated limbs */ } - PREC(x) = p; + MPFR_PREC(x) = p; } mp_prec_t @@ -58,5 +58,5 @@ mpfr_get_prec (x) mpfr_srcptr x; #endif { - return PREC(x); + return MPFR_PREC(x); } @@ -43,17 +43,17 @@ mpfr_set_q (f, q, rnd) num = mpq_numref(q); sign = mpz_cmp_ui(num, 0); if (sign==0) { - SET_ZERO(f); + MPFR_SET_ZERO(f); return; } den = mpq_denref(q); - prec = PREC(f); + prec = MPFR_PREC(f); mpfr_init2(n, mpz_sizeinbase(num, 2)); mpfr_set_z(n, num, GMP_RNDZ); /* result is exact */ mpfr_init2(d, mpz_sizeinbase(den, 2)); mpfr_set_z(d, den, GMP_RNDZ); /* result is exact */ - PREC(f) = prec; + MPFR_PREC(f) = prec; mpfr_div(f, n, d, rnd); mpfr_clear(n); mpfr_clear(d); } @@ -37,31 +37,31 @@ mpfr_set_si(x, i, rnd_mode) { unsigned long xn, cnt; mp_limb_t ai, *xp; - if (i==0) { SET_ZERO(x); return; } - xn = (PREC(x)-1)/BITS_PER_MP_LIMB; + if (i==0) { MPFR_SET_ZERO(x); return; } + xn = (MPFR_PREC(x)-1)/BITS_PER_MP_LIMB; ai = ABS(i); count_leading_zeros(cnt, ai); - xp = MANT(x); + xp = MPFR_MANT(x); xp[xn] = ai << cnt; /* don't forget to put zero in lower limbs */ MPN_ZERO(xp, xn); - EXP(x) = BITS_PER_MP_LIMB - cnt; + MPFR_EXP(x) = BITS_PER_MP_LIMB - cnt; - /* round if PREC(x) smaller than length of i */ - if (PREC(x) < BITS_PER_MP_LIMB-cnt) { - cnt = mpfr_round_raw(xp+xn, xp+xn, BITS_PER_MP_LIMB-cnt, (ai<0), PREC(x), + /* round if MPFR_PREC(x) smaller than length of i */ + if (MPFR_PREC(x) < BITS_PER_MP_LIMB-cnt) { + cnt = mpfr_round_raw(xp+xn, xp+xn, BITS_PER_MP_LIMB-cnt, (ai<0), MPFR_PREC(x), rnd_mode); if (cnt) { /* special case 1.000...000 */ - EXP(x)++; + MPFR_EXP(x)++; xp[xn] = ((mp_limb_t) 1) << (BITS_PER_MP_LIMB-1); } } /* warning: don't change the precision of x! */ - if (i*MPFR_SIGN(x) < 0) CHANGE_SIGN(x); + if (i*MPFR_SIGN(x) < 0) MPFR_CHANGE_SIGN(x); return; } @@ -78,29 +78,29 @@ mpfr_set_ui(x, i, rnd_mode) { unsigned int xn, cnt; mp_limb_t *xp; - if (i==0) { SET_ZERO(x); return; } - xn = (PREC(x)-1)/BITS_PER_MP_LIMB; + if (i==0) { MPFR_SET_ZERO(x); return; } + xn = (MPFR_PREC(x)-1)/BITS_PER_MP_LIMB; count_leading_zeros(cnt, (mp_limb_t) i); - xp = MANT(x); + xp = MPFR_MANT(x); xp[xn] = ((mp_limb_t) i) << cnt; /* don't forget to put zero in lower limbs */ MPN_ZERO(xp, xn); - EXP(x) = BITS_PER_MP_LIMB - cnt; + MPFR_EXP(x) = BITS_PER_MP_LIMB - cnt; - /* round if PREC(x) smaller than length of i */ - if (PREC(x) < BITS_PER_MP_LIMB-cnt) { - cnt = mpfr_round_raw(xp+xn, xp+xn, BITS_PER_MP_LIMB-cnt, 0, PREC(x), + /* round if MPFR_PREC(x) smaller than length of i */ + if (MPFR_PREC(x) < BITS_PER_MP_LIMB-cnt) { + cnt = mpfr_round_raw(xp+xn, xp+xn, BITS_PER_MP_LIMB-cnt, 0, MPFR_PREC(x), rnd_mode); if (cnt) { /* special case 1.000...000 */ - EXP(x)++; + MPFR_EXP(x)++; xp[xn] = ((mp_limb_t) 1) << (BITS_PER_MP_LIMB-1); } } /* warning: don't change the precision of x! */ - if (MPFR_SIGN(x) < 0) CHANGE_SIGN(x); + if (MPFR_SIGN(x) < 0) MPFR_CHANGE_SIGN(x); return; } @@ -102,7 +102,7 @@ mpfr_set_str(x, str, base, rnd_mode) /* the number is mantissa*base^expn */ - q = (PREC(x)/BITS_PER_MP_LIMB)*BITS_PER_MP_LIMB; + q = (MPFR_PREC(x)/BITS_PER_MP_LIMB)*BITS_PER_MP_LIMB; mpfr_init(y); mpfr_init(z); @@ -126,8 +126,8 @@ mpfr_set_str(x, str, base, rnd_mode) /* now y is an approximation of mantissa*base^expn with error at most 2^e*ulp(y) */ - } while (mpfr_can_round(y, q-e, GMP_RNDN, rnd_mode, PREC(x))==0 - && q<=2*PREC(x)); + } while (mpfr_can_round(y, q-e, GMP_RNDN, rnd_mode, MPFR_PREC(x))==0 + && q<=2*MPFR_PREC(x)); mpfr_set(x, y, rnd_mode); diff --git a/set_str_raw.c b/set_str_raw.c index 2629690f4..eb2d2ed1f 100644 --- a/set_str_raw.c +++ b/set_str_raw.c @@ -49,7 +49,7 @@ mpfr_set_str_raw(x, str) long expn = 0, e; char *endstr2; xp = x -> _mp_d; - xsize = 1 + (PREC(x)-1)/BITS_PER_MP_LIMB; + xsize = 1 + (MPFR_PREC(x)-1)/BITS_PER_MP_LIMB; alloc = (strlen(str)+1) * sizeof(char); str0 = str2 = (char *) (*_mp_allocate_func) (alloc); @@ -104,7 +104,7 @@ mpfr_set_str_raw(x, str) if (cnt) mpn_lshift(xp, xp, xsize, cnt); x -> _mp_exp = expn - cnt; - x -> _mp_size = xsize; if (negative) CHANGE_SIGN(x); + x -> _mp_size = xsize; if (negative) MPFR_CHANGE_SIGN(x); (*_mp_free_func) (str0, alloc); @@ -35,27 +35,27 @@ mpfr_set_z (f, z, rnd) mp_rnd_t rnd; #endif { - int fn, zn, k, dif, sign_z, sh; mp_limb_t *fp = MANT(f), *zp, cc, c2; + int fn, zn, k, dif, sign_z, sh; mp_limb_t *fp = MPFR_MANT(f), *zp, cc, c2; sign_z = mpz_cmp_ui(z,0); if (sign_z==0) { - SET_ZERO(f); + MPFR_SET_ZERO(f); return 0; } - fn = 1 + (PREC(f)-1)/BITS_PER_MP_LIMB; + fn = 1 + (MPFR_PREC(f)-1)/BITS_PER_MP_LIMB; zn = ABS(SIZ(z)); dif = zn-fn; zp = PTR(z); count_leading_zeros(k, zp[zn-1]); - EXP(f) = zn*BITS_PER_MP_LIMB-k; - if (MPFR_SIGN(f)*sign_z<0) CHANGE_SIGN(f); + MPFR_EXP(f) = zn*BITS_PER_MP_LIMB-k; + if (MPFR_SIGN(f)*sign_z<0) MPFR_CHANGE_SIGN(f); if (dif>=0) { /* number has to be truncated */ if (k) { mpn_lshift(fp, zp + dif, fn, k); if (dif) *fp += zp[dif-1] >> (BITS_PER_MP_LIMB-k); } else MPN_COPY(fp, zp + dif, fn); - sh = fn*BITS_PER_MP_LIMB-PREC(f); + sh = fn*BITS_PER_MP_LIMB-MPFR_PREC(f); cc = *fp & (((mp_limb_t)1 << sh) - 1); *fp = *fp & ~cc; if (rnd==GMP_RNDN) { @@ -66,22 +66,22 @@ mp_rnd_t rnd_mode; int logn; int tmp_factor; int tmpi; - if (FLAG_NAN(x)) { SET_NAN(sinus); SET_NAN(cosinus); return 1; } - if (!NOTZERO(x)) { + if (MPFR_IS_NAN(x)) { MPFR_SET_NAN(sinus); MPFR_SET_NAN(cosinus); return 1; } + if (!MPFR_NOTZERO(x)) { mpfr_set_ui(sinus, 0, GMP_RNDN); mpfr_set_ui(cosinus, 1, GMP_RNDN); return 0; } prec_x = (int) ceil(log - ((double) (PREC(x)) / (double) BITS_PER_MP_LIMB) + ((double) (MPFR_PREC(x)) / (double) BITS_PER_MP_LIMB) /log(2.0)); logn = (int) ceil(log ((double) prec_x) /log(2.0)); if (logn < 2) logn = 2; - ttt = EXP(x); - mpfr_init2(x_copy,PREC(x)); + ttt = MPFR_EXP(x); + mpfr_init2(x_copy,MPFR_PREC(x)); mpfr_set(x_copy,x,GMP_RNDD); mpz_init(square); /* on fait le shift pour que le nombre soit inferieur a 1 */ @@ -89,9 +89,9 @@ mp_rnd_t rnd_mode; { shift_x = ttt; mpfr_mul_2exp(x_copy,x,-ttt, GMP_RNDN); - ttt = EXP(x_copy); + ttt = MPFR_EXP(x_copy); } - realprec = PREC(sinus)+logn; + realprec = MPFR_PREC(sinus)+logn; while (!good){ Prec = realprec + 2*shift + 2 + shift_x + factor; k = (int) ceil(log @@ -180,22 +180,22 @@ mp_rnd_t rnd_mode; mpfr_mul(tmp_cos, tmp_cos, tmp_cos, GMP_RNDD); mpfr_mul_2exp(tmp_cos, tmp_cos, 1, GMP_RNDD); mpfr_set_ui(tmp, 1,GMP_RNDN); - tmpi = -EXP(tmp_cos); + tmpi = -MPFR_EXP(tmp_cos); mpfr_sub(tmp_cos, tmp_cos, tmp, GMP_RNDD); /* rep\'erer si le nombre de chiffres obtenu est suffisant pour avoir la bonne pr\'ecision. Le probl\`eme : comment faire ? la pr\'ecision s'obtient en comparant (Prec-factor) a la pr\'ecision obtenue r\'eellement, celle-ci - \'etant donn\'ee par Prec + EXP(tmp_cos). - il faut donc comparer EXP(tmp_cos) a factor */ - tmp_factor -= -EXP(tmp_cos) + tmpi; + \'etant donn\'ee par Prec + MPFR_EXP(tmp_cos). + il faut donc comparer MPFR_EXP(tmp_cos) a factor */ + tmp_factor -= -MPFR_EXP(tmp_cos) + tmpi; if (tmp_factor <= 0) { factor += -tmp_factor + 5; goto try_again; } } - if (mpfr_can_round(tmp_sin, realprec, GMP_RNDD, rnd_mode, PREC(sinus))){ + if (mpfr_can_round(tmp_sin, realprec, GMP_RNDD, rnd_mode, MPFR_PREC(sinus))){ mpfr_set(sinus,tmp_sin,rnd_mode); mpfr_set(cosinus,tmp_cos,rnd_mode); good = 1; @@ -49,18 +49,18 @@ mpfr_sqrt (r, u, rnd_mode) char can_round = 0; TMP_DECL (marker); TMP_DECL(marker0); - if (FLAG_NAN(u) || MPFR_SIGN(u) == -1) { SET_NAN(r); return 0; } + if (MPFR_IS_NAN(u) || MPFR_SIGN(u) == -1) { MPFR_SET_NAN(r); return 0; } - prec = PREC(r); + prec = MPFR_PREC(r); - if (!NOTZERO(u)) + if (!MPFR_NOTZERO(u)) { - EXP(r) = 0; - MPN_ZERO(MANT(r), ABSSIZE(r)); + MPFR_EXP(r) = 0; + MPN_ZERO(MPFR_MANT(r), MPFR_ABSSIZE(r)); return 1; } - up = MANT(u); + up = MPFR_MANT(u); #ifdef DEBUG printf("Entering square root : "); @@ -70,9 +70,9 @@ mpfr_sqrt (r, u, rnd_mode) /* Compare the mantissas */ - usize = (PREC(u) - 1)/BITS_PER_MP_LIMB + 1; - rsize = ((PREC(r) + 2 + (EXP(u) & 1))/BITS_PER_MP_LIMB + 1) << 1; - rrsize = (PREC(r) + 2 + (EXP(u) & 1))/BITS_PER_MP_LIMB + 1; + usize = (MPFR_PREC(u) - 1)/BITS_PER_MP_LIMB + 1; + rsize = ((MPFR_PREC(r) + 2 + (MPFR_EXP(u) & 1))/BITS_PER_MP_LIMB + 1) << 1; + rrsize = (MPFR_PREC(r) + 2 + (MPFR_EXP(u) & 1))/BITS_PER_MP_LIMB + 1; /* One extra bit is needed in order to get the square root with enough precision ; take one extra bit for rrsize in order to solve more easily the problem of rounding to nearest. @@ -82,9 +82,9 @@ mpfr_sqrt (r, u, rnd_mode) */ TMP_MARK(marker0); - if (EXP(u) & 1) /* Shift u one bit to the right */ + if (MPFR_EXP(u) & 1) /* Shift u one bit to the right */ { - if (PREC(u) & (BITS_PER_MP_LIMB - 1)) + if (MPFR_PREC(u) & (BITS_PER_MP_LIMB - 1)) { up = TMP_ALLOC(usize*BYTES_PER_MP_LIMB); mpn_rshift(up, u->_mp_d, usize, 1); @@ -92,14 +92,14 @@ mpfr_sqrt (r, u, rnd_mode) else { up = TMP_ALLOC((usize + 1)*BYTES_PER_MP_LIMB); - if (mpn_rshift(up + 1, u->_mp_d, ABSSIZE(u), 1)) + if (mpn_rshift(up + 1, u->_mp_d, MPFR_ABSSIZE(u), 1)) up [0] = ((mp_limb_t) 1) << (BITS_PER_MP_LIMB - 1); else up[0] = 0; usize++; } } - EXP(r) = ((EXP(u) + (EXP(u) & 1)) / 2) ; + MPFR_EXP(r) = ((MPFR_EXP(u) + (MPFR_EXP(u) & 1)) / 2) ; do { @@ -141,7 +141,7 @@ mpfr_sqrt (r, u, rnd_mode) #endif can_round = (mpfr_can_round_raw(rp, rrsize, 1, err, - GMP_RNDZ, rnd_mode, PREC(r))); + GMP_RNDZ, rnd_mode, MPFR_PREC(r))); /* If we used all the limbs of both the dividend and the divisor, then we have the correct RNDZ rounding */ @@ -173,8 +173,8 @@ mpfr_sqrt (r, u, rnd_mode) if (can_round) { - cc = mpfr_round_raw(rp, rp, err, 0, PREC(r), rnd_mode); - rrsize = (PREC(r) - 1)/BITS_PER_MP_LIMB + 1; + cc = mpfr_round_raw(rp, rp, err, 0, MPFR_PREC(r), rnd_mode); + rrsize = (MPFR_PREC(r) - 1)/BITS_PER_MP_LIMB + 1; } else /* Use the return value of sqrtrem to decide of the rounding */ @@ -190,7 +190,7 @@ mpfr_sqrt (r, u, rnd_mode) case GMP_RNDN : /* Not in the situation ...0 111111 */ - rw = (PREC(r) + 1) & (BITS_PER_MP_LIMB - 1); + rw = (MPFR_PREC(r) + 1) & (BITS_PER_MP_LIMB - 1); if (rw) { rw = BITS_PER_MP_LIMB - rw; nw = 0; } else nw = 1; if ((rp[nw] >> rw) & 1 && /* Not 0111111111 */ (q_limb || /* Nonzero remainder */ @@ -202,7 +202,7 @@ mpfr_sqrt (r, u, rnd_mode) case GMP_RNDU : if (q_limb) { - t = PREC(r) & (BITS_PER_MP_LIMB - 1); + t = MPFR_PREC(r) & (BITS_PER_MP_LIMB - 1); if (t) { cc = mpn_add_1(rp, rp, rrsize, 1 << (BITS_PER_MP_LIMB - t)); @@ -220,12 +220,12 @@ mpfr_sqrt (r, u, rnd_mode) fin: rsize = rrsize; - rrsize = (PREC(r) - 1)/BITS_PER_MP_LIMB + 1; + rrsize = (MPFR_PREC(r) - 1)/BITS_PER_MP_LIMB + 1; MPN_COPY(r->_mp_d, rp + rsize - rrsize, rrsize); - if (PREC(r) & (BITS_PER_MP_LIMB - 1)) - MANT(r) [0] &= ~(((mp_limb_t)1 << (BITS_PER_MP_LIMB - - (PREC(r) & (BITS_PER_MP_LIMB - 1)))) - 1) ; + if (MPFR_PREC(r) & (BITS_PER_MP_LIMB - 1)) + MPFR_MANT(r) [0] &= ~(((mp_limb_t)1 << (BITS_PER_MP_LIMB - + (MPFR_PREC(r) & (BITS_PER_MP_LIMB - 1)))) - 1) ; TMP_FREE(marker0); TMP_FREE (marker); return exact; @@ -49,12 +49,12 @@ mpfr_sqrt_ui (r, u, rnd_mode) MON_INIT(up, uu, BITS_PER_MP_LIMB, 1); count_leading_zeros(cnt, (mp_limb_t) u); *up = (mp_limb_t) u << cnt; - EXP(uu) = BITS_PER_MP_LIMB-cnt; + MPFR_EXP(uu) = BITS_PER_MP_LIMB-cnt; exact = mpfr_sqrt(r, uu, rnd_mode); TMP_FREE(marker); } - else SET_ZERO(r); + else MPFR_SET_ZERO(r); return exact; } @@ -81,56 +81,56 @@ mpfr_sub1(a, b, c, rnd_mode, diff_exp) printf("b=%1.20e c=%1.20e\n",mpfr_get_d(b),mpfr_get_d(c)); #endif cancel = mpfr_cmp2(b, c); - /* we have to take into account the first (PREC(a)+cancel) bits from b */ + /* we have to take into account the first (MPFR_PREC(a)+cancel) bits from b */ cancel1 = cancel/BITS_PER_MP_LIMB; cancel2 = cancel%BITS_PER_MP_LIMB; TMP_MARK(marker); - ap = MANT(a); - bp = MANT(b); - cp = MANT(c); + ap = MPFR_MANT(a); + bp = MPFR_MANT(b); + cp = MPFR_MANT(c); if (ap == bp) { - bp = (mp_ptr) TMP_ALLOC(ABSSIZE(b) * BYTES_PER_MP_LIMB); - MPN_COPY (bp, ap, ABSSIZE(b)); + bp = (mp_ptr) TMP_ALLOC(MPFR_ABSSIZE(b) * BYTES_PER_MP_LIMB); + MPN_COPY (bp, ap, MPFR_ABSSIZE(b)); if (ap == cp) { cp = bp; } } else if (ap == cp) { - cp = (mp_ptr) TMP_ALLOC (ABSSIZE(c) * BYTES_PER_MP_LIMB); - MPN_COPY(cp, ap, ABSSIZE(c)); + cp = (mp_ptr) TMP_ALLOC (MPFR_ABSSIZE(c) * BYTES_PER_MP_LIMB); + MPN_COPY(cp, ap, MPFR_ABSSIZE(c)); } - an = (PREC(a)-1)/BITS_PER_MP_LIMB+1; /* number of significant limbs of a */ - sh = an*BITS_PER_MP_LIMB-PREC(a); /* non-significant bits in low limb */ - bn = (PREC(b)-1)/BITS_PER_MP_LIMB+1; /* number of significant limbs of b */ - cn = (PREC(c)-1)/BITS_PER_MP_LIMB + 1; - EXP(a) = EXP(b)-cancel; + an = (MPFR_PREC(a)-1)/BITS_PER_MP_LIMB+1; /* number of significant limbs of a */ + sh = an*BITS_PER_MP_LIMB-MPFR_PREC(a); /* non-significant bits in low limb */ + bn = (MPFR_PREC(b)-1)/BITS_PER_MP_LIMB+1; /* number of significant limbs of b */ + cn = (MPFR_PREC(c)-1)/BITS_PER_MP_LIMB + 1; + MPFR_EXP(a) = MPFR_EXP(b)-cancel; /* adjust sign to that of b */ - if (MPFR_SIGN(a)*MPFR_SIGN(b)<0) CHANGE_SIGN(a); + if (MPFR_SIGN(a)*MPFR_SIGN(b)<0) MPFR_CHANGE_SIGN(a); /* case 1: diff_exp>=prec(a), i.e. c only affects the last bit through rounding */ - dif = PREC(a)-diff_exp; + dif = MPFR_PREC(a)-diff_exp; #ifdef DEBUG - printf("PREC(a)=%d an=%u PREC(b)=%d bn=%u PREC(c)=%d diff_exp=%u dif=%d cancel=%d\n", - PREC(a),an,PREC(b),bn,PREC(c),diff_exp,dif,cancel); + printf("MPFR_PREC(a)=%d an=%u MPFR_PREC(b)=%d bn=%u MPFR_PREC(c)=%d diff_exp=%u dif=%d cancel=%d\n", + MPFR_PREC(a),an,MPFR_PREC(b),bn,MPFR_PREC(c),diff_exp,dif,cancel); #endif - if (dif<=0) { /* diff_exp>=PREC(a): c does not overlap with a */ - /* either PREC(b)<=PREC(a), and we can copy the mantissa of b directly - into that of a, or PREC(b)>PREC(a) and we have to round b-c */ - if (PREC(b)<=PREC(a)+cancel) { + if (dif<=0) { /* diff_exp>=MPFR_PREC(a): c does not overlap with a */ + /* either MPFR_PREC(b)<=MPFR_PREC(a), and we can copy the mantissa of b directly + into that of a, or MPFR_PREC(b)>MPFR_PREC(a) and we have to round b-c */ + if (MPFR_PREC(b)<=MPFR_PREC(a)+cancel) { if (cancel2) mpn_lshift(ap+(an-bn+cancel1), bp, bn-cancel1, cancel2); else MPN_COPY(ap+(an-bn+cancel1), bp, bn-cancel1); /* fill low significant limbs with zero */ MPN_ZERO(ap, an-bn+cancel1); /* now take c into account */ if (rnd_mode==GMP_RNDN) { /* to nearest */ - /* if diff_exp > PREC(a), no change */ - if (diff_exp==PREC(a)) { + /* if diff_exp > MPFR_PREC(a), no change */ + if (diff_exp==MPFR_PREC(a)) { /* if c is not zero, then as it is normalized, we have to subtract one to the lsb of a if c>1/2, or c=1/2 and lsb(a)=1 (round to even) */ - if (NOTZERO(c)) { /* c is not zero */ + if (MPFR_NOTZERO(c)) { /* c is not zero */ /* check whether mant(c)=1/2 or not */ cc = *cp - (ONE<<(BITS_PER_MP_LIMB-1)); if (cc==0) { - bp = cp+(PREC(c)-1)/BITS_PER_MP_LIMB; + bp = cp+(MPFR_PREC(c)-1)/BITS_PER_MP_LIMB; while (cp<bp && cc==0) cc = *++cp; } if (cc || (ap[an-1] & ONE<<sh)) goto sub_one_ulp; @@ -139,23 +139,23 @@ mpfr_sub1(a, b, c, rnd_mode, diff_exp) } else if (ap[an-1]==0) { /* case b=2^n */ ap[an-1] = ONE << (BITS_PER_MP_LIMB-1); - EXP(a)++; + MPFR_EXP(a)++; } } - else if ((ISNONNEG(b) && rnd_mode==GMP_RNDU) || - (ISNEG(b) && rnd_mode==GMP_RNDD)) { + else if ((MPFR_ISNONNEG(b) && rnd_mode==GMP_RNDU) || + (MPFR_ISNEG(b) && rnd_mode==GMP_RNDD)) { /* round up: nothing to do */ if (ap[an-1]==0) { /* case b=2^n */ ap[an-1] = ONE << (BITS_PER_MP_LIMB-1); - EXP(a)++; + MPFR_EXP(a)++; } } else { /* round down: subtract 1 ulp iff c<>0 */ - if (NOTZERO(c)) goto sub_one_ulp; + if (MPFR_NOTZERO(c)) goto sub_one_ulp; } } - else { /* PREC(b)>PREC(a) : we have to round b-c */ + else { /* MPFR_PREC(b)>MPFR_PREC(a) : we have to round b-c */ k=bn-an; /* first copy the 'an' most significant limbs of b to a */ MPN_COPY(ap, bp+k, an); @@ -216,7 +216,7 @@ mpfr_sub1(a, b, c, rnd_mode, diff_exp) if (cout==0) cout=(cc!=0); if (rnd_mode==GMP_RNDU) sign=1; else if (rnd_mode==GMP_RNDD || rnd_mode==GMP_RNDZ) sign=-1; - if (ISNEG(b)) sign=-sign; + if (MPFR_ISNEG(b)) sign=-sign; /* even rounding rule: */ if (rnd_mode==GMP_RNDN && cout==0 && (*ap & (ONE<<sh))) cout=sign; /* round towards infinity if sign=1, towards zero otherwise */ @@ -225,7 +225,7 @@ mpfr_sub1(a, b, c, rnd_mode, diff_exp) } } } - else { /* case 2: diff_exp < PREC(a) : c overlaps with a by dif bits */ + else { /* case 2: diff_exp < MPFR_PREC(a) : c overlaps with a by dif bits */ /* first copy upper part of c into a (after shift) */ int overlap; dif += cancel; @@ -238,7 +238,7 @@ mpfr_sub1(a, b, c, rnd_mode, diff_exp) #ifdef DEBUG printf("cancel=%d dif=%d k=%d cn=%d sh=%d\n",cancel,dif,k,cn,sh); #endif - if (dif<=PREC(c)) { /* c has to be truncated */ + if (dif<=MPFR_PREC(c)) { /* c has to be truncated */ dif = dif % BITS_PER_MP_LIMB; dif = (dif) ? BITS_PER_MP_LIMB-dif-sh : -sh; /* we have to shift by dif bits to the right */ @@ -300,7 +300,7 @@ mpfr_sub1(a, b, c, rnd_mode, diff_exp) /* here overlap=1 iff ulp(c)<ulp(a) */ /* then put high limbs to zero */ /* now add 'an' upper limbs of b in place */ - if (PREC(b)<=PREC(a)+cancel) { int i, imax; + if (MPFR_PREC(b)<=MPFR_PREC(a)+cancel) { int i, imax; overlap += 2; /* invert low limbs */ imax = (int)an-(int)bn+cancel1; @@ -310,7 +310,7 @@ mpfr_sub1(a, b, c, rnd_mode, diff_exp) /* warning: mpn_sub_1 doesn't accept a zero length */ if (i<an) mpn_sub_1(ap+i, ap+i, an-i, ONE-cc); } - else /* PREC(b) > PREC(a): we have to truncate b */ + else /* MPFR_PREC(b) > MPFR_PREC(a): we have to truncate b */ mpn_sub_lshift_n(ap, bp+(bn-an-cancel1), an, cancel2, an); /* remains to do the rounding */ #ifdef DEBUG @@ -320,10 +320,10 @@ mpfr_sub1(a, b, c, rnd_mode, diff_exp) if (rnd_mode==GMP_RNDN) { /* to nearest */ int kc; /* four cases: overlap = - (0) PREC(b) > PREC(a) and diff_exp+PREC(c) <= PREC(a) - (1) PREC(b) > PREC(a) and diff_exp+PREC(c) > PREC(a) - (2) PREC(b) <= PREC(a) and diff_exp+PREC(c) <= PREC(a) - (3) PREC(b) <= PREC(a) and diff_exp+PREC(c) > PREC(a) */ + (0) MPFR_PREC(b) > MPFR_PREC(a) and diff_exp+MPFR_PREC(c) <= MPFR_PREC(a) + (1) MPFR_PREC(b) > MPFR_PREC(a) and diff_exp+MPFR_PREC(c) > MPFR_PREC(a) + (2) MPFR_PREC(b) <= MPFR_PREC(a) and diff_exp+MPFR_PREC(c) <= MPFR_PREC(a) + (3) MPFR_PREC(b) <= MPFR_PREC(a) and diff_exp+MPFR_PREC(c) > MPFR_PREC(a) */ switch (overlap) { case 1: /* both b and c to round */ @@ -373,8 +373,8 @@ mpfr_sub1(a, b, c, rnd_mode, diff_exp) /* otherwise the result is exact: nothing to do */ } } - else if ((ISNONNEG(b) && rnd_mode==GMP_RNDU) || - (ISNEG(b) && rnd_mode==GMP_RNDD)) { + else if ((MPFR_ISNONNEG(b) && rnd_mode==GMP_RNDU) || + (MPFR_ISNEG(b) && rnd_mode==GMP_RNDD)) { cc = *ap & ((ONE<<sh)-1); *ap &= ~cc; /* truncate last bits */ if (cc) goto add_one_ulp; /* will happen most of the time */ @@ -513,12 +513,12 @@ mpfr_sub(a, b, c, rnd_mode) { int diff_exp; - if (FLAG_NAN(b) || FLAG_NAN(c)) { SET_NAN(a); return; } + if (MPFR_IS_NAN(b) || MPFR_IS_NAN(c)) { MPFR_SET_NAN(a); return; } - if (!NOTZERO(b)) { mpfr_neg(a, c, rnd_mode); return; } - if (!NOTZERO(c)) { mpfr_set(a, b, rnd_mode); return; } + if (!MPFR_NOTZERO(b)) { mpfr_neg(a, c, rnd_mode); return; } + if (!MPFR_NOTZERO(c)) { mpfr_set(a, b, rnd_mode); return; } - diff_exp = EXP(b)-EXP(c); + diff_exp = MPFR_EXP(b)-MPFR_EXP(c); if (MPFR_SIGN(b) == MPFR_SIGN(c)) { /* signs are equal, it's a real subtraction */ if (diff_exp<0) { @@ -526,20 +526,20 @@ mpfr_sub(a, b, c, rnd_mode) if (rnd_mode==GMP_RNDU) rnd_mode=GMP_RNDD; else if (rnd_mode==GMP_RNDD) rnd_mode=GMP_RNDU; mpfr_sub1(a, c, b, rnd_mode, -diff_exp); - CHANGE_SIGN(a); + MPFR_CHANGE_SIGN(a); } else if (diff_exp>0) mpfr_sub1(a, b, c, rnd_mode, diff_exp); else { /* diff_exp=0 */ 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); + if (diff_exp==0) MPFR_SET_ZERO(a); else if (diff_exp*MPFR_SIGN(b)>0) mpfr_sub1(a, b, c, rnd_mode, 0); else { /* exchange rounding modes towards +/- infinity */ if (rnd_mode==GMP_RNDU) rnd_mode=GMP_RNDD; else if (rnd_mode==GMP_RNDD) rnd_mode=GMP_RNDU; mpfr_sub1(a, c, b, rnd_mode, 0); - CHANGE_SIGN(a); + MPFR_CHANGE_SIGN(a); } } } @@ -549,7 +549,7 @@ mpfr_sub(a, b, c, rnd_mode) if (rnd_mode==GMP_RNDU) rnd_mode=GMP_RNDD; else if (rnd_mode==GMP_RNDD) rnd_mode=GMP_RNDU; mpfr_add1(a, c, b, rnd_mode, -diff_exp); - CHANGE_SIGN(a); + MPFR_CHANGE_SIGN(a); } else mpfr_add1(a, b, c, rnd_mode, diff_exp); } @@ -47,7 +47,7 @@ mpfr_sub_ui(y, x, u, rnd_mode) MON_INIT(up, uu, BITS_PER_MP_LIMB, 1); count_leading_zeros(cnt, (mp_limb_t) u); *up = (mp_limb_t) u << cnt; - EXP(uu) = BITS_PER_MP_LIMB-cnt; + MPFR_EXP(uu) = BITS_PER_MP_LIMB-cnt; mpfr_sub(y, x, uu, rnd_mode); @@ -82,13 +82,13 @@ FUNC_NAME (r, u) mp_exp_t exp; int signu; long diff; - if (FLAG_NAN(u)) { SET_NAN(r); return; } - if (!NOTZERO(u)) { SET_ZERO(r); return; } + if (MPFR_IS_NAN(u)) { MPFR_SET_NAN(r); return; } + if (!MPFR_NOTZERO(u)) { MPFR_SET_ZERO(r); return; } - signu = SIZE(u); - rp = MANT(r); - exp = EXP(u); - prec = (PREC(r) - 1)/BITS_PER_MP_LIMB + 1; + signu = MPFR_SIZE(u); + rp = MPFR_MANT(r); + exp = MPFR_EXP(u); + prec = (MPFR_PREC(r) - 1)/BITS_PER_MP_LIMB + 1; /* Single out the case where |u| < 1. */ if (exp <= 0) @@ -99,21 +99,21 @@ FUNC_NAME (r, u) rp[prec-1] = (mp_limb_t) 1 << (BITS_PER_MP_LIMB-1); MPN_ZERO(rp, prec-1); /* sign of result is that of u */ - if (MPFR_SIGN(r) * signu < 0) CHANGE_SIGN(r); - EXP(r) = 1; + if (MPFR_SIGN(r) * signu < 0) MPFR_CHANGE_SIGN(r); + MPFR_EXP(r) = 1; return; } #endif - SET_ZERO(r); + MPFR_SET_ZERO(r); return; } - asize = (PREC(u) - 1)/BITS_PER_MP_LIMB + 1; + asize = (MPFR_PREC(u) - 1)/BITS_PER_MP_LIMB + 1; #ifdef _MPFR_FLOOR_OR_CEIL ignored_n = 0; #endif - up = MANT(u); + up = MPFR_MANT(u); if (asize > prec) { @@ -162,6 +162,6 @@ FUNC_NAME (r, u) if (rw) rp[0] &= ~((((mp_limb_t)1)<<rw) - (mp_limb_t)1); - EXP(r) = exp; - if (MPFR_SIGN(r) * signu < 0) CHANGE_SIGN(r); + MPFR_EXP(r) = exp; + if (MPFR_SIGN(r) * signu < 0) MPFR_CHANGE_SIGN(r); } @@ -47,11 +47,11 @@ mpfr_ui_div(y, u, x, rnd_mode) MON_INIT(up, uu, BITS_PER_MP_LIMB, 1); count_leading_zeros(cnt, (mp_limb_t) u); *up = (mp_limb_t) u << cnt; - EXP(uu) = BITS_PER_MP_LIMB-cnt; + MPFR_EXP(uu) = BITS_PER_MP_LIMB-cnt; mpfr_div(y, uu, x, rnd_mode); TMP_FREE(marker); } - else SET_ZERO(y); /* if u=0, then set y to 0 */ + else MPFR_SET_ZERO(y); /* if u=0, then set y to 0 */ } @@ -47,7 +47,7 @@ mpfr_ui_sub(y, u, x, rnd_mode) MON_INIT(up, uu, BITS_PER_MP_LIMB, 1); count_leading_zeros(cnt, (mp_limb_t) u); *up = (mp_limb_t) u << cnt; - EXP(uu) = BITS_PER_MP_LIMB-cnt; + MPFR_EXP(uu) = BITS_PER_MP_LIMB-cnt; mpfr_sub(y, uu, x, rnd_mode); diff --git a/urandomb.c b/urandomb.c index de78d61cc..cd986c1ab 100644 --- a/urandomb.c +++ b/urandomb.c @@ -42,8 +42,8 @@ mpfr_urandomb (rop, rstate) mp_exp_t exp; unsigned long cnt, nbits; - rp = MANT(rop); - nbits = PREC(rop); + rp = MPFR_MANT(rop); + nbits = MPFR_PREC(rop); nlimbs = (nbits + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB; _gmp_rand (rp, rstate, nbits); @@ -71,5 +71,5 @@ mpfr_urandomb (rop, rstate) /* cnt is the number of non significant bits in the low limb */ rp[0] &= ~((((mp_limb_t)1)<<cnt) - 1); - EXP (rop) = exp; + MPFR_EXP (rop) = exp; } @@ -41,7 +41,7 @@ mpfr_zeta(result, op, rnd_mode) /* first version */ if (mpfr_get_d(op) != 2.0 || rnd_mode != GMP_RNDN - || PREC(result) != 53) { + || MPFR_PREC(result) != 53) { fprintf(stderr, "not yet implemented\n"); exit(1); } |