/src/rt/complex.c

http://github.com/AlexeyProkhin/druntime
C | 113 lines | 91 code | 10 blank | 12 comment | 14 complexity | d1047f96e46231aacc71b3f7587b5902 MD5 | raw file
```  1/**
2 * Implementation of complex number support routines.
3 *
6 * Authors:   Walter Bright
7 */
8
9/*          Copyright Digital Mars 2000 - 2010.
11 *    (See accompanying file LICENSE or copy at
13 */
14#include <math.h>
15
16typedef struct Complex
17{
18    long double re;
19    long double im;
20} Complex;
21
22Complex _complex_div(Complex x, Complex y)
23{
24    Complex q;
25    long double r;
26    long double den;
27
28    if (fabs(y.re) < fabs(y.im))
29    {
30        r = y.re / y.im;
31        den = y.im + r * y.re;
32        q.re = (x.re * r + x.im) / den;
33        q.im = (x.im * r - x.re) / den;
34    }
35    else
36    {
37        r = y.im / y.re;
38        den = y.re + r * y.im;
39        q.re = (x.re + r * x.im) / den;
40        q.im = (x.im - r * x.re) / den;
41    }
42    return q;
43}
44
45Complex _complex_mul(Complex x, Complex y)
46{
47    Complex p;
48
49    p.re = x.re * y.re - x.im * y.im;
50    p.im = x.im * y.re + x.re * y.im;
51    return p;
52}
53
54long double _complex_abs(Complex z)
55{
56    long double x,y,ans,temp;
57
58    x = fabs(z.re);
59    y = fabs(z.im);
60    if (x == 0)
61        ans = y;
62    else if (y == 0)
63        ans = x;
64    else if (x > y)
65    {
66        temp = y / x;
67        ans = x * sqrt(1 + temp * temp);
68    }
69    else
70    {
71        temp = x / y;
72        ans = y * sqrt(1 + temp * temp);
73    }
74    return ans;
75}
76
77Complex _complex_sqrt(Complex z)
78{
79    Complex c;
80    long double x,y,w,r;
81
82    if (z.re == 0 && z.im == 0)
83    {
84        c.re = 0;
85        c.im = 0;
86    }
87    else
88    {
89        x = fabs(z.re);
90        y = fabs(z.im);
91        if (x >= y)
92        {
93            r = y / x;
94            w = sqrt(x) * sqrt(0.5 * (1 + sqrt(1 + r * r)));
95        }
96        else
97        {
98            r = x / y;
99            w = sqrt(y) * sqrt(0.5 * (r + sqrt(1 + r * r)));
100        }
101        if (z.re >= 0)
102        {
103            c.re = w;
104            c.im = z.im / (w + w);
105        }
106        else
107        {
108            c.im = (z.im >= 0) ? w : -w;
109            c.re = z.im / (c.im + c.im);
110        }
111    }
112    return c;
113}
```