mp-complex.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 COMPLEX 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 LONG COMPLEX multiplication
33
34 MP_T *cmul_mp (NODE_T * p, MP_T * a, MP_T * b, MP_T * c, MP_T * d, int digs)
35 {
36 ADDR_T pop_sp = A68_SP;
37 int gdigs = FUN_DIGITS (digs);
38 MP_T *la = len_mp (p, a, digs, gdigs);
39 MP_T *lb = len_mp (p, b, digs, gdigs);
40 MP_T *lc = len_mp (p, c, digs, gdigs);
41 MP_T *ld = len_mp (p, d, digs, gdigs);
42 MP_T *ac = nil_mp (p, gdigs);
43 MP_T *bd = nil_mp (p, gdigs);
44 MP_T *ad = nil_mp (p, gdigs);
45 MP_T *bc = nil_mp (p, gdigs);
46 (void) mul_mp (p, ac, la, lc, gdigs);
47 (void) mul_mp (p, bd, lb, ld, gdigs);
48 (void) mul_mp (p, ad, la, ld, gdigs);
49 (void) mul_mp (p, bc, lb, lc, gdigs);
50 (void) sub_mp (p, la, ac, bd, gdigs);
51 (void) add_mp (p, lb, ad, bc, gdigs);
52 (void) shorten_mp (p, a, digs, la, gdigs);
53 (void) shorten_mp (p, b, digs, lb, gdigs);
54 A68_SP = pop_sp;
55 return a;
56 }
57
58 //! @brief LONG COMPLEX division
59
60 MP_T *cdiv_mp (NODE_T * p, MP_T * a, MP_T * b, MP_T * c, MP_T * d, int digs)
61 {
62 ADDR_T pop_sp = A68_SP;
63 if (MP_DIGIT (c, 1) == (MP_T) 0 && MP_DIGIT (d, 1) == (MP_T) 0) {
64 errno = ERANGE;
65 return NaN_MP;
66 }
67 MP_T *q = nil_mp (p, digs);
68 MP_T *r = nil_mp (p, digs);
69 (void) move_mp (q, c, digs);
70 (void) move_mp (r, d, digs);
71 MP_DIGIT (q, 1) = ABS (MP_DIGIT (q, 1));
72 MP_DIGIT (r, 1) = ABS (MP_DIGIT (r, 1));
73 (void) sub_mp (p, q, q, r, digs);
74 if (MP_DIGIT (q, 1) >= 0) {
75 if (div_mp (p, q, d, c, digs) == NaN_MP) {
76 errno = ERANGE;
77 return NaN_MP;
78 }
79 (void) mul_mp (p, r, d, q, digs);
80 (void) add_mp (p, r, r, c, digs);
81 (void) mul_mp (p, c, b, q, digs);
82 (void) add_mp (p, c, c, a, digs);
83 (void) div_mp (p, c, c, r, digs);
84 (void) mul_mp (p, d, a, q, digs);
85 (void) sub_mp (p, d, b, d, digs);
86 (void) div_mp (p, d, d, r, digs);
87 } else {
88 if (div_mp (p, q, c, d, digs) == NaN_MP) {
89 errno = ERANGE;
90 return NaN_MP;
91 }
92 (void) mul_mp (p, r, c, q, digs);
93 (void) add_mp (p, r, r, d, digs);
94 (void) mul_mp (p, c, a, q, digs);
95 (void) add_mp (p, c, c, b, digs);
96 (void) div_mp (p, c, c, r, digs);
97 (void) mul_mp (p, d, b, q, digs);
98 (void) sub_mp (p, d, d, a, digs);
99 (void) div_mp (p, d, d, r, digs);
100 }
101 (void) move_mp (a, c, digs);
102 (void) move_mp (b, d, digs);
103 A68_SP = pop_sp;
104 return a;
105 }
106
107 //! @brief PROC (LONG COMPLEX) LONG COMPLEX csqrt
108
109 MP_T *csqrt_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
110 {
111 ADDR_T pop_sp = A68_SP;
112 int gdigs = FUN_DIGITS (digs);
113 MP_T *re = len_mp (p, r, digs, gdigs);
114 MP_T *im = len_mp (p, i, digs, gdigs);
115 if (IS_ZERO_MP (re) && IS_ZERO_MP (im)) {
116 SET_MP_ZERO (re, gdigs);
117 SET_MP_ZERO (im, gdigs);
118 } else {
119 MP_T *c1 = lit_mp (p, 1, 0, gdigs);
120 MP_T *t = nil_mp (p, gdigs);
121 MP_T *x = nil_mp (p, gdigs);
122 MP_T *y = nil_mp (p, gdigs);
123 MP_T *u = nil_mp (p, gdigs);
124 MP_T *v = nil_mp (p, gdigs);
125 MP_T *w = nil_mp (p, gdigs);
126 SET_MP_ONE (c1, gdigs);
127 (void) move_mp (x, re, gdigs);
128 (void) move_mp (y, im, gdigs);
129 MP_DIGIT (x, 1) = ABS (MP_DIGIT (x, 1));
130 MP_DIGIT (y, 1) = ABS (MP_DIGIT (y, 1));
131 (void) sub_mp (p, w, x, y, gdigs);
132 if (MP_DIGIT (w, 1) >= 0) {
133 (void) div_mp (p, t, y, x, gdigs);
134 (void) mul_mp (p, v, t, t, gdigs);
135 (void) add_mp (p, u, c1, v, gdigs);
136 (void) sqrt_mp (p, v, u, gdigs);
137 (void) add_mp (p, u, c1, v, gdigs);
138 (void) half_mp (p, v, u, gdigs);
139 (void) sqrt_mp (p, u, v, gdigs);
140 (void) sqrt_mp (p, v, x, gdigs);
141 (void) mul_mp (p, w, u, v, gdigs);
142 } else {
143 (void) div_mp (p, t, x, y, gdigs);
144 (void) mul_mp (p, v, t, t, gdigs);
145 (void) add_mp (p, u, c1, v, gdigs);
146 (void) sqrt_mp (p, v, u, gdigs);
147 (void) add_mp (p, u, t, v, gdigs);
148 (void) half_mp (p, v, u, gdigs);
149 (void) sqrt_mp (p, u, v, gdigs);
150 (void) sqrt_mp (p, v, y, gdigs);
151 (void) mul_mp (p, w, u, v, gdigs);
152 }
153 if (MP_DIGIT (re, 1) >= 0) {
154 (void) move_mp (re, w, gdigs);
155 (void) add_mp (p, u, w, w, gdigs);
156 (void) div_mp (p, im, im, u, gdigs);
157 } else {
158 if (MP_DIGIT (im, 1) < 0) {
159 MP_DIGIT (w, 1) = -MP_DIGIT (w, 1);
160 }
161 (void) add_mp (p, v, w, w, gdigs);
162 (void) div_mp (p, re, im, v, gdigs);
163 (void) move_mp (im, w, gdigs);
164 }
165 }
166 (void) shorten_mp (p, r, digs, re, gdigs);
167 (void) shorten_mp (p, i, digs, im, gdigs);
168 A68_SP = pop_sp;
169 return r;
170 }
171
172 //! @brief PROC (LONG COMPLEX) LONG COMPLEX cexp
173
174 MP_T *cexp_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
175 {
176 ADDR_T pop_sp = A68_SP;
177 int gdigs = FUN_DIGITS (digs);
178 MP_T *re = len_mp (p, r, digs, gdigs);
179 MP_T *im = len_mp (p, i, digs, gdigs);
180 MP_T *u = nil_mp (p, gdigs);
181 (void) exp_mp (p, u, re, gdigs);
182 (void) cos_mp (p, re, im, gdigs);
183 (void) sin_mp (p, im, im, gdigs);
184 (void) mul_mp (p, re, re, u, gdigs);
185 (void) mul_mp (p, im, im, u, gdigs);
186 (void) shorten_mp (p, r, digs, re, gdigs);
187 (void) shorten_mp (p, i, digs, im, gdigs);
188 A68_SP = pop_sp;
189 return r;
190 }
191
192 //! @brief PROC (LONG COMPLEX) LONG COMPLEX cln
193
194 MP_T *cln_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
195 {
196 ADDR_T pop_sp = A68_SP;
197 int gdigs = FUN_DIGITS (digs);
198 MP_T *re = len_mp (p, r, digs, gdigs);
199 MP_T *im = len_mp (p, i, digs, gdigs);
200 MP_T *s = nil_mp (p, gdigs);
201 MP_T *t = nil_mp (p, gdigs);
202 MP_T *u = nil_mp (p, gdigs);
203 MP_T *v = nil_mp (p, gdigs);
204 (void) move_mp (u, re, gdigs);
205 (void) move_mp (v, im, gdigs);
206 (void) hypot_mp (p, s, u, v, gdigs);
207 (void) move_mp (u, re, gdigs);
208 (void) move_mp (v, im, gdigs);
209 (void) atan2_mp (p, t, u, v, gdigs);
210 (void) ln_mp (p, re, s, gdigs);
211 (void) move_mp (im, t, gdigs);
212 (void) shorten_mp (p, r, digs, re, gdigs);
213 (void) shorten_mp (p, i, digs, im, gdigs);
214 A68_SP = pop_sp;
215 return r;
216 }
217
218 //! @brief PROC (LONG COMPLEX) LONG COMPLEX csin
219
220 MP_T *csin_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
221 {
222 ADDR_T pop_sp = A68_SP;
223 int gdigs = FUN_DIGITS (digs);
224 MP_T *re = len_mp (p, r, digs, gdigs);
225 MP_T *im = len_mp (p, i, digs, gdigs);
226 MP_T *s = nil_mp (p, gdigs);
227 MP_T *c = nil_mp (p, gdigs);
228 MP_T *sh = nil_mp (p, gdigs);
229 MP_T *ch = nil_mp (p, gdigs);
230 if (IS_ZERO_MP (im)) {
231 (void) sin_mp (p, re, re, gdigs);
232 SET_MP_ZERO (im, gdigs);
233 } else {
234 (void) sin_mp (p, s, re, gdigs);
235 (void) cos_mp (p, c, re, gdigs);
236 (void) hyp_mp (p, sh, ch, im, gdigs);
237 (void) mul_mp (p, re, s, ch, gdigs);
238 (void) mul_mp (p, im, c, sh, gdigs);
239 }
240 (void) shorten_mp (p, r, digs, re, gdigs);
241 (void) shorten_mp (p, i, digs, im, gdigs);
242 A68_SP = pop_sp;
243 return r;
244 }
245
246 //! @brief PROC (LONG COMPLEX) LONG COMPLEX ccos
247
248 MP_T *ccos_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
249 {
250 ADDR_T pop_sp = A68_SP;
251 int gdigs = FUN_DIGITS (digs);
252 MP_T *re = len_mp (p, r, digs, gdigs);
253 MP_T *im = len_mp (p, i, digs, gdigs);
254 MP_T *s = nil_mp (p, gdigs);
255 MP_T *c = nil_mp (p, gdigs);
256 MP_T *sh = nil_mp (p, gdigs);
257 MP_T *ch = nil_mp (p, gdigs);
258 if (IS_ZERO_MP (im)) {
259 (void) cos_mp (p, re, re, gdigs);
260 SET_MP_ZERO (im, gdigs);
261 } else {
262 (void) sin_mp (p, s, re, gdigs);
263 (void) cos_mp (p, c, re, gdigs);
264 (void) hyp_mp (p, sh, ch, im, gdigs);
265 MP_DIGIT (sh, 1) = -MP_DIGIT (sh, 1);
266 (void) mul_mp (p, re, c, ch, gdigs);
267 (void) mul_mp (p, im, s, sh, gdigs);
268 }
269 (void) shorten_mp (p, r, digs, re, gdigs);
270 (void) shorten_mp (p, i, digs, im, gdigs);
271 A68_SP = pop_sp;
272 return r;
273 }
274
275 //! @brief PROC (LONG COMPLEX) LONG COMPLEX ctan
276
277 MP_T *ctan_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
278 {
279 ADDR_T pop_sp = A68_SP;
280 errno = 0;
281 MP_T *s = nil_mp (p, digs);
282 MP_T *t = nil_mp (p, digs);
283 MP_T *u = nil_mp (p, digs);
284 MP_T *v = nil_mp (p, digs);
285 (void) move_mp (u, r, digs);
286 (void) move_mp (v, i, digs);
287 (void) csin_mp (p, u, v, digs);
288 (void) move_mp (s, u, digs);
289 (void) move_mp (t, v, digs);
290 (void) move_mp (u, r, digs);
291 (void) move_mp (v, i, digs);
292 (void) ccos_mp (p, u, v, digs);
293 (void) cdiv_mp (p, s, t, u, v, digs);
294 (void) move_mp (r, s, digs);
295 (void) move_mp (i, t, digs);
296 A68_SP = pop_sp;
297 return r;
298 }
299
300 //! @brief PROC (LONG COMPLEX) LONG COMPLEX casin
301
302 MP_T *casin_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
303 {
304 ADDR_T pop_sp = A68_SP;
305 int gdigs = FUN_DIGITS (digs);
306 MP_T *re = len_mp (p, r, digs, gdigs);
307 MP_T *im = len_mp (p, i, digs, gdigs);
308 BOOL_T negim = MP_DIGIT (im, 1) < 0;
309 if (IS_ZERO_MP (im)) {
310 BOOL_T neg = MP_DIGIT (re, 1) < 0;
311 if (acos_mp (p, im, re, gdigs) == NaN_MP) {
312 errno = 0; // Ignore the acos error
313 MP_DIGIT (re, 1) = ABS (MP_DIGIT (re, 1));
314 (void) acosh_mp (p, im, re, gdigs);
315 }
316 (void) mp_pi (p, re, MP_HALF_PI, gdigs);
317 if (neg) {
318 MP_DIGIT (re, 1) = -MP_DIGIT (re, 1);
319 }
320 } else {
321 MP_T *c1 = lit_mp (p, 1, 0, gdigs);
322 MP_T *u = nil_mp (p, gdigs);
323 MP_T *v = nil_mp (p, gdigs);
324 MP_T *a = nil_mp (p, gdigs);
325 MP_T *b = nil_mp (p, gdigs);
326 // u=sqrt((r+1)^2+i^2), v=sqrt((r-1)^2+i^2).
327 (void) add_mp (p, a, re, c1, gdigs);
328 (void) sub_mp (p, b, re, c1, gdigs);
329 (void) hypot_mp (p, u, a, im, gdigs);
330 (void) hypot_mp (p, v, b, im, gdigs);
331 // a=(u+v)/2, b=(u-v)/2.
332 (void) add_mp (p, a, u, v, gdigs);
333 (void) half_mp (p, a, a, gdigs);
334 (void) sub_mp (p, b, u, v, gdigs);
335 (void) half_mp (p, b, b, gdigs);
336 // r=asin(b), i=ln(a+sqrt(a^2-1)).
337 (void) mul_mp (p, u, a, a, gdigs);
338 (void) sub_mp (p, u, u, c1, gdigs);
339 (void) sqrt_mp (p, u, u, gdigs);
340 (void) add_mp (p, u, a, u, gdigs);
341 (void) ln_mp (p, im, u, gdigs);
342 (void) asin_mp (p, re, b, gdigs);
343 }
344 if (negim) {
345 MP_DIGIT (im, 1) = -MP_DIGIT (im, 1);
346 }
347 (void) shorten_mp (p, r, digs, re, gdigs);
348 (void) shorten_mp (p, i, digs, im, gdigs);
349 A68_SP = pop_sp;
350 return re;
351 }
352
353 //! @brief PROC (LONG COMPLEX) LONG COMPLEX cacos
354
355 MP_T *cacos_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
356 {
357 ADDR_T pop_sp = A68_SP;
358 int gdigs = FUN_DIGITS (digs);
359 MP_T *re = len_mp (p, r, digs, gdigs);
360 MP_T *im = len_mp (p, i, digs, gdigs);
361 BOOL_T negim = MP_DIGIT (im, 1) < 0;
362 if (IS_ZERO_MP (im)) {
363 BOOL_T neg = MP_DIGIT (re, 1) < 0;
364 if (acos_mp (p, im, re, gdigs) == NaN_MP) {
365 errno = 0; // Ignore the acos error
366 MP_DIGIT (re, 1) = ABS (MP_DIGIT (re, 1));
367 (void) acosh_mp (p, im, re, gdigs);
368 MP_DIGIT (im, 1) = -MP_DIGIT (im, 1);
369 }
370 if (neg) {
371 (void) mp_pi (p, re, MP_PI, gdigs);
372 } else {
373 SET_MP_ZERO (re, gdigs);
374 }
375 } else {
376 MP_T *c1 = lit_mp (p, 1, 0, gdigs);
377 MP_T *u = nil_mp (p, gdigs);
378 MP_T *v = nil_mp (p, gdigs);
379 MP_T *a = nil_mp (p, gdigs);
380 MP_T *b = nil_mp (p, gdigs);
381 // u=sqrt((r+1)^2+i^2), v=sqrt((r-1)^2+i^2).
382 (void) add_mp (p, a, re, c1, gdigs);
383 (void) sub_mp (p, b, re, c1, gdigs);
384 (void) hypot_mp (p, u, a, im, gdigs);
385 (void) hypot_mp (p, v, b, im, gdigs);
386 // a=(u+v)/2, b=(u-v)/2.
387 (void) add_mp (p, a, u, v, gdigs);
388 (void) half_mp (p, a, a, gdigs);
389 (void) sub_mp (p, b, u, v, gdigs);
390 (void) half_mp (p, b, b, gdigs);
391 // r=acos(b), i=-ln(a+sqrt(a^2-1)).
392 (void) mul_mp (p, u, a, a, gdigs);
393 (void) sub_mp (p, u, u, c1, gdigs);
394 (void) sqrt_mp (p, u, u, gdigs);
395 (void) add_mp (p, u, a, u, gdigs);
396 (void) ln_mp (p, im, u, gdigs);
397 (void) acos_mp (p, re, b, gdigs);
398 }
399 if (!negim) {
400 MP_DIGIT (im, 1) = -MP_DIGIT (im, 1);
401 }
402 (void) shorten_mp (p, r, digs, re, gdigs);
403 (void) shorten_mp (p, i, digs, im, gdigs);
404 A68_SP = pop_sp;
405 return re;
406 }
407
408 //! @brief PROC (LONG COMPLEX) LONG COMPLEX catan
409
410 MP_T *catan_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
411 {
412 ADDR_T pop_sp = A68_SP;
413 int gdigs = FUN_DIGITS (digs);
414 MP_T *re = len_mp (p, r, digs, gdigs);
415 MP_T *im = len_mp (p, i, digs, gdigs);
416 MP_T *u = nil_mp (p, gdigs);
417 MP_T *v = nil_mp (p, gdigs);
418 if (IS_ZERO_MP (im)) {
419 (void) atan_mp (p, u, re, gdigs);
420 SET_MP_ZERO (v, gdigs);
421 } else {
422 MP_T *c1 = lit_mp (p, 1, 0, gdigs);
423 MP_T *a = nil_mp (p, gdigs);
424 MP_T *b = nil_mp (p, gdigs);
425 // a=sqrt(r^2+(i+1)^2), b=sqrt(r^2+(i-1)^2).
426 (void) add_mp (p, a, im, c1, gdigs);
427 (void) sub_mp (p, b, im, c1, gdigs);
428 (void) hypot_mp (p, u, re, a, gdigs);
429 (void) hypot_mp (p, v, re, b, gdigs);
430 // im=ln(a/b).
431 (void) div_mp (p, u, u, v, gdigs);
432 (void) ln_mp (p, v, u, gdigs);
433 // re=atan(2r/(1-r^2-i^2)).
434 (void) mul_mp (p, a, re, re, gdigs);
435 (void) mul_mp (p, b, im, im, gdigs);
436 (void) add_mp (p, a, a, b, gdigs);
437 (void) sub_mp (p, u, c1, a, gdigs);
438 if (IS_ZERO_MP (u)) {
439 (void) mp_pi (p, u, MP_HALF_PI, gdigs);
440 } else {
441 int neg = MP_DIGIT (u, 1) < 0;
442 (void) add_mp (p, a, re, re, gdigs);
443 (void) div_mp (p, a, a, u, gdigs);
444 (void) atan_mp (p, u, a, gdigs);
445 if (neg) {
446 (void) mp_pi (p, a, MP_PI, gdigs);
447 if (MP_DIGIT (re, 1) < 0) {
448 (void) sub_mp (p, u, u, a, gdigs);
449 } else {
450 (void) add_mp (p, u, u, a, gdigs);
451 }
452 }
453 }
454 }
455 (void) shorten_mp (p, r, digs, u, gdigs);
456 (void) shorten_mp (p, i, digs, v, gdigs);
457 (void) half_mp (p, r, r, digs);
458 (void) half_mp (p, i, i, digs);
459 A68_SP = pop_sp;
460 return re;
461 }
462
463 //! @brief PROC (LONG COMPLEX) LONG COMPLEX csinh
464
465 MP_T *csinh_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
466 {
467 ADDR_T pop_sp = A68_SP;
468 int gdigs = FUN_DIGITS (digs);
469 MP_T *u = len_mp (p, r, digs, gdigs);
470 MP_T *v = len_mp (p, i, digs, gdigs);
471 MP_T *s = nil_mp (p, gdigs);
472 MP_T *t = nil_mp (p, gdigs);
473 // sinh (z) = -i sin (iz)
474 SET_MP_ONE (t, gdigs);
475 (void) cmul_mp (p, u, v, s, t, gdigs);
476 (void) csin_mp (p, u, v, gdigs);
477 SET_MP_ZERO (s, gdigs);
478 SET_MP_MINUS_ONE (t, gdigs);
479 (void) cmul_mp (p, u, v, s, t, gdigs);
480 //
481 (void) shorten_mp (p, r, digs, u, gdigs);
482 (void) shorten_mp (p, i, digs, v, gdigs);
483 A68_SP = pop_sp;
484 return r;
485 }
486
487 //! @brief PROC (LONG COMPLEX) LONG COMPLEX ccosh
488
489 MP_T *ccosh_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
490 {
491 ADDR_T pop_sp = A68_SP;
492 int gdigs = FUN_DIGITS (digs);
493 MP_T *u = len_mp (p, r, digs, gdigs);
494 MP_T *v = len_mp (p, i, digs, gdigs);
495 MP_T *s = nil_mp (p, gdigs);
496 MP_T *t = nil_mp (p, gdigs);
497 // cosh (z) = cos (iz)
498 SET_MP_ZERO (s, digs);
499 SET_MP_ONE (t, gdigs);
500 (void) cmul_mp (p, u, v, s, t, gdigs);
501 (void) ccos_mp (p, u, v, gdigs);
502 //
503 (void) shorten_mp (p, r, digs, u, gdigs);
504 (void) shorten_mp (p, i, digs, v, gdigs);
505 A68_SP = pop_sp;
506 return r;
507 }
508
509 //! @brief PROC (LONG COMPLEX) LONG COMPLEX ctanh
510
511 MP_T *ctanh_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
512 {
513 ADDR_T pop_sp = A68_SP;
514 int gdigs = FUN_DIGITS (digs);
515 MP_T *u = len_mp (p, r, digs, gdigs);
516 MP_T *v = len_mp (p, i, digs, gdigs);
517 MP_T *s = nil_mp (p, gdigs);
518 MP_T *t = nil_mp (p, gdigs);
519 // tanh (z) = -i tan (iz)
520 SET_MP_ONE (t, gdigs);
521 (void) cmul_mp (p, u, v, s, t, gdigs);
522 (void) ctan_mp (p, u, v, gdigs);
523 SET_MP_ZERO (u, gdigs);
524 SET_MP_MINUS_ONE (v, gdigs);
525 (void) cmul_mp (p, u, v, s, t, gdigs);
526 //
527 (void) shorten_mp (p, r, digs, u, gdigs);
528 (void) shorten_mp (p, i, digs, v, gdigs);
529 A68_SP = pop_sp;
530 return r;
531 }
532
533 //! @brief PROC (LONG COMPLEX) LONG COMPLEX casinh
534
535 MP_T *casinh_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
536 {
537 ADDR_T pop_sp = A68_SP;
538 int gdigs = FUN_DIGITS (digs);
539 MP_T *u = len_mp (p, r, digs, gdigs);
540 MP_T *v = len_mp (p, i, digs, gdigs);
541 MP_T *s = nil_mp (p, gdigs);
542 MP_T *t = nil_mp (p, gdigs);
543 // asinh (z) = i asin (-iz)
544 SET_MP_MINUS_ONE (t, gdigs);
545 (void) cmul_mp (p, u, v, s, t, gdigs);
546 (void) casin_mp (p, u, v, gdigs);
547 SET_MP_ZERO (s, gdigs);
548 SET_MP_ONE (t, gdigs);
549 (void) cmul_mp (p, u, v, s, t, gdigs);
550 //
551 (void) shorten_mp (p, r, digs, u, gdigs);
552 (void) shorten_mp (p, i, digs, v, gdigs);
553 A68_SP = pop_sp;
554 return r;
555 }
556
557 //! @brief PROC (LONG COMPLEX) LONG COMPLEX cacosh
558
559 MP_T *cacosh_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
560 {
561 ADDR_T pop_sp = A68_SP;
562 int gdigs = FUN_DIGITS (digs);
563 MP_T *u = len_mp (p, r, digs, gdigs);
564 MP_T *v = len_mp (p, i, digs, gdigs);
565 MP_T *s = nil_mp (p, gdigs);
566 MP_T *t = nil_mp (p, gdigs);
567 // acosh (z) = i * acos (z)
568 (void) cacos_mp (p, u, v, gdigs);
569 SET_MP_ZERO (s, gdigs);
570 SET_MP_ONE (t, gdigs);
571 (void) cmul_mp (p, u, v, s, t, gdigs);
572 //
573 (void) shorten_mp (p, r, digs, u, gdigs);
574 (void) shorten_mp (p, i, digs, v, gdigs);
575 A68_SP = pop_sp;
576 return r;
577 }
578
579 //! @brief PROC (LONG COMPLEX) LONG COMPLEX catanh
580
581 MP_T *catanh_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
582 {
583 ADDR_T pop_sp = A68_SP;
584 int gdigs = FUN_DIGITS (digs);
585 MP_T *u = len_mp (p, r, digs, gdigs);
586 MP_T *v = len_mp (p, i, digs, gdigs);
587 MP_T *s = nil_mp (p, gdigs);
588 MP_T *t = nil_mp (p, gdigs);
589 // atanh (z) = i * atan (-iz)
590 SET_MP_MINUS_ONE (t, gdigs);
591 (void) cmul_mp (p, u, v, s, t, gdigs);
592 (void) catan_mp (p, u, v, gdigs);
593 SET_MP_ZERO (s, gdigs);
594 SET_MP_ONE (t, gdigs);
595 (void) cmul_mp (p, u, v, s, t, gdigs);
596 //
597 (void) shorten_mp (p, r, digs, u, gdigs);
598 (void) shorten_mp (p, i, digs, v, gdigs);
599 A68_SP = pop_sp;
600 return r;
601 }
