sysdeps/ieee754/dbl-64/mpa.h

   1 /*
   2  * IBM Accurate Mathematical Library
   3  * Written by International Business Machines Corp.
   4  * Copyright (C) 2001-2016 Free Software Foundation, Inc.
   5  *
   6  * This program is free software; you can redistribute it and/or modify
   7  * it under the terms of the GNU Lesser General Public License as published by
   8  * the Free Software Foundation; either version 2.1 of the License, or
   9  * (at your option) any later version.
  10  *
  11  * This program is distributed in the hope that it will be useful,
  12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14  * GNU Lesser General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU Lesser General Public License
  17  * along with this program; if not, see <http://www.gnu.org/licenses/>.
  18  */
  19
  20 /************************************************************************/
  21 /*  MODULE_NAME: mpa.h                                                  */
  22 /*                                                                      */
  23 /*  FUNCTIONS:                                                          */
  24 /*               mcr                                                    */
  25 /*               acr                                                    */
  26 /*               cpy                                                    */
  27 /*               mp_dbl                                                 */
  28 /*               dbl_mp                                                 */
  29 /*               add                                                    */
  30 /*               sub                                                    */
  31 /*               mul                                                    */
  32 /*               dvd                                                    */
  33 /*                                                                      */
  34 /* Arithmetic functions for multiple precision numbers.                 */
  35 /* Common types and definition                                          */
  36 /************************************************************************/
  37
  38 #include <mpa-arch.h>
  39
  40 /* The mp_no structure holds the details of a multi-precision floating point
  41    number.
  42
  43    - The radix of the number (R) is 2 ^ 24.
  44
  45    - E: The exponent of the number.
  46
  47    - D[0]: The sign (-1, 1) or 0 if the value is 0.  In the latter case, the
  48      values of the remaining members of the structure are ignored.
  49
  50    - D[1] - D[p]: The mantissa of the number where:
  51
  52         0 <= D[i] < R and
  53         P is the precision of the number and 1 <= p <= 32
  54
  55      D[p+1] ... D[39] have no significance.
  56
  57    - The value of the number is:
  58
  59         D[1] * R ^ (E - 1) + D[2] * R ^ (E - 2) ... D[p] * R ^ (E - p)
  60
  61    */
  62 typedef struct
  63 {
  64   int e;
  65   mantissa_t d[40];
  66 } mp_no;
  67
  68 typedef union
  69 {
  70   int i[2];
  71   double d;
  72 } number;
  73
  74 extern const mp_no __mpone;
  75 extern const mp_no __mptwo;
  76
  77 #define  X   x->d
  78 #define  Y   y->d
  79 #define  Z   z->d
  80 #define  EX  x->e
  81 #define  EY  y->e
  82 #define  EZ  z->e
  83
  84 #ifndef RADIXI
  85 # define  RADIXI    0x1.0p-24           /* 2^-24   */
  86 #endif
  87
  88 #ifndef TWO52
  89 # define  TWO52     0x1.0p52            /* 2^52    */
  90 #endif
  91
  92 #define  TWO5      TWOPOW (5)           /* 2^5     */
  93 #define  TWO8      TWOPOW (8)           /* 2^52    */
  94 #define  TWO10     TWOPOW (10)          /* 2^10    */
  95 #define  TWO18     TWOPOW (18)          /* 2^18    */
  96 #define  TWO19     TWOPOW (19)          /* 2^19    */
  97 #define  TWO23     TWOPOW (23)          /* 2^23    */
  98
  99 #define  HALFRAD   TWO23
 100
 101 #define  TWO57     0x1.0p57             /* 2^57    */
 102 #define  TWO71     0x1.0p71             /* 2^71    */
 103 #define  TWOM1032  0x1.0p-1032          /* 2^-1032 */
 104 #define  TWOM1022  0x1.0p-1022          /* 2^-1022 */
 105
 106 #define  HALF      0x1.0p-1             /* 1/2 */
 107 #define  MHALF     -0x1.0p-1            /* -1/2 */
 108
 109 int __acr (const mp_no *, const mp_no *, int);
 110 void __cpy (const mp_no *, mp_no *, int);
 111 void __mp_dbl (const mp_no *, double *, int);
 112 void __dbl_mp (double, mp_no *, int);
 113 void __add (const mp_no *, const mp_no *, mp_no *, int);
 114 void __sub (const mp_no *, const mp_no *, mp_no *, int);
 115 void __mul (const mp_no *, const mp_no *, mp_no *, int);
 116 void __sqr (const mp_no *, mp_no *, int);
 117 void __dvd (const mp_no *, const mp_no *, mp_no *, int);
 118
 119 extern void __mpatan (mp_no *, mp_no *, int);
 120 extern void __mpatan2 (mp_no *, mp_no *, mp_no *, int);
 121 extern void __mpsqrt (mp_no *, mp_no *, int);
 122 extern void __mpexp (mp_no *, mp_no *, int);
 123 extern void __c32 (mp_no *, mp_no *, mp_no *, int);
 124 extern int __mpranred (double, mp_no *, int);
 125
 126 /* Given a power POW, build a multiprecision number 2^POW.  */
 127 static inline void
 128 __pow_mp (int pow, mp_no *y, int p)
 129 {
 130   int i, rem;
 131
 132   /* The exponent is E such that E is a factor of 2^24.  The remainder (of the
 133      form 2^x) goes entirely into the first digit of the mantissa as it is
 134      always less than 2^24.  */
 135   EY = pow / 24;
 136   rem = pow - EY * 24;
 137   EY++;
 138
 139   /* If the remainder is negative, it means that POW was negative since
 140      |EY * 24| <= |pow|.  Adjust so that REM is positive and still less than
 141      24 because of which, the mantissa digit is less than 2^24.  */
 142   if (rem < 0)
 143     {
 144       EY--;
 145       rem += 24;
 146     }
 147   /* The sign of any 2^x is always positive.  */
 148   Y[0] = 1;
 149   Y[1] = 1 << rem;
 150
 151   /* Everything else is 0.  */
 152   for (i = 2; i <= p; i++)
 153     Y[i] = 0;
 154 }