1 /**********
2 Copyright 1990 Regents of the University of California. All rights reserved.
3 Author: 1985 Wayne A. Christopher, U. C. Berkeley CAD Group
4 **********/
5
6 /*
7 * Routines to do complex mathematical functions. These routines require
8 * the -lm libraries. We sacrifice a lot of space to be able
9 * to avoid having to do a seperate call for every vector element,
10 * but it pays off in time savings. These routines should never
11 * allow FPE's to happen.
12 *
13 * Complex functions are called as follows:
14 * cx_something(data, type, length, &newlength, &newtype),
15 * and return a char * that is cast to complex or double.
16 */
17
18 #include "ngspice/ngspice.h"
19 #include "ngspice/cpdefs.h"
20 #include "ngspice/dvec.h"
21 #include "ngspice/randnumb.h"
22
23 #include "cmath.h"
24 #include "cmath2.h"
25
26
27 static double
cx_max_local(void * data,short int type,int length)28 cx_max_local(void *data, short int type, int length)
29 {
30 double largest = 0.0;
31
32 if (type == VF_COMPLEX) {
33 ngcomplex_t *cc = (ngcomplex_t *) data;
34 int i;
35
36 for (i = 0; i < length; i++)
37 if (largest < cmag(cc[i]))
38 largest = cmag(cc[i]);
39 } else {
40 double *dd = (double *) data;
41 int i;
42
43 for (i = 0; i < length; i++)
44 if (largest < fabs(dd[i]))
45 largest = fabs(dd[i]);
46 }
47 return largest;
48 }
49
50
51 /* Normalize the data so that the magnitude of the greatest value is 1. */
52
53 void *
cx_norm(void * data,short int type,int length,int * newlength,short int * newtype)54 cx_norm(void *data, short int type, int length, int *newlength, short int *newtype)
55 {
56 double largest = 0.0;
57
58 largest = cx_max_local(data, type, length);
59 if (largest == 0.0) {
60 fprintf(cp_err, "Error: can't normalize a 0 vector\n");
61 return (NULL);
62 }
63
64 *newlength = length;
65 if (type == VF_COMPLEX) {
66 ngcomplex_t *c;
67 ngcomplex_t *cc = (ngcomplex_t *) data;
68 int i;
69
70 c = alloc_c(length);
71 *newtype = VF_COMPLEX;
72
73 for (i = 0; i < length; i++) {
74 realpart(c[i]) = realpart(cc[i]) / largest;
75 imagpart(c[i]) = imagpart(cc[i]) / largest;
76 }
77 return ((void *) c);
78 } else {
79 double *d;
80 double *dd = (double *) data;
81 int i;
82
83 d = alloc_d(length);
84 *newtype = VF_REAL;
85
86 for (i = 0; i < length; i++)
87 d[i] = dd[i] / largest;
88 return ((void *) d);
89 }
90 }
91
92
93 void *
cx_uminus(void * data,short int type,int length,int * newlength,short int * newtype)94 cx_uminus(void *data, short int type, int length, int *newlength, short int *newtype)
95 {
96 *newlength = length;
97 if (type == VF_COMPLEX) {
98 ngcomplex_t *c;
99 ngcomplex_t *cc = (ngcomplex_t *) data;
100 int i;
101
102 c = alloc_c(length);
103 *newtype = VF_COMPLEX;
104 for (i = 0; i < length; i++) {
105 realpart(c[i]) = - realpart(cc[i]);
106 imagpart(c[i]) = - imagpart(cc[i]);
107 }
108 return ((void *) c);
109 } else {
110 double *d;
111 double *dd = (double *) data;
112 int i;
113
114 d = alloc_d(length);
115 *newtype = VF_REAL;
116 for (i = 0; i < length; i++)
117 d[i] = - dd[i];
118 return ((void *) d);
119 }
120 }
121
122
123 /* random integers drawn from a uniform distribution
124 * data in: integer numbers, their absolut values are used,
125 * maximum is RAND_MAX (32767)
126 * data out: random integers in interval [0, data[i][
127 * standard library function rand() is used
128 */
129
130 void *
cx_rnd(void * data,short int type,int length,int * newlength,short int * newtype)131 cx_rnd(void *data, short int type, int length, int *newlength, short int *newtype)
132 {
133 *newlength = length;
134 checkseed();
135 if (type == VF_COMPLEX) {
136 ngcomplex_t *c;
137 ngcomplex_t *cc = (ngcomplex_t *) data;
138 int i;
139
140 c = alloc_c(length);
141 *newtype = VF_COMPLEX;
142 for (i = 0; i < length; i++) {
143 int j, k;
144 j = (int)floor(realpart(cc[i]));
145 k = (int)floor(imagpart(cc[i]));
146 realpart(c[i]) = j ? rand() % j : 0;
147 imagpart(c[i]) = k ? rand() % k : 0;
148 }
149 return ((void *) c);
150 } else {
151 double *d;
152 double *dd = (double *) data;
153 int i;
154
155 d = alloc_d(length);
156 *newtype = VF_REAL;
157 for (i = 0; i < length; i++) {
158 int j;
159 j = (int)floor(dd[i]);
160 d[i] = j ? rand() % j : 0;
161 }
162 return ((void *) d);
163 }
164 }
165
166
167 /* random numbers drawn from a uniform distribution
168 * data out: random numbers in interval [-1, 1[
169 */
170
171 void *
cx_sunif(void * data,short int type,int length,int * newlength,short int * newtype)172 cx_sunif(void *data, short int type, int length, int *newlength, short int *newtype)
173 {
174 NG_IGNORE(data);
175
176 *newlength = length;
177 checkseed();
178 if (type == VF_COMPLEX) {
179 ngcomplex_t *c;
180 int i;
181
182 c = alloc_c(length);
183 *newtype = VF_COMPLEX;
184 for (i = 0; i < length; i++) {
185 realpart(c[i]) = drand();
186 imagpart(c[i]) = drand();
187 }
188 return ((void *) c);
189 } else {
190 double *d;
191 int i;
192
193 d = alloc_d(length);
194 *newtype = VF_REAL;
195 for (i = 0; i < length; i++) {
196 d[i] = drand();
197 }
198 return ((void *) d);
199 }
200 }
201
202
203 /* random numbers drawn from a poisson distribution
204 * data in: lambda
205 * data out: random integers according to poisson distribution,
206 * with lambda given by each vector element
207 */
208
209 void *
cx_poisson(void * data,short int type,int length,int * newlength,short int * newtype)210 cx_poisson(void *data, short int type, int length, int *newlength, short int *newtype)
211 {
212 *newlength = length;
213 checkseed();
214 if (type == VF_COMPLEX) {
215 ngcomplex_t *c;
216 ngcomplex_t *cc = (ngcomplex_t *) data;
217 int i;
218
219 c = alloc_c(length);
220 *newtype = VF_COMPLEX;
221 for (i = 0; i < length; i++) {
222 realpart(c[i]) = poisson (realpart(cc[i]));
223 imagpart(c[i]) = poisson (imagpart(cc[i]));
224 }
225 return ((void *) c);
226 } else {
227 double *d;
228 double *dd = (double *) data;
229 int i;
230
231 d = alloc_d(length);
232 *newtype = VF_REAL;
233 for (i = 0; i < length; i++) {
234 d[i] = poisson(dd[i]);
235 }
236 return ((void *) d);
237 }
238 }
239
240
241 /* random numbers drawn from an exponential distribution
242 * data in: Mean values
243 * data out: exponentially distributed random numbers,
244 * with mean given by each vector element
245 */
246
247 void *
cx_exponential(void * data,short int type,int length,int * newlength,short int * newtype)248 cx_exponential(void *data, short int type, int length, int *newlength, short int *newtype)
249 {
250 *newlength = length;
251 checkseed();
252 if (type == VF_COMPLEX) {
253 ngcomplex_t *c;
254 ngcomplex_t *cc = (ngcomplex_t *) data;
255 int i;
256
257 c = alloc_c(length);
258 *newtype = VF_COMPLEX;
259 for (i = 0; i < length; i++) {
260 realpart(c[i]) = exprand(realpart(cc[i]));
261 imagpart(c[i]) = exprand(imagpart(cc[i]));
262 }
263 return ((void *) c);
264 } else {
265 double *d;
266 double *dd = (double *) data;
267 int i;
268
269 d = alloc_d(length);
270 *newtype = VF_REAL;
271 for (i = 0; i < length; i++) {
272 d[i] = exprand(dd[i]);
273 }
274 return ((void *) d);
275 }
276 }
277
278
279 /* random numbers drawn from a Gaussian distribution
280 mean 0, std dev 1
281 */
282
283 void *
cx_sgauss(void * data,short int type,int length,int * newlength,short int * newtype)284 cx_sgauss(void *data, short int type, int length, int *newlength, short int *newtype)
285 {
286 NG_IGNORE(data);
287
288 *newlength = length;
289 checkseed();
290 if (type == VF_COMPLEX) {
291 ngcomplex_t *c;
292 int i;
293
294 c = alloc_c(length);
295 *newtype = VF_COMPLEX;
296 for (i = 0; i < length; i++) {
297 realpart(c[i]) = gauss0();
298 imagpart(c[i]) = gauss0();
299 }
300 return ((void *) c);
301 } else {
302 double *d;
303 int i;
304
305 d = alloc_d(length);
306 *newtype = VF_REAL;
307 for (i = 0; i < length; i++) {
308 d[i] = gauss1();
309 }
310 return ((void *) d);
311 }
312 }
313
314
315 /* Compute the avg of a vector.
316 * Created by A.M.Roldan 2005-05-21
317 */
318
319 void *
cx_avg(void * data,short int type,int length,int * newlength,short int * newtype)320 cx_avg(void *data, short int type, int length, int *newlength, short int *newtype)
321 {
322 double sum_real = 0.0, sum_imag = 0.0;
323 int i;
324
325 if (type == VF_REAL) {
326
327 double *d = alloc_d(length);
328 double *dd = (double *) data;
329
330 *newtype = VF_REAL;
331 *newlength = length;
332
333 for (i = 0; i < length; i++) {
334 sum_real += dd[i];
335 d[i] = sum_real / ((double) i + 1.0);
336 }
337
338 return ((void *) d);
339
340 } else {
341
342 ngcomplex_t *c = alloc_c(length);
343 ngcomplex_t *cc = (ngcomplex_t *) data;
344
345 *newtype = VF_COMPLEX;
346 *newlength = length;
347
348 for (i = 0; i < length; i++) {
349 sum_real += realpart(cc[i]);
350 realpart(c[i]) = sum_real / ((double) i + 1.0);
351
352 sum_imag += imagpart(cc[i]);
353 imagpart(c[i]) = sum_imag / ((double) i + 1.0);
354 }
355
356 return ((void *) c);
357
358 }
359 }
360
361
362 /* Compute the mean of a vector. */
363
cx_mean(void * data,short int type,int length,int * newlength,short int * newtype)364 void *cx_mean(void *data, short int type, int length,
365 int *newlength, short int *newtype)
366 {
367 if (length == 0) {
368 (void) fprintf(cp_err, "mean calculation requires "
369 "at least one element.\n");
370 return NULL;
371 }
372
373 *newlength = 1;
374
375 if (type == VF_REAL) {
376 double *d;
377 double *dd = (double *) data;
378 int i;
379
380 d = alloc_d(1);
381 *newtype = VF_REAL;
382 for (i = 0; i < length; i++)
383 *d += dd[i];
384 *d /= length;
385 return ((void *) d);
386 }
387 else {
388 ngcomplex_t *c;
389 ngcomplex_t *cc = (ngcomplex_t *) data;
390 int i;
391
392 c = alloc_c(1);
393 *newtype = VF_COMPLEX;
394 for (i = 0; i < length; i++) {
395 realpart(*c) += realpart(cc[i]);
396 imagpart(*c) += imagpart(cc[i]);
397 }
398 realpart(*c) /= length;
399 imagpart(*c) /= length;
400 return ((void *) c);
401 }
402 }
403
404
405 /* Compute the standard deviation of all elements of a vector. */
406
cx_stddev(void * data,short int type,int length,int * newlength,short int * newtype)407 void *cx_stddev(void *data, short int type, int length,
408 int *newlength, short int *newtype)
409 {
410 if (length == 0) {
411 (void) fprintf(cp_err, "standard deviation calculation requires "
412 "at least one element.\n");
413 return NULL;
414 }
415
416 *newlength = 1;
417
418 if (type == VF_REAL) {
419 double *p_mean = (double *) cx_mean(data, type, length,
420 newlength, newtype);
421 double mean = *p_mean;
422 double *d, sum = 0.0;
423 double *dd = (double *) data;
424
425 d = alloc_d(1);
426 *newtype = VF_REAL;
427 int i;
428 for (i = 0; i < length; i++) {
429 const double tmp = dd[i] - mean;
430 sum += tmp * tmp;
431 }
432
433 *d = sqrt(sum / ((double) length - 1.0));
434 txfree(p_mean);
435 return (void *) d;
436 }
437 else {
438 ngcomplex_t * const p_cmean = (ngcomplex_t *) cx_mean(data, type, length,
439 newlength, newtype);
440 const ngcomplex_t cmean = *p_cmean;
441 const double rmean = realpart(cmean);
442 const double imean = imagpart(cmean);
443 double *d, sum = 0.0;
444 ngcomplex_t *cc = (ngcomplex_t *) data;
445
446 d = alloc_d(1);
447 *newtype = VF_REAL;
448 int i;
449 for (i = 0; i < length; i++) {
450 const double a = realpart(cc[i]) - rmean;
451 const double b = imagpart(cc[i]) - imean;
452 sum += a * a + b * b;
453 }
454 *d = sqrt(sum / ((double) length - 1.0));
455 txfree(p_cmean);
456 return (void *) d;
457 }
458 } /* end of function cx_stddev */
459
460
461
462 void *
cx_length(void * data,short int type,int length,int * newlength,short int * newtype)463 cx_length(void *data, short int type, int length, int *newlength, short int *newtype)
464 {
465 double *d;
466
467 NG_IGNORE(data);
468 NG_IGNORE(type);
469
470 *newlength = 1;
471 *newtype = VF_REAL;
472 d = alloc_d(1);
473 *d = length;
474 return ((void *) d);
475 }
476
477
478 /* Return a vector from 0 to the magnitude of the argument. Length of the
479 * argument is irrelevent.
480 */
481
482 void *
cx_vector(void * data,short int type,int length,int * newlength,short int * newtype)483 cx_vector(void *data, short int type, int length, int *newlength, short int *newtype)
484 {
485 ngcomplex_t *cc = (ngcomplex_t *) data;
486 double *dd = (double *) data;
487 int i, len;
488 double *d;
489
490 NG_IGNORE(length);
491
492 if (type == VF_REAL)
493 len = (int)fabs(*dd);
494 else
495 len = (int)cmag(*cc);
496 if (len == 0)
497 len = 1;
498 d = alloc_d(len);
499 *newlength = len;
500 *newtype = VF_REAL;
501 for (i = 0; i < len; i++)
502 d[i] = i;
503 return ((void *) d);
504 }
505
506
507 /* Create a vector of the given length composed of all ones. */
508
509 void *
cx_unitvec(void * data,short int type,int length,int * newlength,short int * newtype)510 cx_unitvec(void *data, short int type, int length, int *newlength, short int *newtype)
511 {
512 ngcomplex_t *cc = (ngcomplex_t *) data;
513 double *dd = (double *) data;
514 int i, len;
515 double *d;
516
517 NG_IGNORE(length);
518
519 if (type == VF_REAL)
520 len = (int)fabs(*dd);
521 else
522 len = (int)cmag(*cc);
523 if (len == 0)
524 len = 1;
525 d = alloc_d(len);
526 *newlength = len;
527 *newtype = VF_REAL;
528 for (i = 0; i < len; i++)
529 d[i] = 1;
530 return ((void *) d);
531 }
532
533
534 /* Calling methods for these functions are:
535 * cx_something(data1, data2, datatype1, datatype2, length)
536 *
537 * The length of the two data vectors is always the same, and is the length
538 * of the result. The result type is complex iff one of the args is
539 * complex.
540 */
541
542 void *
cx_plus(void * data1,void * data2,short int datatype1,short int datatype2,int length)543 cx_plus(void *data1, void *data2, short int datatype1, short int datatype2, int length)
544 {
545 double *dd1 = (double *) data1;
546 double *dd2 = (double *) data2;
547 double *d;
548 ngcomplex_t *cc1 = (ngcomplex_t *) data1;
549 ngcomplex_t *cc2 = (ngcomplex_t *) data2;
550 ngcomplex_t *c, c1, c2;
551 int i;
552
553 if ((datatype1 == VF_REAL) && (datatype2 == VF_REAL)) {
554 d = alloc_d(length);
555 for (i = 0; i < length; i++)
556 d[i] = dd1[i] + dd2[i];
557 return ((void *) d);
558 } else {
559 c = alloc_c(length);
560 for (i = 0; i < length; i++) {
561 if (datatype1 == VF_REAL) {
562 realpart(c1) = dd1[i];
563 imagpart(c1) = 0.0;
564 } else {
565 c1 = cc1[i];
566 }
567 if (datatype2 == VF_REAL) {
568 realpart(c2) = dd2[i];
569 imagpart(c2) = 0.0;
570 } else {
571 c2 = cc2[i];
572 }
573 realpart(c[i]) = realpart(c1) + realpart(c2);
574 imagpart(c[i]) = imagpart(c1) + imagpart(c2);
575 }
576 return ((void *) c);
577 }
578 }
579
580
581 void *
cx_minus(void * data1,void * data2,short int datatype1,short int datatype2,int length)582 cx_minus(void *data1, void *data2, short int datatype1, short int datatype2, int length)
583 {
584 double *dd1 = (double *) data1;
585 double *dd2 = (double *) data2;
586 double *d;
587 ngcomplex_t *cc1 = (ngcomplex_t *) data1;
588 ngcomplex_t *cc2 = (ngcomplex_t *) data2;
589 ngcomplex_t *c, c1, c2;
590 int i;
591
592 if ((datatype1 == VF_REAL) && (datatype2 == VF_REAL)) {
593 d = alloc_d(length);
594 for (i = 0; i < length; i++)
595 d[i] = dd1[i] - dd2[i];
596 return ((void *) d);
597 } else {
598 c = alloc_c(length);
599 for (i = 0; i < length; i++) {
600 if (datatype1 == VF_REAL) {
601 realpart(c1) = dd1[i];
602 imagpart(c1) = 0.0;
603 } else {
604 c1 = cc1[i];
605 }
606 if (datatype2 == VF_REAL) {
607 realpart(c2) = dd2[i];
608 imagpart(c2) = 0.0;
609 } else {
610 c2 = cc2[i];
611 }
612 realpart(c[i]) = realpart(c1) - realpart(c2);
613 imagpart(c[i]) = imagpart(c1) - imagpart(c2);
614 }
615 return ((void *) c);
616 }
617 }
618
619
620 void *
cx_times(void * data1,void * data2,short int datatype1,short int datatype2,int length)621 cx_times(void *data1, void *data2, short int datatype1, short int datatype2, int length)
622 {
623 double *dd1 = (double *) data1;
624 double *dd2 = (double *) data2;
625 double *d;
626 ngcomplex_t *cc1 = (ngcomplex_t *) data1;
627 ngcomplex_t *cc2 = (ngcomplex_t *) data2;
628 ngcomplex_t *c, c1, c2;
629 int i;
630
631 if ((datatype1 == VF_REAL) && (datatype2 == VF_REAL)) {
632 d = alloc_d(length);
633 for (i = 0; i < length; i++)
634 d[i] = dd1[i] * dd2[i];
635 return ((void *) d);
636 } else {
637 c = alloc_c(length);
638 for (i = 0; i < length; i++) {
639 if (datatype1 == VF_REAL) {
640 realpart(c1) = dd1[i];
641 imagpart(c1) = 0.0;
642 } else {
643 c1 = cc1[i];
644 }
645 if (datatype2 == VF_REAL) {
646 realpart(c2) = dd2[i];
647 imagpart(c2) = 0.0;
648 } else {
649 c2 = cc2[i];
650 }
651 realpart(c[i]) =
652 realpart(c1) * realpart(c2) - imagpart(c1) * imagpart(c2);
653 imagpart(c[i]) =
654 imagpart(c1) * realpart(c2) + realpart(c1) * imagpart(c2);
655 }
656 return ((void *) c);
657 }
658 }
659
660
cx_mod(void * data1,void * data2,short int datatype1,short int datatype2,int length)661 void *cx_mod(void *data1, void *data2, short int datatype1, short int datatype2,
662 int length)
663 {
664 int xrc = 0;
665 void *rv;
666 double *dd1 = (double *) data1;
667 double *dd2 = (double *) data2;
668
669 if ((datatype1 == VF_REAL) && (datatype2 == VF_REAL)) {
670 double *d;
671 rv = d = alloc_d(length);
672
673 int i;
674 for (i = 0; i < length; i++) {
675 const int r1 = (int) floor(fabs(dd1[i]));
676 rcheck(r1 > 0, "mod");
677 const int r2 = (int)floor(fabs(dd2[i]));
678 rcheck(r2 > 0, "mod");
679 const int r3 = r1 % r2;
680 d[i] = (double) r3;
681 }
682 }
683 else {
684 ngcomplex_t *c, c1, c2;
685 ngcomplex_t *cc1 = (ngcomplex_t *) data1;
686 ngcomplex_t *cc2 = (ngcomplex_t *) data2;
687 rv = c = alloc_c(length);
688
689 int i;
690 for (i = 0; i < length; i++) {
691 if (datatype1 == VF_REAL) {
692 realpart(c1) = dd1[i];
693 imagpart(c1) = 0.0;
694 }
695 else {
696 c1 = cc1[i];
697 }
698 if (datatype2 == VF_REAL) {
699 realpart(c2) = dd2[i];
700 imagpart(c2) = 0.0;
701 } else {
702 c2 = cc2[i];
703 }
704 const int r1 = (int) floor(fabs(realpart(c1)));
705 rcheck(r1 > 0, "mod");
706 const int r2 = (int) floor(fabs(realpart(c2)));
707 rcheck(r2 > 0, "mod");
708 const int i1 = (int) floor(fabs(imagpart(c1)));
709 rcheck(i1 > 0, "mod");
710 const int i2 = (int) floor(fabs(imagpart(c2)));
711 rcheck(i2 > 0, "mod");
712 const int r3 = r1 % r2;
713 const int i3 = i1 % i2;
714 realpart(c[i]) = (double) r3;
715 imagpart(c[i]) = (double) i3;
716 }
717 }
718
719 EXITPOINT:
720 if (xrc != 0) { /* Free resources on error */
721 txfree(rv);
722 rv = NULL;
723 }
724
725 return rv;
726 } /* end of function cx_mod */
727
728
729
730 /* Routoure JM : Compute the max of a vector. */
731
cx_max(void * data,short int type,int length,int * newlength,short int * newtype)732 void *cx_max(void *data, short int type, int length,
733 int *newlength, short int *newtype)
734 {
735 if (length == 0) {
736 (void) fprintf(cp_err, "maximum calculation requires "
737 "at least one element.\n");
738 return NULL;
739 }
740
741 *newlength = 1;
742
743 if (type == VF_REAL) {
744 double largest = 0.0;
745 double *d;
746 double *dd = (double *) data;
747 int i;
748
749 d = alloc_d(1);
750 *newtype = VF_REAL;
751 largest = dd[0];
752 for (i = 1; i < length; i++) {
753 const double tmp = dd[i];
754 if (largest < tmp) {
755 largest = tmp;
756 }
757 }
758 *d = largest;
759 return (void *) d;
760 }
761 else {
762 double largest_real = 0.0;
763 double largest_complex = 0.0;
764 ngcomplex_t *c;
765 ngcomplex_t *cc = (ngcomplex_t *) data;
766 int i;
767
768 c = alloc_c(1);
769 *newtype = VF_COMPLEX;
770 largest_real = realpart(*cc);
771 largest_complex = imagpart(*cc);
772 for (i = 1; i < length; i++) {
773 const double tmpr = realpart(cc[i]);
774 if (largest_real < tmpr) {
775 largest_real = tmpr;
776 }
777 const double tmpi = imagpart(cc[i]);
778 if (largest_complex < tmpi) {
779 largest_complex = tmpi;
780 }
781 }
782 realpart(*c) = largest_real;
783 imagpart(*c) = largest_complex;
784 return (void *) c;
785 }
786 } /* end of function cx_max */
787
788
789
790 /* Routoure JM : Compute the min of a vector. */
791
cx_min(void * data,short int type,int length,int * newlength,short int * newtype)792 void *cx_min(void *data, short int type, int length,
793 int *newlength, short int *newtype)
794 {
795 if (length == 0) {
796 (void) fprintf(cp_err, "minimum calculation requires "
797 "at least one element.\n");
798 return NULL;
799 }
800
801 *newlength = 1;
802
803 if (type == VF_REAL) {
804 double smallest;
805 double *d;
806 double *dd = (double *) data;
807 int i;
808
809 d = alloc_d(1);
810 *newtype = VF_REAL;
811 smallest = dd[0];
812 for (i = 1; i < length; i++) {
813 const double tmp = dd[i];
814 if (smallest > tmp) {
815 smallest = tmp;
816 }
817 }
818 *d = smallest;
819 return (void *) d;
820 }
821 else {
822 double smallest_real;
823 double smallest_complex;
824 ngcomplex_t *c;
825 ngcomplex_t *cc = (ngcomplex_t *) data;
826 int i;
827
828 c = alloc_c(1);
829 *newtype = VF_COMPLEX;
830 smallest_real = realpart(*cc);
831 smallest_complex = imagpart(*cc);
832 for (i = 1; i < length; i++) {
833 const double tmpr = realpart(cc[i]);
834 if (smallest_real > tmpr) {
835 smallest_real = tmpr;
836 }
837 const double tmpi = imagpart(cc[i]);
838 if (smallest_complex > tmpi) {
839 smallest_complex = tmpi;
840 }
841 }
842 realpart(*c) = smallest_real;
843 imagpart(*c) = smallest_complex;
844 return (void *) c;
845 }
846 } /* end of function cx_min */
847
848
849 /* Routoure JM : Compute the differential of a vector. */
850
cx_d(void * data,short int type,int length,int * newlength,short int * newtype)851 void *cx_d(void *data, short int type, int length,
852 int *newlength, short int *newtype)
853 {
854 if (length == 0) {
855 (void) fprintf(cp_err, "differential calculation requires "
856 "at least one element.\n");
857 return NULL;
858 }
859
860 *newlength = length;
861
862 if (type == VF_REAL) {
863 double *d;
864 double *dd = (double *) data;
865 int i;
866
867 d = alloc_d(length);
868 *newtype = VF_REAL;
869 d[0] = dd[1] - dd[0];
870 d[length-1] = dd[length-1] - dd[length-2];
871 for (i = 1; i < length - 1; i++) {
872 d[i] = dd[i+1] - dd[i-1];
873 }
874
875 return (void *) d;
876 }
877 else {
878 ngcomplex_t *c;
879 ngcomplex_t *cc = (ngcomplex_t *) data;
880 int i;
881
882 c = alloc_c(length);
883 *newtype = VF_COMPLEX;
884 realpart(*c) = realpart(cc[1]) - realpart(cc[0]);
885 imagpart(*c) = imagpart(cc[1]) - imagpart(cc[0]);
886 realpart(c[length-1]) = realpart(cc[length-1]) - realpart(cc[length-2]);
887 imagpart(c[length-1]) = imagpart(cc[length-1]) - imagpart(cc[length-2]);
888
889
890 for (i = 1; i < length - 1; i++) {
891 realpart(c[i]) = realpart(cc[i+1]) - realpart(cc[i-1]);
892 imagpart(c[i]) = imagpart(cc[i+1]) - imagpart(cc[i-1]);
893 }
894
895 return (void *) c;
896 }
897 } /* end of function cx_d */
898
899
900 void *
cx_floor(void * data,short int type,int length,int * newlength,short int * newtype)901 cx_floor(void *data, short int type, int length, int *newlength, short int *newtype)
902 {
903 *newlength = length;
904 if (type == VF_COMPLEX) {
905 ngcomplex_t *c;
906 ngcomplex_t *cc = (ngcomplex_t *) data;
907 int i;
908
909 c = alloc_c(length);
910 *newtype = VF_COMPLEX;
911 for (i = 0; i < length; i++) {
912 realpart(c[i]) = floor(realpart(cc[i]));
913 imagpart(c[i]) = floor(imagpart(cc[i]));
914 }
915 return ((void *) c);
916 } else {
917 double *d;
918 double *dd = (double *) data;
919 int i;
920
921 d = alloc_d(length);
922 *newtype = VF_REAL;
923 for (i = 0; i < length; i++)
924 d[i] = floor(dd[i]);
925 return ((void *) d);
926 }
927 }
928
929
930 void *
cx_ceil(void * data,short int type,int length,int * newlength,short int * newtype)931 cx_ceil(void *data, short int type, int length, int *newlength, short int *newtype)
932 {
933 *newlength = length;
934 if (type == VF_COMPLEX) {
935 ngcomplex_t *c;
936 ngcomplex_t *cc = (ngcomplex_t *) data;
937 int i;
938
939 c = alloc_c(length);
940 *newtype = VF_COMPLEX;
941 for (i = 0; i < length; i++) {
942 realpart(c[i]) = ceil(realpart(cc[i]));
943 imagpart(c[i]) = ceil(imagpart(cc[i]));
944 }
945 return ((void *) c);
946 } else {
947 double *d;
948 double *dd = (double *) data;
949 int i;
950
951 d = alloc_d(length);
952 *newtype = VF_REAL;
953 for (i = 0; i < length; i++)
954 d[i] = ceil(dd[i]);
955 return ((void *) d);
956 }
957 }
958
959
960 void *
cx_nint(void * data,short int type,int length,int * newlength,short int * newtype)961 cx_nint(void *data, short int type, int length, int *newlength, short int *newtype)
962 {
963 *newlength = length;
964 if (type == VF_COMPLEX) {
965 ngcomplex_t *c;
966 ngcomplex_t *cc = (ngcomplex_t *) data;
967 int i;
968
969 c = alloc_c(length);
970 *newtype = VF_COMPLEX;
971 for (i = 0; i < length; i++) {
972 realpart(c[i]) = nearbyint(realpart(cc[i]));
973 imagpart(c[i]) = nearbyint(imagpart(cc[i]));
974 }
975 return ((void *) c);
976 } else {
977 double *d;
978 double *dd = (double *) data;
979 int i;
980
981 d = alloc_d(length);
982 *newtype = VF_REAL;
983 for (i = 0; i < length; i++)
984 d[i] = nearbyint(dd[i]);
985 return ((void *) d);
986 }
987 }
988