Initial revision

git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@2 280ebfd0-de03-0410-8827-d642c229c3f4

71 files changed, 4076 insertions, 0 deletions

diff --git a/.pure b/.pure
new file mode 100644
index 000000000..e69de29bb
--- /dev/null
+++ b/.pure
diff --git a/BUGS b/BUGS
new file mode 100644
index 000000000..aae547816
--- /dev/null
+++ b/BUGS
@@ -0,0 +1,2 @@
+Probablement sans doute surement plein.
+Que vaut un arrondi a 0 bits de precision ?...
\ No newline at end of file
diff --git a/Makefile.exp b/Makefile.exp
new file mode 100644
index 000000000..ef2f594df
--- /dev/null
+++ b/Makefile.exp
@@ -0,0 +1,22 @@
+CC=gcc
+INCFLAGS=-I../../../../hanrot/gmp-2.0.2 -I../../../../hanrot/gmp-2.0.2/mpn
+LDFLAGS=../../../../hanrot/gmp-2.0.2/libgmp.a 
+CFLAGS=-g -DExp
+
+essai: libmpfr.a essai.c
+	$(CC) $(INCFLAGS) $(CFLAGS) -o essai essai.c libmpfr.a $(LDFLAGS)
+
+libmpfr.a: init.o set_d.o round.o
+	$(AR) cr libmpfr.a init.o set_d.o round.o
+
+clean: 
+	rm essai *.o *~ *.a
+
+init.o: init.c mpfr.h
+	$(CC) $(CFLAGS) $(INCFLAGS) -c init.c
+
+set_d.o: set_d.c mpfr.h
+	$(CC) $(CFLAGS) $(INCFLAGS) -c set_d.c
+
+round.o: round.c mpfr.h
+	$(CC) $(CFLAGS) $(INCFLAGS) -c round.c
diff --git a/Makefile.msb b/Makefile.msb
new file mode 100644
index 000000000..835212b2a
--- /dev/null
+++ b/Makefile.msb
@@ -0,0 +1,154 @@
+CC=gcc
+INCFLAGS=-I$(GMPDIR)/include
+LDFLAGS=$(GMPDIR)/lib/libgmp.a -lm
+CFLAGS=-g -Wall -DMsb $(special) -D$(dirname)
+
+MPFR_OBJECTS = o.$(dirname)/add.o o.$(dirname)/clear.o o.$(dirname)/cmp.o o.$(dirname)/init.o o.$(dirname)/mul.o o.$(dirname)/print_raw.o o.$(dirname)/rnd_mode.o o.$(dirname)/round.o o.$(dirname)/set_d.o o.$(dirname)/set_prec.o o.$(dirname)/set_dfl_prec.o o.$(dirname)/set_dfl_rnd.o o.$(dirname)/neg.o o.$(dirname)/set_f.o o.$(dirname)/set_si.o o.$(dirname)/set_str_raw.o o.$(dirname)/mul_ui.o o.$(dirname)/mul_2exp.o o.$(dirname)/div_2exp.o o.$(dirname)/sub.o o.$(dirname)/set.o o.$(dirname)/div.o o.$(dirname)/sqrt.o o.$(dirname)/agm.o o.$(dirname)/get_str.o o.$(dirname)/cmp_ui.o
+
+dft target::
+	@echo "Utiliser soit make all, soit ./mmpfr"
+	@echo "./mmpfr peut prendre les options tests ou clean"
+
+libmpfr.a: mpfr.h $(MPFR_OBJECTS)
+	$(AR) cr o.$(dirname)/libmpfr.a $(MPFR_OBJECTS)
+	-chmod g+w o.$(dirname)/*.o o.$(dirname)/libmpfr.a
+	-ranlib o.$(dirname)/libmpfr.a
+
+clean: 
+	-rm o.$(dirname)/*.o *~ o.$(dirname)/*.a core
+
+dist:
+	-rm mpfr.tgz
+	tar cvf mpfr.tar Makefile *.c *.h tests/*.c tests/Makefile
+	gzip -9 mpfr.tar
+
+o.$(dirname)/add.o: mpfr.h add.c
+	$(CC) $(CFLAGS) $(INCFLAGS) -c add.c -o o.$(dirname)/add.o
+
+o.$(dirname)/sub.o: mpfr.h sub.c
+	$(CC) $(CFLAGS) $(INCFLAGS) -c sub.c -o o.$(dirname)/sub.o
+
+o.$(dirname)/set.o: mpfr.h set.c
+	$(CC) $(CFLAGS) $(INCFLAGS) -c set.c -o o.$(dirname)/set.o
+
+o.$(dirname)/neg.o: mpfr.h neg.c
+	$(CC) $(CFLAGS) $(INCFLAGS) -c neg.c -o o.$(dirname)/neg.o
+
+o.$(dirname)/clear.o: mpfr.h clear.c
+	$(CC) $(CFLAGS) $(INCFLAGS) -c clear.c -o o.$(dirname)/clear.o
+
+o.$(dirname)/cmp.o: mpfr.h cmp.c
+	-$(CC) $(CFLAGS) $(INCFLAGS) -c cmp.c -o o.$(dirname)/cmp.o
+
+o.$(dirname)/cmp_ui.o: mpfr.h cmp_ui.c
+	-$(CC) $(CFLAGS) $(INCFLAGS) -c cmp_ui.c -o o.$(dirname)/cmp_ui.o
+
+o.$(dirname)/init.o: mpfr.h init.c
+	$(CC) $(CFLAGS) $(INCFLAGS) -c init.c -o o.$(dirname)/init.o
+
+o.$(dirname)/mul.o: mpfr.h mul.c
+	$(CC) $(CFLAGS) $(INCFLAGS) -c mul.c -o o.$(dirname)/mul.o
+
+o.$(dirname)/print_raw.o: mpfr.h print_raw.c
+	$(CC) $(CFLAGS) $(INCFLAGS) -c print_raw.c -o o.$(dirname)/print_raw.o
+
+o.$(dirname)/rnd_mode.o: mpfr.h rnd_mode.c
+	$(CC) $(CFLAGS) $(INCFLAGS) -c rnd_mode.c -o o.$(dirname)/rnd_mode.o
+
+o.$(dirname)/set_f.o: mpfr.h set_f.c
+	$(CC) $(CFLAGS) $(INCFLAGS) -c set_f.c -o o.$(dirname)/set_f.o
+
+o.$(dirname)/set_si.o: mpfr.h set_si.c
+	$(CC) $(CFLAGS) $(INCFLAGS) -c set_si.c -o o.$(dirname)/set_si.o
+
+o.$(dirname)/round.o: mpfr.h round.c
+	-$(CC) $(CFLAGS) $(INCFLAGS) -c round.c -o o.$(dirname)/round.o
+
+o.$(dirname)/set_d.o: mpfr.h set_d.c
+	$(CC) $(CFLAGS) $(INCFLAGS) -c set_d.c -o o.$(dirname)/set_d.o
+
+o.$(dirname)/set_prec.o: mpfr.h set_prec.c
+	$(CC) $(CFLAGS) $(INCFLAGS) -c set_prec.c -o o.$(dirname)/set_prec.o
+
+o.$(dirname)/set_str_raw.o: mpfr.h set_str_raw.c
+	$(CC) $(CFLAGS) $(INCFLAGS) -c set_str_raw.c -o o.$(dirname)/set_str_raw.o
+
+o.$(dirname)/get_str.o: mpfr.h get_str.c
+	$(CC) $(CFLAGS) $(INCFLAGS) -c get_str.c -o o.$(dirname)/get_str.o
+
+o.$(dirname)/mul_2exp.o: mpfr.h mul_2exp.c
+	$(CC) $(CFLAGS) $(INCFLAGS) -c mul_2exp.c -o o.$(dirname)/mul_2exp.o
+
+o.$(dirname)/div_2exp.o: mpfr.h div_2exp.c
+	$(CC) $(CFLAGS) $(INCFLAGS) -c div_2exp.c -o o.$(dirname)/div_2exp.o
+
+o.$(dirname)/div.o: mpfr.h div.c
+	$(CC) $(CFLAGS) $(INCFLAGS) -c div.c -o o.$(dirname)/div.o
+
+o.$(dirname)/sqrt.o: mpfr.h sqrt.c
+	$(CC) $(CFLAGS) $(INCFLAGS) -c sqrt.c -o o.$(dirname)/sqrt.o
+
+o.$(dirname)/set_dfl_prec.o: mpfr.h set_dfl_prec.c
+	$(CC) $(CFLAGS) $(INCFLAGS) -c set_dfl_prec.c -o o.$(dirname)/set_dfl_prec.o
+
+o.$(dirname)/set_dfl_rnd.o: mpfr.h set_dfl_rnd.c
+	$(CC) $(CFLAGS) $(INCFLAGS) -c set_dfl_rnd.c -o o.$(dirname)/set_dfl_rnd.o
+
+o.$(dirname)/mul_ui.o: mpfr.h mul_ui.c
+	$(CC) $(CFLAGS) $(INCFLAGS) -c mul_ui.c -o o.$(dirname)/mul_ui.o
+
+o.$(dirname)/agm.o: mpfr.h agm.c o.$(dirname)/div.o
+	$(CC) $(CFLAGS) $(INCFLAGS) -c agm.c -o o.$(dirname)/agm.o
+
+tadd: tests/tadd.c o.$(dirname)/add.o libmpfr.a
+	$(CC) $(CFLAGS) $(INCFLAGS) -I. -o tests/tadd tests/tadd.c $(ODIR)/libmpfr.a $(LDFLAGS)
+
+tdiv: tests/tdiv.c o.$(dirname)/div.o libmpfr.a
+	$(CC) $(CFLAGS) $(INCFLAGS) -I. -o tests/tdiv tests/tdiv.c $(ODIR)/libmpfr.a $(LDFLAGS)
+
+tsqrt: tests/tsqrt.c o.$(dirname)/sqrt.o libmpfr.a
+	$(CC) $(CFLAGS) $(INCFLAGS) -I. -o tests/tsqrt tests/tsqrt.c $(ODIR)/libmpfr.a $(LDFLAGS)
+
+tcmp: tests/tcmp.c o.$(dirname)/add.o libmpfr.a
+	$(CC) $(CFLAGS) $(INCFLAGS) -I. -o tests/tcmp tests/tcmp.c $(ODIR)/libmpfr.a $(LDFLAGS)
+
+tcmp_ui: tests/tcmp_ui.c libmpfr.a
+	$(CC) $(CFLAGS) $(INCFLAGS) -I. -o tests/tcmp_ui tests/tcmp_ui.c $(ODIR)/libmpfr.a $(LDFLAGS)
+
+tcmp2: tests/tcmp2.c o.$(dirname)/add.o libmpfr.a
+	$(CC) $(CFLAGS) $(INCFLAGS) -I. -o tests/tcmp2 tests/tcmp2.c $(ODIR)/libmpfr.a $(LDFLAGS)
+
+tmul: tests/tmul.c o.$(dirname)/mul.o mul.c libmpfr.a
+	$(CC) $(CFLAGS) $(INCFLAGS) -I. -o tests/tmul tests/tmul.c $(ODIR)/libmpfr.a $(LDFLAGS)
+
+tround: tests/tround.c round.c $(ODIR)/libmpfr.a
+	$(CC) $(CFLAGS) $(INCFLAGS) -I. -o tests/tround tests/tround.c $(ODIR)/libmpfr.a $(LDFLAGS)
+
+tset_f: tests/tset_f.c set_f.c $(ODIR)/libmpfr.a
+	$(CC) $(CFLAGS) $(INCFLAGS) -I. -o tests/tset_f tests/tset_f.c $(ODIR)/libmpfr.a $(LDFLAGS)
+
+tset_d: tests/tset_d.c set_d.c libmpfr.a
+	$(CC) $(CFLAGS) $(INCFLAGS) -I. -o tests/tset_d tests/tset_d.c $(ODIR)/libmpfr.a $(LDFLAGS)
+
+tset_si: tests/tset_si.c set_si.c $(ODIR)/libmpfr.a
+	$(CC) $(CFLAGS) $(INCFLAGS) -I. -o tests/tset_si tests/tset_si.c $(ODIR)/libmpfr.a $(LDFLAGS)
+
+tmul_ui: tests/tmul_ui.c mul_ui.c $(ODIR)/libmpfr.a
+	$(CC) $(CFLAGS) $(INCFLAGS) -I. -o tests/tmul_ui tests/tmul_ui.c $(ODIR)/libmpfr.a $(LDFLAGS)
+
+tset_str: tests/tset_str.c set_str_raw.c $(ODIR)/libmpfr.a
+	$(CC) $(CFLAGS) $(INCFLAGS) -I. -o tests/tset_str tests/tset_str.c $(ODIR)/libmpfr.a $(LDFLAGS)
+
+tget_str: tests/tget_str.c o.$(dirname)/get_str.o libmpfr.a
+	$(CC) $(CFLAGS) $(INCFLAGS) -I. -o tests/tget_str tests/tget_str.c $(ODIR)/libmpfr.a $(LDFLAGS)
+
+tmul_2exp: tests/tmul_2exp.c mul_2exp.c $(ODIR)/libmpfr.a
+	$(CC) $(CFLAGS) $(INCFLAGS) -I. -o tests/tmul_2exp tests/tmul_2exp.c $(ODIR)/libmpfr.a $(LDFLAGS)
+
+tagm: tests/tagm.c agm.c libmpfr.a
+	$(CC) $(CFLAGS) $(INCFLAGS) -I. -o tests/tagm tests/tagm.c $(ODIR)/libmpfr.a $(LDFLAGS)
+
+
+all: 
+	./mmpfr; 
+
diff --git a/TODO b/TODO
new file mode 100644
index 000000000..e103b7e8c
--- /dev/null
+++ b/TODO
@@ -0,0 +1,2 @@
+- implanter tous les flags, qui sont pour le moment superbement ignores.
+J'ai quand meme laisse deux macros SIZ et SIZE selon qu'on affecte ou non.
diff --git a/add.c b/add.c
new file mode 100644
index 000000000..d0b4d2a3b
--- /dev/null
+++ b/add.c
@@ -0,0 +1,344 @@
+/* (c) PolKA project, Inria Lorraine */
+
+/* written by Paul Zimmermann, February 1999 */
+
+/* to do: - destination may be identical to operands
+*/
+
+
+#include <stdio.h>
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "mpfr.h"
+
+extern mpfr_sub1(mpfr_ptr, mpfr_srcptr, mpfr_srcptr, unsigned char, int);
+
+/* #define DEBUG2 */
+
+#ifdef DEBUG2
+mp_limb_t *ap0, *ap1;
+#define ck(x, dx, i) /* mp_limb_t *x; int dx; */ \
+{ \
+  if (x<ap0 || x+dx-1>ap1) { \
+    printf("error %d\n",i); \
+    printf("%1.20e,%d,%1.20e,%d,%d,%d\n",mpfr_get_d(b),PREC(b),mpfr_get_d(c), \
+	   PREC(c),PREC(a),rnd_mode); \
+    exit(1); } \
+}
+#else
+#define ck(x, dx, i) {}
+#endif
+
+/* signs of b and c are supposed equal,
+   diff_exp is the difference between the exponents of b and c,
+   which is supposed >= 0 */
+void mpfr_add1(a, b, c, rnd_mode, diff_exp) 
+mpfr_ptr a; mpfr_srcptr b, c; unsigned char rnd_mode; int diff_exp;
+{
+  mp_limb_t *ap, *bp, *cp, cc, c2, c3=0; unsigned int an,bn,cn; int sh,dif,k;
+
+#ifdef DEBUG2
+  printf("b=  "); mpfr_print_raw(b); putchar('\n');
+  printf(" =  %1.20e\n",mpfr_get_d(b));
+  printf("c=  "); for (k=0;k<diff_exp;k++) putchar(' '); mpfr_print_raw(c);
+  putchar('\n');
+#endif
+  ap = MANT(a);
+  bp = MANT(b);
+  an = (PREC(a)-1)/mp_bits_per_limb+1; /* number of significant limbs of a */
+#ifdef DEBUG2
+ap0=ap; ap1=ap+an-1;
+#endif
+  sh = an*mp_bits_per_limb-PREC(a); /* non-significant bits in low limb */
+  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);
+  /* case 1: diff_exp>=prec(a), i.e. c only affects the last bit
+     through rounding */
+  dif = PREC(a)-diff_exp;
+#ifdef DEBUG2
+  if (bn & ((mp_limb_t)1<<(mp_bits_per_limb-1))==0) {
+    printf("Error in mpfr_add: b not msb-normalized\n"); exit(1);
+  }
+printf("an=%u prec(a)=%d bn=%u prec(b)=%d prec(c)=%d diff_exp=%u dif=%d\n",
+       an,PREC(a),bn,PREC(b),PREC(c),diff_exp,dif);
+#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)) {
+ck(ap+an-bn,bn,1);
+      MPN_COPY(ap+(an-bn), bp, bn);
+      /* fill low significant limbs with zero */
+      for (bp=ap;bn<an;bn++) { ck(bp,1,2); *bp++=0; }
+      /* 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 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) */
+	  cp=MANT(c);
+	  if (NOTZERO(c)) { /* c is not zero */
+	    /* check whether mant(c)=1/2 or not */
+	    cc = *cp - (1<<(mp_bits_per_limb-1));
+	    if (cc==0) {
+	      bp = cp+(PREC(c)-1)/mp_bits_per_limb;
+	      while (cp<bp && cc==0) cc = *++cp;
+	    }
+ck(ap+an-1, 1, 3);
+	    if (cc || (ap[an-1] & 1<<sh)) goto add_one_ulp;
+	      /* mant(c) != 1/2 or mant(c) = 1/2: add 1 iff lsb(a)=1 */
+	  }
+	}
+      }
+      else if ((ISNONNEG(b) && rnd_mode==GMP_RNDU) || 
+	       (ISNEG(b) && rnd_mode==GMP_RNDD)) {
+	/* round up */
+	if (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 */
+      k=bn-an;
+      /* first copy the 'an' most significant limbs of b to a */
+ck(ap, an, 4);
+      MPN_COPY(ap, bp+k, an);
+      if (rnd_mode==GMP_RNDN) { /* to nearest */
+	/* first check whether the truncated bits from b are 1/2*lsb(a) */
+	if (sh) {
+ck(ap, 1, 5);
+	  cc = *ap & ((1<<sh)-1);
+	  *ap &= ~cc; /* truncate last bits */
+	  cc -= 1<<(sh-1);
+	}
+	else /* no bit to truncate */
+	  cc = bp[--k] - (1<<(mp_bits_per_limb-1));
+	if ((long)cc>0) goto add_one_ulp; /* trunc(b)>1/2*lsb(a) -> round up */
+	else if (cc==0) {
+	  while (k>1 && cc==0) cc=bp[--k];
+	  /* now if the truncated part of b = 1/2*lsb(a), check whether c=0 */
+ck(ap, 1, 6);
+	  if (NOTZERO(c) || (*ap & (1<<sh))) goto add_one_ulp;
+	  /* if trunc(b)+c is exactly 1/2*lsb(a) : round to even lsb */
+	}
+	/* if cc<0 : trunc(b) < 1/2*lsb(a) -> round down, i.e. do nothing */
+      }
+      else if ((ISNONNEG(b) && rnd_mode==GMP_RNDU) || 
+	       (ISNEG(b) && rnd_mode==GMP_RNDD)) {
+	/* first check whether trunc(b)+c is zero or not */
+	if (sh) {
+ck(ap, 1, 7);
+	  cc = *ap & ((1<<sh)-1); *ap &= ~cc; /* truncate last bits */
+	}
+	else cc = bp[--k] - (1<<(mp_bits_per_limb-1));
+	while (cc==0 && k>1) cc=bp[--k];
+	if (cc || NOTZERO(c)) goto add_one_ulp;
+      }
+      /* in the other cases (round to zero, or up/down with sign -/+),
+         nothing to do, since b and c don't overlap, there can't be any
+	 carry */
+    }
+  }
+  else { /* diff_exp < 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)/mp_bits_per_limb + 1; /* only the highest k limbs from c
+					 have to be considered */
+    cp = MANT(c);
+    cn = (PREC(c)-1)/mp_bits_per_limb + 1;
+ck(ap+k, an-k, 8);
+    MPN_ZERO(ap+k, an-k); /* do it now otherwise ap[k] may be destroyed
+			     in case dif<0 */
+#ifdef DEBUG2
+    printf("k=%d cn=%d\n",k,cn);
+#endif
+    if (dif<=PREC(c)) { /* c has to be truncated */
+      dif = dif % mp_bits_per_limb;
+      dif = (dif) ? mp_bits_per_limb-dif-sh : -sh;
+      /* we have to shift by dif bits to the right */
+      if (dif>0) { ck(ap, k, 9); mpn_rshift(ap, cp+(cn-k), k, dif); }
+      else if (dif<0) {
+ck(ap, k+1, 10);
+	ap[k] = mpn_lshift(ap, cp+(cn-k), k, -dif);
+	/* put the non-significant bits in low limb for further rounding */
+	ap[0] += cp[cn-k-1]>>(mp_bits_per_limb+dif);
+      }
+      else { ck(ap, k, 11); MPN_COPY(ap, cp+(cn-k), k); }
+      overlap=1;
+    }
+    else { /* c is not truncated, but we have to fill low limbs with 0 */
+      k = diff_exp/mp_bits_per_limb;
+      overlap = diff_exp%mp_bits_per_limb;
+      /* warning: a shift of zero bit is not allowed */
+      if (overlap) { 
+	cc=mpn_rshift(ap+(an-k-cn), cp, cn, overlap); 
+	if (an-k-cn>0) ap[an-k-cn-1]=cc;
+      }
+      else { ck(ap+(an-k-cn), cn, 13); MPN_COPY(ap+(an-k-cn), cp, cn); }
+      overlap=0;
+    }
+    /* 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)) {
+      overlap += 2;
+ck(ap+(an-bn), bn, 14);
+      cc = mpn_add_n(ap+(an-bn), ap+(an-bn), bp, bn);
+    }
+    else /* PREC(b) > PREC(a): we have to truncate b */
+{ ck(ap, an, 15);      cc = mpn_add_n(ap, ap, bp+(bn-an), an); }
+    if (cc) { /* shift one bit to the right */
+ck(ap, an, 16);
+      c3 = (ap[0]&1) && (PREC(a)%mp_bits_per_limb==0);
+      mpn_rshift(ap, ap, an, 1);
+      ap[an-1] += (mp_limb_t)1<<(mp_bits_per_limb-1);
+      EXP(a)++;
+    }
+    /* remains to do the rounding */
+    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) */
+#ifdef DEBUG2
+printf("overlap=%d\n",overlap);
+#endif
+      switch (overlap)
+	{
+        case 1: /* both b and c to round */
+	  kc = cn-k; /* remains kc limbs from c */
+	  k = bn-an; /* remains k limbs from b */
+	  /* truncate last bits and store the difference with 1/2*ulp in cc */
+ck(ap, 1, 17);
+	  cc = *ap & ((1<<sh)-1);
+	  *ap &= ~cc; /* truncate last bits */
+	  cc -= 1<<(sh-1);
+	  while ((cc==0 || cc==-1) && k!=0 && kc!=0) {
+	    kc--;
+	    cc += mpn_add_1(&c2, bp+(--k), 1,(cp[kc+1]<<(mp_bits_per_limb-dif))
+			    +(cp[kc]>>dif));
+	    if (cc==0 || cc==-1) cc=c2;
+	  }
+	  if ((long)cc>0) goto add_one_ulp;
+	  else if ((long)cc<-1) return; /* the carry can be at most 1 */
+	  else if (kc==0) goto round_b;
+	  /* else round c: go through */
+	case 3: /* only c to round */
+	  bp=cp; k=cn-k; goto to_nearest;
+	case 0: /* only b to round */
+        round_b:
+	  k=bn-an; dif=0; goto to_nearest;
+        /* otherwise the result is exact: nothing to do */
+	}
+    }
+    else if ((ISNONNEG(b) && rnd_mode==GMP_RNDU) || 
+	     (ISNEG(b) && rnd_mode==GMP_RNDD)) {
+ck(ap, 1, 18);
+      cc = *ap & ((1<<sh)-1);
+      *ap &= ~cc; /* truncate last bits */
+      if (cc || c3) goto add_one_ulp; /* will happen most of the time */
+      else { /* same four cases too */
+	int kc = cn-k; /* remains kc limbs from c */
+	switch (overlap)
+	{
+        case 1: /* both b and c to round */
+	  k = bn-an; /* remains k limbs from b */
+	  while (cc==0 && k!=0 && kc!=0) {
+	    kc--;
+	    cc = mpn_add_1(&c2, bp+(--k), 1,(cp[kc+1]<<(mp_bits_per_limb-dif))
+			   + (cp[kc]>>dif));
+	  }
+	  if (cc) goto add_one_ulp;
+	  else if (kc==0) goto round_b2;
+	  /* else round c: go through */
+	case 3: /* only c to round */
+	  while (kc) if (cp[--kc]) goto add_one_ulp;
+	  /* if dif>0 : remains to check last dif bits from c */
+	  if (dif>0 && (cp[0]<<(mp_bits_per_limb-dif))) goto add_one_ulp;
+	  break;
+	case 0: /* only b to round */
+        round_b2:
+	  k=bn-an;
+	  while (k) if (bp[--k]) goto add_one_ulp;
+        /* otherwise the result is exact: nothing to do */
+	}
+      }
+    }
+    /* else nothing to do: round towards zero, i.e. truncate last sh bits */
+    else {
+ck(ap, 1, 19);
+      *ap &= ~((1<<sh)-1);
+    }
+  }
+  goto end_of_add;
+    
+  to_nearest: /* 0 <= sh < mp_bits_per_limb : number of bits of a to truncate
+                 bp[k] : last significant limb from b */
+        if (sh) {
+ck(ap, 1, 20);
+	  cc = *ap & ((1<<sh)-1);
+	  *ap &= ~cc; /* truncate last bits */
+	  c2 = 1<<(sh-1);
+	}
+	else /* no bit to truncate */
+	  { cc = bp[--k]; c2 = 1<<(mp_bits_per_limb-1); }
+#ifdef DEBUG2
+printf("cc=%lu c2=%lu k=%u\n",cc,c2,k);
+#endif
+	if (cc>c2) goto add_one_ulp; /* trunc(b)>1/2*lsb(a) -> round up */
+	else if (cc==c2) {
+	  cc=0; while (k && cc==0) cc=bp[--k];
+#ifdef DEBUG2
+printf("cc=%lu\n",cc);
+#endif
+	  /* special case of rouding c shifted to the right */
+	  if (cc==0 && dif>0) cc=cp[0]<<(mp_bits_per_limb-dif);
+	  /* now if the truncated part of b = 1/2*lsb(a), check whether c=0 */
+ck(ap, 1, 21);
+	  if (cc || (*ap & (1<<sh))) goto add_one_ulp;
+	}
+        goto end_of_add;
+
+  add_one_ulp: /* add one unit in last place to a */
+ck(ap, an, 22);
+    cc = mpn_add_1(ap, ap, an, 1<<sh);
+    if (cc) { printf("carry(3) in mpfr_add\n"); exit(1); }
+
+ end_of_add:
+#ifdef DEBUG2
+printf("b+c="); mpfr_print_raw(a); putchar('\n');
+#endif
+  return;
+}
+
+void mpfr_add(a, b, c, rnd_mode) 
+mpfr_ptr a; mpfr_srcptr b, c; unsigned char rnd_mode;
+{
+  int diff_exp;
+
+  diff_exp = EXP(b)-EXP(c);
+  if (SIGN(b) != SIGN(c)) { /* signs differ, it's a subtraction */
+    if (diff_exp<0) {
+      mpfr_sub1(a, c, b, rnd_mode, -diff_exp);
+    }
+    else if (diff_exp>0) mpfr_sub1(a, b, c, rnd_mode, diff_exp);
+    else { /* diff_exp=0 */
+      CHANGE_SIGN(c);
+      diff_exp = mpfr_cmp(b,c);
+      CHANGE_SIGN(c);
+      /* if b>0 and diff_exp>0 or b<0 and diff_exp<0: abs(b) > abs(c) */
+      if (diff_exp*SIGN(b)>=0) mpfr_sub1(a, b, c, rnd_mode, 0);
+      else {
+	mpfr_sub1(a, c, b, rnd_mode, 0);
+      }
+    }
+  }
+  else /* signs are equal, it's an addition */
+    if (diff_exp<0) mpfr_add1(a, c, b, rnd_mode, -diff_exp);
+    else mpfr_add1(a, b, c, rnd_mode, diff_exp);
+}
+
diff --git a/agm.c b/agm.c
new file mode 100644
index 000000000..3677db915
--- /dev/null
+++ b/agm.c
@@ -0,0 +1,91 @@
+#include <stdio.h>
+#include <math.h>
+#include "gmp.h"
+#include "mpfr.h"
+
+
+void mpfr_agm(mpfr_ptr r, mpfr_srcptr a, mpfr_srcptr b, unsigned char rnd_mode)
+{
+  int i, n, p, q;
+  mpfr_t u, v, tmpu, tmpv, tmp, zero,prec;
+
+
+  /* If a or b is NaN, let's return NaN */
+  if (FLAG_NAN(a) || FLAG_NAN(b)) 
+    { SET_NAN(r); return; }
+
+
+  /* If a or b is negative, let's return NaN */
+  if ((SIGN(a)<0)||(SIGN(b)<0))
+    { SET_NAN(r); return; }
+
+  mpfr_init(zero);
+  mpfr_set_si(zero,0,GMP_RNDZ);
+
+  
+  /* If a or b is 0, let's return 0 set to the precision of the result */
+  if ((mpfr_cmp(a,zero)==0)||(mpfr_cmp(b,zero)==0)) 
+    {  mpfr_set(r,zero,GMP_RNDZ);
+    return;
+    }
+
+ /* precision of the following calculus */
+  q = PREC(r);
+  p = q + 12;
+ 
+  n = floor(log(p)/log(2));
+
+
+  /* Initialisations */
+  mpfr_init2(u,p);
+  mpfr_init2(v,p);
+
+  if (mpfr_cmp(a,b) >= 0) { 
+    mpfr_set(u,b,GMP_RNDZ);
+    mpfr_set(v,a,GMP_RNDU);
+  }
+  else {
+    mpfr_set(u,a,GMP_RNDZ);
+    mpfr_set(v,b,GMP_RNDU);
+  }
+
+  mpfr_init2(tmp,p);
+  mpfr_init2(tmpu,p);
+  mpfr_init2(tmpv,p);
+
+
+  /* Main loop */
+  for(i=0;i<n;i++) {
+    mpfr_mul(tmp,u,v,GMP_RNDZ);
+    mpfr_sqrt(tmpu,tmp,GMP_RNDZ);
+    mpfr_add(tmp,u,v,GMP_RNDU);
+    mpfr_div_2exp(tmpv,tmp,1,GMP_RNDU);
+    
+    mpfr_set(u,tmpu,GMP_RNDZ);
+    mpfr_set(v,tmpv,GMP_RNDU);
+  }
+
+  /* Setting of the result */
+  mpfr_set(r,v,GMP_RNDN);
+
+  /* Exactness of the result */
+  mpfr_init2(prec,p);
+  mpfr_sub(prec,v,u,GMP_RNDU);
+  mpfr_div(tmp,prec,u,GMP_RNDU);
+  mpfr_div_2exp(prec,tmp,q,GMP_RNDU);
+
+  /* printf("entre u et v : %i ulp\n",(int) floor(mpfr_get_d(prec))); */
+  
+
+
+  /* Let's clean */
+  mpfr_clear(u);
+  mpfr_clear(v);
+  mpfr_clear(tmpu);
+  mpfr_clear(tmpv);
+  mpfr_clear(tmp);
+  mpfr_clear(zero);
+  mpfr_clear(prec);
+
+  return ;
+}
diff --git a/clear.c b/clear.c
new file mode 100644
index 000000000..2a09a2de8
--- /dev/null
+++ b/clear.c
@@ -0,0 +1,14 @@
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "mpfr.h"
+
+void
+#if __STDC__
+mpfr_clear (mpfr_ptr m)
+#else
+mpfr_init (m)
+     mpfr_ptr m;
+#endif
+{
+  (*_mp_free_func) (m->_mp_d, ((m->_mp_prec>>3) + 1));
+}
diff --git a/cmp.c b/cmp.c
new file mode 100644
index 000000000..265c86e21
--- /dev/null
+++ b/cmp.c
@@ -0,0 +1,123 @@
+/* (c) PolKA project at Inria Lorraine
+
+written by Paul Zimmermann, February 1999
+
+    returns 0 iff b = c
+            a positive value iff b > c
+            a negative value iff b < c
+
+More precisely, in case b and c are of same sign, the absolute value 
+of the result is one plus the absolute difference between the exponents 
+of b and c, i.e. one plus the number of bits shifts to align b and c
+(this value is useful in mpfr_sub).
+
+*/
+
+#include <stdio.h>
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "longlong.h"
+#include "mpfr.h"
+
+int mpfr_cmp ( mpfr_srcptr b, mpfr_srcptr c )
+{
+   long int s, diff_exp;
+   unsigned long bn, cn;
+   mp_limb_t *bp, *cp;
+
+   s = SIGN(b) * SIGN(c);
+   if (s<0) return(SIGN(b));
+
+   /* now signs are equal */
+
+   diff_exp = EXP(b)-EXP(c);
+   s = (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)/mp_bits_per_limb+1;
+   cn = (PREC(c)-1)/mp_bits_per_limb+1;
+   bp = MANT(b); cp = MANT(c);
+
+   while (bn && cn) {
+     if (bp[--bn] != cp[--cn])
+       return((bp[bn]>cp[cn]) ? s : -s);
+   }
+
+   if (bn) { while (bn) if (bp[--bn]) return(s); }
+   else if (cn) while (cn) if (cp[--cn]) return(-s);
+
+   return 0;
+}
+
+/* returns the number of cancelled bits when one subtracts abs(c) from abs(b). 
+   Assumes b>=c, which implies EXP(b)>=EXP(c).
+   if b=c, returns prec(b).
+*/
+int mpfr_cmp2 ( mpfr_srcptr b, mpfr_srcptr c )
+{
+   long int d, bn, cn, k;
+   mp_limb_t *bp, *cp, t, u=0, cc=0;
+
+   if (NOTZERO(c)==0) return 0;
+   d = EXP(b)-EXP(c);
+   k = 0; /* result can be d or d+1 if d>1, or >= d otherwise */
+#ifdef DEBUG
+   if (d<0) { printf("assumption EXP(b)<EXP(c) violated\n"); exit(1); }
+#endif
+   bn = (PREC(b)-1)/mp_bits_per_limb;
+   cn = (PREC(c)-1)/mp_bits_per_limb;
+   bp = MANT(b); cp = MANT(c);
+   /* subtracts c from b from most significant to less significant limbs,
+      and first determines first non zero limb difference */
+   if (d) {
+     cc = bp[bn--];
+     if (d<mp_bits_per_limb)
+       cc -= cp[cn]>>d; /* cannot be zero since b is normalized */
+   }
+   else { /* d=0 */
+     while (bn>=0 && cn>=0 && (cc=(bp[bn--]-cp[cn--]))==0) {
+       k+=mp_bits_per_limb;
+     }
+     if (cc==0) { /* either bn<0 or cn<0 */
+       while (bn>=0 && (cc=bp[bn--])==0) k+=mp_bits_per_limb;
+     }
+     /* now bn<0 or cc<>0 */
+     if (cc==0 && bn<0) return(PREC(b));
+   }
+   /* the first non-zero limb difference is cc, and the number
+      of cancelled bits in the upper limbs is k */
+   count_leading_zeros(u, cc);
+   k += u;
+   if (cc != (1<<(mp_bits_per_limb-u-1))) return k;
+   /* now cc is an exact power of two */
+   cc = mp_bits_per_limb-u;
+   while (bn>=0) {
+     /* computes limb bn of difference.
+	cc is previous borrow: either 0 or 1.
+     */
+     t = bp[bn--];
+     if (d<mp_bits_per_limb) {
+       if (d) {
+	 u = cp[cn--]<<(mp_bits_per_limb-d);
+	 if (cn>=0) u+=(cp[cn]>>d); else return(k);
+       }
+       else u = cp[cn--];
+       if (t>u) return(k); /* no borrow possible */
+       t -= u;
+       if (t>1) { /* t=0 when diff=0, t=1 when b=0 and c=~0 */
+	 count_leading_zeros(u, t);
+	 if (u!=~(mp_limb_t)0) return(k+cc+u); /* borrow = 1 */
+	 else cc+=mp_bits_per_limb;
+       }
+     }
+     else {
+       d -= mp_bits_per_limb;
+     }
+   }
+   return k;
+}
+
+
diff --git a/cmp_ui.c b/cmp_ui.c
new file mode 100644
index 000000000..d3f769f9f
--- /dev/null
+++ b/cmp_ui.c
@@ -0,0 +1,79 @@
+/* (c) PolKA project at Inria Lorraine */
+
+#include <stdio.h>
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "longlong.h"
+#include "mpfr.h"
+
+/* returns a positive value if b>i,
+           a negative value if b<i,
+           zero if b=i 
+*/
+int mpfr_cmp_ui ( mpfr_srcptr b, unsigned long int i )
+{
+  int e, k, bn; mp_limb_t c, *bp;
+
+  if (SIGN(b)<0) return -1;
+  else if (!NOTZERO(b)) return((i) ? -1 : 0);
+  else { /* b>0 */
+    e = EXP(b); /* 2^(e-1) <= b < 2^e */
+    if (e>mp_bits_per_limb) return 1;
+
+    c = (mp_limb_t) i;
+    count_leading_zeros(k, c);
+    k = mp_bits_per_limb - k; /* 2^(k-1) <= i < 2^k */
+    if (k!=e) return (e-k);
+
+    /* now k=e */
+    c <<= (mp_bits_per_limb-k);
+    bn = (PREC(b)-1)/mp_bits_per_limb;
+    bp = MANT(b) + bn;
+    if (*bp>c) return 1;
+    else if (*bp<c) return -1;
+
+    /* most significant limbs agree, check remaining limbs from b */
+    while (--bn>=0)
+      if (*--bp) return 1;
+    return 0;
+  }
+}
+
+/* returns a positive value if b>i,
+           a negative value if b<i,
+           zero if b=i 
+*/
+int mpfr_cmp_si ( mpfr_srcptr b, long int i )
+{
+  int e, k, bn, si; mp_limb_t c, *bp;
+
+  si = (i<0) ? -1 : 1; /* sign of i */
+  if (SIGN(b)*i<0) return SIGN(b); /* both signs differ */
+  else if (NOTZERO(b)*i==0) { /* one is zero */
+    if (i==0) return ((NOTZERO(b)) ? 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 */
+    if (e>mp_bits_per_limb) return si;
+
+    c = (mp_limb_t) ((i<0) ? -i : i);
+    count_leading_zeros(k, c);
+    k = mp_bits_per_limb - k; /* 2^(k-1) <= i < 2^k */
+    if (k!=e) return (si*(e-k));
+
+    /* now k=e */
+    c <<= (mp_bits_per_limb-k);
+    bn = (PREC(b)-1)/mp_bits_per_limb;
+    bp = MANT(b) + bn;
+    if (*bp>c) return si;
+    else if (*bp<c) return -si;
+
+    /* most significant limbs agree, check remaining limbs from b */
+    while (--bn>=0)
+      if (*--bp) return si;
+    return 0;
+  }
+}
+
diff --git a/div.c b/div.c
new file mode 100644
index 000000000..670c65e8a
--- /dev/null
+++ b/div.c
@@ -0,0 +1,95 @@
+#include <math.h>
+#include "gmp.h"
+#include "mpfr.h"
+
+/* #define DEBUG
+#define DEBUG2 */
+
+/* q <- n/d using Goldschmidt's iteration */
+void mpfr_div(q, n, d, rnd_mode) 
+mpfr_ptr q; mpfr_srcptr n, d; unsigned char rnd_mode;
+{
+  mpfr_t eps, tmp, one; int expd, i, prec, precq, sh, guard, err;
+  mp_limb_t cc;
+
+#ifdef DEBUG
+  printf("enter mpfr_div, prec(q)=%d n=%1.20e prec(n)=%d d=%1.20e prec(d)=%d rnd=%d\n",PREC(q),mpfr_get_d(n),PREC(n),mpfr_get_d(d),PREC(d),rnd_mode); 
+  printf("n="); mpfr_print_raw(n); putchar('\n');
+  printf("d="); mpfr_print_raw(d); putchar('\n');
+#endif
+  if ((MANT(n)[(PREC(n)-1)/mp_bits_per_limb] & 
+      ((mp_limb_t)1<<(mp_bits_per_limb-1)))==0) {
+    printf("Error in mpfr_div: n is not msb-normalized\n"); exit(1);
+  }
+  if (FLAG_NAN(n) || FLAG_NAN(d)) {
+#ifdef DEBUG
+    printf("dividend or divisor is NaN\n");
+#endif
+    SET_NAN(q); return;
+  }
+  prec = precq = PREC(q);
+  for (i=0;i<2;i++)
+    prec = precq + (int) ceil(log(2.0*ceil(log((double)prec)/log(2.0))+7.0)/
+			   log(2.0));
+  err = prec-precq; /* the error is at most 2^err ulp */
+  prec++; /* add one bit otherwise mfpr_can_round will always fail */
+  prec = prec-mp_bits_per_limb;
+  do {
+    prec += mp_bits_per_limb;
+#ifdef DEBUG2
+  printf("PREC(q)=%d prec=%d\n",precq,prec);
+  printf("n=%1.20e d=%1.20e rnd=%d\n",mpfr_get_d(n),mpfr_get_d(d),rnd_mode); 
+#endif
+  mpfr_set_prec(q, prec, GMP_RNDZ);
+  mpfr_set(q, n, GMP_RNDZ);
+  mpfr_init2(eps, prec); mpfr_init2(tmp, prec); mpfr_init2(one, prec);
+  expd = EXP(d);
+  mpfr_set_ui(one, 1, GMP_RNDZ); 
+  mpfr_mul_2exp(eps, one, expd, GMP_RNDZ); /* eps = 2^expd */
+  if (SIGN(d)<0) mpfr_add(tmp, eps, d, GMP_RNDZ);
+  else mpfr_sub(tmp, eps, d, GMP_RNDZ);
+  i=0; while (i<ABSSIZE(tmp) && MANT(tmp)[i]==MANT(d)[i]) i++;
+  if (i==ABSSIZE(tmp)) {
+    /* if tmp=abs(d), then eps=2*d whence d is an exact power of two */
+    PREC(q) = precq;
+    mpfr_set(q, n, rnd_mode);
+    EXP(q) -= expd-1;
+    if (SIGN(d)<0) CHANGE_SIGN(q);
+    return;
+  }
+  mpfr_mul_2exp(eps, tmp, -expd, GMP_RNDZ);
+  for (;;) {
+#ifdef DEBUG2
+    printf("q=%1.20e eps=%1.20e\n",mpfr_get_d(q),mpfr_get_d(eps));
+#endif
+    /* q[n+1] = q[n] * (1 + e[n])
+       Numerically, it's better to compute q + (q*e) than q * (1+e).
+       Rounding towards zero guarantees we get a lower bound of n/d.
+     */
+    mpfr_mul(tmp, q, eps, GMP_RNDZ);
+    mpfr_add(one, q, tmp, GMP_RNDZ);
+    mpfr_set(q, one, GMP_RNDZ);
+    /* e[n+1] = e[n]^2 */
+    mpfr_mul(tmp, eps, eps, GMP_RNDZ);
+    mpfr_set(eps, tmp, GMP_RNDZ);
+    if (-EXP(eps)>prec) break; /* further iterations won't change
+					   q since we round towards zero. */
+  }
+     if (SIGN(d)<0) CHANGE_SIGN(q);
+     cc = mpfr_can_round(q, prec-err, GMP_RNDZ, rnd_mode, precq);
+     if (cc==0) {
+#ifdef DEBUG
+       printf("not enough precision\n");
+#endif
+       if (prec>2*precq) { printf("does not converge\n"); exit(1); }
+     }
+  } while (cc==0);
+  mpfr_round(q, rnd_mode, precq);
+  mpfr_mul_2exp(q, q, -expd, GMP_RNDZ);
+  mpfr_set_prec(q, precq, rnd_mode);
+  mpfr_clear(eps); mpfr_clear(tmp); mpfr_clear(one);
+#ifdef DEBUG
+  printf("q = %1.20e\n", mpfr_get_d(q));
+  printf("n/d=%1.20e\n", mpfr_get_d(n)/mpfr_get_d(d));
+#endif
+}
diff --git a/div_2exp.c b/div_2exp.c
new file mode 100644
index 000000000..1c3f6d79f
--- /dev/null
+++ b/div_2exp.c
@@ -0,0 +1,13 @@
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "mpfr.h"
+
+void
+mpfr_div_2exp(mpfr_ptr y, mpfr_srcptr x, unsigned long int n, unsigned char rnd_mode)
+{
+  /* Important particular case */ 
+  if (y != x) mpfr_set(y, x, rnd_mode);
+  EXP(y) -= n; 
+  return; 
+}
+
diff --git a/from_Torbjorn b/from_Torbjorn
new file mode 100644
index 000000000..147c7b03f
--- /dev/null
+++ b/from_Torbjorn
@@ -0,0 +1,28 @@
+Return-Path: <tege@matematik.su.se>
+X-Address: Department of Mathematics, Stockholm University 
+	      S-106 91  Stockholm
+	      SWEDEN
+X-Phone: int+46 8 162000
+X-Fax:   int+46 8 6126717
+X-Url:   http://www.matematik.su.se
+To: Paul.Zimmermann@loria.fr (Paul Zimmermann)
+Subject: Re: Paul Zimmermann: changing rounding mode on Alpha with gcc 
+In-reply-to: Your message of "Wed, 17 Feb 1999 17:44:33 +0100."
+             <199902171644.RAA03746@leopold.loria.fr> 
+Date: Wed, 17 Feb 1999 19:21:05 +0100
+From: Torbjorn Granlund <tege@matematik.su.se>
+
+  I once thought I had found a bug, but it was only the peculiar rounding
+  to nearest rule when we are just in the middle of two machine number,
+  where we round to the even lsb instead of rounding to infinity as we do
+  in math.
+
+Been there...  ;-) Note that it is simple to get this right with a single
+logical expression.  Perhaps this code,
+
+    roundup = (((m >> (EXP_BITS - 1)) & 1)
+               & (((m & STICKY_MASK) != 0) | ((m >> EXP_BITS) & 1)));
+
+from my P754 library shows the idea.
+
+Tobjörn
diff --git a/get_str.c b/get_str.c
new file mode 100644
index 000000000..60905159a
--- /dev/null
+++ b/get_str.c
@@ -0,0 +1,115 @@
+#include <stdio.h>
+#include <math.h>
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "mpfr.h"
+
+/*
+  Convert op to a string in base 'base' with 'n' digits and writes the 
+  mantissa in 'str', the exponent in 'expptr'.
+  The format is 0.xxxxxxxxEyyyy.
+  The result is rounded wrt 'rnd_mode'.
+ */
+
+/* #define DEBUG */
+
+char *mpfr_get_str(char *str, char *expptr, int base, size_t n,
+		  mpfr_srcptr op, unsigned char rnd_mode)
+{
+  double d; long e, f, q, i, neg, p, err, prec, sh; mpfr_t a, b; mpz_t bz;
+  char *str0; unsigned char rnd1;
+
+#ifdef DEBUG
+  printf("op="); mpfr_print_raw(op); printf(" rnd_mode=%d\n",rnd_mode);
+  printf("  =%1.20e\n",mpfr_get_d(op));
+#endif
+  /* first determines the exponent */
+  e = EXP(op); 
+  EXP(op)=0; d=fabs(mpfr_get_d(op)); EXP(op)=e;
+  /* 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
+     i.e. f = 1 + floor(log(|op|)/log(base))
+     = 1 + floor((log(|m|)+e*log(2))/log(base)) */
+  f = 1 + (int) floor((log(d)+(double)e*log(2.0))/log((double)base));
+#ifdef DEBUG
+  printf("exponent = %d\n",f);
+#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));
+  err = 5;
+  q = prec+err;
+  /* one has to use at least q bits */
+  q = ((q-1)/mp_bits_per_limb)*mp_bits_per_limb;
+  mpfr_init(a); mpfr_init(b);
+  p = n-f; if ((neg=(p<0))) p=-p;
+  rnd1 = (neg) ? GMP_RNDU : GMP_RNDZ; /* if neg we divide by base^p */
+  do {
+    q += mp_bits_per_limb;
+    /* compute base^p with q bits and rounding towards zero */
+    mpfr_set_prec(b, q, GMP_RNDZ);
+    if (p==0) {
+      mpfr_set(b, op, GMP_RNDZ);
+    } 
+    else {
+      mpfr_set_prec(a, q, rnd1); 
+      mpfr_set_ui(a, base, rnd1);
+      for (i=0;(1<<i)<=p;i++);
+      /* now 2^(i-1) <= p < 2^i */
+       for (i-=2; i>=0; i--) {
+	 mpfr_mul(b, a, a, rnd1);
+	 if (p & (1<<i)) mpfr_mul_ui(a, b, base, rnd1);
+	 else mpfr_set(a, b, rnd1);
+       }
+       /* now a is an approximation by default of base^p */
+       if (neg) mpfr_div(b, op, a, GMP_RNDZ);
+       else mpfr_mul(b, op, a, GMP_RNDZ);
+    }
+    if (SIGN(op)<0) CHANGE_SIGN(b);
+  } while (mpfr_can_round(b, q-err, GMP_RNDZ, rnd_mode, prec)==0);
+  if (SIGN(op)<0)
+    switch (rnd_mode) {
+    case GMP_RNDU: rnd_mode=GMP_RNDZ; break;
+    case GMP_RNDD: rnd_mode=GMP_RNDU; break;
+  }
+#ifdef DEBUG
+printf("rnd=%d\n",rnd_mode);
+printf("b="); mpfr_print_raw(b); putchar('\n');
+printf("=%1.20e\n",mpfr_get_d(b));
+#endif
+  prec=EXP(b); 
+  mpfr_round(b, rnd_mode, prec);
+  prec=EXP(b); /* may have chnaged due to rounding */
+#ifdef DEBUG
+printf("b="); mpfr_print_raw(b); putchar('\n');
+printf("prec=%d q=%d b=",prec,q); mpfr_print_raw(b); putchar('\n');
+printf("=%1.20e\n",mpfr_get_d(b));
+#endif
+  /* now the mantissa is the integer part of b */
+  mpz_init(bz); q=1+(prec-1)/mp_bits_per_limb; _mpz_realloc(bz, q);
+  sh = prec%mp_bits_per_limb;
+  if (sh) mpn_rshift(PTR(bz), MANT(b), q, mp_bits_per_limb-sh);
+  else MPN_COPY(PTR(bz), MANT(b), q);
+  bz->_mp_size=q;
+#ifdef DEBUG
+printf("bz="); mpz_out_str(stdout,10,bz); putchar('\n');
+printf("b="); mpfr_print_raw(b); putchar('\n');
+#endif
+  /* computes the number of characters needed */
+  q = ((SIGN(op)<0) ? 1 : 0) + 2 + n + 2 +
+    + (int) ceil(log((double)fabs(f))/log(10.0));
+  if (f<10 || f>-10) q++;
+  if (str==NULL) str0=str=malloc(q);
+  if (SIGN(op)<0) *str++='-';
+  if (n>1) *str++ = '.';
+  mpz_get_str(str, base, bz); /* n digits of mantissa */
+  if (strlen(str)==n+1) f++; /* possible due to rounding */
+  str[n++] = 'e'; 
+  f--; /* replaces 0.xxx*b^f by x.xx*b^(f-1) */
+  str[n++] = (f>=0) ? '+' : '-'; /* is there a rule for f=0 ? */
+  if (f<10 && f>-10) str[n++]='0';
+  sprintf(str+n, "%ld", (f<0) ? -f : f);
+  if (str[-1]=='.') { str[-1]=str[0]; str[0]='.'; }
+  mpfr_clear(a); mpfr_clear(b); mpz_clear(bz);
+  return str0;
+}
diff --git a/init.c b/init.c
new file mode 100644
index 000000000..fd5d31cba
--- /dev/null
+++ b/init.c
@@ -0,0 +1,28 @@
+#include <stdio.h>
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "mpfr.h"
+
+void
+#if __STDC__
+mpfr_init2 (mpfr_t x, unsigned long int p)
+#else
+mpfr_init2 (x, p)
+     mpfr_t x;
+     unsigned long int p;
+#endif
+{
+  unsigned long xsize; 
+
+  if (p==0) {
+    printf("*** cannot initialize mpfr with precision 0\n"); exit(1);
+  }
+
+  xsize = (p - 1)/BITS_PER_MP_LIMB + 1; 
+
+  x -> _mp_prec = p;
+  x -> _mp_d = (mp_ptr) (*_mp_allocate_func) 
+    (xsize * BYTES_PER_MP_LIMB);
+  x -> _mp_size = xsize; 
+  x -> _mp_exp = 0; 
+}
diff --git a/init_set.h b/init_set.h
new file mode 100644
index 000000000..f779fee12
--- /dev/null
+++ b/init_set.h
@@ -0,0 +1,35 @@
+#define mpfr_init_set_si(x, i, p, rnd) {\
+ mpfr_init2((x), (p)); \
+ mpfr_set_si((x), (i), (rnd)); \
+}
+
+#define mpfr_init_set_ui(x, i, p, rnd) {\
+ mpfr_init2((x), (p)); \
+ mpfr_set_ui((x), (i), (rnd)); \
+}
+
+#define mpfr_init_set_d(x, d, p, rnd) {\
+ mpfr_init2((x), (p)); \
+ mpfr_set_d((x), (d), (rnd)); \
+}
+
+#define mpfr_init_set(x, y, p, rnd) {\
+ mpfr_init2((x), (p)); \
+ mpfr_set((x), (y), (rnd)); \
+}
+
+#define mpfr_init_set_f(x, y, p, rnd) {\
+ mpfr_init2((x), (p)); \
+ mpfr_set_f((x), (y), (rnd)); \
+}
+
+#define mpfr_init_set_str(x, y, p, rnd) {\
+ mpfr_init2((x), (p)); \
+ mpfr_set_str((x), (y), (rnd)); \
+}
+
+#define mpfr_init_set_str_raw(x, y, p, rnd) {\
+ mpfr_init2((x), (p)); \
+ mpfr_set_str_raw((x), (y), (rnd)); \
+}
+
diff --git a/mmpfr b/mmpfr
new file mode 100755
index 000000000..cfeacfdf4
--- /dev/null
+++ b/mmpfr
@@ -0,0 +1,34 @@
+#!/bin/sh
+
+SYST=`uname -a | awk '{print $1$3}'`
+
+special=
+
+case $SYST in
+    SunOS5*)
+        dirname=solaris;;
+    Linux*)
+        special="-DLIBC211"; # for __setfpucw
+        dirname=linux;;
+    OSF*)
+        special="-mfp-rounding-mode=d -mieee-with-inexact";
+        dirname=alpha;;
+    SunOS4*)
+	dirname=sunos;;
+	*)
+        echo gmp n\'est pas compile pour le systeme $SYST. 
+        exit 0;;
+esac
+
+GMPDIR=/users/polka/logiciels/gmp/$dirname
+ODIR=/users/polka/logiciels/mpfr/C/msb_norm/o.$dirname
+export GMPDIR ODIR dirname special
+if [ "$1" = "tests" ]
+then 
+    cd tests; make all
+elif [ "$1" = "" ]
+then
+    make libmpfr.a
+else
+    make $*
+fi
diff --git a/mpfr.h b/mpfr.h
new file mode 100644
index 000000000..5da9b29b0
--- /dev/null
+++ b/mpfr.h
@@ -0,0 +1,134 @@
+/* Definition of rounding modes */
+
+#define GMP_RNDN 0
+#define GMP_RNDZ 1
+#define GMP_RNDU 2
+#define GMP_RNDD 3
+
+/* Definitions of types and their semantics */
+
+typedef struct {  
+  unsigned long int _mp_prec; /* WARNING : for the mpfr type, the precision */
+                              /* should be understood as the number of BITS,*/
+			      /* not the number of mp_limb_t's. This means  */
+			      /* that the corresponding number of allocated
+				 limbs is >= ceil(_mp_prec/BITS_PER_MP_LIMB) */
+  mp_size_t _mp_size;         /* abs(_mp_size) is the number of allocated 
+				 limbs the field _mp_d points to.
+				 The sign is that of _mp_size.
+				 The number 0 is such that _mp_d[k-1]=0
+				 where k = ceil(_mp_prec/BITS_PER_MP_LIMB) */
+  mp_exp_t _mp_exp; 
+  mp_limb_t *_mp_d;
+}
+__mpfr_struct; 
+
+/*
+   The number represented is
+
+    sign(_mp_size)*(_mp_d[k-1]/B+_mp_d[k-2]/B^2+...+_mp_d[0]/B^k)*2^_mp_exp
+
+   where k=ceil(_mp_prec/BITS_PER_MP_LIMB) and B=2^BITS_PER_MP_LIMB.
+
+   For the msb (most significant bit) normalized representation, we must have
+   _mp_d[k-1]>=B/2, unless the number is zero (in that case its sign is still
+   given by sign(_mp_size)).
+
+   We must also have the last k*BITS_PER_MP_LIMB-_mp_prec bits set to zero.
+*/
+
+typedef __mpfr_struct mpfr_t[1]; 
+typedef __mpfr_struct *mpfr_ptr; 
+typedef __gmp_const __mpfr_struct *mpfr_srcptr;
+
+
+/* Prototypes */
+
+#ifndef _PROTO
+#if defined (__STDC__) || defined (__cplusplus)
+#define _PROTO(x) x
+#else
+#define _PROTO(x) ()
+#endif
+#endif
+
+/* bit 31 of _mp_size is used for sign,
+   bit 30 of _mp_size is used for Nan flag,
+   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 SIGN(x) (((x)->_mp_size >> 31) ? -1 : 1)
+#define ISNONNEG(x) (SIGN(x)>=0)
+#define ISNEG(x) (SIGN(x)==-1)
+#define CHANGE_SIGN(x) (SIZE(x) = SIZE(x) ^ (1<<31))
+#define PREC(x) ((x)->_mp_prec)
+#define NOTZERO(x) (ABSSIZE(x))
+
+/* 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, (x)->_mp_prec>>3 + 1, (q)>>3 + 1); \
+    (x)->_mp_prec = q; }
+
+void mpfr_init2 _PROTO ((mpfr_ptr, unsigned long int));
+int mpfr_round_raw _PROTO ((mp_limb_t *, mp_limb_t *, char, long,
+			     unsigned long));
+void mpfr_round _PROTO ((mpfr_ptr, char, unsigned long)); 
+int mpfr_can_round _PROTO ((mpfr_ptr, unsigned int, unsigned char, 
+			    unsigned char, unsigned long));
+void mpfr_set_d _PROTO ((mpfr_ptr, double, unsigned char)); 
+double mpfr_get_d _PROTO ((mpfr_ptr)); 
+void mpfr_set_f _PROTO ((mpfr_ptr, mpf_srcptr, unsigned long, char)); 
+void mpfr_set_si _PROTO ((mpfr_ptr, long, unsigned char)); 
+void mpfr_set_ui _PROTO ((mpfr_ptr, unsigned long, unsigned char)); 
+void mpfr_print_raw _PROTO ((mpfr_srcptr)); 
+void mpfr_clear _PROTO ((mpfr_ptr)); 
+void mpfr_set_str_raw _PROTO ((mpfr_ptr, char *)); 
+void mpfr_get_str_raw _PROTO ((char *, mpfr_srcptr));
+char* mpfr_get_str _PROTO ((char *, char *, int, size_t, mpfr_srcptr, unsigned char));
+void mpfr_mul _PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, unsigned char));
+void mpfr_div _PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, unsigned char));
+void mpfr_agm _PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, unsigned char));
+void mpfr_sqrt _PROTO ((mpfr_ptr, mpfr_srcptr, unsigned char));
+void mpfr_add _PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, unsigned char));
+void mpfr_sub _PROTO ((mpfr_ptr, mpfr_srcptr, mpfr_srcptr, unsigned char));
+void mpfr_set _PROTO ((mpfr_ptr, mpfr_srcptr, unsigned char));
+void mpfr_set_machine_rnd_mode _PROTO ((unsigned char));
+int mpfr_cmp _PROTO ((mpfr_srcptr, mpfr_srcptr));
+int mpfr_cmp_ui _PROTO ((mpfr_srcptr, unsigned long int));
+int mpfr_cmp_si _PROTO ((mpfr_srcptr, long int));
+int mpfr_cmp2 _PROTO ((mpfr_srcptr, mpfr_srcptr));
+void mpfr_mul_2exp _PROTO((mpfr_ptr, mpfr_srcptr, unsigned long int,unsigned char));
+void mpfr_div_2exp _PROTO((mpfr_ptr, mpfr_srcptr, unsigned long int,unsigned char));
+void mpfr_set_prec _PROTO((mpfr_ptr, unsigned long int,unsigned char));
+extern mp_size_t __gmp_default_fp_bit_precision;
+extern char __gmp_default_rounding_mode;
+
+#define mpfr_init(x) mpfr_init2(x, __gmp_default_fp_bit_precision)
+
+#if (BITS_PER_MP_LIMB==32)
+#define LOG_MP_BITS_PER_LIMB 5
+#define MPFR_LIMBS_PER_DOUBLE 2
+#elif (BITS_PER_MP_LIMB==64)
+#define LOG_MP_BITS_PER_LIMB 6
+#define MPFR_LIMBS_PER_DOUBLE 1
+#endif
+
+#define mpfr_init_set_si(x, i, p, rnd) \
+ mpfr_init2((x), (p)); mpfr_set_si((x), (i), (rnd)); 
+#define mpfr_init_set_ui(x, i, p, rnd) \
+ mpfr_init2((x), (p)); mpfr_set_ui((x), (i), (rnd));
+#define mpfr_init_set_d(x, d, p, rnd) \
+ mpfr_init2((x), (p)); mpfr_set_d((x), (d), (rnd));
+#define mpfr_init_set(x, y, p, rnd) \
+ mpfr_init2((x), (p)); mpfr_set((x), (y), (rnd)); 
+#define mpfr_init_set_f(x, y, p, rnd) \
+ mpfr_init2((x), (p)); mpfr_set_f((x), (y), (rnd));
+#define mpfr_init_set_str(x, y, p, rnd) \
+ mpfr_init2((x), (p)); mpfr_set_str((x), (y), (rnd)); 
+#define mpfr_init_set_str_raw(x, y, p, rnd) \
+ mpfr_init2((x), (p)); mpfr_set_str_raw((x), (y), (rnd)); 
diff --git a/mul.c b/mul.c
new file mode 100644
index 000000000..1c800ad73
--- /dev/null
+++ b/mul.c
@@ -0,0 +1,161 @@
+/* written by Paul Zimmermann, November 1998-January 1999 */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "mpfr.h"
+
+/* #define DEBUG */ /* check possible errors */
+/* #define DEBUG2 */ /* print informations */
+
+/* avoids two divisions */
+#define GETBIT(a, n) (a)[(n)>>LOG_MP_BITS_PER_LIMB]>>((n)&((1<<LOG_MP_BITS_PER_LIMB)-1))
+
+/* and here comes the code. To do:
+- take sign(b) and sign(c) into account, and write sign(a)
+- does not use all bits of b and c when PREC(b)>PREC(a) or PREC(c)>PREC(a)
+- use fast mpn_mul multiplication
+*/
+
+void mpfr_mul(a, b, c, rnd_mode) 
+mpfr_ptr a; mpfr_srcptr b, c; unsigned char rnd_mode;
+{
+  unsigned int bn, cn, an, k; int sh; unsigned char b1;
+  mp_limb_t *ap, *bp, *cp, cc;
+  long int sign_product;
+
+  /* the multiplication works as follows:
+     1. multiply all significant limbs from the mantissa of b and c in a 
+        temporary space (ap)
+     2. rounds the result up to PREC(a) bits
+     3. copies the result into the destination
+  */
+
+#ifdef DEBUG2
+  printf("b="); mpfr_print_raw(b); putchar('\n');
+  printf("c="); mpfr_print_raw(c); putchar('\n');
+#endif
+  sign_product = SIGN(b) * 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 */
+  k = bn+cn; /* effective nb of limbs used by b*c */
+  ap = (mp_limb_t*) alloca(k*BYTES_PER_MP_LIMB);
+  bp = MANT(b); cp = MANT(c);
+
+  /* step 1: multiplies two mantissa */
+  if (bn>=cn) mpn_mul(ap, bp, bn, cp, cn);
+  else mpn_mul(ap, cp, cn, bp, bn);
+
+  /* now ap[0]..ap[k-1] contains the product of both mantissa,
+     with ap[k-1]>=2^(mp_bits_per_limb-2) */
+  an = (PREC(a)-1)/mp_bits_per_limb+1; /* number of significant limbs of a */
+  sh = k*mp_bits_per_limb-PREC(a); /* number of bits to truncate */
+  b1 = GETBIT(ap+k-1,mp_bits_per_limb-1); /* msb from the product */
+  sh = sh+b1-1; /* sh-- if b1=0 */
+#ifdef DEBUG2
+printf("an=%u bn=%u cn=%u b1=%u sh=%u bp[0]=%lu cp[0]=%lu\n",an,bn,cn,b1,sh,*bp,*cp);
+#endif
+  /* if sh<=0, then no rounding is needed */
+  if (sh<=0) { /* copy most significant k*mp_bits_per_limb bits */
+    bp = a->_mp_d; cp = bp+an-k;
+    if (b1) {
+      MPN_COPY(cp, ap, k); 
+      EXP(a) = EXP(b) + EXP(c);
+    }
+    else { 
+      mpn_lshift(cp, ap, k, 1); 
+      EXP(a) = EXP(b) + EXP(c) - 1;
+    }
+    /* set to zero low bits */
+    while (cp>bp) *--cp=0;
+  }
+  else { /* k*mp_bits_per_limb + (b1-1) > PREC(a) */
+     /* sh>0 is the number of low significant bits to truncate.
+        (i) down mode: just put them to zero
+        (ii) nearest mode: add 2^(sh-1) and truncate
+        (iii) up mode: add 2^sh-1 and truncate
+     */
+     if (rnd_mode==GMP_RNDN) { /* nearest mode */
+       /* if middle of two representable numbers, we must
+          set the lower bit to zero according to the IEEE norm */
+       sh--; /* to keep one guard bit */
+       cc=0; /* 0 iff middle */
+       while (cc==0 && sh>=mp_bits_per_limb) {
+	 cc = *ap++; k--; sh -= mp_bits_per_limb;
+       }
+       ap += (sh/mp_bits_per_limb);
+       k -= (sh/mp_bits_per_limb);
+       sh %= mp_bits_per_limb; /* now 0<=sh<mp_bits_per_limb */
+       if (cc==0)
+	 /* middle iff last sh bits are zero and bit (sh+1) is 1 */
+	 cc = (*ap<<(mp_bits_per_limb-sh-1)) - ((mp_limb_t)1<<(mp_bits_per_limb-1));
+       if (cc)
+	 cc = mpn_add_1(ap, ap, k, (mp_limb_t)1<<sh);
+       else /* put lower (sh+1) bits to zero */
+	 *ap = *ap & ~(((mp_limb_t)2<<sh)-1);
+       if (cc!=0) { /* may happen especially when PREC(a) is small */
+	 mpn_rshift(ap, ap, k, 1);
+	 ap[k-1] += (mp_limb_t)1<<(mp_bits_per_limb-1);
+	 b1++; /* so that EXP(a) will be incremented by 1 below */
+       }
+       else sh++;
+       /* 0 <= sh <= mp_bits_per_limb */
+     }
+     else if ((sign_product>=0 && rnd_mode==GMP_RNDU) ||
+         (sign_product<0 && rnd_mode==GMP_RNDD)) { /* up mode */
+       cc = mpn_sub_1(ap, ap, k, 1); /* subtract 1 */
+       ap += (sh/mp_bits_per_limb);
+       k -= (sh/mp_bits_per_limb);
+       sh %= mp_bits_per_limb; /* now 0<=sh<mp_bits_per_limb */
+       cc = mpn_add_1(ap, ap, k, (mp_limb_t)1<<sh)-cc;
+       if (cc!=0) { /* may happen especially when PREC(a) is small */
+	 mpn_rshift(ap, ap, k, 1);
+	 ap[k-1] += (mp_limb_t)1<<(mp_bits_per_limb-1);
+	 b1++; /* so that EXP(a) will be incremented by 1 below */
+	 if (sh==0) { sh=mp_bits_per_limb; ap--; k++; }
+	 sh--;
+       }
+     }
+     else { /* down mode */
+       ap += (sh/mp_bits_per_limb);
+       k -= (sh/mp_bits_per_limb);
+       sh %= mp_bits_per_limb;
+     }
+     /* now truncate last sh bits in ap[0] */
+#ifdef DEBUG2
+printf("sh=%u *ap=%u\n",sh,*ap);
+#endif
+     if (sh==mp_bits_per_limb) { ap++; k--; sh=0; }
+     else *ap = *ap & (~(mp_limb_t)0 << sh);
+#ifdef DEBUG2
+printf("*ap=%u\n",*ap);
+#endif
+     /* step 3: copies the result into the destination */
+     cp=a->_mp_d; 
+     EXP(a) = EXP(b) + EXP(c) + b1 - 1;
+     if (b1) {
+#ifdef DEBUG
+       if (k!=an) /* normally this should not happen */
+	 { printf("k=%u != an=%u\n",k,an); exit(1); }
+#endif
+       MPN_COPY(cp, ap, k);
+     }
+     else { 
+       if (k!=an) { /* only when k=an+1 and sh=mp_bits_per_limb-1 */
+#ifdef DEBUG
+	 if (k!=an+1 || sh!=mp_bits_per_limb-1) 
+	   { printf("k=%u != an+1=%u or sh=%u != mp_bits_per_limb-1\n",k,an+1,
+		    sh); exit(1); }
+#endif
+	 mpn_lshift(cp, ap+1, an, 1);
+	 *cp += ap[0]>>sh;
+       }
+       else mpn_lshift(cp, ap, an, 1); 
+     }
+   }
+  if (sign_product<0) CHANGE_SIGN(a);
+#ifdef DEBUG2
+printf("b*c="); mpfr_print_raw(a); putchar('\n');
+#endif
+}
+
+     
diff --git a/mul_2exp.c b/mul_2exp.c
new file mode 100644
index 000000000..6c595463d
--- /dev/null
+++ b/mul_2exp.c
@@ -0,0 +1,13 @@
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "mpfr.h"
+
+void
+mpfr_mul_2exp(mpfr_ptr y, mpfr_srcptr x, unsigned long int n, unsigned char rnd_mode)
+{
+  /* Important particular case */ 
+  if (y != x) mpfr_set(y, x, rnd_mode);
+  EXP(y) += n;
+  return;
+}
+
diff --git a/mul_ui.c b/mul_ui.c
new file mode 100644
index 000000000..5976a4633
--- /dev/null
+++ b/mul_ui.c
@@ -0,0 +1,42 @@
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "longlong.h"
+#include "mpfr.h"
+
+void
+mpfr_mul_ui(mpfr_ptr y, mpfr_ptr x, unsigned long u, char RND_MODE)
+     /* on suppose SIZ(y)=SIZ(x) */
+{
+  mp_limb_t carry = 0, *my, *old_my; unsigned long c; 
+  unsigned long xsize, ysize, cnt; 
+
+  my = MANT(y); 
+  ysize = ABSSIZE(y); xsize = ABSSIZE(x); 
+
+  /*  
+      if (ysize < xsize) 
+      {
+      old_my = my; 
+      my = _mp_allocate_func(xsize*BYTES_PER_MP_LIMB);
+      _mp_free_func(old_my, ysize*BYTES_PER_MP_LIMB); 
+      }
+  */
+  carry = mpn_mul_1(my, MANT(x), xsize, u); 
+  count_leading_zeros(cnt, carry); 
+      
+  c = mpfr_round_raw(my, my, RND_MODE, ysize, PREC(y)-BITS_PER_MP_LIMB+cnt);
+  
+  /* If cnt = 1111111111111 and c = 1 we shall get depressed */
+  if (c && (carry == (1UL << (BITS_PER_MP_LIMB - cnt)) - 1))
+    {
+      cnt--; 
+      mpn_rshift(my, my, ysize, BITS_PER_MP_LIMB - cnt); 
+      my[ysize - 1] |= 1 << (BITS_PER_MP_LIMB - 1);
+    }
+  else
+    {
+      mpn_rshift(my, my, ysize, BITS_PER_MP_LIMB - cnt); 
+      my[ysize - 1] |= (carry << cnt); 
+    }
+  EXP(y) = EXP(x) + BITS_PER_MP_LIMB - cnt; 
+}
diff --git a/neg.c b/neg.c
new file mode 100644
index 000000000..6a85f3971
--- /dev/null
+++ b/neg.c
@@ -0,0 +1,13 @@
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "mpfr.h"
+
+void mpfr_neg(mpfr_t a, mpfr_t b, unsigned char rnd_mode)
+{
+  CHANGE_SIGN(b);
+  if (a != b) {
+    mpfr_set(a, b, rnd_mode);
+    CHANGE_SIGN(b);
+  }
+  return; 
+}
diff --git a/o.solaris/.pure b/o.solaris/.pure
new file mode 100644
index 000000000..e69de29bb
--- /dev/null
+++ b/o.solaris/.pure
diff --git a/print_raw.c b/print_raw.c
new file mode 100644
index 000000000..afc2590d6
--- /dev/null
+++ b/print_raw.c
@@ -0,0 +1,54 @@
+#include <stdio.h>
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "mpfr.h"
+
+void
+mpfr_get_str_raw(char *digit_ptr, mpfr_srcptr x)
+{
+  mp_limb_t *mx, wd, t; long ex, sx, k, l, p;
+
+  mx = MANT(x); 
+  ex = EXP(x); 
+  p = PREC(x); 
+
+  if (SIGN(x) < 0) { *digit_ptr = '-'; digit_ptr++; }
+  sprintf(digit_ptr, "0."); digit_ptr += 2; 
+
+  sx = 1+(p-1)/mp_bits_per_limb; /* number of significant limbs */
+  for (k = sx - 1; k >= 0 ; k--)
+    { 
+      wd = mx[k]; 
+      t = 1UL << (BITS_PER_MP_LIMB - 1); 
+      for (l = BITS_PER_MP_LIMB - 1; l>=0; l--)
+	{
+	  if (wd & t) 
+	    { *digit_ptr = '1'; digit_ptr++; } 
+	  else 
+	    { *digit_ptr = '0'; digit_ptr++; }
+	  t >>= 1; 
+	  if (--p==0) { *digit_ptr = '['; digit_ptr++; }
+	}
+    }
+  sprintf(digit_ptr, "]E%ld", ex); 
+}
+ 
+void
+mpfr_print_raw(mpfr_srcptr x)
+{
+  char *str; 
+
+  /* 3 char for sign + 0 + binary point
+     + 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
+   */
+  str = (char *) malloc((17 + ABSSIZE(x) * BITS_PER_MP_LIMB)*sizeof(char)); 
+  mpfr_get_str_raw(str, x);
+
+  printf("%s", str); 
+  free(str); 
+}
+
diff --git a/rnd_mode.c b/rnd_mode.c
new file mode 100644
index 000000000..ad48ce95c
--- /dev/null
+++ b/rnd_mode.c
@@ -0,0 +1,77 @@
+/* functions to set/get the machine rounding mode,
+   this should use the interface from the ISO C9X draft
+   (file <fenv.h>) */
+
+#include <stdio.h>
+#include "gmp.h"
+#include "mpfr.h"
+
+#ifdef IRIX64
+#include <sys/fpu.h>
+#define TOZERO swapRM(ROUND_TO_ZERO)
+#define TOINFP swapRM(ROUND_TO_PLUS_INFINITY)
+#define TONEAREST swapRM(ROUND_TO_NEAREST)
+#define TOINFM swapRM(ROUND_TO_MINUS_INFINITY)
+#elif (defined (solaris) || defined (sun4))
+#include <ieeefp.h>
+#define TOZERO fpsetround(FP_RZ)
+#define TOINFP fpsetround(FP_RP)
+#define TONEAREST fpsetround(FP_RN)
+#define TOINFM fpsetround(FP_RM)
+#elif alpha
+#ifdef __GNUC__
+/* GCC patched include files forget to define those... */
+#define FP_RND_RZ       0
+#define FP_RND_RN       1
+#define FP_RND_RP       2
+#define FP_RND_RM       3
+#endif
+#include <float.h>
+#define TOZERO write_rnd(FP_RND_RZ)
+#define TOINFP write_rnd(FP_RND_RP)
+#define TONEAREST write_rnd(FP_RND_RN)
+#define TOINFM write_rnd(FP_RND_RM)
+#elif AIX
+#include <float.h>
+#define TOZERO fp_swap_rnd(FP_RND_RZ)
+#define TOINFP fp_swap_rnd(FP_RND_RP)
+#define TONEAREST fp_swap_rnd(FP_RND_RN)
+#define TOINFM fp_swap_rnd(FP_RND_RM)
+#elif sunos
+#include <floatingpoint.h>
+char *out;
+#define TOZERO ieee_flags("set","direction","tozero",&out)
+#define TOINFP ieee_flags("set","direction","positive",&out)
+#define TONEAREST ieee_flags("set","direction","nearest",&out)
+#define TOINFM ieee_flags("set","direction","negative",&out)
+#elif hp700
+#define TOZERO fpsetround(FP_RZ)
+#define TOINFP fpsetround(FP_RP)
+#define TONEAREST fpsetround(FP_RN)
+#define TOINFM fpsetround(FP_RM)
+#elif (defined (__i386__) || defined (__i486__) || defined (linux))
+#include <fpu_control.h>
+#ifdef LIBC211
+#define __setfpucw(cw) __asm__ ("fldcw %0" : : "m" (cw))
+#endif
+/* be careful to put bits 9-8 (Precision control) 
+   to 10 (round to double precision) instead of the
+   default 11 (round to extended precision) */
+#define TOZERO __setfpucw(0x1e7f)
+#define TOINFP __setfpucw(0x1a7f)
+#define TOINFM __setfpucw(0x167f)
+#define TONEAREST __setfpucw(0x127f)
+#endif
+
+/* sets the machine rounding mode to the value rnd_mode */
+void mpfr_set_machine_rnd_mode(unsigned char rnd_mode)
+{
+  switch (rnd_mode) {
+  case GMP_RNDN: TONEAREST; break;
+  case GMP_RNDZ: TOZERO; break;
+  case GMP_RNDU: TOINFP; break;
+  case GMP_RNDD: TOINFM; break;
+  default: fprintf(stderr,"invalid rounding mode\n"); exit(1);
+  }
+}
+
diff --git a/round.c b/round.c
new file mode 100644
index 000000000..556677a5d
--- /dev/null
+++ b/round.c
@@ -0,0 +1,197 @@
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "mpfr.h"
+#include <stdio.h>
+
+#ifdef Exp
+#include "longlong.h"
+#endif
+
+/* returns 0 if round(sign*xp[0..xn-1], prec, rnd) = 
+   round(sign*xp[0..xn-1], prec, GMP_RNDZ), 1 otherwise,
+   where sign=1 if neg=0, sign=-1 otherwise.
+
+   Does *not* modify anything.
+*/
+int mpfr_round_raw2(mp_limb_t *xp, unsigned long xn, char neg, char rnd, 
+		    unsigned long prec)
+{
+  unsigned long mask, nw; long wd; char rw; short l;
+
+  nw = prec / BITS_PER_MP_LIMB; rw = prec & (BITS_PER_MP_LIMB - 1); 
+  if (rw) nw++; 
+  if (rnd==GMP_RNDZ || xn<nw || (rnd==GMP_RNDU && neg)
+      || (rnd==GMP_RNDD && neg==0)) return 0;
+  mask = ~((1UL<<(BITS_PER_MP_LIMB - rw)) - 1);
+  switch (rnd)
+    {
+    case GMP_RNDU:
+    case GMP_RNDD:
+      if (xp[xn-nw] & ~mask) return 1;
+      for (l=nw+1;l<=xn;l++)
+	if (xp[xn-l]) break;
+      return (l<=xn);
+    case GMP_RNDN:
+    /* First check if we are just halfway between two representable numbers */
+      wd = xn - nw;
+      if (!rw)
+	{
+	  wd--; 
+	  if ((xp[wd] & (~mask)) == (1UL << (BITS_PER_MP_LIMB - rw - 1)))
+	      do wd--; while (wd > 0 && !xp[wd]);
+
+	  if (wd)
+	      return (xp[xn - nw - 1] & (1 << (BITS_PER_MP_LIMB - 1)));
+	  else
+	      return xp[xn - nw] & 1;
+	}
+      else
+      if (rw + 1 < BITS_PER_MP_LIMB)
+	{
+	  if ((xp[wd] & (~mask)) == (1UL << (BITS_PER_MP_LIMB - rw - 1)))
+	      do { wd--; } while (wd >= 0 && !xp[wd]);
+	  else return (xp[wd] & (1UL << (BITS_PER_MP_LIMB - rw - 1)));
+	  
+	  /* first limb was in the middle, and others down to wd+1 were 0 */
+	  if (wd>=0) return 1;
+	  else
+	      return ((xp[xn - nw] & mask) >> (BITS_PER_MP_LIMB - rw)) & 1;
+	}
+      else
+	/* Modified PZ, 27 May 1999:
+	   rw, i.e. the number of bits stored in xp[xn-nw], is 
+	   BITS_PER_MP_LIMB-1, i.e. there is exactly one non significant bit. 
+	   We are just halway iff xp[wd] has its low significant bit 
+	   set and all limbs xp[0]...xp[wd-1] are zero */
+	{
+	  if (xp[wd] & 1)
+	      do wd--; while (wd >= 0 && !xp[wd]);
+	  return xp[xn-nw] & ((wd<0) ? 2 : 1);
+	}
+    default: return 0;
+    }
+}
+
+/* puts in y the value of xp (which has xsize limbs) rounded to precision
+   xprec and direction RND_MODE */
+int
+mpfr_round_raw(mp_limb_t *y, mp_limb_t *xp, char RND_MODE, 
+	       long xsize, unsigned long xprec)
+{
+  unsigned long nw, mask;
+  char negative = 0, rw, carry = 0;
+
+  if (xsize == 0) { return 0; }
+  /* warning: -xsize is not correct since some bits are used for Nan, exact */
+  if (xsize < 0) { negative = 1; xsize = xsize ^ (1<<31); }
+
+#ifdef Exp
+  count_leading_zeros(flag, xp[xsize-1]); 
+  xprec += flag; 
+#endif
+
+  nw = xprec / BITS_PER_MP_LIMB; rw = xprec & (BITS_PER_MP_LIMB - 1); 
+  if (rw) nw++; 
+  /* number of words needed to represent x */
+  /* precision is larger than the size of x. No rounding is necessary. 
+     Just add zeroes at the end */
+  if (xsize < nw) { 
+    MPN_COPY(y + nw - xsize, xp, xsize); 
+    MPN_ZERO(y, nw - xsize); /* PZ 27 May 99 */
+    return 0; 
+  }
+  
+  mask = ~((1UL<<(BITS_PER_MP_LIMB - rw)) - 1); 
+
+  if (mpfr_round_raw2(xp, xsize, negative, RND_MODE, xprec))
+    carry = mpn_add_1(y, xp + xsize - nw, nw,
+                          1UL << (BITS_PER_MP_LIMB - rw));
+  else MPN_COPY(y, xp + xsize - nw, nw);
+
+  y[0] &= mask;
+  return carry; 
+}
+
+void
+mpfr_round(mpfr_t x, char RND_MODE, unsigned long prec)
+{
+  mp_limb_t *tmp; int carry; unsigned long nw; 
+
+  nw = prec / BITS_PER_MP_LIMB; 
+  if (prec & (BITS_PER_MP_LIMB - 1)) nw++;
+  tmp = (mp_ptr) (*_mp_allocate_func) (nw * BYTES_PER_MP_LIMB);
+  carry = mpfr_round_raw(tmp, x->_mp_d, RND_MODE, x->_mp_size, prec); 
+
+  if (carry)
+    {      
+      mpn_rshift(tmp, tmp, nw, 1); 
+      tmp [nw-1] |= (1UL << (BITS_PER_MP_LIMB - 1)); 
+      EXP(x)++; 
+    }
+
+  (*_mp_free_func) (x->_mp_d, ABSSIZE(x) * BYTES_PER_MP_LIMB);
+  if (SIGN(x)<0) { SIZE(x) = nw; CHANGE_SIGN(x); }
+  else SIZE(x) = nw;
+  PREC(x) = prec; 
+  MANT(x) = tmp; 
+}
+
+/* hypotheses : BITS_PER_MP_LIMB est une puissance de 2 
+                dans un test, la premiere partie du && est evaluee d'abord */
+
+
+/* assuming b is an approximation of x in direction rnd1 
+   with error at most 2^(EXP(b)-err), returns 1 if one is 
+   able to round exactly x to precision prec with direction rnd2,
+   and 0 otherwise.
+
+   Side effects: none.
+*/
+int mpfr_can_round(b, err, rnd1, rnd2, prec) mpfr_t b; unsigned int err;
+unsigned char rnd1, rnd2; unsigned long prec;
+{
+  int k, k1, l, l1, bn; mp_limb_t cc, cc2; char neg;
+
+  if (err<=prec) return 0;
+  bn = 1+(PREC(b)-1)/mp_bits_per_limb; /* number of sign. limbs of b */
+  neg = (SIGN(b)<0);
+  k = err/mp_bits_per_limb;
+  l = err % mp_bits_per_limb; if (l) l = mp_bits_per_limb-l;
+  /* the error corresponds to bit l in limb k */
+  k1 = prec/mp_bits_per_limb;
+  l1 = prec%mp_bits_per_limb; if (l1) l1 = mp_bits_per_limb-l1;
+  /* the last significant bit is bit l1 in limb k1 */
+  if (rnd1==GMP_RNDU && neg) rnd1=GMP_RNDZ;
+  switch (rnd1) {
+  case GMP_RNDZ: /* b <= x <= b+2^(EXP(b)-err) */
+    /* first round b */
+    cc = (MANT(b)[bn-k1-1]>>l1) & 1;
+    cc += mpfr_round_raw2(MANT(b), bn, neg, rnd2, prec);
+    /* now round b+2^(EXP(b)-err) */
+    mpn_add_1(MANT(b)+bn-k-1, MANT(b)+bn-k-1, k+1, (mp_limb_t)1<<l);
+    cc2 = (MANT(b)[bn-k1-1]>>l1) & 1;
+    /* as mpfr_round_raw2 returns a nonnegative value, if cc2>cc
+       then we already know we can't round */
+    if (cc2<=cc) cc2+=mpfr_round_raw2(MANT(b), bn, neg, rnd2, prec);
+    /* if parity of cc and cc2 equals, then one is able to round */
+    /* reset b to original value */
+    mpn_sub_1(MANT(b)+bn-k-1, MANT(b)+bn-k-1, k+1, (mp_limb_t)1<<l);
+    return (cc == cc2);
+  case GMP_RNDU: /* b-2^(EXP(b)-err) <= x <= b */
+    /* first round b */
+    cc = (MANT(b)[bn-k1-1]>>l1) & 1;
+    cc += mpfr_round_raw2(MANT(b), bn, neg, rnd2, prec);
+    /* now round b-2^(EXP(b)-err) */
+    cc2 = mpn_sub_1(MANT(b)+bn-k-1, MANT(b)+bn-k-1, k+1, (mp_limb_t)1<<l);
+    if (cc2) return 0;
+    cc2 = (MANT(b)[bn-k1-1]>>l1) & 1;
+    cc2 += mpfr_round_raw2(MANT(b), bn, neg, rnd2, prec);
+    /* if parity of cc and cc2 equals, then one is able to round */
+    /* reset b to original value */
+    mpn_add_1(MANT(b)+bn-k-1, MANT(b)+bn-k-1, k+1, (mp_limb_t)1<<l);
+    return (cc == cc2);
+  default:
+    printf("rnd1=%d not yet implemented in mpfr_round2\n",rnd1);
+    exit(1);
+  }
+}
diff --git a/set.c b/set.c
new file mode 100644
index 000000000..aa99ac9c8
--- /dev/null
+++ b/set.c
@@ -0,0 +1,11 @@
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "mpfr.h"
+
+void mpfr_set(a, b, rnd_mode) 
+mpfr_ptr a; mpfr_srcptr b; unsigned char rnd_mode;
+{
+  mpfr_round_raw(MANT(a), b->_mp_d, rnd_mode, b->_mp_size, PREC(a));
+  EXP(a) = EXP(b);
+  if (SIGN(a) != SIGN(b)) CHANGE_SIGN(a);
+}
diff --git a/set_d.c b/set_d.c
new file mode 100644
index 000000000..64bc05aeb
--- /dev/null
+++ b/set_d.c
@@ -0,0 +1,276 @@
+#include <math.h> /* for isnan */
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "longlong.h"
+#include "mpfr.h"
+
+#define NaN sqrt(-1) /* ensures a machine-independent NaN */
+
+/* Included from gmp-2.0.2, patched to support denorms */
+
+#ifdef XDEBUG
+#undef _GMP_IEEE_FLOATS
+#endif
+
+#ifndef _GMP_IEEE_FLOATS
+#define _GMP_IEEE_FLOATS 0
+#endif
+
+#define MP_BASE_AS_DOUBLE (2.0 * ((mp_limb_t) 1 << (BITS_PER_MP_LIMB - 1)))
+int
+#if __STDC__
+__mpfr_extract_double (mp_ptr rp, double d, int e)
+#else
+__mpfr_extract_double (rp, d)
+     mp_ptr rp;
+     double d;
+     int e;
+#endif
+     /* e=0 iff rp has only one limb */
+{
+  long exp;
+  mp_limb_t manh, manl;
+
+  /* BUGS
+
+     1. Should handle Inf and NaN in IEEE specific code.
+     2. Handle Inf and NaN also in default code, to avoid hangs.
+     3. Generalize to handle all BITS_PER_MP_LIMB >= 32.
+     4. This lits is incomplete and misspelled.
+   */
+
+  if (d == 0.0)
+    {
+      rp[0] = 0;
+#if BITS_PER_MP_LIMB == 32
+      if (e) rp[1] = 0;
+#endif
+      return 0;
+    }
+
+#if _GMP_IEEE_FLOATS
+  {
+    union ieee_double_extract x;
+    x.d = d;
+
+    exp = x.s.exp;
+    if (exp) 
+      {
+#if BITS_PER_MP_LIMB == 64
+	manl = (((mp_limb_t) 1 << 63)
+		| ((mp_limb_t) x.s.manh << 43) | ((mp_limb_t) x.s.manl << 11));
+#else
+	manh = ((mp_limb_t) 1 << 31) | (x.s.manh << 11) | (x.s.manl >> 21);
+	manl = x.s.manl << 11;      
+#endif
+      }
+    else
+      {
+#if BITS_PER_MP_LIMB == 64
+	manl = ((mp_limb_t) x.s.manh << 43) | ((mp_limb_t) x.s.manl << 11);
+#else
+    manh = (x.s.manh << 11) | (x.s.manl >> 21);
+	manl = x.s.manl << 11;      
+#endif
+      }
+  }
+#else
+  {
+    /* Unknown (or known to be non-IEEE) double format.  */
+    exp = 0;
+    if (d >= 1.0)
+      {
+        if (d * 0.5 == d)
+          abort ();
+
+        while (d >= 32768.0)
+          {
+            d *= (1.0 / 65536.0);
+            exp += 16;
+          }
+        while (d >= 1.0)
+          {
+            d *= 0.5;
+            exp += 1;
+          }
+      }
+    else if (d < 0.5)
+      {
+        while (d < (1.0 / 65536.0))
+          {
+            d *=  65536.0;
+            exp -= 16;
+          }
+        while (d < 0.5)
+          {
+            d *= 2.0;
+            exp -= 1;
+          }
+      }
+
+    d *= MP_BASE_AS_DOUBLE;
+#if BITS_PER_MP_LIMB == 64
+    manl = d;
+#else
+    manh = d;
+    manl = (d - manh) * MP_BASE_AS_DOUBLE;
+#endif
+
+    exp += 1022;
+  }
+#endif
+
+  if (exp) exp = (unsigned) exp - 1022; else exp = -1021; 
+
+#if BITS_PER_MP_LIMB == 64
+      rp[0] = manl;
+#else
+      if (e) {
+	rp[1] = manh;
+	rp[0] = manl;
+      }
+      else {
+	rp[0] = manh;
+      }
+#endif
+
+  return exp;
+}
+
+/* End of part included from gmp-2.0.2 */
+/* Part included from gmp temporary releases */
+double
+#if __STDC__
+__mpfr_scale2 (double d, int exp)
+#else
+__mpfr_scale2 (d, exp)
+     double d;
+     int exp;
+#endif
+{
+#if _GMP_IEEE_FLOATS
+  {
+    union ieee_double_extract x;
+    x.d = d;
+    exp += x.s.exp;
+    x.s.exp = exp;
+    if (exp >= 2047)
+      {
+        /* Return +-infinity */
+        x.s.exp = 2047;
+        x.s.manl = x.s.manh = 0;
+      }
+    else if (exp < 1)
+      {
+        x.s.exp = 1;            /* smallest exponent (biased) */
+        /* Divide result by 2 until we have scaled it to the right IEEE
+           denormalized number, but stop if it becomes zero.  */
+        while (exp < 1 && x.d != 0)
+          {
+            x.d *= 0.5;
+            exp++;
+          }
+      }
+    return x.d;
+  }
+#else
+  {
+    double factor, r;
+
+    factor = 2.0;
+    if (exp < 0)
+      {
+        factor = 0.5;
+        exp = -exp;
+      }
+    r = d;
+    if (exp != 0)
+      {
+        if ((exp & 1) != 0)
+          r *= factor;
+        exp >>= 1;
+        while (exp != 0)
+          {
+            factor *= factor;
+            if ((exp & 1) != 0)
+              r *= factor;
+            exp >>= 1;
+          }
+      }
+    return r;
+  }
+#endif
+}
+
+
+/* End of part included from gmp */
+
+void
+mpfr_set_d(mpfr_t r, double d, unsigned char rnd_mode)
+{
+  int negative, sizer; unsigned int cnt;
+
+  if (d == 0)
+    {      
+      EXP(r) = 0;
+      return;
+    }
+  else if (isnan(d)) { SET_NAN(r); return; }
+
+  negative = d < 0;
+  d = ABS (d);
+
+  sizer = MPFR_LIMBS_PER_DOUBLE; if (ABSSIZE(r)<sizer) sizer=ABSSIZE(r);
+  /* warning: __mpfr_extract_double requires at least two limbs */
+  EXP(r) = __mpfr_extract_double (MANT(r), d, (sizer>=2) );
+  
+  count_leading_zeros(cnt, MANT(r)[sizer-1]);
+  if (cnt) mpn_lshift(MANT(r), MANT(r), sizer, cnt); 
+  
+  EXP(r) -= cnt; 
+  SIZE(r) = sizer; if (negative) CHANGE_SIGN(r);
+
+  mpfr_round(r, rnd_mode, PREC(r)); 
+  return; 
+}
+
+
+double
+mpfr_get_d(mpfr_t src)
+{
+  double res;
+  mp_size_t size, i, n_limbs_to_use;
+  mp_ptr qp;
+  int negative;
+
+  if (FLAG_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);
+  negative = (SIGN(src)==-1);
+
+  /* Warning: don't compute the abs(res) and set the sign afterwards,
+     otherwise the current machine rounding mode will not be taken
+     correctly into account. */
+  /* res = (negative) ? -(double)qp[size - 1] : qp[size - 1]; */
+  res = 0.0;
+  /* Warning: an arbitrary number of limbs may be required for an exact 
+     rounding. The following code is correct but not optimal since one
+     may be able to decide without considering all limbs. */
+  /* n_limbs_to_use = MIN (MPFR_LIMBS_PER_DOUBLE, size); */
+  n_limbs_to_use = size;
+  /* Accumulate the limbs from less significant to most significant
+     otherwise due to rounding we may accumulate several ulps,
+     especially in rounding towards -/+infinity. */
+  for (i = n_limbs_to_use; i>=1; i--)
+    res = res / MP_BASE_AS_DOUBLE +
+      ((negative) ? -(double)qp[size - i] : qp[size - i]);
+  res = __mpfr_scale2 (res, EXP(src) - BITS_PER_MP_LIMB); 
+
+  return res;
+}
+
diff --git a/set_dfl_prec.c b/set_dfl_prec.c
new file mode 100644
index 000000000..aa7b7b668
--- /dev/null
+++ b/set_dfl_prec.c
@@ -0,0 +1,16 @@
+#include "gmp.h"
+#include "gmp-impl.h"
+
+/* default is 53 bits */
+mp_size_t __gmp_default_fp_bit_precision = 53;
+
+void
+#if __STDC__
+mpfr_set_default_prec (unsigned long int prec_in_bits)
+#else
+mpfr_set_default_prec (prec_in_bits)
+     unsigned long int prec_in_bits;
+#endif
+{
+  __gmp_default_fp_bit_precision = prec_in_bits;
+}
diff --git a/set_dfl_rnd.c b/set_dfl_rnd.c
new file mode 100644
index 000000000..efef4d388
--- /dev/null
+++ b/set_dfl_rnd.c
@@ -0,0 +1,16 @@
+#include "gmp.h"
+#include "gmp-impl.h"
+
+char __gmp_default_rounding_mode = 0;
+
+void
+#if __STDC__
+mpfr_set_default_rounding_mode (char rnd_mode)
+#else
+mpfr_set_default_rounding_mode (rnd_mode)
+     char rnd_mode;
+#endif
+{
+  __gmp_default_rounding_mode = rnd_mode;
+}
+
diff --git a/set_f.c b/set_f.c
new file mode 100644
index 000000000..e352d6a35
--- /dev/null
+++ b/set_f.c
@@ -0,0 +1,36 @@
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "longlong.h"
+#include "mpfr.h"
+
+void mpfr_set_f(mpfr_ptr y, mpf_srcptr x, unsigned long prec, char rnd_mode)
+{
+  mp_limb_t *my, *mx; unsigned long cnt, sx, sy; 
+
+  sx = SIZ (x); sy = ABSSIZE(y);
+  my = MANT(y); mx = MANT(x); 
+
+  count_leading_zeros(cnt, mx[sx - 1]);  
+
+  if (sy < sx)
+    {
+      if (cnt) 
+	{ 
+	  mpn_lshift(my, mx + 1, sy, cnt); 
+	  my [0] |= mx[0] >> (BITS_PER_MP_LIMB - cnt); 
+	  EXP(y) = EXP(x)* BITS_PER_MP_LIMB - cnt; 
+	}
+      else { MPN_COPY(my, mx, sy); EXP(y) = EXP(x) * BITS_PER_MP_LIMB; }
+    }
+  else
+    {
+      if (cnt) 
+	{ 
+	  mpn_lshift(my, mx, sx, cnt); 
+	  EXP(y) = EXP(x)* BITS_PER_MP_LIMB - cnt; 
+	}
+      else { MPN_COPY(my, mx, sy); EXP(y) = EXP(x) * BITS_PER_MP_LIMB; }
+    }
+  
+  mpfr_round(y, rnd_mode, prec); 
+}
diff --git a/set_prec.c b/set_prec.c
new file mode 100644
index 000000000..a044a5344
--- /dev/null
+++ b/set_prec.c
@@ -0,0 +1,46 @@
+#include <stdio.h>
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "mpfr.h"
+
+void
+#if __STDC__
+mpfr_set_prec (mpfr_t x, unsigned long int p, unsigned char rnd_mode)
+#else
+mpfr_set_prec (x, p, rnd_mode)
+     mpfr_t x;
+     unsigned long int p;
+     unsigned char rnd_mode;
+#endif
+{
+  unsigned long xsize,oldp,oldsize; mp_limb_t *old;
+
+  if (p==0) {
+    printf("*** cannot set precision to 0 bits\n"); exit(1);
+  }
+
+  oldp = x -> _mp_prec;
+  oldsize = (oldp-1)/BITS_PER_MP_LIMB + 1;
+  if (SIGN(x)<0) oldsize = oldsize ^ (1<<31);
+  xsize = (p - 1)/BITS_PER_MP_LIMB + 1; /* new limb size */
+
+  old = x -> _mp_d; 
+  x -> _mp_d = (mp_ptr) (*_mp_allocate_func) 
+    (xsize * BYTES_PER_MP_LIMB);
+  x -> _mp_prec = p;
+  mpfr_round_raw(x -> _mp_d, old, rnd_mode, oldsize, p);
+  SIZE(x) = (SIGN(x)>0) ? xsize : (xsize ^ (1<<31));
+
+  (*_mp_free_func) (old, 1 + ((oldp-1)>>3));
+}
+
+unsigned long int
+#if __STDC__
+mpfr_get_prec (mpfr_t x)
+#else
+mpfr_set_prec (x)
+     mpfr_t x;
+#endif
+{
+  return x -> _mp_prec;
+}
diff --git a/set_si.c b/set_si.c
new file mode 100644
index 000000000..3dccefc10
--- /dev/null
+++ b/set_si.c
@@ -0,0 +1,43 @@
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "longlong.h"
+#include "mpfr.h"
+
+void
+mpfr_set_si(mpfr_ptr x, long i, unsigned char rnd_mode)
+{
+  unsigned long xsize, ai, cnt; 
+
+  xsize = ABSSIZE(x);
+  ai = ABS(i); 
+
+  count_leading_zeros(cnt, ai); 
+
+  x -> _mp_d[xsize - 1] = ai << cnt;
+  /* don't forget to put zero in lower limbs */
+  MPN_ZERO(MANT(x), xsize-1);
+  x -> _mp_exp = BITS_PER_MP_LIMB - cnt;
+  /* warning: don't change the precision of x! */
+  x -> _mp_size = xsize; if (i < 0) CHANGE_SIGN(x);
+
+  return; 
+}
+
+void
+mpfr_set_ui(mpfr_ptr x, unsigned long i, unsigned char rnd_mode)
+{
+  unsigned int xsize, cnt; 
+
+  xsize = ABSSIZE(x);
+  count_leading_zeros(cnt, i); 
+
+  x -> _mp_d[xsize - 1] = i << cnt; 
+  /* don't forget to put zero in lower limbs */
+  MPN_ZERO(MANT(x), xsize-1);
+  x -> _mp_exp = BITS_PER_MP_LIMB - cnt;
+  /* warning: don't change the precision of x! */
+  x -> _mp_size = xsize; 
+
+  return; 
+}
+
diff --git a/set_str.c b/set_str.c
new file mode 100644
index 000000000..8b1378917
--- /dev/null
+++ b/set_str.c
@@ -0,0 +1 @@
+
diff --git a/set_str_raw.c b/set_str_raw.c
new file mode 100644
index 000000000..bef006126
--- /dev/null
+++ b/set_str_raw.c
@@ -0,0 +1,82 @@
+#include <stdio.h>
+#include <stdlib.h>
+#ifdef HAS_STRING_H
+#include <string.h>
+#else
+#include <strings.h>
+#endif
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "longlong.h"
+#include "mpfr.h"
+
+/* Currently the number should be of the form +/- xxxx.xxxxxxEyy, with 
+   decimal exponent. The mantissa of x is supposed to be large enough 
+   to hold all the bits of str. */
+
+void
+mpfr_set_str_raw(mpfr_ptr x, char *str)
+{
+  char *str2, *str0, negative = 0; 
+  unsigned long j, l, k = 0, xsize, cnt; mp_limb_t *xp; 
+  long expn = 0; char *endstr2;
+
+  xp = x -> _mp_d; 
+  xsize = 1 + (PREC(x)-1)/BITS_PER_MP_LIMB;
+  str0 = str2 = (char *) malloc(strlen(str)*sizeof(char)); 
+
+  if (*str == '-') { negative = 1; str++; }
+
+  while (*str == '0') { str++; } 
+
+  while (*str == '0' || *str == '1')
+    { *(str2++) = *(str++); k++; }
+
+  if (*str == '.') 
+    {
+      str++;
+      while (*str == '0' || *str == '1')
+	{ *(str2++) = *(str++); }
+
+      if (*str == '[') { while (*str != ']') str++; }
+    }
+    
+  if (*str == '[') { while (*str != ']') str++; }
+
+  if (*str == 'e' || *str == 'E') 
+    {
+      expn = k + atoi(++str); 
+      if (expn < k)
+	{
+	  fprintf(stderr, "Warning : possible overflow in exponent in Str -> mpfr\n"); 
+	}
+    }
+  endstr2 = str2;
+  l = (strlen(str0) - 1) / BITS_PER_MP_LIMB + 1; str2 = str0; 
+
+  /* str2[0]..endstr2[-1] contains the mantissa */
+  for (k = 1; k <= l; k++)
+    {
+      j = 0; 
+      xp[xsize - k] = 0; 
+      while (str2<endstr2 && j < BITS_PER_MP_LIMB)
+	{
+	  xp[xsize - k] = (xp[xsize - k] << 1) + (*str2 - '0'); 
+	  str2++; j++; 
+	}
+      xp[xsize - k] <<= (BITS_PER_MP_LIMB - j); 
+    }
+
+  for (; k <= xsize; k++) { xp[xsize - k] = 0; }
+
+  count_leading_zeros(cnt, xp[xsize - 1]); 
+  if (cnt) mpn_lshift(xp, xp, xsize, cnt); 
+
+  x -> _mp_exp = expn - cnt; 
+  x -> _mp_size = (negative ? -xsize : xsize); 
+
+  free(str0); 
+  
+  /* May change to take into account : syntax errors, overflow in exponent, 
+     string truncated because of size of x */
+}
diff --git a/sqrt.c b/sqrt.c
new file mode 100644
index 000000000..c6f10900c
--- /dev/null
+++ b/sqrt.c
@@ -0,0 +1,73 @@
+#include <stdio.h>
+#include <math.h>
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "mpfr.h"
+
+/* #define DEBUG */
+
+void print_double(x) mpfr_t x;
+{
+  int e, i; double d;
+
+  d = mpfr_get_d(x);
+  e = (int) ceil(log(d)/log(2.0));
+  /* d <= 2^e */
+  e -= 53;
+  if (e>0) for (i=0;i<e;i++) d /= 2.0;
+  else for (i=0;i<-e;i++) d *= 2.0;
+  printf("%1.0f*2^(%d)",d,e);
+}
+
+void mpfr_sqrt(mpfr_ptr x, mpfr_srcptr a, unsigned char rnd_mode)
+{
+  int p, q, err, i, e, n; mpfr_t t, u;
+
+#ifdef DEBUG
+  printf("enter mpfr_sqrt, a=%1.20e, rnd=%d\n",mpfr_get_d(a), rnd_mode);
+  printf("a="); print_double(a); putchar('\n');
+#endif
+  /* use Newton's iteration x[n+1] = 1/2*(x[n]+a/x[n]),
+     the error e[n] = x[n]-sqrt(a) satisfies e[n+1] <= e[n]/2/sqrt(a) */
+  if (SIGN(a)<0) { SET_NAN(x); return; }
+  e = EXP(a)/2; if (2*e<EXP(a)) e++;
+#ifdef DEBUG
+  printf("e=%d\n",e);
+#endif
+  /* now 2^(2*e-2) <= a <= 2^(2*e) i.e. 1/4 <= a/2^(2e) <= 1 */
+  q = p = PREC(x);
+  for (i=0; i<3; i++)
+    q = p + (int) ceil(log(4.0*ceil(log((double)q)/log(2.0))+2.0)/log(2.0));
+  err = q-p; /* the error is at most 2^err ulp */
+  q = ((q-1)/mp_bits_per_limb)*mp_bits_per_limb; /* adjust to entire limb */
+  mpfr_init(t); mpfr_init(u);
+  do {
+    q += mp_bits_per_limb;
+#ifdef DEBUG
+    printf("prec=%d q=%d err=%d\n",p,q,err);
+#endif
+    mpfr_set_prec(t, q, GMP_RNDU); 
+    mpfr_set_prec(x, q, GMP_RNDU); 
+    mpfr_set_prec(u, q, GMP_RNDU);
+    mpfr_set_ui(x, 1, GMP_RNDU);
+    EXP(x) += e;
+    n = (int) ceil(log((double) q)/log(2.0));
+    for (i=0; i<n; i++) {
+      mpfr_div(t, a, x, GMP_RNDU);
+      mpfr_add(u, x, t, GMP_RNDU);
+      mpfr_div_2exp(x, u, 1, GMP_RNDU);
+#ifdef DEBUG
+      printf("t="); print_double(t); putchar('\n');
+      printf("u="); print_double(u); putchar('\n');
+      printf("x="); print_double(x); putchar('\n');
+      printf("i=%d t=%1.20e u=%1.20e x=%1.20e\n",i,mpfr_get_d(t),mpfr_get_d(u),
+	     mpfr_get_d(x));
+      printf("t="); mpfr_print_raw(t); putchar('\n');
+      printf("u="); mpfr_print_raw(u); putchar('\n');
+      printf("x="); mpfr_print_raw(x); putchar('\n');
+#endif
+    }
+  } while (mpfr_can_round(x, q-err, GMP_RNDU, rnd_mode, p)==-1);
+  mpfr_round(x, rnd_mode, p);
+  mpfr_clear(t); mpfr_clear(u);
+}
diff --git a/sub.c b/sub.c
new file mode 100644
index 000000000..46e52427a
--- /dev/null
+++ b/sub.c
@@ -0,0 +1,412 @@
+#include <stdio.h>
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "mpfr.h"
+
+/* #define DEBUG2 */
+
+#ifdef DEBUG2
+mp_limb_t *ap0, *ap1;
+#define ck(x, dx, i) /* mp_limb_t *x; int dx; */ \
+{ \
+  if ((x)<ap0 || (x)+(dx)-1>ap1) { \
+    printf("error %d\n",(i)); \
+    printf("%1.20e,%d,%1.20e,%d,%d,%d\n",mpfr_get_d(b),PREC(b),mpfr_get_d(c), \
+	   PREC(c),PREC(a),rnd_mode); \
+    exit(1); } \
+}
+#else
+#define ck(x, dx, i) {}
+#endif
+
+void mpfr_sub(a, b, c, rnd_mode) 
+mpfr_ptr a; mpfr_srcptr b, c; unsigned char rnd_mode;
+{
+  CHANGE_SIGN(c);
+  mpfr_add(a, b, c, rnd_mode);
+  CHANGE_SIGN(c);
+}
+
+/* put in ap[0]..ap[an-1] the value of bp[0]..bp[n-1] shifted by sh bits
+   to the left minus ap[0]..ap[n-1], with 0 <= sh < mp_bits_per_limb, and
+   returns the borrow.
+*/
+mp_limb_t mpn_sub_lshift_n (ap, bp, n, sh, an) mp_limb_t *ap, *bp; int n,sh,an;
+{
+  mp_limb_t c, bh;
+
+  /* shift b to the left */
+  if (sh) bh = mpn_lshift(bp, bp, n, sh);
+  c = (n<an) ? mpn_sub_n(ap, bp, ap, n) : mpn_sub_n(ap, bp+(n-an), ap, an);
+  /* shift back b to the right */
+  if (sh) {
+    mpn_rshift(bp, bp, n, sh);
+    bp[n-1] += bh<<(mp_bits_per_limb-sh);
+  }
+  return c;
+}
+
+/* signs of b and c differ, abs(b)>=abs(c), diff_exp>=0 */
+void mpfr_sub1(a, b, c, rnd_mode, diff_exp) 
+mpfr_t a, b, c; unsigned char rnd_mode; int diff_exp;
+{
+  mp_limb_t *ap, *bp, *cp, cc, c2; unsigned int an,bn,cn; 
+  int sh,dif,k,cancel,cancel1,cancel2;
+
+#ifdef DEBUG2
+  printf("b=  "); if (SIGN(b)>=0) putchar(' ');
+  mpfr_print_raw(b); putchar('\n');
+  printf("c=  "); if (SIGN(c)>=0) putchar(' ');
+  for (k=0;k<diff_exp;k++) putchar(' '); mpfr_print_raw(c);
+  putchar('\n');
+  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 */
+  cancel1 = cancel/mp_bits_per_limb; cancel2 = cancel%mp_bits_per_limb;
+  ap = MANT(a);
+  bp = MANT(b);
+  cp = MANT(c);
+  an = (PREC(a)-1)/mp_bits_per_limb+1; /* number of significant limbs of a */
+#ifdef DEBUG2
+ap0=ap; ap1=ap+an-1;
+#endif
+  sh = an*mp_bits_per_limb-PREC(a); /* non-significant bits in low limb */
+  bn = (PREC(b)-1)/mp_bits_per_limb+1; /* number of significant limbs of b */
+  cn = (PREC(c)-1)/mp_bits_per_limb + 1;
+  EXP(a) = EXP(b)-cancel;
+  /* adjust sign to that of b */
+  if (SIGN(a)*SIGN(b)<0) CHANGE_SIGN(a);
+  /* case 1: diff_exp>=prec(a), i.e. c only affects the last bit
+     through rounding */
+  dif = PREC(a)-diff_exp;
+#ifdef DEBUG2
+  printf("PREC(a)=%d an=%u PREC(b)=%d bn=%u PREC(c)=%d diff_exp=%u dif=%d\n",
+	 PREC(a),an,PREC(b),bn,PREC(c),diff_exp,dif);
+#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 (cancel2) { ck(ap+(an-bn+cancel1), bn-cancel1, 1);
+	mpn_lshift(ap+(an-bn+cancel1), bp, bn-cancel1, cancel2); }
+      else { ck(ap+(an-bn+cancel1), bn-cancel1, 2);
+	MPN_COPY(ap+(an-bn+cancel1), bp, bn-cancel1); }
+      /* fill low significant limbs with zero */
+ck(ap, an-bn+cancel1, 3);
+      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 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 */
+	    /* check whether mant(c)=1/2 or not */
+	    cc = *cp - (1<<(mp_bits_per_limb-1));
+	    if (cc==0) {
+	      bp = cp+(PREC(c)-1)/mp_bits_per_limb;
+	      while (cp<bp && cc==0) cc = *++cp;
+	    }
+ck(ap+an-1, 1, 4);
+	    if (cc || (ap[an-1] & 1<<sh)) goto sub_one_ulp;
+	      /* mant(c) > 1/2 or mant(c) = 1/2: subtract 1 iff lsb(a)=1 */
+	  }
+	}
+      }
+      else if ((ISNONNEG(b) && rnd_mode==GMP_RNDU) || 
+	       (ISNEG(b) && rnd_mode==GMP_RNDD)) {
+	/* round up: nothing to do */
+      }
+      else {
+	/* round down: subtract 1 ulp iff c<>0 */
+	if (NOTZERO(c)) goto sub_one_ulp;
+      }
+    }
+    else { /* PREC(b)>PREC(a) : we have to round b-c */
+      k=bn-an;
+      /* first copy the 'an' most significant limbs of b to a */
+ck(ap, an, 5);
+      MPN_COPY(ap, bp+k, an);
+      if (rnd_mode==GMP_RNDN) { /* to nearest */
+	/* first check whether the truncated bits from b are 1/2*lsb(a) */
+	if (sh) {
+ck(ap, 1, 6);
+	  cc = *ap & (((mp_limb_t)1<<sh)-1);
+	  *ap &= ~cc; /* truncate last bits */
+	  cc -= (mp_limb_t)1<<(sh-1);
+	}
+	else /* no bit to truncate */
+	  cc = bp[--k] - (1<<(mp_bits_per_limb-1));
+	if ((long)cc>0) { /* suppose sizeof(long)=sizeof(mp_limb_t) */
+	  goto add_one_ulp; /* trunc(b)>1/2*lsb(a) -> round up */
+	}
+	else if (cc==0) {
+	  while (k>1 && cc==0) cc=bp[--k];
+	  /* now if the truncated part of b = 1/2*lsb(a), check whether c=0 */
+ck(ap, 1, 7);
+	  if (NOTZERO(c) || (*ap & (1<<sh))) goto sub_one_ulp;
+	  /* if trunc(b)-c is exactly 1/2*lsb(a) : round to even lsb */
+	}
+	/* if cc<0 : trunc(b) < 1/2*lsb(a) -> round down, i.e. do nothing */
+      }
+      else { /* round towards infinity or zero */
+	if (sh) {
+ck(ap, 1, 8);
+	  cc = *ap & ((1<<sh)-1);
+	  *ap &= ~cc; /* truncate last bits */
+	}
+	else cc=0;
+	cn--;
+        c2 = (dif>-sh) ? cp[cn]>>(mp_bits_per_limb-dif-sh) : 0;
+	while (cc==c2 && (k || cn)) {
+	  if (k) cc = bp[--k];
+	  if (cn) {
+	    c2 = cp[cn]<<(dif+sh);
+	    if (--cn) c2 += cp[cn]>>(mp_bits_per_limb-dif-sh);
+	  }
+	}
+	dif = ((ISNONNEG(b) && rnd_mode==GMP_RNDU) || 
+	       (ISNEG(b) && rnd_mode==GMP_RNDD));
+	/* round towards infinity if dif=1, towards zero otherwise */
+	if (dif && cc>c2) goto add_one_ulp;
+	else if (dif==0 && cc<c2) goto sub_one_ulp;
+      }
+    }
+  }
+  else { /* case 2: diff_exp < PREC(a) : c overlaps with a by dif bits */
+    /* first copy upper part of c into a (after shift) */
+    int overlap;
+    dif += cancel;
+    k = (dif-1)/mp_bits_per_limb + 1; /* only the highest k limbs from c
+					 have to be considered */
+    if (k<an) {
+ck(ap+k, an-k, 9);
+      MPN_ZERO(ap+k, an-k); /* do it now otherwise ap[k] may be 
+				       destroyed in case dif<0 */
+    }
+#ifdef DEBUG2
+    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 */
+      dif = dif % mp_bits_per_limb;
+      dif = (dif) ? mp_bits_per_limb-dif-sh : -sh;
+      /* we have to shift by dif bits to the right */
+      if (dif>0) {
+ck(ap, (k<=an) ? k : an, 10);
+	mpn_rshift(ap, cp+(cn-k), (k<=an) ? k : an, dif);
+        if (k>an) ap[an-1] += cp[cn-k+an]<<(mp_bits_per_limb-dif);
+      }
+      else if (dif<0) {
+ck(ap, k, 11);
+	cc = mpn_lshift(ap, cp+(cn-k), k, -dif);
+	if (k<an) { ck(ap+k, 1, 12); ap[k]=cc; }
+	/* put the non-significant bits in low limb for further rounding */
+	ap[0] += cp[cn-k-1]>>(mp_bits_per_limb+dif);
+      }
+      else { ck(ap, k, 13); MPN_COPY(ap, cp+(cn-k), k); }
+      overlap=1;
+    }
+    else { /* c is not truncated, but we have to fill low limbs with 0 */
+ck(ap, k-cn, 14);
+      MPN_ZERO(ap, k-cn);
+      overlap = cancel-diff_exp;
+#ifdef DEBUG2
+      printf("0:a="); mpfr_print_raw(a); putchar('\n');
+      printf("overlap=%d\n",overlap);
+#endif
+      if (overlap>=0) {
+	cn -= overlap/mp_bits_per_limb;
+	/* warning: a shift of zero with mpn_lshift is not allowed */
+	if (overlap) {
+	  if (an<cn) {
+ck(ap, an, 15);
+	    mpn_lshift(ap, cp+(cn-an), an, overlap%mp_bits_per_limb);
+	    ap[0] += cp[cn-an-1]>>(mp_bits_per_limb-overlap);
+	  }
+	  else {
+ck(ap+(an-cn), cn, 15);
+	    mpn_lshift(ap+(an-cn), cp, cn, overlap%mp_bits_per_limb);
+	  }
+	}
+	else { ck(ap+(an-cn), cn, 15); MPN_COPY(ap+(an-cn), cp, cn); }
+      }
+      else { /* shift to the right by -overlap bits */
+	overlap = -overlap;
+	k = overlap/mp_bits_per_limb;
+	cc=mpn_rshift(ap+(an-k-cn),cp,cn,overlap%mp_bits_per_limb);
+	if (an>k+cn) ap[an-k-cn-1]=cc;
+      }
+      overlap=0;
+    }
+#ifdef DEBUG2
+      printf("1:a="); mpfr_print_raw(a); putchar('\n');
+#endif
+    /* 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;
+      overlap += 2;
+      /* invert low limbs */
+      imax = (int)an-(int)bn+cancel1;
+      for (i=0;i<imax;i++) { ck(ap+i, 1, 17); ap[i] = ~ap[i]; }
+if (i) ck(ap, i, 18);
+      cc = (i) ? mpn_add_1(ap, ap, i, 1) : 1;
+ck(ap+i, ((bn-cancel1)<an) ? (bn-cancel1) : an, 19);
+      mpn_sub_lshift_n(ap+i, bp, bn-cancel1, cancel2, an);
+ck(ap+i, an-i, 20);
+      mpn_sub_1(ap+i, ap+i, an-i, 1-cc);
+    }
+    else { /* PREC(b) > PREC(a): we have to truncate b */
+ck(ap, an, 21);
+      mpn_sub_lshift_n(ap, bp+(bn-an-cancel1), an, cancel2, an);
+    }
+    /* remains to do the rounding */
+#ifdef DEBUG2
+      printf("2:a="); mpfr_print_raw(a); putchar('\n');
+      printf("overlap=%d\n",overlap);
+#endif
+    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) */
+      switch (overlap)
+	{
+        case 1: /* both b and c to round */
+	  kc = cn-k; /* remains kc limbs from c */
+	  k = bn-an; /* remains k limbs from b */
+	  /* truncate last bits and store the difference with 1/2*ulp in cc */
+	  cc = *ap & ((1<<sh)-1);
+	  *ap &= ~cc; /* truncate last bits */
+	  cc -= 1<<(sh-1);
+	  while ((cc==0 || cc==-1) && k!=0 && kc!=0) {
+	    kc--;
+	    cc -= mpn_sub_1(&c2, bp+(--k), 1, (cp[kc]>>dif) +
+			    (cp[kc+1]<<(mp_bits_per_limb-dif)));
+	    if (cc==0 || cc==-1) cc=c2;
+	  }
+	  if ((long)cc>0) goto add_one_ulp;
+	  else if ((long)cc<-1) return; /* the carry can be at most 1 */
+	  else if (kc==0) goto round_b;
+	  /* else round c: go through */
+	case 3: /* only c to round */
+	  bp=cp; k=cn-k; goto to_nearest;
+	case 0: /* only b to round */
+        round_b:
+	  k=bn-an; dif=0; goto to_nearest;
+        /* otherwise the result is exact: nothing to do */
+	}
+    }
+    else if ((ISNONNEG(b) && rnd_mode==GMP_RNDU) || 
+	     (ISNEG(b) && rnd_mode==GMP_RNDD)) {
+      cc = *ap & ((1<<sh)-1);
+      *ap &= ~cc; /* truncate last bits */
+      if (cc) goto add_one_ulp; /* will happen most of the time */
+      else { /* same four cases too */
+	int kc = cn-k; /* remains kc limbs from c */
+	switch (overlap)
+	{
+        case 1: /* both b and c to round */
+	  k = bn-an; /* remains k limbs from b */
+          dif = diff_exp % mp_bits_per_limb;
+	  while (cc==0 && k!=0 && kc!=0) {
+	    kc--;
+	    cc = bp[--k] - (cp[kc]>>dif);
+	    if (dif) cc -= (cp[kc+1]<<(mp_bits_per_limb-dif));
+	  }
+	  if (cc) goto add_one_ulp;
+	  else if (kc==0) goto round_b2;
+	  /* else round c: go through */
+	case 3: /* only c to round: nothing to do */
+	  /* while (kc) if (cp[--kc]) goto add_one_ulp; */
+	  /* if dif>0 : remains to check last dif bits from c */
+	  /* if (dif>0 && (cp[0]<<(mp_bits_per_limb-dif))) goto add_one_ulp; */
+	  break;
+	case 0: /* only b to round */
+        round_b2:
+	  k=bn-an;
+	  while (k) if (bp[--k]) goto add_one_ulp;
+        /* otherwise the result is exact: nothing to do */
+	}
+      }
+    }
+    /* else round to zero: remove 1 ulp if neglected bits from b-c are < 0 */
+    else {
+      cc = *ap & ((1<<sh)-1); 
+      *ap &= ~cc;
+      if (cc==0) { /* otherwise neglected difference cannot be < 0 */
+	/* take into account bp[0]..bp[bn-cancel1-1] shifted by cancel2 to left
+	   and cp[0]..cp[cn-k-1] shifted by dif bits to right */
+	switch (overlap) {
+	case 0:
+	case 2:
+	  break; /* c is not truncated ==> no borrow */
+	case 1: /* both b and c are truncated */
+	  break;
+	case 3: /* only c is truncated */
+	  cn -= k; /* take into account cp[0]..cp[cn-1] shifted by dif bits
+		      to the right */
+	  cc = (dif>0) ? cp[cn]<<(mp_bits_per_limb-dif) : 0;
+	  while (cc==0 && cn>0) cc = cp[--cn];
+          if (cc) goto sub_one_ulp;
+	  break;
+	}
+      }
+    }
+  }
+  goto end_of_sub;
+    
+  to_nearest: /* 0 <= sh < mp_bits_per_limb : number of bits of a to truncate
+                 bp[k] : last significant limb from b */
+#ifdef DEBUG2
+mpfr_print_raw(a); putchar('\n');
+#endif
+        if (sh) {
+	  cc = *ap & ((1<<sh)-1);
+	  *ap &= ~cc; /* truncate last bits */
+	  c2 = 1<<(sh-1);
+	}
+	else /* no bit to truncate */
+	  { cc = bp[--k]; c2 = 1<<(mp_bits_per_limb-1); }
+#ifdef DEBUG2
+	printf("cc=%lu c2=%lu k=%u\n",cc,c2,k);
+#endif
+	if (cc>c2) goto add_one_ulp; /* trunc(b)>1/2*lsb(a) -> round up */
+	else if (cc==c2) {
+	  cc=0; while (k && cc==0) cc=bp[--k];
+#ifdef DEBUG2
+	  printf("cc=%lu\n",cc);
+#endif
+	  /* special case of rouding c shifted to the right */
+	  if (cc==0 && dif>0) cc=bp[0]<<(mp_bits_per_limb-dif);
+	  /* now if the truncated part of b = 1/2*lsb(a), check whether c=0 */
+          if (bp!=cp) { 
+	    if (cc || (*ap & (1<<sh))) goto add_one_ulp;
+	  }
+	  else {
+	    /* subtract: if cc>0, do nothing */
+	    if (cc==0 && (*ap & (1<<sh))) goto add_one_ulp;
+	  }
+	}
+        goto end_of_sub;
+
+ sub_one_ulp:
+    cc = mpn_sub_1(ap, ap, an, 1<<sh);
+    if (cc) { printf("carry(2) in mpfr_sub1\n"); exit(1); }
+  goto end_of_sub;
+
+  add_one_ulp: /* add one unit in last place to a */
+    cc = mpn_add_1(ap, ap, an, 1<<sh);
+    if (cc) { printf("carry(3) in mpfr_sub1\n"); exit(1); }
+
+ end_of_sub:
+#ifdef DEBUG2
+printf("b-c="); if (SIGN(a)>0) putchar(' '); mpfr_print_raw(a); putchar('\n');
+#endif
+  return;
+}
+
diff --git a/tests/Makefile b/tests/Makefile
new file mode 100644
index 000000000..4e4c837bd
--- /dev/null
+++ b/tests/Makefile
@@ -0,0 +1,71 @@
+CC=gcc
+INCFLAGS=-I$(GMPDIR)/include -I..
+LDFLAGS=libmpfr.a $(GMPDIR)/lib/libgmp.a 
+CFLAGS=-g -DMsb
+
+all_files=tmul tround tset_f tcmp tcmp2 tdiv tmul_2exp tmul_ui tset_d tset_i tset_str tadd tsqrt tagm
+
+all: $(all_files)
+	for f in $(all_files) ; do \
+	   echo Testing $$f; \
+	   ./$$f; \
+	done
+
+../$(ODIR)/libmpfr.a: ../mpfr.h
+	cd .. ; ./mmpfr
+
+tadd: tadd.c ../add.c
+	cd .. ; ./mmpfr tadd
+	-chmod g+w tadd
+
+tcmp: tcmp.c ../cmp.c
+	cd .. ; ./mmpfr tcmp
+	-chmod g+w tcmp
+
+tcmp2: tcmp2.c ../cmp.c
+	cd .. ; ./mmpfr tcmp2
+	-chmod g+w tcmp2
+
+tmul: tmul.c ../mul.c
+	cd .. ; ./mmpfr tmul
+	-chmod g+w tmul
+
+tdiv: tdiv.c ../div.c
+	cd .. ; ./mmpfr tdiv
+	-chmod g+w tdiv
+
+tround: tround.c ../round.c
+	cd .. ; ./mmpfr tround
+	-chmod g+w tround
+
+tset_f: tset_f.c ../set_f.c
+	cd .. ; ./mmpfr tset_f
+	-chmod g+w tset_f
+
+tset_d: tset_d.c ../set_d.c
+	cd .. ; ./mmpfr tset_d
+	-chmod g+w tset_d
+
+tset_si: tset_si.c ../set_si.c
+	cd .. ; ./mmpfr tset_si
+	-chmod g+w tset_si
+
+tset_str: tset_str.c ../set_str_raw.c
+	cd .. ; ./mmpfr tset_str
+	-chmod g+w tset_str
+
+tmul_ui: tmul_ui.c ../mul_ui.c
+	cd .. ; ./mmpfr tmul_ui
+	-chmod g+w tmul_ui
+
+tmul_2exp: tmul_2exp.c ../mul_2exp.c
+	cd .. ; ./mmpfr tmul_2exp
+	-chmod g+w tmul_2exp
+
+tagm: tagm.c ../agm.c
+	cd .. ; ./mmpfr tagm
+	-chmod g+w tagm
+
+clean: 
+	-rm *~ $(all_files)
+
diff --git a/tests/mon_fichier b/tests/mon_fichier
new file mode 100644
index 000000000..bbc5844f7
--- /dev/null
+++ b/tests/mon_fichier
@@ -0,0 +1,9 @@
+check for a=7.14111721922249727293e+141, b=1.19068418426567829913e+142
+mpfr_agm failed for a=7.14111721922249727292e+141, b=1.19068418426567829912e+142, rnd_mode=3
+expected result is 9.37191344337123761296e+141, got 9.37191344337124167449e+141 (3 ulp)
+
+check for a=1.10542733675580811212e+116, b=2.98141268622731060995e+116
+mpfr_agm failed for a=1.10542733675580811213e+116, b=2.98141268622731060995e+116, rnd_mode=0
+expected result is 1.92773215447859987843e+116, got 1.92773215447859952847e+116 (1 ulp)
+
+fin
diff --git a/tests/tadd b/tests/tadd
new file mode 100755
index 000000000..5c07eb3e1
--- /dev/null
+++ b/tests/tadd
diff --git a/tests/tadd.c b/tests/tadd.c
new file mode 100644
index 000000000..6957181c2
--- /dev/null
+++ b/tests/tadd.c
@@ -0,0 +1,238 @@
+/* 10^6 additions on a PII-400:
+precision       *    mpf_add mpfr_add(RNDZ/RNDN/RNDU)   maple   mupad
+53 bits       0.003       0.45       0.64/0.59/0.63
+100 bits                  0.52       0.76/0.75/0.77
+225                       0.54       1.14
+500 bits                  0.72       1.58/1.59/1.65
+1000 bits                 1.10       2.34
+2017                      1.87       3.31
+5025                      4.17       4.87
+10017                     7.69       7.52
+20017                     15.0       13.4
+50017                     57.8       37.8
+100017                    124.       105.
+*/  
+
+/* #define DEBUG */
+
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include "gmp.h"
+#include "mpfr.h"
+
+#define ABS(x) (((x)>0) ? (x) : (-x))
+
+double drand()
+{
+  double d; long int *i;
+
+  i = (long int*) &d;
+  i[0] = lrand48();
+  i[1] = lrand48();
+  if (lrand48()%2) d=-d; /* generates negative numbers */
+  return d;
+}
+
+/* checks that x+y gives the same results in double
+   and with mpfr with 53 bits of precision */
+void check(double x, double y, unsigned int rnd_mode, unsigned int px, 
+unsigned int py, unsigned int pz, double res)
+{
+  double z1,z2; mpfr_t xx,yy,zz;
+
+  /*  printf("x=%1.20e, y=%1.20e, rnd_mode=%u px=%u py=%u pz=%u\n",x,y,rnd_mode,
+      px, py, pz); */
+  mpfr_init2(xx, px);
+  mpfr_init2(yy, py);
+  mpfr_init2(zz, pz);
+  mpfr_set_d(xx, x, rnd_mode);
+  mpfr_set_d(yy, y, rnd_mode);
+  mpfr_add(zz, xx, yy, rnd_mode);
+  mpfr_set_machine_rnd_mode(rnd_mode);
+  z1 = (res==0.0) ? x+y : res;
+  z2 = mpfr_get_d(zz);
+  /* printf("x+y=%1.20e\n",z2); */
+  if (px==53 && py==53 && pz==53) res=1.0;
+  if (res!=0.0 && z1!=z2) {
+    printf("expected sum is %1.20e, got %1.20e\n",z1,z2);
+    printf("mpfr_add failed for x=%1.20e y=%1.20e with rnd_mode=%u\n",x,y,rnd_mode);
+    exit(1);
+  }
+  mpfr_clear(xx); mpfr_clear(yy); mpfr_clear(zz);
+}
+
+/* returns the number of ulp's between a and b */
+int ulp(a,b) double a,b;
+{
+  double eps=1.1102230246251565404e-16; /* 2^(-53) */
+  b = (a-b)/a; if (b<0) b = -b;
+  return (int) floor(b/eps);
+}
+
+void check2(x,px,y,py,pz,rnd_mode) double x,y; int px,py,pz,rnd_mode;
+{
+  mpfr_t xx, yy, zz; double z,z2; int u;
+
+#ifdef DEBUG
+printf("x=%1.20e,%d y=%1.20e,%d pz=%d,rnd=%d\n",x,px,y,py,pz,rnd_mode);
+#endif
+  mpfr_init2(xx,px); mpfr_init2(yy,py); mpfr_init2(zz,pz);
+  mpfr_set_d(xx, x, rnd_mode);
+  mpfr_set_d(yy, y, rnd_mode);
+  mpfr_add(zz, xx, yy, rnd_mode);
+  mpfr_set_machine_rnd_mode(rnd_mode);
+  z = x+y; z2=mpfr_get_d(zz); u=ulp(z,z2);
+#ifdef DEBUG
+  printf("x+y=%1.20e,%d (%d ulp) rnd_mode=%d\n",z2,pz,u,rnd_mode);
+  mpfr_set_d(zz, z2, rnd_mode); 
+  printf("i.e."); mpfr_print_raw(zz); putchar('\n');
+#endif
+    /* one ulp difference is possible due to composed rounding */
+  if (px>=53 && py>=53 && pz>=53 && ABS(u)>1) { 
+    printf("x=%1.20e,%d y=%1.20e,%d pz=%d,rnd=%d\n",x,px,y,py,pz,rnd_mode);
+    printf("got %1.20e\n",z2);
+    printf("result should be %1.20e (diff=%d ulp)\n",z,u);
+    mpfr_set_d(zz, z, rnd_mode);
+    printf("i.e."); mpfr_print_raw(zz); putchar('\n');
+    exit(1); }
+  mpfr_clear(xx); mpfr_clear(yy); mpfr_clear(zz);
+}
+
+check64()
+{
+  mpfr_t x, t, u;
+  mpfr_init2(x, 65); mpfr_init2(t, 65); mpfr_init2(u, 65);
+  mpfr_set_str_raw(x, "0.10011010101000110101010000000011001001001110001011101011111011101E623");
+  mpfr_set_str_raw(t, "0.10011010101000110101010000000011001001001110001011101011111011100E623");
+  mpfr_sub(u, x, t, GMP_RNDU);
+  if (mpfr_get_d(u) != 9.4349060620538533806e167) { /* 2^558 */
+    printf("Error (1) in mpfr_sub\n"); exit(1);
+  }
+  mpfr_init2(x, 64); mpfr_init2(t, 64); mpfr_init2(u, 64);
+  mpfr_set_str_raw(x, "0.1000011110101111011110111111000011101011101111101101101100000100E-220");
+  mpfr_set_str_raw(t, "0.1000011110101111011110111111000011101011101111101101010011111101E-220");
+  mpfr_add(u, x, t, GMP_RNDU);
+  if ((MANT(u)[0] & 1) != 1) { 
+    printf("error in mpfr_add with rnd_mode=GMP_RNDU\n");
+    printf("b=  "); mpfr_print_raw(x); putchar('\n');
+    printf("c=  "); mpfr_print_raw(t); putchar('\n');
+    printf("b+c="); mpfr_print_raw(u); putchar('\n');
+    exit(1);
+  }
+  mpfr_clear(x); mpfr_clear(t); mpfr_clear(u);
+}
+
+main(argc,argv) int argc; char *argv[];
+{
+  double x,y; int i,prec,rnd_mode,px,py,pz;
+
+  check64();
+  check(1.16809465359248765399e+196, 7.92883212101990665259e+196,
+	GMP_RNDU, 53, 53, 53, 0.0);
+  check(3.14553393112021279444e-67, 3.14553401015952024126e-67,
+	GMP_RNDU, 53, 53, 53, 0.0);
+  printf("Checking random precisions\n");
+  srand(getpid());
+  check2(2.26531902208967707071e+168,99,-2.67795218510613988524e+168,67,94,2);
+  check2(-1.21510626304662318398e+145,70,1.21367733647758957118e+145,65,61,3);
+  check2(2.73028857032080744543e+155,83,-1.16446121423113355603e+163,59,125,1);
+  check2(-4.38589520019641698848e+78,155,-1.09923643769309483415e+72,15,159,3);
+  check2(-1.49963910666191123860e+265,76,-2.30915090591874527520e-191,8,75,1);
+  check2(5.17945380930936917508e+112,119,1.11369077158813567738e+108,15,150,1);
+  check2(-2.66910493504493276454e-52,117,1.61188644159592323415e-52,61,68,1);
+  check2(-1.87427265794105342764e-57,175,1.76570844587489516446e+190,2,115,1);
+  check2(6.85523243386777784171e+107,187,-2.78148588123699111146e+48,87,178,3);
+  check2(-1.15706375390780417299e-135,94,-1.07455137477117851576e-129,66,111,2);
+  check2(-1.15706375390780417299e-135,94,-1.07455137477117851576e-129,66,111,3);
+  check2(-3.31624349995221499866e-22,107,-8.20150212714204839621e+156,79,99,3);
+  check2(-1.15706375390780417299e-135,94,-1.07455137477117851576e-129,66,111,3);
+  check2(-1.08007920352320089721e+150,63,1.77607317509426332389e+73,64,64,0);
+  check2(4.49465557237618783128e+53,108,-2.45103927353799477871e+48,60,105,0);
+  check2(3.25471707846623300604e-160,81,-7.93846654265839958715e-274,58,54,0);
+  check2(-8.88471912490080158206e+253,79,-7.84488427404526918825e+124,95,53,3);
+  check2(-2.18548638152863831959e-125,61,-1.22788940592412363564e+226,71,54,0);
+  check2(-7.94156823309993162569e+77,74,-5.26820160805275124882e+80,99,101,3);
+  check2(-3.85170653452493859064e+189,62,2.18827389706660314090e+158,94,106,3);
+  check2(1.07966151149311101950e+46,88,1.13198076934147457400e+46,67,53,0);
+  check2(3.36768223223409657622e+209,55,-9.61624007357265441884e+219,113,53,0);
+  check2(-6.47376909368979326475e+159,111,5.11127211034490340501e+159,99,62,3);
+  check2(-4.95229483271607845549e+220,110,-6.06992115033276044773e+213,109,55,0);
+  check2(-6.47376909368979326475e+159,74,5.11127211034490340501e+159,111,75,2);
+  check2(2.26531902208967707070e+168,99,-2.67795218510613988525e+168,67,94,2);
+  check2(-2.28886326552077479586e-188,67,3.41419438647157839320e-177,60,110,2);
+  check2(-2.66910493504493276454e-52,117,1.61188644159592323415e-52,61,68,1);
+  check2(2.90983392714730768886e+50,101,2.31299792168440591870e+50,74,105,1);
+  check2(2.72046257722708717791e+243,97,-1.62158447436486437113e+243,83,96,0);
+  /* checks random precisions */
+  for (i=0;i<1000000;i++) {
+#ifdef DEBUG
+printf("\nTest i=%d\n",i);
+#endif
+    px = 53 + (rand() % 64); 
+    py = 53 + (rand() % 64); 
+    pz = 53 + (rand() % 64); 
+    rnd_mode = rand() % 4;
+    do { x = drand(); } while (isnan(x));
+    do { y = drand(); } while (isnan(y));
+    check2(x,px,y,py,pz,rnd_mode);
+  }
+  printf("Checking double precision (53 bits)\n");
+  prec = (argc<2) ? 53 : atoi(argv[1]);
+  rnd_mode = (argc<3) ? -1 : atoi(argv[2]);
+  check(-8.22183238641455905806e-19, 7.42227178769761587878e-19,
+	GMP_RNDD, 53, 53, 53, 0.0);
+  check(5.82106394662028628236e+234, -5.21514064202368477230e+89,
+	GMP_RNDD, 53, 53, 53, 0.0);
+  check(5.72931679569871602371e+122, -5.72886070363264321230e+122,
+	GMP_RNDN, 53, 53, 53, 0.0);
+  check(-5.09937369394650450820e+238, 2.70203299854862982387e+250,
+	GMP_RNDD, 53, 53, 53, 0.0);
+  check(-2.96695924472363684394e+27, 1.22842938251111500000e+16,
+	GMP_RNDD, 53, 53, 53, 0.0);
+  check(1.74693641655743793422e-227, -7.71776956366861843469e-229,
+	GMP_RNDN, 53, 53, 53, 0.0);
+  check(-1.03432206392780011159e-125, 1.30127034799251347548e-133,
+	GMP_RNDN, 53, 53, 53, 0.0);
+  check(1.05824655795525779205e+71, -1.06022698059744327881e+71,
+	GMP_RNDZ, 53, 53, 53, 0.0);
+  check(-5.84204911040921732219e+240, 7.26658169050749590763e+240,
+	GMP_RNDD, 53, 53, 53, 0.0);
+  /* the following check double overflow */
+  /*  check(6.27557402141211962228e+307, 1.32141396570101687757e+308,
+      GMP_RNDZ, 53, 53, 53, 0.0); */
+  check(1.00944884131046636376e+221, 2.33809162651471520268e+215,
+	GMP_RNDN, 53, 53, 53, 0.0);
+  check(4.29232078932667367325e-278, 1.07735250473897938332e-281,
+	GMP_RNDU, 53, 53, 53, 0.0);
+  check(5.27584773801377058681e-80, 8.91207657803547196421e-91,
+	GMP_RNDN, 53, 53, 53, 0.0);
+  check(2.99280481918991653800e+272, 5.34637717585790933424e+271,
+	GMP_RNDN, 53, 53, 53, 0.0);
+  check(4.67302514390488041733e-184, 2.18321376145645689945e-190,
+	GMP_RNDN, 53, 53, 53, 0.0);
+  check(5.57294120336300389254e+71, 2.60596167942024924040e+65,
+	GMP_RNDZ, 53, 53, 53, 0.0);
+  x=6151626677899716.0; for (i=0;i<30;i++) x = 2.0*x;
+  check(x, 4938448004894539.0, GMP_RNDU, 53, 53, 53, 0.0);
+  check(1.23056185051606761523e-190, 1.64589756643433857138e-181,
+	GMP_RNDU, 53, 53, 53, 0.0);
+  check(2.93231171510175981584e-280, 3.26266919161341483877e-273,
+	GMP_RNDU, 53, 53, 53, 0.0);
+  check(5.76707395945001907217e-58, 4.74752971449827687074e-51,
+	GMP_RNDD, 53, 53, 53, 0.0);
+  check(277363943109.0, 11.0, GMP_RNDN, 53, 53, 53, 0.0);
+  /* test denormalized numbers too */
+  check(8.06294740693074521573e-310, 6.95250701071929654575e-310,
+	GMP_RNDU, 53, 53, 53, 0.0);
+  /* compares to results with double precision using machine arithmetic */
+  for (i=0;i<1000000;i++) {
+    x = drand(); 
+    y = drand();
+    if (ABS(x)>2.2e-307 && ABS(y)>2.2e-307 && x+y<1.7e+308 && x+y>-1.7e308) 
+      /* avoid denormalized numbers and overflows */
+      check(x, y, (rnd_mode==-1) ? lrand48()%4 : rnd_mode, 
+	    prec, prec, prec, 0.0);
+  } 
+}
+
diff --git a/tests/tagm b/tests/tagm
new file mode 100755
index 000000000..e54c62ada
--- /dev/null
+++ b/tests/tagm
diff --git a/tests/tagm.c b/tests/tagm.c
new file mode 100644
index 000000000..6a89b534e
--- /dev/null
+++ b/tests/tagm.c
@@ -0,0 +1,162 @@
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <sys/types.h>
+#include <unistd.h>
+#include "gmp.h"
+#include "mpfr.h"
+
+
+double drand()
+{
+  double d; long int *i;
+
+  i = (long int*) &d;
+  do {
+    i[0] = lrand48();
+    i[1] = lrand48();
+  /*if (lrand48()%2) d=-d; */ /* generates negative numbers */
+                              /* useless here */
+  } while ((d<1e-153)||(d>1e153));    /* to avoid underflow or overflow
+					 in double calculus in sqrt(u*v) */
+
+  return d;
+}
+
+
+/* returns the number of ulp's between a and b */
+int ulp(a,b) double a,b;
+{
+  double eps=1.1102230246251565404e-16; /* 2^(-53) */
+  b = (a-b)/a; if (b<0) b = -b;
+  return (int) floor(b/eps);
+}
+
+
+double max(double a,double b) {
+  if (a>=b)
+    return a;
+  return b;
+}
+
+double min(double a,double b) {
+   if (b>=a)
+    return a;
+  return b;
+}
+
+
+
+double dagm(double a, double b) { 
+  double u,v,tmpu,tmpv;
+  
+  if ((isnan(a))||(isnan(b)))
+    return a+b;
+
+  tmpv=max(a,b);
+  tmpu=min(a,b);
+
+  while (!(((tmpu==u)&&(tmpv==v))||(ulp(u,v)==0))) {
+    u=tmpu;
+    v=tmpv;
+    tmpu=sqrt(u*v);
+    tmpv=(u+v)/2.0;
+  }
+  /*  printf("difference : %i ulp\n",ulp(u,v)); */
+	 return u;
+}
+
+
+
+void check(double a, double b, unsigned char rnd_mode)
+{
+  mpfr_t ta, tb, tres;
+  double res1, res2;
+
+  mpfr_init2(ta, 53);
+  mpfr_init2(tb, 53);
+  mpfr_init2(tres, 53);
+  
+  mpfr_set_d(ta, a, rnd_mode);
+  mpfr_set_d(tb, b, rnd_mode);
+
+  mpfr_agm(tres, ta, tb, rnd_mode);
+  mpfr_set_machine_rnd_mode(rnd_mode);
+  
+  res1=dagm(a,b);
+  res2 = mpfr_get_d(tres);
+
+  if (res1!=res2 && (!isnan(res1) || !isnan(res2))) {
+    printf("mpfr_agm failed for a=%1.20e, b=%1.20e, rnd_mode=%d\n",a,b,rnd_mode);
+    printf("expected result is %1.20e, got %1.20e (%d ulp)\n",res1,res2,
+	   ulp(res2,res1));
+  }
+  else
+    printf("GOAL !!!\n");
+  mpfr_clear(ta); mpfr_clear(tb); mpfr_clear(tres); 
+  
+}
+/*
+  void main()
+{
+  int i;
+  double a,b,gd,pt;
+
+  check(2,1,GMP_RNDN);
+  check(6,4,GMP_RNDN);
+  check(62,61,GMP_RNDN);
+  check(0.5,1,GMP_RNDN);
+  check(1,2,GMP_RNDN);
+  check(234375765,234375000,GMP_RNDN);
+  check(8,1,GMP_RNDU);
+  check(1,44,GMP_RNDU); 
+  
+  srand(getpid());
+
+  for (i=0;i<100;i++) {
+   a = drand(); 
+   b = drand();
+
+   gd=max(a,b);
+   pt=min(a,b);
+
+   if (((gd-pt)/pt)<2) {
+     printf("check for a=%1.20e, b=%1.20e\n",a,b);
+     check(a, b, rand() % 4);
+     printf("\n");
+     }
+    else {
+     printf("suit pas hypotheses !\n");
+     check(a, b, rand() % 4);
+     } 
+   }
+     printf("fin\n");
+     }
+*/
+
+void main() {
+  mpfr_t a,b,res;
+  int p;
+  double op1,op2;
+  char * strptr,* expptr;
+  printf("debut\n");
+  p=53;
+  op1=1;
+  op2=2;
+
+  mpfr_init2(a,p);
+  mpfr_init2(b,p);
+  mpfr_init2(res,p);
+
+  mpfr_set_d(a,op1,GMP_RNDN);
+  mpfr_set_d(b,op2,GMP_RNDN);
+  mpfr_agm(res,a,b,GMP_RNDN);
+
+  printf("avant\n");
+  strptr= mpfr_get_str(NULL,expptr,10,0,res,GMP_RNDN);
+
+  printf("mag entre %f et %f a la precision %i :\n",op1,op2,p);
+  printf("%sE%s\n",strptr,expptr);
+}
+
+  
diff --git a/tests/tcmp b/tests/tcmp
new file mode 100755
index 000000000..faba0dc31
--- /dev/null
+++ b/tests/tcmp
diff --git a/tests/tcmp.c b/tests/tcmp.c
new file mode 100644
index 000000000..baf36ca6e
--- /dev/null
+++ b/tests/tcmp.c
@@ -0,0 +1,57 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "gmp.h"
+#include "longlong.h"
+#include "mpfr.h"
+
+double drand()
+{
+  double d; long int *i;
+
+  i = (long int*) &d;
+  i[0] = lrand48();
+  i[1] = lrand48();
+  if (lrand48()%2) d=-d; /* generates negative numbers */
+  return d;
+}
+
+void main()
+{
+  double x,y; mpfr_t xx,yy; int i,c;
+
+  mpfr_init2(xx, 65); mpfr_init2(yy, 65);
+  mpfr_set_str_raw(xx, "0.10011010101000110101010000000011001001001110001011101011111011101E623");
+  mpfr_set_str_raw(yy, "0.10011010101000110101010000000011001001001110001011101011111011100E623");
+  if (mpfr_cmp2(xx,yy)!=64) { printf("Error (1) in mpfr_cmp\n"); exit(1); }
+  mpfr_set_str_raw(xx, "0.10100010001110110111000010001000010011111101000100011101000011100");
+  mpfr_set_str_raw(yy, "0.10100010001110110111000010001000010011111101000100011101000011011");
+  if (mpfr_cmp2(xx,yy)!=64) { printf("Error (1) in mpfr_cmp\n"); exit(1); }
+  mpfr_init2(xx,53);
+  mpfr_init2(yy,200);
+  mpfr_set_d(xx, 1.0, 0);
+  mpfr_set_d(yy, 1.0, 0);
+  if (mpfr_cmp(xx,yy)!=0) {
+    printf("Error in mpfr_cmp: 1.0 != 1.0\n"); exit(1);
+  }
+  mpfr_set_prec(yy, 31, GMP_RNDZ);
+  mpfr_set_d(xx, 1.0000000002, 0);
+  mpfr_set_d(yy, 1.0, 0);
+  if (!(mpfr_cmp(xx,yy)>0)) {
+    printf("Error in mpfr_cmp: not 1.0000000002 > 1.0\n"); exit(1);
+  }
+  mpfr_set_prec(yy, 53, GMP_RNDZ);
+  for (i=0;i<1000000;) {    x=drand(); y=drand();
+    if (!isnan(x) && !isnan(y)) {
+      i++;
+      mpfr_set_d(xx, x, 0);
+      mpfr_set_d(yy, y, 0);
+      c = mpfr_cmp(xx,yy);
+      if ((c>0 && x<=y) || (c==0 && x!=y) || (c<0 && x>=y)) {
+	printf("Error in mpfr_cmp with x=%1.20e, y=%1.20e mpfr_cmp(x,y)=%d\n",
+	       x,y,c); exit(1);
+      }
+    }
+  }
+  mpfr_clear(xx); mpfr_clear(yy);
+}
diff --git a/tests/tcmp2 b/tests/tcmp2
new file mode 100755
index 000000000..1cb3e6073
--- /dev/null
+++ b/tests/tcmp2
diff --git a/tests/tcmp2.c b/tests/tcmp2.c
new file mode 100644
index 000000000..c4c080b97
--- /dev/null
+++ b/tests/tcmp2.c
@@ -0,0 +1,52 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "gmp.h"
+#include "longlong.h"
+#include "mpfr.h"
+
+double drand()
+{
+  double d; long int *i;
+
+  i = (long int*) &d;
+  i[0] = lrand48();
+  i[1] = lrand48();
+  return d;
+}
+
+void tcmp2(x, y, i) double x, y; int i;
+{
+  mpfr_t xx,yy; int j;
+
+  if (i==-1) i = (int) floor(log(x)/log(2.0)) - (int) floor(log(x-y)/log(2.0));
+  mpfr_init2(xx, 53); mpfr_init2(yy, 53);
+  mpfr_set_d(xx, x, 0);
+  mpfr_set_d(yy, y, 0);
+  if ((j=mpfr_cmp2(xx, yy)) != i) {
+    printf("Error in mpfr_cmp2: x=%1.16e y=%1.16e mpfr_cmp(x,y)=%d instead of %d\n",x,y,j,i); 
+    exit(1);
+  }
+  mpfr_clear(xx); mpfr_clear(yy);
+}
+
+void main()
+{
+  int i,j; double x=1.0, y, z;
+
+  tcmp2(1.06022698059744327881e+71, 1.05824655795525779205e+71, -1);
+  tcmp2(1.0, 1.0, 53);
+  for (i=0;i<54;i++) {
+    tcmp2(1.0, 1.0-x, i);
+    x /= 2.0;
+  }
+  for (j=0;j<1000000;j++) {
+    x = drand(); if (x<0) x = -x;
+    y = drand(); if (y<0) y = -y;
+    if (!isnan(x) && !isnan(y)) {
+      if (x<y) { z=x; x=y; y=z; }
+      tcmp2(x, y, -1);
+    }
+  }
+}
+
diff --git a/tests/tcmp_ui b/tests/tcmp_ui
new file mode 100755
index 000000000..9ad3d0cb8
--- /dev/null
+++ b/tests/tcmp_ui
diff --git a/tests/tcmp_ui.c b/tests/tcmp_ui.c
new file mode 100644
index 000000000..98642f7ca
--- /dev/null
+++ b/tests/tcmp_ui.c
@@ -0,0 +1,35 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "gmp.h"
+#include "longlong.h"
+#include "mpfr.h"
+
+main()
+{
+  mpfr_t x; unsigned long i; long s;
+  
+  mpfr_init(x);
+
+  mpfr_set_ui(x, 3, GMP_RNDZ);
+  if (mpfr_cmp_ui(x, i=3)!=0) {
+    printf("Error in mpfr_cmp_ui(%1.20f,%d)\n",mpfr_get_d(x), i); exit(1);
+  }
+  if (mpfr_cmp_ui(x, i=2)<=0) {
+    printf("Error in mpfr_cmp_ui(%1.20f,%d)\n",mpfr_get_d(x), i); exit(1);
+  }
+  if (mpfr_cmp_ui(x, i=4)>=0) {
+    printf("Error in mpfr_cmp_ui(%1.20f,%d)\n",mpfr_get_d(x), i); exit(1);
+  }
+
+  mpfr_set_si(x, -3, GMP_RNDZ);
+  if (mpfr_cmp_si(x, s=-3)!=0) {
+    printf("Error in mpfr_cmp_si(%1.20f,%d)\n",mpfr_get_d(x), s); exit(1);
+  }
+  if (mpfr_cmp_si(x, i=-4)<=0) {
+    printf("Error in mpfr_cmp_si(%1.20f,%d)\n",mpfr_get_d(x), s); exit(1);
+  }
+  if (mpfr_cmp_si(x, i=1)>=0) {
+    printf("Error in mpfr_cmp_si(%1.20f,%d)\n",mpfr_get_d(x), s); exit(1);
+  }
+}
diff --git a/tests/tdiv b/tests/tdiv
new file mode 100755
index 000000000..55fe33689
--- /dev/null
+++ b/tests/tdiv
diff --git a/tests/tdiv.c b/tests/tdiv.c
new file mode 100644
index 000000000..971970089
--- /dev/null
+++ b/tests/tdiv.c
@@ -0,0 +1,78 @@
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "mpfr.h"
+
+/* #define DEBUG */
+
+double drand()
+{
+  double d; long int *i;
+
+  i = (long int*) &d;
+  i[0] = lrand48();
+  i[1] = lrand48();
+  if (lrand48()%2) d=-d; /* generates negative numbers */
+  return d;
+}
+
+/* returns the number of ulp's between a and b */
+int ulp(a,b) double a,b;
+{
+  double eps=1.1102230246251565404e-16; /* 2^(-53) */
+  b = (a-b)/a; if (b<0) b = -b;
+  return (int) floor(b/eps);
+}
+
+#define check(n,d,r) check4(n,d,r,53)
+
+void check4(N, D, rnd_mode, p) double N, D; unsigned char rnd_mode; int p;
+{
+  mpfr_t q, n, d; double Q,Q2;
+
+#ifdef DEBUG
+  printf("N=%1.20e D=%1.20e rnd_mode=%d\n",N,D,rnd_mode);
+#endif
+  mpfr_init2(q, p); mpfr_init2(n, p); mpfr_init2(d, p);
+  mpfr_set_d(n, N, rnd_mode);
+  mpfr_set_d(d, D, rnd_mode);
+  mpfr_div(q, n, d, rnd_mode);
+  mpfr_set_machine_rnd_mode(rnd_mode);
+  Q = N/D;
+  Q2 = mpfr_get_d(q);
+#ifdef DEBUG
+    printf("expected quotient is %1.20e, got %1.20e (%d ulp)\n",Q,Q2,
+	   ulp(Q2,Q));
+    mpfr_print_raw(q); putchar('\n');
+#endif
+  if (Q!=Q2 && (!isnan(Q) || !isnan(Q2))) {
+    printf("mpfr_div failed for n=%1.20e, d=%1.20e, rnd_mode=%d\n",N,D,rnd_mode);
+    printf("expected quotient is %1.20e, got %1.20e (%d ulp)\n",Q,Q2,
+	   ulp(Q2,Q));
+    exit(1);
+  }
+  mpfr_clear(q); mpfr_clear(n); mpfr_clear(d);  
+}
+
+void main()
+{
+  int i; double n, d, e; 
+
+  check4(2.44394909079968374564e-150, 2.10263340267725788209e+187, 2, 65);
+exit(1);
+  /* the following tests when d is an exact power of two */
+  check(9.89438396044940256501e-134, 5.93472984109987421717e-67, 2);
+  check(9.89438396044940256501e-134, -5.93472984109987421717e-67, 2);
+  check(-4.53063926135729747564e-308, 7.02293374921793516813e-84, 3);
+  check(6.25089225176473806123e-01, -2.35527154824420243364e-230, 3);
+  check(6.52308934689126000000e+15, -1.62063546601505417497e+273, 0);
+  check(1.04636807108079349236e-189, 3.72295730823253012954e-292, 1);
+  for (i=0;i<100000;i++) {
+    do { n = drand(); d = drand(); e = fabs(n)/fabs(d); }
+    /* smallest normalized is 2^(-1022), largest is 2^(1023)*(2-2^(-52)) */
+    while (e>=1.7976931348623157081e308 || e<2.225073858507201383e-308);
+    check(n, d, rand() % 4);
+  }
+}
diff --git a/tests/tget_str b/tests/tget_str
new file mode 100755
index 000000000..4ac65d743
--- /dev/null
+++ b/tests/tget_str
diff --git a/tests/tget_str.c b/tests/tget_str.c
new file mode 100644
index 000000000..6093bce67
--- /dev/null
+++ b/tests/tget_str.c
@@ -0,0 +1,77 @@
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "gmp.h"
+#include "mpfr.h"
+#include <time.h>
+
+void print_double(d) double d;
+{
+  int e, i;
+
+  e = (int) ceil(log(fabs(d))/log(2.0));
+  /* d <= 2^e */
+  e -= 53;
+  if (e>0) for (i=0;i<e;i++) d /= 2.0;
+  else for (i=0;i<-e;i++) d *= 2.0;
+  printf("%1.0f*2^(%d)",d,e);
+}
+
+double drand()
+{
+  double d; long int *i;
+
+  i = (long int*) &d;
+  i[0] = lrand48();
+  i[1] = lrand48();
+  if (lrand48()%2) d=-d; /* generates negative numbers */
+  return d;
+}
+
+check(d, rnd) double d; unsigned char rnd;
+{
+  mpfr_t x; char *str, str2[30]; int l, l2;
+
+  mpfr_init2(x, 53);
+  mpfr_set_d(x, d, 53, 0, rnd);
+  str = mpfr_get_str(NULL, NULL, 10, 5, x, rnd);
+  mpfr_set_machine_rnd_mode(rnd);
+  sprintf(str2, "%1.4e", d);
+  l2 = strlen(str2);
+  l = strlen(str); 
+  if (l!=l2) printf("l=%d l2=%d\n",l,l2);
+  if (str2[l2-3]=='-' && str2[l2-2]=='0' && str2[l2-1]=='0')
+    str2[l2-3]='+'; /* rule for sign of exponent 0 ? */
+  if (strcmp(str, str2)) {
+    printf("Error in mpfr_get_str for d=%s=",str2);
+    print_double(d);
+    printf("\ngot %s\n", str);
+    exit(1);
+  }
+  mpfr_clear(x);
+  free(str);
+}
+
+int
+main(int argc, char **argv)
+{
+  int i; double d;
+
+  srand(getpid());
+  /* printf seems to round towards nearest in all cases, at least with gcc */
+  check(4.059650008e-83, 0);
+  check(-6.606499965302424244461355e233, 0);
+  check(-7.4, 0);
+  check(0.997, 0);
+  check(-4.53063926135729747564e-308, 0);
+  check(2.14478198760196000000e+16, 0);
+  check(7.02293374921793516813e-84, 0);
+  for (i=0;i<100000;i++) {
+    do { d = drand(); } while (isnan(d));
+    check(d, 0);
+  }
+}
+
+
+
diff --git a/tests/tmul b/tests/tmul
new file mode 100755
index 000000000..4e43f999a
--- /dev/null
+++ b/tests/tmul
diff --git a/tests/tmul.c b/tests/tmul.c
new file mode 100644
index 000000000..2f561ca58
--- /dev/null
+++ b/tests/tmul.c
@@ -0,0 +1,107 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include "gmp.h"
+#include "mpfr.h"
+
+#define MINNORM 2.2250738585072013831e-308 /* 2^(-1022), smallest normalized */
+
+/* 10^6 multiplications on a PII-400:
+precision       *    mpf_mul mpfr_mul(RNDZ/RNDN/RNDU)   maple   mupad
+53 bits       0.018    0.89       1.44/1.58/1.50       15.2[16]   17.1
+100 bits               1.66       2.01/2.38/2.27       20.0[30]   18.7
+200 bits               3.45       Seg. fault           27.6[60]   20.2
+225                    4.12       4.33
+500 bits              12.53       12.39                81.5[151]  29.3
+1000 bits             40.26       38.35               190.5[301]  58.6
+2000 bits                                                        169.3
+2017                  123.        119.
+5025 bits             552         544                   1860      918
+*/
+
+double drand()
+{
+  double d; long int *i;
+
+  i = (long int*) &d;
+  i[0] = lrand48();
+  i[1] = lrand48();
+  if (lrand48()%2) d=-d; /* generates negative numbers */
+  return d;
+}
+
+/* checks that x*y gives the same results in double
+   and with mpfr with 53 bits of precision */
+int check(double x, double y, unsigned int rnd_mode, unsigned int px, 
+unsigned int py, unsigned int pz, double res)
+{
+  double z1,z2,z3; mpfr_t xx,yy,zz; int i;
+  mpf_t xxx,yyy,zzz;
+
+  /* printf("x=%1.20e, y=%1.20e, rnd_mode=%u px=%u py=%u pz=%u\n",x,y,rnd_mode,
+     px, py, pz); */
+  mpfr_init2(xx, px);
+  mpfr_init2(yy, py);
+  mpfr_init2(zz, pz);
+  mpf_init2(xxx,px); mpf_init2(yyy,py); mpf_init2(zzz,pz);
+  mpf_set_d(xxx, x); mpf_set_d(yyy, y);
+  mpfr_set_d(xx, x, rnd_mode);
+  mpfr_set_d(yy, y, rnd_mode);
+for (i=0;i<1;i++)  mpfr_mul(zz, xx, yy, rnd_mode);
+  mpf_mul(zzz, xxx, yyy);
+  mpfr_set_machine_rnd_mode(rnd_mode);
+  z1 = (res==0.0) ? x*y : res;
+  z2 = mpfr_get_d(zz);
+  z3 = mpf_get_d(zzz);
+  if (px==53 && py==53 && pz==53) res=1.0;
+  if (res!=0.0 && z1!=z2 && (z1>=MINNORM || z1<=-MINNORM)) {
+    printf("expected product is %1.20e, got %1.20e\n",z1,z2);
+    printf("mpfr_mul failed for x=%1.20e y=%1.20e with rnd_mode=%u\n",x,y,rnd_mode);
+mpfr_print_raw(zz); putchar('\n');
+    exit(1);
+  }
+  mpfr_clear(xx); mpfr_clear(yy); mpfr_clear(zz);
+  mpf_clear(xxx); mpf_clear(yyy); mpf_clear(zzz);
+}
+
+main(argc,argv) int argc; char *argv[];
+{
+  double x,y,z; int i,prec,rnd_mode;
+
+  prec = (argc<2) ? 53 : atoi(argv[1]);
+  rnd_mode = (argc<3) ? -1 : atoi(argv[2]);
+  check(2.71331408349172961467e-08, -6.72658901114033715233e-165,
+	GMP_RNDZ, 53, 53, 53, 0.0);
+  x=0.31869277231188065; y=0.88642843322303122;
+  check(x, y, GMP_RNDZ, 53, 53, 53, 0.0);
+  x=8.47622108205396074254e-01; y=3.24039313247872939883e-01;
+  check(x, y, GMP_RNDU, 28, 45, 1, 0.5);
+  x=2.63978122803639081440e-01; 
+  y=5736014.0/8388608.0; /* 6.83786096444222835089e-01; */
+  check(x, y, GMP_RNDN, 34, 23, 31, 0.180504585267044603);
+  x=9.84891017624509146344e-01; /* rounded to 1.0 with prec=6 */
+  x=1.0;
+  y=1.18351709358762491320e-01;
+  check(x, y, GMP_RNDU, 6, 41, 36, 0.1183517093595583);
+  /* the following checks that rounding to nearest sets the last
+     bit to zero in case of equal distance */
+  check(67108865.0, 134217729.0, GMP_RNDN, 53, 53, 53, 0.0);
+  x=1.37399642157394197284e-01; y=2.28877275604219221350e-01;
+  check(x, y, GMP_RNDN, 49, 15, 32, 0.0314472340833162888);
+  x=4.03160720978664954828e-01; y=5.85483042917246621073e-01;
+  check(x, y, GMP_RNDZ, 51, 22, 32, 0.2360436821472831);
+  x=3.90798504668055102229e-14; y=9.85394674650308388664e-04;
+  check(x, y, GMP_RNDN, 46, 22, 12, 0.385027296503914762e-16);
+  x=4.58687081072827851358e-01; y=2.20543551472118792844e-01;
+  check(x, y, GMP_RNDN, 49, 3, 1, 0.125);
+  for (i=0;i<1000000;) {
+    x = drand();
+    y = drand();
+    z = x*y; if (z<0) z=-z;
+    if (z<1e+308 && z>1e-308) /* don't test overflow/underflow for now */
+      { i++;
+      check(x, y, (rnd_mode==-1) ? lrand48()%4 : rnd_mode, 
+	    prec, prec, prec, 0.0);
+      }
+  } 
+}
+
diff --git a/tests/tmul_2exp b/tests/tmul_2exp
new file mode 100755
index 000000000..849710aa3
--- /dev/null
+++ b/tests/tmul_2exp
diff --git a/tests/tmul_2exp.c b/tests/tmul_2exp.c
new file mode 100644
index 000000000..3f6055727
--- /dev/null
+++ b/tests/tmul_2exp.c
@@ -0,0 +1,25 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+#include "gmp.h"
+#include "mpfr.h"
+
+/* checks that x*y gives the same results in double
+   and with mpfr with 53 bits of precision */
+
+void main(argc,argv) int argc; char *argv[];
+{
+  double x, z; mpfr_t w;
+
+  mpfr_init2(w, 53); 
+
+  srand48(time(NULL)); 
+  x = drand48(); 
+  mpfr_set_d(w, x, 53, 1, 0);
+  mpfr_mul_2exp(w, w, 10, GMP_RNDZ); 
+  if (x != (z = mpfr_get_d(w)/1024))
+    {
+      fprintf(stderr, "%lf != %lf\n", x, z); 
+    };
+}
+
diff --git a/tests/tmul_ui b/tests/tmul_ui
new file mode 100755
index 000000000..473d7e0b1
--- /dev/null
+++ b/tests/tmul_ui
diff --git a/tests/tmul_ui.c b/tests/tmul_ui.c
new file mode 100644
index 000000000..e173f4e87
--- /dev/null
+++ b/tests/tmul_ui.c
@@ -0,0 +1,23 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include "gmp.h"
+#include "mpfr.h"
+#include "time.h"
+
+int
+main(int argc, char **argv)
+{
+  mpfr_t x;
+  
+  mpfr_init2(x, 53);
+  mpfr_set_d(x, 1.0/3.0, GMP_RNDZ); 
+  mpfr_print_raw(x); putchar('\n'); 
+  
+  mpfr_mul_ui(x, x, 3, GMP_RNDU); 
+  mpfr_print_raw(x); putchar('\n'); 
+
+  printf("%f\n", mpfr_get_d(x)); 
+
+  mpfr_clear(x); 
+  return(0);
+}
diff --git a/tests/tround b/tests/tround
new file mode 100755
index 000000000..10fd312dd
--- /dev/null
+++ b/tests/tround
diff --git a/tests/tround.c b/tests/tround.c
new file mode 100644
index 000000000..1a743781e
--- /dev/null
+++ b/tests/tround.c
@@ -0,0 +1,22 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include "gmp.h"
+#include "mpfr.h"
+
+main()
+{
+   mpfr_t x;
+
+   /* checks that rounds to nearest sets the last
+     bit to zero in case of equal distance */
+   mpfr_init2(x, 2);
+   mpfr_set_d(x, 5.0, 0);
+   if (mpfr_get_d(x) != 4.0) { printf("Error in tround: got %1.1f instead of 4.0\n",mpfr_get_d(x)); }
+
+   mpfr_set_d(x, 0.00098539467465030839, 0);
+
+   mpfr_set_d(x, 9.84891017624509146344e-01, GMP_RNDU); 
+   if (mpfr_get_d(x) != 1.0) { printf("Error in tround: got %f instead of 1.0\n",mpfr_get_d(x)); exit(1); }
+
+   mpfr_clear(x); 
+}
diff --git a/tests/tset_d b/tests/tset_d
new file mode 100755
index 000000000..f639320b0
--- /dev/null
+++ b/tests/tset_d
diff --git a/tests/tset_d.c b/tests/tset_d.c
new file mode 100644
index 000000000..54478c3f4
--- /dev/null
+++ b/tests/tset_d.c
@@ -0,0 +1,60 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include "gmp.h"
+#include "mpfr.h"
+#include <time.h>
+#include <math.h>
+
+double drand()
+{
+  double d; long int *i;
+
+  i = (long int*) &d;
+  i[0] = lrand48();
+  i[1] = lrand48();
+  return d;
+}
+
+int
+main(int argc, char **argv)
+{
+  mpfr_t x,y,z; unsigned long k,n; double d; 
+
+  mpfr_init2(z, 32);
+  mpfr_set_d(z, 1.0, 0);
+  if (mpfr_get_d(z) != 1.0) {
+    mpfr_print_raw(z); putchar('\n');
+    printf("Error: 1.0 != 1.0\n"); exit(1);
+  }
+  mpfr_init2(x, 53); mpfr_init2(y, 53);
+  mpfr_set_d(x, d=-1.08007920352320089721e+150, 0);
+  if (mpfr_get_d(x) != d) {
+    mpfr_print_raw(x); putchar('\n');
+    printf("Error: get_d o set_d <> identity for d = %1.20e %1.20e\n",d,
+	   mpfr_get_d(x)); exit(1);
+  }
+  srand48(time(NULL)); 
+  mpfr_set_d(x, 8.06294740693074521573e-310, 0); 
+  d = -6.72658901114033715233e-165;
+  mpfr_set_d(x, d, 0);
+  if (d != mpfr_get_d(x)) {
+    mpfr_print_raw(x); putchar('\n');
+    printf("Error: get_d o set_d <> identity for d = %1.20e %1.20e\n",d,
+	   mpfr_get_d(x)); exit(1);
+  }
+  n = (argc==1) ? 1000000 : atoi(argv[1]);
+  for (k = 1; k <= n; k++)
+    {      
+      do { d = drand(); } while (isnan(d)); /* does not yet work for NaN */
+      mpfr_set_d(x, d, 0); 
+      if (d != mpfr_get_d(x)) 
+	{ 
+	  fprintf(stderr, 
+		  "Mismatch on : %1.18g != %1.18g\n", d, mpfr_get_d(x)); 
+	  mpfr_print_raw(x); putchar('\n');
+	} 
+    }
+
+  mpfr_clear(x); 
+  return 0; 
+}
diff --git a/tests/tset_f b/tests/tset_f
new file mode 100755
index 000000000..1c0d8e334
--- /dev/null
+++ b/tests/tset_f
diff --git a/tests/tset_f.c b/tests/tset_f.c
new file mode 100644
index 000000000..8a8dc0a50
--- /dev/null
+++ b/tests/tset_f.c
@@ -0,0 +1,33 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include "gmp.h"
+#include "mpfr.h"
+#include "time.h"
+
+int
+main()
+{
+  mpfr_t x; mpf_t y; mpf_t z; unsigned long k, pr; 
+  
+  mpfr_init2(x, 100);
+  mpf_init(y); 
+
+  srandom(time(NULL)); 
+  mpf_random2(y, 10, 0); 
+  mpfr_set_f(x, y, 53, rand() & 3); 
+
+  mpf_out_str(stdout, 10, 0, y); 
+  printf("\n%1.21g\n", mpfr_get_d(x)); 
+
+  mpf_clear(y); mpfr_clear(x); 
+
+  for (k = 1; k <= 100000; k++)
+    {
+      pr = 1 + (rand()&255); 
+      mpf_init2(z, pr);
+      mpf_random2(z, z->_mp_prec, 0);
+      mpfr_init2(x, pr);
+      mpfr_set_f(x, z, pr, 0);
+    }
+  return(0);
+}
diff --git a/tests/tset_i b/tests/tset_i
new file mode 100755
index 000000000..4863b90fa
--- /dev/null
+++ b/tests/tset_i
diff --git a/tests/tset_si.c b/tests/tset_si.c
new file mode 100644
index 000000000..f3d3bcef9
--- /dev/null
+++ b/tests/tset_si.c
@@ -0,0 +1,37 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include "gmp.h"
+#include "mpfr.h"
+#include "time.h"
+
+int
+main(int argc, char **argv)
+{
+  mpfr_t x; long k, z, d; unsigned long zl, dl; 
+  
+  mpfr_init2(x, 100);
+
+  srandom(time(NULL)); 
+
+  for (k = 1; k <= atoi(argv[1]); k++)
+    {
+      z = random() - (1 << 30);      
+      mpfr_set_si(x, z, GMP_RNDZ); 
+      d = (int)mpfr_get_d(x);
+      if (d != z)
+	printf("Expected %ld got %ld\n", z, d); 
+      
+    }
+
+  for (k = 1; k <= atoi(argv[1]); k++)
+    {
+      zl = random();
+      mpfr_set_ui(x, zl, GMP_RNDZ); 
+      dl = (unsigned int) mpfr_get_d(x);
+      if (dl != zl)
+	printf("Expected %lu got %lu\n", zl, dl); 
+    }
+
+  mpfr_clear(x); 
+  return(0);
+}
diff --git a/tests/tset_str b/tests/tset_str
new file mode 100755
index 000000000..7d81d195d
--- /dev/null
+++ b/tests/tset_str
diff --git a/tests/tset_str.c b/tests/tset_str.c
new file mode 100644
index 000000000..e37544fe0
--- /dev/null
+++ b/tests/tset_str.c
@@ -0,0 +1,45 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include "gmp.h"
+#include "mpfr.h"
+#include <time.h>
+
+int
+main(int argc, char **argv)
+{
+  mpfr_t x; unsigned long k, bd, nc; char *str, *str2; 
+
+  srandom(time(NULL)); 
+  
+  if (argc > 1) { nc = atoi(argv[1]); } else { nc = 53; }
+  if (nc < 24) { nc = 24; }
+
+  bd = random()&8;
+  
+  str2 = str = (char *) malloc (nc*sizeof(char)); 
+
+  if (bd) 
+    {
+      for(k = 1; k <= bd; k++) 
+	{ *(str2++) = (random() & 1) + '0'; }
+    }
+  else { *(str2++) = '0'; }
+
+  *(str2++) = '.'; 
+
+  for(k = 1; k < nc - 17 - bd; k++)
+    {
+      *(str2++) = '0' + (random() & 1); 
+    }
+
+  *(str2++) = 'e'; 
+  sprintf(str2, "%d", random() - (1 << 30)); 
+
+  printf("%s\n", str); 
+  mpfr_init2(x, nc + 10); 
+  mpfr_set_str_raw(x, str); 
+  mpfr_print_raw(x); printf("\n"); 
+
+  mpfr_clear(x); free(str); 
+  return 0; 
+}
diff --git a/tests/tsqrt b/tests/tsqrt
new file mode 100755
index 000000000..8c802dce5
--- /dev/null
+++ b/tests/tsqrt
diff --git a/tests/tsqrt.c b/tests/tsqrt.c
new file mode 100644
index 000000000..c3af019c3
--- /dev/null
+++ b/tests/tsqrt.c
@@ -0,0 +1,63 @@
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "mpfr.h"
+
+double drand()
+{
+  double d; long int *i;
+
+  i = (long int*) &d;
+  i[0] = lrand48();
+  i[1] = lrand48();
+  if (lrand48()%2) d=-d; /* generates negative numbers */
+  return d;
+}
+
+/* returns the number of ulp's between a and b */
+int ulp(a,b) double a,b;
+{
+  double eps=1.1102230246251565404e-16; /* 2^(-53) */
+  b = (a-b)/a; if (b<0) b = -b;
+  return (int) floor(b/eps);
+}
+
+void check(a, rnd_mode) double a; unsigned char rnd_mode;
+{
+  mpfr_t q, n; double Q,Q2;
+
+#ifdef DEBUG
+  printf("a=%1.20e rnd_mode=%d\n",a,rnd_mode);
+#endif
+  mpfr_init2(q, 53); mpfr_init2(n, 53);
+  mpfr_set_d(n, a, 53, 0, rnd_mode);
+  mpfr_set_machine_rnd_mode(rnd_mode);
+  mpfr_sqrt(q, n, rnd_mode);
+  Q = sqrt(a);
+  Q2 = mpfr_get_d(q);
+#ifdef DEBUG
+    printf("expected sqrt is %1.20e, got %1.20e (%d ulp)\n",Q,Q2,
+	   ulp(Q2,Q));
+    mpfr_print_raw(q); putchar('\n');
+#endif
+  if (Q!=Q2 && (!isnan(Q) || !isnan(Q2))) {
+    printf("mpfr_sqrt failed for a=%1.20e, rnd_mode=%d\n",a,rnd_mode);
+    printf("expected sqrt is %1.20e, got %1.20e (%d ulp)\n",Q,Q2,
+	   ulp(Q2,Q));
+    exit(1);
+  }
+  mpfr_clear(q); mpfr_clear(n);
+}
+
+void main()
+{
+  int i; double a;
+
+  check(9.89438396044940256501e-134, 2);
+  for (i=0;i<1000000;i++) {
+    do { a = drand(); } while (isnan(a));
+    check(a, rand() % 4);
+  }
+}