diff options
Diffstat (limited to 'REORG.TODO/sysdeps/ieee754/dbl-64/mpa.h')
| -rw-r--r-- | REORG.TODO/sysdeps/ieee754/dbl-64/mpa.h | 154 |
1 files changed, 154 insertions, 0 deletions
diff --git a/REORG.TODO/sysdeps/ieee754/dbl-64/mpa.h b/REORG.TODO/sysdeps/ieee754/dbl-64/mpa.h new file mode 100644 index 0000000000..a665e6b8f7 --- /dev/null +++ b/REORG.TODO/sysdeps/ieee754/dbl-64/mpa.h @@ -0,0 +1,154 @@ +/* + * IBM Accurate Mathematical Library + * Written by International Business Machines Corp. + * Copyright (C) 2001-2017 Free Software Foundation, Inc. + * + * 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, see <http://www.gnu.org/licenses/>. + */ + +/************************************************************************/ +/* MODULE_NAME: mpa.h */ +/* */ +/* FUNCTIONS: */ +/* mcr */ +/* acr */ +/* cpy */ +/* mp_dbl */ +/* dbl_mp */ +/* add */ +/* sub */ +/* mul */ +/* dvd */ +/* */ +/* Arithmetic functions for multiple precision numbers. */ +/* Common types and definition */ +/************************************************************************/ + +#include <mpa-arch.h> + +/* The mp_no structure holds the details of a multi-precision floating point + number. + + - The radix of the number (R) is 2 ^ 24. + + - E: The exponent of the number. + + - D[0]: The sign (-1, 1) or 0 if the value is 0. In the latter case, the + values of the remaining members of the structure are ignored. + + - D[1] - D[p]: The mantissa of the number where: + + 0 <= D[i] < R and + P is the precision of the number and 1 <= p <= 32 + + D[p+1] ... D[39] have no significance. + + - The value of the number is: + + D[1] * R ^ (E - 1) + D[2] * R ^ (E - 2) ... D[p] * R ^ (E - p) + + */ +typedef struct +{ + int e; + mantissa_t d[40]; +} mp_no; + +typedef union +{ + int i[2]; + double d; +} number; + +extern const mp_no __mpone; +extern const mp_no __mptwo; + +#define X x->d +#define Y y->d +#define Z z->d +#define EX x->e +#define EY y->e +#define EZ z->e + +#ifndef RADIXI +# define RADIXI 0x1.0p-24 /* 2^-24 */ +#endif + +#ifndef TWO52 +# define TWO52 0x1.0p52 /* 2^52 */ +#endif + +#define TWO5 TWOPOW (5) /* 2^5 */ +#define TWO8 TWOPOW (8) /* 2^52 */ +#define TWO10 TWOPOW (10) /* 2^10 */ +#define TWO18 TWOPOW (18) /* 2^18 */ +#define TWO19 TWOPOW (19) /* 2^19 */ +#define TWO23 TWOPOW (23) /* 2^23 */ + +#define HALFRAD TWO23 + +#define TWO57 0x1.0p57 /* 2^57 */ +#define TWO71 0x1.0p71 /* 2^71 */ +#define TWOM1032 0x1.0p-1032 /* 2^-1032 */ +#define TWOM1022 0x1.0p-1022 /* 2^-1022 */ + +#define HALF 0x1.0p-1 /* 1/2 */ +#define MHALF -0x1.0p-1 /* -1/2 */ + +int __acr (const mp_no *, const mp_no *, int); +void __cpy (const mp_no *, mp_no *, int); +void __mp_dbl (const mp_no *, double *, int); +void __dbl_mp (double, mp_no *, int); +void __add (const mp_no *, const mp_no *, mp_no *, int); +void __sub (const mp_no *, const mp_no *, mp_no *, int); +void __mul (const mp_no *, const mp_no *, mp_no *, int); +void __sqr (const mp_no *, mp_no *, int); +void __dvd (const mp_no *, const mp_no *, mp_no *, int); + +extern void __mpatan (mp_no *, mp_no *, int); +extern void __mpatan2 (mp_no *, mp_no *, mp_no *, int); +extern void __mpsqrt (mp_no *, mp_no *, int); +extern void __mpexp (mp_no *, mp_no *, int); +extern void __c32 (mp_no *, mp_no *, mp_no *, int); +extern int __mpranred (double, mp_no *, int); + +/* Given a power POW, build a multiprecision number 2^POW. */ +static inline void +__pow_mp (int pow, mp_no *y, int p) +{ + int i, rem; + + /* The exponent is E such that E is a factor of 2^24. The remainder (of the + form 2^x) goes entirely into the first digit of the mantissa as it is + always less than 2^24. */ + EY = pow / 24; + rem = pow - EY * 24; + EY++; + + /* If the remainder is negative, it means that POW was negative since + |EY * 24| <= |pow|. Adjust so that REM is positive and still less than + 24 because of which, the mantissa digit is less than 2^24. */ + if (rem < 0) + { + EY--; + rem += 24; + } + /* The sign of any 2^x is always positive. */ + Y[0] = 1; + Y[1] = 1 << rem; + + /* Everything else is 0. */ + for (i = 2; i <= p; i++) + Y[i] = 0; +} |