mpa.c

glibc 2.9,最新版的C语言库函数

/*
 * IBM Accurate Mathematical Library
 * written by International Business Machines Corp.
 * Copyright (C) 2001, 2006 Free Software Foundation
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published by
 * the Free Software Foundation; either version 2.1 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
 */
/************************************************************************/
/*  MODULE_NAME: mpa.c                                                  */
/*                                                                      */
/*  FUNCTIONS:                                                          */
/*               mcr                                                    */
/*               acr                                                    */
/*               cr                                                     */
/*               cpy                                                    */
/*               cpymn                                                  */
/*               norm                                                   */
/*               denorm                                                 */
/*               mp_dbl                                                 */
/*               dbl_mp                                                 */
/*               add_magnitudes                                         */
/*               sub_magnitudes                                         */
/*               add                                                    */
/*               sub                                                    */
/*               mul                                                    */
/*               inv                                                    */
/*               dvd                                                    */
/*                                                                      */
/* Arithmetic functions for multiple precision numbers.                 */
/* Relative errors are bounded                                          */
/************************************************************************/


#include "endian.h"
#include "mpa.h"
#include "mpa2.h"
#include <sys/param.h>	/* For MIN() */
/* mcr() compares the sizes of the mantissas of two multiple precision  */
/* numbers. Mantissas are compared regardless of the signs of the       */
/* numbers, even if x->d[0] or y->d[0] are zero. Exponents are also     */
/* disregarded.                                                         */
static int mcr(const mp_no *x, const mp_no *y, int p) {
  long i;
  long p2 = p;
  for (i=1; i<=p2; i++) {
    if      (X[i] == Y[i])  continue;
    else if (X[i] >  Y[i])  return  1;
    else                    return -1; }
  return 0;
}



/* acr() compares the absolute values of two multiple precision numbers */
int __acr(const mp_no *x, const mp_no *y, int p) {
  long i;

  if      (X[0] == ZERO) {
    if    (Y[0] == ZERO) i= 0;
    else                 i=-1;
  }
  else if (Y[0] == ZERO) i= 1;
  else {
    if      (EX >  EY)   i= 1;
    else if (EX <  EY)   i=-1;
    else                 i= mcr(x,y,p);
  }

  return i;
}


/* cr90 compares the values of two multiple precision numbers           */
int  __cr(const mp_no *x, const mp_no *y, int p) {
  int i;

  if      (X[0] > Y[0])  i= 1;
  else if (X[0] < Y[0])  i=-1;
  else if (X[0] < ZERO ) i= __acr(y,x,p);
  else                   i= __acr(x,y,p);

  return i;
}


/* Copy a multiple precision number. Set *y=*x. x=y is permissible.      */
void __cpy(const mp_no *x, mp_no *y, int p) {
  long i;

  EY = EX;
  for (i=0; i <= p; i++)    Y[i] = X[i];

  return;
}


/* Copy a multiple precision number x of precision m into a */
/* multiple precision number y of precision n. In case n>m, */
/* the digits of y beyond the m'th are set to zero. In case */
/* n<m, the digits of x beyond the n'th are ignored.        */
/* x=y is permissible.                                      */

void __cpymn(const mp_no *x, int m, mp_no *y, int n) {

  long i,k;
  long n2 = n;
  long m2 = m;

  EY = EX;     k=MIN(m2,n2);
  for (i=0; i <= k; i++)    Y[i] = X[i];
  for (   ; i <= n2; i++)    Y[i] = ZERO;

  return;
}

/* Convert a multiple precision number *x into a double precision */
/* number *y, normalized case  (|x| >= 2**(-1022))) */
static void norm(const mp_no *x, double *y, int p)
{
  #define R  radixi.d
  long i;
#if 0
  int k;
#endif
  double a,c,u,v,z[5];
  if (p<5) {
    if      (p==1) c = X[1];
    else if (p==2) c = X[1] + R* X[2];
    else if (p==3) c = X[1] + R*(X[2]  +   R* X[3]);
    else if (p==4) c =(X[1] + R* X[2]) + R*R*(X[3] + R*X[4]);
  }
  else {
    for (a=ONE, z[1]=X[1]; z[1] < TWO23; )
        {a *= TWO;   z[1] *= TWO; }

    for (i=2; i<5; i++) {
      z[i] = X[i]*a;
      u = (z[i] + CUTTER)-CUTTER;
      if  (u > z[i])  u -= RADIX;
      z[i] -= u;
      z[i-1] += u*RADIXI;
    }

    u = (z[3] + TWO71) - TWO71;
    if (u > z[3])   u -= TWO19;
    v = z[3]-u;

    if (v == TWO18) {
      if (z[4] == ZERO) {
        for (i=5; i <= p; i++) {
          if (X[i] == ZERO)   continue;
          else                {z[3] += ONE;   break; }
        }
      }
      else              z[3] += ONE;
    }

    c = (z[1] + R *(z[2] + R * z[3]))/a;
  }

  c *= X[0];

  for (i=1; i<EX; i++)   c *= RADIX;
  for (i=1; i>EX; i--)   c *= RADIXI;

  *y = c;
  return;
#undef R
}

/* Convert a multiple precision number *x into a double precision */
/* number *y, denormalized case  (|x| < 2**(-1022))) */
static void denorm(const mp_no *x, double *y, int p)
{
  long i,k;
  long p2 = p;
  double c,u,z[5];
#if 0
  double a,v;
#endif

#define R  radixi.d
  if (EX<-44 || (EX==-44 && X[1]<TWO5))
     { *y=ZERO; return; }

  if      (p2==1) {
    if      (EX==-42) {z[1]=X[1]+TWO10;  z[2]=ZERO;  z[3]=ZERO;  k=3;}
    else if (EX==-43) {z[1]=     TWO10;  z[2]=X[1];  z[3]=ZERO;  k=2;}
    else              {z[1]=     TWO10;  z[2]=ZERO;  z[3]=X[1];  k=1;}
  }
  else if (p2==2) {
    if      (EX==-42) {z[1]=X[1]+TWO10;  z[2]=X[2];  z[3]=ZERO;  k=3;}
    else if (EX==-43) {z[1]=     TWO10;  z[2]=X[1];  z[3]=X[2];  k=2;}
    else              {z[1]=     TWO10;  z[2]=ZERO;  z[3]=X[1];  k=1;}
  }
  else {
    if      (EX==-42) {z[1]=X[1]+TWO10;  z[2]=X[2];  k=3;}
    else if (EX==-43) {z[1]=     TWO10;  z[2]=X[1];  k=2;}
    else              {z[1]=     TWO10;  z[2]=ZERO;  k=1;}
    z[3] = X[k];
  }

  u = (z[3] + TWO57) - TWO57;
  if  (u > z[3])   u -= TWO5;

  if (u==z[3]) {
    for (i=k+1; i <= p2; i++) {
      if (X[i] == ZERO)   continue;
      else {z[3] += ONE;   break; }
    }
  }

  c = X[0]*((z[1] + R*(z[2] + R*z[3])) - TWO10);

  *y = c*TWOM1032;
  return;

#undef R
}

/* Convert a multiple precision number *x into a double precision number *y. */
/* The result is correctly rounded to the nearest/even. *x is left unchanged */

void __mp_dbl(const mp_no *x, double *y, int p) {
#if 0
  int i,k;
  double a,c,u,v,z[5];
#endif

  if (X[0] == ZERO)  {*y = ZERO;  return; }

  if      (EX> -42)                 norm(x,y,p);
  else if (EX==-42 && X[1]>=TWO10)  norm(x,y,p);
  else                              denorm(x,y,p);
}


/* dbl_mp() converts a double precision number x into a multiple precision  */
/* number *y. If the precision p is too small the result is truncated. x is */
/* left unchanged.                                                          */

void __dbl_mp(double x, mp_no *y, int p) {

  long i,n;
  long p2 = p;
  double u;

  /* Sign */
  if      (x == ZERO)  {Y[0] = ZERO;  return; }
  else if (x >  ZERO)   Y[0] = ONE;
  else                 {Y[0] = MONE;  x=-x;   }

  /* Exponent */
  for (EY=ONE; x >= RADIX; EY += ONE)   x *= RADIXI;
  for (      ; x <  ONE;   EY -= ONE)   x *= RADIX;

  /* Digits */
  n=MIN(p2,4);
  for (i=1; i<=n; i++) {
    u = (x + TWO52) - TWO52;
    if (u>x)   u -= ONE;
    Y[i] = u;     x -= u;    x *= RADIX; }