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