]>
Commit | Line | Data |
---|---|---|
e4d82761 UD |
1 | /* |
2 | * IBM Accurate Mathematical Library | |
c6c6dd48 | 3 | * Written by International Business Machines Corp. |
b168057a | 4 | * Copyright (C) 2001-2015 Free Software Foundation, Inc. |
e4d82761 UD |
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 | |
cc7375ce | 8 | * the Free Software Foundation; either version 2.1 of the License, or |
e4d82761 | 9 | * (at your option) any later version. |
50944bca | 10 | * |
e4d82761 UD |
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 | |
c6c6dd48 | 14 | * GNU Lesser General Public License for more details. |
e4d82761 UD |
15 | * |
16 | * You should have received a copy of the GNU Lesser General Public License | |
59ba27a6 | 17 | * along with this program; if not, see <http://www.gnu.org/licenses/>. |
e4d82761 | 18 | */ |
c6c6dd48 | 19 | |
e4d82761 UD |
20 | /************************************************************************/ |
21 | /* MODULE_NAME: mpa.h */ | |
22 | /* */ | |
23 | /* FUNCTIONS: */ | |
24 | /* mcr */ | |
25 | /* acr */ | |
e4d82761 | 26 | /* cpy */ |
e4d82761 UD |
27 | /* mp_dbl */ |
28 | /* dbl_mp */ | |
29 | /* add */ | |
30 | /* sub */ | |
31 | /* mul */ | |
e4d82761 UD |
32 | /* dvd */ |
33 | /* */ | |
34 | /* Arithmetic functions for multiple precision numbers. */ | |
35 | /* Common types and definition */ | |
36 | /************************************************************************/ | |
37 | ||
6f2e90e7 | 38 | #include <mpa-arch.h> |
e4d82761 | 39 | |
a9e48ab4 SP |
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 | */ | |
6420d207 SP |
62 | typedef struct |
63 | { | |
a9e48ab4 | 64 | int e; |
6f2e90e7 | 65 | mantissa_t d[40]; |
a9e48ab4 | 66 | } mp_no; |
e4d82761 | 67 | |
6420d207 SP |
68 | typedef union |
69 | { | |
70 | int i[2]; | |
71 | double d; | |
72 | } number; | |
e4d82761 | 73 | |
107a5bf0 JM |
74 | extern const mp_no __mpone; |
75 | extern const mp_no __mptwo; | |
b76eb5f0 | 76 | |
e4d82761 UD |
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 | ||
6f2e90e7 SP |
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 | |
f93a8d15 | 91 | |
6f2e90e7 SP |
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 | ||
e7906a47 SP |
99 | #define HALFRAD TWO23 |
100 | ||
f93a8d15 SP |
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 */ | |
f93a8d15 | 108 | |
6420d207 SP |
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); | |
d6752ccd | 116 | void __sqr (const mp_no *, mp_no *, int); |
6420d207 | 117 | void __dvd (const mp_no *, const mp_no *, mp_no *, int); |
1d052247 AJ |
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); | |
302913e1 | 122 | extern void __mpexp (mp_no *, mp_no *, int); |
1d052247 AJ |
123 | extern void __c32 (mp_no *, mp_no *, mp_no *, int); |
124 | extern int __mpranred (double, mp_no *, int); | |
caa99d06 SP |
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. */ | |
c2d94018 | 148 | Y[0] = 1; |
caa99d06 SP |
149 | Y[1] = 1 << rem; |
150 | ||
a64d7e0e | 151 | /* Everything else is 0. */ |
caa99d06 | 152 | for (i = 2; i <= p; i++) |
a64d7e0e | 153 | Y[i] = 0; |
caa99d06 | 154 | } |