Mercurial > projects > ldc

comparison druntime/src/compiler/dmd/complex.c @ 759:d3eb054172f9

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Added copy of druntime from DMD 2.020 modified for LDC.
author Tomas Lindquist Olsen <tomas.l.olsen@gmail.com>
date Tue, 11 Nov 2008 01:52:37 +0100
parents
children
comparison
equal deleted inserted replaced
758:f04dde6e882c 759:d3eb054172f9
1 /*
2 * Placed into the public domain.
3 * Written by Walter Bright
4 * www.digitalmars.com
5 */
6
7
8 #include <math.h>
9
10 typedef struct Complex
11 {
12 long double re;
13 long double im;
14 } Complex;
15
16 Complex _complex_div(Complex x, Complex y)
17 {
18 Complex q;
19 long double r;
20 long double den;
21
22 if (fabs(y.re) < fabs(y.im))
23 {
24 r = y.re / y.im;
25 den = y.im + r * y.re;
26 q.re = (x.re * r + x.im) / den;
27 q.im = (x.im * r - x.re) / den;
28 }
29 else
30 {
31 r = y.im / y.re;
32 den = y.re + r * y.im;
33 q.re = (x.re + r * x.im) / den;
34 q.im = (x.im - r * x.re) / den;
35 }
36 return q;
37 }
38
39 Complex _complex_mul(Complex x, Complex y)
40 {
41 Complex p;
42
43 p.re = x.re * y.re - x.im * y.im;
44 p.im = x.im * y.re + x.re * y.im;
45 return p;
46 }
47
48 long double _complex_abs(Complex z)
49 {
50 long double x,y,ans,temp;
51
52 x = fabs(z.re);
53 y = fabs(z.im);
54 if (x == 0)
55 ans = y;
56 else if (y == 0)
57 ans = x;
58 else if (x > y)
59 {
60 temp = y / x;
61 ans = x * sqrt(1 + temp * temp);
62 }
63 else
64 {
65 temp = x / y;
66 ans = y * sqrt(1 + temp * temp);
67 }
68 return ans;
69 }
70
71 Complex _complex_sqrt(Complex z)
72 {
73 Complex c;
74 long double x,y,w,r;
75
76 if (z.re == 0 && z.im == 0)
77 {
78 c.re = 0;
79 c.im = 0;
80 }
81 else
82 {
83 x = fabs(z.re);
84 y = fabs(z.im);
85 if (x >= y)
86 {
87 r = y / x;
88 w = sqrt(x) * sqrt(0.5 * (1 + sqrt(1 + r * r)));
89 }
90 else
91 {
92 r = x / y;
93 w = sqrt(y) * sqrt(0.5 * (r + sqrt(1 + r * r)));
94 }
95 if (z.re >= 0)
96 {
97 c.re = w;
98 c.im = z.im / (w + w);
99 }
100 else
101 {
102 c.im = (z.im >= 0) ? w : -w;
103 c.re = z.im / (c.im + c.im);
104 }
105 }
106 return c;
107 }