GRASS GIS 8 Programmer's Manual  8.2.2dev(2023)-3d2c704037
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 #include <stdlib.h>
9 #include "ccmath.h"
10 void chousv(Cpx * a, double *d, double *dp, int n)
11 {
12  double sc, x, y;
13 
14  Cpx cc, u, *qs;
15 
16  int i, j, k, m, e;
17 
18  Cpx *qw, *pc, *p, *q;
19 
20  qs = (Cpx *) calloc(2 * n, sizeof(Cpx));
21  q = qs + n;
22  for (j = 0, pc = a; j < n - 2; ++j, pc += n + 1, ++q) {
23  m = n - j - 1;
24  for (i = 1, sc = 0.; i <= m; ++i)
25  sc += pc[i].re * pc[i].re + pc[i].im * pc[i].im;
26  if (sc > 0.) {
27  sc = sqrt(sc);
28  p = pc + 1;
29  y = sc + (x = sqrt(p->re * p->re + p->im * p->im));
30  if (x > 0.) {
31  cc.re = p->re / x;
32  cc.im = p->im / x;
33  }
34  else {
35  cc.re = 1.;
36  cc.im = 0.;
37  }
38  q->re = -cc.re;
39  q->im = -cc.im;
40  x = 1. / sqrt(2. * sc * y);
41  y *= x;
42  for (i = 0, qw = pc + 1; i < m; ++i) {
43  qs[i].re = qs[i].im = 0.;
44  if (i) {
45  qw[i].re *= x;
46  qw[i].im *= -x;
47  }
48  else {
49  qw[0].re = y * cc.re;
50  qw[0].im = -y * cc.im;
51  }
52  }
53  for (i = 0, e = j + 2, p = pc + n + 1, y = 0.; i < m;
54  ++i, p += e++) {
55  qs[i].re += (u.re = qw[i].re) * p->re - (u.im =
56  qw[i].im) * p->im;
57  qs[i].im += u.re * p->im + u.im * p->re;
58  ++p;
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;
64  }
65  y += u.re * qs[i].re + u.im * qs[i].im;
66  }
67  for (i = 0; i < m; ++i) {
68  qs[i].re -= y * qw[i].re;
69  qs[i].re += qs[i].re;
70  qs[i].im -= y * qw[i].im;
71  qs[i].im += qs[i].im;
72  }
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;
79  }
80  }
81  }
82  d[j] = pc->re;
83  dp[j] = sc;
84  }
85  d[j] = pc->re;
86  cc = *(pc + 1);
87  d[j + 1] = (pc += n + 1)->re;
88  dp[j] = sc = sqrt(cc.re * cc.re + cc.im * cc.im);
89  q->re = cc.re /= sc;
90  q->im = cc.im /= sc;
91  for (i = 0, m = n + n, p = pc; i < m; ++i, --p)
92  p->re = p->im = 0.;
93  pc->re = 1.;
94  (pc -= n + 1)->re = 1.;
95  qw = pc - n;
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;
101  }
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;
105  }
106  }
107  for (i = 0, p = qw + m - 1; i < n; ++i, --p)
108  p->re = p->im = 0.;
109  (pc -= n + 1)->re = 1.;
110  }
111  for (j = 1, p = a + n + 1, q = qs + n, u.re = 1., u.im = 0.; j < n;
112  ++j, ++p, ++q) {
113  sc = u.re * q->re - u.im * q->im;
114  u.im = u.im * q->re + u.re * q->im;
115  u.re = sc;
116  for (i = 1; i < n; ++i, ++p) {
117  sc = u.re * p->re - u.im * p->im;
118  p->im = u.re * p->im + u.im * p->re;
119  p->re = sc;
120  }
121  }
122  free(qs);
123 }
void free(void *)
#define x
void chousv(Cpx *a, double *d, double *dp, int n)
Definition: chousv.c:10
double re
Definition: ccmath.h:38
double im
Definition: ccmath.h:38
Definition: la.h:54