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