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