diff --git a/std/math.d b/std/math.d index 948d6c2db4e..1606ebaa452 100644 --- a/std/math.d +++ b/std/math.d @@ -94,6 +94,14 @@ version(LDC) } } +version(LDC) +{ + version(Win64) + { + version = LDC_Win64; + } +} + version(DigitalMars) { version = INLINE_YL2X; // x87 has opcodes for these @@ -431,23 +439,25 @@ T floorImpl(T)(T x) @trusted pure nothrow @nogc public: -// Values obtained from Wolfram Alpha. 116 bits ought to be enough for anybody. -// Wolfram Alpha LLC. 2011. Wolfram|Alpha. http://www.wolframalpha.com/input/?i=e+in+base+16 (access July 6, 2011). -enum real E = 0x1.5bf0a8b1457695355fb8ac404e7a8p+1L; /** e = 2.718281... */ -enum real LOG2T = 0x1.a934f0979a3715fc9257edfe9b5fbp+1L; /** $(SUB log, 2)10 = 3.321928... */ -enum real LOG2E = 0x1.71547652b82fe1777d0ffda0d23a8p+0L; /** $(SUB log, 2)e = 1.442695... */ -enum real LOG2 = 0x1.34413509f79fef311f12b35816f92p-2L; /** $(SUB log, 10)2 = 0.301029... */ -enum real LOG10E = 0x1.bcb7b1526e50e32a6ab7555f5a67cp-2L; /** $(SUB log, 10)e = 0.434294... */ -enum real LN2 = 0x1.62e42fefa39ef35793c7673007e5fp-1L; /** ln 2 = 0.693147... */ -enum real LN10 = 0x1.26bb1bbb5551582dd4adac5705a61p+1L; /** ln 10 = 2.302585... */ -enum real PI = 0x1.921fb54442d18469898cc51701b84p+1L; /** $(_PI) = 3.141592... */ -enum real PI_2 = PI/2; /** $(PI) / 2 = 1.570796... */ -enum real PI_4 = PI/4; /** $(PI) / 4 = 0.785398... */ -enum real M_1_PI = 0x1.45f306dc9c882a53f84eafa3ea69cp-2L; /** 1 / $(PI) = 0.318309... */ -enum real M_2_PI = 2*M_1_PI; /** 2 / $(PI) = 0.636619... */ -enum real M_2_SQRTPI = 0x1.20dd750429b6d11ae3a914fed7fd8p+0L; /** 2 / $(SQRT)$(PI) = 1.128379... */ -enum real SQRT2 = 0x1.6a09e667f3bcc908b2fb1366ea958p+0L; /** $(SQRT)2 = 1.414213... */ -enum real SQRT1_2 = SQRT2/2; /** $(SQRT)$(HALF) = 0.707106... */ +// Values obtained from Wolfram Alpha. +// http://www.wolframalpha.com/input/?i=pi +// http://www.wolframalpha.com/input/?i=e in decimal +// etc +enum real E = 2.718281828459045235360287; /** e = 2.718281... */ +enum real LOG2T = 3.321928094887362347870319; /** $(SUB log, 2)10 = 3.321928... */ +enum real LOG2E = 1.442695040888963407359924; /** $(SUB log, 2)e = 1.442695... */ +enum real LOG2 = 0.301029995663981195213738; /** $(SUB log, 10)2 = 0.301029... */ +enum real LOG10E = 0.434294481903251827651128; /** $(SUB log, 10)e = 0.434294... */ +enum real LN2 = 0.693147180559945309417232; /** ln 2 = 0.693147... */ +enum real LN10 = 2.302585092994045684017991; /** ln 10 = 2.302585... */ +enum real PI = 3.141592653589793238462643; /** $(_PI) = 3.141592... */ +enum real PI_2 = 1.570796326794896619231321; /** $(PI) / 2 = 1.570796... */ +enum real PI_4 = 0.785398163397448309615660; /** $(PI) / 4 = 0.785398... */ +enum real M_1_PI = 0.318309886183790671537767; /** 1 / $(PI) = 0.318309... */ +enum real M_2_PI = 0.636619772367581343075535; /** 2 / $(PI) = 0.636619... */ +enum real M_2_SQRTPI = 1.128379167095512573896158; /** 2 / $(SQRT)$(PI) = 1.128379... */ +enum real SQRT2 = 1.414213562373095048801688; /** $(SQRT)2 = 1.414213... */ +enum real SQRT1_2 = 0.707106781186547524400844; /** $(SQRT)$(HALF) = 0.707106... */ // Note: Make sure the magic numbers in compiler backend for x87 match these. @@ -792,33 +802,33 @@ Lret: {} unittest { static real[2][] vals = // angle,tan - [ - [ 0, 0], - [ .5, .5463024898], - [ 1, 1.557407725], - [ 1.5, 14.10141995], - [ 2, -2.185039863], - [ 2.5,-.7470222972], - [ 3, -.1425465431], - [ 3.5, .3745856402], - [ 4, 1.157821282], - [ 4.5, 4.637332055], - [ 5, -3.380515006], - [ 5.5,-.9955840522], - [ 6, -.2910061914], - [ 6.5, .2202772003], - [ 10, .6483608275], - - // special angles - [ PI_4, 1], - //[ PI_2, real.infinity], // PI_2 is not _exactly_ pi/2. - [ 3*PI_4, -1], - [ PI, 0], - [ 5*PI_4, 1], - //[ 3*PI_2, -real.infinity], - [ 7*PI_4, -1], - [ 2*PI, 0], - ]; + [ + [ 0, 0], + [ .5, .5463024898], + [ 1, 1.557407725], + [ 1.5, 14.10141995], + [ 2, -2.185039863], + [ 2.5,-.7470222972], + [ 3, -.1425465431], + [ 3.5, .3745856402], + [ 4, 1.157821282], + [ 4.5, 4.637332055], + [ 5, -3.380515006], + [ 5.5,-.9955840522], + [ 6, -.2910061914], + [ 6.5, .2202772003], + [ 10, .6483608275], + + // special angles + [ PI_4, 1], + //[ PI_2, real.infinity], // PI_2 is not _exactly_ pi/2. + [ 3*PI_4, -1], + [ PI, 0], + [ 5*PI_4, 1], + //[ 3*PI_2, -real.infinity], + [ 7*PI_4, -1], + [ 2*PI, 0], + ]; int i; for (i = 0; i < vals.length; i++) @@ -1017,11 +1027,18 @@ real atan2(real y, real x) @trusted pure nothrow @nogc { version (Win64) { - asm { + asm + { naked; - fld real ptr [RDX]; // y - fld real ptr [RCX]; // x - fpatan; + sub RSP,8; // space for one 64 bit value (no redzone allowed on win64) + movq [RSP], XMM1; // move xmm1 into memory + fld double ptr [RSP]; // load floating point stack + movq [RSP], XMM0; // move xmm0 into memory + fld double ptr [RSP]; // load floating point stack + fpatan; // calculate the result + fst double ptr [RSP]; // move result into memory + movq XMM0, [RSP]; // store the result in xmm0 + add RSP,8; // restore stack ret; } } @@ -1245,9 +1262,13 @@ public: real acosh(real x) @safe pure nothrow @nogc { if (x > 1/real.epsilon) + { return LN2 + log(x); + } else + { return log(x + sqrt(x*x - 1)); + } } /// ditto @@ -1391,9 +1412,11 @@ extern (C) real rndtonl(real x); */ version(LDC) { - real sqrt(real x) @safe pure nothrow @nogc { return llvm_sqrt(x); } - double sqrt(double x) @safe pure nothrow @nogc { return llvm_sqrt(x); } - float sqrt(float x) @safe pure nothrow @nogc { return llvm_sqrt(x); } + // http://llvm.org/docs/LangRef.html#llvm-sqrt-intrinsic + // sqrt(x) when x is less than zero is undefined + real sqrt(real x) @safe pure nothrow @nogc { return x<0?NAN:llvm_sqrt(x); } + double sqrt(double x) @safe pure nothrow @nogc { return x<0?NAN:llvm_sqrt(x); } + float sqrt(float x) @safe pure nothrow @nogc { return x<0?NAN:llvm_sqrt(x); } } else { @@ -1408,6 +1431,11 @@ real sqrt(real x) @nogc @safe pure nothrow; /* intrinsic */ } +unittest +{ + assert(isNaN(sqrt(-1.0))); +} + unittest { //ctfe @@ -2062,73 +2090,80 @@ unittest } } -unittest +version(LDC_Win64) { - FloatingPointControl ctrl; - if(FloatingPointControl.hasExceptionTraps) - ctrl.disableExceptions(FloatingPointControl.allExceptions); - ctrl.rounding = FloatingPointControl.roundToNearest; - - // @@BUG@@: Non-immutable array literals are ridiculous. - // Note that these are only valid for 80-bit reals: overflow will be different for 64-bit reals. - static const real [2][] exptestpoints = - [ // x, exp(x) - [1.0L, E ], - [0.5L, 0x1.A612_98E1_E069_BC97p+0L ], - [3.0L, E*E*E ], - [0x1.1p13L, 0x1.29aeffefc8ec645p+12557L ], // near overflow - [-0x1.18p13L, 0x1.5e4bf54b4806db9p-12927L ], // near underflow - [-0x1.625p13L, 0x1.a6bd68a39d11f35cp-16358L], - [-0x1p30L, 0 ], // underflow - subnormal - [-0x1.62DAFp13L, 0x1.96c53d30277021dp-16383L ], - [-0x1.643p13L, 0x1p-16444L ], - [-0x1.645p13L, 0 ], // underflow to zero - [0x1p80L, real.infinity ], // far overflow - [real.infinity, real.infinity ], - [0x1.7p13L, real.infinity ] // close overflow - ]; - real x; - IeeeFlags f; - for (int i=0; i