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