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