Mercurial > projects > ldc
comparison druntime/src/compiler/dmd/cmath2.d @ 759:d3eb054172f9
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 } |