root/trunk/phobos/std/math.d

Revision 2293, 126.1 kB (checked in by andrei, 3 years ago)

unconditionalized unittest for round

  • Property svn:eol-style set to native
Line 
1 // Written in the D programming language.
2
3 /**
4  * Elementary mathematical functions
5  *
6  * Contains the elementary mathematical functions (powers, roots,
7  * and trignometric functions), and low-level floating-point operations.
8  * Mathematical special functions are available in std.mathspecial.
9  *
10  * The functionality closely follows the IEEE754-2008 standard for
11  * floating-point arithmetic, including the use of camelCase names rather
12  * than C99-style lower case names. All of these functions behave correctly
13  * when presented with an infinity or NaN.
14  *
15  * Unlike C, there is no global 'errno' variable. Consequently, almost all of
16  * these functions are pure nothrow.
17  *
18  * Status:
19  * The gamma and error functions have been superceded by improved versions in
20  * std.mathspecial. They will be officially deprecated in std.math in DMD2.055.
21  * The semantics and names of feqrel and approxEqual will be revised.
22  *
23  * Source: $(PHOBOSSRC std/_math.d)
24  * Macros:
25  *      WIKI = Phobos/StdMath
26  *
27  *      TABLE_SV = <table border=1 cellpadding=4 cellspacing=0>
28  *              <caption>Special Values</caption>
29  *              $0</table>
30  *      SVH = $(TR $(TH $1) $(TH $2))
31  *      SV  = $(TR $(TD $1) $(TD $2))
32  *
33  *      NAN = $(RED NAN)
34  *      SUP = <span style="vertical-align:super;font-size:smaller">$0</span>
35  *      GAMMA = &#915;
36  *      THETA = &theta;
37  *      INTEGRAL = &#8747;
38  *      INTEGRATE = $(BIG &#8747;<sub>$(SMALL $1)</sub><sup>$2</sup>)
39  *      POWER = $1<sup>$2</sup>
40  *      SUB = $1<sub>$2</sub>
41  *      BIGSUM = $(BIG &Sigma; <sup>$2</sup><sub>$(SMALL $1)</sub>)
42  *      CHOOSE = $(BIG &#40;) <sup>$(SMALL $1)</sup><sub>$(SMALL $2)</sub> $(BIG &#41;)
43  *      PLUSMN = &plusmn;
44  *      INFIN = &infin;
45  *      PLUSMNINF = &plusmn;&infin;
46  *      PI = &pi;
47  *      LT = &lt;
48  *      GT = &gt;
49  *      SQRT = &radic;
50  *      HALF = &frac12;
51  *
52  * Copyright: Copyright Digital Mars 2000 - 2009.
53  * License:   <a href="http://www.boost.org/LICENSE_1_0.txt">Boost License 1.0</a>.
54  * Authors:   $(WEB digitalmars.com, Walter Bright),
55  *                        Don Clugston
56  */
57 module std.math;
58
59 import core.stdc.math;
60 import std.range, std.traits;
61
62 version(unittest) {
63     import std.typetuple;
64 }
65
66 version(LDC) {
67     import ldc.intrinsics;
68 }
69
70 version(DigitalMars){
71     version = INLINE_YL2X;        // x87 has opcodes for these
72 }
73
74 version (X86){
75     version = X86_Any;
76 }
77
78 version (X86_64){
79     version = X86_Any;
80 }
81
82 version(D_InlineAsm_X86){
83     version = InlineAsm_X86_Any;
84 }
85 else version(D_InlineAsm_X86_64){
86     version = InlineAsm_X86_Any;
87 }
88
89
90
91
92 private:
93 /*
94  * The following IEEE 'real' formats are currently supported:
95  * 64 bit Big-endian  'double' (eg PowerPC)
96  * 128 bit Big-endian 'quadruple' (eg SPARC)
97  * 64 bit Little-endian 'double' (eg x86-SSE2)
98  * 80 bit Little-endian, with implied bit 'real80' (eg x87, Itanium).
99  * 128 bit Little-endian 'quadruple' (not implemented on any known processor!)
100  *
101  * Non-IEEE 128 bit Big-endian 'doubledouble' (eg PowerPC) has partial support
102  */
103 version(LittleEndian) {
104     static assert(real.mant_dig == 53 || real.mant_dig==64
105                || real.mant_dig == 113,
106       "Only 64-bit, 80-bit, and 128-bit reals"
107       " are supported for LittleEndian CPUs");
108 } else {
109     static assert(real.mant_dig == 53 || real.mant_dig==106
110                || real.mant_dig == 113,
111     "Only 64-bit and 128-bit reals are supported for BigEndian CPUs."
112     " double-double reals have partial support");
113 }
114
115 // Constants used for extracting the components of the representation.
116 // They supplement the built-in floating point properties.
117 template floatTraits(T) {
118     // EXPMASK is a ushort mask to select the exponent portion (without sign)
119     // EXPPOS_SHORT is the index of the exponent when represented as a ushort array.
120     // SIGNPOS_BYTE is the index of the sign when represented as a ubyte array.
121     // RECIP_EPSILON is the value such that (smallest_denormal) * RECIP_EPSILON == T.min_normal
122     enum T RECIP_EPSILON = (1/T.epsilon);
123     static if (T.mant_dig == 24)
124     { // float
125         enum ushort EXPMASK = 0x7F80;
126         enum ushort EXPBIAS = 0x3F00;
127         enum uint EXPMASK_INT = 0x7F80_0000;
128         enum uint MANTISSAMASK_INT = 0x007F_FFFF;
129         version(LittleEndian) {
130             enum EXPPOS_SHORT = 1;
131         } else {
132             enum EXPPOS_SHORT = 0;
133         }
134     }
135     else static if (T.mant_dig == 53) // double, or real==double
136     {
137         enum ushort EXPMASK = 0x7FF0;
138         enum ushort EXPBIAS = 0x3FE0;
139         enum uint EXPMASK_INT = 0x7FF0_0000;
140         enum uint MANTISSAMASK_INT = 0x000F_FFFF; // for the MSB only
141         version(LittleEndian) {
142             enum EXPPOS_SHORT = 3;
143             enum SIGNPOS_BYTE = 7;
144         } else {
145             enum EXPPOS_SHORT = 0;
146             enum SIGNPOS_BYTE = 0;
147         }
148     }
149     else static if (T.mant_dig == 64) // real80
150     {
151         enum ushort EXPMASK = 0x7FFF;
152         enum ushort EXPBIAS = 0x3FFE;
153         version(LittleEndian)
154         {
155             enum EXPPOS_SHORT = 4;
156             enum SIGNPOS_BYTE = 9;
157         }
158         else
159         {
160             enum EXPPOS_SHORT = 0;
161             enum SIGNPOS_BYTE = 0;
162         }
163     } else static if (T.mant_dig == 113){ // quadruple
164         enum ushort EXPMASK = 0x7FFF;
165         version(LittleEndian)
166         {
167             enum EXPPOS_SHORT = 7;
168             enum SIGNPOS_BYTE = 15;
169         }
170         else
171         {
172             enum EXPPOS_SHORT = 0;
173             enum SIGNPOS_BYTE = 0;
174         }
175     } else static if (T.mant_dig == 106) { // doubledouble
176         enum ushort EXPMASK = 0x7FF0;
177         // the exponent byte is not unique
178         version(LittleEndian)
179         {
180             enum EXPPOS_SHORT = 7; // [3] is also an exp short
181             enum SIGNPOS_BYTE = 15;
182         }
183         else
184         {
185             enum EXPPOS_SHORT = 0; // [4] is also an exp short
186             enum SIGNPOS_BYTE = 0;
187         }
188     }
189 }
190
191 // These apply to all floating-point types
192 version(LittleEndian)
193 {
194     enum MANTISSA_LSB = 0;
195     enum MANTISSA_MSB = 1;
196 }
197 else
198 {
199     enum MANTISSA_LSB = 1;
200     enum MANTISSA_MSB = 0;
201 }
202
203 public:
204
205 enum real E =          2.7182818284590452354L;  /** e */ // 0x1.5BF0A8B1_45769535_5FF5p+1L
206 enum real LOG2T =      0x1.a934f0979a3715fcp+1; /** $(SUB log, 2)10 */ // 3.32193 fldl2t
207 enum real LOG2E =      0x1.71547652b82fe178p+0; /** $(SUB log, 2)e */ // 1.4427 fldl2e
208 enum real LOG2 =       0x1.34413509f79fef32p-2; /** $(SUB log, 10)2 */ // 0.30103 fldlg2
209 enum real LOG10E =     0.43429448190325182765;  /** $(SUB log, 10)e */
210 enum real LN2 =        0x1.62e42fefa39ef358p-1; /** ln 2 */  // 0.693147 fldln2
211 enum real LN10 =       2.30258509299404568402;  /** ln 10 */
212 enum real PI =         0x1.921fb54442d1846ap+1; /** $(_PI) */ // 3.14159 fldpi
213 enum real PI_2 =       1.57079632679489661923;  /** $(PI) / 2 */
214 enum real PI_4 =       0.78539816339744830962;  /** $(PI) / 4 */
215 enum real M_1_PI =     0.31830988618379067154;  /** 1 / $(PI) */
216 enum real M_2_PI =     0.63661977236758134308;  /** 2 / $(PI) */
217 enum real M_2_SQRTPI = 1.12837916709551257390;  /** 2 / $(SQRT)$(PI) */
218 enum real SQRT2 =      1.41421356237309504880;  /** $(SQRT)2 */
219 enum real SQRT1_2 =    0.70710678118654752440;  /** $(SQRT)$(HALF) */
220
221 /*
222         Octal versions:
223         PI/64800        0.00001 45530 36176 77347 02143 15351 61441 26767
224         PI/180          0.01073 72152 11224 72344 25603 54276 63351 22056
225         PI/8            0.31103 75524 21026 43021 51423 06305 05600 67016
226         SQRT(1/PI)      0.44067 27240 41233 33210 65616 51051 77327 77303
227         2/PI            0.50574 60333 44710 40522 47741 16537 21752 32335
228         PI/4            0.62207 73250 42055 06043 23046 14612 13401 56034
229         SQRT(2/PI)      0.63041 05147 52066 24106 41762 63612 00272 56161
230
231         PI              3.11037 55242 10264 30215 14230 63050 56006 70163
232         LOG2            0.23210 11520 47674 77674 61076 11263 26013 37111
233  */
234
235 /***********************************
236  * Calculates the absolute value
237  *
238  * For complex numbers, abs(z) = sqrt( $(POWER z.re, 2) + $(POWER z.im, 2) )
239  * = hypot(z.re, z.im).
240  */
241 Num abs(Num)(Num x) @safe pure nothrow
242     if (is(typeof(Num.init >= 0)) && is(typeof(-Num.init)) &&
243             !(is(Num* : const(ifloat*)) || is(Num* : const(idouble*))
244                     || is(Num* : const(ireal*))))
245 {
246     static if (isFloatingPoint!(Num))
247         return fabs(x);
248     else
249         return x>=0 ? x : -x;
250 }
251
252 auto abs(Num)(Num z) @safe pure nothrow
253     if (is(Num* : const(cfloat*)) || is(Num* : const(cdouble*))
254             || is(Num* : const(creal*)))
255 {
256     return hypot(z.re, z.im);
257 }
258
259 /** ditto */
260 real abs(Num)(Num y) @safe pure nothrow
261     if (is(Num* : const(ifloat*)) || is(Num* : const(idouble*))
262             || is(Num* : const(ireal*)))
263 {
264     return fabs(y.im);
265 }
266
267
268 unittest
269 {
270     assert(isIdentical(abs(-0.0L), 0.0L));
271     assert(isNaN(abs(real.nan)));
272     assert(abs(-real.infinity) == real.infinity);
273     assert(abs(-3.2Li) == 3.2L);
274     assert(abs(71.6Li) == 71.6L);
275     assert(abs(-56) == 56);
276     assert(abs(2321312L)  == 2321312L);
277     assert(abs(-1+1i) == sqrt(2.0));
278 }
279
280 /***********************************
281  * Complex conjugate
282  *
283  *  conj(x + iy) = x - iy
284  *
285  * Note that z * conj(z) = $(POWER z.re, 2) - $(POWER z.im, 2)
286  * is always a real number
287  */
288 creal conj(creal z) @safe pure nothrow
289 {
290     return z.re - z.im*1i;
291 }
292
293 /** ditto */
294 ireal conj(ireal y) @safe pure nothrow
295 {
296     return -y;
297 }
298
299
300 unittest
301 {
302     assert(conj(7 + 3i) == 7-3i);
303     ireal z = -3.2Li;
304     assert(conj(z) == -z);
305 }
306
307 /***********************************
308  * Returns cosine of x. x is in radians.
309  *
310  *      $(TABLE_SV
311  *      $(TR $(TH x)                 $(TH cos(x)) $(TH invalid?))
312  *      $(TR $(TD $(NAN))            $(TD $(NAN)) $(TD yes)     )
313  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(NAN)) $(TD yes)     )
314  *      )
315  * Bugs:
316  *      Results are undefined if |x| >= $(POWER 2,64).
317  */
318
319 real cos(real x) @safe pure nothrow;       /* intrinsic */
320
321 /***********************************
322  * Returns sine of x. x is in radians.
323  *
324  *      $(TABLE_SV
325  *      $(TR $(TH x)               $(TH sin(x))      $(TH invalid?))
326  *      $(TR $(TD $(NAN))          $(TD $(NAN))      $(TD yes))
327  *      $(TR $(TD $(PLUSMN)0.0)    $(TD $(PLUSMN)0.0) $(TD no))
328  *      $(TR $(TD $(PLUSMNINF))    $(TD $(NAN))      $(TD yes))
329  *      )
330  * Bugs:
331  *      Results are undefined if |x| >= $(POWER 2,64).
332  */
333
334 real sin(real x) @safe pure nothrow;       /* intrinsic */
335
336
337 /***********************************
338  *  sine, complex and imaginary
339  *
340  *  sin(z) = sin(z.re)*cosh(z.im) + cos(z.re)*sinh(z.im)i
341  *
342  * If both sin($(THETA)) and cos($(THETA)) are required,
343  * it is most efficient to use expi($(THETA)).
344  */
345 creal sin(creal z) @safe pure nothrow
346 {
347     creal cs = expi(z.re);
348     creal csh = coshisinh(z.im);
349     return cs.im * csh.re + cs.re * csh.im * 1i;
350 }
351
352 /** ditto */
353 ireal sin(ireal y) @safe pure nothrow
354 {
355     return cosh(y.im)*1i;
356 }
357
358 unittest
359 {
360   assert(sin(0.0+0.0i) == 0.0);
361   assert(sin(2.0+0.0i) == sin(2.0L) );
362 }
363
364 /***********************************
365  *  cosine, complex and imaginary
366  *
367  *  cos(z) = cos(z.re)*cosh(z.im) - sin(z.re)*sinh(z.im)i
368  */
369 creal cos(creal z) @safe pure nothrow
370 {
371     creal cs = expi(z.re);
372     creal csh = coshisinh(z.im);
373     return cs.re * csh.re - cs.im * csh.im * 1i;
374 }
375
376 /** ditto */
377 real cos(ireal y) @safe pure nothrow
378 {
379     return cosh(y.im);
380 }
381
382 unittest{
383     assert(cos(0.0+0.0i)==1.0);
384     assert(cos(1.3L+0.0i)==cos(1.3L));
385     assert(cos(5.2Li)== cosh(5.2L));
386 }
387
388 /****************************************************************************
389  * Returns tangent of x. x is in radians.
390  *
391  *      $(TABLE_SV
392  *      $(TR $(TH x)             $(TH tan(x))       $(TH invalid?))
393  *      $(TR $(TD $(NAN))        $(TD $(NAN))       $(TD yes))
394  *      $(TR $(TD $(PLUSMN)0.0)  $(TD $(PLUSMN)0.0) $(TD no))
395  *      $(TR $(TD $(PLUSMNINF))  $(TD $(NAN))       $(TD yes))
396  *      )
397  */
398
399 real tan(real x) @trusted pure nothrow
400 {
401     version(D_InlineAsm_X86) {
402         asm
403         {
404             fld     x[EBP]                  ; // load theta
405         }
406     } else version(D_InlineAsm_X86_64) {
407         asm {
408             fld     x[RBP]                  ; // load theta
409         }
410     }
411     version(InlineAsm_X86_Any) {
412     asm
413     {
414         fxam                            ; // test for oddball values
415         fstsw   AX                      ;
416         sahf                            ;
417         jc      trigerr                 ; // x is NAN, infinity, or empty
418                                           // 387's can handle denormals
419 SC18:   fptan                           ;
420         fstp    ST(0)                   ; // dump X, which is always 1
421         fstsw   AX                      ;
422         sahf                            ;
423         jnp     Lret                    ; // C2 = 1 (x is out of range)
424
425         // Do argument reduction to bring x into range
426         fldpi                           ;
427         fxch                            ;
428 SC17:   fprem1                          ;
429         fstsw   AX                      ;
430         sahf                            ;
431         jp      SC17                    ;
432         fstp    ST(1)                   ; // remove pi from stack
433         jmp     SC18                    ;
434
435 trigerr:
436         jnp     Lret                    ; // if theta is NAN, return theta
437         fstp    ST(0)                   ; // dump theta
438     }
439     return real.nan;
440
441 Lret:
442     ;
443     } else {
444         return core.stdc.math.tanl(x);
445     }
446 }
447
448 unittest
449 {
450     static real vals[][2] =     // angle,tan
451         [
452          [   0,   0],
453          [   .5,  .5463024898],
454          [   1,   1.557407725],
455          [   1.5, 14.10141995],
456          [   2,  -2.185039863],
457          [   2.5,-.7470222972],
458          [   3,  -.1425465431],
459          [   3.5, .3745856402],
460          [   4,   1.157821282],
461          [   4.5, 4.637332055],
462          [   5,  -3.380515006],
463          [   5.5,-.9955840522],
464          [   6,  -.2910061914],
465          [   6.5, .2202772003],
466          [   10,  .6483608275],
467
468          // special angles
469          [   PI_4,   1],
470          //[   PI_2,   real.infinity], // PI_2 is not _exactly_ pi/2.
471          [   3*PI_4, -1],
472          [   PI,     0],
473          [   5*PI_4, 1],
474          //[   3*PI_2, -real.infinity],
475          [   7*PI_4, -1],
476          [   2*PI,   0],
477          ];
478     int i;
479
480     for (i = 0; i < vals.length; i++)
481     {
482         real x = vals[i][0];
483         real r = vals[i][1];
484         real t = tan(x);
485
486         //printf("tan(%Lg) = %Lg, should be %Lg\n", x, t, r);
487         if (!isIdentical(r, t)) assert(fabs(r-t) <= .0000001);
488
489         x = -x;
490         r = -r;
491         t = tan(x);
492         //printf("tan(%Lg) = %Lg, should be %Lg\n", x, t, r);
493         if (!isIdentical(r, t) && !(r!<>=0 && t!<>=0)) assert(fabs(r-t) <= .0000001);
494     }
495     // overflow
496     assert(isNaN(tan(real.infinity)));
497     assert(isNaN(tan(-real.infinity)));
498     // NaN propagation
499     assert(isIdentical( tan(NaN(0x0123L)), NaN(0x0123L) ));
500 }
501
502 /***************
503  * Calculates the arc cosine of x,
504  * returning a value ranging from 0 to $(PI).
505  *
506  *      $(TABLE_SV
507  *      $(TR $(TH x)         $(TH acos(x)) $(TH invalid?))
508  *      $(TR $(TD $(GT)1.0)  $(TD $(NAN))  $(TD yes))
509  *      $(TR $(TD $(LT)-1.0) $(TD $(NAN))  $(TD yes))
510  *      $(TR $(TD $(NAN))    $(TD $(NAN))  $(TD yes))
511  *  )
512  */
513 real acos(real x) @safe pure nothrow
514 {
515     return atan2(sqrt(1-x*x), x);
516 }
517
518 /// ditto
519 double acos(double x) @safe pure nothrow { return acos(cast(real)x); }
520 /// ditto
521 float acos(float x) @safe pure nothrow  { return acos(cast(real)x); }
522
523 /***************
524  * Calculates the arc sine of x,
525  * returning a value ranging from -$(PI)/2 to $(PI)/2.
526  *
527  *      $(TABLE_SV
528  *      $(TR $(TH x)            $(TH asin(x))      $(TH invalid?))
529  *      $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no))
530  *      $(TR $(TD $(GT)1.0)     $(TD $(NAN))       $(TD yes))
531  *      $(TR $(TD $(LT)-1.0)    $(TD $(NAN))       $(TD yes))
532  *  )
533  */
534 real asin(real x) @safe pure nothrow
535 {
536     return atan2(x, sqrt(1-x*x));
537 }
538 /// ditto
539 double asin(double x) @safe pure nothrow { return asin(cast(real)x); }
540 /// ditto
541 float asin(float x) @safe pure nothrow  { return asin(cast(real)x); }
542
543 /***************
544  * Calculates the arc tangent of x,
545  * returning a value ranging from -$(PI)/2 to $(PI)/2.
546  *
547  *      $(TABLE_SV
548  *      $(TR $(TH x)                 $(TH atan(x))      $(TH invalid?))
549  *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(PLUSMN)0.0) $(TD no))
550  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(NAN))       $(TD yes))
551  *  )
552  */
553 real atan(real x) @safe pure nothrow { return atan2(x, 1.0L); }
554 /// ditto
555 double atan(double x) @safe pure nothrow { return atan(cast(real)x); }
556 /// ditto
557 float atan(float x)  @safe pure nothrow { return atan(cast(real)x); }
558
559 /***************
560  * Calculates the arc tangent of y / x,
561  * returning a value ranging from -$(PI) to $(PI).
562  *
563  *      $(TABLE_SV
564  *      $(TR $(TH y)                 $(TH x)            $(TH atan(y, x)))
565  *      $(TR $(TD $(NAN))            $(TD anything)     $(TD $(NAN)) )
566  *      $(TR $(TD anything)          $(TD $(NAN))       $(TD $(NAN)) )
567  *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(GT)0.0)     $(TD $(PLUSMN)0.0) )
568  *      $(TR $(TD $(PLUSMN)0.0)      $(TD +0.0)         $(TD $(PLUSMN)0.0) )
569  *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(LT)0.0)     $(TD $(PLUSMN)$(PI)))
570  *      $(TR $(TD $(PLUSMN)0.0)      $(TD -0.0)         $(TD $(PLUSMN)$(PI)))
571  *      $(TR $(TD $(GT)0.0)          $(TD $(PLUSMN)0.0) $(TD $(PI)/2) )
572  *      $(TR $(TD $(LT)0.0)          $(TD $(PLUSMN)0.0) $(TD -$(PI)/2) )
573  *      $(TR $(TD $(GT)0.0)          $(TD $(INFIN))     $(TD $(PLUSMN)0.0) )
574  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD anything)     $(TD $(PLUSMN)$(PI)/2))
575  *      $(TR $(TD $(GT)0.0)          $(TD -$(INFIN))    $(TD $(PLUSMN)$(PI)) )
576  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(INFIN))     $(TD $(PLUSMN)$(PI)/4))
577  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD -$(INFIN))    $(TD $(PLUSMN)3$(PI)/4))
578  *      )
579  */
580 real atan2(real y, real x) @trusted pure nothrow
581 {
582     version(InlineAsm_X86_Any)
583     {
584         asm {
585             fld y;
586             fld x;
587             fpatan;
588         }
589     }
590     else
591     {
592         return core.stdc.math.atan2l(y,x);
593     }
594 }
595
596 /// ditto
597 double atan2(double y, double x) @safe pure nothrow
598 {
599     return atan2(cast(real)y, cast(real)x);
600 }
601
602 /// ditto
603 float atan2(float y, float x) @safe pure nothrow
604 {
605     return atan2(cast(real)y, cast(real)x);
606 }
607
608 /***********************************
609  * Calculates the hyperbolic cosine of x.
610  *
611  *      $(TABLE_SV
612  *      $(TR $(TH x)                 $(TH cosh(x))      $(TH invalid?))
613  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)0.0) $(TD no) )
614  *      )
615  */
616 real cosh(real x) @safe pure nothrow
617 {
618     //  cosh = (exp(x)+exp(-x))/2.
619     // The naive implementation works correctly.
620     real y = exp(x);
621     return (y + 1.0/y) * 0.5;
622 }
623 /// ditto
624 double cosh(double x) @safe pure nothrow { return cosh(cast(real)x); }
625 /// ditto
626 float cosh(float x) @safe pure nothrow  { return cosh(cast(real)x); }
627
628
629 /***********************************
630  * Calculates the hyperbolic sine of x.
631  *
632  *      $(TABLE_SV
633  *      $(TR $(TH x)                 $(TH sinh(x))           $(TH invalid?))
634  *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(PLUSMN)0.0)      $(TD no))
635  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)$(INFIN)) $(TD no))
636  *      )
637  */
638 real sinh(real x) @safe pure nothrow
639 {
640     //  sinh(x) =  (exp(x)-exp(-x))/2;
641     // Very large arguments could cause an overflow, but
642     // the maximum value of x for which exp(x) + exp(-x)) != exp(x)
643     // is x = 0.5 * (real.mant_dig) * LN2. // = 22.1807 for real80.
644     if (fabs(x) > real.mant_dig * LN2) {
645         return copysign(0.5 * exp(fabs(x)), x);
646     }
647     real y = expm1(x);
648     return 0.5 * y / (y+1) * (y+2);
649 }
650 /// ditto
651 double sinh(double x) @safe pure nothrow { return sinh(cast(real)x); }
652 /// ditto
653 float sinh(float x) @safe pure nothrow  { return sinh(cast(real)x); }
654
655
656 /***********************************
657  * Calculates the hyperbolic tangent of x.
658  *
659  *      $(TABLE_SV
660  *      $(TR $(TH x)                 $(TH tanh(x))      $(TH invalid?))
661  *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(PLUSMN)0.0) $(TD no) )
662  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)1.0) $(TD no))
663  *      )
664  */
665 real tanh(real x) @safe pure nothrow
666 {
667     //  tanh(x) = (exp(x) - exp(-x))/(exp(x)+exp(-x))
668     if (fabs(x) > real.mant_dig * LN2) {
669         return copysign(1, x);
670     }
671     real y = expm1(2*x);
672     return y / (y + 2);
673 }
674 /// ditto
675 double tanh(double x) @safe pure nothrow { return tanh(cast(real)x); }
676 /// ditto
677 float tanh(float x) @safe pure nothrow { return tanh(cast(real)x); }
678
679 private:
680 /* Returns cosh(x) + I * sinh(x)
681  * Only one call to exp() is performed.
682  */
683 creal coshisinh(real x) @safe pure nothrow
684 {
685     // See comments for cosh, sinh.
686     if (fabs(x) > real.mant_dig * LN2) {
687         real y = exp(fabs(x));
688         return y * 0.5 + 0.5i * copysign(y, x);
689     } else {
690         real y = expm1(x);
691         return (y + 1.0 + 1.0/(y + 1.0)) * 0.5 + 0.5i * y / (y+1) * (y+2);
692     }
693 }
694
695 unittest {
696     creal c = coshisinh(3.0L);
697     assert(c.re == cosh(3.0L));
698     assert(c.im == sinh(3.0L));
699 }
700
701 public:
702
703 /***********************************
704  * Calculates the inverse hyperbolic cosine of x.
705  *
706  *  Mathematically, acosh(x) = log(x + sqrt( x*x - 1))
707  *
708  * $(TABLE_DOMRG
709  *  $(DOMAIN 1..$(INFIN))
710  *  $(RANGE  1..log(real.max), $(INFIN)) )
711  *      $(TABLE_SV
712  *    $(SVH  x,     acosh(x) )
713  *    $(SV  $(NAN), $(NAN) )
714  *    $(SV  $(LT)1,     $(NAN) )
715  *    $(SV  1,      0       )
716  *    $(SV  +$(INFIN),+$(INFIN))
717  *  )
718  */
719 real acosh(real x) @safe pure nothrow
720 {
721     if (x > 1/real.epsilon)
722         return LN2 + log(x);
723     else
724         return log(x + sqrt(x*x - 1));
725 }
726 /// ditto
727 double acosh(double x) @safe pure nothrow { return acosh(cast(real)x); }
728 /// ditto
729 float acosh(float x) @safe pure nothrow  { return acosh(cast(real)x); }
730
731
732 unittest
733 {
734     assert(isNaN(acosh(0.9)));
735     assert(isNaN(acosh(real.nan)));
736     assert(acosh(1.0)==0.0);
737     assert(acosh(real.infinity) == real.infinity);
738 }
739
740 /***********************************
741  * Calculates the inverse hyperbolic sine of x.
742  *
743  *  Mathematically,
744  *  ---------------
745  *  asinh(x) =  log( x + sqrt( x*x + 1 )) // if x >= +0
746  *  asinh(x) = -log(-x + sqrt( x*x + 1 )) // if x <= -0
747  *  -------------
748  *
749  *    $(TABLE_SV
750  *    $(SVH x,                asinh(x)       )
751  *    $(SV  $(NAN),           $(NAN)         )
752  *    $(SV  $(PLUSMN)0,       $(PLUSMN)0      )
753  *    $(SV  $(PLUSMN)$(INFIN),$(PLUSMN)$(INFIN))
754  *    )
755  */
756 real asinh(real x) @safe pure nothrow
757 {
758     return (fabs(x) > 1 / real.epsilon)
759        // beyond this point, x*x + 1 == x*x
760        ?  copysign(LN2 + log(fabs(x)), x)
761        // sqrt(x*x + 1) ==  1 + x * x / ( 1 + sqrt(x*x + 1) )
762        : copysign(log1p(fabs(x) + x*x / (1 + sqrt(x*x + 1)) ), x);
763 }
764 /// ditto
765 double asinh(double x) @safe pure nothrow { return asinh(cast(real)x); }
766 /// ditto
767 float asinh(float x) @safe pure nothrow { return asinh(cast(real)x); }
768
769 unittest
770 {
771     assert(isIdentical(asinh(0.0), 0.0));
772     assert(isIdentical(asinh(-0.0), -0.0));
773     assert(asinh(real.infinity) == real.infinity);
774     assert(asinh(-real.infinity) == -real.infinity);
775     assert(isNaN(asinh(real.nan)));
776 }
777
778 /***********************************
779  * Calculates the inverse hyperbolic tangent of x,
780  * returning a value from ranging from -1 to 1.
781  *
782  * Mathematically, atanh(x) = log( (1+x)/(1-x) ) / 2
783  *
784  *
785  * $(TABLE_DOMRG
786  *  $(DOMAIN -$(INFIN)..$(INFIN))
787  *  $(RANGE  -1..1) )
788  * $(TABLE_SV
789  *    $(SVH  x,     acosh(x) )
790  *    $(SV  $(NAN), $(NAN) )
791  *    $(SV  $(PLUSMN)0, $(PLUSMN)0)
792  *    $(SV  -$(INFIN), -0)
793  * )
794  */
795 real atanh(real x) @safe pure nothrow
796 {
797     // log( (1+x)/(1-x) ) == log ( 1 + (2*x)/(1-x) )
798     return  0.5 * log1p( 2 * x / (1 - x) );
799 }
800 /// ditto
801 double atanh(double x) @safe pure nothrow { return atanh(cast(real)x); }
802 /// ditto
803 float atanh(float x) @safe pure nothrow { return atanh(cast(real)x); }
804
805
806 unittest
807 {
808     assert(isIdentical(atanh(0.0), 0.0));
809     assert(isIdentical(atanh(-0.0),-0.0));
810     assert(isNaN(atanh(real.nan)));
811     assert(isNaN(atanh(-real.infinity)));
812 }
813
814 /*****************************************
815  * Returns x rounded to a long value using the current rounding mode.
816  * If the integer value of x is
817  * greater than long.max, the result is
818  * indeterminate.
819  */
820 long rndtol(real x) @safe pure nothrow;    /* intrinsic */
821
822
823 /*****************************************
824  * Returns x rounded to a long value using the FE_TONEAREST rounding mode.
825  * If the integer value of x is
826  * greater than long.max, the result is
827  * indeterminate.
828  */
829 extern (C) real rndtonl(real x);
830
831 /***************************************
832  * Compute square root of x.
833  *
834  *      $(TABLE_SV
835  *      $(TR $(TH x)         $(TH sqrt(x))   $(TH invalid?))
836  *      $(TR $(TD -0.0)      $(TD -0.0)      $(TD no))
837  *      $(TR $(TD $(LT)0.0)  $(TD $(NAN))    $(TD yes))
838  *      $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no))
839  *      )
840  */
841
842 @safe pure nothrow
843 {
844     float sqrt(float x);    /* intrinsic */
845     double sqrt(double x);  /* intrinsic */ /// ditto
846     real sqrt(real x);      /* intrinsic */ /// ditto
847 }
848
849 @trusted pure nothrow {  // Should be @safe.  See bugs 4628, 4630.
850     // Create explicit overloads for integer sqrts.  No ddoc for these because
851     // hopefully a more elegant solution will eventually be found, so we don't
852     // want people relying too heavily on the minutiae of this, for example,
853     // by taking the address of sqrt(int) or something.
854     real sqrt(byte x) { return sqrt(cast(real) x); }
855     real sqrt(ubyte x) { return sqrt(cast(real) x); }
856     real sqrt(short x) { return sqrt(cast(real) x); }
857     real sqrt(ushort x) { return sqrt(cast(real) x); }
858     real sqrt(int x) { return sqrt(cast(real) x); }
859     real sqrt(uint x) { return sqrt(cast(real) x); }
860     real sqrt(long x) { return sqrt(cast(real) x); }
861     real sqrt(ulong x) { return sqrt(cast(real) x); }
862 }
863
864 unittest {
865     alias TypeTuple!(byte, ubyte, short, ushort,
866                      int, uint, long, ulong, float, double, real) Numerics;
867     foreach(T; Numerics) {
868         immutable T two = 2;
869         assert(approxEqual(sqrt(two), SQRT2),
870             "sqrt unittest failed on type " ~ T.stringof);
871     }
872 }
873
874 creal sqrt(creal z) @safe pure nothrow
875 {
876     creal c;
877     real x,y,w,r;
878
879     if (z == 0)
880     {
881         c = 0 + 0i;
882     }
883     else
884     {
885         real z_re = z.re;
886         real z_im = z.im;
887
888         x = fabs(z_re);
889         y = fabs(z_im);
890         if (x >= y)
891         {
892             r = y / x;
893             w = sqrt(x) * sqrt(0.5 * (1 + sqrt(1 + r * r)));
894         }
895         else
896         {
897             r = x / y;
898             w = sqrt(y) * sqrt(0.5 * (r + sqrt(1 + r * r)));
899         }
900
901         if (z_re >= 0)
902         {
903             c = w + (z_im / (w + w)) * 1.0i;
904         }
905         else
906         {
907             if (z_im < 0)
908                 w = -w;
909             c = z_im / (w + w) + w * 1.0i;
910         }
911     }
912     return c;
913 }
914
915 /**
916  * Calculates e$(SUP x).
917  *
918  *  $(TABLE_SV
919  *    $(TR $(TH x)             $(TH e$(SUP x)) )
920  *    $(TR $(TD +$(INFIN))     $(TD +$(INFIN)) )
921  *    $(TR $(TD -$(INFIN))     $(TD +0.0)      )
922  *    $(TR $(TD $(NAN))        $(TD $(NAN))    )
923  *  )
924  */
925 real exp(real x) @safe pure nothrow
926 {
927     version(D_InlineAsm_X86)
928     {
929         //  e^^x = 2^^(LOG2E*x)
930         // (This is valid because the overflow & underflow limits for exp
931         // and exp2 are so similar).
932         return exp2(LOG2E*x);
933     }
934     else version(D_InlineAsm_X86_64)
935     {
936         //  e^^x = 2^^(LOG2E*x)
937         // (This is valid because the overflow & underflow limits for exp
938         // and exp2 are so similar).
939         return exp2(LOG2E*x);
940     } else {
941         return core.stdc.math.exp(x);
942     }
943 }
944 /// ditto
945 double exp(double x) @safe pure nothrow  { return exp(cast(real)x); }
946 /// ditto
947 float exp(float x)  @safe pure nothrow   { return exp(cast(real)x); }
948
949
950 /**
951  * Calculates the value of the natural logarithm base (e)
952  * raised to the power of x, minus 1.
953  *
954  * For very small x, expm1(x) is more accurate
955  * than exp(x)-1.
956  *
957  *  $(TABLE_SV
958  *    $(TR $(TH x)             $(TH e$(SUP x)-1)  )
959  *    $(TR $(TD $(PLUSMN)0.0)  $(TD $(PLUSMN)0.0) )
960  *    $(TR $(TD +$(INFIN))     $(TD +$(INFIN))    )
961  *    $(TR $(TD -$(INFIN))     $(TD -1.0)         )
962  *    $(TR $(TD $(NAN))        $(TD $(NAN))       )
963  *  )
964  */
965 real expm1(real x) @trusted pure nothrow
966 {
967     version(D_InlineAsm_X86) {
968         enum { PARAMSIZE = (real.sizeof+3)&(0xFFFF_FFFC) } // always a multiple of 4
969         asm {
970             /*  expm1() for x87 80-bit reals, IEEE754-2008 conformant.
971              * Author: Don Clugston.
972              *
973              *    expm1(x) = 2^^(rndint(y))* 2^^(y-rndint(y)) - 1 where y = LN2*x.
974              *    = 2rndy * 2ym1 + 2rndy - 1, where 2rndy = 2^^(rndint(y))
975              *     and 2ym1 = (2^^(y-rndint(y))-1).
976              *    If 2rndy  < 0.5*real.epsilon, result is -1.
977              *    Implementation is otherwise the same as for exp2()
978              */
979             naked;
980             fld real ptr [ESP+4] ; // x
981             mov AX, [ESP+4+8]; // AX = exponent and sign
982             sub ESP, 12+8; // Create scratch space on the stack
983             // [ESP,ESP+2] = scratchint
984             // [ESP+4..+6, +8..+10, +10] = scratchreal
985             // set scratchreal mantissa = 1.0
986             mov dword ptr [ESP+8], 0;
987             mov dword ptr [ESP+8+4], 0x80000000;
988             and AX, 0x7FFF; // drop sign bit
989             cmp AX, 0x401D; // avoid InvalidException in fist
990             jae L_extreme;
991             fldl2e;
992             fmulp ST(1), ST; // y = x*log2(e)
993             fist dword ptr [ESP]; // scratchint = rndint(y)
994             fisub dword ptr [ESP]; // y - rndint(y)
995             // and now set scratchreal exponent
996             mov EAX, [ESP];
997             add EAX, 0x3fff;
998             jle short L_largenegative;
999             cmp EAX,0x8000;
1000             jge short L_largepositive;
1001             mov [ESP+8+8],AX;
1002             f2xm1; // 2ym1 = 2^^(y-rndint(y)) -1
1003             fld real ptr [ESP+8] ; // 2rndy = 2^^rndint(y)
1004             fmul ST(1), ST;  // ST=2rndy, ST(1)=2rndy*2ym1
1005             fld1;
1006             fsubp ST(1), ST; // ST = 2rndy-1, ST(1) = 2rndy * 2ym1 - 1
1007             faddp ST(1), ST; // ST = 2rndy * 2ym1 + 2rndy - 1
1008             add ESP,12+8;
1009             ret PARAMSIZE;
1010
1011 L_extreme:  // Extreme exponent. X is very large positive, very
1012             // large negative, infinity, or NaN.
1013             fxam;
1014             fstsw AX;
1015             test AX, 0x0400; // NaN_or_zero, but we already know x!=0
1016             jz L_was_nan;  // if x is NaN, returns x
1017             test AX, 0x0200;
1018             jnz L_largenegative;
1019 L_largepositive:
1020             // Set scratchreal = real.max.
1021             // squaring it will create infinity, and set overflow flag.
1022             mov word  ptr [ESP+8+8], 0x7FFE;
1023             fstp ST(0), ST;
1024             fld real ptr [ESP+8];  // load scratchreal
1025             fmul ST(0), ST;        // square it, to create havoc!
1026 L_was_nan:
1027             add ESP,12+8;
1028             ret PARAMSIZE;
1029 L_largenegative:
1030             fstp ST(0), ST;
1031             fld1;
1032             fchs; // return -1. Underflow flag is not set.
1033             add ESP,12+8;
1034             ret PARAMSIZE;
1035         }
1036     } else version(D_InlineAsm_X86_64) {
1037         asm
1038         {
1039             /*  expm1() for x87 80-bit reals, IEEE754-2008 conformant.
1040              * Author: Don Clugston.
1041              *
1042              *    expm1(x) = 2^(rndint(y))* 2^(y-rndint(y)) - 1 where y = LN2*x.
1043              *    = 2rndy * 2ym1 + 2rndy - 1, where 2rndy = 2^(rndint(y))
1044              *     and 2ym1 = (2^(y-rndint(y))-1).
1045              *    If 2rndy  < 0.5*real.epsilon, result is -1.
1046              *    Implementation is otherwise the same as for exp2()
1047              */
1048             naked;
1049             fld real ptr [RSP+8] ; // x
1050             mov AX, [RSP+8+8]; // AX = exponent and sign
1051             sub RSP, 24;       // Create scratch space on the stack
1052             // [RSP,RSP+2] = scratchint
1053             // [RSP+4..+6, +8..+10, +10] = scratchreal
1054             // set scratchreal mantissa = 1.0
1055             mov dword ptr [RSP+8], 0;
1056             mov dword ptr [RSP+8+4], 0x80000000;
1057             and AX, 0x7FFF; // drop sign bit
1058             cmp AX, 0x401D; // avoid InvalidException in fist
1059             jae L_extreme;
1060             fldl2e;
1061             fmul ; // y = x*log2(e)
1062             fist dword ptr [RSP]; // scratchint = rndint(y)
1063             fisub dword ptr [RSP]; // y - rndint(y)
1064             // and now set scratchreal exponent
1065             mov EAX, [RSP];
1066             add EAX, 0x3fff;
1067             jle short L_largenegative;
1068             cmp EAX,0x8000;
1069             jge short L_largepositive;
1070             mov [RSP+8+8],AX;
1071             f2xm1; // 2^(y-rndint(y)) -1
1072             fld real ptr [RSP+8] ; // 2^rndint(y)
1073             fmul ST(1), ST;
1074             fld1;
1075             fsubp ST(1), ST;
1076             fadd;
1077             add RSP,24;
1078             ret;
1079            
1080 L_extreme: // Extreme exponent. X is very large positive, very
1081             // large negative, infinity, or NaN.
1082             fxam;
1083             fstsw AX;
1084             test AX, 0x0400; // NaN_or_zero, but we already know x!=0
1085             jz L_was_nan;  // if x is NaN, returns x
1086             test AX, 0x0200;
1087             jnz L_largenegative;
1088 L_largepositive: 
1089             // Set scratchreal = real.max.
1090             // squaring it will create infinity, and set overflow flag.
1091             mov word  ptr [RSP+8+8], 0x7FFE;
1092             fstp ST(0), ST;
1093             fld real ptr [RSP+8];  // load scratchreal
1094             fmul ST(0), ST;        // square it, to create havoc!
1095 L_was_nan: 
1096             add RSP,24;
1097             ret;
1098            
1099 L_largenegative: 
1100             fstp ST(0), ST;
1101             fld1;
1102             fchs; // return -1. Underflow flag is not set.
1103             add RSP,24;
1104             ret;
1105         }
1106     } else {
1107         return core.stdc.math.expm1(x);
1108     }
1109 }
1110
1111
1112
1113 /**
1114  * Calculates 2$(SUP x).
1115  *
1116  *  $(TABLE_SV
1117  *    $(TR $(TH x)             $(TH exp2(x))   )
1118  *    $(TR $(TD +$(INFIN))     $(TD +$(INFIN)) )
1119  *    $(TR $(TD -$(INFIN))     $(TD +0.0)      )
1120  *    $(TR $(TD $(NAN))        $(TD $(NAN))    )
1121  *  )
1122  */
1123 real exp2(real x) @trusted pure nothrow
1124 {
1125     version(D_InlineAsm_X86) {
1126         enum { PARAMSIZE = (real.sizeof+3)&(0xFFFF_FFFC) } // always a multiple of 4
1127         asm {
1128             /*  exp2() for x87 80-bit reals, IEEE754-2008 conformant.
1129              * Author: Don Clugston.
1130              *
1131              * exp2(x) = 2^^(rndint(x))* 2^^(y-rndint(x))
1132              * The trick for high performance is to avoid the fscale(28cycles on core2),
1133              * frndint(19 cycles), leaving f2xm1(19 cycles) as the only slow instruction.
1134              *
1135              * We can do frndint by using fist. BUT we can't use it for huge numbers,
1136              * because it will set the Invalid Operation flag if overflow or NaN occurs.
1137              * Fortunately, whenever this happens the result would be zero or infinity.
1138              *
1139              * We can perform fscale by directly poking into the exponent. BUT this doesn't
1140              * work for the (very rare) cases where the result is subnormal. So we fall back
1141              * to the slow method in that case.
1142              */
1143             naked;
1144             fld real ptr [ESP+4] ; // x
1145             mov AX, [ESP+4+8]; // AX = exponent and sign
1146             sub ESP, 12+8; // Create scratch space on the stack
1147             // [ESP,ESP+2] = scratchint
1148             // [ESP+4..+6, +8..+10, +10] = scratchreal
1149             // set scratchreal mantissa = 1.0
1150             mov dword ptr [ESP+8], 0;
1151             mov dword ptr [ESP+8+4], 0x80000000;
1152             and AX, 0x7FFF; // drop sign bit
1153             cmp AX, 0x401D; // avoid InvalidException in fist
1154             jae L_extreme;
1155             fist dword ptr [ESP]; // scratchint = rndint(x)
1156             fisub dword ptr [ESP]; // x - rndint(x)
1157             // and now set scratchreal exponent
1158             mov EAX, [ESP];
1159             add EAX, 0x3fff;
1160             jle short L_subnormal;
1161             cmp EAX,0x8000;
1162             jge short L_overflow;
1163             mov [ESP+8+8],AX;
1164 L_normal:
1165             f2xm1;
1166             fld1;
1167             faddp ST(1), ST; // 2^^(x-rndint(x))
1168             fld real ptr [ESP+8] ; // 2^^rndint(x)
1169             add ESP,12+8;
1170             fmulp ST(1), ST;
1171             ret PARAMSIZE;
1172
1173 L_subnormal:
1174             // Result will be subnormal.
1175             // In this rare case, the simple poking method doesn't work.
1176             // The speed doesn't matter, so use the slow fscale method.
1177             fild dword ptr [ESP];  // scratchint
1178             fld1;
1179             fscale;
1180             fstp real ptr [ESP+8]; // scratchreal = 2^^scratchint
1181             fstp ST(0),ST;         // drop scratchint
1182             jmp L_normal;
1183
1184 L_extreme:  // Extreme exponent. X is very large positive, very
1185             // large negative, infinity, or NaN.
1186             fxam;
1187             fstsw AX;
1188             test AX, 0x0400; // NaN_or_zero, but we already know x!=0
1189             jz L_was_nan;  // if x is NaN, returns x
1190             // set scratchreal = real.min_normal
1191             // squaring it will return 0, setting underflow flag
1192             mov word  ptr [ESP+8+8], 1;
1193             test AX, 0x0200;
1194             jnz L_waslargenegative;
1195 L_overflow:
1196             // Set scratchreal = real.max.
1197             // squaring it will create infinity, and set overflow flag.
1198             mov word  ptr [ESP+8+8], 0x7FFE;
1199 L_waslargenegative:
1200             fstp ST(0), ST;
1201             fld real ptr [ESP+8];  // load scratchreal
1202             fmul ST(0), ST;        // square it, to create havoc!
1203 L_was_nan:
1204             add ESP,12+8;
1205             ret PARAMSIZE;
1206         }
1207     } else version(D_InlineAsm_X86_64) {
1208         asm {
1209             /*  exp2() for x87 80-bit reals, IEEE754-2008 conformant.
1210              * Author: Don Clugston.
1211              *
1212              * exp2(x) = 2^(rndint(x))* 2^(y-rndint(x))
1213              * The trick for high performance is to avoid the fscale(28cycles on core2),
1214              * frndint(19 cycles), leaving f2xm1(19 cycles) as the only slow instruction.
1215              *
1216              * We can do frndint by using fist. BUT we can't use it for huge numbers,
1217              * because it will set the Invalid Operation flag is overflow or NaN occurs.
1218              * Fortunately, whenever this happens the result would be zero or infinity.
1219              *
1220              * We can perform fscale by directly poking into the exponent. BUT this doesn't
1221              * work for the (very rare) cases where the result is subnormal. So we fall back
1222              * to the slow method in that case.
1223              */
1224             naked;
1225             fld real ptr [RSP+8] ; // x
1226             mov AX, [RSP+8+8]; // AX = exponent and sign
1227             sub RSP, 24; // Create scratch space on the stack
1228             // [RSP,RSP+2] = scratchint
1229             // [RSP+4..+6, +8..+10, +10] = scratchreal
1230             // set scratchreal mantissa = 1.0
1231             mov dword ptr [RSP+8], 0;
1232             mov dword ptr [RSP+8+4], 0x80000000;
1233             and AX, 0x7FFF; // drop sign bit
1234             cmp AX, 0x401D; // avoid InvalidException in fist
1235             jae L_extreme;
1236             fist dword ptr [RSP]; // scratchint = rndint(x)
1237             fisub dword ptr [RSP]; // x - rndint(x)
1238             // and now set scratchreal exponent
1239             mov EAX, [RSP];
1240             add EAX, 0x3fff;
1241             jle short L_subnormal;
1242             cmp EAX,0x8000;
1243             jge short L_overflow;
1244             mov [RSP+8+8],AX;
1245 L_normal: 
1246             f2xm1;
1247             fld1;
1248             fadd; // 2^(x-rndint(x))
1249             fld real ptr [RSP+8] ; // 2^rndint(x)
1250             add RSP,24;
1251             fmulp ST(1), ST;
1252             ret;
1253      
1254 L_subnormal: 
1255             // Result will be subnormal.
1256             // In this rare case, the simple poking method doesn't work.
1257             // The speed doesn't matter, so use the slow fscale method.
1258             fild dword ptr [RSP];  // scratchint
1259             fld1;
1260             fscale;
1261             fstp real ptr [RSP+8]; // scratchreal = 2^scratchint
1262             fstp ST(0),ST;         // drop scratchint
1263             jmp L_normal;
1264      
1265 L_extreme: // Extreme exponent. X is very large positive, very
1266             // large negative, infinity, or NaN.
1267             fxam;
1268             fstsw AX;
1269             test AX, 0x0400; // NaN_or_zero, but we already know x!=0
1270             jz L_was_nan;  // if x is NaN, returns x
1271             // set scratchreal = real.min
1272             // squaring it will return 0, setting underflow flag
1273             mov word  ptr [RSP+8+8], 1;
1274             test AX, 0x0200;
1275             jnz L_waslargenegative;
1276 L_overflow: 
1277             // Set scratchreal = real.max.
1278             // squaring it will create infinity, and set overflow flag.
1279             mov word  ptr [RSP+8+8], 0x7FFE;
1280 L_waslargenegative: 
1281             fstp ST(0), ST;
1282             fld real ptr [RSP+8];  // load scratchreal
1283             fmul ST(0), ST;        // square it, to create havoc!
1284 L_was_nan: 
1285             add RSP,24;
1286             ret;
1287         }
1288     } else {
1289         return core.stdc.math.exp2(x);
1290     }
1291 }
1292
1293 unittest{
1294     assert(exp2(0.5L)== SQRT2);
1295     assert(exp2(8.0L) == 256.0);
1296     assert(exp2(-9.0L)== 1.0L/512.0);
1297     assert(exp(3.0L) == E*E*E);
1298 }
1299
1300 unittest
1301 {
1302     FloatingPointControl ctrl;
1303     ctrl.disableExceptions(FloatingPointControl.allExceptions);
1304     ctrl.rounding = FloatingPointControl.roundToNearest;
1305
1306     // @@BUG@@: Non-immutable array literals are ridiculous.
1307     // Note that these are only valid for 80-bit reals: overflow will be different for 64-bit reals.
1308     static const real [2][] exptestpoints =
1309     [ // x,            exp(x)
1310         [1.0L,           E                           ],
1311         [0.5L,           0x1.A612_98E1_E069_BC97p+0L ],
1312         [3.0L,           E*E*E                       ],
1313         [0x1.1p13L,      0x1.29aeffefc8ec645p+12557L ], // near overflow
1314         [-0x1.18p13L,    0x1.5e4bf54b4806db9p-12927L ], // near underflow
1315         [-0x1.625p13L,   0x1.a6bd68a39d11f35cp-16358L],
1316         [-0x1p30L,       0                           ], // underflow - subnormal
1317         [-0x1.62DAFp13L, 0x1.96c53d30277021dp-16383L ],
1318         [-0x1.643p13L,   0x1p-16444L                 ],
1319         [-0x1.645p13L,   0                           ], // underflow to zero
1320         [0x1p80L,        real.infinity               ], // far overflow
1321         [real.infinity,  real.infinity               ],
1322         [0x1.7p13L,      real.infinity               ]  // close overflow
1323     ];
1324     real x;
1325     IeeeFlags f;
1326     for (int i=0; i<exptestpoints.length;++i) {
1327         resetIeeeFlags();
1328         x = exp(exptestpoints[i][0]);
1329         f = ieeeFlags();
1330         assert(x == exptestpoints[i][1]);
1331         // Check the overflow bit
1332         assert(f.overflow() == (fabs(x) == real.infinity));
1333         // Check the underflow bit
1334         assert(f.underflow() == (fabs(x) < real.min_normal));
1335         // Invalid and div by zero shouldn't be affected.
1336         assert(!f.invalid);
1337         assert(!f.divByZero);
1338     }
1339     // Ideally, exp(0) would not set the inexact flag.
1340     // Unfortunately, fldl2e sets it!
1341     // So it's not realistic to avoid setting it.
1342     assert(exp(0.0L) == 1.0);
1343
1344     // NaN propagation. Doesn't set flags, bcos was already NaN.
1345     resetIeeeFlags();
1346     x = exp(real.nan);
1347     f = ieeeFlags();
1348     assert(isIdentical(x,real.nan));
1349     assert(f.flags == 0);
1350
1351     resetIeeeFlags();
1352     x = exp(-real.nan);
1353     f = ieeeFlags;
1354     assert(isIdentical(x, -real.nan));
1355     assert(f.flags == 0);
1356
1357     x = exp(NaN(0x123));
1358     assert(isIdentical(x, NaN(0x123)));
1359
1360     // High resolution test
1361     assert(exp(0.5L) == 0x1.A612_98E1_E069_BC97_2DFE_FAB6D_33Fp+0L);
1362
1363 }
1364
1365
1366 /**
1367  * Calculate cos(y) + i sin(y).
1368  *
1369  * On many CPUs (such as x86), this is a very efficient operation;
1370  * almost twice as fast as calculating sin(y) and cos(y) separately,
1371  * and is the preferred method when both are required.
1372  */
1373 creal expi(real y) @trusted pure nothrow
1374 {
1375     version(InlineAsm_X86_Any)
1376     {
1377         asm
1378         {
1379             fld y;
1380             fsincos;
1381             fxch ST(1), ST(0);
1382         }
1383     }
1384     else
1385     {
1386         return cos(y) + sin(y)*1i;
1387     }
1388 }
1389
1390 unittest
1391 {
1392     assert(expi(1.3e5L) == cos(1.3e5L) + sin(1.3e5L) * 1i);
1393     assert(expi(0.0L) == 1L + 0.0Li);
1394 }
1395
1396 /*********************************************************************
1397  * Separate floating point value into significand and exponent.
1398  *
1399  * Returns:
1400  *      Calculate and return $(I x) and $(I exp) such that
1401  *      value =$(I x)*2$(SUP exp) and
1402  *      .5 $(LT)= |$(I x)| $(LT) 1.0
1403  *
1404  *      $(I x) has same sign as value.
1405  *
1406  *      $(TABLE_SV
1407  *      $(TR $(TH value)           $(TH returns)         $(TH exp))
1408  *      $(TR $(TD $(PLUSMN)0.0)    $(TD $(PLUSMN)0.0)    $(TD 0))
1409  *      $(TR $(TD +$(INFIN))       $(TD +$(INFIN))       $(TD int.max))
1410  *      $(TR $(TD -$(INFIN))       $(TD -$(INFIN))       $(TD int.min))
1411  *      $(TR $(TD $(PLUSMN)$(NAN)) $(TD $(PLUSMN)$(NAN)) $(TD int.min))
1412  *      )
1413  */
1414 real frexp(real value, out int exp) @trusted pure nothrow
1415 {
1416     ushort* vu = cast(ushort*)&value;
1417     long* vl = cast(long*)&value;
1418     uint ex;
1419     alias floatTraits!(real) F;
1420
1421     ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
1422     static if (real.mant_dig == 64) { // real80
1423         if (ex) { // If exponent is non-zero
1424             if (ex == F.EXPMASK) {   // infinity or NaN
1425                 if (*vl &  0x7FFF_FFFF_FFFF_FFFF) {  // NaN
1426                     *vl |= 0xC000_0000_0000_0000;  // convert NaNS to NaNQ
1427                     exp = int.min;
1428                 } else if (vu[F.EXPPOS_SHORT] & 0x8000) {   // negative infinity
1429                     exp = int.min;
1430                 } else {   // positive infinity
1431                     exp = int.max;
1432                 }
1433             } else {
1434                 exp = ex - F.EXPBIAS;
1435                 vu[F.EXPPOS_SHORT] = (0x8000 & vu[F.EXPPOS_SHORT]) | 0x3FFE;
1436             }
1437         } else if (!*vl) {
1438             // value is +-0.0
1439             exp = 0;
1440         } else {
1441             // denormal
1442             value *= F.RECIP_EPSILON;
1443             ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
1444             exp = ex - F.EXPBIAS - real.mant_dig + 1;
1445             vu[F.EXPPOS_SHORT] = (0x8000 & vu[F.EXPPOS_SHORT]) | 0x3FFE;
1446         }
1447     } else static if (real.mant_dig == 113) { // quadruple
1448         if (ex) { // If exponent is non-zero
1449             if (ex == F.EXPMASK) {   // infinity or NaN
1450                 if (vl[MANTISSA_LSB] |
1451                     ( vl[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) {  // NaN
1452                     // convert NaNS to NaNQ
1453                     vl[MANTISSA_MSB] |= 0x0000_8000_0000_0000;
1454                     exp = int.min;
1455                 } else if (vu[F.EXPPOS_SHORT] & 0x8000) {   // negative infinity
1456                     exp = int.min;
1457                 } else {   // positive infinity
1458                     exp = int.max;
1459                 }
1460             } else {
1461                 exp = ex - F.EXPBIAS;
1462                 vu[F.EXPPOS_SHORT] =
1463                     cast(ushort)((0x8000 & vu[F.EXPPOS_SHORT]) | 0x3FFE);
1464             }
1465         } else if ((vl[MANTISSA_LSB]
1466                        |(vl[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) == 0) {
1467             // value is +-0.0
1468             exp = 0;
1469         } else {
1470             // denormal
1471             value *= F.RECIP_EPSILON;
1472             ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
1473             exp = ex - F.EXPBIAS - real.mant_dig + 1;
1474             vu[F.EXPPOS_SHORT] =
1475                 cast(ushort)((0x8000 & vu[F.EXPPOS_SHORT]) | 0x3FFE);
1476         }
1477     } else static if (real.mant_dig==53) { // real is double
1478         if (ex) { // If exponent is non-zero
1479             if (ex == F.EXPMASK) {   // infinity or NaN
1480                 if (*vl == 0x7FF0_0000_0000_0000) {  // positive infinity
1481                     exp = int.max;
1482                 } else if (*vl == 0xFFF0_0000_0000_0000) { // negative infinity
1483                     exp = int.min;
1484                 } else { // NaN
1485                     *vl |= 0x0008_0000_0000_0000;  // convert NaNS to NaNQ
1486                     exp = int.min;
1487                 }
1488             } else {
1489                 exp = (ex - F.EXPBIAS) >> 4;
1490                 vu[F.EXPPOS_SHORT] = cast(ushort)((0x8000 & vu[F.EXPPOS_SHORT]) | 0x3FE0);
1491             }
1492         } else if (!(*vl & 0x7FFF_FFFF_FFFF_FFFF)) {
1493             // value is +-0.0
1494             exp = 0;
1495         } else {
1496             // denormal
1497             value *= F.RECIP_EPSILON;
1498             ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
1499             exp = ((ex - F.EXPBIAS)>> 4) - real.mant_dig + 1;
1500             vu[F.EXPPOS_SHORT] =
1501                 cast(ushort)((0x8000 & vu[F.EXPPOS_SHORT]) | 0x3FE0);
1502         }
1503     } else { //static if(real.mant_dig==106) // doubledouble
1504         assert (0, "frexp not implemented");
1505     }
1506     return value;
1507 }
1508
1509
1510 unittest
1511 {
1512     static real vals[][3] =     // x,frexp,exp
1513         [
1514          [0.0,   0.0,    0],
1515          [-0.0,  -0.0,   0],
1516          [1.0,   .5,     1],
1517          [-1.0,  -.5,    1],
1518          [2.0,   .5,     2],
1519          [double.min_normal/2.0, .5, -1022],
1520          [real.infinity,real.infinity,int.max],
1521          [-real.infinity,-real.infinity,int.min],
1522          [real.nan,real.nan,int.min],
1523          [-real.nan,-real.nan,int.min],
1524          ];
1525
1526     int i;
1527
1528     for (i = 0; i < vals.length; i++) {
1529         real x = vals[i][0];
1530         real e = vals[i][1];
1531         int exp = cast(int)vals[i][2];
1532         int eptr;
1533         real v = frexp(x, eptr);
1534         assert(isIdentical(e, v));
1535         assert(exp == eptr);
1536
1537     }
1538     static if (real.mant_dig == 64) {
1539         static real extendedvals[][3] = [ // x,frexp,exp
1540                                           [0x1.a5f1c2eb3fe4efp+73L, 0x1.A5F1C2EB3FE4EFp-1L,   74],    // normal
1541                                           [0x1.fa01712e8f0471ap-1064L,  0x1.fa01712e8f0471ap-1L,     -1063],
1542                                           [real.min_normal,  .5,     -16381],
1543                                           [real.min_normal/2.0L, .5,     -16382]    // denormal
1544                                            ];
1545
1546         for (i = 0; i < extendedvals.length; i++) {
1547             real x = extendedvals[i][0];
1548             real e = extendedvals[i][1];
1549             int exp = cast(int)extendedvals[i][2];
1550             int eptr;
1551             real v = frexp(x, eptr);
1552             assert(isIdentical(e, v));
1553             assert(exp == eptr);
1554
1555         }
1556     }
1557 }
1558
1559 /******************************************
1560  * Extracts the exponent of x as a signed integral value.
1561  *
1562  * If x is not a special value, the result is the same as
1563  * $(D cast(int)logb(x)).
1564  *
1565  *      $(TABLE_SV
1566  *      $(TR $(TH x)                $(TH ilogb(x))     $(TH Range error?))
1567  *      $(TR $(TD 0)                 $(TD FP_ILOGB0)   $(TD yes))
1568  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD int.max)     $(TD no))
1569  *      $(TR $(TD $(NAN))            $(TD FP_ILOGBNAN) $(TD no))
1570  *      )
1571  */
1572 int ilogb(real x)  @trusted nothrow    { return core.stdc.math.ilogbl(x); }
1573
1574 alias core.stdc.math.FP_ILOGB0   FP_ILOGB0;
1575 alias core.stdc.math.FP_ILOGBNAN FP_ILOGBNAN;
1576
1577
1578 /*******************************************
1579  * Compute n * 2$(SUP exp)
1580  * References: frexp
1581  */
1582
1583 real ldexp(real n, int exp) @safe pure nothrow;    /* intrinsic */
1584
1585 unittest {
1586     assert(ldexp(1, -16384) == 0x1p-16384L);
1587     assert(ldexp(1, -16382) == 0x1p-16382L);
1588     int x;
1589     real n = frexp(0x1p-16384L, x);
1590     assert(n==0.5L);
1591     assert(x==-16383);
1592     assert(ldexp(n, x)==0x1p-16384L);
1593
1594 }
1595
1596 /**************************************
1597  * Calculate the natural logarithm of x.
1598  *
1599  *    $(TABLE_SV
1600  *    $(TR $(TH x)            $(TH log(x))    $(TH divide by 0?) $(TH invalid?))
1601  *    $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes)          $(TD no))
1602  *    $(TR $(TD $(LT)0.0)     $(TD $(NAN))    $(TD no)           $(TD yes))
1603  *    $(TR $(TD +$(INFIN))    $(TD +$(INFIN)) $(TD no)           $(TD no))
1604  *    )
1605  */
1606
1607 real log(real x) @safe pure nothrow
1608 {
1609     version (INLINE_YL2X)
1610         return yl2x(x, LN2);
1611     else
1612         return core.stdc.math.logl(x);
1613 }
1614
1615 unittest
1616 {
1617     assert(log(E) == 1);
1618 }
1619
1620 /**************************************
1621  * Calculate the base-10 logarithm of x.
1622  *
1623  *      $(TABLE_SV
1624  *      $(TR $(TH x)            $(TH log10(x))  $(TH divide by 0?) $(TH invalid?))
1625  *      $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes)          $(TD no))
1626  *      $(TR $(TD $(LT)0.0)     $(TD $(NAN))    $(TD no)           $(TD yes))
1627  *      $(TR $(TD +$(INFIN))    $(TD +$(INFIN)) $(TD no)           $(TD no))
1628  *      )
1629  */
1630
1631 real log10(real x) @safe pure nothrow
1632 {
1633     version (INLINE_YL2X)
1634         return yl2x(x, LOG2);
1635     else
1636         return core.stdc.math.log10l(x);
1637 }
1638
1639 unittest
1640 {
1641     //printf("%Lg\n", log10(1000) - 3);
1642     assert(fabs(log10(1000) - 3) < .000001);
1643 }
1644
1645 /******************************************
1646  *      Calculates the natural logarithm of 1 + x.
1647  *
1648  *      For very small x, log1p(x) will be more accurate than
1649  *      log(1 + x).
1650  *
1651  *  $(TABLE_SV
1652  *  $(TR $(TH x)            $(TH log1p(x))     $(TH divide by 0?) $(TH invalid?))
1653  *  $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no)           $(TD no))
1654  *  $(TR $(TD -1.0)         $(TD -$(INFIN))    $(TD yes)          $(TD no))
1655  *  $(TR $(TD $(LT)-1.0)    $(TD $(NAN))       $(TD no)           $(TD yes))
1656  *  $(TR $(TD +$(INFIN))    $(TD -$(INFIN))    $(TD no)           $(TD no))
1657  *  )
1658  */
1659
1660 real log1p(real x) @safe pure nothrow
1661 {
1662     version(INLINE_YL2X)
1663     {
1664         // On x87, yl2xp1 is valid if and only if -0.5 <= lg(x) <= 0.5,
1665         //    ie if -0.29<=x<=0.414
1666         return (fabs(x) <= 0.25)  ? yl2xp1(x, LN2) : yl2x(x+1, LN2);
1667     }
1668     else
1669     {
1670         return core.stdc.math.log1pl(x);
1671     }
1672 }
1673
1674 /***************************************
1675  * Calculates the base-2 logarithm of x:
1676  * $(SUB log, 2)x
1677  *
1678  *  $(TABLE_SV
1679  *  $(TR $(TH x)            $(TH log2(x))   $(TH divide by 0?) $(TH invalid?))
1680  *  $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes)          $(TD no) )
1681  *  $(TR $(TD $(LT)0.0)     $(TD $(NAN))    $(TD no)           $(TD yes) )
1682  *  $(TR $(TD +$(INFIN))    $(TD +$(INFIN)) $(TD no)           $(TD no) )
1683  *  )
1684  */
1685 real log2(real x) @safe pure nothrow
1686 {
1687     version (INLINE_YL2X)
1688         return yl2x(x, 1);
1689     else
1690         return core.stdc.math.log2l(x);
1691 }
1692
1693 /*****************************************
1694  * Extracts the exponent of x as a signed integral value.
1695  *
1696  * If x is subnormal, it is treated as if it were normalized.
1697  * For a positive, finite x:
1698  *
1699  * 1 $(LT)= $(I x) * FLT_RADIX$(SUP -logb(x)) $(LT) FLT_RADIX
1700  *
1701  *      $(TABLE_SV
1702  *      $(TR $(TH x)                 $(TH logb(x))   $(TH divide by 0?) )
1703  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD +$(INFIN)) $(TD no))
1704  *      $(TR $(TD $(PLUSMN)0.0)      $(TD -$(INFIN)) $(TD yes) )
1705  *      )
1706  */
1707 real logb(real x) @trusted nothrow    { return core.stdc.math.logbl(x); }
1708
1709 /************************************
1710  * Calculates the remainder from the calculation x/y.
1711  * Returns:
1712  * The value of x - i * y, where i is the number of times that y can
1713  * be completely subtracted from x. The result has the same sign as x.
1714  *
1715  * $(TABLE_SV
1716  *  $(TR $(TH x)              $(TH y)             $(TH modf(x, y))   $(TH invalid?))
1717  *  $(TR $(TD $(PLUSMN)0.0)   $(TD not 0.0)       $(TD $(PLUSMN)0.0) $(TD no))
1718  *  $(TR $(TD $(PLUSMNINF))   $(TD anything)      $(TD $(NAN))       $(TD yes))
1719  *  $(TR $(TD anything)       $(TD $(PLUSMN)0.0)  $(TD $(NAN))       $(TD yes))
1720  *  $(TR $(TD !=$(PLUSMNINF)) $(TD $(PLUSMNINF))  $(TD x)            $(TD no))
1721  * )
1722  */
1723 real modf(real x, ref real y) @trusted nothrow { return core.stdc.math.modfl(x,&y); }
1724
1725 /*************************************
1726  * Efficiently calculates x * 2$(SUP n).
1727  *
1728  * scalbn handles underflow and overflow in
1729  * the same fashion as the basic arithmetic operators.
1730  *
1731  *      $(TABLE_SV
1732  *      $(TR $(TH x)                 $(TH scalb(x)))
1733  *      $(TR $(TD $(PLUSMNINF))      $(TD $(PLUSMNINF)) )
1734  *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(PLUSMN)0.0) )
1735  *      )
1736  */
1737 real scalbn(real x, int n) @trusted nothrow
1738 {
1739     version(InlineAsm_X86_Any) {
1740         // scalbnl is not supported on DMD-Windows, so use asm.
1741         asm {
1742             fild n;
1743             fld x;
1744             fscale;
1745             fstp ST(1), ST;
1746         }
1747     } else {
1748         return core.stdc.math.scalbnl(x, n);
1749     }
1750 }
1751
1752 unittest {
1753     assert(scalbn(-real.infinity, 5) == -real.infinity);
1754 }
1755
1756 /***************
1757  * Calculates the cube root of x.
1758  *
1759  *      $(TABLE_SV
1760  *      $(TR $(TH $(I x))            $(TH cbrt(x))           $(TH invalid?))
1761  *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(PLUSMN)0.0)      $(TD no) )
1762  *      $(TR $(TD $(NAN))            $(TD $(NAN))            $(TD yes) )
1763  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)$(INFIN)) $(TD no) )
1764  *      )
1765  */
1766 real cbrt(real x) @trusted nothrow    { return core.stdc.math.cbrtl(x); }
1767
1768
1769 /*******************************
1770  * Returns |x|
1771  *
1772  *      $(TABLE_SV
1773  *      $(TR $(TH x)                 $(TH fabs(x)))
1774  *      $(TR $(TD $(PLUSMN)0.0)      $(TD +0.0) )
1775  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD +$(INFIN)) )
1776  *      )
1777  */
1778 real fabs(real x) @safe pure nothrow;      /* intrinsic */
1779
1780
1781 /***********************************************************************
1782  * Calculates the length of the
1783  * hypotenuse of a right-angled triangle with sides of length x and y.
1784  * The hypotenuse is the value of the square root of
1785  * the sums of the squares of x and y:
1786  *
1787  *      sqrt($(POW x, 2) + $(POW y, 2))
1788  *
1789  * Note that hypot(x, y), hypot(y, x) and
1790  * hypot(x, -y) are equivalent.
1791  *
1792  *  $(TABLE_SV
1793  *  $(TR $(TH x)            $(TH y)            $(TH hypot(x, y)) $(TH invalid?))
1794  *  $(TR $(TD x)            $(TD $(PLUSMN)0.0) $(TD |x|)         $(TD no))
1795  *  $(TR $(TD $(PLUSMNINF)) $(TD y)            $(TD +$(INFIN))   $(TD no))
1796  *  $(TR $(TD $(PLUSMNINF)) $(TD $(NAN))       $(TD +$(INFIN))   $(TD no))
1797  *  )
1798  */
1799
1800 real hypot(real x, real y) @safe pure nothrow
1801 {
1802     // Scale x and y to avoid underflow and overflow.
1803     // If one is huge and the other tiny, return the larger.
1804     // If both are huge, avoid overflow by scaling by 1/sqrt(real.max/2).
1805     // If both are tiny, avoid underflow by scaling by sqrt(real.min_normal*real.epsilon).
1806
1807     enum real SQRTMIN = 0.5*sqrt(real.min_normal); // This is a power of 2.
1808     enum real SQRTMAX = 1.0L/SQRTMIN; // 2^^((max_exp)/2) = nextUp(sqrt(real.max))
1809
1810     static assert(2*(SQRTMAX/2)*(SQRTMAX/2) <= real.max);
1811     static assert(real.min_normal*real.max>2 && real.min_normal*real.max<=4); // Proves that sqrt(real.max) ~~  0.5/sqrt(real.min_normal)
1812
1813     real u = fabs(x);
1814     real v = fabs(y);
1815     if (!(u >= v))  // check for NaN as well.
1816     {
1817         v = u;
1818         u = fabs(y);
1819         if (u == real.infinity) return u; // hypot(inf, nan) == inf
1820         if (v == real.infinity) return v; // hypot(nan, inf) == inf
1821     }
1822     // Now u >= v, or else one is NaN.
1823     if (v >= SQRTMAX*0.5)
1824     {
1825             // hypot(huge, huge) -- avoid overflow
1826         u *= SQRTMIN*0.5;
1827         v *= SQRTMIN*0.5;
1828         return sqrt(u*u + v*v) * SQRTMAX * 2.0;
1829     }
1830     if (u <= SQRTMIN)
1831     {
1832         // hypot (tiny, tiny) -- avoid underflow
1833         // This is only necessary to avoid setting the underflow
1834         // flag.
1835         u *= SQRTMAX / real.epsilon;
1836         v *= SQRTMAX / real.epsilon;
1837         return sqrt(u*u + v*v) * SQRTMIN * real.epsilon;
1838     }
1839     if (u * real.epsilon > v)
1840     {
1841         // hypot (huge, tiny) = huge
1842         return u;
1843     }
1844
1845     // both are in the normal range
1846     return sqrt(u*u + v*v);
1847 }
1848
1849 unittest
1850 {
1851     static real vals[][3] =     // x,y,hypot
1852         [
1853             [ 0.0,     0.0,   0.0],
1854             [ 0.0,    -0.0,   0.0],
1855             [ -0.0,   -0.0,   0.0],
1856             [ 3.0,     4.0,   5.0],
1857             [ -300,   -400,   500],
1858             [0.0,      7.0,   7.0],
1859             [9.0,   9*real.epsilon,   9.0],
1860             [88/(64*sqrt(real.min_normal)), 105/(64*sqrt(real.min_normal)), 137/(64*sqrt(real.min_normal))],
1861             [88/(128*sqrt(real.min_normal)), 105/(128*sqrt(real.min_normal)), 137/(128*sqrt(real.min_normal))],
1862             [3*real.min_normal*real.epsilon, 4*real.min_normal*real.epsilon, 5*real.min_normal*real.epsilon],
1863             [ real.min_normal, real.min_normal, sqrt(2.0L)*real.min_normal],
1864             [ real.max/sqrt(2.0L), real.max/sqrt(2.0L), real.max],
1865             [ real.infinity, real.nan, real.infinity],
1866             [ real.nan, real.infinity, real.infinity],
1867             [ real.nan, real.nan, real.nan],
1868             [ real.nan, real.max, real.nan],
1869             [ real.max, real.nan, real.nan],
1870         ];
1871         for (int i = 0; i < vals.length; i++)
1872         {
1873             real x = vals[i][0];
1874             real y = vals[i][1];
1875             real z = vals[i][2];
1876             real h = hypot(x, y);
1877             assert(isIdentical(z, h));
1878         }
1879 }
1880
1881 /**********************************
1882  * Returns the error function of x.
1883  *
1884  * <img src="erf.gif" alt="error function">
1885  */
1886 real erf(real x)  @trusted nothrow   { return core.stdc.math.erfl(x); }
1887
1888 /**********************************
1889  * Returns the complementary error function of x, which is 1 - erf(x).
1890  *
1891  * <img src="erfc.gif" alt="complementary error function">
1892  */
1893 real erfc(real x)  @trusted nothrow  { return core.stdc.math.erfcl(x); }
1894
1895 /***********************************
1896  * Natural logarithm of gamma function.
1897  *
1898  * Returns the base e (2.718...) logarithm of the absolute
1899  * value of the gamma function of the argument.
1900  *
1901  * For reals, lgamma is equivalent to log(fabs(gamma(x))).
1902  *
1903  *      $(TABLE_SV
1904  *      $(TR $(TH x)                 $(TH lgamma(x)) $(TH invalid?))
1905  *      $(TR $(TD $(NAN))            $(TD $(NAN))    $(TD yes))
1906  *      $(TR $(TD integer $(LT)= 0)      $(TD +$(INFIN)) $(TD yes))
1907  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD +$(INFIN)) $(TD no))
1908  *      )
1909  */
1910 real lgamma(real x) @trusted nothrow
1911 {
1912     return core.stdc.math.lgammal(x);
1913
1914     // Use etc.gamma.lgamma for those C systems that are missing it
1915 }
1916
1917 /***********************************
1918  *  The Gamma function, $(GAMMA)(x)
1919  *
1920  *  $(GAMMA)(x) is a generalisation of the factorial function
1921  *  to real and complex numbers.
1922  *  Like x!, $(GAMMA)(x+1) = x*$(GAMMA)(x).
1923  *
1924  *  Mathematically, if z.re > 0 then
1925  *   $(GAMMA)(z) = $(INTEGRATE 0, $(INFIN)) $(POWER t, z-1)$(POWER e, -t) dt
1926  *
1927  *    $(TABLE_SV
1928  *      $(TR $(TH x)              $(TH $(GAMMA)(x))       $(TH invalid?))
1929  *      $(TR $(TD $(NAN))         $(TD $(NAN))            $(TD yes))
1930  *      $(TR $(TD $(PLUSMN)0.0)   $(TD $(PLUSMNINF))      $(TD yes))
1931  *      $(TR $(TD integer $(GT)0) $(TD (x-1)!)            $(TD no))
1932  *      $(TR $(TD integer $(LT)0) $(TD $(NAN))            $(TD yes))
1933  *      $(TR $(TD +$(INFIN))      $(TD +$(INFIN))         $(TD no))
1934  *      $(TR $(TD -$(INFIN))      $(TD $(NAN))            $(TD yes))
1935  *    )
1936  *
1937  *  References:
1938  *      $(LINK http://en.wikipedia.org/wiki/Gamma_function),
1939  *      $(LINK http://www.netlib.org/cephes/ldoubdoc.html#gamma)
1940  */
1941 real tgamma(real x) @trusted nothrow
1942 {
1943     return core.stdc.math.tgammal(x);
1944
1945     // Use etc.gamma.tgamma for those C systems that are missing it
1946 }
1947
1948 /**************************************
1949  * Returns the value of x rounded upward to the next integer
1950  * (toward positive infinity).
1951  */
1952 real ceil(real x)  @trusted nothrow    { return core.stdc.math.ceill(x); }
1953
1954 /**************************************
1955  * Returns the value of x rounded downward to the next integer
1956  * (toward negative infinity).
1957  */
1958 real floor(real x) @trusted nothrow    { return core.stdc.math.floorl(x); }
1959
1960 /******************************************
1961  * Rounds x to the nearest integer value, using the current rounding
1962  * mode.
1963  *
1964  * Unlike the rint functions, nearbyint does not raise the
1965  * FE_INEXACT exception.
1966  */
1967 real nearbyint(real x) @trusted nothrow { return core.stdc.math.nearbyintl(x); }
1968
1969 /**********************************
1970  * Rounds x to the nearest integer value, using the current rounding
1971  * mode.
1972  * If the return value is not equal to x, the FE_INEXACT
1973  * exception is raised.
1974  * $(B nearbyint) performs
1975  * the same operation, but does not set the FE_INEXACT exception.
1976  */
1977 real rint(real x) @safe pure nothrow;      /* intrinsic */
1978
1979 /***************************************
1980  * Rounds x to the nearest integer value, using the current rounding
1981  * mode.
1982  *
1983  * This is generally the fastest method to convert a floating-point number
1984  * to an integer. Note that the results from this function
1985  * depend on the rounding mode, if the fractional part of x is exactly 0.5.
1986  * If using the default rounding mode (ties round to even integers)
1987  * lrint(4.5) == 4, lrint(5.5)==6.
1988  */
1989 long lrint(real x) @trusted pure nothrow
1990 {
1991     version(InlineAsm_X86_Any)
1992     {
1993         long n;
1994         asm
1995         {
1996             fld x;
1997             fistp n;
1998         }
1999         return n;
2000     } else {
2001         return core.stdc.math.llrintl(x);
2002     }
2003 }
2004
2005 /*******************************************
2006  * Return the value of x rounded to the nearest integer.
2007  * If the fractional part of x is exactly 0.5, the return value is rounded to
2008  * the even integer.
2009  */
2010 real round(real x) @trusted nothrow { return core.stdc.math.roundl(x); }
2011
2012 /**********************************************
2013  * Return the value of x rounded to the nearest integer.
2014  *
2015  * If the fractional part of x is exactly 0.5, the return value is rounded
2016  * away from zero.
2017  */
2018 long lround(real x) @trusted nothrow
2019 {
2020     version (Posix)
2021         return core.stdc.math.llroundl(x);
2022     else
2023         assert (0, "lround not implemented");
2024 }
2025
2026 version(Posix)
2027 {
2028     unittest
2029     {
2030         assert(lround(0.49) == 0);
2031         assert(lround(0.5) == 1);
2032         assert(lround(1.5) == 2);
2033     }
2034 }
2035
2036 /****************************************************
2037  * Returns the integer portion of x, dropping the fractional portion.
2038  *
2039  * This is also known as "chop" rounding.
2040  */
2041 real trunc(real x) @trusted nothrow { return core.stdc.math.truncl(x); }
2042
2043 /****************************************************
2044  * Calculate the remainder x REM y, following IEC 60559.
2045  *
2046  * REM is the value of x - y * n, where n is the integer nearest the exact
2047  * value of x / y.
2048  * If |n - x / y| == 0.5, n is even.
2049  * If the result is zero, it has the same sign as x.
2050  * Otherwise, the sign of the result is the sign of x / y.
2051  * Precision mode has no effect on the remainder functions.
2052  *
2053  * remquo returns n in the parameter n.
2054  *
2055  * $(TABLE_SV
2056  *  $(TR $(TH x)               $(TH y)            $(TH remainder(x, y)) $(TH n)   $(TH invalid?))
2057  *  $(TR $(TD $(PLUSMN)0.0)    $(TD not 0.0)      $(TD $(PLUSMN)0.0)    $(TD 0.0) $(TD no))
2058  *  $(TR $(TD $(PLUSMNINF))    $(TD anything)     $(TD $(NAN))          $(TD ?)   $(TD yes))
2059  *  $(TR $(TD anything)        $(TD $(PLUSMN)0.0) $(TD $(NAN))          $(TD ?)   $(TD yes))
2060  *  $(TR $(TD != $(PLUSMNINF)) $(TD $(PLUSMNINF)) $(TD x)               $(TD ?)   $(TD no))
2061  * )
2062  *
2063  * Note: remquo not supported on windows
2064  */
2065 real remainder(real x, real y) @trusted nothrow { return core.stdc.math.remainderl(x, y); }
2066
2067 real remquo(real x, real y, out int n) @trusted nothrow  /// ditto
2068 {
2069     version (Posix)
2070         return core.stdc.math.remquol(x, y, &n);
2071     else
2072         assert (0, "remquo not implemented");
2073 }
2074
2075 /** IEEE exception status flags ('sticky bits')
2076
2077  These flags indicate that an exceptional floating-point condition has occurred.
2078  They indicate that a NaN or an infinity has been generated, that a result
2079  is inexact, or that a signalling NaN has been encountered. If floating-point
2080  exceptions are enabled (unmasked), a hardware exception will be generated
2081  instead of setting these flags.
2082
2083  Example:
2084  ----
2085     real a=3.5;
2086     // Set all the flags to zero
2087     resetIeeeFlags();
2088     assert(!ieeeFlags.divByZero);
2089     // Perform a division by zero.
2090     a/=0.0L;
2091     assert(a==real.infinity);
2092     assert(ieeeFlags.divByZero);
2093     // Create a NaN
2094     a*=0.0L;
2095     assert(ieeeFlags.invalid);
2096     assert(isNaN(a));
2097
2098     // Check that calling func() has no effect on the
2099     // status flags.
2100     IeeeFlags f = ieeeFlags;
2101     func();
2102     assert(ieeeFlags == f);
2103
2104  ----
2105  */
2106 struct IeeeFlags
2107 {
2108 private:
2109     // The x87 FPU status register is 16 bits.
2110     // The Pentium SSE2 status register is 32 bits.
2111     uint flags;
2112     version (X86_Any) {
2113         // Applies to both x87 status word (16 bits) and SSE2 status word(32 bits).
2114         enum : int {
2115             INEXACT_MASK   = 0x20,
2116             UNDERFLOW_MASK = 0x10,
2117             OVERFLOW_MASK  = 0x08,
2118             DIVBYZERO_MASK = 0x04,
2119             INVALID_MASK   = 0x01
2120         }
2121         // Don't bother about denormals, they are not supported on most CPUs.
2122         //  DENORMAL_MASK = 0x02;
2123     } else version (PPC) {
2124         // PowerPC FPSCR is a 32-bit register.
2125         enum : int {
2126             INEXACT_MASK   = 0x600,
2127             UNDERFLOW_MASK = 0x010,
2128             OVERFLOW_MASK  = 0x008,
2129             DIVBYZERO_MASK = 0x020,
2130             INVALID_MASK   = 0xF80 // PowerPC has five types of invalid exceptions.
2131         }
2132     } else version(SPARC) { // SPARC FSR is a 32bit register
2133              //(64 bits for Sparc 7 & 8, but high 32 bits are uninteresting).
2134         enum : int {
2135             INEXACT_MASK   = 0x020,
2136             UNDERFLOW_MASK = 0x080,
2137             OVERFLOW_MASK  = 0x100,
2138             DIVBYZERO_MASK = 0x040,
2139             INVALID_MASK   = 0x200
2140         }
2141     } else
2142         static assert(0, "Not implemented");
2143 private:
2144     static uint getIeeeFlags()
2145     {
2146         version(D_InlineAsm_X86) {
2147             asm {
2148                  fstsw AX;
2149                  // NOTE: If compiler supports SSE2, need to OR the result with
2150                  // the SSE2 status register.
2151                  // Clear all irrelevant bits
2152                  and EAX, 0x03D;
2153             }
2154         } else version(D_InlineAsm_X86_64) {
2155             asm {
2156                  fstsw AX;
2157                  // NOTE: If compiler supports SSE2, need to OR the result with
2158                  // the SSE2 status register.
2159                  // Clear all irrelevant bits
2160                  and RAX, 0x03D;
2161             }
2162         } else version (SPARC) {
2163            /*
2164                int retval;
2165                asm { st %fsr, retval; }
2166                return retval;
2167             */
2168            assert(0, "Not yet supported");
2169         } else
2170             assert(0, "Not yet supported");
2171     }
2172     static void resetIeeeFlags()
2173     {
2174         version(InlineAsm_X86_Any) {
2175             asm {
2176                 fnclex;
2177             }
2178         } else {
2179             /* SPARC:
2180               int tmpval;
2181               asm { st %fsr, tmpval; }
2182               tmpval &=0xFFFF_FC00;
2183               asm { ld tmpval, %fsr; }
2184             */
2185            assert(0, "Not yet supported");
2186         }
2187     }
2188 public:
2189      /// The result cannot be represented exactly, so rounding occured.
2190      /// (example: x = sin(0.1); )
2191      bool inexact() { return (flags & INEXACT_MASK) != 0; }
2192      /// A zero was generated by underflow (example: x = real.min*real.epsilon/2;)
2193      bool underflow() { return (flags & UNDERFLOW_MASK) != 0; }
2194      /// An infinity was generated by overflow (example: x = real.max*2;)
2195      bool overflow() { return (flags & OVERFLOW_MASK) != 0; }
2196      /// An infinity was generated by division by zero (example: x = 3/0.0; )
2197      bool divByZero() { return (flags & DIVBYZERO_MASK) != 0; }
2198      /// A machine NaN was generated. (example: x = real.infinity * 0.0; )
2199      bool invalid() { return (flags & INVALID_MASK) != 0; }
2200 }
2201
2202
2203 /// Set all of the floating-point status flags to false.
2204 void resetIeeeFlags() { IeeeFlags.resetIeeeFlags; }
2205
2206 /// Return a snapshot of the current state of the floating-point status flags.
2207 IeeeFlags ieeeFlags()
2208 {
2209    return IeeeFlags(IeeeFlags.getIeeeFlags());
2210 }
2211
2212 /** Control the Floating point hardware
2213
2214   Change the IEEE754 floating-point rounding mode and the floating-point
2215   hardware exceptions.
2216
2217   By default, the rounding mode is roundToNearest and all hardware exceptions
2218   are disabled. For most applications, debugging is easier if the $(I division
2219   by zero), $(I overflow), and $(I invalid operation) exceptions are enabled.
2220   These three are combined into a $(I severeExceptions) value for convenience.
2221   Note in particular that if $(I invalidException) is enabled, a hardware trap
2222   will be generated whenever an uninitialized floating-point variable is used.
2223
2224   All changes are temporary. The previous state is restored at the
2225   end of the scope.
2226
2227
2228 Example:
2229  ----
2230   {
2231     // Enable hardware exceptions for division by zero, overflow to infinity,
2232     // invalid operations, and uninitialized floating-point variables.
2233
2234     FloatingPointControl fpctrl;
2235     fpctrl.enableExceptions(FloatingPointControl.severeExceptions);
2236
2237     double y = x*3.0; // will generate a hardware exception, if x is uninitialized.
2238     //
2239     fpctrl.rounding = FloatingPointControl.roundUp;
2240
2241     // The hardware exceptions will be disabled when leaving this scope.
2242     // The original rounding mode will also be restored.
2243   }
2244
2245  ----
2246
2247  */
2248 struct FloatingPointControl
2249 {
2250     alias uint RoundingMode;
2251
2252     /** IEEE rounding modes.
2253      * The default mode is roundToNearest.
2254      */
2255     enum : RoundingMode
2256     {
2257         roundToNearest = 0x0000,
2258         roundDown      = 0x0400,
2259         roundUp        = 0x0800,
2260         roundToZero    = 0x0C00
2261     };
2262     /** IEEE hardware exceptions.
2263      *  By default, all exceptions are masked (disabled).
2264      */
2265     enum : uint
2266     {
2267         inexactException   = 0x20,
2268         underflowException = 0x10,
2269         overflowException  = 0x08,
2270         divByZeroException = 0x04,
2271         invalidException   = 0x01,
2272         /// Severe = The overflow, division by zero, and invalid exceptions.
2273         severeExceptions   = overflowException | divByZeroException
2274                              | invalidException,
2275         allExceptions      = severeExceptions | underflowException
2276                              | inexactException,
2277     };
2278 private:
2279     enum ushort EXCEPTION_MASK = 0x3F;
2280     enum ushort ROUNDING_MASK = 0xC00;
2281 public:
2282     /// Enable (unmask) specific hardware exceptions. Multiple exceptions may be ORed together.
2283     void enableExceptions(uint exceptions)
2284     {
2285         initialize();
2286         setControlState(getControlState() & ~(exceptions & EXCEPTION_MASK));
2287     }
2288     /// Disable (mask) specific hardware exceptions. Multiple exceptions may be ORed together.
2289     void disableExceptions(uint exceptions)
2290     {
2291         initialize();
2292         setControlState(getControlState() | (exceptions & EXCEPTION_MASK));
2293     }
2294     //// Change the floating-point hardware rounding mode
2295     void rounding(RoundingMode newMode)
2296     {
2297         ushort old = getControlState();
2298         setControlState((old & ~ROUNDING_MASK) | (newMode & ROUNDING_MASK));
2299     }
2300     /// Return the exceptions which are currently enabled (unmasked)
2301     static uint enabledExceptions()
2302     {
2303         return (getControlState() & EXCEPTION_MASK) ^ EXCEPTION_MASK;
2304     }
2305     /// Return the currently active rounding mode
2306     static RoundingMode rounding()
2307     {
2308         return cast(RoundingMode)(getControlState() & ROUNDING_MASK);
2309     }
2310     ///  Clear all pending exceptions, then restore the original exception state and rounding mode.
2311     ~this()
2312     {
2313         clearExceptions();
2314         setControlState(savedState);
2315     }
2316 private:
2317     ushort savedState;
2318
2319     bool initialized=false;
2320     void initialize()
2321     {
2322         // BUG: This works around the absence of this() constructors.
2323         if (initialized) return;
2324         clearExceptions();
2325         savedState = getControlState();
2326         initialized=true;
2327     }
2328     // Clear all pending exceptions
2329     static void clearExceptions()
2330     {
2331         version (InlineAsm_X86_Any)
2332         {
2333             asm
2334             {
2335                 fclex;
2336             }
2337         }
2338         else
2339             assert(0, "Not yet supported");
2340     }
2341     // Read from the control register
2342     static ushort getControlState()
2343     {
2344         version (D_InlineAsm_X86)
2345         {
2346             short cont;
2347             asm
2348             {
2349                 xor EAX, EAX;
2350                 fstcw cont;
2351             }
2352             return cont;
2353         }
2354         else
2355         version (D_InlineAsm_X86_64)
2356         {
2357             short cont;
2358             asm
2359             {
2360                 xor RAX, RAX;
2361                 fstcw cont;
2362             }
2363             return cont;
2364         }
2365         else
2366             assert(0, "Not yet supported");
2367     }
2368     // Set the control register
2369     static void setControlState(ushort newState)
2370     {
2371         version (InlineAsm_X86_Any)
2372         {
2373             asm
2374             {
2375                  fclex;
2376                  fldcw newState;
2377             }
2378         }
2379         else
2380             assert(0, "Not yet supported");
2381     }
2382 }
2383
2384 unittest
2385 {
2386    {
2387         FloatingPointControl ctrl;
2388         ctrl.enableExceptions(FloatingPointControl.divByZeroException
2389                            | FloatingPointControl.overflowException);
2390         assert(ctrl.enabledExceptions() ==
2391             (FloatingPointControl.divByZeroException
2392           | FloatingPointControl.overflowException));
2393
2394         ctrl.rounding = FloatingPointControl.roundUp;
2395         assert(FloatingPointControl.rounding == FloatingPointControl.roundUp);
2396     }
2397     assert(FloatingPointControl.rounding
2398        == FloatingPointControl.roundToNearest);
2399     assert(FloatingPointControl.enabledExceptions() ==0);
2400 }
2401
2402
2403 /*********************************
2404  * Returns !=0 if e is a NaN.
2405  */
2406
2407 bool isNaN(real x) @trusted pure nothrow
2408 {
2409     alias floatTraits!(real) F;
2410     static if (real.mant_dig==53) { // double
2411         ulong*  p = cast(ulong *)&x;
2412         return ((*p & 0x7FF0_0000_0000_0000) == 0x7FF0_0000_0000_0000)
2413         && *p & 0x000F_FFFF_FFFF_FFFF;
2414     } else static if (real.mant_dig==64) {     // real80
2415         ushort e = F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT];
2416         ulong*  ps = cast(ulong *)&x;
2417         return e == F.EXPMASK &&
2418         *ps & 0x7FFF_FFFF_FFFF_FFFF; // not infinity
2419     } else static if (real.mant_dig==113) {  // quadruple
2420         ushort e = F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT];
2421         ulong*  ps = cast(ulong *)&x;
2422         return e == F.EXPMASK &&
2423         (ps[MANTISSA_LSB] | (ps[MANTISSA_MSB]& 0x0000_FFFF_FFFF_FFFF))!=0;
2424     } else {
2425         return x!=x;
2426     }
2427 }
2428
2429
2430 unittest
2431 {
2432     assert(isNaN(float.nan));
2433     assert(isNaN(-double.nan));
2434     assert(isNaN(real.nan));
2435
2436     assert(!isNaN(53.6));
2437     assert(!isNaN(float.infinity));
2438 }
2439
2440 /*********************************
2441  * Returns !=0 if e is finite (not infinite or $(NAN)).
2442  */
2443
2444 int isFinite(real e) @trusted pure nothrow
2445 {
2446     alias floatTraits!(real) F;
2447     ushort* pe = cast(ushort *)&e;
2448     return (pe[F.EXPPOS_SHORT] & F.EXPMASK) != F.EXPMASK;
2449 }
2450
2451 unittest
2452 {
2453     assert(isFinite(1.23));
2454     assert(!isFinite(double.infinity));
2455     assert(!isFinite(float.nan));
2456 }
2457
2458
2459 /*********************************
2460  * Returns !=0 if x is normalized (not zero, subnormal, infinite, or $(NAN)).
2461  */
2462
2463 /* Need one for each format because subnormal floats might
2464  * be converted to normal reals.
2465  */
2466
2467 int isNormal(X)(X x) @trusted pure nothrow
2468 {
2469     alias floatTraits!(X) F;
2470
2471     static if(real.mant_dig==106) { // doubledouble
2472         // doubledouble is normal if the least significant part is normal.
2473         return isNormal((cast(double*)&x)[MANTISSA_LSB]);
2474     } else {
2475         ushort e = F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT];
2476         return (e != F.EXPMASK && e!=0);
2477     }
2478 }
2479
2480
2481 unittest
2482 {
2483     float f = 3;
2484     double d = 500;
2485     real e = 10e+48;
2486
2487     assert(isNormal(f));
2488     assert(isNormal(d));
2489     assert(isNormal(e));
2490     f = d = e = 0;
2491     assert(!isNormal(f));
2492     assert(!isNormal(d));
2493     assert(!isNormal(e));
2494     assert(!isNormal(real.infinity));
2495     assert(isNormal(-real.max));
2496     assert(!isNormal(real.min_normal/4));
2497
2498 }
2499
2500 /*********************************
2501  * Is number subnormal? (Also called "denormal".)
2502  * Subnormals have a 0 exponent and a 0 most significant mantissa bit.
2503  */
2504
2505 /* Need one for each format because subnormal floats might
2506  * be converted to normal reals.
2507  */
2508
2509 int isSubnormal(float f) @trusted pure nothrow
2510 {
2511     uint *p = cast(uint *)&f;
2512     return (*p & 0x7F80_0000) == 0 && *p & 0x007F_FFFF;
2513 }
2514
2515 unittest
2516 {
2517     float f = 3.0;
2518
2519     for (f = 1.0; !isSubnormal(f); f /= 2)
2520         assert(f != 0);
2521 }
2522
2523 /// ditto
2524
2525 int isSubnormal(double d) @trusted pure nothrow
2526 {
2527     uint *p = cast(uint *)&d;
2528     return (p[MANTISSA_MSB] & 0x7FF0_0000) == 0
2529         && (p[MANTISSA_LSB] || p[MANTISSA_MSB] & 0x000F_FFFF);
2530 }
2531
2532 unittest
2533 {
2534     double f;
2535
2536     for (f = 1; !isSubnormal(f); f /= 2)
2537         assert(f != 0);
2538 }
2539
2540 /// ditto
2541
2542 int isSubnormal(real x) @trusted pure nothrow
2543 {
2544     alias floatTraits!(real) F;
2545     static if (real.mant_dig == 53) { // double
2546         return isSubnormal(cast(double)x);
2547     } else static if (real.mant_dig == 113) { // quadruple
2548         ushort e = F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT];
2549         long*   ps = cast(long *)&x;
2550         return (e == 0 &&
2551           (((ps[MANTISSA_LSB]|(ps[MANTISSA_MSB]& 0x0000_FFFF_FFFF_FFFF))) !=0));
2552     } else static if (real.mant_dig==64) { // real80
2553         ushort* pe = cast(ushort *)&x;
2554         long*   ps = cast(long *)&x;
2555
2556         return (pe[F.EXPPOS_SHORT] & F.EXPMASK) == 0 && *ps > 0;
2557     } else { // double double
2558         return isSubnormal((cast(double*)&x)[MANTISSA_MSB]);
2559     }
2560 }
2561
2562 unittest
2563 {
2564     real f;
2565
2566     for (f = 1; !isSubnormal(f); f /= 2)
2567         assert(f != 0);
2568 }
2569
2570 /*********************************
2571  * Return !=0 if e is $(PLUSMN)$(INFIN).
2572  */
2573
2574 bool isInfinity(real x) @trusted pure nothrow
2575 {
2576     alias floatTraits!(real) F;
2577     static if (real.mant_dig == 53) { // double
2578         return ((*cast(ulong *)&x) & 0x7FFF_FFFF_FFFF_FFFF)
2579             == 0x7FF8_0000_0000_0000;
2580     } else static if(real.mant_dig == 106) { //doubledouble
2581         return (((cast(ulong *)&x)[MANTISSA_MSB]) & 0x7FFF_FFFF_FFFF_FFFF)
2582             == 0x7FF8_0000_0000_0000;
2583     } else static if (real.mant_dig == 113) { // quadruple
2584         long*   ps = cast(long *)&x;
2585         return (ps[MANTISSA_LSB] == 0)
2586             && (ps[MANTISSA_MSB] & 0x7FFF_FFFF_FFFF_FFFF) == 0x7FFF_0000_0000_0000;
2587     } else { // real80
2588         ushort e = cast(ushort)(F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT]);
2589         ulong*  ps = cast(ulong *)&x;
2590         // On Motorola 68K, infinity can have hidden bit=1 or 0. On x86, it is always 1.
2591         return e == F.EXPMASK && (*ps & 0x7FFF_FFFF_FFFF_FFFF) == 0;
2592     }
2593 }
2594
2595 unittest
2596 {
2597     assert(isInfinity(float.infinity));
2598     assert(!isInfinity(float.nan));
2599     assert(isInfinity(double.infinity));
2600     assert(isInfinity(-real.infinity));
2601
2602     assert(isInfinity(-1.0 / 0.0));
2603 }
2604
2605 /*********************************
2606  * Is the binary representation of x identical to y?
2607  *
2608  * Same as ==, except that positive and negative zero are not identical,
2609  * and two $(NAN)s are identical if they have the same 'payload'.
2610  */
2611
2612 bool isIdentical(real x, real y) @trusted pure nothrow
2613 {
2614     // We're doing a bitwise comparison so the endianness is irrelevant.
2615     long*   pxs = cast(long *)&x;
2616     long*   pys = cast(long *)&y;
2617     static if (real.mant_dig == 53)
2618     { //double
2619         return pxs[0] == pys[0];
2620     }
2621     else static if (real.mant_dig == 113 || real.mant_dig==106)
2622     {
2623         // quadruple or doubledouble
2624         return pxs[0] == pys[0] && pxs[1] == pys[1];
2625     }
2626     else
2627     { // real80
2628         ushort* pxe = cast(ushort *)&x;
2629         ushort* pye = cast(ushort *)&y;
2630         return pxe[4] == pye[4] && pxs[0] == pys[0];
2631     }
2632 }
2633
2634 /*********************************
2635  * Return 1 if sign bit of e is set, 0 if not.
2636  */
2637
2638 int signbit(real x) @trusted pure nothrow
2639 {
2640     return ((cast(ubyte *)&x)[floatTraits!(real).SIGNPOS_BYTE] & 0x80) != 0;
2641 }
2642
2643 unittest
2644 {
2645     debug (math) printf("math.signbit.unittest\n");
2646     assert(!signbit(float.nan));
2647     assert(signbit(-float.nan));
2648     assert(!signbit(168.1234));
2649     assert(signbit(-168.1234));
2650     assert(!signbit(0.0));
2651     assert(signbit(-0.0));
2652     assert(signbit(-double.max));
2653     assert(!signbit(double.max));
2654 }
2655
2656 /*********************************
2657  * Return a value composed of to with from's sign bit.
2658  */
2659
2660 real copysign(real to, real from) @trusted pure nothrow
2661 {
2662     ubyte* pto   = cast(ubyte *)&to;
2663     const ubyte* pfrom = cast(ubyte *)&from;
2664
2665     alias floatTraits!(real) F;
2666     pto[F.SIGNPOS_BYTE] &= 0x7F;
2667     pto[F.SIGNPOS_BYTE] |= pfrom[F.SIGNPOS_BYTE] & 0x80;
2668     return to;
2669 }
2670
2671 unittest
2672 {
2673     real e;
2674
2675     e = copysign(21, 23.8);
2676     assert(e == 21);
2677
2678     e = copysign(-21, 23.8);
2679     assert(e == 21);
2680
2681     e = copysign(21, -23.8);
2682     assert(e == -21);
2683
2684     e = copysign(-21, -23.8);
2685     assert(e == -21);
2686
2687     e = copysign(real.nan, -23.8);
2688     assert(isNaN(e) && signbit(e));
2689 }
2690
2691 /*********************************
2692 Returns $(D -1) if $(D x < 0), $(D x) if $(D x == 0), $(D 1) if
2693 $(D x > 0), and $(NAN) if x==$(NAN).
2694  */
2695 F sgn(F)(F x) @safe pure nothrow
2696 {
2697     // @@@TODO@@@: make this faster
2698     return x > 0 ? 1 : x < 0 ? -1 : x;
2699 }
2700
2701 unittest
2702 {
2703     debug (math) printf("math.sgn.unittest\n");
2704     assert(sgn(168.1234) == 1);
2705     assert(sgn(-168.1234) == -1);
2706     assert(sgn(0.0) == 0);
2707     assert(sgn(-0.0) == 0);
2708 }
2709
2710 // Functions for NaN payloads
2711 /*
2712  * A 'payload' can be stored in the significand of a $(NAN). One bit is required
2713  * to distinguish between a quiet and a signalling $(NAN). This leaves 22 bits
2714  * of payload for a float; 51 bits for a double; 62 bits for an 80-bit real;
2715  * and 111 bits for a 128-bit quad.
2716 */
2717 /**
2718  * Create a quiet $(NAN), storing an integer inside the payload.
2719  *
2720  * For floats, the largest possible payload is 0x3F_FFFF.
2721  * For doubles, it is 0x3_FFFF_FFFF_FFFF.
2722  * For 80-bit or 128-bit reals, it is 0x3FFF_FFFF_FFFF_FFFF.
2723  */
2724 real NaN(ulong payload) @trusted pure nothrow
2725 {
2726     static if (real.mant_dig == 64) { //real80
2727         ulong v = 3; // implied bit = 1, quiet bit = 1
2728     } else {
2729         ulong v = 2; // no implied bit. quiet bit = 1
2730     }
2731
2732     ulong a = payload;
2733
2734     // 22 Float bits
2735     ulong w = a & 0x3F_FFFF;
2736     a -= w;
2737
2738     v <<=22;
2739     v |= w;
2740     a >>=22;
2741
2742     // 29 Double bits
2743     v <<=29;
2744     w = a & 0xFFF_FFFF;
2745     v |= w;
2746     a -= w;
2747     a >>=29;
2748
2749     static if (real.mant_dig == 53) { // double
2750         v |=0x7FF0_0000_0000_0000;
2751         real x;
2752         * cast(ulong *)(&x) = v;
2753         return x;
2754     } else {
2755         v <<=11;
2756         a &= 0x7FF;
2757         v |= a;
2758         real x = real.nan;
2759         // Extended real bits
2760         static if (real.mant_dig==113) { //quadruple
2761             v<<=1; // there's no implicit bit
2762             version(LittleEndian) {
2763                 *cast(ulong*)(6+cast(ubyte*)(&x)) = v;
2764             } else {
2765                 *cast(ulong*)(2+cast(ubyte*)(&x)) = v;
2766             }
2767         } else { // real80
2768             * cast(ulong *)(&x) = v;
2769         }
2770         return x;
2771     }
2772 }
2773
2774 /**
2775  * Extract an integral payload from a $(NAN).
2776  *
2777  * Returns:
2778  * the integer payload as a ulong.
2779  *
2780  * For floats, the largest possible payload is 0x3F_FFFF.
2781  * For doubles, it is 0x3_FFFF_FFFF_FFFF.
2782  * For 80-bit or 128-bit reals, it is 0x3FFF_FFFF_FFFF_FFFF.
2783  */
2784 ulong getNaNPayload(real x) @trusted pure nothrow
2785 {
2786     //  assert(isNaN(x));
2787     static if (real.mant_dig == 53) {
2788         ulong m = *cast(ulong *)(&x);
2789         // Make it look like an 80-bit significand.
2790         // Skip exponent, and quiet bit
2791         m &= 0x0007_FFFF_FFFF_FFFF;
2792         m <<= 10;
2793     } else static if (real.mant_dig==113) { // quadruple
2794         version(LittleEndian) {
2795             ulong m = *cast(ulong*)(6+cast(ubyte*)(&x));
2796         } else {
2797             ulong m = *cast(ulong*)(2+cast(ubyte*)(&x));
2798         }
2799         m>>=1; // there's no implicit bit
2800     } else {
2801         ulong m = *cast(ulong *)(&x);
2802     }
2803     // ignore implicit bit and quiet bit
2804     ulong f = m & 0x3FFF_FF00_0000_0000L;
2805     ulong w = f >>> 40;
2806             w |= (m & 0x00FF_FFFF_F800L) << (22 - 11);
2807             w |= (m & 0x7FF) << 51;
2808             return w;
2809 }
2810
2811 debug(UnitTest) {
2812     unittest {
2813         real nan4 = NaN(0x789_ABCD_EF12_3456);
2814         static if (real.mant_dig == 64 || real.mant_dig==113) {
2815             assert (getNaNPayload(nan4) == 0x789_ABCD_EF12_3456);
2816         } else {
2817             assert (getNaNPayload(nan4) == 0x1_ABCD_EF12_3456);
2818         }
2819         double nan5 = nan4;
2820         assert (getNaNPayload(nan5) == 0x1_ABCD_EF12_3456);
2821         float nan6 = nan4;
2822         assert (getNaNPayload(nan6) == 0x12_3456);
2823         nan4 = NaN(0xFABCD);
2824         assert (getNaNPayload(nan4) == 0xFABCD);
2825         nan6 = nan4;
2826         assert (getNaNPayload(nan6) == 0xFABCD);
2827         nan5 = NaN(0x100_0000_0000_3456);
2828         assert(getNaNPayload(nan5) == 0x0000_0000_3456);
2829     }
2830 }
2831
2832 /**
2833  * Calculate the next largest floating point value after x.
2834  *
2835  * Return the least number greater than x that is representable as a real;
2836  * thus, it gives the next point on the IEEE number line.
2837  *
2838  *  $(TABLE_SV
2839  *    $(SVH x,            nextUp(x)   )
2840  *    $(SV  -$(INFIN),    -real.max   )
2841  *    $(SV  $(PLUSMN)0.0, real.min_normal*real.epsilon )
2842  *    $(SV  real.max,     $(INFIN) )
2843  *    $(SV  $(INFIN),     $(INFIN) )
2844  *    $(SV  $(NAN),       $(NAN)   )
2845  * )
2846  */
2847 real nextUp(real x) @trusted pure nothrow
2848 {
2849     alias floatTraits!(real) F;
2850     static if (real.mant_dig == 53) { // double
2851         return nextUp(cast(double)x);
2852     } else static if(real.mant_dig==113) {  // quadruple
2853         ushort e = F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT];
2854         if (e == F.EXPMASK) { // NaN or Infinity
2855             if (x == -real.infinity) return -real.max;
2856             return x; // +Inf and NaN are unchanged.
2857         }
2858         ulong*   ps = cast(ulong *)&e;
2859         if (ps[MANTISSA_LSB] & 0x8000_0000_0000_0000)  { // Negative number
2860             if (ps[MANTISSA_LSB] == 0
2861                 && ps[MANTISSA_MSB] == 0x8000_0000_0000_0000) {
2862                 // it was negative zero, change to smallest subnormal
2863                 ps[MANTISSA_LSB] = 0x0000_0000_0000_0001;
2864                 ps[MANTISSA_MSB] = 0;
2865                 return x;
2866             }
2867             --*ps;
2868             if (ps[MANTISSA_LSB]==0) --ps[MANTISSA_MSB];
2869         } else { // Positive number
2870             ++ps[MANTISSA_LSB];
2871             if (ps[MANTISSA_LSB]==0) ++ps[MANTISSA_MSB];
2872         }
2873         return x;
2874
2875     } else static if(real.mant_dig==64){ // real80
2876         // For 80-bit reals, the "implied bit" is a nuisance...
2877         ushort *pe = cast(ushort *)&x;
2878         ulong  *ps = cast(ulong  *)&x;
2879
2880         if ((pe[F.EXPPOS_SHORT] & F.EXPMASK) == F.EXPMASK) {
2881             // First, deal with NANs and infinity
2882             if (x == -real.infinity) return -real.max;
2883             return x; // +Inf and NaN are unchanged.
2884         }
2885         if (pe[F.EXPPOS_SHORT] & 0x8000)  {
2886             // Negative number -- need to decrease the significand
2887             --*ps;
2888             // Need to mask with 0x7FFF... so subnormals are treated correctly.
2889             if ((*ps & 0x7FFF_FFFF_FFFF_FFFF) == 0x7FFF_FFFF_FFFF_FFFF) {
2890                 if (pe[F.EXPPOS_SHORT] == 0x8000) { // it was negative zero
2891                     *ps = 1;
2892                     pe[F.EXPPOS_SHORT] = 0; // smallest subnormal.
2893                     return x;
2894                 }
2895                 --pe[F.EXPPOS_SHORT];
2896                 if (pe[F.EXPPOS_SHORT] == 0x8000) {
2897                     return x; // it's become a subnormal, implied bit stays low.
2898                 }
2899                 *ps = 0xFFFF_FFFF_FFFF_FFFF; // set the implied bit
2900                 return x;
2901             }
2902             return x;
2903         } else {
2904             // Positive number -- need to increase the significand.
2905             // Works automatically for positive zero.
2906             ++*ps;
2907             if ((*ps & 0x7FFF_FFFF_FFFF_FFFF) == 0) {
2908                 // change in exponent
2909                 ++pe[F.EXPPOS_SHORT];
2910                 *ps = 0x8000_0000_0000_0000; // set the high bit
2911             }
2912         }
2913         return x;
2914     } // doubledouble is not supported
2915 }
2916
2917 /** ditto */
2918 double nextUp(double x) @trusted pure nothrow
2919 {
2920     ulong *ps = cast(ulong *)&x;
2921
2922     if ((*ps & 0x7FF0_0000_0000_0000) == 0x7FF0_0000_0000_0000) {
2923         // First, deal with NANs and infinity
2924         if (x == -x.infinity) return -x.max;
2925         return x; // +INF and NAN are unchanged.
2926     }
2927     if (*ps & 0x8000_0000_0000_0000)  { // Negative number
2928         if (*ps == 0x8000_0000_0000_0000) { // it was negative zero
2929             *ps = 0x0000_0000_0000_0001; // change to smallest subnormal
2930             return x;
2931         }
2932         --*ps;
2933     } else { // Positive number
2934         ++*ps;
2935     }
2936     return x;
2937 }
2938
2939 /** ditto */
2940 float nextUp(float x) @trusted pure nothrow
2941 {
2942     uint *ps = cast(uint *)&x;
2943
2944     if ((*ps & 0x7F80_0000) == 0x7F80_0000) {
2945         // First, deal with NANs and infinity
2946         if (x == -x.infinity) return -x.max;
2947         return x; // +INF and NAN are unchanged.
2948     }
2949     if (*ps & 0x8000_0000)  { // Negative number
2950         if (*ps == 0x8000_0000) { // it was negative zero
2951             *ps = 0x0000_0001; // change to smallest subnormal
2952             return x;
2953         }
2954         --*ps;
2955     } else { // Positive number
2956         ++*ps;
2957     }
2958     return x;
2959 }
2960
2961 /**
2962  * Calculate the next smallest floating point value before x.
2963  *
2964  * Return the greatest number less than x that is representable as a real;
2965  * thus, it gives the previous point on the IEEE number line.
2966  *
2967  *  $(TABLE_SV
2968  *    $(SVH x,            nextDown(x)   )
2969  *    $(SV  $(INFIN),     real.max  )
2970  *    $(SV  $(PLUSMN)0.0, -real.min_normal*real.epsilon )
2971  *    $(SV  -real.max,    -$(INFIN) )
2972  *    $(SV  -$(INFIN),    -$(INFIN) )
2973  *    $(SV  $(NAN),       $(NAN)    )
2974  * )
2975  */
2976 real nextDown(real x) @safe pure nothrow
2977 {
2978     return -nextUp(-x);
2979 }
2980
2981 /** ditto */
2982 double nextDown(double x) @safe pure nothrow
2983 {
2984     return -nextUp(-x);
2985 }
2986
2987 /** ditto */
2988 float nextDown(float x) @safe pure nothrow
2989 {
2990     return -nextUp(-x);
2991 }
2992
2993 unittest {
2994     assert( nextDown(1.0 + real.epsilon) == 1.0);
2995 }
2996
2997 unittest {
2998     static if (real.mant_dig == 64) {
2999
3000         // Tests for 80-bit reals
3001         assert(isIdentical(nextUp(NaN(0xABC)), NaN(0xABC)));
3002         // negative numbers
3003         assert( nextUp(-real.infinity) == -real.max );
3004         assert( nextUp(-1.0L-real.epsilon) == -1.0 );
3005         assert( nextUp(-2.0L) == -2.0 + real.epsilon);
3006         // denormals and zero
3007         assert( nextUp(-real.min_normal) == -real.min_normal*(1-real.epsilon) );
3008         assert( nextUp(-real.min_normal*(1-real.epsilon)) == -real.min_normal*(1-2*real.epsilon) );
3009         assert( isIdentical(-0.0L, nextUp(-real.min_normal*real.epsilon)) );
3010         assert( nextUp(-0.0L) == real.min_normal*real.epsilon );
3011         assert( nextUp(0.0L) == real.min_normal*real.epsilon );
3012         assert( nextUp(real.min_normal*(1-real.epsilon)) == real.min_normal );
3013         assert( nextUp(real.min_normal) == real.min_normal*(1+real.epsilon) );
3014         // positive numbers
3015         assert( nextUp(1.0L) == 1.0 + real.epsilon );
3016         assert( nextUp(2.0L-real.epsilon) == 2.0 );
3017         assert( nextUp(real.max) == real.infinity );
3018         assert( nextUp(real.infinity)==real.infinity );
3019     }
3020
3021     double n = NaN(0xABC);
3022     assert(isIdentical(nextUp(n), n));
3023     // negative numbers
3024     assert( nextUp(-double.infinity) == -double.max );
3025     assert( nextUp(-1-double.epsilon) == -1.0 );
3026     assert( nextUp(-2.0) == -2.0 + double.epsilon);
3027     // denormals and zero
3028
3029     assert( nextUp(-double.min_normal) == -double.min_normal*(1-double.epsilon) );
3030     assert( nextUp(-double.min_normal*(1-double.epsilon)) == -double.min_normal*(1-2*double.epsilon) );
3031     assert( isIdentical(-0.0, nextUp(-double.min_normal*double.epsilon)) );
3032     assert( nextUp(0.0) == double.min_normal*double.epsilon );
3033     assert( nextUp(-0.0) == double.min_normal*double.epsilon );
3034     assert( nextUp(double.min_normal*(1-double.epsilon)) == double.min_normal );
3035     assert( nextUp(double.min_normal) == double.min_normal*(1+double.epsilon) );
3036     // positive numbers
3037     assert( nextUp(1.0) == 1.0 + double.epsilon );
3038     assert( nextUp(2.0-double.epsilon) == 2.0 );
3039     assert( nextUp(double.max) == double.infinity );
3040
3041     float fn = NaN(0xABC);
3042     assert(isIdentical(nextUp(fn), fn));
3043     float f = -float.min_normal*(1-float.epsilon);
3044     float f1 = -float.min_normal;
3045     assert( nextUp(f1) ==  f);
3046     f = 1.0f+float.epsilon;
3047     f1 = 1.0f;
3048     assert( nextUp(f1) == f );
3049     f1 = -0.0f;
3050     assert( nextUp(f1) == float.min_normal*float.epsilon);
3051     assert( nextUp(float.infinity)==float.infinity );
3052
3053     assert(nextDown(1.0L+real.epsilon)==1.0);
3054     assert(nextDown(1.0+double.epsilon)==1.0);
3055     f = 1.0f+float.epsilon;
3056     assert(nextDown(f)==1.0);
3057     assert(nextafter(1.0+real.epsilon, -real.infinity)==1.0);
3058 }
3059
3060
3061
3062 /******************************************
3063  * Calculates the next representable value after x in the direction of y.
3064  *
3065  * If y > x, the result will be the next largest floating-point value;
3066  * if y < x, the result will be the next smallest value.
3067  * If x == y, the result is y.
3068  *
3069  * Remarks:
3070  * This function is not generally very useful; it's almost always better to use
3071  * the faster functions nextUp() or nextDown() instead.
3072  *
3073  * The FE_INEXACT and FE_OVERFLOW exceptions will be raised if x is finite and
3074  * the function result is infinite. The FE_INEXACT and FE_UNDERFLOW
3075  * exceptions will be raised if the function value is subnormal, and x is
3076  * not equal to y.
3077  */
3078 T nextafter(T)(T x, T y) @safe pure nothrow
3079 {
3080     if (x==y) return y;
3081     return ((y>x) ? nextUp(x) :  nextDown(x));
3082 }
3083
3084 unittest
3085 {
3086     float a = 1;
3087     assert(is(typeof(nextafter(a, a)) == float));
3088     assert(nextafter(a, a.infinity) > a);
3089
3090     double b = 2;
3091     assert(is(typeof(nextafter(b, b)) == double));
3092     assert(nextafter(b, b.infinity) > b);
3093
3094     real c = 3;
3095     assert(is(typeof(nextafter(c, c)) == real));
3096     assert(nextafter(c, c.infinity) > c);
3097 }
3098
3099 //real nexttoward(real x, real y) { return core.stdc.math.nexttowardl(x, y); }
3100
3101 /*******************************************
3102  * Returns the positive difference between x and y.
3103  * Returns:
3104  *      $(TABLE_SV
3105  *      $(TR $(TH x, y)       $(TH fdim(x, y)))
3106  *      $(TR $(TD x $(GT) y)  $(TD x - y))
3107  *      $(TR $(TD x $(LT)= y) $(TD +0.0))
3108  *      )
3109  */
3110 real fdim(real x, real y) @safe pure nothrow { return (x > y) ? x - y : +0.0; }
3111
3112 /****************************************
3113  * Returns the larger of x and y.
3114  */
3115 real fmax(real x, real y) @safe pure nothrow { return x > y ? x : y; }
3116
3117 /****************************************
3118  * Returns the smaller of x and y.
3119  */
3120 real fmin(real x, real y) @safe pure nothrow { return x < y ? x : y; }
3121
3122 /**************************************
3123  * Returns (x * y) + z, rounding only once according to the
3124  * current rounding mode.
3125  *
3126  * BUGS: Not currently implemented - rounds twice.
3127  */
3128 real fma(real x, real y, real z) @safe pure nothrow { return (x * y) + z; }
3129
3130 /*******************************************************************
3131  * Compute the value of x $(SUP n), where n is an integer
3132  */
3133 Unqual!F pow(F, G)(F x, G n) @trusted pure nothrow
3134     if (isFloatingPoint!(F) && isIntegral!(G))
3135 {
3136     real p = 1.0, v = void;
3137     Unsigned!(Unqual!G) m = n;
3138     if (n < 0)
3139     {
3140         switch (n)
3141         {
3142         case -1:
3143             return 1 / x;
3144         case -2:
3145             return 1 / (x * x);
3146         default:
3147         }
3148
3149         m = -n;
3150         v = p / x;
3151     }
3152     else
3153     {
3154         switch (n)
3155         {
3156         case 0:
3157             return 1.0;
3158         case 1:
3159             return x;
3160         case 2:
3161             return x * x;
3162         default:
3163         }
3164
3165         v = x;
3166     }
3167
3168     while (1)
3169     {
3170         if (m & 1)
3171             p *= v;
3172         m >>= 1;
3173         if (!m)
3174             break;
3175         v *= v;
3176     }
3177     return p;
3178 }
3179
3180 unittest
3181 {
3182     // Make sure it instantiates and works properly on immutable values and
3183     // with various integer and float types.
3184     immutable real x = 46;
3185     immutable float xf = x;
3186     immutable double xd = x;
3187     immutable uint one = 1;
3188     immutable ushort two = 2;
3189     immutable ubyte three = 3;
3190     immutable ulong eight = 8;
3191
3192     immutable int neg1 = -1;
3193     immutable short neg2 = -2;
3194     immutable byte neg3 = -3;
3195     immutable long neg8 = -8;
3196
3197
3198     assert(pow(x,0) == 1.0);
3199     assert(pow(xd,one) == x);
3200     assert(pow(xf,two) == x * x);
3201     assert(pow(x,three) == x * x * x);
3202     assert(pow(x,eight) == (x * x) * (x * x) * (x * x) * (x * x));
3203
3204     assert(pow(x, neg1) == 1 / x);
3205     assert(pow(xd, neg2) == 1 / (x * x));
3206     assert(pow(x, neg3) == 1 / (x * x * x));
3207     assert(pow(xf, neg8) == 1 / ((x * x) * (x * x) * (x * x) * (x * x)));
3208 }
3209
3210 /** Compute the value of an integer x, raised to the power of a positive
3211  * integer n.
3212  *
3213  *  If both x and n are 0, the result is 1.
3214  *  If n is negative, an integer divide error will occur at runtime,
3215  * regardless of the value of x.
3216  */
3217
3218 typeof(Unqual!(F).init * Unqual!(G).init) pow(F, G)(F x, G n) @trusted pure nothrow
3219 if (isIntegral!(F) && isIntegral!(G))
3220 {
3221     if (n<0) return x/0; // Only support positive powers
3222     typeof(return) p, v = void;
3223     Unqual!G m = n;
3224
3225     switch (m)
3226     {
3227     case 0:
3228         p = 1;
3229         break;
3230
3231     case 1:
3232         p = x;
3233         break;
3234
3235     case 2:
3236         p = x * x;
3237         break;
3238
3239     default:
3240         v = x;
3241         p = 1;
3242         while (1){
3243             if (m & 1)
3244                 p *= v;
3245             m >>= 1;
3246             if (!m)
3247                 break;
3248             v *= v;
3249         }
3250         break;
3251     }
3252     return p;
3253 }
3254
3255 unittest
3256 {
3257     immutable int one = 1;
3258     immutable byte two = 2;
3259     immutable ubyte three = 3;
3260     immutable short four = 4;
3261     immutable long ten = 10;
3262
3263     assert(pow(two, three) == 8);
3264     assert(pow(two, ten) == 1024);
3265     assert(pow(one, ten) == 1);
3266     assert(pow(ten, four) == 10_000);
3267     assert(pow(four, 10) == 1_048_576);
3268     assert(pow(three, four) == 81);
3269
3270 }
3271
3272 /**Computes integer to floating point powers.*/
3273 real pow(I, F)(I x, F y) @trusted pure nothrow
3274     if(isIntegral!I && isFloatingPoint!F)
3275 {
3276     return pow(cast(real) x, cast(Unqual!F) y);
3277 }
3278
3279 /*********************************************
3280  * Calculates x$(SUP y).
3281  *
3282  * $(TABLE_SV
3283  * $(TR $(TH x) $(TH y) $(TH pow(x, y))
3284  *      $(TH div 0) $(TH invalid?))
3285  * $(TR $(TD anything)      $(TD $(PLUSMN)0.0)                $(TD 1.0)
3286  *      $(TD no)        $(TD no) )
3287  * $(TR $(TD |x| $(GT) 1)    $(TD +$(INFIN))                  $(TD +$(INFIN))
3288  *      $(TD no)        $(TD no) )
3289  * $(TR $(TD |x| $(LT) 1)    $(TD +$(INFIN))                  $(TD +0.0)
3290  *      $(TD no)        $(TD no) )
3291  * $(TR $(TD |x| $(GT) 1)    $(TD -$(INFIN))                  $(TD +0.0)
3292  *      $(TD no)        $(TD no) )
3293  * $(TR $(TD |x| $(LT) 1)    $(TD -$(INFIN))                  $(TD +$(INFIN))
3294  *      $(TD no)        $(TD no) )
3295  * $(TR $(TD +$(INFIN))      $(TD $(GT) 0.0)                  $(TD +$(INFIN))
3296  *      $(TD no)        $(TD no) )
3297  * $(TR $(TD +$(INFIN))      $(TD $(LT) 0.0)                  $(TD +0.0)
3298  *      $(TD no)        $(TD no) )
3299  * $(TR $(TD -$(INFIN))      $(TD odd integer $(GT) 0.0)      $(TD -$(INFIN))
3300  *      $(TD no)        $(TD no) )
3301  * $(TR $(TD -$(INFIN))      $(TD $(GT) 0.0, not odd integer) $(TD +$(INFIN))
3302  *      $(TD no)        $(TD no))
3303  * $(TR $(TD -$(INFIN))      $(TD odd integer $(LT) 0.0)      $(TD -0.0)
3304  *      $(TD no)        $(TD no) )
3305  * $(TR $(TD -$(INFIN))      $(TD $(LT) 0.0, not odd integer) $(TD +0.0)
3306  *      $(TD no)        $(TD no) )
3307  * $(TR $(TD $(PLUSMN)1.0)   $(TD $(PLUSMN)$(INFIN))          $(TD $(NAN))
3308  *      $(TD no)        $(TD yes) )
3309  * $(TR $(TD $(LT) 0.0)      $(TD finite, nonintegral)        $(TD $(NAN))
3310  *      $(TD no)        $(TD yes))
3311  * $(TR $(TD $(PLUSMN)0.0)   $(TD odd integer $(LT) 0.0)      $(TD $(PLUSMNINF))
3312  *      $(TD yes)       $(TD no) )
3313  * $(TR $(TD $(PLUSMN)0.0)   $(TD $(LT) 0.0, not odd integer) $(TD +$(INFIN))
3314  *      $(TD yes)       $(TD no))
3315  * $(TR $(TD $(PLUSMN)0.0)   $(TD odd integer $(GT) 0.0)      $(TD $(PLUSMN)0.0)
3316  *      $(TD no)        $(TD no) )
3317  * $(TR $(TD $(PLUSMN)0.0)   $(TD $(GT) 0.0, not odd integer) $(TD +0.0)
3318  *      $(TD no)        $(TD no) )
3319  * )
3320  */
3321
3322 Unqual!(Largest!(F, G)) pow(F, G)(F x, G y) @trusted pure nothrow
3323     if (isFloatingPoint!(F) && isFloatingPoint!(G))
3324 {
3325     alias typeof(return) Float;
3326
3327     static real impl(real x, real y) pure nothrow
3328     {
3329         if (isNaN(y))
3330             return y;
3331
3332         if (y == 0)
3333             return 1;           // even if x is $(NAN)
3334         if (isNaN(x) && y != 0)
3335             return x;
3336         if (isInfinity(y))
3337         {
3338             if (fabs(x) > 1)
3339             {
3340                 if (signbit(y))
3341                     return +0.0;
3342                 else
3343                     return F.infinity;
3344             }
3345             else if (fabs(x) == 1)
3346             {
3347                 return y * 0; // generate NaN.
3348             }
3349             else // < 1
3350             {
3351                 if (signbit(y))
3352                     return F.infinity;
3353                 else
3354                     return +0.0;
3355             }
3356         }
3357         if (isInfinity(x))
3358         {
3359             if (signbit(x))
3360             {   long i;
3361
3362                 i = cast(long)y;
3363                 if (y > 0)
3364                 {
3365                     if (i == y && i & 1)
3366                         return -F.infinity;
3367                     else
3368                         return F.infinity;
3369                 }
3370                 else if (y < 0)
3371                 {
3372                     if (i == y && i & 1)
3373                         return -0.0;
3374                     else
3375                         return +0.0;
3376                 }
3377             }
3378             else
3379             {
3380                 if (y > 0)
3381                     return F.infinity;
3382                 else if (y < 0)
3383                     return +0.0;
3384             }
3385         }
3386
3387         if (x == 0.0)
3388         {
3389             if (signbit(x))
3390             {   long i;
3391
3392                 i = cast(long)y;
3393                 if (y > 0)
3394                 {
3395                     if (i == y && i & 1)
3396                         return -0.0;
3397                     else
3398                         return +0.0;
3399                 }
3400                 else if (y < 0)
3401                 {
3402                     if (i == y && i & 1)
3403                         return -F.infinity;
3404                     else
3405                         return F.infinity;
3406                 }
3407             }
3408             else
3409             {
3410                 if (y > 0)
3411                     return +0.0;
3412                 else if (y < 0)
3413                     return F.infinity;
3414             }
3415         }
3416         double sign = 1.0;
3417         if (x < 0) {
3418             // Result is real only if y is an integer
3419             // Check for a non-zero fractional part
3420             if (y > -1.0 / real.epsilon && y < 1.0 / real.epsilon)
3421             {
3422                 long w = cast(long)y;
3423                 if (w != y)
3424                     return sqrt(x); // Complex result -- create a NaN
3425                 if (w & 1) sign = -1.0;
3426             }
3427             x = -x;
3428         }
3429         version(INLINE_YL2X) {       
3430             // If x > 0, x ^^ y == 2 ^^ ( y * log2(x) )
3431             // TODO: This is not accurate in practice. A fast and accurate
3432             // (though complicated) method is described in:
3433             // "An efficient rounding boundary test for pow(x, y)
3434             // in double precision", C.Q. Lauter and V. Lefèvre, INRIA (2007).
3435             return sign * exp2( yl2x(x, y) );
3436         } else {
3437             return sign * core.stdc.math.powl(x, y);
3438         }
3439     }
3440     return impl(x, y);
3441 }
3442
3443 unittest
3444 {
3445     // Test all the special values.  These unittests can be run on Windows
3446     // by temporarily changing the version(linux) to version(all).
3447     immutable float zero = 0;
3448     immutable real one = 1;
3449     immutable double two = 2;
3450     immutable float three = 3;
3451     immutable float fnan = float.nan;
3452     immutable double dnan = double.nan;
3453     immutable real rnan = real.nan;
3454     immutable dinf = double.infinity;
3455     immutable rninf = -real.infinity;
3456
3457     assert(pow(fnan, zero) == 1);
3458     assert(pow(dnan, zero) == 1);
3459     assert(pow(rnan, zero) == 1);
3460
3461     assert(pow(two, dinf) == double.infinity);
3462     assert(isIdentical(pow(0.2f, dinf), +0.0));
3463     assert(pow(0.99999999L, rninf) == real.infinity);
3464     assert(isIdentical(pow(1.000000001, rninf), +0.0));
3465     assert(pow(dinf, 0.001) == dinf);
3466     assert(isIdentical(pow(dinf, -0.001), +0.0));
3467     assert(pow(rninf, 3.0L) == rninf);
3468     assert(pow(rninf, 2.0L) == real.infinity);
3469     assert(isIdentical(pow(rninf, -3.0), -0.0));
3470     assert(isIdentical(pow(rninf, -2.0), +0.0));
3471
3472     // @@@BUG@@@ somewhere
3473     version(OSX) {} else assert(isNaN(pow(one, dinf)));
3474     version(OSX) {} else assert(isNaN(pow(-one, dinf)));
3475     assert(isNaN(pow(-0.2, PI)));
3476     // boundary cases. Note that epsilon == 2^^-n for some n,
3477     // so 1/epsilon == 2^^n is always even.
3478     assert(pow(-1.0L, 1/real.epsilon - 1.0L) == -1.0L);
3479     assert(pow(-1.0L, 1/real.epsilon) == 1.0L);
3480     assert(isNaN(pow(-1.0L, 1/real.epsilon-0.5L)));
3481     assert(isNaN(pow(-1.0L, -1/real.epsilon+0.5L)));
3482    
3483     assert(pow(0.0, -3.0) == double.infinity);
3484     assert(pow(-0.0, -3.0) == -double.infinity);
3485     assert(pow(0.0, -PI) == double.infinity);
3486     assert(pow(-0.0, -PI) == double.infinity);
3487     assert(isIdentical(pow(0.0, 5.0), 0.0));
3488     assert(isIdentical(pow(-0.0, 5.0), -0.0));
3489     assert(isIdentical(pow(0.0, 6.0), 0.0));
3490     assert(isIdentical(pow(-0.0, 6.0), 0.0));
3491
3492     // Now, actual numbers.
3493     assert(approxEqual(pow(two, three), 8.0));
3494     assert(approxEqual(pow(two, -2.5), 0.1767767));
3495
3496     // Test integer to float power.
3497     immutable uint twoI = 2;
3498     assert(approxEqual(pow(twoI, three), 8.0));
3499 }
3500
3501 /**************************************
3502  * To what precision is x equal to y?
3503  *
3504  * Returns: the number of mantissa bits which are equal in x and y.
3505  * eg, 0x1.F8p+60 and 0x1.F1p+60 are equal to 5 bits of precision.
3506  *
3507  *      $(TABLE_SV
3508  *      $(TR $(TH x)      $(TH y)          $(TH feqrel(x, y)))
3509  *      $(TR $(TD x)      $(TD x)          $(TD real.mant_dig))
3510  *      $(TR $(TD x)      $(TD $(GT)= 2*x) $(TD 0))
3511  *      $(TR $(TD x)      $(TD $(LT)= x/2) $(TD 0))
3512  *      $(TR $(TD $(NAN)) $(TD any)        $(TD 0))
3513  *      $(TR $(TD any)    $(TD $(NAN))     $(TD 0))
3514  *      )
3515  */
3516 int feqrel(X)(X x, X y) @trusted pure nothrow
3517     if (isFloatingPoint!(X))
3518 {
3519     /* Public Domain. Author: Don Clugston, 18 Aug 2005.
3520      */
3521     static if (X.mant_dig == 106) { // doubledouble.
3522         if (cast(double*)(&x)[MANTISSA_MSB] == cast(double*)(&y)[MANTISSA_MSB]) {
3523             return double.mant_dig
3524             + feqrel(cast(double*)(&x)[MANTISSA_LSB],
3525                     cast(double*)(&y)[MANTISSA_LSB]);
3526         } else {
3527             return feqrel(cast(double*)(&x)[MANTISSA_MSB],
3528                     cast(double*)(&y)[MANTISSA_MSB]);
3529         }
3530     } else static if (X.mant_dig==64 || X.mant_dig==113 || X.mant_dig==53) {
3531
3532         if (x == y) return X.mant_dig; // ensure diff!=0, cope with INF.
3533
3534         X diff = fabs(x - y);
3535
3536         ushort *pa = cast(ushort *)(&x);
3537         ushort *pb = cast(ushort *)(&y);
3538         ushort *pd = cast(ushort *)(&diff);
3539
3540         alias floatTraits!(X) F;
3541
3542         // The difference in abs(exponent) between x or y and abs(x-y)
3543         // is equal to the number of significand bits of x which are
3544         // equal to y. If negative, x and y have different exponents.
3545         // If positive, x and y are equal to 'bitsdiff' bits.
3546         // AND with 0x7FFF to form the absolute value.
3547         // To avoid out-by-1 errors, we subtract 1 so it rounds down
3548         // if the exponents were different. This means 'bitsdiff' is
3549         // always 1 lower than we want, except that if bitsdiff==0,
3550         // they could have 0 or 1 bits in common.
3551
3552         static if (X.mant_dig==64 || X.mant_dig==113) { // real80 or quadruple
3553             int bitsdiff = ( ((pa[F.EXPPOS_SHORT] & F.EXPMASK)
3554                               + (pb[F.EXPPOS_SHORT] & F.EXPMASK) - 1) >> 1)
3555                               - pd[F.EXPPOS_SHORT];
3556         } else static if (X.mant_dig==53) { // double
3557             int bitsdiff = (( ((pa[F.EXPPOS_SHORT]&0x7FF0)
3558                                + (pb[F.EXPPOS_SHORT]&0x7FF0)-0x10)>>1)
3559                                - (pd[F.EXPPOS_SHORT]&0x7FF0))>>4;
3560         }
3561         if (pd[F.EXPPOS_SHORT] == 0)
3562         {   // Difference is denormal
3563             // For denormals, we need to add the number of zeros that
3564             // lie at the start of diff's significand.
3565             // We do this by multiplying by 2^^real.mant_dig
3566             diff *= F.RECIP_EPSILON;
3567             return bitsdiff + X.mant_dig - pd[F.EXPPOS_SHORT];
3568         }
3569
3570         if (bitsdiff > 0)
3571             return bitsdiff + 1; // add the 1 we subtracted before
3572
3573         // Avoid out-by-1 errors when factor is almost 2.
3574         static if (X.mant_dig==64 || X.mant_dig==113) { // real80 or quadruple
3575             return (bitsdiff == 0) ? (pa[F.EXPPOS_SHORT] == pb[F.EXPPOS_SHORT]) : 0;
3576         } else static if (X.mant_dig==53) { // double
3577             if (bitsdiff == 0
3578                 && !((pa[F.EXPPOS_SHORT] ^ pb[F.EXPPOS_SHORT])& F.EXPMASK)) {
3579                 return 1;
3580             } else return 0;
3581         }
3582     }
3583 }
3584
3585 unittest
3586 {
3587    // Exact equality
3588    assert(feqrel(real.max,real.max)==real.mant_dig);
3589    assert(feqrel(0.0L,0.0L)==real.mant_dig);
3590    assert(feqrel(7.1824L,7.1824L)==real.mant_dig);
3591    assert(feqrel(real.infinity,real.infinity)==real.mant_dig);
3592
3593    // a few bits away from exact equality
3594    real w=1;
3595    for (int i=1; i<real.mant_dig-1; ++i) {
3596       assert(feqrel(1+w*real.epsilon,1.0L)==real.mant_dig-i);
3597       assert(feqrel(1-w*real.epsilon,1.0L)==real.mant_dig-i);
3598       assert(feqrel(1.0L,1+(w-1)*real.epsilon)==real.mant_dig-i+1);
3599       w*=2;
3600    }
3601    assert(feqrel(1.5+real.epsilon,1.5L)==real.mant_dig-1);
3602    assert(feqrel(1.5-real.epsilon,1.5L)==real.mant_dig-1);
3603    assert(feqrel(1.5-real.epsilon,1.5+real.epsilon)==real.mant_dig-2);
3604
3605    assert(feqrel(real.min_normal/8,real.min_normal/17)==3);;
3606
3607    // Numbers that are close
3608    assert(feqrel(0x1.Bp+84, 0x1.B8p+84)==5);
3609    assert(feqrel(0x1.8p+10, 0x1.Cp+10)==2);
3610    assert(feqrel(1.5*(1-real.epsilon), 1.0L)==2);
3611    assert(feqrel(1.5, 1.0)==1);
3612    assert(feqrel(2*(1-real.epsilon), 1.0L)==1);
3613
3614    // Factors of 2
3615    assert(feqrel(real.max,real.infinity)==0);
3616    assert(feqrel(2*(1-real.epsilon), 1.0L)==1);
3617    assert(feqrel(1.0, 2.0)==0);
3618    assert(feqrel(4.0, 1.0)==0);
3619
3620    // Extreme inequality
3621    assert(feqrel(real.nan,real.nan)==0);
3622    assert(feqrel(0.0L,-real.nan)==0);
3623    assert(feqrel(real.nan,real.infinity)==0);
3624    assert(feqrel(real.infinity,-real.infinity)==0);
3625    assert(feqrel(-real.max,real.infinity)==0);
3626    assert(feqrel(real.max,-real.max)==0);
3627 }
3628
3629 package: // Not public yet
3630 /* Return the value that lies halfway between x and y on the IEEE number line.
3631  *
3632  * Formally, the result is the arithmetic mean of the binary significands of x
3633  * and y, multiplied by the geometric mean of the binary exponents of x and y.
3634  * x and y must have the same sign, and must not be NaN.
3635  * Note: this function is useful for ensuring O(log n) behaviour in algorithms
3636  * involving a 'binary chop'.
3637  *
3638  * Special cases:
3639  * If x and y are within a factor of 2, (ie, feqrel(x, y) > 0), the return value
3640  * is the arithmetic mean (x + y) / 2.
3641  * If x and y are even powers of 2, the return value is the geometric mean,
3642  *   ieeeMean(x, y) = sqrt(x * y).
3643  *
3644  */
3645 T ieeeMean(T)(T x, T y)  @trusted pure nothrow
3646 in {
3647     // both x and y must have the same sign, and must not be NaN.
3648     assert(signbit(x) == signbit(y));
3649     assert(x<>=0 && y<>=0);
3650 }
3651 body {
3652     // Runtime behaviour for contract violation:
3653     // If signs are opposite, or one is a NaN, return 0.
3654     if (!((x>=0 && y>=0) || (x<=0 && y<=0))) return 0.0;
3655
3656     // The implementation is simple: cast x and y to integers,
3657     // average them (avoiding overflow), and cast the result back to a floating-point number.
3658
3659     alias floatTraits!(real) F;
3660     T u;
3661     static if (T.mant_dig==64) { // real80
3662         // There's slight additional complexity because they are actually
3663         // 79-bit reals...
3664         ushort *ue = cast(ushort *)&u;
3665         ulong *ul = cast(ulong *)&u;
3666         ushort *xe = cast(ushort *)&x;
3667         ulong *xl = cast(ulong *)&x;
3668         ushort *ye = cast(ushort *)&y;
3669         ulong *yl = cast(ulong *)&y;
3670         // Ignore the useless implicit bit. (Bonus: this prevents overflows)
3671         ulong m = ((*xl) & 0x7FFF_FFFF_FFFF_FFFFL) + ((*yl) & 0x7FFF_FFFF_FFFF_FFFFL);
3672
3673         // @@@ BUG? @@@
3674         // Cast shouldn't be here
3675         ushort e = cast(ushort) ((xe[F.EXPPOS_SHORT] & F.EXPMASK)
3676                                  + (ye[F.EXPPOS_SHORT] & F.EXPMASK));
3677         if (m & 0x8000_0000_0000_0000L) {
3678             ++e;
3679             m &= 0x7FFF_FFFF_FFFF_FFFFL;
3680         }
3681         // Now do a multi-byte right shift
3682         uint c = e & 1; // carry
3683         e >>= 1;
3684     m >>>= 1;
3685     if (c) m |= 0x4000_0000_0000_0000L; // shift carry into significand
3686     if (e) *ul = m | 0x8000_0000_0000_0000L; // set implicit bit...
3687     else *ul = m; // ... unless exponent is 0 (denormal or zero).
3688     ue[4]= e | (xe[F.EXPPOS_SHORT]& 0x8000); // restore sign bit
3689     } else static if(T.mant_dig == 113) { //quadruple
3690         // This would be trivial if 'ucent' were implemented...
3691         ulong *ul = cast(ulong *)&u;
3692         ulong *xl = cast(ulong *)&x;
3693         ulong *yl = cast(ulong *)&y;
3694         // Multi-byte add, then multi-byte right shift.
3695         ulong mh = ((xl[MANTISSA_MSB] & 0x7FFF_FFFF_FFFF_FFFFL)
3696                     + (yl[MANTISSA_MSB] & 0x7FFF_FFFF_FFFF_FFFFL));
3697         // Discard the lowest bit (to avoid overflow)
3698         ulong ml = (xl[MANTISSA_LSB]>>>1) + (yl[MANTISSA_LSB]>>>1);
3699         // add the lowest bit back in, if necessary.
3700         if (xl[MANTISSA_LSB] & yl[MANTISSA_LSB] & 1) {
3701             ++ml;
3702             if (ml==0) ++mh;
3703         }
3704         mh >>>=1;
3705         ul[MANTISSA_MSB] = mh | (xl[MANTISSA_MSB] & 0x8000_0000_0000_0000);
3706         ul[MANTISSA_LSB] = ml;
3707     } else static if (T.mant_dig == double.mant_dig) {
3708         ulong *ul = cast(ulong *)&u;
3709         ulong *xl = cast(ulong *)&x;
3710         ulong *yl = cast(ulong *)&y;
3711         ulong m = (((*xl) & 0x7FFF_FFFF_FFFF_FFFFL)
3712                    + ((*yl) & 0x7FFF_FFFF_FFFF_FFFFL)) >>> 1;
3713                    m |= ((*xl) & 0x8000_0000_0000_0000L);
3714                    *ul = m;
3715     } else static if (T.mant_dig == float.mant_dig) {
3716         uint *ul = cast(uint *)&u;
3717         uint *xl = cast(uint *)&x;
3718         uint *yl = cast(uint *)&y;
3719         uint m = (((*xl) & 0x7FFF_FFFF) + ((*yl) & 0x7FFF_FFFF)) >>> 1;
3720         m |= ((*xl) & 0x8000_0000);
3721         *ul = m;
3722     } else {
3723         assert(0, "Not implemented");
3724     }
3725     return u;
3726 }
3727
3728 unittest {
3729     assert(ieeeMean(-0.0,-1e-20)<0);
3730     assert(ieeeMean(0.0,1e-20)>0);
3731
3732     assert(ieeeMean(1.0L,4.0L)==2L);
3733     assert(ieeeMean(2.0*1.013,8.0*1.013)==4*1.013);
3734     assert(ieeeMean(-1.0L,-4.0L)==-2L);
3735     assert(ieeeMean(-1.0,-4.0)==-2);
3736     assert(ieeeMean(-1.0f,-4.0f)==-2f);
3737     assert(ieeeMean(-1.0,-2.0)==-1.5);
3738     assert(ieeeMean(-1*(1+8*real.epsilon),-2*(1+8*real.epsilon))
3739                  ==-1.5*(1+5*real.epsilon));
3740     assert(ieeeMean(0x1p60,0x1p-10)==0x1p25);
3741     static if (real.mant_dig==64) { // x87, 80-bit reals
3742       assert(ieeeMean(1.0L,real.infinity)==0x1p8192L);
3743       assert(ieeeMean(0.0L,real.infinity)==1.5);
3744     }
3745     assert(ieeeMean(0.5*real.min_normal*(1-4*real.epsilon),0.5*real.min_normal)
3746            == 0.5*real.min_normal*(1-2*real.epsilon));
3747 }
3748
3749 public:
3750
3751
3752 /***********************************
3753  * Evaluate polynomial A(x) = $(SUB a, 0) + $(SUB a, 1)x + $(SUB a, 2)$(POWER x,2)
3754  *                          + $(SUB a,3)$(POWER x,3); ...
3755  *
3756  * Uses Horner's rule A(x) = $(SUB a, 0) + x($(SUB a, 1) + x($(SUB a, 2)
3757  *                         + x($(SUB a, 3) + ...)))
3758  * Params:
3759  *      A =     array of coefficients $(SUB a, 0), $(SUB a, 1), etc.
3760  */
3761 real poly(real x, const real[] A) @trusted pure nothrow
3762 in
3763 {
3764     assert(A.length > 0);
3765 }
3766 body
3767 {
3768     version (D_InlineAsm_X86)
3769     {
3770         version (Windows)
3771         {
3772         // BUG: This code assumes a frame pointer in EBP.
3773             asm // assembler by W. Bright
3774             {
3775                 // EDX = (A.length - 1) * real.sizeof
3776                 mov     ECX,A[EBP]              ; // ECX = A.length
3777                 dec     ECX                     ;
3778                 lea     EDX,[ECX][ECX*8]        ;
3779                 add     EDX,ECX                 ;
3780                 add     EDX,A+4[EBP]            ;
3781                 fld     real ptr [EDX]          ; // ST0 = coeff[ECX]
3782                 jecxz   return_ST               ;
3783                 fld     x[EBP]                  ; // ST0 = x
3784                 fxch    ST(1)                   ; // ST1 = x, ST0 = r
3785                 align   4                       ;
3786         L2:     fmul    ST,ST(1)                ; // r *= x
3787                 fld     real ptr -10[EDX]       ;
3788                 sub     EDX,10                  ; // deg--
3789                 faddp   ST(1),ST                ;
3790                 dec     ECX                     ;
3791                 jne     L2                      ;
3792                 fxch    ST(1)                   ; // ST1 = r, ST0 = x
3793                 fstp    ST(0)                   ; // dump x
3794                 align   4                       ;
3795         return_ST:                              ;
3796                 ;
3797             }
3798         }
3799         else version (linux)
3800         {
3801             asm // assembler by W. Bright
3802             {
3803                 // EDX = (A.length - 1) * real.sizeof
3804                 mov     ECX,A[EBP]              ; // ECX = A.length
3805                 dec     ECX                     ;
3806                 lea     EDX,[ECX*8]             ;
3807                 lea     EDX,[EDX][ECX*4]        ;
3808                 add     EDX,A+4[EBP]            ;
3809                 fld     real ptr [EDX]          ; // ST0 = coeff[ECX]
3810                 jecxz   return_ST               ;
3811                 fld     x[EBP]                  ; // ST0 = x
3812                 fxch    ST(1)                   ; // ST1 = x, ST0 = r
3813                 align   4                       ;
3814         L2:     fmul    ST,ST(1)                ; // r *= x
3815                 fld     real ptr -12[EDX]       ;
3816                 sub     EDX,12                  ; // deg--
3817                 faddp   ST(1),ST                ;
3818                 dec     ECX                     ;
3819                 jne     L2                      ;
3820                 fxch    ST(1)                   ; // ST1 = r, ST0 = x
3821                 fstp    ST(0)                   ; // dump x
3822                 align   4                       ;
3823         return_ST:                              ;
3824                 ;
3825             }
3826         }
3827         else version (OSX)
3828         {
3829             asm // assembler by W. Bright
3830             {
3831                 // EDX = (A.length - 1) * real.sizeof
3832                 mov     ECX,A[EBP]              ; // ECX = A.length
3833                 dec     ECX                     ;
3834                 lea     EDX,[ECX*8]             ;
3835                 add     EDX,EDX                 ;
3836                 add     EDX,A+4[EBP]            ;
3837                 fld     real ptr [EDX]          ; // ST0 = coeff[ECX]
3838                 jecxz   return_ST               ;
3839                 fld     x[EBP]                  ; // ST0 = x
3840                 fxch    ST(1)                   ; // ST1 = x, ST0 = r
3841                 align   4                       ;
3842         L2:     fmul    ST,ST(1)                ; // r *= x
3843                 fld     real ptr -16[EDX]       ;
3844                 sub     EDX,16                  ; // deg--
3845                 faddp   ST(1),ST                ;
3846                 dec     ECX                     ;
3847                 jne     L2                      ;
3848                 fxch    ST(1)                   ; // ST1 = r, ST0 = x
3849                 fstp    ST(0)                   ; // dump x
3850                 align   4                       ;
3851         return_ST:                              ;
3852                 ;
3853             }
3854         }
3855         else version (FreeBSD)
3856         {
3857             asm // assembler by W. Bright
3858             {
3859                 // EDX = (A.length - 1) * real.sizeof
3860                 mov     ECX,A[EBP]              ; // ECX = A.length
3861                 dec     ECX                     ;
3862                 lea     EDX,[ECX*8]             ;
3863                 lea     EDX,[EDX][ECX*4]        ;
3864                 add     EDX,A+4[EBP]            ;
3865                 fld     real ptr [EDX]          ; // ST0 = coeff[ECX]
3866                 jecxz   return_ST               ;
3867                 fld     x[EBP]                  ; // ST0 = x
3868                 fxch    ST(1)                   ; // ST1 = x, ST0 = r
3869                 align   4                       ;
3870         L2:     fmul    ST,ST(1)                ; // r *= x
3871                 fld     real ptr -12[EDX]       ;
3872                 sub     EDX,12                  ; // deg--
3873                 faddp   ST(1),ST                ;
3874                 dec     ECX                     ;
3875                 jne     L2                      ;
3876                 fxch    ST(1)                   ; // ST1 = r, ST0 = x
3877                 fstp    ST(0)                   ; // dump x
3878                 align   4                       ;
3879         return_ST:                              ;
3880                 ;
3881             }
3882         }
3883         else
3884         {
3885             static assert(0);
3886         }
3887     }
3888     else
3889     {
3890         sizediff_t i = A.length - 1;
3891         real r = A[i];
3892         while (--i >= 0)
3893         {
3894             r *= x;
3895             r += A[i];
3896         }
3897         return r;
3898     }
3899 }
3900
3901 unittest
3902 {
3903     debug (math) printf("math.poly.unittest\n");
3904     real x = 3.1;
3905     static real pp[] = [56.1, 32.7, 6];
3906
3907     assert( poly(x, pp) == (56.1L + (32.7L + 6L * x) * x) );
3908 }
3909
3910 /**
3911    Computes whether $(D lhs) is approximately equal to $(D rhs)
3912    admitting a maximum relative difference $(D maxRelDiff) and a
3913    maximum absolute difference $(D maxAbsDiff).
3914
3915    If the two inputs are ranges, $(D approxEqual) returns true if and
3916    only if the ranges have the same number of elements and if $(D
3917    approxEqual) evaluates to $(D true) for each pair of elements.
3918  */
3919 bool approxEqual(T, U, V)(T lhs, U rhs, V maxRelDiff, V maxAbsDiff = 1e-5)
3920 {
3921     static if (isInputRange!T)
3922     {
3923         static if (isInputRange!U)
3924         {
3925             // Two ranges
3926             for (;; lhs.popFront, rhs.popFront)
3927             {
3928                 if (lhs.empty) return rhs.empty;
3929                 if (rhs.empty) return lhs.empty;
3930                 if (!approxEqual(lhs.front, rhs.front, maxRelDiff, maxAbsDiff))
3931                     return false;
3932             }
3933         }
3934         else
3935         {
3936             // lhs is range, rhs is number
3937             for (; !lhs.empty; lhs.popFront)
3938             {
3939                 if (!approxEqual(lhs.front, rhs, maxRelDiff, maxAbsDiff))
3940                     return false;
3941             }
3942             return true;
3943         }
3944     }
3945     else
3946     {
3947         static if (isInputRange!U)
3948         {
3949             // lhs is number, rhs is array
3950             return approxEqual(rhs, lhs, maxRelDiff, maxAbsDiff);
3951         }
3952         else
3953         {
3954             // two numbers
3955             //static assert(is(T : real) && is(U : real));
3956             if (rhs == 0)
3957             {
3958                 return fabs(lhs) <= maxAbsDiff;
3959             }
3960             static if (is(typeof(lhs.infinity)) && is(typeof(rhs.infinity)))
3961             {
3962                 if (lhs == lhs.infinity && rhs == rhs.infinity ||
3963                     lhs == -lhs.infinity && rhs == -rhs.infinity) return true;
3964             }
3965             return fabs((lhs - rhs) / rhs) <= maxRelDiff
3966                 || maxAbsDiff != 0 && fabs(lhs - rhs) <= maxAbsDiff;
3967         }
3968     }
3969 }
3970
3971 /**
3972    Returns $(D approxEqual(lhs, rhs, 1e-2, 1e-5)).
3973  */
3974 bool approxEqual(T, U)(T lhs, U rhs)
3975 {
3976     return approxEqual(lhs, rhs, 1e-2, 1e-5);
3977 }
3978
3979 unittest
3980 {
3981     assert(approxEqual(1.0, 1.0099));
3982     assert(!approxEqual(1.0, 1.011));
3983     float[] arr1 = [ 1.0, 2.0, 3.0 ];
3984     double[] arr2 = [ 1.001, 1.999, 3 ];
3985     assert(approxEqual(arr1, arr2));
3986
3987     real num = real.infinity;
3988     assert(num == real.infinity);  // Passes.
3989     assert(approxEqual(num, real.infinity));  // Fails.
3990     num = -real.infinity;
3991     assert(num == -real.infinity);  // Passes.
3992     assert(approxEqual(num, -real.infinity));  // Fails.
3993 }
3994
3995 // Included for backwards compatibility with Phobos1
3996 alias isNaN isnan;
3997 alias isFinite isfinite;
3998 alias isNormal isnormal;
3999 alias isSubnormal issubnormal;
4000 alias isInfinity isinf;
4001
4002 /* **********************************
4003  * Building block functions, they
4004  * translate to a single x87 instruction.
4005  */
4006
4007 real yl2x(real x, real y)   @safe pure nothrow;       // y * log2(x)
4008 real yl2xp1(real x, real y) @safe pure nothrow;       // y * log2(x + 1)
4009
4010 unittest
4011 {
4012     version (INLINE_YL2X)
4013     {
4014         assert(yl2x(1024, 1) == 10);
4015         assert(yl2xp1(1023, 1) == 10);
4016     }
4017 }
4018
4019 unittest
4020 {
4021     real num = real.infinity;
4022     assert(num == real.infinity);  // Passes.
4023     assert(approxEqual(num, real.infinity));  // Fails.
4024 }
Note: See TracBrowser for help on using the browser.