libgcc-math/dbl-64/mpa.c

   1
   2 /*
   3  * IBM Accurate Mathematical Library
   4  * written by International Business Machines Corp.
   5  * Copyright (C) 2001 Free Software Foundation
   6  *
   7  * This program is free software; you can redistribute it and/or modify
   8  * it under the terms of the GNU Lesser General Public License as published by
   9  * the Free Software Foundation; either version 2.1 of the License, or
  10  * (at your option) any later version.
  11  *
  12  * This program is distributed in the hope that it will be useful,
  13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15  * GNU Lesser General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU Lesser General Public License
  18  * along with this program; if not, write to the Free Software
  19  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
  20  */
  21 /************************************************************************/
  22 /*  MODULE_NAME: mpa.c                                                  */
  23 /*                                                                      */
  24 /*  FUNCTIONS:                                                          */
  25 /*               mcr                                                    */
  26 /*               acr                                                    */
  27 /*               cr                                                     */
  28 /*               cpy                                                    */
  29 /*               cpymn                                                  */
  30 /*               norm                                                   */
  31 /*               denorm                                                 */
  32 /*               mp_dbl                                                 */
  33 /*               dbl_mp                                                 */
  34 /*               add_magnitudes                                         */
  35 /*               sub_magnitudes                                         */
  36 /*               add                                                    */
  37 /*               sub                                                    */
  38 /*               mul                                                    */
  39 /*               inv                                                    */
  40 /*               dvd                                                    */
  41 /*                                                                      */
  42 /* Arithmetic functions for multiple precision numbers.                 */
  43 /* Relative errors are bounded                                          */
  44 /************************************************************************/
  45
  46
  47 #include "endian.h"
  48 #include "mpa.h"
  49 #include "mpa2.h"
  50
  51 /* mcr() compares the sizes of the mantissas of two multiple precision  */
  52 /* numbers. Mantissas are compared regardless of the signs of the       */
  53 /* numbers, even if x->d[0] or y->d[0] are zero. Exponents are also     */
  54 /* disregarded.                                                         */
  55 static int mcr(const mp_no *x, const mp_no *y, int p) {
  56   int i;
  57   for (i=1; i<=p; i++) {
  58     if      (X[i] == Y[i])  continue;
  59     else if (X[i] >  Y[i])  return  1;
  60     else                    return -1; }
  61   return 0;
  62 }
  63
  64
  65
  66 /* acr() compares the absolute values of two multiple precision numbers */
  67 int __acr(const mp_no *x, const mp_no *y, int p) {
  68   int i;
  69
  70   if      (X[0] == ZERO) {
  71     if    (Y[0] == ZERO) i= 0;
  72     else                 i=-1;
  73   }
  74   else if (Y[0] == ZERO) i= 1;
  75   else {
  76     if      (EX >  EY)   i= 1;
  77     else if (EX <  EY)   i=-1;
  78     else                 i= mcr(x,y,p);
  79   }
  80
  81   return i;
  82 }
  83
  84
  85 /* cr90 compares the values of two multiple precision numbers           */
  86 int  __cr(const mp_no *x, const mp_no *y, int p) {
  87   int i;
  88
  89   if      (X[0] > Y[0])  i= 1;
  90   else if (X[0] < Y[0])  i=-1;
  91   else if (X[0] < ZERO ) i= __acr(y,x,p);
  92   else                   i= __acr(x,y,p);
  93
  94   return i;
  95 }
  96
  97
  98 /* Copy a multiple precision number. Set *y=*x. x=y is permissible.      */
  99 void __cpy(const mp_no *x, mp_no *y, int p) {
 100   int i;
 101
 102   EY = EX;
 103   for (i=0; i <= p; i++)    Y[i] = X[i];
 104
 105   return;
 106 }
 107
 108
 109 /* Copy a multiple precision number x of precision m into a */
 110 /* multiple precision number y of precision n. In case n>m, */
 111 /* the digits of y beyond the m'th are set to zero. In case */
 112 /* n<m, the digits of x beyond the n'th are ignored.        */
 113 /* x=y is permissible.                                      */
 114
 115 void __cpymn(const mp_no *x, int m, mp_no *y, int n) {
 116
 117   int i,k;
 118
 119   EY = EX;     k=MIN(m,n);
 120   for (i=0; i <= k; i++)    Y[i] = X[i];
 121   for (   ; i <= n; i++)    Y[i] = ZERO;
 122
 123   return;
 124 }
 125
 126 /* Convert a multiple precision number *x into a double precision */
 127 /* number *y, normalized case  (|x| >= 2**(-1022))) */
 128 static void norm(const mp_no *x, double *y, int p)
 129 {
 130   #define R  radixi.d
 131   int i;
 132 #if 0
 133   int k;
 134 #endif
 135   double a,c=c,u,v,z[5];
 136   if (p<5) {
 137     if      (p==1) c = X[1];
 138     else if (p==2) c = X[1] + R* X[2];
 139     else if (p==3) c = X[1] + R*(X[2]  +   R* X[3]);
 140     else if (p==4) c =(X[1] + R* X[2]) + R*R*(X[3] + R*X[4]);
 141   }
 142   else {
 143     for (a=ONE, z[1]=X[1]; z[1] < TWO23; )
 144         {a *= TWO;   z[1] *= TWO; }
 145
 146     for (i=2; i<5; i++) {
 147       z[i] = X[i]*a;
 148       u = (z[i] + CUTTER)-CUTTER;
 149       if  (u > z[i])  u -= RADIX;
 150       z[i] -= u;
 151       z[i-1] += u*RADIXI;
 152     }
 153
 154     u = (z[3] + TWO71) - TWO71;
 155     if (u > z[3])   u -= TWO19;
 156     v = z[3]-u;
 157
 158     if (v == TWO18) {
 159       if (z[4] == ZERO) {
 160         for (i=5; i <= p; i++) {
 161           if (X[i] == ZERO)   continue;
 162           else                {z[3] += ONE;   break; }
 163         }
 164       }
 165       else              z[3] += ONE;
 166     }
 167
 168     c = (z[1] + R *(z[2] + R * z[3]))/a;
 169   }
 170
 171   c *= X[0];
 172
 173   for (i=1; i<EX; i++)   c *= RADIX;
 174   for (i=1; i>EX; i--)   c *= RADIXI;
 175
 176   *y = c;
 177   return;
 178 #undef R
 179 }
 180
 181 /* Convert a multiple precision number *x into a double precision */
 182 /* number *y, denormalized case  (|x| < 2**(-1022))) */
 183 static void denorm(const mp_no *x, double *y, int p)
 184 {
 185   int i,k;
 186   double c,u,z[5];
 187 #if 0
 188   double a,v;
 189 #endif
 190
 191 #define R  radixi.d
 192   if (EX<-44 || (EX==-44 && X[1]<TWO5))
 193      { *y=ZERO; return; }
 194
 195   if      (p==1) {
 196     if      (EX==-42) {z[1]=X[1]+TWO10;  z[2]=ZERO;  z[3]=ZERO;  k=3;}
 197     else if (EX==-43) {z[1]=     TWO10;  z[2]=X[1];  z[3]=ZERO;  k=2;}
 198     else              {z[1]=     TWO10;  z[2]=ZERO;  z[3]=X[1];  k=1;}
 199   }
 200   else if (p==2) {
 201     if      (EX==-42) {z[1]=X[1]+TWO10;  z[2]=X[2];  z[3]=ZERO;  k=3;}
 202     else if (EX==-43) {z[1]=     TWO10;  z[2]=X[1];  z[3]=X[2];  k=2;}
 203     else              {z[1]=     TWO10;  z[2]=ZERO;  z[3]=X[1];  k=1;}
 204   }
 205   else {
 206     if      (EX==-42) {z[1]=X[1]+TWO10;  z[2]=X[2];  k=3;}
 207     else if (EX==-43) {z[1]=     TWO10;  z[2]=X[1];  k=2;}
 208     else              {z[1]=     TWO10;  z[2]=ZERO;  k=1;}
 209     z[3] = X[k];
 210   }
 211
 212   u = (z[3] + TWO57) - TWO57;
 213   if  (u > z[3])   u -= TWO5;
 214
 215   if (u==z[3]) {
 216     for (i=k+1; i <= p; i++) {
 217       if (X[i] == ZERO)   continue;
 218       else {z[3] += ONE;   break; }
 219     }
 220   }
 221
 222   c = X[0]*((z[1] + R*(z[2] + R*z[3])) - TWO10);
 223
 224   *y = c*TWOM1032;
 225   return;
 226
 227 #undef R
 228 }
 229
 230 /* Convert a multiple precision number *x into a double precision number *y. */
 231 /* The result is correctly rounded to the nearest/even. *x is left unchanged */
 232
 233 void __mp_dbl(const mp_no *x, double *y, int p) {
 234 #if 0
 235   int i,k;
 236   double a,c,u,v,z[5];
 237 #endif
 238
 239   if (X[0] == ZERO)  {*y = ZERO;  return; }
 240
 241   if      (EX> -42)                 norm(x,y,p);
 242   else if (EX==-42 && X[1]>=TWO10)  norm(x,y,p);
 243   else                              denorm(x,y,p);
 244 }
 245
 246
 247 /* dbl_mp() converts a double precision number x into a multiple precision  */
 248 /* number *y. If the precision p is too small the result is truncated. x is */
 249 /* left unchanged.                                                          */
 250
 251 void __dbl_mp(double x, mp_no *y, int p) {
 252
 253   int i,n;
 254   double u;
 255
 256   /* Sign */
 257   if      (x == ZERO)  {Y[0] = ZERO;  return; }
 258   else if (x >  ZERO)   Y[0] = ONE;
 259   else                 {Y[0] = MONE;  x=-x;   }
 260
 261   /* Exponent */
 262   for (EY=ONE; x >= RADIX; EY += ONE)   x *= RADIXI;
 263   for (      ; x <  ONE;   EY -= ONE)   x *= RADIX;
 264
 265   /* Digits */
 266   n=MIN(p,4);
 267   for (i=1; i<=n; i++) {
 268     u = (x + TWO52) - TWO52;
 269     if (u>x)   u -= ONE;
 270     Y[i] = u;     x -= u;    x *= RADIX; }
 271   for (   ; i<=p; i++)     Y[i] = ZERO;
 272   return;
 273 }
 274
 275
 276 /*  add_magnitudes() adds the magnitudes of *x & *y assuming that           */
 277 /*  abs(*x) >= abs(*y) > 0.                                                 */
 278 /* The sign of the sum *z is undefined. x&y may overlap but not x&z or y&z. */
 279 /* No guard digit is used. The result equals the exact sum, truncated.      */
 280 /* *x & *y are left unchanged.                                              */
 281
 282 static void add_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
 283
 284   int i,j,k;
 285
 286   EZ = EX;
 287
 288   i=p;    j=p+ EY - EX;    k=p+1;
 289
 290   if (j<1)
 291      {__cpy(x,z,p);  return; }
 292   else   Z[k] = ZERO;
 293
 294   for (; j>0; i--,j--) {
 295     Z[k] += X[i] + Y[j];
 296     if (Z[k] >= RADIX) {
 297       Z[k]  -= RADIX;
 298       Z[--k] = ONE; }
 299     else
 300       Z[--k] = ZERO;
 301   }
 302
 303   for (; i>0; i--) {
 304     Z[k] += X[i];
 305     if (Z[k] >= RADIX) {
 306       Z[k]  -= RADIX;
 307       Z[--k] = ONE; }
 308     else
 309       Z[--k] = ZERO;
 310   }
 311
 312   if (Z[1] == ZERO) {
 313     for (i=1; i<=p; i++)    Z[i] = Z[i+1]; }
 314   else   EZ += ONE;
 315 }
 316
 317
 318 /*  sub_magnitudes() subtracts the magnitudes of *x & *y assuming that      */
 319 /*  abs(*x) > abs(*y) > 0.                                                  */
 320 /* The sign of the difference *z is undefined. x&y may overlap but not x&z  */
 321 /* or y&z. One guard digit is used. The error is less than one ulp.         */
 322 /* *x & *y are left unchanged.                                              */
 323
 324 static void sub_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
 325
 326   int i,j,k;
 327
 328   EZ = EX;
 329
 330   if (EX == EY) {
 331     i=j=k=p;
 332     Z[k] = Z[k+1] = ZERO; }
 333   else {
 334     j= EX - EY;
 335     if (j > p)  {__cpy(x,z,p);  return; }
 336     else {
 337       i=p;   j=p+1-j;   k=p;
 338       if (Y[j] > ZERO) {
 339         Z[k+1] = RADIX - Y[j--];
 340         Z[k]   = MONE; }
 341       else {
 342         Z[k+1] = ZERO;
 343         Z[k]   = ZERO;   j--;}
 344     }
 345   }
 346
 347   for (; j>0; i--,j--) {
 348     Z[k] += (X[i] - Y[j]);
 349     if (Z[k] < ZERO) {
 350       Z[k]  += RADIX;
 351       Z[--k] = MONE; }
 352     else
 353       Z[--k] = ZERO;
 354   }
 355
 356   for (; i>0; i--) {
 357     Z[k] += X[i];
 358     if (Z[k] < ZERO) {
 359       Z[k]  += RADIX;
 360       Z[--k] = MONE; }
 361     else
 362       Z[--k] = ZERO;
 363   }
 364
 365   for (i=1; Z[i] == ZERO; i++) ;
 366   EZ = EZ - i + 1;
 367   for (k=1; i <= p+1; )
 368     Z[k++] = Z[i++];
 369   for (; k <= p; )
 370     Z[k++] = ZERO;
 371
 372   return;
 373 }
 374
 375
 376 /* Add two multiple precision numbers. Set *z = *x + *y. x&y may overlap  */
 377 /* but not x&z or y&z. One guard digit is used. The error is less than    */
 378 /* one ulp. *x & *y are left unchanged.                                   */
 379
 380 void __add(const mp_no *x, const mp_no *y, mp_no *z, int p) {
 381
 382   int n;
 383
 384   if      (X[0] == ZERO)     {__cpy(y,z,p);  return; }
 385   else if (Y[0] == ZERO)     {__cpy(x,z,p);  return; }
 386
 387   if (X[0] == Y[0])   {
 388     if (__acr(x,y,p) > 0)      {add_magnitudes(x,y,z,p);  Z[0] = X[0]; }
 389     else                     {add_magnitudes(y,x,z,p);  Z[0] = Y[0]; }
 390   }
 391   else                       {
 392     if ((n=__acr(x,y,p)) == 1) {sub_magnitudes(x,y,z,p);  Z[0] = X[0]; }
 393     else if (n == -1)        {sub_magnitudes(y,x,z,p);  Z[0] = Y[0]; }
 394     else                      Z[0] = ZERO;
 395   }
 396   return;
 397 }
 398
 399
 400 /* Subtract two multiple precision numbers. *z is set to *x - *y. x&y may */
 401 /* overlap but not x&z or y&z. One guard digit is used. The error is      */
 402 /* less than one ulp. *x & *y are left unchanged.                         */
 403
 404 void __sub(const mp_no *x, const mp_no *y, mp_no *z, int p) {
 405
 406   int n;
 407
 408   if      (X[0] == ZERO)     {__cpy(y,z,p);  Z[0] = -Z[0];  return; }
 409   else if (Y[0] == ZERO)     {__cpy(x,z,p);                 return; }
 410
 411   if (X[0] != Y[0])    {
 412     if (__acr(x,y,p) > 0)      {add_magnitudes(x,y,z,p);  Z[0] =  X[0]; }
 413     else                     {add_magnitudes(y,x,z,p);  Z[0] = -Y[0]; }
 414   }
 415   else                       {
 416     if ((n=__acr(x,y,p)) == 1) {sub_magnitudes(x,y,z,p);  Z[0] =  X[0]; }
 417     else if (n == -1)        {sub_magnitudes(y,x,z,p);  Z[0] = -Y[0]; }
 418     else                      Z[0] = ZERO;
 419   }
 420   return;
 421 }
 422
 423
 424 /* Multiply two multiple precision numbers. *z is set to *x * *y. x&y      */
 425 /* may overlap but not x&z or y&z. In case p=1,2,3 the exact result is     */
 426 /* truncated to p digits. In case p>3 the error is bounded by 1.001 ulp.   */
 427 /* *x & *y are left unchanged.                                             */
 428
 429 void __mul(const mp_no *x, const mp_no *y, mp_no *z, int p) {
 430
 431   int i, i1, i2, j, k, k2;
 432   double u;
 433
 434                       /* Is z=0? */
 435   if (X[0]*Y[0]==ZERO)
 436      { Z[0]=ZERO;  return; }
 437
 438                        /* Multiply, add and carry */
 439   k2 = (p<3) ? p+p : p+3;
 440   Z[k2]=ZERO;
 441   for (k=k2; k>1; ) {
 442     if (k > p)  {i1=k-p; i2=p+1; }
 443     else        {i1=1;   i2=k;   }
 444     for (i=i1,j=i2-1; i<i2; i++,j--)  Z[k] += X[i]*Y[j];
 445
 446     u = (Z[k] + CUTTER)-CUTTER;
 447     if  (u > Z[k])  u -= RADIX;
 448     Z[k]  -= u;
 449     Z[--k] = u*RADIXI;
 450   }
 451
 452                  /* Is there a carry beyond the most significant digit? */
 453   if (Z[1] == ZERO) {
 454     for (i=1; i<=p; i++)  Z[i]=Z[i+1];
 455     EZ = EX + EY - 1; }
 456   else
 457     EZ = EX + EY;
 458
 459   Z[0] = X[0] * Y[0];
 460   return;
 461 }
 462
 463
 464 /* Invert a multiple precision number. Set *y = 1 / *x.                     */
 465 /* Relative error bound = 1.001*r**(1-p) for p=2, 1.063*r**(1-p) for p=3,   */
 466 /* 2.001*r**(1-p) for p>3.                                                  */
 467 /* *x=0 is not permissible. *x is left unchanged.                           */
 468
 469 void __inv(const mp_no *x, mp_no *y, int p) {
 470   int i;
 471 #if 0
 472   int l;
 473 #endif
 474   double t;
 475   mp_no z,w;
 476   static const int np1[] = {0,0,0,0,1,2,2,2,2,3,3,3,3,3,3,3,3,3,
 477                             4,4,4,4,4,4,4,4,4,4,4,4,4,4,4};
 478   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,
 479                          0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
 480                          0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
 481                          0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}};
 482
 483   __cpy(x,&z,p);  z.e=0;  __mp_dbl(&z,&t,p);
 484   t=ONE/t;   __dbl_mp(t,y,p);    EY -= EX;
 485
 486   for (i=0; i<np1[p]; i++) {
 487     __cpy(y,&w,p);
 488     __mul(x,&w,y,p);
 489     __sub(&mptwo,y,&z,p);
 490     __mul(&w,&z,y,p);
 491   }
 492   return;
 493 }
 494
 495
 496 /* Divide one multiple precision number by another.Set *z = *x / *y. *x & *y */
 497 /* are left unchanged. x&y may overlap but not x&z or y&z.                   */
 498 /* Relative error bound = 2.001*r**(1-p) for p=2, 2.063*r**(1-p) for p=3     */
 499 /* and 3.001*r**(1-p) for p>3. *y=0 is not permissible.                      */
 500
 501 void __dvd(const mp_no *x, const mp_no *y, mp_no *z, int p) {
 502
 503   mp_no w;
 504
 505   if (X[0] == ZERO)    Z[0] = ZERO;
 506   else                {__inv(y,&w,p);   __mul(x,&w,z,p);}
 507   return;
 508 }