1/*
2 * IBM Accurate Mathematical Library
3 * Written by International Business Machines Corp.
4 * Copyright (C) 2001-2018 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 */
62typedef struct
63{
64 int e;
65 mantissa_t d[40];
66} mp_no;
67
68typedef union
69{
70 int i[2];
71 double d;
72} number;
73
74extern const mp_no __mpone;
75extern 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
109int __acr (const mp_no *, const mp_no *, int);
110void __cpy (const mp_no *, mp_no *, int);
111void __mp_dbl (const mp_no *, double *, int);
112void __dbl_mp (double, mp_no *, int);
113void __add (const mp_no *, const mp_no *, mp_no *, int);
114void __sub (const mp_no *, const mp_no *, mp_no *, int);
115void __mul (const mp_no *, const mp_no *, mp_no *, int);
116void __sqr (const mp_no *, mp_no *, int);
117void __dvd (const mp_no *, const mp_no *, mp_no *, int);
118
119extern void __mpatan (mp_no *, mp_no *, int);
120extern void __mpatan2 (mp_no *, mp_no *, mp_no *, int);
121extern void __mpsqrt (mp_no *, mp_no *, int);
122extern void __mpexp (mp_no *, mp_no *, int);
123extern void __c32 (mp_no *, mp_no *, mp_no *, int);
124extern int __mpranred (double, mp_no *, int);
125
126/* Given a power POW, build a multiprecision number 2^POW. */
127static 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}
155