mp-math.c
1 //! @file mp-math.c
2 //! @author J. Marcel van der Veer
3 //!
4 //! @section Copyright
5 //!
6 //! This file is part of Algol68G - an Algol 68 compiler-interpreter.
7 //! Copyright 2001-2023 J. Marcel van der Veer [algol68g@xs4all.nl].
8 //!
9 //! @section License
10 //!
11 //! This program is free software; you can redistribute it and/or modify it
12 //! under the terms of the GNU General Public License as published by the
13 //! Free Software Foundation; either version 3 of the License, or
14 //! (at your option) any later version.
15 //!
16 //! This program is distributed in the hope that it will be useful, but
17 //! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18 //! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
19 //! more details. You should have received a copy of the GNU General Public
20 //! License along with this program. If not, see [http://www.gnu.org/licenses/].
21
22 //! @section Synopsis
23 //!
24 //! [LONG] LONG REAL math functions.
25
26 #include "a68g.h"
27 #include "a68g-genie.h"
28 #include "a68g-prelude.h"
29 #include "a68g-double.h"
30 #include "a68g-mp.h"
31
32 //! @brief Test on |"z"| > 0.001 for argument reduction in "sin" and "exp".
33
34 static inline BOOL_T eps_mp (MP_T * z, int digs)
35 {
36 if (MP_DIGIT (z, 1) == 0) {
37 return A68_FALSE;
38 } else if (MP_EXPONENT (z) > -1) {
39 return A68_TRUE;
40 } else if (MP_EXPONENT (z) < -1) {
41 return A68_FALSE;
42 } else {
43 // More or less optimised for LONG and default LONG LONG precisions.
44 return (BOOL_T) (digs <= 10 ? ABS (MP_DIGIT (z, 1)) > 0.01 * MP_RADIX : ABS (MP_DIGIT (z, 1)) > 0.001 * MP_RADIX);
45 }
46 }
47
48 //! @brief PROC (LONG REAL) LONG REAL sqrt
49
50 MP_T *sqrt_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
51 {
52 ADDR_T pop_sp = A68_SP;
53 if (MP_DIGIT (x, 1) == 0) {
54 A68_SP = pop_sp;
55 SET_MP_ZERO (z, digs);
56 return z;
57 }
58 if (MP_DIGIT (x, 1) < 0) {
59 A68_SP = pop_sp;
60 errno = EDOM;
61 return NaN_MP;
62 }
63 int gdigs = FUN_DIGITS (digs), hdigs;
64 BOOL_T reciprocal = A68_FALSE;
65 MP_T *z_g = nil_mp (p, gdigs);
66 MP_T *x_g = len_mp (p, x, digs, gdigs);
67 MP_T *tmp = nil_mp (p, gdigs);
68 // Scaling for small x; sqrt (x) = 1 / sqrt (1 / x).
69 if ((reciprocal = (BOOL_T) (MP_EXPONENT (x_g) < 0)) == A68_TRUE) {
70 (void) rec_mp (p, x_g, x_g, gdigs);
71 }
72 if (ABS (MP_EXPONENT (x_g)) >= 2) {
73 // For extreme arguments we want accurate results as well.
74 int expo = (int) MP_EXPONENT (x_g);
75 MP_EXPONENT (x_g) = (MP_T) (expo % 2);
76 (void) sqrt_mp (p, z_g, x_g, gdigs);
77 MP_EXPONENT (z_g) += (MP_T) (expo / 2);
78 } else {
79 // Argument is in range. Estimate the root as double.
80 #if (A68_LEVEL >= 3)
81 DOUBLE_T x_d = mp_to_double_real (p, x_g, gdigs);
82 (void) double_real_to_mp (p, z_g, sqrtq (x_d), gdigs);
83 #else
84 REAL_T x_d = mp_to_real (p, x_g, gdigs);
85 (void) real_to_mp (p, z_g, sqrt (x_d), gdigs);
86 #endif
87 // Newton's method: x<n+1> = (x<n> + a / x<n>) / 2.
88 int decimals = DOUBLE_ACCURACY;
89 do {
90 decimals <<= 1;
91 hdigs = MINIMUM (1 + decimals / LOG_MP_RADIX, gdigs);
92 (void) div_mp (p, tmp, x_g, z_g, hdigs);
93 (void) add_mp (p, tmp, z_g, tmp, hdigs);
94 (void) half_mp (p, z_g, tmp, hdigs);
95 } while (decimals < 2 * gdigs * LOG_MP_RADIX);
96 }
97 if (reciprocal) {
98 (void) rec_mp (p, z_g, z_g, digs);
99 }
100 (void) shorten_mp (p, z, digs, z_g, gdigs);
101 // Exit.
102 A68_SP = pop_sp;
103 return z;
104 }
105
106 //! @brief PROC (LONG REAL) LONG REAL curt
107
108 MP_T *curt_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
109 {
110 ADDR_T pop_sp = A68_SP;
111 if (MP_DIGIT (x, 1) == 0) {
112 A68_SP = pop_sp;
113 SET_MP_ZERO (z, digs);
114 return z;
115 }
116 BOOL_T change_sign = A68_FALSE;
117 if (MP_DIGIT (x, 1) < 0) {
118 change_sign = A68_TRUE;
119 MP_DIGIT (x, 1) = -MP_DIGIT (x, 1);
120 }
121 int gdigs = FUN_DIGITS (digs), hdigs;
122 BOOL_T reciprocal = A68_FALSE;
123 MP_T *z_g = nil_mp (p, gdigs);
124 MP_T *x_g = len_mp (p, x, digs, gdigs);
125 MP_T *tmp = nil_mp (p, gdigs);
126 // Scaling for small x; curt (x) = 1 / curt (1 / x).
127 if ((reciprocal = (BOOL_T) (MP_EXPONENT (x_g) < 0)) == A68_TRUE) {
128 (void) rec_mp (p, x_g, x_g, gdigs);
129 }
130 if (ABS (MP_EXPONENT (x_g)) >= 3) {
131 // For extreme arguments we want accurate results as well.
132 int expo = (int) MP_EXPONENT (x_g);
133 MP_EXPONENT (x_g) = (MP_T) (expo % 3);
134 (void) curt_mp (p, z_g, x_g, gdigs);
135 MP_EXPONENT (z_g) += (MP_T) (expo / 3);
136 } else {
137 // Argument is in range. Estimate the root as double.
138 int decimals;
139 #if (A68_LEVEL >= 3)
140 DOUBLE_T x_d = mp_to_double_real (p, x_g, gdigs);
141 (void) double_real_to_mp (p, z_g, cbrtq (x_d), gdigs);
142 #else
143 REAL_T x_d = mp_to_real (p, x_g, gdigs);
144 (void) real_to_mp (p, z_g, cbrt (x_d), gdigs);
145 #endif
146 // Newton's method: x<n+1> = (2 x<n> + a / x<n> ^ 2) / 3.
147 decimals = DOUBLE_ACCURACY;
148 do {
149 decimals <<= 1;
150 hdigs = MINIMUM (1 + decimals / LOG_MP_RADIX, gdigs);
151 (void) mul_mp (p, tmp, z_g, z_g, hdigs);
152 (void) div_mp (p, tmp, x_g, tmp, hdigs);
153 (void) add_mp (p, tmp, z_g, tmp, hdigs);
154 (void) add_mp (p, tmp, z_g, tmp, hdigs);
155 (void) div_mp_digit (p, z_g, tmp, (MP_T) 3, hdigs);
156 } while (decimals < gdigs * LOG_MP_RADIX);
157 }
158 if (reciprocal) {
159 (void) rec_mp (p, z_g, z_g, digs);
160 }
161 (void) shorten_mp (p, z, digs, z_g, gdigs);
162 // Exit.
163 A68_SP = pop_sp;
164 if (change_sign) {
165 MP_DIGIT (z, 1) = -MP_DIGIT (z, 1);
166 }
167 return z;
168 }
169
170 //! @brief PROC (LONG REAL) LONG REAL hypot
171
172 MP_T *hypot_mp (NODE_T * p, MP_T * z, MP_T * x, MP_T * y, int digs)
173 {
174 // sqrt (x^2 + y^2).
175 ADDR_T pop_sp = A68_SP;
176 MP_T *t = nil_mp (p, digs);
177 MP_T *u = nil_mp (p, digs);
178 MP_T *v = nil_mp (p, digs);
179 (void) move_mp (u, x, digs);
180 (void) move_mp (v, y, digs);
181 MP_DIGIT (u, 1) = ABS (MP_DIGIT (u, 1));
182 MP_DIGIT (v, 1) = ABS (MP_DIGIT (v, 1));
183 if (IS_ZERO_MP (u)) {
184 (void) move_mp (z, v, digs);
185 } else if (IS_ZERO_MP (v)) {
186 (void) move_mp (z, u, digs);
187 } else {
188 SET_MP_ONE (t, digs);
189 (void) sub_mp (p, z, u, v, digs);
190 if (MP_DIGIT (z, 1) > 0) {
191 (void) div_mp (p, z, v, u, digs);
192 (void) mul_mp (p, z, z, z, digs);
193 (void) add_mp (p, z, t, z, digs);
194 (void) sqrt_mp (p, z, z, digs);
195 (void) mul_mp (p, z, u, z, digs);
196 } else {
197 (void) div_mp (p, z, u, v, digs);
198 (void) mul_mp (p, z, z, z, digs);
199 (void) add_mp (p, z, t, z, digs);
200 (void) sqrt_mp (p, z, z, digs);
201 (void) mul_mp (p, z, v, z, digs);
202 }
203 }
204 A68_SP = pop_sp;
205 return z;
206 }
207
208 //! @brief PROC (LONG REAL) LONG REAL exp
209
210 MP_T *exp_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
211 {
212 // Argument is reduced by using exp (z / (2 ** n)) ** (2 ** n) = exp(z).
213 int m, n, gdigs = FUN_DIGITS (digs);
214 ADDR_T pop_sp = A68_SP;
215 BOOL_T iterate;
216 if (MP_DIGIT (x, 1) == 0) {
217 SET_MP_ONE (z, digs);
218 return z;
219 }
220 MP_T *x_g = len_mp (p, x, digs, gdigs);
221 MP_T *sum = nil_mp (p, gdigs);
222 MP_T *a68_pow = nil_mp (p, gdigs);
223 MP_T *fac = nil_mp (p, gdigs);
224 MP_T *tmp = nil_mp (p, gdigs);
225 m = 0;
226 // Scale x down.
227 while (eps_mp (x_g, gdigs)) {
228 m++;
229 (void) half_mp (p, x_g, x_g, gdigs);
230 }
231 // Calculate Taylor sum exp (z) = 1 + z / 1 ! + z ** 2 / 2 ! + ..
232 SET_MP_ONE (sum, gdigs);
233 (void) add_mp (p, sum, sum, x_g, gdigs);
234 (void) mul_mp (p, a68_pow, x_g, x_g, gdigs);
235 (void) half_mp (p, tmp, a68_pow, gdigs);
236 (void) add_mp (p, sum, sum, tmp, gdigs);
237 (void) mul_mp (p, a68_pow, a68_pow, x_g, gdigs);
238 (void) div_mp_digit (p, tmp, a68_pow, (MP_T) 6, gdigs);
239 (void) add_mp (p, sum, sum, tmp, gdigs);
240 (void) mul_mp (p, a68_pow, a68_pow, x_g, gdigs);
241 (void) div_mp_digit (p, tmp, a68_pow, (MP_T) 24, gdigs);
242 (void) add_mp (p, sum, sum, tmp, gdigs);
243 (void) mul_mp (p, a68_pow, a68_pow, x_g, gdigs);
244 (void) div_mp_digit (p, tmp, a68_pow, (MP_T) 120, gdigs);
245 (void) add_mp (p, sum, sum, tmp, gdigs);
246 (void) mul_mp (p, a68_pow, a68_pow, x_g, gdigs);
247 (void) div_mp_digit (p, tmp, a68_pow, (MP_T) 720, gdigs);
248 (void) add_mp (p, sum, sum, tmp, gdigs);
249 (void) mul_mp (p, a68_pow, a68_pow, x_g, gdigs);
250 (void) div_mp_digit (p, tmp, a68_pow, (MP_T) 5040, gdigs);
251 (void) add_mp (p, sum, sum, tmp, gdigs);
252 (void) mul_mp (p, a68_pow, a68_pow, x_g, gdigs);
253 (void) div_mp_digit (p, tmp, a68_pow, (MP_T) 40320, gdigs);
254 (void) add_mp (p, sum, sum, tmp, gdigs);
255 (void) mul_mp (p, a68_pow, a68_pow, x_g, gdigs);
256 (void) div_mp_digit (p, tmp, a68_pow, (MP_T) 362880, gdigs);
257 (void) add_mp (p, sum, sum, tmp, gdigs);
258 (void) mul_mp (p, a68_pow, a68_pow, x_g, gdigs);
259 (void) set_mp (fac, (MP_T) (MP_T) 3628800, 0, gdigs);
260 n = 10;
261 iterate = (BOOL_T) (MP_DIGIT (a68_pow, 1) != 0);
262 while (iterate) {
263 (void) div_mp (p, tmp, a68_pow, fac, gdigs);
264 if (MP_EXPONENT (tmp) <= (MP_EXPONENT (sum) - gdigs)) {
265 iterate = A68_FALSE;
266 } else {
267 (void) add_mp (p, sum, sum, tmp, gdigs);
268 (void) mul_mp (p, a68_pow, a68_pow, x_g, gdigs);
269 n++;
270 (void) mul_mp_digit (p, fac, fac, (MP_T) n, gdigs);
271 }
272 }
273 // Square exp (x) up.
274 while (m--) {
275 (void) mul_mp (p, sum, sum, sum, gdigs);
276 }
277 (void) shorten_mp (p, z, digs, sum, gdigs);
278 // Exit.
279 A68_SP = pop_sp;
280 return z;
281 }
282
283 //! @brief PROC (LONG REAL) LONG REAL exp (x) - 1
284 // assuming "x" to be close to 0.
285
286 MP_T *expm1_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
287 {
288 int n, gdigs = FUN_DIGITS (digs);
289 ADDR_T pop_sp = A68_SP;
290 BOOL_T iterate;
291 if (MP_DIGIT (x, 1) == 0) {
292 SET_MP_ONE (z, digs);
293 return z;
294 }
295 MP_T *x_g = len_mp (p, x, digs, gdigs);
296 MP_T *sum = nil_mp (p, gdigs);
297 MP_T *a68_pow = nil_mp (p, gdigs);
298 MP_T *fac = nil_mp (p, gdigs);
299 MP_T *tmp = nil_mp (p, gdigs);
300 // Calculate Taylor sum expm1 (z) = z / 1 ! + z ** 2 / 2 ! + ...
301 (void) add_mp (p, sum, sum, x_g, gdigs);
302 (void) mul_mp (p, a68_pow, x_g, x_g, gdigs);
303 (void) half_mp (p, tmp, a68_pow, gdigs);
304 (void) add_mp (p, sum, sum, tmp, gdigs);
305 (void) mul_mp (p, a68_pow, a68_pow, x_g, gdigs);
306 (void) div_mp_digit (p, tmp, a68_pow, (MP_T) 6, gdigs);
307 (void) add_mp (p, sum, sum, tmp, gdigs);
308 (void) mul_mp (p, a68_pow, a68_pow, x_g, gdigs);
309 (void) div_mp_digit (p, tmp, a68_pow, (MP_T) 24, gdigs);
310 (void) add_mp (p, sum, sum, tmp, gdigs);
311 (void) mul_mp (p, a68_pow, a68_pow, x_g, gdigs);
312 (void) div_mp_digit (p, tmp, a68_pow, (MP_T) 120, gdigs);
313 (void) add_mp (p, sum, sum, tmp, gdigs);
314 (void) mul_mp (p, a68_pow, a68_pow, x_g, gdigs);
315 (void) div_mp_digit (p, tmp, a68_pow, (MP_T) 720, gdigs);
316 (void) add_mp (p, sum, sum, tmp, gdigs);
317 (void) mul_mp (p, a68_pow, a68_pow, x_g, gdigs);
318 (void) div_mp_digit (p, tmp, a68_pow, (MP_T) 5040, gdigs);
319 (void) add_mp (p, sum, sum, tmp, gdigs);
320 (void) mul_mp (p, a68_pow, a68_pow, x_g, gdigs);
321 (void) div_mp_digit (p, tmp, a68_pow, (MP_T) 40320, gdigs);
322 (void) add_mp (p, sum, sum, tmp, gdigs);
323 (void) mul_mp (p, a68_pow, a68_pow, x_g, gdigs);
324 (void) div_mp_digit (p, tmp, a68_pow, (MP_T) 362880, gdigs);
325 (void) add_mp (p, sum, sum, tmp, gdigs);
326 (void) mul_mp (p, a68_pow, a68_pow, x_g, gdigs);
327 (void) set_mp (fac, (MP_T) (MP_T) 3628800, 0, gdigs);
328 n = 10;
329 iterate = (BOOL_T) (MP_DIGIT (a68_pow, 1) != 0);
330 while (iterate) {
331 (void) div_mp (p, tmp, a68_pow, fac, gdigs);
332 if (MP_EXPONENT (tmp) <= (MP_EXPONENT (sum) - gdigs)) {
333 iterate = A68_FALSE;
334 } else {
335 (void) add_mp (p, sum, sum, tmp, gdigs);
336 (void) mul_mp (p, a68_pow, a68_pow, x_g, gdigs);
337 n++;
338 (void) mul_mp_digit (p, fac, fac, (MP_T) n, gdigs);
339 }
340 }
341 (void) shorten_mp (p, z, digs, sum, gdigs);
342 // Exit.
343 A68_SP = pop_sp;
344 return z;
345 }
346
347 //! @brief ln scale
348
349 MP_T *mp_ln_scale (NODE_T * p, MP_T * z, int digs)
350 {
351 ADDR_T pop_sp = A68_SP;
352 int gdigs = FUN_DIGITS (digs);
353 MP_T *z_g = nil_mp (p, gdigs);
354 // First see if we can restore a previous calculation.
355 if (gdigs <= A68_MP (mp_ln_scale_size)) {
356 (void) move_mp (z_g, A68_MP (mp_ln_scale), gdigs);
357 } else {
358 // No luck with the kept value, we generate a longer one.
359 (void) set_mp (z_g, (MP_T) 1, 1, gdigs);
360 (void) ln_mp (p, z_g, z_g, gdigs);
361 A68_MP (mp_ln_scale) = (MP_T *) get_heap_space ((unt) SIZE_MP (gdigs));
362 (void) move_mp (A68_MP (mp_ln_scale), z_g, gdigs);
363 A68_MP (mp_ln_scale_size) = gdigs;
364 }
365 (void) shorten_mp (p, z, digs, z_g, gdigs);
366 A68_SP = pop_sp;
367 return z;
368 }
369
370 //! @brief ln 10
371
372 MP_T *mp_ln_10 (NODE_T * p, MP_T * z, int digs)
373 {
374 ADDR_T pop_sp = A68_SP;
375 int gdigs = FUN_DIGITS (digs);
376 MP_T *z_g = nil_mp (p, gdigs);
377 // First see if we can restore a previous calculation.
378 if (gdigs <= A68_MP (mp_ln_10_size)) {
379 (void) move_mp (z_g, A68_MP (mp_ln_10), gdigs);
380 } else {
381 // No luck with the kept value, we generate a longer one.
382 (void) set_mp (z_g, (MP_T) 10, 0, gdigs);
383 (void) ln_mp (p, z_g, z_g, gdigs);
384 A68_MP (mp_ln_10) = (MP_T *) get_heap_space ((unt) SIZE_MP (gdigs));
385 (void) move_mp (A68_MP (mp_ln_10), z_g, gdigs);
386 A68_MP (mp_ln_10_size) = gdigs;
387 }
388 (void) shorten_mp (p, z, digs, z_g, gdigs);
389 A68_SP = pop_sp;
390 return z;
391 }
392
393 //! @brief PROC (LONG REAL) LONG REAL ln
394
395 MP_T *ln_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
396 {
397 // Depending on the argument we choose either Taylor or Newton.
398 ADDR_T pop_sp = A68_SP;
399 int gdigs = FUN_DIGITS (digs);
400 BOOL_T negative, scale;
401 MP_T expo = 0;
402 if (MP_DIGIT (x, 1) <= 0) {
403 errno = EDOM;
404 return NaN_MP;
405 }
406 MP_T *x_g = len_mp (p, x, digs, gdigs);
407 MP_T *z_g = nil_mp (p, gdigs);
408 // We use ln (1 / x) = - ln (x).
409 negative = (BOOL_T) (MP_EXPONENT (x_g) < 0);
410 if (negative) {
411 (void) rec_mp (p, x_g, x_g, digs);
412 }
413 // We want correct results for extreme arguments. We scale when "x_g" exceeds
414 // "MP_RADIX ** +- 2", using ln (x * MP_RADIX ** n) = ln (x) + n * ln (MP_RADIX).
415 scale = (BOOL_T) (ABS (MP_EXPONENT (x_g)) >= 2);
416 if (scale) {
417 expo = MP_EXPONENT (x_g);
418 MP_EXPONENT (x_g) = (MP_T) 0;
419 }
420 if (MP_EXPONENT (x_g) == 0 && MP_DIGIT (x_g, 1) == 1 && MP_DIGIT (x_g, 2) == 0) {
421 // Taylor sum for x close to unity.
422 // ln (x) = (x - 1) - (x - 1) ** 2 / 2 + (x - 1) ** 3 / 3 - ...
423 // This is faster for small x and avoids cancellation.
424 int n = 2;
425 BOOL_T iterate;
426 MP_T *tmp = nil_mp (p, gdigs);
427 MP_T *a68_pow = nil_mp (p, gdigs);
428 (void) minus_one_mp (p, x_g, x_g, gdigs);
429 (void) mul_mp (p, a68_pow, x_g, x_g, gdigs);
430 (void) move_mp (z_g, x_g, gdigs);
431 iterate = (BOOL_T) (MP_DIGIT (a68_pow, 1) != 0);
432 while (iterate) {
433 (void) div_mp_digit (p, tmp, a68_pow, (MP_T) n, gdigs);
434 if (MP_EXPONENT (tmp) <= (MP_EXPONENT (z_g) - gdigs)) {
435 iterate = A68_FALSE;
436 } else {
437 MP_DIGIT (tmp, 1) = (EVEN (n) ? -MP_DIGIT (tmp, 1) : MP_DIGIT (tmp, 1));
438 (void) add_mp (p, z_g, z_g, tmp, gdigs);
439 (void) mul_mp (p, a68_pow, a68_pow, x_g, gdigs);
440 n++;
441 }
442 }
443 } else {
444 // Newton's method: x<n+1> = x<n> - 1 + a / exp (x<n>).
445 MP_T *tmp = nil_mp (p, gdigs);
446 // Construct an estimate.
447 #if (A68_LEVEL >= 3)
448 (void) double_real_to_mp (p, z_g, logq (mp_to_double_real (p, x_g, gdigs)), gdigs);
449 #else
450 (void) real_to_mp (p, z_g, log (mp_to_real (p, x_g, gdigs)), gdigs);
451 #endif
452 int decimals = DOUBLE_ACCURACY;
453 do {
454 int hdigs;
455 decimals <<= 1;
456 hdigs = MINIMUM (1 + decimals / LOG_MP_RADIX, gdigs);
457 (void) exp_mp (p, tmp, z_g, hdigs);
458 (void) div_mp (p, tmp, x_g, tmp, hdigs);
459 (void) minus_one_mp (p, z_g, z_g, hdigs);
460 (void) add_mp (p, z_g, z_g, tmp, hdigs);
461 } while (decimals < gdigs * LOG_MP_RADIX);
462 }
463 // Inverse scaling.
464 if (scale) {
465 // ln (x * MP_RADIX ** n) = ln (x) + n * ln (MP_RADIX).
466 MP_T *ln_base = nil_mp (p, gdigs);
467 (void) mp_ln_scale (p, ln_base, gdigs);
468 (void) mul_mp_digit (p, ln_base, ln_base, (MP_T) expo, gdigs);
469 (void) add_mp (p, z_g, z_g, ln_base, gdigs);
470 }
471 if (negative) {
472 MP_DIGIT (z_g, 1) = -MP_DIGIT (z_g, 1);
473 }
474 (void) shorten_mp (p, z, digs, z_g, gdigs);
475 // Exit.
476 A68_SP = pop_sp;
477 return z;
478 }
479
480 //! @brief PROC (LONG REAL) LONG REAL log
481
482 MP_T *log_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
483 {
484 ADDR_T pop_sp = A68_SP;
485 MP_T *ln_10 = nil_mp (p, digs);
486 if (ln_mp (p, z, x, digs) == NaN_MP) {
487 errno = EDOM;
488 return NaN_MP;
489 }
490 (void) mp_ln_10 (p, ln_10, digs);
491 (void) div_mp (p, z, z, ln_10, digs);
492 A68_SP = pop_sp;
493 return z;
494 }
495
496 //! @brief sinh ("z") and cosh ("z")
497
498 MP_T *hyp_mp (NODE_T * p, MP_T * sh, MP_T * ch, MP_T * z, int digs)
499 {
500 ADDR_T pop_sp = A68_SP;
501 MP_T *x_g = nil_mp (p, digs);
502 MP_T *y_g = nil_mp (p, digs);
503 MP_T *z_g = nil_mp (p, digs);
504 (void) move_mp (z_g, z, digs);
505 (void) exp_mp (p, x_g, z_g, digs);
506 (void) rec_mp (p, y_g, x_g, digs);
507 (void) add_mp (p, ch, x_g, y_g, digs);
508 // Avoid cancellation for sinh.
509 if ((MP_DIGIT (x_g, 1) == 1 && MP_DIGIT (x_g, 2) == 0) || (MP_DIGIT (y_g, 1) == 1 && MP_DIGIT (y_g, 2) == 0)) {
510 (void) expm1_mp (p, x_g, z_g, digs);
511 MP_DIGIT (z_g, 1) = -MP_DIGIT (z_g, 1);
512 (void) expm1_mp (p, y_g, z_g, digs);
513 }
514 (void) sub_mp (p, sh, x_g, y_g, digs);
515 (void) half_mp (p, sh, sh, digs);
516 (void) half_mp (p, ch, ch, digs);
517 A68_SP = pop_sp;
518 return z;
519 }
520
521 //! @brief PROC (LONG REAL) LONG REAL sinh
522
523 MP_T *sinh_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
524 {
525 ADDR_T pop_sp = A68_SP;
526 int gdigs = FUN_DIGITS (digs);
527 MP_T *x_g = len_mp (p, x, digs, gdigs);
528 MP_T *y_g = nil_mp (p, gdigs);
529 MP_T *z_g = nil_mp (p, gdigs);
530 (void) hyp_mp (p, z_g, y_g, x_g, gdigs);
531 (void) shorten_mp (p, z, digs, z_g, gdigs);
532 A68_SP = pop_sp;
533 return z;
534 }
535
536 //! @brief PROC (LONG REAL) LONG REAL asinh
537
538 MP_T *asinh_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
539 {
540 if (IS_ZERO_MP (x)) {
541 SET_MP_ZERO (z, digs);
542 return z;
543 } else {
544 ADDR_T pop_sp = A68_SP;
545 int gdigs;
546 if (MP_EXPONENT (x) >= -1) {
547 gdigs = FUN_DIGITS (digs);
548 } else {
549 // Extra precision when x^2+1 gets close to 1.
550 gdigs = 2 * FUN_DIGITS (digs);
551 }
552 MP_T *x_g = len_mp (p, x, digs, gdigs);
553 MP_T *y_g = nil_mp (p, gdigs);
554 MP_T *z_g = nil_mp (p, gdigs);
555 (void) mul_mp (p, z_g, x_g, x_g, gdigs);
556 SET_MP_ONE (y_g, gdigs);
557 (void) add_mp (p, y_g, z_g, y_g, gdigs);
558 (void) sqrt_mp (p, y_g, y_g, gdigs);
559 (void) add_mp (p, y_g, y_g, x_g, gdigs);
560 (void) ln_mp (p, z_g, y_g, gdigs);
561 if (IS_ZERO_MP (z_g)) {
562 (void) move_mp (z, x, digs);
563 } else {
564 (void) shorten_mp (p, z, digs, z_g, gdigs);
565 }
566 A68_SP = pop_sp;
567 return z;
568 }
569 }
570
571 //! @brief PROC (LONG REAL) LONG REAL cosh
572
573 MP_T *cosh_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
574 {
575 ADDR_T pop_sp = A68_SP;
576 int gdigs = FUN_DIGITS (digs);
577 MP_T *x_g = len_mp (p, x, digs, gdigs);
578 MP_T *y_g = nil_mp (p, gdigs);
579 MP_T *z_g = nil_mp (p, gdigs);
580 (void) hyp_mp (p, y_g, z_g, x_g, gdigs);
581 (void) shorten_mp (p, z, digs, z_g, gdigs);
582 A68_SP = pop_sp;
583 return z;
584 }
585
586 //! @brief PROC (LONG REAL) LONG REAL acosh
587
588 MP_T *acosh_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
589 {
590 ADDR_T pop_sp = A68_SP;
591 int gdigs;
592 if (MP_DIGIT (x, 1) == 1 && MP_DIGIT (x, 2) == 0) {
593 // Extra precision when x^2-1 gets close to 0.
594 gdigs = 2 * FUN_DIGITS (digs);
595 } else {
596 gdigs = FUN_DIGITS (digs);
597 }
598 MP_T *x_g = len_mp (p, x, digs, gdigs);
599 MP_T *y_g = nil_mp (p, gdigs);
600 MP_T *z_g = nil_mp (p, gdigs);
601 (void) mul_mp (p, z_g, x_g, x_g, gdigs);
602 SET_MP_ONE (y_g, gdigs);
603 (void) sub_mp (p, y_g, z_g, y_g, gdigs);
604 (void) sqrt_mp (p, y_g, y_g, gdigs);
605 (void) add_mp (p, y_g, y_g, x_g, gdigs);
606 (void) ln_mp (p, z_g, y_g, gdigs);
607 (void) shorten_mp (p, z, digs, z_g, gdigs);
608 A68_SP = pop_sp;
609 return z;
610 }
611
612 //! @brief PROC (LONG REAL) LONG REAL tanh
613
614 MP_T *tanh_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
615 {
616 ADDR_T pop_sp = A68_SP;
617 int gdigs = FUN_DIGITS (digs);
618 MP_T *x_g = len_mp (p, x, digs, gdigs);
619 MP_T *y_g = nil_mp (p, gdigs);
620 MP_T *z_g = nil_mp (p, gdigs);
621 (void) hyp_mp (p, y_g, z_g, x_g, gdigs);
622 (void) div_mp (p, z_g, y_g, z_g, gdigs);
623 (void) shorten_mp (p, z, digs, z_g, gdigs);
624 A68_SP = pop_sp;
625 return z;
626 }
627
628 //! @brief PROC (LONG REAL) LONG REAL atanh
629
630 MP_T *atanh_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
631 {
632 ADDR_T pop_sp = A68_SP;
633 int gdigs = FUN_DIGITS (digs);
634 MP_T *x_g = len_mp (p, x, digs, gdigs);
635 MP_T *y_g = nil_mp (p, gdigs);
636 MP_T *z_g = nil_mp (p, gdigs);
637 SET_MP_ONE (y_g, gdigs);
638 (void) add_mp (p, z_g, y_g, x_g, gdigs);
639 (void) sub_mp (p, y_g, y_g, x_g, gdigs);
640 (void) div_mp (p, y_g, z_g, y_g, gdigs);
641 (void) ln_mp (p, z_g, y_g, gdigs);
642 (void) half_mp (p, z_g, z_g, gdigs);
643 (void) shorten_mp (p, z, digs, z_g, gdigs);
644 A68_SP = pop_sp;
645 return z;
646 }
647
648 //! @brief PROC (LONG REAL) LONG REAL sin
649
650 MP_T *sin_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
651 {
652 // Use triple-angle relation to reduce argument.
653 ADDR_T pop_sp = A68_SP;
654 int m, n, gdigs = FUN_DIGITS (digs);
655 BOOL_T flip, negative, iterate, even;
656 // We will use "pi".
657 MP_T *pi = nil_mp (p, gdigs);
658 MP_T *tpi = nil_mp (p, gdigs);
659 MP_T *hpi = nil_mp (p, gdigs);
660 (void) mp_pi (p, pi, MP_PI, gdigs);
661 (void) mp_pi (p, tpi, MP_TWO_PI, gdigs);
662 (void) mp_pi (p, hpi, MP_HALF_PI, gdigs);
663 // Argument reduction (1): sin (x) = sin (x mod 2 pi).
664 MP_T *x_g = len_mp (p, x, digs, gdigs);
665 (void) mod_mp (p, x_g, x_g, tpi, gdigs);
666 // Argument reduction (2): sin (-x) = sin (x)
667 // sin (x) = - sin (x - pi); pi < x <= 2 pi
668 // sin (x) = sin (pi - x); pi / 2 < x <= pi
669 negative = (BOOL_T) (MP_DIGIT (x_g, 1) < 0);
670 if (negative) {
671 MP_DIGIT (x_g, 1) = -MP_DIGIT (x_g, 1);
672 }
673 MP_T *tmp = nil_mp (p, gdigs);
674 (void) sub_mp (p, tmp, x_g, pi, gdigs);
675 flip = (BOOL_T) (MP_DIGIT (tmp, 1) > 0);
676 if (flip) { // x > pi
677 (void) sub_mp (p, x_g, x_g, pi, gdigs);
678 }
679 (void) sub_mp (p, tmp, x_g, hpi, gdigs);
680 if (MP_DIGIT (tmp, 1) > 0) { // x > pi / 2
681 (void) sub_mp (p, x_g, pi, x_g, gdigs);
682 }
683 // Argument reduction (3): (follows from De Moivre's theorem)
684 // sin (3x) = sin (x) * (3 - 4 sin ^ 2 (x))
685 m = 0;
686 while (eps_mp (x_g, gdigs)) {
687 m++;
688 (void) div_mp_digit (p, x_g, x_g, (MP_T) 3, gdigs);
689 }
690 // Taylor sum.
691 MP_T *sqr = nil_mp (p, gdigs);
692 MP_T *a68_pow = nil_mp (p, gdigs);
693 MP_T *fac = nil_mp (p, gdigs);
694 MP_T *z_g = nil_mp (p, gdigs);
695 (void) mul_mp (p, sqr, x_g, x_g, gdigs);
696 (void) mul_mp (p, a68_pow, sqr, x_g, gdigs);
697 (void) move_mp (z_g, x_g, gdigs);
698 (void) div_mp_digit (p, tmp, a68_pow, (MP_T) 6, gdigs);
699 (void) sub_mp (p, z_g, z_g, tmp, gdigs);
700 (void) mul_mp (p, a68_pow, a68_pow, sqr, gdigs);
701 (void) div_mp_digit (p, tmp, a68_pow, (MP_T) 120, gdigs);
702 (void) add_mp (p, z_g, z_g, tmp, gdigs);
703 (void) mul_mp (p, a68_pow, a68_pow, sqr, gdigs);
704 (void) div_mp_digit (p, tmp, a68_pow, (MP_T) 5040, gdigs);
705 (void) sub_mp (p, z_g, z_g, tmp, gdigs);
706 (void) mul_mp (p, a68_pow, a68_pow, sqr, gdigs);
707 (void) set_mp (fac, (MP_T) 362880, 0, gdigs);
708 n = 9;
709 even = A68_TRUE;
710 iterate = (BOOL_T) (MP_DIGIT (a68_pow, 1) != 0);
711 while (iterate) {
712 (void) div_mp (p, tmp, a68_pow, fac, gdigs);
713 if (MP_EXPONENT (tmp) <= (MP_EXPONENT (z_g) - gdigs)) {
714 iterate = A68_FALSE;
715 } else {
716 if (even) {
717 (void) add_mp (p, z_g, z_g, tmp, gdigs);
718 even = A68_FALSE;
719 } else {
720 (void) sub_mp (p, z_g, z_g, tmp, gdigs);
721 even = A68_TRUE;
722 }
723 (void) mul_mp (p, a68_pow, a68_pow, sqr, gdigs);
724 (void) mul_mp_digit (p, fac, fac, (MP_T) (++n), gdigs);
725 (void) mul_mp_digit (p, fac, fac, (MP_T) (++n), gdigs);
726 }
727 }
728 // Inverse scaling using sin (3x) = sin (x) * (3 - 4 sin ** 2 (x)).
729 // Use existing mp's for intermediates.
730 (void) set_mp (fac, (MP_T) 3, 0, gdigs);
731 while (m--) {
732 (void) mul_mp (p, a68_pow, z_g, z_g, gdigs);
733 (void) mul_mp_digit (p, a68_pow, a68_pow, (MP_T) 4, gdigs);
734 (void) sub_mp (p, a68_pow, fac, a68_pow, gdigs);
735 (void) mul_mp (p, z_g, a68_pow, z_g, gdigs);
736 }
737 (void) shorten_mp (p, z, digs, z_g, gdigs);
738 if (negative ^ flip) {
739 MP_DIGIT (z, 1) = -MP_DIGIT (z, 1);
740 }
741 // Exit.
742 A68_SP = pop_sp;
743 return z;
744 }
745
746 //! @brief PROC (LONG REAL) LONG REAL cos
747
748 MP_T *cos_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
749 {
750 // Use cos (x) = sin (pi / 2 - x).
751 // Compute x mod 2 pi before subtracting to avoid cancellation.
752 ADDR_T pop_sp = A68_SP;
753 int gdigs = FUN_DIGITS (digs);
754 MP_T *hpi = nil_mp (p, gdigs);
755 MP_T *tpi = nil_mp (p, gdigs);
756 MP_T *x_g = len_mp (p, x, digs, gdigs);
757 MP_T *y = nil_mp (p, digs);
758 (void) mp_pi (p, hpi, MP_HALF_PI, gdigs);
759 (void) mp_pi (p, tpi, MP_TWO_PI, gdigs);
760 (void) mod_mp (p, x_g, x_g, tpi, gdigs);
761 (void) sub_mp (p, x_g, hpi, x_g, gdigs);
762 (void) shorten_mp (p, y, digs, x_g, gdigs);
763 (void) sin_mp (p, z, y, digs);
764 A68_SP = pop_sp;
765 return z;
766 }
767
768 //! @brief PROC (LONG REAL) LONG REAL tan
769
770 MP_T *tan_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
771 {
772 // Use tan (x) = sin (x) / sqrt (1 - sin ^ 2 (x)).
773 ADDR_T pop_sp = A68_SP;
774 int gdigs = FUN_DIGITS (digs);
775 BOOL_T negate;
776 MP_T *pi = nil_mp (p, gdigs);
777 MP_T *hpi = nil_mp (p, gdigs);
778 (void) mp_pi (p, pi, MP_PI, gdigs);
779 (void) mp_pi (p, hpi, MP_HALF_PI, gdigs);
780 MP_T *x_g = len_mp (p, x, digs, gdigs);
781 MP_T *y_g = nil_mp (p, gdigs);
782 MP_T *sns = nil_mp (p, digs);
783 MP_T *cns = nil_mp (p, digs);
784 // Argument mod pi.
785 (void) mod_mp (p, x_g, x_g, pi, gdigs);
786 if (MP_DIGIT (x_g, 1) >= 0) {
787 (void) sub_mp (p, y_g, x_g, hpi, gdigs);
788 negate = (BOOL_T) (MP_DIGIT (y_g, 1) > 0);
789 } else {
790 (void) add_mp (p, y_g, x_g, hpi, gdigs);
791 negate = (BOOL_T) (MP_DIGIT (y_g, 1) < 0);
792 }
793 (void) shorten_mp (p, x, digs, x_g, gdigs);
794 // tan(x) = sin(x) / sqrt (1 - sin ** 2 (x)).
795 (void) sin_mp (p, sns, x, digs);
796 (void) mul_mp (p, cns, sns, sns, digs);
797 (void) one_minus_mp (p, cns, cns, digs);
798 (void) sqrt_mp (p, cns, cns, digs);
799 if (div_mp (p, z, sns, cns, digs) == NaN_MP) {
800 errno = EDOM;
801 A68_SP = pop_sp;
802 return NaN_MP;
803 }
804 A68_SP = pop_sp;
805 if (negate) {
806 MP_DIGIT (z, 1) = -MP_DIGIT (z, 1);
807 }
808 return z;
809 }
810
811 //! @brief PROC (LONG REAL) LONG REAL arcsin
812
813 MP_T *asin_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
814 {
815 ADDR_T pop_sp = A68_SP;
816 int gdigs = FUN_DIGITS (digs);
817 MP_T *y = nil_mp (p, digs);
818 MP_T *x_g = len_mp (p, x, digs, gdigs);
819 MP_T *z_g = nil_mp (p, gdigs);
820 (void) mul_mp (p, z_g, x_g, x_g, gdigs);
821 (void) one_minus_mp (p, z_g, z_g, gdigs);
822 if (sqrt_mp (p, z_g, z_g, digs) == NaN_MP) {
823 errno = EDOM;
824 A68_SP = pop_sp;
825 return NaN_MP;
826 }
827 if (MP_DIGIT (z_g, 1) == 0) {
828 (void) mp_pi (p, z, MP_HALF_PI, digs);
829 MP_DIGIT (z, 1) = (MP_DIGIT (x_g, 1) >= 0 ? MP_DIGIT (z, 1) : -MP_DIGIT (z, 1));
830 A68_SP = pop_sp;
831 return z;
832 }
833 if (div_mp (p, x_g, x_g, z_g, gdigs) == NaN_MP) {
834 errno = EDOM;
835 A68_SP = pop_sp;
836 return NaN_MP;
837 }
838 (void) shorten_mp (p, y, digs, x_g, gdigs);
839 (void) atan_mp (p, z, y, digs);
840 A68_SP = pop_sp;
841 return z;
842 }
843
844 //! @brief PROC (LONG REAL) LONG REAL arccos
845
846 MP_T *acos_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
847 {
848 ADDR_T pop_sp = A68_SP;
849 int gdigs = FUN_DIGITS (digs);
850 BOOL_T negative = (BOOL_T) (MP_DIGIT (x, 1) < 0);
851 if (MP_DIGIT (x, 1) == 0) {
852 (void) mp_pi (p, z, MP_HALF_PI, digs);
853 A68_SP = pop_sp;
854 return z;
855 }
856 MP_T *y = nil_mp (p, digs);
857 MP_T *x_g = len_mp (p, x, digs, gdigs);
858 MP_T *z_g = nil_mp (p, gdigs);
859 (void) mul_mp (p, z_g, x_g, x_g, gdigs);
860 (void) one_minus_mp (p, z_g, z_g, gdigs);
861 if (sqrt_mp (p, z_g, z_g, digs) == NaN_MP) {
862 errno = EDOM;
863 A68_SP = pop_sp;
864 return NaN_MP;
865 }
866 if (div_mp (p, x_g, z_g, x_g, gdigs) == NaN_MP) {
867 errno = EDOM;
868 A68_SP = pop_sp;
869 return NaN_MP;
870 }
871 (void) shorten_mp (p, y, digs, x_g, gdigs);
872 (void) atan_mp (p, z, y, digs);
873 if (negative) {
874 (void) mp_pi (p, y, MP_PI, digs);
875 (void) add_mp (p, z, z, y, digs);
876 }
877 A68_SP = pop_sp;
878 return z;
879 }
880
881 //! @brief PROC (LONG REAL) LONG REAL arctan
882
883 MP_T *atan_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
884 {
885 // Depending on the argument we choose either Taylor or Newton.
886 ADDR_T pop_sp = A68_SP;
887 if (MP_DIGIT (x, 1) == 0) {
888 A68_SP = pop_sp;
889 SET_MP_ZERO (z, digs);
890 return z;
891 }
892 int gdigs = FUN_DIGITS (digs);
893 MP_T *x_g = len_mp (p, x, digs, gdigs);
894 MP_T *z_g = nil_mp (p, gdigs);
895 BOOL_T negative = (BOOL_T) (MP_DIGIT (x_g, 1) < 0);
896 if (negative) {
897 MP_DIGIT (x_g, 1) = -MP_DIGIT (x_g, 1);
898 }
899 // For larger arguments we use atan(x) = pi/2 - atan(1/x).
900 BOOL_T flip = (BOOL_T) (((MP_EXPONENT (x_g) > 0) || (MP_EXPONENT (x_g) == 0 && MP_DIGIT (x_g, 1) > 1)) && (MP_DIGIT (x_g, 1) != 0));
901 if (flip) {
902 (void) rec_mp (p, x_g, x_g, gdigs);
903 }
904 if (MP_EXPONENT (x_g) < -1 || (MP_EXPONENT (x_g) == -1 && MP_DIGIT (x_g, 1) < MP_RADIX / 100)) {
905 // Taylor sum for x close to zero.
906 // arctan (x) = x - x ** 3 / 3 + x ** 5 / 5 - x ** 7 / 7 + ..
907 // This is faster for small x and avoids cancellation
908 int n = 3;
909 BOOL_T iterate, even;
910 MP_T *tmp = nil_mp (p, gdigs);
911 MP_T *a68_pow = nil_mp (p, gdigs);
912 MP_T *sqr = nil_mp (p, gdigs);
913 (void) mul_mp (p, sqr, x_g, x_g, gdigs);
914 (void) mul_mp (p, a68_pow, sqr, x_g, gdigs);
915 (void) move_mp (z_g, x_g, gdigs);
916 even = A68_FALSE;
917 iterate = (BOOL_T) (MP_DIGIT (a68_pow, 1) != 0);
918 while (iterate) {
919 (void) div_mp_digit (p, tmp, a68_pow, (MP_T) n, gdigs);
920 if (MP_EXPONENT (tmp) <= (MP_EXPONENT (z_g) - gdigs)) {
921 iterate = A68_FALSE;
922 } else {
923 if (even) {
924 (void) add_mp (p, z_g, z_g, tmp, gdigs);
925 even = A68_FALSE;
926 } else {
927 (void) sub_mp (p, z_g, z_g, tmp, gdigs);
928 even = A68_TRUE;
929 }
930 (void) mul_mp (p, a68_pow, a68_pow, sqr, gdigs);
931 n += 2;
932 }
933 }
934 } else {
935 // Newton's method: x<n+1> = x<n> - cos (x<n>) * (sin (x<n>) - a cos (x<n>)).
936 int decimals, hdigs;
937 MP_T *tmp = nil_mp (p, gdigs);
938 MP_T *sns = nil_mp (p, gdigs);
939 MP_T *cns = nil_mp (p, gdigs);
940 // Construct an estimate.
941 #if (A68_LEVEL >= 3)
942 (void) double_real_to_mp (p, z_g, atanq (mp_to_double_real (p, x_g, gdigs)), gdigs);
943 #else
944 (void) real_to_mp (p, z_g, atan (mp_to_real (p, x_g, gdigs)), gdigs);
945 #endif
946 decimals = DOUBLE_ACCURACY;
947 do {
948 decimals <<= 1;
949 hdigs = MINIMUM (1 + decimals / LOG_MP_RADIX, gdigs);
950 (void) sin_mp (p, sns, z_g, hdigs);
951 (void) mul_mp (p, tmp, sns, sns, hdigs);
952 (void) one_minus_mp (p, tmp, tmp, hdigs);
953 (void) sqrt_mp (p, cns, tmp, hdigs);
954 (void) mul_mp (p, tmp, x_g, cns, hdigs);
955 (void) sub_mp (p, tmp, sns, tmp, hdigs);
956 (void) mul_mp (p, tmp, tmp, cns, hdigs);
957 (void) sub_mp (p, z_g, z_g, tmp, hdigs);
958 } while (decimals < gdigs * LOG_MP_RADIX);
959 }
960 if (flip) {
961 MP_T *hpi = nil_mp (p, gdigs);
962 (void) sub_mp (p, z_g, mp_pi (p, hpi, MP_HALF_PI, gdigs), z_g, gdigs);
963 }
964 (void) shorten_mp (p, z, digs, z_g, gdigs);
965 MP_DIGIT (z, 1) = (negative ? -MP_DIGIT (z, 1) : MP_DIGIT (z, 1));
966 // Exit.
967 A68_SP = pop_sp;
968 return z;
969 }
970
971 //! @brief PROC (LONG REAL, LONG REAL) LONG REAL atan2
972
973 MP_T *atan2_mp (NODE_T * p, MP_T * z, MP_T * x, MP_T * y, int digs)
974 {
975 ADDR_T pop_sp = A68_SP;
976 MP_T *t = nil_mp (p, digs);
977 if (MP_DIGIT (x, 1) == 0 && MP_DIGIT (y, 1) == 0) {
978 errno = EDOM;
979 A68_SP = pop_sp;
980 return NaN_MP;
981 } else {
982 BOOL_T flip = (BOOL_T) (MP_DIGIT (y, 1) < 0);
983 MP_DIGIT (y, 1) = ABS (MP_DIGIT (y, 1));
984 if (IS_ZERO_MP (x)) {
985 (void) mp_pi (p, z, MP_HALF_PI, digs);
986 } else {
987 BOOL_T flop = (BOOL_T) (MP_DIGIT (x, 1) <= 0);
988 MP_DIGIT (x, 1) = ABS (MP_DIGIT (x, 1));
989 (void) div_mp (p, z, y, x, digs);
990 (void) atan_mp (p, z, z, digs);
991 if (flop) {
992 (void) mp_pi (p, t, MP_PI, digs);
993 (void) sub_mp (p, z, t, z, digs);
994 }
995 }
996 if (flip) {
997 MP_DIGIT (z, 1) = -MP_DIGIT (z, 1);
998 }
999 }
1000 A68_SP = pop_sp;
1001 return z;
1002 }
1003
1004 //! @brief PROC (LONG REAL) LONG REAL csc
1005
1006 MP_T *csc_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1007 {
1008 sin_mp (p, z, x, digs);
1009 if (rec_mp (p, z, z, digs) == NaN_MP) {
1010 return NaN_MP;
1011 }
1012 return z;
1013 }
1014
1015 //! @brief PROC (LONG REAL) LONG REAL acsc
1016
1017 MP_T *acsc_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1018 {
1019 if (rec_mp (p, z, x, digs) == NaN_MP) {
1020 return NaN_MP;
1021 }
1022 return asin_mp (p, z, z, digs);
1023 }
1024
1025 //! @brief PROC (LONG REAL) LONG REAL sec
1026
1027 MP_T *sec_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1028 {
1029 cos_mp (p, z, x, digs);
1030 if (rec_mp (p, z, z, digs) == NaN_MP) {
1031 return NaN_MP;
1032 }
1033 return z;
1034 }
1035
1036 //! @brief PROC (LONG REAL) LONG REAL asec
1037
1038 MP_T *asec_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1039 {
1040 if (rec_mp (p, z, x, digs) == NaN_MP) {
1041 return NaN_MP;
1042 }
1043 return acos_mp (p, z, z, digs);
1044 }
1045
1046 //! @brief PROC (LONG REAL) LONG REAL cot
1047
1048 MP_T *cot_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1049 {
1050 // Use tan (x) = sin (x) / sqrt (1 - sin ^ 2 (x)).
1051 ADDR_T pop_sp = A68_SP;
1052 int gdigs = FUN_DIGITS (digs);
1053 BOOL_T negate;
1054 MP_T *pi = nil_mp (p, gdigs);
1055 MP_T *hpi = nil_mp (p, gdigs);
1056 (void) mp_pi (p, pi, MP_PI, gdigs);
1057 (void) mp_pi (p, hpi, MP_HALF_PI, gdigs);
1058 MP_T *x_g = len_mp (p, x, digs, gdigs);
1059 MP_T *y_g = nil_mp (p, gdigs);
1060 MP_T *sns = nil_mp (p, digs);
1061 MP_T *cns = nil_mp (p, digs);
1062 // Argument mod pi.
1063 (void) mod_mp (p, x_g, x_g, pi, gdigs);
1064 if (MP_DIGIT (x_g, 1) >= 0) {
1065 (void) sub_mp (p, y_g, x_g, hpi, gdigs);
1066 negate = (BOOL_T) (MP_DIGIT (y_g, 1) > 0);
1067 } else {
1068 (void) add_mp (p, y_g, x_g, hpi, gdigs);
1069 negate = (BOOL_T) (MP_DIGIT (y_g, 1) < 0);
1070 }
1071 (void) shorten_mp (p, x, digs, x_g, gdigs);
1072 // tan(x) = sin(x) / sqrt (1 - sin ** 2 (x)).
1073 (void) sin_mp (p, sns, x, digs);
1074 (void) mul_mp (p, cns, sns, sns, digs);
1075 (void) one_minus_mp (p, cns, cns, digs);
1076 (void) sqrt_mp (p, cns, cns, digs);
1077 if (div_mp (p, z, cns, sns, digs) == NaN_MP) {
1078 errno = EDOM;
1079 A68_SP = pop_sp;
1080 return NaN_MP;
1081 }
1082 A68_SP = pop_sp;
1083 if (negate) {
1084 MP_DIGIT (z, 1) = -MP_DIGIT (z, 1);
1085 }
1086 return z;
1087 }
1088
1089 //! @brief PROC (LONG REAL) LONG REAL arccot
1090
1091 MP_T *acot_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1092 {
1093 ADDR_T pop_sp = A68_SP;
1094 MP_T *f = nil_mp (p, digs);
1095 if (rec_mp (p, f, x, digs) == NaN_MP) {
1096 errno = EDOM;
1097 A68_SP = pop_sp;
1098 return NaN_MP;
1099 } else {
1100 (void) atan_mp (p, z, f, digs);
1101 A68_SP = pop_sp;
1102 return z;
1103 }
1104 }
1105
1106 //! @brief PROC (LONG REAL) LONG REAL sindg
1107
1108 MP_T *sindg_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1109 {
1110 ADDR_T pop_sp = A68_SP;
1111 MP_T *f = nil_mp (p, digs);
1112 MP_T *g = nil_mp (p, digs);
1113 (void) mp_pi (p, f, MP_PI_OVER_180, digs);
1114 (void) mul_mp (p, g, x, f, digs);
1115 (void) sin_mp (p, z, g, digs);
1116 A68_SP = pop_sp;
1117 return z;
1118 }
1119
1120 //! @brief PROC (LONG REAL) LONG REAL cosdg
1121
1122 MP_T *cosdg_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1123 {
1124 ADDR_T pop_sp = A68_SP;
1125 MP_T *f = nil_mp (p, digs);
1126 MP_T *g = nil_mp (p, digs);
1127 (void) mp_pi (p, f, MP_PI_OVER_180, digs);
1128 (void) mul_mp (p, g, x, f, digs);
1129 (void) cos_mp (p, z, g, digs);
1130 A68_SP = pop_sp;
1131 return z;
1132 }
1133
1134 //! @brief PROC (LONG REAL) LONG REAL tandg
1135
1136 MP_T *tandg_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1137 {
1138 ADDR_T pop_sp = A68_SP;
1139 MP_T *f = nil_mp (p, digs);
1140 MP_T *g = nil_mp (p, digs);
1141 (void) mp_pi (p, f, MP_PI_OVER_180, digs);
1142 (void) mul_mp (p, g, x, f, digs);
1143 if (tan_mp (p, z, g, digs) == NaN_MP) {
1144 errno = EDOM;
1145 A68_SP = pop_sp;
1146 return NaN_MP;
1147 } else {
1148 A68_SP = pop_sp;
1149 return z;
1150 }
1151 }
1152
1153 //! @brief PROC (LONG REAL) LONG REAL cotdg
1154
1155 MP_T *cotdg_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1156 {
1157 ADDR_T pop_sp = A68_SP;
1158 MP_T *f = nil_mp (p, digs);
1159 MP_T *g = nil_mp (p, digs);
1160 (void) mp_pi (p, f, MP_PI_OVER_180, digs);
1161 (void) mul_mp (p, g, x, f, digs);
1162 if (cot_mp (p, z, g, digs) == NaN_MP) {
1163 errno = EDOM;
1164 A68_SP = pop_sp;
1165 return NaN_MP;
1166 } else {
1167 A68_SP = pop_sp;
1168 return z;
1169 }
1170 }
1171
1172 //! @brief PROC (LONG REAL) LONG REAL arcsindg
1173
1174 MP_T *asindg_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1175 {
1176 ADDR_T pop_sp = A68_SP;
1177 MP_T *f = nil_mp (p, digs);
1178 MP_T *g = nil_mp (p, digs);
1179 (void) asin_mp (p, f, x, digs);
1180 (void) mp_pi (p, g, MP_180_OVER_PI, digs);
1181 (void) mul_mp (p, z, f, g, digs);
1182 A68_SP = pop_sp;
1183 return z;
1184 }
1185
1186 //! @brief PROC (LONG REAL) LONG REAL arccosdg
1187
1188 MP_T *acosdg_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1189 {
1190 ADDR_T pop_sp = A68_SP;
1191 MP_T *f = nil_mp (p, digs);
1192 MP_T *g = nil_mp (p, digs);
1193 (void) acos_mp (p, f, x, digs);
1194 (void) mp_pi (p, g, MP_180_OVER_PI, digs);
1195 (void) mul_mp (p, z, f, g, digs);
1196 A68_SP = pop_sp;
1197 return z;
1198 }
1199
1200 //! @brief PROC (LONG REAL) LONG REAL arctandg
1201
1202 MP_T *atandg_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1203 {
1204 ADDR_T pop_sp = A68_SP;
1205 MP_T *f = nil_mp (p, digs);
1206 MP_T *g = nil_mp (p, digs);
1207 (void) atan_mp (p, f, x, digs);
1208 (void) mp_pi (p, g, MP_180_OVER_PI, digs);
1209 (void) mul_mp (p, z, f, g, digs);
1210 A68_SP = pop_sp;
1211 return z;
1212 }
1213
1214 //! @brief PROC (LONG REAL) LONG REAL arccotdg
1215
1216 MP_T *acotdg_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1217 {
1218 ADDR_T pop_sp = A68_SP;
1219 MP_T *f = nil_mp (p, digs);
1220 MP_T *g = nil_mp (p, digs);
1221 if (acot_mp (p, f, x, digs) == NaN_MP) {
1222 errno = EDOM;
1223 A68_SP = pop_sp;
1224 return NaN_MP;
1225 } else {
1226 (void) mp_pi (p, g, MP_180_OVER_PI, digs);
1227 (void) mul_mp (p, z, f, g, digs);
1228 A68_SP = pop_sp;
1229 return z;
1230 }
1231 }
1232
1233 //! @brief PROC (LONG REAL, LONG REAL) LONG REAL atan2dg
1234
1235 MP_T *atan2dg_mp (NODE_T * p, MP_T * z, MP_T * x, MP_T * y, int digs)
1236 {
1237 ADDR_T pop_sp = A68_SP;
1238 MP_T *f = nil_mp (p, digs);
1239 MP_T *g = nil_mp (p, digs);
1240 if (atan2_mp (p, f, x, y, digs) == NaN_MP) {
1241 errno = EDOM;
1242 A68_SP = pop_sp;
1243 return NaN_MP;
1244 } else {
1245 (void) mp_pi (p, g, MP_180_OVER_PI, digs);
1246 (void) mul_mp (p, z, f, g, digs);
1247 A68_SP = pop_sp;
1248 return z;
1249 }
1250 }
1251
1252 //! @brief PROC (LONG REAL) LONG REAL sinpi
1253
1254 MP_T *sinpi_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1255 {
1256 ADDR_T pop_sp = A68_SP;
1257 MP_T *f = nil_mp (p, digs);
1258 MP_T *g = lit_mp (p, 1, 0, digs);
1259 (void) mod_mp (p, f, x, g, digs);
1260 if (IS_ZERO_MP (f)) {
1261 SET_MP_ZERO (z, digs);
1262 A68_SP = pop_sp;
1263 return z;
1264 }
1265 (void) mp_pi (p, f, MP_PI, digs);
1266 (void) mul_mp (p, g, x, f, digs);
1267 (void) sin_mp (p, z, g, digs);
1268 A68_SP = pop_sp;
1269 return z;
1270 }
1271
1272 //! @brief PROC (LONG REAL) LONG REAL acospi
1273
1274 MP_T *cospi_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1275 {
1276 ADDR_T pop_sp = A68_SP;
1277 MP_T *f = nil_mp (p, digs);
1278 MP_T *g = lit_mp (p, 1, 0, digs);
1279 (void) mod_mp (p, f, x, g, digs);
1280 abs_mp (p, f, f, digs);
1281 SET_MP_HALF (g, digs);
1282 (void) sub_mp (p, g, f, g, digs);
1283 if (IS_ZERO_MP (g)) {
1284 SET_MP_ZERO (z, digs);
1285 A68_SP = pop_sp;
1286 return z;
1287 }
1288 (void) mp_pi (p, f, MP_PI, digs);
1289 (void) mul_mp (p, g, x, f, digs);
1290 (void) cos_mp (p, z, g, digs);
1291 A68_SP = pop_sp;
1292 return z;
1293 }
1294
1295 //! @brief PROC (LONG REAL) LONG REAL tanpi
1296
1297 MP_T *tanpi_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1298 {
1299 ADDR_T pop_sp = A68_SP;
1300 MP_T *f = nil_mp (p, digs);
1301 MP_T *g = lit_mp (p, 1, 0, digs);
1302 MP_T *h = nil_mp (p, digs);
1303 MP_T *half = nil_mp (p, digs);
1304 SET_MP_ONE (g, digs);
1305 (void) mod_mp (p, f, x, g, digs);
1306 if (IS_ZERO_MP (f)) {
1307 SET_MP_ZERO (z, digs);
1308 A68_SP = pop_sp;
1309 return z;
1310 }
1311 SET_MP_MINUS_HALF (half, digs);
1312 (void) sub_mp (p, h, f, half, digs);
1313 if (MP_DIGIT (h, 1) < 0) {
1314 (void) plus_one_mp (p, f, f, digs);
1315 } else {
1316 SET_MP_HALF (half, digs);
1317 (void) sub_mp (p, h, f, half, digs);
1318 if (MP_DIGIT (h, 1) < 0) {
1319 (void) minus_one_mp (p, f, f, digs);
1320 }
1321 }
1322 BOOL_T neg = MP_DIGIT (f, 1) < 0;
1323 (void) abs_mp (p, f, f, digs);
1324 SET_MP_QUART (g, digs);
1325 (void) sub_mp (p, g, f, g, digs);
1326 if (IS_ZERO_MP (g)) {
1327 if (neg) {
1328 SET_MP_MINUS_ONE (z, digs);
1329 } else {
1330 SET_MP_ONE (z, digs);
1331 }
1332 A68_SP = pop_sp;
1333 return z;
1334 }
1335 (void) mp_pi (p, f, MP_PI, digs);
1336 (void) mul_mp (p, g, x, f, digs);
1337 if (tan_mp (p, z, g, digs) == NaN_MP) {
1338 errno = EDOM;
1339 A68_SP = pop_sp;
1340 return NaN_MP;
1341 } else {
1342 A68_SP = pop_sp;
1343 return z;
1344 }
1345 }
1346
1347 //! @brief PROC (LONG REAL) LONG REAL cotpi
1348
1349 MP_T *cotpi_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1350 {
1351 ADDR_T pop_sp = A68_SP;
1352 MP_T *f = nil_mp (p, digs);
1353 MP_T *g = lit_mp (p, 1, 0, digs);
1354 MP_T *h = nil_mp (p, digs);
1355 MP_T *half = nil_mp (p, digs);
1356 (void) mod_mp (p, f, x, g, digs);
1357 if (IS_ZERO_MP (f)) {
1358 errno = EDOM;
1359 A68_SP = pop_sp;
1360 return NaN_MP;
1361 }
1362 SET_MP_MINUS_HALF (half, digs);
1363 (void) sub_mp (p, h, f, half, digs);
1364 if (MP_DIGIT (h, 1) < 0) {
1365 (void) plus_one_mp (p, f, f, digs);
1366 } else {
1367 SET_MP_HALF (half, digs);
1368 (void) sub_mp (p, h, f, half, digs);
1369 if (MP_DIGIT (h, 1) < 0) {
1370 (void) minus_one_mp (p, f, f, digs);
1371 }
1372 }
1373 BOOL_T neg = MP_DIGIT (f, 1) < 0;
1374 (void) abs_mp (p, f, f, digs);
1375 SET_MP_QUART (g, digs);
1376 (void) sub_mp (p, g, f, g, digs);
1377 if (IS_ZERO_MP (g)) {
1378 if (neg) {
1379 SET_MP_MINUS_ONE (z, digs);
1380 } else {
1381 SET_MP_ONE (z, digs);
1382 }
1383 A68_SP = pop_sp;
1384 return z;
1385 }
1386 (void) mp_pi (p, f, MP_PI, digs);
1387 (void) mul_mp (p, g, x, f, digs);
1388 if (cot_mp (p, z, g, digs) == NaN_MP) {
1389 errno = EDOM;
1390 A68_SP = pop_sp;
1391 return NaN_MP;
1392 } else {
1393 A68_SP = pop_sp;
1394 return z;
1395 }
1396 }
