diff --git a/arraycat.d b/arraycat.d index 07b93b914..50b725706 100644 --- a/arraycat.d +++ b/arraycat.d @@ -62,7 +62,7 @@ byte[] _d_arraycatn(uint size, uint n, ...) byte[] _d_arraycopy(uint size, byte[] from, byte[] to) { - //printf("f = %p,%d, t = %p,%d\n", (void*)from, from.length, (void*)to, to.length); + //printf("f = %p,%d, t = %p,%d, size = %d\n", (void*)from, from.length, (void*)to, to.length, size); if (to.length != from.length) { diff --git a/cmath.d b/cmath.d index 15b36533e..68fdbb542 100644 --- a/cmath.d +++ b/cmath.d @@ -1,206 +1,52 @@ - -// Copyright (C) 2001-2003 by Digital Mars +// cmath.d +// Written by Walter Bright +// Copyright (c) 2001-2003 Digital Mars // All Rights Reserved // www.digitalmars.com -// Runtime support for complex arithmetic code generation -import math; +module math; -extern (C): +real fabs(real); +real sqrt(real); -/**************************** - * 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() +creal sqrt(creal z) { - // 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 + creal c; + real x,y,w,r; - 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)) + if (z == 0) { - // 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 - */ - -creal _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; + c = 0; } 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; + { real z_re = z.re; + real z_im = z.im; - // non-zero / zero - if (denom == 0.0 && (!isnan(x_re) || !isnan(x_im))) + x = fabs(z_re); + y = fabs(z_im); + if (x >= y) { - q_re = copysignl(INFINITY, y_re) * x_re; - q_im = copysignl(INFINITY, y_re) * x_im; + r = y / x; + w = sqrt(x) * sqrt(0.5 * (1 + sqrt(1 + r * r))); } - // infinite / finite - else if ((isinf(x_re) || isinf(x_im)) && isfinite(y_re) && isfinite(y_im)) + else { - 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); + r = x / y; + w = sqrt(y) * sqrt(0.5 * (r + sqrt(1 + r * r))); } - // finite / infinite - else if (isinf(logbw) && isfinite(x_re) && isfinite(x_im)) + + if (z_re >= 0) { - 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); + c = w + (z_im / (w + w)) * 1.0i; + } + else + { + if (z_im < 0) + w = -w; + c = z_im / (w + w) + w * 1.0i; } } -+/ - return q_re + q_im * 1.0i; + return c; } -/**************************** - * 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 ; - } -} diff --git a/cmath2.d b/cmath2.d new file mode 100644 index 000000000..15b36533e --- /dev/null +++ b/cmath2.d @@ -0,0 +1,206 @@ + +// Copyright (C) 2001-2003 by Digital Mars +// All Rights Reserved +// www.digitalmars.com + +// Runtime support for complex arithmetic code generation + +import 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 + */ + +creal _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; +} + +/**************************** + * 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 ; + } +} diff --git a/intrinsic.d b/intrinsic.d index 999428fac..9296a3473 100644 --- a/intrinsic.d +++ b/intrinsic.d @@ -8,6 +8,8 @@ /* These functions are built-in intrinsics to the compiler. */ +module math; + int bsf(uint v); int bsr(uint v); int bt(uint *p, uint bitnum); @@ -28,5 +30,6 @@ real fabs(real); real rint(real); long rndtol(real); real sin(real); -real sqrt(real); real ldexp(real, int); +real sqrt(real); + diff --git a/linux.mak b/linux.mak index ee1c0b0cf..f3dd4d520 100644 --- a/linux.mak +++ b/linux.mak @@ -47,12 +47,12 @@ unittest.o : unittest.d OBJS= assert.o deh2.o switch.o complex.o gcstats.o \ critical.o object.o monitor.o arraycat.o invariant.o \ dmain2.o outofmemory.o achar.o aaA.o adi.o file.o \ - compiler.o system.o moduleinit.o \ + compiler.o system.o moduleinit.o cmath.o \ cast.o syserror.o path.o string.o memset.o math.o \ outbuffer.o ctype.o regexp.o random.o linux.o \ stream.o switcherr.o array.o gc.o adi.o \ qsort.o thread.o obj.o \ - crc32.o conv.o arraycast.o errno.o alloca.o cmath.o \ + crc32.o conv.o arraycast.o errno.o alloca.o cmath2.o \ ti_wchar.o ti_uint.o ti_short.o ti_ushort.o \ ti_byte.o ti_ubyte.o ti_long.o ti_ulong.o ti_ptr.o \ ti_float.o ti_double.o ti_real.o ti_delegate.o \ @@ -72,7 +72,7 @@ SRC= mars.h switch.d complex.c critical.c minit.asm \ moduleinit.d cast.d math.d qsort.d \ outbuffer.d unittest.d stream.d ctype.d random.d adi.d \ math2.d thread.d obj.d iunknown.d intrinsic.d time.d memset.d \ - array.d switcherr.d arraycast.d errno.c alloca.d cmath.d \ + array.d switcherr.d arraycast.d errno.c alloca.d cmath.d cmath2.d \ ti_wchar.d ti_uint.d ti_short.d ti_ushort.d \ ti_byte.d ti_ubyte.d ti_long.d ti_ulong.d ti_ptr.d \ ti_float.d ti_double.d ti_real.d ti_delegate.d \ @@ -115,6 +115,9 @@ cast.o : cast.d cmath.o : cmath.d $(DMD) -c $(DFLAGS) cmath.d +cmath2.o : cmath2.d + $(DMD) -c $(DFLAGS) cmath2.d + compiler.o : compiler.d $(DMD) -c $(DFLAGS) compiler.d diff --git a/math.d b/math.d index 33b00e720..c4ce059e1 100644 --- a/math.d +++ b/math.d @@ -7,7 +7,21 @@ module math; import c.stdio; -import intrinsic; + +/* Intrinsics */ + +real cos(real); +real fabs(real); +real rint(real); +long rndtol(real); +real sin(real); +real ldexp(real, int); + +float sqrt(float); +double sqrt(double); +real sqrt(real); +creal sqrt(creal); + extern (C) { diff --git a/math2.d b/math2.d index 9c82646fc..64cfceece 100644 --- a/math2.d +++ b/math2.d @@ -12,7 +12,6 @@ */ module math2; -import intrinsic; private import math, string, c.stdlib, c.stdio; //debug=math2; @@ -48,7 +47,7 @@ long abs(long n) real abs(real n) { // just let the compiler handle it - return intrinsic.fabs(n); + return math.fabs(n); } /********************************* @@ -487,7 +486,7 @@ unittest real acos(real x) { - return atan2(intrinsic.sqrt(1 - x * x), x); + return atan2(math.sqrt(1 - x * x), x); } unittest @@ -502,7 +501,7 @@ unittest real asin(real x) { - return atan2(x, intrinsic.sqrt(1 - x * x)); + return atan2(x, math.sqrt(1 - x * x)); } unittest @@ -528,7 +527,7 @@ real atan(real x) unittest { - assert(feq(atan(intrinsic.sqrt(3)), PI / 3)); + assert(feq(atan(math.sqrt(3.0L)), PI / 3)); } @@ -549,7 +548,7 @@ real atan2(real y, real x) unittest { - assert(feq(atan2(1, intrinsic.sqrt(3)), PI / 6)); + assert(feq(atan2(1, math.sqrt(3.0L)), PI / 6)); } @@ -573,7 +572,7 @@ unittest real asec(real x) { - return intrinsic.cos(1.0 / x); + return math.cos(1.0 / x); } @@ -583,7 +582,7 @@ real asec(real x) real acosec(real x) { - return intrinsic.sin(1.0 / x); + return math.sin(1.0 / x); } /************************************* @@ -604,7 +603,7 @@ real tan(real x) unittest { - assert(feq(tan(PI / 3), intrinsic.sqrt(3))); + assert(feq(tan(PI / 3), math.sqrt(3))); } +/ @@ -625,7 +624,7 @@ real cot(real x) unittest { - assert(feq(cot(PI / 6), intrinsic.sqrt(3))); + assert(feq(cot(PI / 6), math.sqrt(3.0L))); } /************************************* @@ -1057,7 +1056,7 @@ real acosh(real x) else if (x > 1.0e10) return log(2) + log(x); else - return log(x + intrinsic.sqrt((x - 1) * (x + 1))); + return log(x + math.sqrt((x - 1) * (x + 1))); } unittest @@ -1082,8 +1081,8 @@ real asinh(real x) { real z = x * x; return x > 0 ? - log1p(x + z / (1.0 + intrinsic.sqrt(1.0 + z))) : - -log1p(-x + z / (1.0 + intrinsic.sqrt(1.0 + z))); + log1p(x + z / (1.0 + math.sqrt(1.0 + z))) : + -log1p(-x + z / (1.0 + math.sqrt(1.0 + z))); } } diff --git a/unittest.d b/unittest.d index 35e47eddf..ae5ce04a8 100644 --- a/unittest.d +++ b/unittest.d @@ -45,6 +45,12 @@ int main(char[][] args) isValidDchar((dchar)0); // utf uri.ascii2hex(0); // uri +{ + creal c = 3.0 + 4.0i; + c = sqrt(c); + printf("re = %Lg, im = %Lg\n", c.re, c.im); +} + printf("hello world\n"); printf("args.length = %d\n", args.length); for (int i = 0; i < args.length; i++) diff --git a/win32.mak b/win32.mak index ebc358b3f..16c29df07 100644 --- a/win32.mak +++ b/win32.mak @@ -52,7 +52,7 @@ unittest.exe : unittest.d phobos.lib OBJS= assert.obj deh.obj switch.obj complex.obj gcstats.obj \ critical.obj object.obj monitor.obj arraycat.obj invariant.obj \ dmain2.obj outofmemory.obj achar.obj aaA.obj adi.obj file.obj \ - compiler.obj system.obj moduleinit.obj \ + compiler.obj system.obj moduleinit.obj cmath.obj \ cast.obj syserror.obj path.obj string.obj memset.obj math.obj \ outbuffer.obj ctype.obj regexp.obj random.obj windows.obj \ stream.obj switcherr.obj com.obj array.obj gc.obj adi.obj \ @@ -68,7 +68,7 @@ OBJS= assert.obj deh.obj switch.obj complex.obj gcstats.obj \ HDR=mars.h -SRC= switch.d complex.c critical.c errno.c alloca.d cmath.d \ +SRC= switch.d complex.c critical.c errno.c alloca.d cmath2.d \ minit.asm linux.d deh2.d date.d linuxextern.d llmath.d SRC2=deh.c object.d gc.d math.d c\stdio.d c\stdlib.d time.d monitor.c arraycat.d \ @@ -92,7 +92,7 @@ SRC7=ti_wchar.d ti_uint.d ti_short.d ti_ushort.d \ ti_creal.d ti_ireal.d ti_cfloat.d ti_ifloat.d \ ti_cdouble.d ti_idouble.d -SRC8=crc32.d stdint.d conv.d gcstats.d utf.d uri.d +SRC8=crc32.d stdint.d conv.d gcstats.d utf.d uri.d cmath.d phobos.lib : $(OBJS) minit.obj gc2\dmgc.lib win32.mak lib -c phobos.lib $(OBJS) minit.obj gc2\dmgc.lib @@ -103,6 +103,7 @@ adi.obj : adi.d arraycat.obj : arraycat.d assert.obj : assert.d cast.obj : cast.d +cmath.obj : cmath.d compiler.obj : compiler.d complex.obj : mars.h complex.c critical.obj : mars.h critical.c