1 /*
2 * commath - extended precision complex arithmetic primitive routines
3 *
4 * Copyright (C) 1999-2007,2021 David I. Bell
5 *
6 * Calc is open software; you can redistribute it and/or modify it under
7 * the terms of the version 2.1 of the GNU Lesser General Public License
8 * as published by the Free Software Foundation.
9 *
10 * Calc is distributed in the hope that it will be useful, but WITHOUT
11 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
12 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General
13 * Public License for more details.
14 *
15 * A copy of version 2.1 of the GNU Lesser General Public License is
16 * distributed with calc under the filename COPYING-LGPL. You should have
17 * received a copy with calc; if not, write to Free Software Foundation, Inc.
18 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
19 *
20 * Under source code control: 1990/02/15 01:48:10
21 * File existed as early as: before 1990
22 *
23 * Share and enjoy! :-) http://www.isthe.com/chongo/tech/comp/calc/
24 */
25
26
27 #include "cmath.h"
28
29
30 #include "banned.h" /* include after system header <> includes */
31
32
33 COMPLEX _czero_ = { &_qzero_, &_qzero_, 1 };
34 COMPLEX _cone_ = { &_qone_, &_qzero_, 1 };
35 COMPLEX _conei_ = { &_qzero_, &_qone_, 1 };
36
37 STATIC COMPLEX _cnegone_ = { &_qnegone_, &_qzero_, 1 };
38
39
40 /*
41 * Add two complex numbers.
42 */
43 COMPLEX *
c_add(COMPLEX * c1,COMPLEX * c2)44 c_add(COMPLEX *c1, COMPLEX *c2)
45 {
46 COMPLEX *r;
47
48 if (ciszero(c1))
49 return clink(c2);
50 if (ciszero(c2))
51 return clink(c1);
52 r = comalloc();
53 if (!qiszero(c1->real) || !qiszero(c2->real)) {
54 qfree(r->real);
55 r->real = qqadd(c1->real, c2->real);
56 }
57 if (!qiszero(c1->imag) || !qiszero(c2->imag)) {
58 qfree(r->imag);
59 r->imag = qqadd(c1->imag, c2->imag);
60 }
61 return r;
62 }
63
64
65 /*
66 * Subtract two complex numbers.
67 */
68 COMPLEX *
c_sub(COMPLEX * c1,COMPLEX * c2)69 c_sub(COMPLEX *c1, COMPLEX *c2)
70 {
71 COMPLEX *r;
72
73 if ((c1->real == c2->real) && (c1->imag == c2->imag))
74 return clink(&_czero_);
75 if (ciszero(c2))
76 return clink(c1);
77 r = comalloc();
78 if (!qiszero(c1->real) || !qiszero(c2->real)) {
79 qfree(r->real);
80 r->real = qsub(c1->real, c2->real);
81 }
82 if (!qiszero(c1->imag) || !qiszero(c2->imag)) {
83 qfree(r->imag);
84 r->imag = qsub(c1->imag, c2->imag);
85 }
86 return r;
87 }
88
89
90 /*
91 * Multiply two complex numbers.
92 * This saves one multiplication over the obvious algorithm by
93 * trading it for several extra additions, as follows. Let
94 * q1 = (a + b) * (c + d)
95 * q2 = a * c
96 * q3 = b * d
97 * Then (a+bi) * (c+di) = (q2 - q3) + (q1 - q2 - q3)i.
98 */
99 COMPLEX *
c_mul(COMPLEX * c1,COMPLEX * c2)100 c_mul(COMPLEX *c1, COMPLEX *c2)
101 {
102 COMPLEX *r;
103 NUMBER *q1, *q2, *q3, *q4;
104
105 if (ciszero(c1) || ciszero(c2))
106 return clink(&_czero_);
107 if (cisone(c1))
108 return clink(c2);
109 if (cisone(c2))
110 return clink(c1);
111 if (cisreal(c2))
112 return c_mulq(c1, c2->real);
113 if (cisreal(c1))
114 return c_mulq(c2, c1->real);
115 /*
116 * Need to do the full calculation.
117 */
118 r = comalloc();
119 q2 = qqadd(c1->real, c1->imag);
120 q3 = qqadd(c2->real, c2->imag);
121 q1 = qmul(q2, q3);
122 qfree(q2);
123 qfree(q3);
124 q2 = qmul(c1->real, c2->real);
125 q3 = qmul(c1->imag, c2->imag);
126 q4 = qqadd(q2, q3);
127 qfree(r->real);
128 r->real = qsub(q2, q3);
129 qfree(r->imag);
130 r->imag = qsub(q1, q4);
131 qfree(q1);
132 qfree(q2);
133 qfree(q3);
134 qfree(q4);
135 return r;
136 }
137
138
139 /*
140 * Square a complex number.
141 */
142 COMPLEX *
c_square(COMPLEX * c)143 c_square(COMPLEX *c)
144 {
145 COMPLEX *r;
146 NUMBER *q1, *q2;
147
148 if (ciszero(c))
149 return clink(&_czero_);
150 if (cisrunit(c))
151 return clink(&_cone_);
152 if (cisiunit(c))
153 return clink(&_cnegone_);
154 r = comalloc();
155 if (cisreal(c)) {
156 qfree(r->real);
157 r->real = qsquare(c->real);
158 return r;
159 }
160 if (cisimag(c)) {
161 qfree(r->real);
162 q1 = qsquare(c->imag);
163 r->real = qneg(q1);
164 qfree(q1);
165 return r;
166 }
167 q1 = qsquare(c->real);
168 q2 = qsquare(c->imag);
169 qfree(r->real);
170 r->real = qsub(q1, q2);
171 qfree(q1);
172 qfree(q2);
173 qfree(r->imag);
174 q1 = qmul(c->real, c->imag);
175 r->imag = qscale(q1, 1L);
176 qfree(q1);
177 return r;
178 }
179
180
181 /*
182 * Divide two complex numbers.
183 */
184 COMPLEX *
c_div(COMPLEX * c1,COMPLEX * c2)185 c_div(COMPLEX *c1, COMPLEX *c2)
186 {
187 COMPLEX *r;
188 NUMBER *q1, *q2, *q3, *den;
189
190 if (ciszero(c2)) {
191 math_error("Division by zero");
192 /*NOTREACHED*/
193 }
194 if ((c1->real == c2->real) && (c1->imag == c2->imag))
195 return clink(&_cone_);
196 r = comalloc();
197 if (cisreal(c1) && cisreal(c2)) {
198 qfree(r->real);
199 r->real = qqdiv(c1->real, c2->real);
200 return r;
201 }
202 if (cisimag(c1) && cisimag(c2)) {
203 qfree(r->real);
204 r->real = qqdiv(c1->imag, c2->imag);
205 return r;
206 }
207 if (cisimag(c1) && cisreal(c2)) {
208 qfree(r->imag);
209 r->imag = qqdiv(c1->imag, c2->real);
210 return r;
211 }
212 if (cisreal(c1) && cisimag(c2)) {
213 qfree(r->imag);
214 q1 = qqdiv(c1->real, c2->imag);
215 r->imag = qneg(q1);
216 qfree(q1);
217 return r;
218 }
219 if (cisreal(c2)) {
220 qfree(r->real);
221 qfree(r->imag);
222 r->real = qqdiv(c1->real, c2->real);
223 r->imag = qqdiv(c1->imag, c2->real);
224 return r;
225 }
226 q1 = qsquare(c2->real);
227 q2 = qsquare(c2->imag);
228 den = qqadd(q1, q2);
229 qfree(q1);
230 qfree(q2);
231 q1 = qmul(c1->real, c2->real);
232 q2 = qmul(c1->imag, c2->imag);
233 q3 = qqadd(q1, q2);
234 qfree(q1);
235 qfree(q2);
236 qfree(r->real);
237 r->real = qqdiv(q3, den);
238 qfree(q3);
239 q1 = qmul(c1->real, c2->imag);
240 q2 = qmul(c1->imag, c2->real);
241 q3 = qsub(q2, q1);
242 qfree(q1);
243 qfree(q2);
244 qfree(r->imag);
245 r->imag = qqdiv(q3, den);
246 qfree(q3);
247 qfree(den);
248 return r;
249 }
250
251
252 /*
253 * Invert a complex number.
254 */
255 COMPLEX *
c_inv(COMPLEX * c)256 c_inv(COMPLEX *c)
257 {
258 COMPLEX *r;
259 NUMBER *q1, *q2, *den;
260
261 if (ciszero(c)) {
262 math_error("Inverting zero");
263 /*NOTREACHED*/
264 }
265 r = comalloc();
266 if (cisreal(c)) {
267 qfree(r->real);
268 r->real = qinv(c->real);
269 return r;
270 }
271 if (cisimag(c)) {
272 q1 = qinv(c->imag);
273 qfree(r->imag);
274 r->imag = qneg(q1);
275 qfree(q1);
276 return r;
277 }
278 q1 = qsquare(c->real);
279 q2 = qsquare(c->imag);
280 den = qqadd(q1, q2);
281 qfree(q1);
282 qfree(q2);
283 qfree(r->real);
284 r->real = qqdiv(c->real, den);
285 q1 = qqdiv(c->imag, den);
286 qfree(r->imag);
287 r->imag = qneg(q1);
288 qfree(q1);
289 qfree(den);
290 return r;
291 }
292
293
294 /*
295 * Negate a complex number.
296 */
297 COMPLEX *
c_neg(COMPLEX * c)298 c_neg(COMPLEX *c)
299 {
300 COMPLEX *r;
301
302 if (ciszero(c))
303 return clink(&_czero_);
304 r = comalloc();
305 if (!qiszero(c->real)) {
306 qfree(r->real);
307 r->real = qneg(c->real);
308 }
309 if (!qiszero(c->imag)) {
310 qfree(r->imag);
311 r->imag = qneg(c->imag);
312 }
313 return r;
314 }
315
316
317 /*
318 * Take the integer part of a complex number.
319 * This means take the integer part of both components.
320 */
321 COMPLEX *
c_int(COMPLEX * c)322 c_int(COMPLEX *c)
323 {
324 COMPLEX *r;
325
326 if (cisint(c))
327 return clink(c);
328 r = comalloc();
329 qfree(r->real);
330 r->real = qint(c->real);
331 qfree(r->imag);
332 r->imag = qint(c->imag);
333 return r;
334 }
335
336
337 /*
338 * Take the fractional part of a complex number.
339 * This means take the fractional part of both components.
340 */
341 COMPLEX *
c_frac(COMPLEX * c)342 c_frac(COMPLEX *c)
343 {
344 COMPLEX *r;
345
346 if (cisint(c))
347 return clink(&_czero_);
348 r = comalloc();
349 qfree(r->real);
350 r->real = qfrac(c->real);
351 qfree(r->imag);
352 r->imag = qfrac(c->imag);
353 return r;
354 }
355
356
357 /*
358 * Take the conjugate of a complex number.
359 * This negates the complex part.
360 */
361 COMPLEX *
c_conj(COMPLEX * c)362 c_conj(COMPLEX *c)
363 {
364 COMPLEX *r;
365
366 if (cisreal(c))
367 return clink(c);
368 r = comalloc();
369 if (!qiszero(c->real)) {
370 qfree(r->real);
371 r->real = qlink(c->real);
372 }
373 qfree(r->imag);
374 r->imag = qneg(c->imag);
375 return r;
376 }
377
378
379 /*
380 * Return the real part of a complex number.
381 */
382 COMPLEX *
c_real(COMPLEX * c)383 c_real(COMPLEX *c)
384 {
385 COMPLEX *r;
386
387 if (cisreal(c))
388 return clink(c);
389 r = comalloc();
390 if (!qiszero(c->real)) {
391 qfree(r->real);
392 r->real = qlink(c->real);
393 }
394 return r;
395 }
396
397
398 /*
399 * Return the imaginary part of a complex number as a real.
400 */
401 COMPLEX *
c_imag(COMPLEX * c)402 c_imag(COMPLEX *c)
403 {
404 COMPLEX *r;
405
406 if (cisreal(c))
407 return clink(&_czero_);
408 r = comalloc();
409 qfree(r->real);
410 r->real = qlink(c->imag);
411 return r;
412 }
413
414
415 /*
416 * Add a real number to a complex number.
417 */
418 COMPLEX *
c_addq(COMPLEX * c,NUMBER * q)419 c_addq(COMPLEX *c, NUMBER *q)
420 {
421 COMPLEX *r;
422
423 if (qiszero(q))
424 return clink(c);
425 r = comalloc();
426 qfree(r->real);
427 qfree(r->imag);
428 r->real = qqadd(c->real, q);
429 r->imag = qlink(c->imag);
430 return r;
431 }
432
433
434 /*
435 * Subtract a real number from a complex number.
436 */
437 COMPLEX *
c_subq(COMPLEX * c,NUMBER * q)438 c_subq(COMPLEX *c, NUMBER *q)
439 {
440 COMPLEX *r;
441
442 if (qiszero(q))
443 return clink(c);
444 r = comalloc();
445 qfree(r->real);
446 qfree(r->imag);
447 r->real = qsub(c->real, q);
448 r->imag = qlink(c->imag);
449 return r;
450 }
451
452
453 /*
454 * Shift the components of a complex number left by the specified
455 * number of bits. Negative values shift to the right.
456 */
457 COMPLEX *
c_shift(COMPLEX * c,long n)458 c_shift(COMPLEX *c, long n)
459 {
460 COMPLEX *r;
461
462 if (ciszero(c) || (n == 0))
463 return clink(c);
464 r = comalloc();
465 qfree(r->real);
466 qfree(r->imag);
467 r->real = qshift(c->real, n);
468 r->imag = qshift(c->imag, n);
469 return r;
470 }
471
472
473 /*
474 * Scale a complex number by a power of two.
475 */
476 COMPLEX *
c_scale(COMPLEX * c,long n)477 c_scale(COMPLEX *c, long n)
478 {
479 COMPLEX *r;
480
481 if (ciszero(c) || (n == 0))
482 return clink(c);
483 r = comalloc();
484 qfree(r->real);
485 qfree(r->imag);
486 r->real = qscale(c->real, n);
487 r->imag = qscale(c->imag, n);
488 return r;
489 }
490
491
492 /*
493 * Multiply a complex number by a real number.
494 */
495 COMPLEX *
c_mulq(COMPLEX * c,NUMBER * q)496 c_mulq(COMPLEX *c, NUMBER *q)
497 {
498 COMPLEX *r;
499
500 if (qiszero(q))
501 return clink(&_czero_);
502 if (qisone(q))
503 return clink(c);
504 if (qisnegone(q))
505 return c_neg(c);
506 r = comalloc();
507 qfree(r->real);
508 qfree(r->imag);
509 r->real = qmul(c->real, q);
510 r->imag = qmul(c->imag, q);
511 return r;
512 }
513
514
515 /*
516 * Divide a complex number by a real number.
517 */
518 COMPLEX *
c_divq(COMPLEX * c,NUMBER * q)519 c_divq(COMPLEX *c, NUMBER *q)
520 {
521 COMPLEX *r;
522
523 if (qiszero(q)) {
524 math_error("Division by zero");
525 /*NOTREACHED*/
526 }
527 if (qisone(q))
528 return clink(c);
529 if (qisnegone(q))
530 return c_neg(c);
531 r = comalloc();
532 qfree(r->real);
533 qfree(r->imag);
534 r->real = qqdiv(c->real, q);
535 r->imag = qqdiv(c->imag, q);
536 return r;
537 }
538
539
540
541
542 /*
543 * Construct a complex number given the real and imaginary components.
544 */
545 COMPLEX *
qqtoc(NUMBER * q1,NUMBER * q2)546 qqtoc(NUMBER *q1, NUMBER *q2)
547 {
548 COMPLEX *r;
549
550 if (qiszero(q1) && qiszero(q2))
551 return clink(&_czero_);
552 r = comalloc();
553 qfree(r->real);
554 qfree(r->imag);
555 r->real = qlink(q1);
556 r->imag = qlink(q2);
557 return r;
558 }
559
560
561 /*
562 * Compare two complex numbers for equality, returning FALSE if they are equal,
563 * and TRUE if they differ.
564 */
565 BOOL
c_cmp(COMPLEX * c1,COMPLEX * c2)566 c_cmp(COMPLEX *c1, COMPLEX *c2)
567 {
568 BOOL i;
569
570 i = qcmp(c1->real, c2->real);
571 if (!i)
572 i = qcmp(c1->imag, c2->imag);
573 return i;
574 }
575
576
577 /*
578 * Compare two complex numbers and return a complex number with real and
579 * imaginary parts -1, 0 or 1 indicating relative values of the real and
580 * imaginary parts of the two numbers.
581 */
582 COMPLEX *
c_rel(COMPLEX * c1,COMPLEX * c2)583 c_rel(COMPLEX *c1, COMPLEX *c2)
584 {
585 COMPLEX *c;
586
587 c = comalloc();
588 qfree(c->real);
589 qfree(c->imag);
590 c->real = itoq((long) qrel(c1->real, c2->real));
591 c->imag = itoq((long) qrel(c1->imag, c2->imag));
592
593 return c;
594 }
595
596
597 /*
598 * Allocate a new complex number.
599 */
600 COMPLEX *
comalloc(void)601 comalloc(void)
602 {
603 COMPLEX *r;
604
605 r = (COMPLEX *) malloc(sizeof(COMPLEX));
606 if (r == NULL) {
607 math_error("Cannot allocate complex number");
608 /*NOTREACHED*/
609 }
610 r->links = 1;
611 r->real = qlink(&_qzero_);
612 r->imag = qlink(&_qzero_);
613 return r;
614 }
615
616
617 /*
618 * Free a complex number.
619 */
620 void
comfree(COMPLEX * c)621 comfree(COMPLEX *c)
622 {
623 if (--(c->links) > 0)
624 return;
625 qfree(c->real);
626 qfree(c->imag);
627 free(c);
628 }
629