1 /***************************************************************************
2 JSPICE3 adaptation of Spice3e2 - Copyright (c) Stephen R. Whiteley 1992
3 Copyright 1990 Regents of the University of California. All rights reserved.
4 Authors: 1985 Wayne A. Christopher
5 1992 Stephen R. Whiteley
6 ****************************************************************************/
7
8 /*
9 * Routines to do complex mathematical functions. These routines require
10 * the -lm libraries. We sacrifice a lot of space to be able
11 * to avoid having to do a seperate call for every vector element,
12 * but it pays off in time savings. These routines should never
13 * allow FPE's to happen.
14 *
15 * Complex functions are called as follows:
16 * cx_something(data, type, length, &newlength, &newtype),
17 * and return a char * that is cast to complex or double.
18 */
19
20 #include "spice.h"
21 #include "ftedefs.h"
22 #include "ftecmath.h"
23
24
25 char *
cx_tan(data,type,length,newlength,newtype)26 cx_tan(data, type, length, newlength, newtype)
27
28 char *data;
29 short type;
30 int length;
31 int *newlength;
32 short *newtype;
33 {
34 double *d, u, v;
35 complex *c;
36 double *dd = (double *) data;
37 complex *cc = (complex *) data;
38 int i;
39
40 if (type == VF_REAL) {
41 d = alloc_d(length);
42 *newtype = VF_REAL;
43 }
44 else {
45 c = alloc_c(length);
46 *newtype = VF_COMPLEX;
47 }
48 *newlength = length;
49 if (type == VF_COMPLEX) {
50 for (i = 0; i < length; i++) {
51 rcheck(cos(degtorad(realpart(&cc[i])))*cosh(degtorad(imagpart(&cc[i]))),"tan");
52 rcheck(sin(degtorad(realpart(&cc[i])))*sinh(degtorad(imagpart(&cc[i]))),"tan");
53 u = degtorad(realpart(&cc[i]));
54 v = degtorad(imagpart(&cc[i]));
55 /* The Lattice C compiler won't take multi-line macros, and
56 * CMS won't take >80 column lines....
57 */
58 #define xx1 sin(u) * cosh(v)
59 #define xx2 cos(u) * sinh(v)
60 #define xx3 cos(u) * cosh(v)
61 #define xx4 sin(u) * sinh(v)
62 cdiv(xx1, xx2, xx3, xx4, realpart(&c[i]), imagpart(&c[i]));
63 }
64 return ((char *) c);
65 }
66 else {
67 for (i = 0; i < length; i++) {
68 rcheck(cos(degtorad(dd[i])) != 0, "tan");
69 d[i] = sin(degtorad(dd[i])) / cos(degtorad(dd[i]));
70 }
71 return ((char *) d);
72 }
73 }
74
75
76 char *
cx_atan(data,type,length,newlength,newtype)77 cx_atan(data, type, length, newlength, newtype)
78
79 char *data;
80 short type;
81 int length;
82 int *newlength;
83 short *newtype;
84 {
85 double *d;
86 double *dd = (double *) data;
87 complex *cc = (complex *) data;
88 int i;
89
90 d = alloc_d(length);
91 *newtype = VF_REAL;
92 *newlength = length;
93 if (type == VF_COMPLEX) {
94 for (i = 0; i < length; i++)
95 d[i] = radtodeg(atan(realpart(&cc[i])));
96 }
97 else {
98 for (i = 0; i < length; i++)
99 d[i] = radtodeg(atan(dd[i]));
100 }
101 return ((char *) d);
102 }
103
104 /* Normalize the data so that the magnitude of the greatest value is 1. */
105
106 char *
cx_norm(data,type,length,newlength,newtype)107 cx_norm(data, type, length, newlength, newtype)
108
109 char *data;
110 short type;
111 int length;
112 int *newlength;
113 short *newtype;
114 {
115 double *d, largest = 0.0;
116 complex *c;
117 double *dd = (double *) data;
118 complex *cc = (complex *) data;
119 int i;
120
121 if (type == VF_REAL) {
122 d = alloc_d(length);
123 *newtype = VF_REAL;
124 }
125 else {
126 c = alloc_c(length);
127 *newtype = VF_COMPLEX;
128 }
129 *newlength = length;
130 if (type == VF_COMPLEX) {
131 for (i = 0; i < length; i++)
132 if (cmag(&cc[i]) > largest)
133 largest = cmag(&cc[i]);
134 }
135 else {
136 for (i = 0; i < length; i++)
137 if (FTEcabs(dd[i]) > largest)
138 largest = FTEcabs(dd[i]);
139 }
140 if (largest == 0.0) {
141 fprintf(cp_err, "Error: can't normalize a 0 vector\n");
142 return (NULL);
143 }
144 if (type == VF_COMPLEX) {
145 for (i = 0; i < length; i++) {
146 realpart(&c[i]) = realpart(&cc[i]) / largest;
147 imagpart(&c[i]) = imagpart(&cc[i]) / largest;
148 }
149 return ((char *) c);
150 }
151 else {
152 for (i = 0; i < length; i++)
153 d[i] = dd[i] / largest;
154 return ((char *) d);
155 }
156 }
157
158
159 char *
cx_uminus(data,type,length,newlength,newtype)160 cx_uminus(data, type, length, newlength, newtype)
161
162 char *data;
163 short type;
164 int length;
165 int *newlength;
166 short *newtype;
167 {
168 double *d;
169 complex *c;
170 double *dd = (double *) data;
171 complex *cc = (complex *) data;
172 int i;
173
174 if (type == VF_REAL) {
175 d = alloc_d(length);
176 *newtype = VF_REAL;
177 }
178 else {
179 c = alloc_c(length);
180 *newtype = VF_COMPLEX;
181 }
182 *newlength = length;
183 if (type == VF_COMPLEX) {
184 for (i = 0; i < length; i++) {
185 realpart(&c[i]) = - realpart(&cc[i]);
186 imagpart(&c[i]) = - imagpart(&cc[i]);
187 }
188 return ((char *) c);
189 }
190 else {
191 for (i = 0; i < length; i++)
192 d[i] = - dd[i];
193 return ((char *) d);
194 }
195 }
196
197 /* Compute the mean of a vector. */
198
199 char *
cx_mean(data,type,length,newlength,newtype)200 cx_mean(data, type, length, newlength, newtype)
201
202 char *data;
203 short type;
204 int length;
205 int *newlength;
206 short *newtype;
207 {
208 complex *c;
209 double *d;
210 complex *cc = (complex *) data;
211 double *dd = (double *) data;
212 int i;
213
214 *newlength = 1;
215 rcheck(length > 0, "mean");
216 if (type == VF_REAL) {
217 d = alloc_d(1);
218 *newtype = VF_REAL;
219 for (i = 0; i < length; i++)
220 *d += dd[i];
221 *d /= length;
222 return ((char *) d);
223 }
224 else {
225 c = alloc_c(1);
226 *newtype = VF_COMPLEX;
227 for (i = 0; i < length; i++) {
228 realpart(c) += realpart(cc + i);
229 imagpart(c) += imagpart(cc + i);
230 }
231 realpart(c) /= length;
232 imagpart(c) /= length;
233 return ((char *) c);
234 }
235 }
236
237 /* ARGSUSED */
238 char *
cx_length(data,type,length,newlength,newtype)239 cx_length(data, type, length, newlength, newtype)
240
241 char *data;
242 short type;
243 int length;
244 int *newlength;
245 short *newtype;
246 {
247 double *d;
248
249 *newlength = 1;
250 *newtype = VF_REAL;
251 d = alloc_d(1);
252 *d = length;
253 return ((char *) d);
254 }
255
256 /* Return a vector from 0 to the magnitude of the argument. Length of the
257 * argument is irrelevent.
258 */
259
260 /* ARGSUSED */
261 char *
cx_vector(data,type,length,newlength,newtype)262 cx_vector(data, type, length, newlength, newtype)
263
264 char *data;
265 short type;
266 int length;
267 int *newlength;
268 short *newtype;
269 {
270 complex *cc = (complex *) data;
271 double *dd = (double *) data;
272 int i, len;
273 double *d;
274
275 if (type == VF_REAL)
276 len = FTEcabs(*dd);
277 else
278 len = cmag(cc);
279 if (len == 0)
280 len = 1;
281 d = alloc_d(len);
282 *newlength = len;
283 *newtype = VF_REAL;
284 for (i = 0; i < len; i++)
285 d[i] = i;
286 return ((char *) d);
287 }
288
289 /* Create a vector of the given length composed of all ones. */
290
291 /* ARGSUSED */
292 char *
cx_unitvec(data,type,length,newlength,newtype)293 cx_unitvec(data, type, length, newlength, newtype)
294
295 char *data;
296 short type;
297 int length;
298 int *newlength;
299 short *newtype;
300 {
301 complex *cc = (complex *) data;
302 double *dd = (double *) data;
303 int i, len;
304 double *d;
305
306 if (type == VF_REAL)
307 len = FTEcabs(*dd);
308 else
309 len = cmag(cc);
310 if (len == 0)
311 len = 1;
312 d = alloc_d(len);
313 *newlength = len;
314 *newtype = VF_REAL;
315 for (i = 0; i < len; i++)
316 d[i] = 1;
317 return ((char *) d);
318 }
319
320 /* Calling methods for these functions are:
321 * cx_something(data1, data2, datatype1, datatype2, length)
322 * The length of the two data vectors is always the same, and is the length
323 * of the result. The result type is complex iff one of the args is
324 * complex.
325 */
326
327 char *
cx_plus(data1,data2,datatype1,datatype2,length)328 cx_plus(data1, data2, datatype1, datatype2, length)
329
330 char *data1, *data2;
331 short datatype1, datatype2;
332 {
333 double *dd1 = (double *) data1;
334 double *dd2 = (double *) data2;
335 double *d;
336 complex *cc1 = (complex *) data1;
337 complex *cc2 = (complex *) data2;
338 complex *c, c1, c2;
339 int i;
340
341 if ((datatype1 == VF_REAL) && (datatype2 == VF_REAL)) {
342 d = alloc_d(length);
343 for (i = 0; i < length; i++)
344 d[i] = dd1[i] + dd2[i];
345 return ((char *) d);
346 }
347 else {
348 c = alloc_c(length);
349 for (i = 0; i < length; i++) {
350 if (datatype1 == VF_REAL) {
351 realpart(&c1) = dd1[i];
352 imagpart(&c1) = 0.0;
353 } else {
354 realpart(&c1) = realpart(&cc1[i]);
355 imagpart(&c1) = imagpart(&cc1[i]);
356 }
357 if (datatype2 == VF_REAL) {
358 realpart(&c2) = dd2[i];
359 imagpart(&c2) = 0.0;
360 } else {
361 realpart(&c2) = realpart(&cc2[i]);
362 imagpart(&c2) = imagpart(&cc2[i]);
363 }
364 realpart(&c[i]) = realpart(&c1) + realpart(&c2);
365 imagpart(&c[i]) = imagpart(&c1) + imagpart(&c2);
366 }
367 return ((char *) c);
368 }
369 }
370
371
372 char *
cx_minus(data1,data2,datatype1,datatype2,length)373 cx_minus(data1, data2, datatype1, datatype2, length)
374
375 char *data1, *data2;
376 short datatype1, datatype2;
377 {
378 double *dd1 = (double *) data1;
379 double *dd2 = (double *) data2;
380 double *d;
381 complex *cc1 = (complex *) data1;
382 complex *cc2 = (complex *) data2;
383 complex *c, c1, c2;
384 int i;
385
386 if ((datatype1 == VF_REAL) && (datatype2 == VF_REAL)) {
387 d = alloc_d(length);
388 for (i = 0; i < length; i++)
389 d[i] = dd1[i] - dd2[i];
390 return ((char *) d);
391 }
392 else {
393 c = alloc_c(length);
394 for (i = 0; i < length; i++) {
395 if (datatype1 == VF_REAL) {
396 realpart(&c1) = dd1[i];
397 imagpart(&c1) = 0.0;
398 }
399 else {
400 realpart(&c1) = realpart(&cc1[i]);
401 imagpart(&c1) = imagpart(&cc1[i]);
402 }
403 if (datatype2 == VF_REAL) {
404 realpart(&c2) = dd2[i];
405 imagpart(&c2) = 0.0;
406 }
407 else {
408 realpart(&c2) = realpart(&cc2[i]);
409 imagpart(&c2) = imagpart(&cc2[i]);
410 }
411 realpart(&c[i]) = realpart(&c1) - realpart(&c2);
412 imagpart(&c[i]) = imagpart(&c1) - imagpart(&c2);
413 }
414 return ((char *) c);
415 }
416 }
417
418
419 char *
cx_times(data1,data2,datatype1,datatype2,length)420 cx_times(data1, data2, datatype1, datatype2, length)
421
422 char *data1, *data2;
423 short datatype1, datatype2;
424 {
425 double *dd1 = (double *) data1;
426 double *dd2 = (double *) data2;
427 double *d;
428 complex *cc1 = (complex *) data1;
429 complex *cc2 = (complex *) data2;
430 complex *c, c1, c2;
431 int i;
432
433 if ((datatype1 == VF_REAL) && (datatype2 == VF_REAL)) {
434 d = alloc_d(length);
435 for (i = 0; i < length; i++)
436 d[i] = dd1[i] * dd2[i];
437 return ((char *) d);
438 }
439 else {
440 c = alloc_c(length);
441 for (i = 0; i < length; i++) {
442 if (datatype1 == VF_REAL) {
443 realpart(&c1) = dd1[i];
444 imagpart(&c1) = 0.0;
445 } else {
446 realpart(&c1) = realpart(&cc1[i]);
447 imagpart(&c1) = imagpart(&cc1[i]);
448 }
449 if (datatype2 == VF_REAL) {
450 realpart(&c2) = dd2[i];
451 imagpart(&c2) = 0.0;
452 } else {
453 realpart(&c2) = realpart(&cc2[i]);
454 imagpart(&c2) = imagpart(&cc2[i]);
455 }
456 realpart(&c[i]) = realpart(&c1) * realpart(&c2)
457 - imagpart(&c1) * imagpart(&c2);
458 imagpart(&c[i]) = imagpart(&c1) * realpart(&c2)
459 + realpart(&c1) * imagpart(&c2);
460 }
461 return ((char *) c);
462 }
463 }
464
465
466 char *
cx_mod(data1,data2,datatype1,datatype2,length)467 cx_mod(data1, data2, datatype1, datatype2, length)
468
469 char *data1, *data2;
470 short datatype1, datatype2;
471 {
472 double *dd1 = (double *) data1;
473 double *dd2 = (double *) data2;
474 double *d;
475 complex *cc1 = (complex *) data1;
476 complex *cc2 = (complex *) data2;
477 complex *c, c1, c2;
478 int i, r1, r2, i1, i2, r3, i3;
479
480 if ((datatype1 == VF_REAL) && (datatype2 == VF_REAL)) {
481 d = alloc_d(length);
482 for (i = 0; i < length; i++) {
483 r1 = floor(FTEcabs(dd1[i]));
484 rcheck(r1 > 0, "mod");
485 r2 = floor(FTEcabs(dd2[i]));
486 rcheck(r2 > 0, "mod");
487 r3 = r1 % r2;
488 d[i] = (double) r3;
489 }
490 return ((char *) d);
491 }
492 else {
493 c = alloc_c(length);
494 for (i = 0; i < length; i++) {
495 if (datatype1 == VF_REAL) {
496 realpart(&c1) = dd1[i];
497 imagpart(&c1) = 0.0;
498 } else {
499 realpart(&c1) = realpart(&cc1[i]);
500 imagpart(&c1) = imagpart(&cc1[i]);
501 }
502 if (datatype2 == VF_REAL) {
503 realpart(&c2) = dd2[i];
504 imagpart(&c2) = 0.0;
505 } else {
506 realpart(&c2) = realpart(&cc2[i]);
507 imagpart(&c2) = imagpart(&cc2[i]);
508 }
509 r1 = floor(FTEcabs(realpart(&c1)));
510 rcheck(r1 > 0, "mod");
511 r2 = floor(FTEcabs(realpart(&c2)));
512 rcheck(r2 > 0, "mod");
513 i1 = floor(FTEcabs(imagpart(&c1)));
514 rcheck(i1 > 0, "mod");
515 i2 = floor(FTEcabs(imagpart(&c2)));
516 rcheck(i2 > 0, "mod");
517 r3 = r1 % r2;
518 i3 = i1 % i2;
519 realpart(&c[i]) = (double) r3;
520 imagpart(&c[i]) = (double) i3;
521 }
522 return ((char *) c);
523 }
524 }
525
526
527 #ifdef notdef
528
529 /* this is stupid */
530 char *
cx_rnd(data,type,length,newlength,newtype)531 cx_rnd(data, type, length, newlength, newtype)
532
533 char *data;
534 short type;
535 int length;
536 int *newlength;
537 short *newtype;
538 {
539 double *d;
540 complex *c;
541 int i, j, k;
542 double *dd = (double *) data;
543 complex *cc = (complex *) data;
544
545 if (type == VF_REAL) {
546 d = alloc_d(length);
547 *newtype = VF_REAL;
548 } else {
549 c = alloc_c(length);
550 *newtype = VF_COMPLEX;
551 }
552 *newlength = length;
553 if (type == VF_COMPLEX) {
554 for (i = 0; i < length; i++) {
555 j = floor(realpart(&cc[i]));
556 k = floor(imagpart(&cc[i]));
557 realpart(&c[i]) = j ? random() % j : 0;
558 imagpart(&c[i]) = k ? random() % k : 0;
559 }
560 return ((char *) c);
561 }else {
562 for (i = 0; i < length; i++) {
563 j = floor(dd[i]);
564 d[i] = j ? random() % j : 0;
565 }
566 return ((char *) d);
567 }
568 }
569 #endif
570
571
572 char *
cx_rnd(data,type,length,newlength,newtype)573 cx_rnd(data, type, length, newlength, newtype)
574
575 /* Return a vector of random values between zero and the number given
576 * in the input vector. If the input vector is complex, so is the
577 * vector produced.
578 */
579
580 char *data;
581 short type;
582 int length;
583 int *newlength;
584 short *newtype;
585 {
586 double *d;
587 complex *c;
588 int i, j, k;
589 double *dd = (double *) data;
590 complex *cc = (complex *) data;
591
592 if (type == VF_REAL) {
593 d = alloc_d(length);
594 *newtype = VF_REAL;
595 }
596 else {
597 c = alloc_c(length);
598 *newtype = VF_COMPLEX;
599 }
600 *newlength = length;
601 if (type == VF_COMPLEX) {
602 for (i = 0; i < length; i++) {
603 realpart(&c[i]) = xrandom()*realpart(&cc[i]);
604 imagpart(&c[i]) = xrandom()*imagpart(&cc[i]);
605 }
606 return ((char *) c);
607 }
608 else {
609 for (i = 0; i < length; i++) {
610 d[i] = xrandom()*dd[i];
611 }
612 return ((char *) d);
613 }
614 }
615
616
617 char *
cx_gauss(data,type,length,newlength,newtype)618 cx_gauss(data, type, length, newlength, newtype)
619
620 /* Return a real vector with normal distributed coefficients.
621 * If the given vector is complex, terms are interpreted as (sd,mean).
622 * Otherwise the returned coefficients have sd given by the term
623 * in the real vector given, with zero mean.
624 */
625
626 char *data;
627 short type;
628 int length;
629 int *newlength;
630 short *newtype;
631 {
632 double *d;
633 complex *c;
634 int i, j, k;
635 double *dd = (double *) data;
636 complex *cc = (complex *) data;
637
638 d = alloc_d(length);
639 *newtype = VF_REAL;
640 *newlength = length;
641
642 if (type == VF_COMPLEX) {
643 for (i = 0; i < length; i++) {
644 d[i] = xgauss()*realpart(&cc[i]) + imagpart(&cc[i]);
645 }
646 }
647 else {
648 for (i = 0; i < length; i++) {
649 d[i] = xgauss()*dd[i];
650 }
651 }
652 return ((char *) d);
653 }
654