14 Cpx *
ps, *p, *q, *pa, *pd;
18 double s,
t, tq = 0., zr = 1.e-15;
20 q0 = (
Cpx *) calloc(n,
sizeof(
Cpx));
23 for (j = 0; j < n; ++j, ++pa, pd += n + 1) {
25 for (i = 0, p = pa, q = q0; i < n; ++i, p += n)
27 for (i = 1; i < n; ++i) {
30 for (k = 0, p = pa + i * n - j, q = q0; k < lc; ++k, ++q, ++p) {
31 z.
re += p->re * q->re - p->im * q->im;
32 z.
im += p->im * q->re + p->re * q->im;
37 for (i = 0, p = pa, q = q0; i < n; ++i, p += n)
40 s = fabs(pd->re) + fabs(pd->im);
42 for (k = j + 1, ps = pd; k < n; ++k) {
44 if ((t = fabs(ps->
re) + fabs(ps->
im)) > s) {
60 for (k = 0; k < n; ++k) {
66 t = pd->
re * pd->re + pd->im * pd->im;
69 for (k = j + 1, ps = pd + n; k < n; ++k, ps += n) {
70 z.
re = ps->re * h.
re - ps->im * h.
im;
71 z.
im = ps->im * h.
re + ps->re * h.
im;
75 for (j = 1, ps = b + 1; j < n; ++j, ++
ps) {
76 for (k = 0, p = a + n * j, q = b, z.
re = z.
im = 0.; k < j; ++k) {
85 for (j = n - 1, --ps, pd = a + n * n - 1; j >= 0; --j, pd -= n + 1) {
86 for (k = j + 1, p = pd + 1, q = b + j + 1, z.
re = z.
im = 0.; k < n;
95 t = pd->re * pd->re + pd->im * pd->im;
96 ps->
re = (h.
re * pd->re + h.
im * pd->im) / t;
97 ps->
im = (h.
im * pd->re - h.
re * pd->im) / t;
int csolv(Cpx *a, Cpx *b, int n)