GRASS GIS 8 Programmer's Manual  8.5.0dev(2024)-bea8435a9e
chouse.c
Go to the documentation of this file.
1 /* chouse.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 
12 void chouse(Cpx *a, double *d, double *dp, int n)
13 {
14  double sc, x, y;
15 
16  Cpx cc, u, *q0;
17 
18  int i, j, k, m, e;
19 
20  Cpx *qw, *pc, *p;
21 
22  q0 = (Cpx *)calloc(2 * n, sizeof(Cpx));
23  for (i = 0, p = q0 + n, pc = a; i < n; ++i, pc += n + 1)
24  *p++ = *pc;
25  for (j = 0, pc = a; j < n - 2; ++j, pc += n + 1) {
26  m = n - j - 1;
27  for (i = 1, sc = 0.; i <= m; ++i)
28  sc += pc[i].re * pc[i].re + pc[i].im * pc[i].im;
29  if (sc > 0.) {
30  sc = sqrt(sc);
31  p = pc + 1;
32  y = sc + (x = sqrt(p->re * p->re + p->im * p->im));
33  if (x > 0.) {
34  cc.re = p->re / x;
35  cc.im = p->im / x;
36  }
37  else {
38  cc.re = 1.;
39  cc.im = 0.;
40  }
41  x = 1. / sqrt(2. * sc * y);
42  y *= x;
43  for (i = 0, qw = pc + 1; i < m; ++i) {
44  q0[i].re = q0[i].im = 0.;
45  if (i) {
46  qw[i].re *= x;
47  qw[i].im *= -x;
48  }
49  else {
50  qw[0].re = y * cc.re;
51  qw[0].im = -y * cc.im;
52  }
53  }
54  for (i = 0, e = j + 2, p = pc + n + 1, y = 0.; i < m;
55  ++i, p += e++) {
56  q0[i].re +=
57  (u.re = qw[i].re) * p->re - (u.im = qw[i].im) * p->im;
58  q0[i].im += u.re * p->im + u.im * p->re;
59  ++p;
60  for (k = i + 1; k < m; ++k, ++p) {
61  q0[i].re += qw[k].re * p->re - qw[k].im * p->im;
62  q0[i].im += qw[k].im * p->re + qw[k].re * p->im;
63  q0[k].re += u.re * p->re + u.im * p->im;
64  q0[k].im += u.im * p->re - u.re * p->im;
65  }
66  y += u.re * q0[i].re + u.im * q0[i].im;
67  }
68  for (i = 0; i < m; ++i) {
69  q0[i].re -= y * qw[i].re;
70  q0[i].re += q0[i].re;
71  q0[i].im -= y * qw[i].im;
72  q0[i].im += q0[i].im;
73  }
74  for (i = 0, e = j + 2, p = pc + n + 1; i < m; ++i, p += e++) {
75  for (k = i; k < m; ++k, ++p) {
76  p->re -= qw[i].re * q0[k].re + qw[i].im * q0[k].im +
77  q0[i].re * qw[k].re + q0[i].im * qw[k].im;
78  p->im -= qw[i].im * q0[k].re - qw[i].re * q0[k].im +
79  q0[i].im * qw[k].re - q0[i].re * qw[k].im;
80  }
81  }
82  }
83  d[j] = pc->re;
84  dp[j] = sc;
85  }
86  d[j] = pc->re;
87  d[j + 1] = (pc + n + 1)->re;
88  u = *(pc + 1);
89  dp[j] = sqrt(u.re * u.re + u.im * u.im);
90  for (j = 0, pc = a, qw = q0 + n; j < n; ++j, pc += n + 1) {
91  *pc = qw[j];
92  for (i = 1, p = pc + n; i < n - j; ++i, p += n) {
93  pc[i].re = p->re;
94  pc[i].im = -p->im;
95  }
96  }
97  free(q0);
98 }
void chouse(Cpx *a, double *d, double *dp, int n)
Definition: chouse.c:12
void free(void *)
Definition: ccmath.h:38
double re
Definition: ccmath.h:39
double im
Definition: ccmath.h:39
#define x