Mercurial > projects > ldc

comparison druntime/src/compiler/ldc/complex.c @ 1458:e0b2d67cfe7c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Added druntime (this should be removed once it works).
author Robert Clipsham <robert@octarineparrot.com>
date Tue, 02 Jun 2009 17:43:06 +0100
parents
children
comparison
equal deleted inserted replaced
1456:7b218ec1044f 1458:e0b2d67cfe7c
1 /**
2 * Implementation of complex number support routines.
3 *
4 * Copyright: Copyright Digital Mars 2000 - 2009.
5 * License: <a href="http://www.boost.org/LICENSE_1_0.txt>Boost License 1.0</a>.
6 * Authors: Walter Bright
7 *
8 * Copyright Digital Mars 2000 - 2009.
9 * Distributed under the Boost Software License, Version 1.0.
10 * (See accompanying file LICENSE_1_0.txt or copy at
11 * http://www.boost.org/LICENSE_1_0.txt)
12 */
13 #include <math.h>
14
15 typedef struct Complex
16 {
17 long double re;
18 long double im;
19 } Complex;
20
21 Complex _complex_div(Complex x, Complex y)
22 {
23 Complex q;
24 long double r;
25 long double den;
26
27 if (fabs(y.re) < fabs(y.im))
28 {
29 r = y.re / y.im;
30 den = y.im + r * y.re;
31 q.re = (x.re * r + x.im) / den;
32 q.im = (x.im * r - x.re) / den;
33 }
34 else
35 {
36 r = y.im / y.re;
37 den = y.re + r * y.im;
38 q.re = (x.re + r * x.im) / den;
39 q.im = (x.im - r * x.re) / den;
40 }
41 return q;
42 }
43
44 Complex _complex_mul(Complex x, Complex y)
45 {
46 Complex p;
47
48 p.re = x.re * y.re - x.im * y.im;
49 p.im = x.im * y.re + x.re * y.im;
50 return p;
51 }
52
53 long double _complex_abs(Complex z)
54 {
55 long double x,y,ans,temp;
56
57 x = fabs(z.re);
58 y = fabs(z.im);
59 if (x == 0)
60 ans = y;
61 else if (y == 0)
62 ans = x;
63 else if (x > y)
64 {
65 temp = y / x;
66 ans = x * sqrt(1 + temp * temp);
67 }
68 else
69 {
70 temp = x / y;
71 ans = y * sqrt(1 + temp * temp);
72 }
73 return ans;
74 }
75
76 Complex _complex_sqrt(Complex z)
77 {
78 Complex c;
79 long double x,y,w,r;
80
81 if (z.re == 0 && z.im == 0)
82 {
83 c.re = 0;
84 c.im = 0;
85 }
86 else
87 {
88 x = fabs(z.re);
89 y = fabs(z.im);
90 if (x >= y)
91 {
92 r = y / x;
93 w = sqrt(x) * sqrt(0.5 * (1 + sqrt(1 + r * r)));
94 }
95 else
96 {
97 r = x / y;
98 w = sqrt(y) * sqrt(0.5 * (r + sqrt(1 + r * r)));
99 }
100 if (z.re >= 0)
101 {
102 c.re = w;
103 c.im = z.im / (w + w);
104 }
105 else
106 {
107 c.im = (z.im >= 0) ? w : -w;
108 c.re = z.im / (c.im + c.im);
109 }
110 }
111 return c;
112 }