Back to index

glibc  2.9
mpa.c
Go to the documentation of this file.
00001 
00002 /*
00003  * IBM Accurate Mathematical Library
00004  * written by International Business Machines Corp.
00005  * Copyright (C) 2001 Free Software Foundation
00006  *
00007  * This program is free software; you can redistribute it and/or modify
00008  * it under the terms of the GNU Lesser General Public License as published by
00009  * the Free Software Foundation; either version 2.1 of the License, or
00010  * (at your option) any later version.
00011  *
00012  * This program is distributed in the hope that it will be useful,
00013  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00014  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015  * GNU Lesser General Public License for more details.
00016  *
00017  * You should have received a copy of the GNU Lesser General Public License
00018  * along with this program; if not, write to the Free Software
00019  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
00020  */
00021 /************************************************************************/
00022 /*  MODULE_NAME: mpa.c                                                  */
00023 /*                                                                      */
00024 /*  FUNCTIONS:                                                          */
00025 /*               mcr                                                    */
00026 /*               acr                                                    */
00027 /*               cr                                                     */
00028 /*               cpy                                                    */
00029 /*               cpymn                                                  */
00030 /*               norm                                                   */
00031 /*               denorm                                                 */
00032 /*               mp_dbl                                                 */
00033 /*               dbl_mp                                                 */
00034 /*               add_magnitudes                                         */
00035 /*               sub_magnitudes                                         */
00036 /*               add                                                    */
00037 /*               sub                                                    */
00038 /*               mul                                                    */
00039 /*               inv                                                    */
00040 /*               dvd                                                    */
00041 /*                                                                      */
00042 /* Arithmetic functions for multiple precision numbers.                 */
00043 /* Relative errors are bounded                                          */
00044 /************************************************************************/
00045 
00046 
00047 #include "endian.h"
00048 #include "mpa.h"
00049 #include "mpa2.h"
00050 #include <sys/param.h>      /* For MIN() */
00051 /* mcr() compares the sizes of the mantissas of two multiple precision  */
00052 /* numbers. Mantissas are compared regardless of the signs of the       */
00053 /* numbers, even if x->d[0] or y->d[0] are zero. Exponents are also     */
00054 /* disregarded.                                                         */
00055 static int mcr(const mp_no *x, const mp_no *y, int p) {
00056   int i;
00057   for (i=1; i<=p; i++) {
00058     if      (X[i] == Y[i])  continue;
00059     else if (X[i] >  Y[i])  return  1;
00060     else                    return -1; }
00061   return 0;
00062 }
00063 
00064 
00065 
00066 /* acr() compares the absolute values of two multiple precision numbers */
00067 int __acr(const mp_no *x, const mp_no *y, int p) {
00068   int i;
00069 
00070   if      (X[0] == ZERO) {
00071     if    (Y[0] == ZERO) i= 0;
00072     else                 i=-1;
00073   }
00074   else if (Y[0] == ZERO) i= 1;
00075   else {
00076     if      (EX >  EY)   i= 1;
00077     else if (EX <  EY)   i=-1;
00078     else                 i= mcr(x,y,p);
00079   }
00080 
00081   return i;
00082 }
00083 
00084 
00085 /* cr90 compares the values of two multiple precision numbers           */
00086 int  __cr(const mp_no *x, const mp_no *y, int p) {
00087   int i;
00088 
00089   if      (X[0] > Y[0])  i= 1;
00090   else if (X[0] < Y[0])  i=-1;
00091   else if (X[0] < ZERO ) i= __acr(y,x,p);
00092   else                   i= __acr(x,y,p);
00093 
00094   return i;
00095 }
00096 
00097 
00098 /* Copy a multiple precision number. Set *y=*x. x=y is permissible.      */
00099 void __cpy(const mp_no *x, mp_no *y, int p) {
00100   int i;
00101 
00102   EY = EX;
00103   for (i=0; i <= p; i++)    Y[i] = X[i];
00104 
00105   return;
00106 }
00107 
00108 
00109 /* Copy a multiple precision number x of precision m into a */
00110 /* multiple precision number y of precision n. In case n>m, */
00111 /* the digits of y beyond the m'th are set to zero. In case */
00112 /* n<m, the digits of x beyond the n'th are ignored.        */
00113 /* x=y is permissible.                                      */
00114 
00115 void __cpymn(const mp_no *x, int m, mp_no *y, int n) {
00116 
00117   int i,k;
00118 
00119   EY = EX;     k=MIN(m,n);
00120   for (i=0; i <= k; i++)    Y[i] = X[i];
00121   for (   ; i <= n; i++)    Y[i] = ZERO;
00122 
00123   return;
00124 }
00125 
00126 /* Convert a multiple precision number *x into a double precision */
00127 /* number *y, normalized case  (|x| >= 2**(-1022))) */
00128 static void norm(const mp_no *x, double *y, int p)
00129 {
00130   #define R  radixi.d
00131   int i;
00132 #if 0
00133   int k;
00134 #endif
00135   double a,c,u,v,z[5];
00136   if (p<5) {
00137     if      (p==1) c = X[1];
00138     else if (p==2) c = X[1] + R* X[2];
00139     else if (p==3) c = X[1] + R*(X[2]  +   R* X[3]);
00140     else if (p==4) c =(X[1] + R* X[2]) + R*R*(X[3] + R*X[4]);
00141   }
00142   else {
00143     for (a=ONE, z[1]=X[1]; z[1] < TWO23; )
00144         {a *= TWO;   z[1] *= TWO; }
00145 
00146     for (i=2; i<5; i++) {
00147       z[i] = X[i]*a;
00148       u = (z[i] + CUTTER)-CUTTER;
00149       if  (u > z[i])  u -= RADIX;
00150       z[i] -= u;
00151       z[i-1] += u*RADIXI;
00152     }
00153 
00154     u = (z[3] + TWO71) - TWO71;
00155     if (u > z[3])   u -= TWO19;
00156     v = z[3]-u;
00157 
00158     if (v == TWO18) {
00159       if (z[4] == ZERO) {
00160         for (i=5; i <= p; i++) {
00161           if (X[i] == ZERO)   continue;
00162           else                {z[3] += ONE;   break; }
00163         }
00164       }
00165       else              z[3] += ONE;
00166     }
00167 
00168     c = (z[1] + R *(z[2] + R * z[3]))/a;
00169   }
00170 
00171   c *= X[0];
00172 
00173   for (i=1; i<EX; i++)   c *= RADIX;
00174   for (i=1; i>EX; i--)   c *= RADIXI;
00175 
00176   *y = c;
00177   return;
00178 #undef R
00179 }
00180 
00181 /* Convert a multiple precision number *x into a double precision */
00182 /* number *y, denormalized case  (|x| < 2**(-1022))) */
00183 static void denorm(const mp_no *x, double *y, int p)
00184 {
00185   int i,k;
00186   double c,u,z[5];
00187 #if 0
00188   double a,v;
00189 #endif
00190 
00191 #define R  radixi.d
00192   if (EX<-44 || (EX==-44 && X[1]<TWO5))
00193      { *y=ZERO; return; }
00194 
00195   if      (p==1) {
00196     if      (EX==-42) {z[1]=X[1]+TWO10;  z[2]=ZERO;  z[3]=ZERO;  k=3;}
00197     else if (EX==-43) {z[1]=     TWO10;  z[2]=X[1];  z[3]=ZERO;  k=2;}
00198     else              {z[1]=     TWO10;  z[2]=ZERO;  z[3]=X[1];  k=1;}
00199   }
00200   else if (p==2) {
00201     if      (EX==-42) {z[1]=X[1]+TWO10;  z[2]=X[2];  z[3]=ZERO;  k=3;}
00202     else if (EX==-43) {z[1]=     TWO10;  z[2]=X[1];  z[3]=X[2];  k=2;}
00203     else              {z[1]=     TWO10;  z[2]=ZERO;  z[3]=X[1];  k=1;}
00204   }
00205   else {
00206     if      (EX==-42) {z[1]=X[1]+TWO10;  z[2]=X[2];  k=3;}
00207     else if (EX==-43) {z[1]=     TWO10;  z[2]=X[1];  k=2;}
00208     else              {z[1]=     TWO10;  z[2]=ZERO;  k=1;}
00209     z[3] = X[k];
00210   }
00211 
00212   u = (z[3] + TWO57) - TWO57;
00213   if  (u > z[3])   u -= TWO5;
00214 
00215   if (u==z[3]) {
00216     for (i=k+1; i <= p; i++) {
00217       if (X[i] == ZERO)   continue;
00218       else {z[3] += ONE;   break; }
00219     }
00220   }
00221 
00222   c = X[0]*((z[1] + R*(z[2] + R*z[3])) - TWO10);
00223 
00224   *y = c*TWOM1032;
00225   return;
00226 
00227 #undef R
00228 }
00229 
00230 /* Convert a multiple precision number *x into a double precision number *y. */
00231 /* The result is correctly rounded to the nearest/even. *x is left unchanged */
00232 
00233 void __mp_dbl(const mp_no *x, double *y, int p) {
00234 #if 0
00235   int i,k;
00236   double a,c,u,v,z[5];
00237 #endif
00238 
00239   if (X[0] == ZERO)  {*y = ZERO;  return; }
00240 
00241   if      (EX> -42)                 norm(x,y,p);
00242   else if (EX==-42 && X[1]>=TWO10)  norm(x,y,p);
00243   else                              denorm(x,y,p);
00244 }
00245 
00246 
00247 /* dbl_mp() converts a double precision number x into a multiple precision  */
00248 /* number *y. If the precision p is too small the result is truncated. x is */
00249 /* left unchanged.                                                          */
00250 
00251 void __dbl_mp(double x, mp_no *y, int p) {
00252 
00253   int i,n;
00254   double u;
00255 
00256   /* Sign */
00257   if      (x == ZERO)  {Y[0] = ZERO;  return; }
00258   else if (x >  ZERO)   Y[0] = ONE;
00259   else                 {Y[0] = MONE;  x=-x;   }
00260 
00261   /* Exponent */
00262   for (EY=ONE; x >= RADIX; EY += ONE)   x *= RADIXI;
00263   for (      ; x <  ONE;   EY -= ONE)   x *= RADIX;
00264 
00265   /* Digits */
00266   n=MIN(p,4);
00267   for (i=1; i<=n; i++) {
00268     u = (x + TWO52) - TWO52;
00269     if (u>x)   u -= ONE;
00270     Y[i] = u;     x -= u;    x *= RADIX; }
00271   for (   ; i<=p; i++)     Y[i] = ZERO;
00272   return;
00273 }
00274 
00275 
00276 /*  add_magnitudes() adds the magnitudes of *x & *y assuming that           */
00277 /*  abs(*x) >= abs(*y) > 0.                                                 */
00278 /* The sign of the sum *z is undefined. x&y may overlap but not x&z or y&z. */
00279 /* No guard digit is used. The result equals the exact sum, truncated.      */
00280 /* *x & *y are left unchanged.                                              */
00281 
00282 static void add_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
00283 
00284   int i,j,k;
00285 
00286   EZ = EX;
00287 
00288   i=p;    j=p+ EY - EX;    k=p+1;
00289 
00290   if (j<1)
00291      {__cpy(x,z,p);  return; }
00292   else   Z[k] = ZERO;
00293 
00294   for (; j>0; i--,j--) {
00295     Z[k] += X[i] + Y[j];
00296     if (Z[k] >= RADIX) {
00297       Z[k]  -= RADIX;
00298       Z[--k] = ONE; }
00299     else
00300       Z[--k] = ZERO;
00301   }
00302 
00303   for (; i>0; i--) {
00304     Z[k] += X[i];
00305     if (Z[k] >= RADIX) {
00306       Z[k]  -= RADIX;
00307       Z[--k] = ONE; }
00308     else
00309       Z[--k] = ZERO;
00310   }
00311 
00312   if (Z[1] == ZERO) {
00313     for (i=1; i<=p; i++)    Z[i] = Z[i+1]; }
00314   else   EZ += ONE;
00315 }
00316 
00317 
00318 /*  sub_magnitudes() subtracts the magnitudes of *x & *y assuming that      */
00319 /*  abs(*x) > abs(*y) > 0.                                                  */
00320 /* The sign of the difference *z is undefined. x&y may overlap but not x&z  */
00321 /* or y&z. One guard digit is used. The error is less than one ulp.         */
00322 /* *x & *y are left unchanged.                                              */
00323 
00324 static void sub_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
00325 
00326   int i,j,k;
00327 
00328   EZ = EX;
00329 
00330   if (EX == EY) {
00331     i=j=k=p;
00332     Z[k] = Z[k+1] = ZERO; }
00333   else {
00334     j= EX - EY;
00335     if (j > p)  {__cpy(x,z,p);  return; }
00336     else {
00337       i=p;   j=p+1-j;   k=p;
00338       if (Y[j] > ZERO) {
00339         Z[k+1] = RADIX - Y[j--];
00340         Z[k]   = MONE; }
00341       else {
00342         Z[k+1] = ZERO;
00343         Z[k]   = ZERO;   j--;}
00344     }
00345   }
00346 
00347   for (; j>0; i--,j--) {
00348     Z[k] += (X[i] - Y[j]);
00349     if (Z[k] < ZERO) {
00350       Z[k]  += RADIX;
00351       Z[--k] = MONE; }
00352     else
00353       Z[--k] = ZERO;
00354   }
00355 
00356   for (; i>0; i--) {
00357     Z[k] += X[i];
00358     if (Z[k] < ZERO) {
00359       Z[k]  += RADIX;
00360       Z[--k] = MONE; }
00361     else
00362       Z[--k] = ZERO;
00363   }
00364 
00365   for (i=1; Z[i] == ZERO; i++) ;
00366   EZ = EZ - i + 1;
00367   for (k=1; i <= p+1; )
00368     Z[k++] = Z[i++];
00369   for (; k <= p; )
00370     Z[k++] = ZERO;
00371 
00372   return;
00373 }
00374 
00375 
00376 /* Add two multiple precision numbers. Set *z = *x + *y. x&y may overlap  */
00377 /* but not x&z or y&z. One guard digit is used. The error is less than    */
00378 /* one ulp. *x & *y are left unchanged.                                   */
00379 
00380 void __add(const mp_no *x, const mp_no *y, mp_no *z, int p) {
00381 
00382   int n;
00383 
00384   if      (X[0] == ZERO)     {__cpy(y,z,p);  return; }
00385   else if (Y[0] == ZERO)     {__cpy(x,z,p);  return; }
00386 
00387   if (X[0] == Y[0])   {
00388     if (__acr(x,y,p) > 0)      {add_magnitudes(x,y,z,p);  Z[0] = X[0]; }
00389     else                     {add_magnitudes(y,x,z,p);  Z[0] = Y[0]; }
00390   }
00391   else                       {
00392     if ((n=__acr(x,y,p)) == 1) {sub_magnitudes(x,y,z,p);  Z[0] = X[0]; }
00393     else if (n == -1)        {sub_magnitudes(y,x,z,p);  Z[0] = Y[0]; }
00394     else                      Z[0] = ZERO;
00395   }
00396   return;
00397 }
00398 
00399 
00400 /* Subtract two multiple precision numbers. *z is set to *x - *y. x&y may */
00401 /* overlap but not x&z or y&z. One guard digit is used. The error is      */
00402 /* less than one ulp. *x & *y are left unchanged.                         */
00403 
00404 void __sub(const mp_no *x, const mp_no *y, mp_no *z, int p) {
00405 
00406   int n;
00407 
00408   if      (X[0] == ZERO)     {__cpy(y,z,p);  Z[0] = -Z[0];  return; }
00409   else if (Y[0] == ZERO)     {__cpy(x,z,p);                 return; }
00410 
00411   if (X[0] != Y[0])    {
00412     if (__acr(x,y,p) > 0)      {add_magnitudes(x,y,z,p);  Z[0] =  X[0]; }
00413     else                     {add_magnitudes(y,x,z,p);  Z[0] = -Y[0]; }
00414   }
00415   else                       {
00416     if ((n=__acr(x,y,p)) == 1) {sub_magnitudes(x,y,z,p);  Z[0] =  X[0]; }
00417     else if (n == -1)        {sub_magnitudes(y,x,z,p);  Z[0] = -Y[0]; }
00418     else                      Z[0] = ZERO;
00419   }
00420   return;
00421 }
00422 
00423 
00424 /* Multiply two multiple precision numbers. *z is set to *x * *y. x&y      */
00425 /* may overlap but not x&z or y&z. In case p=1,2,3 the exact result is     */
00426 /* truncated to p digits. In case p>3 the error is bounded by 1.001 ulp.   */
00427 /* *x & *y are left unchanged.                                             */
00428 
00429 void __mul(const mp_no *x, const mp_no *y, mp_no *z, int p) {
00430 
00431   int i, i1, i2, j, k, k2;
00432   double u;
00433 
00434                       /* Is z=0? */
00435   if (X[0]*Y[0]==ZERO)
00436      { Z[0]=ZERO;  return; }
00437 
00438                        /* Multiply, add and carry */
00439   k2 = (p<3) ? p+p : p+3;
00440   Z[k2]=ZERO;
00441   for (k=k2; k>1; ) {
00442     if (k > p)  {i1=k-p; i2=p+1; }
00443     else        {i1=1;   i2=k;   }
00444     for (i=i1,j=i2-1; i<i2; i++,j--)  Z[k] += X[i]*Y[j];
00445 
00446     u = (Z[k] + CUTTER)-CUTTER;
00447     if  (u > Z[k])  u -= RADIX;
00448     Z[k]  -= u;
00449     Z[--k] = u*RADIXI;
00450   }
00451 
00452                  /* Is there a carry beyond the most significant digit? */
00453   if (Z[1] == ZERO) {
00454     for (i=1; i<=p; i++)  Z[i]=Z[i+1];
00455     EZ = EX + EY - 1; }
00456   else
00457     EZ = EX + EY;
00458 
00459   Z[0] = X[0] * Y[0];
00460   return;
00461 }
00462 
00463 
00464 /* Invert a multiple precision number. Set *y = 1 / *x.                     */
00465 /* Relative error bound = 1.001*r**(1-p) for p=2, 1.063*r**(1-p) for p=3,   */
00466 /* 2.001*r**(1-p) for p>3.                                                  */
00467 /* *x=0 is not permissible. *x is left unchanged.                           */
00468 
00469 void __inv(const mp_no *x, mp_no *y, int p) {
00470   int i;
00471 #if 0
00472   int l;
00473 #endif
00474   double t;
00475   mp_no z,w;
00476   static const int np1[] = {0,0,0,0,1,2,2,2,2,3,3,3,3,3,3,3,3,3,
00477                             4,4,4,4,4,4,4,4,4,4,4,4,4,4,4};
00478   const mp_no mptwo = {1,{1.0,2.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
00479                          0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
00480                          0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
00481                          0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}};
00482 
00483   __cpy(x,&z,p);  z.e=0;  __mp_dbl(&z,&t,p);
00484   t=ONE/t;   __dbl_mp(t,y,p);    EY -= EX;
00485 
00486   for (i=0; i<np1[p]; i++) {
00487     __cpy(y,&w,p);
00488     __mul(x,&w,y,p);
00489     __sub(&mptwo,y,&z,p);
00490     __mul(&w,&z,y,p);
00491   }
00492   return;
00493 }
00494 
00495 
00496 /* Divide one multiple precision number by another.Set *z = *x / *y. *x & *y */
00497 /* are left unchanged. x&y may overlap but not x&z or y&z.                   */
00498 /* Relative error bound = 2.001*r**(1-p) for p=2, 2.063*r**(1-p) for p=3     */
00499 /* and 3.001*r**(1-p) for p>3. *y=0 is not permissible.                      */
00500 
00501 void __dvd(const mp_no *x, const mp_no *y, mp_no *z, int p) {
00502 
00503   mp_no w;
00504 
00505   if (X[0] == ZERO)    Z[0] = ZERO;
00506   else                {__inv(y,&w,p);   __mul(x,&w,z,p);}
00507   return;
00508 }