added libtomfloat-0.02

changes.txt

@@ -1,3 +1,17 @@
-May 5th, 2004 v0.01 -- wrote base of LTF
+June 21st, 2004
+v0.02 -- Added missing objects to the makefile [oops]
+      -- fixed up mpf_add and mpf_sub to be more reliable about the precision
+      -- added required limited ln function mpf_ln_l (domain 0 < x <= 2)
+         It's still incomplete as it converges slowly (and therefore yields incorrect results)
+      -- Added mpf_ln and mpf_atan
+      -- Added short-circuits to sin, cos, invsqrt, atan, ln and sqrt [huge speedup]
+      -- Optimized mpf_sqrt and mpf_invsqrt by using quick estimates.  Fixed 
+         circular dependency as well (mpf_sqrt requires mpf_invsqrt but not vice versa)
+      ++ Note: No further releases are planned for a while.  I encourage interested
+         parties to fork this code base and extend it!
+
+May 5th, 2004 
+v0.01 -- wrote base of LTF
 
-Anytime       v0.00 -- no LTF existed.
+Anytime       
+v0.00 -- no LTF existed.

demos/ex1.c

@@ -10,7 +10,9 @@ void draw(mp_float *a)
 int main(void)
 {
   mp_float a, b, c, d, e;
-  mpf_init_multi(96, &a, &b, &c, &d, &e, NULL);
+  int err;
+
+  mpf_init_multi(100, &a, &b, &c, &d, &e, NULL);
 
   mpf_const_d(&a, 1); draw(&a);
   mpf_const_d(&b, 2); draw(&b);
@@ -21,6 +23,12 @@ int main(void)
   mpf_sub(&b, &c, &e); printf("2 - 3            ==");  draw(&e);
   mpf_mul(&b, &c, &e); printf("2 * 3            == "); draw(&e);
   mpf_div(&b, &c, &e); printf("2 / 3            == "); draw(&e);
+  mpf_add_d(&b, 3, &e); printf("2 + 3            == "); draw(&e);
+  mpf_sub_d(&b, 3, &e); printf("2 - 3            ==");  draw(&e);
+  mpf_mul_d(&b, 3, &e); printf("2 * 3            == "); draw(&e);
+  mpf_div_d	(&b, 3, &e); printf("2 / 3            == "); draw(&e);
+  mpf_const_d(&e, 0); mpf_add_d(&e, 1, &e); printf("0 + 1            == "); draw(&e);
+  mpf_const_d(&e, 0); mpf_sub_d(&e, 1, &e); printf("0 - 1            == "); draw(&e);
   printf("\n");
   mpf_invsqrt(&d, &e); printf("1/sqrt(4) == 1/2 == "); draw(&e);
   mpf_invsqrt(&c, &e); printf("1/sqrt(3)        == "); draw(&e);
@@ -29,6 +37,8 @@ int main(void)
   mpf_inv(&c, &e);     printf("1/3              == "); draw(&e);
   mpf_inv(&d, &e);     printf("1/4              == "); draw(&e);
   printf("\n");
+  mpf_const_pi(&e);    printf("Pi               == "); draw(&e);
+  printf("\n");
   mpf_const_e(&e);     printf("e                == "); draw(&e);
   mpf_exp(&c, &e);     printf("e^3              == "); draw(&e);
   mpf_sqrt(&e, &e);    printf("sqrt(e^3)        == "); draw(&e);
@@ -46,7 +56,21 @@ int main(void)
   mpf_tan(&b, &e);     printf("tan(2)           == "); draw(&e);
   mpf_tan(&c, &e);     printf("tan(3)           == "); draw(&e);
   mpf_tan(&d, &e);     printf("tan(4)           == "); draw(&e);
+  mpf_inv(&a, &e);     mpf_atan(&e, &e);    printf("atan(1/1)        == "); draw(&e);
+  mpf_inv(&b, &e);     mpf_atan(&e, &e);    printf("atan(1/2)        == "); draw(&e);
+  mpf_inv(&c, &e);     mpf_atan(&e, &e);    printf("atan(1/3)        == "); draw(&e);
+  mpf_inv(&d, &e);     mpf_atan(&e, &e);    printf("atan(1/4)        == "); draw(&e);
   printf("\n");
+  #define lntest(x) if ((err = mpf_const_ln_d(&e, x)) != MP_OKAY) { printf("Failed ln(%3d), %d\n", x, err); } else { printf("ln(%3d)          == ", x); draw(&e); };
+  lntest(0);
+  lntest(1);
+  lntest(2);
+  lntest(4);
+  lntest(8);
+  lntest(17);
+  lntest(1000);
+  lntest(100000);
+  lntest(250000);
   return 0; 
 }
 

float.tex

@@ -49,7 +49,7 @@
 \begin{document}
 \frontmatter
 \pagestyle{empty}
-\title{LibTomFloat User Manual \\ v0.01}
+\title{LibTomFloat User Manual \\ v0.02}
 \author{Tom St Denis \\ [email protected]}
 \maketitle
 This text and the library are hereby placed in the public domain.  This book has been formatted for B5 [176x250] paper using the \LaTeX{} {\em book} 

makefile

@@ -6,7 +6,7 @@ default: libtomfloat.a
 
 CFLAGS += -Os -Wall -W -I./
 
-VERSION=0.01
+VERSION=0.02
 
 #default files to install
 LIBNAME=libtomfloat.a
@@ -24,7 +24,7 @@ DATAPATH=/usr/share/doc/libtomfloat/pdf
 OBJECTS = \
 mpf_init.o mpf_clear.o mpf_init_multi.o mpf_clear_multi.o mpf_init_copy.o \
 \
-mpf_copy.o mpf_exch.o \
+mpf_copy.o mpf_exch.o mpf_abs.o mpf_neg.o \
 \
 mpf_cmp.o mpf_cmp_d.o \
 \
@@ -70,7 +70,7 @@ install: libtomfloat.a
 
 clean:
 	rm -f $(OBJECTS) libtomfloat.a *~ demos/*.o demos/*~ ex1
-	rm -f float.aux float.dvi float.log float.idx float.lof float.out float.toc
+	rm -f float.aux float.dvi float.log float.idx float.lof float.out float.toc float.ilg float.ind float.pdf
 
 zipup: clean manual
 	cd .. ; rm -rf ltf* libtomfloat-$(VERSION) ; mkdir libtomfloat-$(VERSION) ; \

mpf_add.c

@@ -18,6 +18,20 @@ int  mpf_add(mp_float *a, mp_float *b, mp_float *c)
    mp_float tmp, *other;
    long diff;
 
+   if (mpf_iszero(a)) {
+      diff = c->radix;
+      if ((err = mpf_copy(b, c)) != MP_OKAY) {
+         return err;
+      }
+      return mpf_normalize_to(c, diff);
+   } else if (mpf_iszero(b)) {
+      diff = c->radix;
+      if ((err = mpf_copy(a, c)) != MP_OKAY) {
+         return err;
+      }
+      return mpf_normalize_to(c, diff);
+   }
+
    if (a->exp < b->exp) {
       /* tmp == a normalize to b's exp */
       if ((err = mpf_init_copy(a, &tmp)) != MP_OKAY) {

mpf_atan.c

@@ -12,6 +12,85 @@
  */
 #include <tomfloat.h>
 
+/* y = y - (tan(y) - x)/(tan(y+0.1)-tan(x)) */
+
 int  mpf_atan(mp_float *a, mp_float *b)
 {
+   mp_float oldval, tmp, tmpx, res, sqr;
+   int      oddeven, ires, err, itts;
+   long     n;
+
+   /* ensure -1 <= a <= 1 */
+   if ((err = mpf_cmp_d(a, -1, &ires)) != MP_OKAY) {
+      return err; 
+   }
+   if (ires == MP_LT) {
+      return MP_VAL;
+   }
+
+   if ((err = mpf_cmp_d(a, 1, &ires)) != MP_OKAY) {
+      return err; 
+   }
+   if (ires == MP_GT) {
+      return MP_VAL;
+   }
+
+   /* easy out if a == 0 */
+   if (mpf_iszero(a) == MP_YES) {
+      return mpf_const_d(b, 1);
+   }
+
+   /* now a != 0 */
+
+   /* initialize temps */
+   if ((err = mpf_init_multi(b->radix, &oldval, &tmpx, &tmp, &res, &sqr, NULL)) != MP_OKAY) {
+      return err;
+   }
+
+   /* initlialize temps */
+   /* res = 0 */
+   /* tmpx = 1/a */
+   if ((err = mpf_inv(a, &tmpx)) != MP_OKAY)                                            { goto __ERR; }
+
+   /* sqr = a^2 */
+   if ((err = mpf_sqr(a, &sqr)) != MP_OKAY)                                             { goto __ERR; }
+
+   /* this is the denom counter.  Goes up by two per pass */
+   n       = 1;
+
+   /* we alternate between adding and subtracting */
+   oddeven = 0;
+
+   /* get number of iterations */
+   itts = mpf_iterations(b);
+
+   while (itts-- > 0) {
+       if ((err = mpf_copy(&res, &oldval)) != MP_OKAY)                                  { goto __ERR; }
+
+       /* compute 1/(2n-1) */
+       if ((err = mpf_const_d(&tmp, (2*n++ - 1))) != MP_OKAY)                           { goto __ERR; }
+       if ((err = mpf_inv(&tmp, &tmp)) != MP_OKAY)                                      { goto __ERR; }
+
+       /* now multiply a into tmpx twice */
+       if ((err = mpf_mul(&tmpx, &sqr, &tmpx)) != MP_OKAY)                              { goto __ERR; }
+
+       /* now multiply the two */
+       if ((err = mpf_mul(&tmpx, &tmp, &tmp)) != MP_OKAY)                               { goto __ERR; }
+
+       /* now depending on if this is even or odd we add/sub */
+       oddeven ^= 1;
+       if (oddeven  == 1) {
+          if ((err = mpf_add(&res, &tmp, &res)) != MP_OKAY)                             { goto __ERR; }
+       } else {
+          if ((err = mpf_sub(&res, &tmp, &res)) != MP_OKAY)                             { goto __ERR; }
+       }
+
+       if (mpf_cmp(&oldval, &res) == MP_EQ) {
+          break;
+       }
+   }
+   mpf_exch(&res, b);
+__ERR: mpf_clear_multi(&oldval, &tmpx, &tmp, &res, &sqr, NULL);
+   return err;
+
 }

mpf_const_ln_d.c

@@ -15,6 +15,20 @@
 int  mpf_const_ln_d(mp_float *a, long b)
 {
    int err;
+
+   /* test input */
+   if (b < 0) {
+      return MP_VAL;
+   }
+
+   if (b == 0) {
+      return mpf_const_d(a, 1);
+   }
+
+   if (b == 1) {
+      return mpf_const_d(a, 0);
+   }
+
    if ((err = mpf_const_d(a, b)) != MP_OKAY) {
       return err;
    }

mpf_cos.c

@@ -15,11 +15,11 @@
 /* using cos x == \sum_{n=0}^{\infty} ((-1)^n/(2n)!) * x^2n */
 int  mpf_cos(mp_float *a, mp_float *b)
 {
-   mp_float tmpovern, tmp, tmpx, res, sqr;
+   mp_float oldval, tmpovern, tmp, tmpx, res, sqr;
    int      oddeven, err, itts;
    long     n;
    /* initialize temps */
-   if ((err = mpf_init_multi(b->radix, &tmpx, &tmpovern, &tmp, &res, &sqr, NULL)) != MP_OKAY) {
+   if ((err = mpf_init_multi(b->radix, &oldval, &tmpx, &tmpovern, &tmp, &res, &sqr, NULL)) != MP_OKAY) {
       return err;
    }
 
@@ -40,6 +40,7 @@ int  mpf_cos(mp_float *a, mp_float *b)
    itts = mpf_iterations(b);
 
    while (itts-- > 0) {
+       if ((err = mpf_copy(&res, &oldval)) != MP_OKAY)                                  { goto __ERR; }
        /* compute 1/(2n)! from 1/(2(n-1))! by multiplying by (1/n)(1/(n+1)) */
        if ((err = mpf_const_d(&tmp, ++n)) != MP_OKAY)                                   { goto __ERR; }
        if ((err = mpf_inv(&tmp, &tmp)) != MP_OKAY)                                      { goto __ERR; }
@@ -62,8 +63,12 @@ int  mpf_cos(mp_float *a, mp_float *b)
        } else {
           if ((err = mpf_sub(&res, &tmp, &res)) != MP_OKAY)                             { goto __ERR; }
        }
+
+       if (mpf_cmp(&res, &oldval) == MP_EQ) {
+          break;
+       }
    }
    mpf_exch(&res, b);
-__ERR: mpf_clear_multi(&tmpx, &tmpovern, &tmp, &res, &sqr, NULL);
+__ERR: mpf_clear_multi(&oldval, &tmpx, &tmpovern, &tmp, &res, &sqr, NULL);
    return err;
 }

mpf_exp.c

@@ -16,12 +16,12 @@
 /* compute b = e^a using e^x == \sum_{n=0}^{\infty} {1 \over n!}x^n */
 int  mpf_exp(mp_float *a, mp_float *b)
 {
-   mp_float  tmpx, tmpovern, tmp, res; 
+   mp_float  oldval, tmpx, tmpovern, tmp, res; 
    int       err, itts;
    long      n;
 
    /* initialize temps */
-   if ((err = mpf_init_multi(b->radix, &tmpx, &tmpovern, &tmp, &res, NULL)) != MP_OKAY) {
+   if ((err = mpf_init_multi(b->radix, &oldval, &tmpx, &tmpovern, &tmp, &res, NULL)) != MP_OKAY) {
       return err;
    }
 
@@ -36,8 +36,9 @@ int  mpf_exp(mp_float *a, mp_float *b)
    itts = mpf_iterations(b);
 
    while (itts-- > 0) {
+       if ((err = mpf_copy(&res, &oldval)) != MP_OKAY)                                  { goto __ERR; }
+
        /* compute 1/n! as 1/(n-1)! * 1/n */
-// hack: this won't be portable for n>127
        if ((err = mpf_const_d(&tmp, n++)) != MP_OKAY)                                   { goto __ERR; }
        if ((err = mpf_inv(&tmp, &tmp)) != MP_OKAY)                                      { goto __ERR; }
        if ((err = mpf_mul(&tmp, &tmpovern, &tmpovern)) != MP_OKAY)                      { goto __ERR; }
@@ -48,10 +49,14 @@ int  mpf_exp(mp_float *a, mp_float *b)
        /* multiply and sum them */
        if ((err = mpf_mul(&tmpovern, &tmpx, &tmp)) != MP_OKAY)                          { goto __ERR; }
        if ((err = mpf_add(&tmp, &res, &res)) != MP_OKAY)                                { goto __ERR; }
+
+       if (mpf_cmp(&oldval, &res) == MP_EQ) {
+          break;
+       }
    }
 
    mpf_exch(&res, b);
-__ERR: mpf_clear_multi(&tmpx, &tmpovern, &tmp, &res, NULL);
+__ERR: mpf_clear_multi(&oldval, &tmpx, &tmpovern, &tmp, &res, NULL);
    return err;
 }
 

mpf_invsqrt.c

@@ -15,7 +15,7 @@
 /* using newtons method we have 1/sqrt(x) = Y_{n+1} = y_n * ((3 - xy^2_n)/2) */
 int  mpf_invsqrt(mp_float *a, mp_float *b)
 {
-   mp_float tmp1, tmp2, const_3;
+   mp_float oldval, tmp1, tmp2, const_3;
    int err, itts;
 
    /* ensure a is not zero or negative */
@@ -27,7 +27,7 @@ int  mpf_invsqrt(mp_float *a, mp_float *b)
    itts = mpf_iterations(b);
 
    /* init temps */
-   if ((err = mpf_init_multi(b->radix, &tmp1, &tmp2, &const_3, NULL)) != MP_OKAY) {
+   if ((err = mpf_init_multi(b->radix, &oldval, &tmp1, &tmp2, &const_3, NULL)) != MP_OKAY) {
       return err;
    }
 
@@ -36,13 +36,13 @@ int  mpf_invsqrt(mp_float *a, mp_float *b)
 
    /* tmp1 == reasonable guess at sqrt */
    if ((err = mpf_copy(a, &tmp1)) != MP_OKAY)                                     { goto __ERR; }
-
-   /* negate exponent and halve */
-   tmp1.radix = b->radix;
-   if ((err = mp_sqrt(&(tmp1.mantissa), &(tmp1.mantissa))) != MP_OKAY)            { goto __ERR; }
-   if ((err = mpf_normalize(&tmp1)) != MP_OKAY)                                   { goto __ERR; }
+   mp_rshd(&(tmp1.mantissa), tmp1.mantissa.used>>1);
+   if ((err = mpf_normalize_to(&tmp1, b->radix)) != MP_OKAY)                      { goto __ERR; }
 
    while (itts-- > 0) {
+       /* grap copy of tmp1 for early out */
+       if ((err = mpf_copy(&tmp1, &oldval)) != MP_OKAY)                           { goto __ERR; }
+
        /* first tmp2 = y^2 == tmp1^2 */
        if ((err = mpf_sqr(&tmp1, &tmp2)) != MP_OKAY)                              { goto __ERR; }
        /* multiply by x, tmp1 * a */
@@ -53,10 +53,15 @@ int  mpf_invsqrt(mp_float *a, mp_float *b)
        if ((err = mpf_div_2(&tmp2, &tmp2)) != MP_OKAY)                            { goto __ERR; }
        /* multiply by y_n and feedback */
        if ((err = mpf_mul(&tmp1, &tmp2, &tmp1)) != MP_OKAY)                       { goto __ERR; }
+
+       /* early out if stable */
+       if (mpf_cmp(&oldval, &tmp1) == MP_EQ) {
+          break;
+       }
    }
 
    mpf_exch(&tmp1, b);
-__ERR: mpf_clear_multi(&tmp1, &tmp2, &const_3, NULL);
+__ERR: mpf_clear_multi(&oldval, &tmp1, &tmp2, &const_3, NULL);
    return err;
 }