20 qs = (
Cpx *) calloc(2 * n,
sizeof(
Cpx));
22 for (j = 0, pc = a; j < n - 2; ++j, pc += n + 1, ++q) {
24 for (i = 1, sc = 0.; i <= m; ++i)
25 sc += pc[i].re * pc[i].re + pc[i].im * pc[i].im;
29 y = sc + (x = sqrt(p->
re * p->
re + p->
im * p->
im));
40 x = 1. / sqrt(2. * sc * y);
42 for (i = 0, qw = pc + 1; i < m; ++i) {
43 qs[i].
re = qs[i].
im = 0.;
50 qw[0].
im = -y * cc.
im;
53 for (i = 0, e = j + 2, p = pc + n + 1, y = 0.; i < m;
55 qs[i].
re += (u.
re = qw[i].
re) * p->re - (u.
im =
57 qs[i].
im += u.
re * p->im + u.
im * p->re;
59 for (k = i + 1; k < m; ++k, ++p) {
60 qs[i].
re += qw[k].
re * p->re - qw[k].
im * p->im;
61 qs[i].
im += qw[k].
im * p->re + qw[k].
re * p->im;
62 qs[k].
re += u.
re * p->re + u.
im * p->im;
63 qs[k].
im += u.
im * p->re - u.
re * p->im;
65 y += u.
re * qs[i].
re + u.
im * qs[i].
im;
67 for (i = 0; i < m; ++i) {
68 qs[i].
re -= y * qw[i].
re;
70 qs[i].
im -= y * qw[i].
im;
73 for (i = 0, e = j + 2, p = pc + n + 1; i < m; ++i, p += e++) {
74 for (k = i; k < m; ++k, ++p) {
75 p->re -= qw[i].
re * qs[k].
re + qw[i].
im * qs[k].
im 76 + qs[i].
re * qw[k].
re + qs[i].
im * qw[k].
im;
77 p->im -= qw[i].
im * qs[k].
re - qw[i].
re * qs[k].
im 78 + qs[i].
im * qw[k].
re - qs[i].
re * qw[k].
im;
87 d[j + 1] = (pc += n + 1)->re;
88 dp[j] = sc = sqrt(cc.
re * cc.
re + cc.
im * cc.
im);
91 for (i = 0, m = n + n, p = pc; i < m; ++i, --p)
94 (pc -= n + 1)->re = 1.;
96 for (m = 2; m < n; ++m, qw -= n + 1) {
97 for (j = 0, p = pc, pc->
re = 1.; j < m; ++j, p += n) {
98 for (i = 0, q = p, u.
re = u.
im = 0.; i < m; ++i, ++q) {
99 u.
re += qw[i].re * q->
re - qw[i].im * q->
im;
100 u.
im += qw[i].re * q->
im + qw[i].im * q->
re;
102 for (i = 0, q = p, u.
re += u.
re, u.
im += u.
im; i < m; ++i, ++q) {
103 q->
re -= u.
re * qw[i].re + u.
im * qw[i].im;
104 q->
im -= u.
im * qw[i].re - u.
re * qw[i].im;
107 for (i = 0, p = qw + m - 1; i < n; ++i, --p)
109 (pc -= n + 1)->re = 1.;
111 for (j = 1, p = a + n + 1, q = qs + n, u.
re = 1., u.
im = 0.; j < n;
116 for (i = 1; i < n; ++i, ++p) {
void chousv(Cpx *a, double *d, double *dp, int n)