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, 2006 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   long i;
00057   long p2 = p;
00058   for (i=1; i<=p2; i++) {
00059     if      (X[i] == Y[i])  continue;
00060     else if (X[i] >  Y[i])  return  1;
00061     else                    return -1; }
00062   return 0;
00063 }
00064 
00065 
00066 
00067 /* acr() compares the absolute values of two multiple precision numbers */
00068 int __acr(const mp_no *x, const mp_no *y, int p) {
00069   long i;
00070 
00071   if      (X[0] == ZERO) {
00072     if    (Y[0] == ZERO) i= 0;
00073     else                 i=-1;
00074   }
00075   else if (Y[0] == ZERO) i= 1;
00076   else {
00077     if      (EX >  EY)   i= 1;
00078     else if (EX <  EY)   i=-1;
00079     else                 i= mcr(x,y,p);
00080   }
00081 
00082   return i;
00083 }
00084 
00085 
00086 /* cr90 compares the values of two multiple precision numbers           */
00087 int  __cr(const mp_no *x, const mp_no *y, int p) {
00088   int i;
00089 
00090   if      (X[0] > Y[0])  i= 1;
00091   else if (X[0] < Y[0])  i=-1;
00092   else if (X[0] < ZERO ) i= __acr(y,x,p);
00093   else                   i= __acr(x,y,p);
00094 
00095   return i;
00096 }
00097 
00098 
00099 /* Copy a multiple precision number. Set *y=*x. x=y is permissible.      */
00100 void __cpy(const mp_no *x, mp_no *y, int p) {
00101   long i;
00102 
00103   EY = EX;
00104   for (i=0; i <= p; i++)    Y[i] = X[i];
00105 
00106   return;
00107 }
00108 
00109 
00110 /* Copy a multiple precision number x of precision m into a */
00111 /* multiple precision number y of precision n. In case n>m, */
00112 /* the digits of y beyond the m'th are set to zero. In case */
00113 /* n<m, the digits of x beyond the n'th are ignored.        */
00114 /* x=y is permissible.                                      */
00115 
00116 void __cpymn(const mp_no *x, int m, mp_no *y, int n) {
00117 
00118   long i,k;
00119   long n2 = n;
00120   long m2 = m;
00121 
00122   EY = EX;     k=MIN(m2,n2);
00123   for (i=0; i <= k; i++)    Y[i] = X[i];
00124   for (   ; i <= n2; i++)    Y[i] = ZERO;
00125 
00126   return;
00127 }
00128 
00129 /* Convert a multiple precision number *x into a double precision */
00130 /* number *y, normalized case  (|x| >= 2**(-1022))) */
00131 static void norm(const mp_no *x, double *y, int p)
00132 {
00133   #define R  radixi.d
00134   long i;
00135 #if 0
00136   int k;
00137 #endif
00138   double a,c,u,v,z[5];
00139   if (p<5) {
00140     if      (p==1) c = X[1];
00141     else if (p==2) c = X[1] + R* X[2];
00142     else if (p==3) c = X[1] + R*(X[2]  +   R* X[3]);
00143     else if (p==4) c =(X[1] + R* X[2]) + R*R*(X[3] + R*X[4]);
00144   }
00145   else {
00146     for (a=ONE, z[1]=X[1]; z[1] < TWO23; )
00147         {a *= TWO;   z[1] *= TWO; }
00148 
00149     for (i=2; i<5; i++) {
00150       z[i] = X[i]*a;
00151       u = (z[i] + CUTTER)-CUTTER;
00152       if  (u > z[i])  u -= RADIX;
00153       z[i] -= u;
00154       z[i-1] += u*RADIXI;
00155     }
00156 
00157     u = (z[3] + TWO71) - TWO71;
00158     if (u > z[3])   u -= TWO19;
00159     v = z[3]-u;
00160 
00161     if (v == TWO18) {
00162       if (z[4] == ZERO) {
00163         for (i=5; i <= p; i++) {
00164           if (X[i] == ZERO)   continue;
00165           else                {z[3] += ONE;   break; }
00166         }
00167       }
00168       else              z[3] += ONE;
00169     }
00170 
00171     c = (z[1] + R *(z[2] + R * z[3]))/a;
00172   }
00173 
00174   c *= X[0];
00175 
00176   for (i=1; i<EX; i++)   c *= RADIX;
00177   for (i=1; i>EX; i--)   c *= RADIXI;
00178 
00179   *y = c;
00180   return;
00181 #undef R
00182 }
00183 
00184 /* Convert a multiple precision number *x into a double precision */
00185 /* number *y, denormalized case  (|x| < 2**(-1022))) */
00186 static void denorm(const mp_no *x, double *y, int p)
00187 {
00188   long i,k;
00189   long p2 = p;
00190   double c,u,z[5];
00191 #if 0
00192   double a,v;
00193 #endif
00194 
00195 #define R  radixi.d
00196   if (EX<-44 || (EX==-44 && X[1]<TWO5))
00197      { *y=ZERO; return; }
00198 
00199   if      (p2==1) {
00200     if      (EX==-42) {z[1]=X[1]+TWO10;  z[2]=ZERO;  z[3]=ZERO;  k=3;}
00201     else if (EX==-43) {z[1]=     TWO10;  z[2]=X[1];  z[3]=ZERO;  k=2;}
00202     else              {z[1]=     TWO10;  z[2]=ZERO;  z[3]=X[1];  k=1;}
00203   }
00204   else if (p2==2) {
00205     if      (EX==-42) {z[1]=X[1]+TWO10;  z[2]=X[2];  z[3]=ZERO;  k=3;}
00206     else if (EX==-43) {z[1]=     TWO10;  z[2]=X[1];  z[3]=X[2];  k=2;}
00207     else              {z[1]=     TWO10;  z[2]=ZERO;  z[3]=X[1];  k=1;}
00208   }
00209   else {
00210     if      (EX==-42) {z[1]=X[1]+TWO10;  z[2]=X[2];  k=3;}
00211     else if (EX==-43) {z[1]=     TWO10;  z[2]=X[1];  k=2;}
00212     else              {z[1]=     TWO10;  z[2]=ZERO;  k=1;}
00213     z[3] = X[k];
00214   }
00215 
00216   u = (z[3] + TWO57) - TWO57;
00217   if  (u > z[3])   u -= TWO5;
00218 
00219   if (u==z[3]) {
00220     for (i=k+1; i <= p2; i++) {
00221       if (X[i] == ZERO)   continue;
00222       else {z[3] += ONE;   break; }
00223     }
00224   }
00225 
00226   c = X[0]*((z[1] + R*(z[2] + R*z[3])) - TWO10);
00227 
00228   *y = c*TWOM1032;
00229   return;
00230 
00231 #undef R
00232 }
00233 
00234 /* Convert a multiple precision number *x into a double precision number *y. */
00235 /* The result is correctly rounded to the nearest/even. *x is left unchanged */
00236 
00237 void __mp_dbl(const mp_no *x, double *y, int p) {
00238 #if 0
00239   int i,k;
00240   double a,c,u,v,z[5];
00241 #endif
00242 
00243   if (X[0] == ZERO)  {*y = ZERO;  return; }
00244 
00245   if      (EX> -42)                 norm(x,y,p);
00246   else if (EX==-42 && X[1]>=TWO10)  norm(x,y,p);
00247   else                              denorm(x,y,p);
00248 }
00249 
00250 
00251 /* dbl_mp() converts a double precision number x into a multiple precision  */
00252 /* number *y. If the precision p is too small the result is truncated. x is */
00253 /* left unchanged.                                                          */
00254 
00255 void __dbl_mp(double x, mp_no *y, int p) {
00256 
00257   long i,n;
00258   long p2 = p;
00259   double u;
00260 
00261   /* Sign */
00262   if      (x == ZERO)  {Y[0] = ZERO;  return; }
00263   else if (x >  ZERO)   Y[0] = ONE;
00264   else                 {Y[0] = MONE;  x=-x;   }
00265 
00266   /* Exponent */
00267   for (EY=ONE; x >= RADIX; EY += ONE)   x *= RADIXI;
00268   for (      ; x <  ONE;   EY -= ONE)   x *= RADIX;
00269 
00270   /* Digits */
00271   n=MIN(p2,4);
00272   for (i=1; i<=n; i++) {
00273     u = (x + TWO52) - TWO52;
00274     if (u>x)   u -= ONE;
00275     Y[i] = u;     x -= u;    x *= RADIX; }
00276   for (   ; i<=p2; i++)     Y[i] = ZERO;
00277   return;
00278 }
00279 
00280 
00281 /*  add_magnitudes() adds the magnitudes of *x & *y assuming that           */
00282 /*  abs(*x) >= abs(*y) > 0.                                                 */
00283 /* The sign of the sum *z is undefined. x&y may overlap but not x&z or y&z. */
00284 /* No guard digit is used. The result equals the exact sum, truncated.      */
00285 /* *x & *y are left unchanged.                                              */
00286 
00287 static void add_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
00288 
00289   long i,j,k;
00290   long p2 = p;
00291 
00292   EZ = EX;
00293 
00294   i=p2;    j=p2+ EY - EX;    k=p2+1;
00295 
00296   if (j<1)
00297      {__cpy(x,z,p);  return; }
00298   else   Z[k] = ZERO;
00299 
00300   for (; j>0; i--,j--) {
00301     Z[k] += X[i] + Y[j];
00302     if (Z[k] >= RADIX) {
00303       Z[k]  -= RADIX;
00304       Z[--k] = ONE; }
00305     else
00306       Z[--k] = ZERO;
00307   }
00308 
00309   for (; i>0; i--) {
00310     Z[k] += X[i];
00311     if (Z[k] >= RADIX) {
00312       Z[k]  -= RADIX;
00313       Z[--k] = ONE; }
00314     else
00315       Z[--k] = ZERO;
00316   }
00317 
00318   if (Z[1] == ZERO) {
00319     for (i=1; i<=p2; i++)    Z[i] = Z[i+1]; }
00320   else   EZ += ONE;
00321 }
00322 
00323 
00324 /*  sub_magnitudes() subtracts the magnitudes of *x & *y assuming that      */
00325 /*  abs(*x) > abs(*y) > 0.                                                  */
00326 /* The sign of the difference *z is undefined. x&y may overlap but not x&z  */
00327 /* or y&z. One guard digit is used. The error is less than one ulp.         */
00328 /* *x & *y are left unchanged.                                              */
00329 
00330 static void sub_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
00331 
00332   long i,j,k;
00333   long p2 = p;
00334 
00335   EZ = EX;
00336 
00337   if (EX == EY) {
00338     i=j=k=p2;
00339     Z[k] = Z[k+1] = ZERO; }
00340   else {
00341     j= EX - EY;
00342     if (j > p2)  {__cpy(x,z,p);  return; }
00343     else {
00344       i=p2;   j=p2+1-j;   k=p2;
00345       if (Y[j] > ZERO) {
00346         Z[k+1] = RADIX - Y[j--];
00347         Z[k]   = MONE; }
00348       else {
00349         Z[k+1] = ZERO;
00350         Z[k]   = ZERO;   j--;}
00351     }
00352   }
00353 
00354   for (; j>0; i--,j--) {
00355     Z[k] += (X[i] - Y[j]);
00356     if (Z[k] < ZERO) {
00357       Z[k]  += RADIX;
00358       Z[--k] = MONE; }
00359     else
00360       Z[--k] = ZERO;
00361   }
00362 
00363   for (; i>0; i--) {
00364     Z[k] += X[i];
00365     if (Z[k] < ZERO) {
00366       Z[k]  += RADIX;
00367       Z[--k] = MONE; }
00368     else
00369       Z[--k] = ZERO;
00370   }
00371 
00372   for (i=1; Z[i] == ZERO; i++) ;
00373   EZ = EZ - i + 1;
00374   for (k=1; i <= p2+1; )
00375     Z[k++] = Z[i++];
00376   for (; k <= p2; )
00377     Z[k++] = ZERO;
00378 
00379   return;
00380 }
00381 
00382 
00383 /* Add two multiple precision numbers. Set *z = *x + *y. x&y may overlap  */
00384 /* but not x&z or y&z. One guard digit is used. The error is less than    */
00385 /* one ulp. *x & *y are left unchanged.                                   */
00386 
00387 void __add(const mp_no *x, const mp_no *y, mp_no *z, int p) {
00388 
00389   int n;
00390 
00391   if      (X[0] == ZERO)     {__cpy(y,z,p);  return; }
00392   else if (Y[0] == ZERO)     {__cpy(x,z,p);  return; }
00393 
00394   if (X[0] == Y[0])   {
00395     if (__acr(x,y,p) > 0)      {add_magnitudes(x,y,z,p);  Z[0] = X[0]; }
00396     else                     {add_magnitudes(y,x,z,p);  Z[0] = Y[0]; }
00397   }
00398   else                       {
00399     if ((n=__acr(x,y,p)) == 1) {sub_magnitudes(x,y,z,p);  Z[0] = X[0]; }
00400     else if (n == -1)        {sub_magnitudes(y,x,z,p);  Z[0] = Y[0]; }
00401     else                      Z[0] = ZERO;
00402   }
00403   return;
00404 }
00405 
00406 
00407 /* Subtract two multiple precision numbers. *z is set to *x - *y. x&y may */
00408 /* overlap but not x&z or y&z. One guard digit is used. The error is      */
00409 /* less than one ulp. *x & *y are left unchanged.                         */
00410 
00411 void __sub(const mp_no *x, const mp_no *y, mp_no *z, int p) {
00412 
00413   int n;
00414 
00415   if      (X[0] == ZERO)     {__cpy(y,z,p);  Z[0] = -Z[0];  return; }
00416   else if (Y[0] == ZERO)     {__cpy(x,z,p);                 return; }
00417 
00418   if (X[0] != Y[0])    {
00419     if (__acr(x,y,p) > 0)      {add_magnitudes(x,y,z,p);  Z[0] =  X[0]; }
00420     else                     {add_magnitudes(y,x,z,p);  Z[0] = -Y[0]; }
00421   }
00422   else                       {
00423     if ((n=__acr(x,y,p)) == 1) {sub_magnitudes(x,y,z,p);  Z[0] =  X[0]; }
00424     else if (n == -1)        {sub_magnitudes(y,x,z,p);  Z[0] = -Y[0]; }
00425     else                      Z[0] = ZERO;
00426   }
00427   return;
00428 }
00429 
00430 
00431 /* Multiply two multiple precision numbers. *z is set to *x * *y. x&y      */
00432 /* may overlap but not x&z or y&z. In case p=1,2,3 the exact result is     */
00433 /* truncated to p digits. In case p>3 the error is bounded by 1.001 ulp.   */
00434 /* *x & *y are left unchanged.                                             */
00435 
00436 void __mul(const mp_no *x, const mp_no *y, mp_no *z, int p) {
00437 
00438   long i, i1, i2, j, k, k2;
00439   long p2 = p;
00440   double u, zk, zk2;
00441 
00442                       /* Is z=0? */
00443   if (X[0]*Y[0]==ZERO)
00444      { Z[0]=ZERO;  return; }
00445 
00446                        /* Multiply, add and carry */
00447   k2 = (p2<3) ? p2+p2 : p2+3;
00448   zk = Z[k2]=ZERO;
00449   for (k=k2; k>1; ) {
00450     if (k > p2)  {i1=k-p2; i2=p2+1; }
00451     else        {i1=1;   i2=k;   }
00452 #if 1
00453     /* rearange this inner loop to allow the fmadd instructions to be
00454        independent and execute in parallel on processors that have
00455        dual symetrical FP pipelines.  */
00456     if (i1 < (i2-1))
00457     {
00458        /* make sure we have at least 2 iterations */
00459        if (((i2 - i1) & 1L) == 1L)
00460        {
00461                 /* Handle the odd iterations case.  */
00462               zk2 = x->d[i2-1]*y->d[i1];
00463        }
00464        else
00465               zk2 = zero.d;
00466        /* Do two multiply/adds per loop iteration, using independent
00467            accumulators; zk and zk2.  */
00468        for (i=i1,j=i2-1; i<i2-1; i+=2,j-=2) 
00469        {
00470               zk += x->d[i]*y->d[j];
00471               zk2 += x->d[i+1]*y->d[j-1];
00472        }
00473        zk += zk2; /* final sum.  */
00474     }
00475     else
00476     {
00477         /* Special case when iterations is 1.  */
00478        zk += x->d[i1]*y->d[i1];
00479     }
00480 #else
00481     /* The orginal code.  */
00482     for (i=i1,j=i2-1; i<i2; i++,j--)  zk += X[i]*Y[j];
00483 #endif
00484 
00485     u = (zk + CUTTER)-CUTTER;
00486     if  (u > zk)  u -= RADIX;
00487     Z[k]  = zk - u;
00488     zk = u*RADIXI;
00489     --k;
00490   }
00491   Z[k] = zk;
00492 
00493                  /* Is there a carry beyond the most significant digit? */
00494   if (Z[1] == ZERO) {
00495     for (i=1; i<=p2; i++)  Z[i]=Z[i+1];
00496     EZ = EX + EY - 1; }
00497   else
00498     EZ = EX + EY;
00499 
00500   Z[0] = X[0] * Y[0];
00501   return;
00502 }
00503 
00504 
00505 /* Invert a multiple precision number. Set *y = 1 / *x.                     */
00506 /* Relative error bound = 1.001*r**(1-p) for p=2, 1.063*r**(1-p) for p=3,   */
00507 /* 2.001*r**(1-p) for p>3.                                                  */
00508 /* *x=0 is not permissible. *x is left unchanged.                           */
00509 
00510 void __inv(const mp_no *x, mp_no *y, int p) {
00511   long i;
00512 #if 0
00513   int l;
00514 #endif
00515   double t;
00516   mp_no z,w;
00517   static const int np1[] = {0,0,0,0,1,2,2,2,2,3,3,3,3,3,3,3,3,3,
00518                             4,4,4,4,4,4,4,4,4,4,4,4,4,4,4};
00519   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,
00520                          0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
00521                          0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
00522                          0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}};
00523 
00524   __cpy(x,&z,p);  z.e=0;  __mp_dbl(&z,&t,p);
00525   t=ONE/t;   __dbl_mp(t,y,p);    EY -= EX;
00526 
00527   for (i=0; i<np1[p]; i++) {
00528     __cpy(y,&w,p);
00529     __mul(x,&w,y,p);
00530     __sub(&mptwo,y,&z,p);
00531     __mul(&w,&z,y,p);
00532   }
00533   return;
00534 }
00535 
00536 
00537 /* Divide one multiple precision number by another.Set *z = *x / *y. *x & *y */
00538 /* are left unchanged. x&y may overlap but not x&z or y&z.                   */
00539 /* Relative error bound = 2.001*r**(1-p) for p=2, 2.063*r**(1-p) for p=3     */
00540 /* and 3.001*r**(1-p) for p>3. *y=0 is not permissible.                      */
00541 
00542 void __dvd(const mp_no *x, const mp_no *y, mp_no *z, int p) {
00543 
00544   mp_no w;
00545 
00546   if (X[0] == ZERO)    Z[0] = ZERO;
00547   else                {__inv(y,&w,p);   __mul(x,&w,z,p);}
00548   return;
00549 }