← Back to branch summary

~ubuntu-branches/ubuntu/vivid/grass/vivid-proposed

« back to all changes in this revision

Viewing changes to lib/external/ccmath/chousv.c

  • Committer: Package Import Robot
  • Author(s): Bas Couwenberg
  • Date: 2015-02-20 23:12:08 UTC
  • mfrom: (8.2.6 experimental)
  • Revision ID: package-import@ubuntu.com-20150220231208-1u6qvqm84v430b10
Tags: 7.0.0-1~exp1
* New upstream release.
* Update python-ctypes-ternary.patch to use if/else instead of and/or.
* Drop check4dev patch, rely on upstream check.
* Add build dependency on libpq-dev to grass-dev for libpq-fe.h.
* Drop patches applied upstream, refresh remaining patches.
* Update symlinks for images switched from jpg to png.

Show diffs side-by-side

added added

removed removed

Lines of Context:
lib/external/ccmath/chousv.c
 
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
}