Mercurial > projects > ldc

comparison druntime/src/compiler/dmd/cmath2.d @ 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 // Copyright (C) 2001-2003 by Digital Mars
3 // All Rights Reserved
4 // www.digitalmars.com
5
6 // Runtime support for complex arithmetic code generation
7 // (for linux)
8
9 /*
10 * Modified by Sean Kelly for use with the D Runtime Project
11 */
12
13 module rt.cmath2;
14
15 private import stdc.math;
16
17 extern (C):
18
19 /****************************
20 * Multiply two complex floating point numbers, x and y.
21 * Input:
22 * x.re ST3
23 * x.im ST2
24 * y.re ST1
25 * y.im ST0
26 * Output:
27 * ST1 real part
28 * ST0 imaginary part
29 */
30
31 void _Cmul()
32 {
33 // p.re = x.re * y.re - x.im * y.im;
34 // p.im = x.im * y.re + x.re * y.im;
35 asm
36 { naked ;
37 fld ST(1) ; // x.re
38 fmul ST,ST(4) ; // ST0 = x.re * y.re
39
40 fld ST(1) ; // y.im
41 fmul ST,ST(4) ; // ST0 = x.im * y.im
42
43 fsubp ST(1),ST ; // ST0 = x.re * y.re - x.im * y.im
44
45 fld ST(3) ; // x.im
46 fmul ST,ST(3) ; // ST0 = x.im * y.re
47
48 fld ST(5) ; // x.re
49 fmul ST,ST(3) ; // ST0 = x.re * y.im
50
51 faddp ST(1),ST ; // ST0 = x.im * y.re + x.re * y.im
52
53 fxch ST(4),ST ;
54 fstp ST(0) ;
55 fxch ST(4),ST ;
56 fstp ST(0) ;
57 fstp ST(0) ;
58 fstp ST(0) ;
59
60 ret ;
61 }
62 /+
63 if (isnan(x) && isnan(y))
64 {
65 // Recover infinities that computed as NaN+ iNaN ...
66 int recalc = 0;
67 if ( isinf( a) || isinf( b) )
68 { // z is infinite
69 // "Box" the infinity and change NaNs in the other factor to 0
70 a = copysignl( isinf( a) ? 1.0 : 0.0, a);
71 b = copysignl( isinf( b) ? 1.0 : 0.0, b);
72 if (isnan( c)) c = copysignl( 0.0, c);
73 if (isnan( d)) d = copysignl( 0.0, d);
74 recalc = 1;
75 }
76 if (isinf(c) || isinf(d))
77 { // w is infinite
78 // "Box" the infinity and change NaNs in the other factor to 0
79 c = copysignl( isinf( c) ? 1.0 : 0.0, c);
80 d = copysignl( isinf( d) ? 1.0 : 0.0, d);
81 if (isnan( a)) a = copysignl( 0.0, a);
82 if (isnan( b)) b = copysignl( 0.0, b);
83 recalc = 1;
84 }
85 if (!recalc && (isinf(ac) || isinf(bd) ||
86 isinf(ad) || isinf(bc)))
87 {
88 // Recover infinities from overflow by changing NaNs to 0 ...
89 if (isnan( a)) a = copysignl( 0.0, a);
90 if (isnan( b)) b = copysignl( 0.0, b);
91 if (isnan( c)) c = copysignl( 0.0, c);
92 if (isnan( d)) d = copysignl( 0.0, d);
93 recalc = 1;
94 }
95 if (recalc)
96 {
97 x = INFINITY * (a * c - b * d);
98 y = INFINITY * (a * d + b * c);
99 }
100 }
101 +/
102 }
103
104 /****************************
105 * Divide two complex floating point numbers, x / y.
106 * Input:
107 * x.re ST3
108 * x.im ST2
109 * y.re ST1
110 * y.im ST0
111 * Output:
112 * ST1 real part
113 * ST0 imaginary part
114 */
115
116 void _Cdiv()
117 {
118 real x_re, x_im;
119 real y_re, y_im;
120 real q_re, q_im;
121 real r;
122 real den;
123
124 asm
125 {
126 fstp y_im ;
127 fstp y_re ;
128 fstp x_im ;
129 fstp x_re ;
130 }
131
132 if (fabs(y_re) < fabs(y_im))
133 {
134 r = y_re / y_im;
135 den = y_im + r * y_re;
136 q_re = (x_re * r + x_im) / den;
137 q_im = (x_im * r - x_re) / den;
138 }
139 else
140 {
141 r = y_im / y_re;
142 den = y_re + r * y_im;
143 q_re = (x_re + r * x_im) / den;
144 q_im = (x_im - r * x_re) / den;
145 }
146 //printf("q.re = %g, q.im = %g\n", (double)q_re, (double)q_im);
147 /+
148 if (isnan(q_re) && isnan(q_im))
149 {
150 real denom = y_re * y_re + y_im * y_im;
151
152 // non-zero / zero
153 if (denom == 0.0 && (!isnan(x_re) || !isnan(x_im)))
154 {
155 q_re = copysignl(INFINITY, y_re) * x_re;
156 q_im = copysignl(INFINITY, y_re) * x_im;
157 }
158 // infinite / finite
159 else if ((isinf(x_re) || isinf(x_im)) && isfinite(y_re) && isfinite(y_im))
160 {
161 x_re = copysignl(isinf(x_re) ? 1.0 : 0.0, x_re);
162 x_im = copysignl(isinf(x_im) ? 1.0 : 0.0, x_im);
163 q_re = INFINITY * (x_re * y_re + x_im * y_im);
164 q_im = INFINITY * (x_im * y_re - x_re * y_im);
165 }
166 // finite / infinite
167 else if (isinf(logbw) && isfinite(x_re) && isfinite(x_im))
168 {
169 y_re = copysignl(isinf(y_re) ? 1.0 : 0.0, y_re);
170 y_im = copysignl(isinf(y_im) ? 1.0 : 0.0, y_im);
171 q_re = 0.0 * (x_re * y_re + x_im * y_im);
172 q_im = 0.0 * (x_im * y_re - x_re * y_im);
173 }
174 }
175 return q_re + q_im * 1.0i;
176 +/
177 asm
178 {
179 fld q_re;
180 fld q_im;
181 }
182 }
183
184 /****************************
185 * Compare two complex floating point numbers, x and y.
186 * Input:
187 * x.re ST3
188 * x.im ST2
189 * y.re ST1
190 * y.im ST0
191 * Output:
192 * 8087 stack is cleared
193 * flags set
194 */
195
196 void _Ccmp()
197 {
198 asm
199 { naked ;
200 fucomp ST(2) ; // compare x.im and y.im
201 fstsw AX ;
202 sahf ;
203 jne L1 ;
204 jp L1 ; // jmp if NAN
205 fucomp ST(2) ; // compare x.re and y.re
206 fstsw AX ;
207 sahf ;
208 fstp ST(0) ; // pop
209 fstp ST(0) ; // pop
210 ret ;
211
212 L1:
213 fstp ST(0) ; // pop
214 fstp ST(0) ; // pop
215 fstp ST(0) ; // pop
216 ret ;
217 }
218 }