Mercurial > projects > ldc
diff 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 |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/druntime/src/compiler/dmd/cmath2.d Tue Nov 11 01:52:37 2008 +0100 @@ -0,0 +1,218 @@ + +// Copyright (C) 2001-2003 by Digital Mars +// All Rights Reserved +// www.digitalmars.com + +// Runtime support for complex arithmetic code generation +// (for linux) + +/* + * Modified by Sean Kelly for use with the D Runtime Project + */ + +module rt.cmath2; + +private import stdc.math; + +extern (C): + +/**************************** + * Multiply two complex floating point numbers, x and y. + * Input: + * x.re ST3 + * x.im ST2 + * y.re ST1 + * y.im ST0 + * Output: + * ST1 real part + * ST0 imaginary part + */ + +void _Cmul() +{ + // p.re = x.re * y.re - x.im * y.im; + // p.im = x.im * y.re + x.re * y.im; + asm + { naked ; + fld ST(1) ; // x.re + fmul ST,ST(4) ; // ST0 = x.re * y.re + + fld ST(1) ; // y.im + fmul ST,ST(4) ; // ST0 = x.im * y.im + + fsubp ST(1),ST ; // ST0 = x.re * y.re - x.im * y.im + + fld ST(3) ; // x.im + fmul ST,ST(3) ; // ST0 = x.im * y.re + + fld ST(5) ; // x.re + fmul ST,ST(3) ; // ST0 = x.re * y.im + + faddp ST(1),ST ; // ST0 = x.im * y.re + x.re * y.im + + fxch ST(4),ST ; + fstp ST(0) ; + fxch ST(4),ST ; + fstp ST(0) ; + fstp ST(0) ; + fstp ST(0) ; + + ret ; + } +/+ + if (isnan(x) && isnan(y)) + { + // Recover infinities that computed as NaN+ iNaN ... + int recalc = 0; + if ( isinf( a) || isinf( b) ) + { // z is infinite + // "Box" the infinity and change NaNs in the other factor to 0 + a = copysignl( isinf( a) ? 1.0 : 0.0, a); + b = copysignl( isinf( b) ? 1.0 : 0.0, b); + if (isnan( c)) c = copysignl( 0.0, c); + if (isnan( d)) d = copysignl( 0.0, d); + recalc = 1; + } + if (isinf(c) || isinf(d)) + { // w is infinite + // "Box" the infinity and change NaNs in the other factor to 0 + c = copysignl( isinf( c) ? 1.0 : 0.0, c); + d = copysignl( isinf( d) ? 1.0 : 0.0, d); + if (isnan( a)) a = copysignl( 0.0, a); + if (isnan( b)) b = copysignl( 0.0, b); + recalc = 1; + } + if (!recalc && (isinf(ac) || isinf(bd) || + isinf(ad) || isinf(bc))) + { + // Recover infinities from overflow by changing NaNs to 0 ... + if (isnan( a)) a = copysignl( 0.0, a); + if (isnan( b)) b = copysignl( 0.0, b); + if (isnan( c)) c = copysignl( 0.0, c); + if (isnan( d)) d = copysignl( 0.0, d); + recalc = 1; + } + if (recalc) + { + x = INFINITY * (a * c - b * d); + y = INFINITY * (a * d + b * c); + } + } ++/ +} + +/**************************** + * Divide two complex floating point numbers, x / y. + * Input: + * x.re ST3 + * x.im ST2 + * y.re ST1 + * y.im ST0 + * Output: + * ST1 real part + * ST0 imaginary part + */ + +void _Cdiv() +{ + real x_re, x_im; + real y_re, y_im; + real q_re, q_im; + real r; + real den; + + asm + { + fstp y_im ; + fstp y_re ; + fstp x_im ; + fstp x_re ; + } + + if (fabs(y_re) < fabs(y_im)) + { + r = y_re / y_im; + den = y_im + r * y_re; + q_re = (x_re * r + x_im) / den; + q_im = (x_im * r - x_re) / den; + } + else + { + r = y_im / y_re; + den = y_re + r * y_im; + q_re = (x_re + r * x_im) / den; + q_im = (x_im - r * x_re) / den; + } +//printf("q.re = %g, q.im = %g\n", (double)q_re, (double)q_im); +/+ + if (isnan(q_re) && isnan(q_im)) + { + real denom = y_re * y_re + y_im * y_im; + + // non-zero / zero + if (denom == 0.0 && (!isnan(x_re) || !isnan(x_im))) + { + q_re = copysignl(INFINITY, y_re) * x_re; + q_im = copysignl(INFINITY, y_re) * x_im; + } + // infinite / finite + else if ((isinf(x_re) || isinf(x_im)) && isfinite(y_re) && isfinite(y_im)) + { + x_re = copysignl(isinf(x_re) ? 1.0 : 0.0, x_re); + x_im = copysignl(isinf(x_im) ? 1.0 : 0.0, x_im); + q_re = INFINITY * (x_re * y_re + x_im * y_im); + q_im = INFINITY * (x_im * y_re - x_re * y_im); + } + // finite / infinite + else if (isinf(logbw) && isfinite(x_re) && isfinite(x_im)) + { + y_re = copysignl(isinf(y_re) ? 1.0 : 0.0, y_re); + y_im = copysignl(isinf(y_im) ? 1.0 : 0.0, y_im); + q_re = 0.0 * (x_re * y_re + x_im * y_im); + q_im = 0.0 * (x_im * y_re - x_re * y_im); + } + } + return q_re + q_im * 1.0i; ++/ + asm + { + fld q_re; + fld q_im; + } +} + +/**************************** + * Compare two complex floating point numbers, x and y. + * Input: + * x.re ST3 + * x.im ST2 + * y.re ST1 + * y.im ST0 + * Output: + * 8087 stack is cleared + * flags set + */ + +void _Ccmp() +{ + asm + { naked ; + fucomp ST(2) ; // compare x.im and y.im + fstsw AX ; + sahf ; + jne L1 ; + jp L1 ; // jmp if NAN + fucomp ST(2) ; // compare x.re and y.re + fstsw AX ; + sahf ; + fstp ST(0) ; // pop + fstp ST(0) ; // pop + ret ; + + L1: + fstp ST(0) ; // pop + fstp ST(0) ; // pop + fstp ST(0) ; // pop + ret ; + } +}