now conforms to the specification

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

1 files changed, 41 insertions, 31 deletions

diff --git a/get_str.c b/get_str.c
index 60905159a..6a1af0452 100644
--- a/get_str.c
+++ b/get_str.c
@@ -2,22 +2,33 @@
 #include <math.h>
 #include "gmp.h"
 #include "gmp-impl.h"
+#include "longlong.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'.
+
+  For op = 3.1416 we get str = "31416" and expptr=1.
  */
 
 /* #define DEBUG */
 
-char *mpfr_get_str(char *str, char *expptr, int base, size_t n,
+char *mpfr_get_str(char *str, mp_exp_t *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;
+  double d; long e, q, i, neg, p, err, prec, sh; mpfr_t a, b; mpz_t bz;
+  char *str0; unsigned char rnd1; int f, pow2;
+
+  if (base<2 || 36<base) {
+    fprintf(stderr, "Error: too small or too large base in mpfr_get_str: %d\n",
+	    base);
+    exit(1);
+  }
+  count_leading_zeros(pow2, base); pow2 = mp_bits_per_limb - pow2 - 1;
+  if (base != (1<<pow2)) pow2=0; 
+  /* if pow2 <> 0, then base = 2^pow2 */
 
 #ifdef DEBUG
   printf("op="); mpfr_print_raw(op); printf(" rnd_mode=%d\n",rnd_mode);
@@ -34,6 +45,9 @@ char *mpfr_get_str(char *str, char *expptr, int base, size_t n,
 #ifdef DEBUG
   printf("exponent = %d\n",f);
 #endif
+  if (n==0)
+    n = (int) ceil((double)PREC(op)*log(2.0)/log((double)base)+
+		   log(4.0*fabs((double)((f==0) ? 1 : f)))/log(2.0));
   /* 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));
@@ -46,27 +60,26 @@ char *mpfr_get_str(char *str, char *expptr, int base, size_t n,
   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);
+    if (pow2) {
+      if (neg) mpfr_div_2exp(b, op, pow2*p, GMP_RNDZ);
+      else mpfr_mul_2exp(b, op, pow2*p, 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);
+       /* 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_ui_pow_ui(a, base, p, 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);
        }
-       /* 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);
+  } while (pow2==0 && 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;
@@ -96,20 +109,17 @@ 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);
+  q = ((SIGN(op)<0) ? 1 : 0) + n + 1;
+  if (str==NULL) str0=str=malloc(q); 
   if (SIGN(op)<0) *str++='-';
-  if (n>1) *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]='.'; }
+  /*  str[n++] = 'e'; */
+  /* str[n++] = (f>=0) ? '+' : '-'; */ /* is there a rule for f=0 ? */
+  /*  if (f<10 && f>-10) str[n++]='0'; */
+  *expptr = f;
+  /*  if (str[-1]=='.') { str[-1]=str[0]; str[0]='.'; } */
   mpfr_clear(a); mpfr_clear(b); mpz_clear(bz);
   return str0;
 }