GRASS GIS 8 Programmer's Manual 8.3.2(2024)-exported
Loading...
Searching...
No Matches
chousv.c
Go to the documentation of this file.
1/* chousv.c CCMATH mathematics library source code.
2 *
3 * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
4 * This code may be redistributed under the terms of the GNU library
5 * public license (LGPL). ( See the lgpl.license file for details.)
6 * ------------------------------------------------------------------------
7 */
8
9#include <stdlib.h>
10#include "ccmath.h"
11
12void chousv(Cpx *a, double *d, double *dp, int n)
13{
14 double sc, x, y;
15
16 Cpx cc, u, *qs;
17
18 int i, j, k, m, e;
19
20 Cpx *qw, *pc, *p, *q;
21
22 qs = (Cpx *)calloc(2 * n, sizeof(Cpx));
23 q = qs + n;
24 for (j = 0, pc = a; j < n - 2; ++j, pc += n + 1, ++q) {
25 m = n - j - 1;
26 for (i = 1, sc = 0.; i <= m; ++i)
27 sc += pc[i].re * pc[i].re + pc[i].im * pc[i].im;
28 if (sc > 0.) {
29 sc = sqrt(sc);
30 p = pc + 1;
31 y = sc + (x = sqrt(p->re * p->re + p->im * p->im));
32 if (x > 0.) {
33 cc.re = p->re / x;
34 cc.im = p->im / x;
35 }
36 else {
37 cc.re = 1.;
38 cc.im = 0.;
39 }
40 q->re = -cc.re;
41 q->im = -cc.im;
42 x = 1. / sqrt(2. * sc * y);
43 y *= x;
44 for (i = 0, qw = pc + 1; i < m; ++i) {
45 qs[i].re = qs[i].im = 0.;
46 if (i) {
47 qw[i].re *= x;
48 qw[i].im *= -x;
49 }
50 else {
51 qw[0].re = y * cc.re;
52 qw[0].im = -y * cc.im;
53 }
54 }
55 for (i = 0, e = j + 2, p = pc + n + 1, y = 0.; i < m;
56 ++i, p += e++) {
57 qs[i].re +=
58 (u.re = qw[i].re) * p->re - (u.im = qw[i].im) * p->im;
59 qs[i].im += u.re * p->im + u.im * p->re;
60 ++p;
61 for (k = i + 1; k < m; ++k, ++p) {
62 qs[i].re += qw[k].re * p->re - qw[k].im * p->im;
63 qs[i].im += qw[k].im * p->re + qw[k].re * p->im;
64 qs[k].re += u.re * p->re + u.im * p->im;
65 qs[k].im += u.im * p->re - u.re * p->im;
66 }
67 y += u.re * qs[i].re + u.im * qs[i].im;
68 }
69 for (i = 0; i < m; ++i) {
70 qs[i].re -= y * qw[i].re;
71 qs[i].re += qs[i].re;
72 qs[i].im -= y * qw[i].im;
73 qs[i].im += qs[i].im;
74 }
75 for (i = 0, e = j + 2, p = pc + n + 1; i < m; ++i, p += e++) {
76 for (k = i; k < m; ++k, ++p) {
77 p->re -= qw[i].re * qs[k].re + qw[i].im * qs[k].im +
78 qs[i].re * qw[k].re + qs[i].im * qw[k].im;
79 p->im -= qw[i].im * qs[k].re - qw[i].re * qs[k].im +
80 qs[i].im * qw[k].re - qs[i].re * qw[k].im;
81 }
82 }
83 }
84 d[j] = pc->re;
85 dp[j] = sc;
86 }
87 d[j] = pc->re;
88 cc = *(pc + 1);
89 d[j + 1] = (pc += n + 1)->re;
90 dp[j] = sc = sqrt(cc.re * cc.re + cc.im * cc.im);
91 q->re = cc.re /= sc;
92 q->im = cc.im /= sc;
93 for (i = 0, m = n + n, p = pc; i < m; ++i, --p)
94 p->re = p->im = 0.;
95 pc->re = 1.;
96 (pc -= n + 1)->re = 1.;
97 qw = pc - n;
98 for (m = 2; m < n; ++m, qw -= n + 1) {
99 for (j = 0, p = pc, pc->re = 1.; j < m; ++j, p += n) {
100 for (i = 0, q = p, u.re = u.im = 0.; i < m; ++i, ++q) {
101 u.re += qw[i].re * q->re - qw[i].im * q->im;
102 u.im += qw[i].re * q->im + qw[i].im * q->re;
103 }
104 for (i = 0, q = p, u.re += u.re, u.im += u.im; i < m; ++i, ++q) {
105 q->re -= u.re * qw[i].re + u.im * qw[i].im;
106 q->im -= u.im * qw[i].re - u.re * qw[i].im;
107 }
108 }
109 for (i = 0, p = qw + m - 1; i < n; ++i, --p)
110 p->re = p->im = 0.;
111 (pc -= n + 1)->re = 1.;
112 }
113 for (j = 1, p = a + n + 1, q = qs + n, u.re = 1., u.im = 0.; j < n;
114 ++j, ++p, ++q) {
115 sc = u.re * q->re - u.im * q->im;
116 u.im = u.im * q->re + u.re * q->im;
117 u.re = sc;
118 for (i = 1; i < n; ++i, ++p) {
119 sc = u.re * p->re - u.im * p->im;
120 p->im = u.re * p->im + u.im * p->re;
121 p->re = sc;
122 }
123 }
124 free(qs);
125}
void chousv(Cpx *a, double *d, double *dp, int n)
Definition chousv.c:12
double re
Definition ccmath.h:39
double im
Definition ccmath.h:39
#define x