1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
|
/* cpmul.c
*
* Multiply two polynomials with complex coefficients
*
*
*
* SYNOPSIS:
*
* typedef struct
* {
* double r;
* double i;
* }cmplx;
*
* cmplx a[], b[], c[];
* int da, db, dc;
*
* cpmul( a, da, b, db, c, &dc );
*
*
*
* DESCRIPTION:
*
* The two argument polynomials are multiplied together, and
* their product is placed in c.
*
* Each polynomial is represented by its coefficients stored
* as an array of complex number structures (see the typedef).
* The degree of a is da, which must be passed to the routine
* as an argument; similarly the degree db of b is an argument.
* Array a has da + 1 elements and array b has db + 1 elements.
* Array c must have storage allocated for at least da + db + 1
* elements. The value da + db is returned in dc; this is
* the degree of the product polynomial.
*
* Polynomial coefficients are stored in ascending order; i.e.,
* a(x) = a[0]*x**0 + a[1]*x**1 + ... + a[da]*x**da.
*
*
* If desired, c may be the same as either a or b, in which
* case the input argument array is replaced by the product
* array (but only up to terms of degree da + db).
*
*/
/* cpmul */
typedef struct {
double r;
double i;
} cmplx;
void cpmul(cmplx *, int, cmplx *, int, cmplx *, int *);
void cpmul(a, da, b, db, c, dc)
cmplx *a, *b, *c;
int da, db;
int *dc;
{
int i, j, k;
cmplx y;
register cmplx *pa, *pb, *pc;
if (da > db) { /* Know which polynomial has higher degree */
i = da; /* Swapping is OK because args are on the stack */
da = db;
db = i;
pa = a;
a = b;
b = pa;
}
k = da + db;
*dc = k; /* Output the degree of the product */
pc = &c[db + 1];
for (i = db + 1; i <= k; i++) { /* Clear high order terms of output */
pc->r = 0;
pc->i = 0;
pc++;
}
/* To permit replacement of input, work backward from highest degree */
pb = &b[db];
for (j = 0; j <= db; j++) {
pa = &a[da];
pc = &c[k - j];
for (i = 0; i < da; i++) {
y.r = pa->r * pb->r - pa->i * pb->i; /* cmpx multiply */
y.i = pa->r * pb->i + pa->i * pb->r;
pc->r += y.r; /* accumulate partial product */
pc->i += y.i;
pa--;
pc--;
}
y.r = pa->r * pb->r - pa->i * pb->i; /* replace last term, */
y.i = pa->r * pb->i + pa->i * pb->r; /* ...do not accumulate */
pc->r = y.r;
pc->i = y.i;
pb--;
}
}
|
