1 /**
2  * Elementary Mathematical Functions
3  *
4  * Copyright:
5  *     Portions Copyright (C) 2001-2005 Digital Mars.
6  *     Some parts copyright (c) 2009-2016 dunnhumby Germany GmbH.
7  *     All rights reserved.
8  *
9  * License:
10  *     Tango Dual License: 3-Clause BSD License / Academic Free License v3.0.
11  *     See LICENSE_TANGO.txt for details.
12  *
13  * Authors: Walter Bright, Don Clugston, Sean Kelly
14  *
15  */
16 
17 /**
18  * Macros:
19  *  NAN = $(RED NAN)
20  *  TEXTNAN = $(RED NAN:$1 )
21  *  SUP = <span style="vertical-align:super;font-size:smaller">$0</span>
22  *  GAMMA =  &#915;
23  *  INTEGRAL = &#8747;
24  *  INTEGRATE = $(BIG &#8747;<sub>$(SMALL $1)</sub><sup>$2</sup>)
25  *  POWER = $1<sup>$2</sup>
26  *  BIGSUM = $(BIG &Sigma; <sup>$2</sup><sub>$(SMALL $1)</sub>)
27  *  CHOOSE = $(BIG &#40;) <sup>$(SMALL $1)</sup><sub>$(SMALL $2)</sub> $(BIG &#41;)
28  *  PLUSMN = &plusmn;
29  *  INFIN = &infin;
30  *  PLUSMNINF = &plusmn;&infin;
31  *  PI = &pi;
32  *  LT = &lt;
33  *  GT = &gt;
34  *  SQRT = &radix;
35  *  HALF = &frac12;
36  *  TABLE_SV = <table border=1 cellpadding=4 cellspacing=0>
37  *      <caption>Special Values</caption>
38  *      $0</table>
39  *  SVH = $(TR $(TH $1) $(TH $2))
40  *  SV  = $(TR $(TD $1) $(TD $2))
41  *  TABLE_DOMRG = <table border=1 cellpadding=4 cellspacing=0>$0</table>
42  *  DOMAIN = $(TR $(TD Domain) $(TD $0))
43  *  RANGE  = $(TR $(TD Range) $(TD $0))
44  */
45 
46 module ocean.math.Math;
47 
48 import ocean.meta.types.Qualifiers;
49 
50 static import core.stdc.math;
51 import ocean.math.IEEE;
52 import ocean.core.Verify;
53 
54 version (unittest) import ocean.core.Test;
55 
56 version(TangoNoAsm) {
57 
58 } else version(D_InlineAsm_X86) {
59     version = Naked_D_InlineAsm_X86;
60 } else version(D_InlineAsm_X86_64) {
61     version = Naked_D_InlineAsm_X86_64;
62 }
63 
64 /*
65  * Constants
66  */
67 
68 static immutable real E =          2.7182818284590452354L;  /** e */ // 3.32193 fldl2t 0x1.5BF0A8B1_45769535_5FF5p+1L
69 static immutable real LOG2T =      0x1.a934f0979a3715fcp+1; /** $(SUB log, 2)10 */ // 1.4427 fldl2e
70 static immutable real LOG2E =      0x1.71547652b82fe178p+0; /** $(SUB log, 2)e */ // 0.30103 fldlg2
71 static immutable real LOG2 =       0x1.34413509f79fef32p-2; /** $(SUB log, 10)2 */
72 static immutable real LOG10E =     0.43429448190325182765;  /** $(SUB log, 10)e */
73 static immutable real LN2 =        0x1.62e42fefa39ef358p-1; /** ln 2 */  // 0.693147 fldln2
74 static immutable real LN10 =       2.30258509299404568402;  /** ln 10 */
75 static immutable real PI =         0x1.921fb54442d1846ap+1; /** $(_PI) */ // 3.14159 fldpi
76 static immutable real PI_2 =       1.57079632679489661923;  /** $(PI) / 2 */
77 static immutable real PI_4 =       0.78539816339744830962;  /** $(PI) / 4 */
78 static immutable real M_1_PI =     0.31830988618379067154;  /** 1 / $(PI) */
79 static immutable real M_2_PI =     0.63661977236758134308;  /** 2 / $(PI) */
80 static immutable real M_2_SQRTPI = 1.12837916709551257390;  /** 2 / $(SQRT)$(PI) */
81 static immutable real SQRT2 =      1.41421356237309504880;  /** $(SQRT)2 */
82 static immutable real SQRT1_2 =    0.70710678118654752440;  /** $(SQRT)$(HALF) */
83 
84 //const real SQRTPI  = 1.77245385090551602729816748334114518279754945612238L; /** &radic;&pi; */
85 //const real SQRT2PI = 2.50662827463100050242E0L; /** &radic;(2 &pi;) */
86 //const real SQRTE   = 1.64872127070012814684865078781416357L; /** &radic;(e) */
87 
88 static immutable real MAXLOG = 0x1.62e42fefa39ef358p+13L;  /** log(real.max) */
89 static immutable real MINLOG = -0x1.6436716d5406e6d8p+13L; /** log(real.min*real.epsilon) */
90 static immutable real EULERGAMMA = 0.57721_56649_01532_86060_65120_90082_40243_10421_59335_93992L; /** Euler-Mascheroni constant 0.57721566.. */
91 
92 /*
93  * Primitives
94  */
95 
96 /**
97  * Calculates the absolute value
98  *
99  * For complex numbers, abs(z) = sqrt( $(POWER z.re, 2) + $(POWER z.im, 2) )
100  * = hypot(z.re, z.im).
101  */
102 real abs(real x)
103 {
104     return ocean.math.IEEE.fabs(x);
105 }
106 
107 /** ditto */
108 long abs(long x)
109 {
110     return x>=0 ? x : -x;
111 }
112 
113 /** ditto */
114 int abs(int x)
115 {
116     return x>=0 ? x : -x;
117 }
118 
119 /** ditto */
120 real abs(creal z)
121 {
122     return hypot(z.re, z.im);
123 }
124 
125 /** ditto */
126 real abs(ireal y)
127 {
128     return ocean.math.IEEE.fabs(y.im);
129 }
130 
131 unittest
132 {
133     test(isIdentical(0.0L,abs(-0.0L)));
134     test(isNaN(abs(real.nan)));
135     test(abs(-real.infinity) == real.infinity);
136     test(abs(-3.2Li) == 3.2L);
137     test(abs(71.6Li) == 71.6L);
138     test(abs(-56) == 56);
139     test(abs(2321312L)  == 2321312L);
140     test(abs(-1.0L+1.0Li) == sqrt(2.0L));
141 }
142 
143 /**
144  * Complex conjugate
145  *
146  *  conj(x + iy) = x - iy
147  *
148  * Note that z * conj(z) = $(POWER z.re, 2) + $(POWER z.im, 2)
149  * is always a real number
150  */
151 creal conj(creal z)
152 {
153     return z.re - z.im*1i;
154 }
155 
156 /** ditto */
157 ireal conj(ireal y)
158 {
159     return -y;
160 }
161 
162 unittest
163 {
164     test(conj(7 + 3i) == 7-3i);
165     ireal z = -3.2Li;
166     test(conj(z) == -z);
167 }
168 
169 private {
170     // Return the type which would be returned by a max or min operation
171 template minmaxtype(T...){
172     static if(T.length == 1) alias T[0] minmaxtype;
173     else static if(T.length > 2)
174         alias minmaxtype!(minmaxtype!(T[0..2]), T[2..$]) minmaxtype;
175     else alias typeof (T[1].init > T[0].init ? T[1].init : T[0].init) minmaxtype;
176 }
177 }
178 
179 unittest
180 {
181     static assert (is(minmaxtype!(int, long) == long));
182 }
183 
184 /** Return the minimum of the supplied arguments.
185  *
186  * Note: If the arguments are floating-point numbers, and at least one is a NaN,
187  * the result is undefined.
188  */
189 minmaxtype!(T) min(T...)(T arg){
190     static if(arg.length == 1) return arg[0];
191     else static if(arg.length == 2) return arg[1] < arg[0] ? arg[1] : arg[0];
192     static if(arg.length > 2) return min(arg[1] < arg[0] ? arg[1] : arg[0], arg[2..$]);
193 }
194 
195 /** Return the maximum of the supplied arguments.
196  *
197  * Note: If the arguments are floating-point numbers, and at least one is a NaN,
198  * the result is undefined.
199  */
200 minmaxtype!(T) max(T...)(T arg){
201     static if(arg.length == 1) return arg[0];
202     else static if(arg.length == 2) return arg[1] > arg[0] ? arg[1] : arg[0];
203     static if(arg.length > 2) return max(arg[1] > arg[0] ? arg[1] : arg[0], arg[2..$]);
204 }
205 unittest
206 {
207     test(max('e', 'f')=='f');
208     test(min(3.5, 3.8)==3.5);
209     // check implicit conversion to integer.
210     test(min(3.5, 18)==3.5);
211 
212 }
213 
214 /** Returns the minimum number of x and y, favouring numbers over NaNs.
215  *
216  * If both x and y are numbers, the minimum is returned.
217  * If both parameters are NaN, either will be returned.
218  * If one parameter is a NaN and the other is a number, the number is
219  * returned (this behaviour is mandated by IEEE 754R, and is useful
220  * for determining the range of a function).
221  */
222 real minNum(real x, real y) {
223     if (x<=y || isNaN(y)) return x; else return y;
224 }
225 
226 /** Returns the maximum number of x and y, favouring numbers over NaNs.
227  *
228  * If both x and y are numbers, the maximum is returned.
229  * If both parameters are NaN, either will be returned.
230  * If one parameter is a NaN and the other is a number, the number is
231  * returned (this behaviour is mandated by IEEE 754-2008, and is useful
232  * for determining the range of a function).
233  */
234 real maxNum(real x, real y) {
235     if (x>=y || isNaN(y)) return x; else return y;
236 }
237 
238 /** Returns the minimum of x and y, favouring NaNs over numbers
239  *
240  * If both x and y are numbers, the minimum is returned.
241  * If both parameters are NaN, either will be returned.
242  * If one parameter is a NaN and the other is a number, the NaN is returned.
243  */
244 real minNaN(real x, real y) {
245     return (x<=y || isNaN(x))? x : y;
246 }
247 
248 /** Returns the maximum of x and y, favouring NaNs over numbers
249  *
250  * If both x and y are numbers, the maximum is returned.
251  * If both parameters are NaN, either will be returned.
252  * If one parameter is a NaN and the other is a number, the NaN is returned.
253  */
254 real maxNaN(real x, real y) {
255     return (x>=y || isNaN(x))? x : y;
256 }
257 
258 unittest
259 {
260     test(maxNum(NaN(0xABC), 56.1L)== 56.1L);
261     test(isIdentical(maxNaN(NaN(1389), 56.1L), NaN(1389)));
262     test(maxNum(28.0, NaN(0xABC))== 28.0);
263     test(minNum(1e12, NaN(0xABC))== 1e12);
264     test(isIdentical(minNaN(1e12, NaN(23454)), NaN(23454)));
265     test(isIdentical(minNum(NaN(489), NaN(23)), NaN(489)));
266 }
267 
268 /*
269  * Trig Functions
270  */
271 
272 /***********************************
273  * Returns cosine of x. x is in radians.
274  *
275  *      $(TABLE_SV
276  *      $(TR $(TH x)                 $(TH cos(x)) $(TH invalid?))
277  *      $(TR $(TD $(NAN))            $(TD $(NAN)) $(TD yes)     )
278  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(NAN)) $(TD yes)     )
279  *      )
280  * Bugs:
281  *      Results are undefined if |x| >= $(POWER 2,64).
282  */
283 
284 real cos(real x) /* intrinsic */
285 {
286     version(D_InlineAsm_X86)
287     {
288         asm
289         {
290             fld x;
291             fcos;
292         }
293     }
294     else
295     {
296         return core.stdc.math.cosl(x);
297     }
298 }
299 
300 unittest {
301     // NaN payloads
302     test(isIdentical(cos(NaN(314)), NaN(314)));
303 }
304 
305 /***********************************
306  * Returns sine of x. x is in radians.
307  *
308  *      $(TABLE_SV
309  *      $(TR $(TH x)               $(TH sin(x))      $(TH invalid?))
310  *      $(TR $(TD $(NAN))          $(TD $(NAN))      $(TD yes))
311  *      $(TR $(TD $(PLUSMN)0.0)    $(TD $(PLUSMN)0.0) $(TD no))
312  *      $(TR $(TD $(PLUSMNINF))    $(TD $(NAN))      $(TD yes))
313  *      )
314  * Bugs:
315  *      Results are undefined if |x| >= $(POWER 2,64).
316  */
317 real sin(real x) /* intrinsic */
318 {
319     version(D_InlineAsm_X86)
320     {
321         asm
322         {
323             fld x;
324             fsin;
325         }
326     }
327     else
328     {
329         return core.stdc.math.sinl(x);
330     }
331 }
332 
333 unittest {
334     // NaN payloads
335     test(isIdentical(sin(NaN(314)), NaN(314)));
336 }
337 
338 /**
339  * Returns tangent of x. x is in radians.
340  *
341  *	$(TABLE_SV
342  *	$(TR $(TH x)               $(TH tan(x))       $(TH invalid?))
343  *	$(TR $(TD $(NAN))          $(TD $(NAN))       $(TD yes))
344  *	$(TR $(TD $(PLUSMN)0.0)    $(TD $(PLUSMN)0.0) $(TD no))
345  *	$(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(NAN))     $(TD yes))
346  *	)
347  */
348 real tan(real x)
349 {
350     asm
351     {
352         fld x[EBP]      ; // load theta
353         fxam            ; // test for oddball values
354         fstsw   AX      ;
355         sahf            ;
356         jc  trigerr     ; // x is NAN, infinity, or empty
357                               // 387's can handle denormals
358 SC18:   fptan           ;
359         fstp    ST(0)   ; // dump X, which is always 1
360         fstsw   AX      ;
361         sahf            ;
362         jnp Lret        ; // C2 = 1 (x is out of range)
363 
364         // Do argument reduction to bring x into range
365         fldpi           ;
366         fxch            ;
367 SC17:   fprem1          ;
368         fstsw   AX      ;
369         sahf            ;
370         jp  SC17        ;
371         fstp    ST(1)   ; // remove pi from stack
372         jmp SC18        ;
373 
374 trigerr:
375         jnp Lret        ; // if x is NaN, return x.
376         fstp    ST(0)   ; // dump x, which will be infinity
377     }
378     return NaN(TANGO_NAN.TAN_DOMAIN);
379 Lret:
380     ;
381 }
382 
383 unittest
384 {
385     static real[2][] vals =     // angle,tan
386     [
387             [   0,   0],
388             [   .5,  .5463024898],
389             [   1,   1.557407725],
390             [   1.5, 14.10141995],
391             [   2,  -2.185039863],
392             [   2.5,-.7470222972],
393             [   3,  -.1425465431],
394             [   3.5, .3745856402],
395             [   4,   1.157821282],
396             [   4.5, 4.637332055],
397             [   5,  -3.380515006],
398             [   5.5,-.9955840522],
399             [   6,  -.2910061914],
400             [   6.5, .2202772003],
401             [   10,  .6483608275],
402 
403             // special angles
404             [   PI_4,   1],
405             //[   PI_2,   real.infinity], // PI_2 is not _exactly_ pi/2.
406             [   3*PI_4, -1],
407             [   PI,     0],
408             [   5*PI_4, 1],
409             //[   3*PI_2, -real.infinity],
410             [   7*PI_4, -1],
411             [   2*PI,   0],
412     ];
413     int i;
414 
415     for (i = 0; i < vals.length; i++)
416     {
417         real x = vals[i][0];
418         real r = vals[i][1];
419         real t = tan(x);
420 
421         //printf("tan(%Lg) = %Lg, should be %Lg\n", x, t, r);
422         if (!isIdentical(r, t)) test(fabs(r-t) <= .0000001);
423 
424         x = -x;
425         r = -r;
426         t = tan(x);
427         //printf("tan(%Lg) = %Lg, should be %Lg\n", x, t, r);
428         if (!isIdentical(r, t) && !(isNaN(r) && isNaN(t)))
429             test(fabs(r-t) <= .0000001);
430     }
431     // overflow
432     test(isNaN(tan(real.infinity)));
433     test(isNaN(tan(-real.infinity)));
434     // NaN propagation
435     test(isIdentical( tan(NaN(0x0123L)), NaN(0x0123L) ));
436 }
437 
438 /*****************************************
439  * Sine, cosine, and arctangent of multiple of &pi;
440  *
441  * Accuracy is preserved for large values of x.
442  */
443 real cosPi(real x)
444 {
445     return cos((x%2.0)*PI);
446 }
447 
448 /** ditto */
449 real sinPi(real x)
450 {
451     return sin((x%2.0)*PI);
452 }
453 
454 /** ditto */
455 real atanPi(real x)
456 {
457     return PI * atan(x); // BUG: Fix this.
458 }
459 
460 unittest {
461     test(isIdentical(sinPi(0.0), 0.0));
462     test(isIdentical(sinPi(-0.0), -0.0));
463     test(isIdentical(atanPi(0.0), 0.0));
464     test(isIdentical(atanPi(-0.0), -0.0));
465 }
466 
467 /***********************************
468  *  sine, complex and imaginary
469  *
470  *  sin(z) = sin(z.re)*cosh(z.im) + cos(z.re)*sinh(z.im)i
471  *
472  * If both sin(&theta;) and cos(&theta;) are required,
473  * it is most efficient to use expi(&theta).
474  */
475 creal sin(creal z)
476 {
477   creal cs = expi(z.re);
478   return cs.im * cosh(z.im) + cs.re * sinh(z.im) * 1i;
479 }
480 
481 /** ditto */
482 ireal sin(ireal y)
483 {
484   return cosh(y.im)*1i;
485 }
486 
487 unittest
488 {
489   test(sin(0.0+0.0i) == 0.0);
490   test(sin(2.0+0.0i) == sin(2.0L) );
491 }
492 
493 /***********************************
494  *  cosine, complex and imaginary
495  *
496  *  cos(z) = cos(z.re)*cosh(z.im) + sin(z.re)*sinh(z.im)i
497  */
498 creal cos(creal z)
499 {
500   creal cs = expi(z.re);
501   return cs.re * cosh(z.im) - cs.im * sinh(z.im) * 1i;
502 }
503 
504 /** ditto */
505 real cos(ireal y)
506 {
507   return cosh(y.im);
508 }
509 
510 unittest
511 {
512   test(cos(0.0+0.0i)==1.0);
513   test(cos(1.3L+0.0i)==cos(1.3L));
514   test(cos(5.2Li)== cosh(5.2L));
515 }
516 
517 /***************
518  * Calculates the arc cosine of x,
519  * returning a value ranging from 0 to $(PI).
520  *
521  *      $(TABLE_SV
522  *      $(TR $(TH x)         $(TH acos(x)) $(TH invalid?))
523  *      $(TR $(TD $(GT)1.0)  $(TD $(NAN))  $(TD yes))
524  *      $(TR $(TD $(LT)-1.0) $(TD $(NAN))  $(TD yes))
525  *      $(TR $(TD $(NAN))    $(TD $(NAN))  $(TD yes))
526  *      )
527  */
528 real acos(real x)
529 {
530     return core.stdc.math.acosl(x);
531 }
532 
533 unittest {
534     // NaN payloads
535     version(darwin){}
536     else {
537         test(isIdentical(acos(NaN(254)), NaN(254)));
538     }
539 }
540 
541 /***************
542  * Calculates the arc sine of x,
543  * returning a value ranging from -$(PI)/2 to $(PI)/2.
544  *
545  *      $(TABLE_SV
546  *      $(TR $(TH x)            $(TH asin(x))      $(TH invalid?))
547  *      $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no))
548  *      $(TR $(TD $(GT)1.0)     $(TD $(NAN))       $(TD yes))
549  *      $(TR $(TD $(LT)-1.0)    $(TD $(NAN))       $(TD yes))
550  *      )
551  */
552 real asin(real x)
553 {
554     return core.stdc.math.asinl(x);
555 }
556 
557 unittest {
558     // NaN payloads
559     version(darwin){}
560     else{
561         test(isIdentical(asin(NaN(7249)), NaN(7249)));
562     }
563 }
564 
565 /***************
566  * Calculates the arc tangent of x,
567  * returning a value ranging from -$(PI)/2 to $(PI)/2.
568  *
569  *      $(TABLE_SV
570  *      $(TR $(TH x)                 $(TH atan(x))      $(TH invalid?))
571  *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(PLUSMN)0.0) $(TD no))
572  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(NAN))       $(TD yes))
573  *      )
574  */
575 real atan(real x)
576 {
577     return core.stdc.math.atanl(x);
578 }
579 
580 unittest
581 {
582     // NaN payloads
583     test(isIdentical(atan(NaN(9876)), NaN(9876)));
584 }
585 
586 /***************
587  * Calculates the arc tangent of y / x,
588  * returning a value ranging from -$(PI) to $(PI).
589  *
590  *      $(TABLE_SV
591  *      $(TR $(TH y)                 $(TH x)            $(TH atan(y, x)))
592  *      $(TR $(TD $(NAN))            $(TD anything)     $(TD $(NAN)) )
593  *      $(TR $(TD anything)          $(TD $(NAN))       $(TD $(NAN)) )
594  *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(GT)0.0)     $(TD $(PLUSMN)0.0) )
595  *      $(TR $(TD $(PLUSMN)0.0)      $(TD +0.0)         $(TD $(PLUSMN)0.0) )
596  *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(LT)0.0)     $(TD $(PLUSMN)$(PI)))
597  *      $(TR $(TD $(PLUSMN)0.0)      $(TD -0.0)         $(TD $(PLUSMN)$(PI)))
598  *      $(TR $(TD $(GT)0.0)          $(TD $(PLUSMN)0.0) $(TD $(PI)/2) )
599  *      $(TR $(TD $(LT)0.0)          $(TD $(PLUSMN)0.0) $(TD -$(PI)/2) )
600  *      $(TR $(TD $(GT)0.0)          $(TD $(INFIN))     $(TD $(PLUSMN)0.0) )
601  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD anything)     $(TD $(PLUSMN)$(PI)/2))
602  *      $(TR $(TD $(GT)0.0)          $(TD -$(INFIN))    $(TD $(PLUSMN)$(PI)) )
603  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(INFIN))     $(TD $(PLUSMN)$(PI)/4))
604  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD -$(INFIN))    $(TD $(PLUSMN)3$(PI)/4))
605  *      )
606  */
607 real atan2(real y, real x)
608 {
609     return core.stdc.math.atan2l(y,x);
610 }
611 
612 unittest
613 {
614     // NaN payloads
615     test(isIdentical(atan2(5.3, NaN(9876)), NaN(9876)));
616     test(isIdentical(atan2(NaN(9876), 2.18), NaN(9876)));
617 }
618 
619 /***********************************
620  * Complex inverse sine
621  *
622  * asin(z) = -i log( sqrt(1-$(POWER z, 2)) + iz)
623  * where both log and sqrt are complex.
624  */
625 creal asin(creal z)
626 {
627     return -log(sqrt(1-z*z) + z*1i)*1i;
628 }
629 
630 unittest
631 {
632    test(asin(sin(0+0i)) == 0 + 0i);
633 }
634 
635 /***********************************
636  * Complex inverse cosine
637  *
638  * acos(z) = $(PI)/2 - asin(z)
639  */
640 creal acos(creal z)
641 {
642     return PI_2 - asin(z);
643 }
644 
645 
646 /***********************************
647  * Calculates the hyperbolic cosine of x.
648  *
649  *      $(TABLE_SV
650  *      $(TR $(TH x)                 $(TH cosh(x))      $(TH invalid?))
651  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)0.0) $(TD no) )
652  *      )
653  */
654 real cosh(real x)
655 {
656     //  cosh = (exp(x)+exp(-x))/2.
657     // The naive implementation works correctly.
658     real y = exp(x);
659     return (y + 1.0/y) * 0.5;
660 }
661 
662 unittest
663 {
664     // NaN payloads
665     test(isIdentical(cosh(NaN(432)), NaN(432)));
666 }
667 
668 /***********************************
669  * Calculates the hyperbolic sine of x.
670  *
671  *      $(TABLE_SV
672  *      $(TR $(TH x)                 $(TH sinh(x))           $(TH invalid?))
673  *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(PLUSMN)0.0)      $(TD no))
674  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)$(INFIN)) $(TD no))
675  *      )
676  */
677 real sinh(real x)
678 {
679     //  sinh(x) =  (exp(x)-exp(-x))/2;
680     // Very large arguments could cause an overflow, but
681     // the maximum value of x for which exp(x) + exp(-x)) != exp(x)
682     // is x = 0.5 * (real.mant_dig) * LN2. // = 22.1807 for real80.
683     if (fabs(x) > real.mant_dig * LN2) {
684         return copysign(0.5*exp(fabs(x)), x);
685     }
686     real y = expm1(x);
687     return 0.5 * y / (y+1) * (y+2);
688 }
689 
690 unittest
691 {
692     // NaN payloads
693     test(isIdentical(sinh(NaN(0xABC)), NaN(0xABC)));
694 }
695 
696 /***********************************
697  * Calculates the hyperbolic tangent of x.
698  *
699  *      $(TABLE_SV
700  *      $(TR $(TH x)                 $(TH tanh(x))      $(TH invalid?))
701  *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(PLUSMN)0.0) $(TD no) )
702  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)1.0) $(TD no))
703  *      )
704  */
705 real tanh(real x)
706 {
707     //  tanh(x) = (exp(x) - exp(-x))/(exp(x)+exp(-x))
708     if (fabs(x)> real.mant_dig * LN2){
709         return copysign(1, x);
710     }
711     real y = expm1(2*x);
712     return y/(y + 2);
713 }
714 
715 unittest
716 {
717     // NaN payloads
718     test(isIdentical(tanh(NaN(0xABC)), NaN(0xABC)));
719 }
720 
721 /***********************************
722  *  hyperbolic sine, complex and imaginary
723  *
724  *  sinh(z) = cos(z.im)*sinh(z.re) + sin(z.im)*cosh(z.re)i
725  */
726 creal sinh(creal z)
727 {
728   creal cs = expi(z.im);
729   return cs.re * sinh(z.re) + cs.im * cosh(z.re) * 1i;
730 }
731 
732 /** ditto */
733 ireal sinh(ireal y)
734 {
735   return sin(y.im)*1i;
736 }
737 
738 unittest
739 {
740   test(sinh(4.2L + 0i)==sinh(4.2L));
741 }
742 
743 /***********************************
744  *  hyperbolic cosine, complex and imaginary
745  *
746  *  cosh(z) = cos(z.im)*cosh(z.re) + sin(z.im)*sinh(z.re)i
747  */
748 creal cosh(creal z)
749 {
750   creal cs = expi(z.im);
751   return cs.re * cosh(z.re) + cs.im * sinh(z.re) * 1i;
752 }
753 
754 /** ditto */
755 real cosh(ireal y)
756 {
757   return cos(y.im);
758 }
759 
760 unittest
761 {
762   test(cosh(8.3L + 0i)==cosh(8.3L));
763 }
764 
765 
766 /***********************************
767  * Calculates the inverse hyperbolic cosine of x.
768  *
769  *  Mathematically, acosh(x) = log(x + sqrt( x*x - 1))
770  *
771  *    $(TABLE_SV
772  *    $(SVH  x,     acosh(x) )
773  *    $(SV  $(NAN), $(NAN) )
774  *    $(SV  $(LT)1,     $(NAN) )
775  *    $(SV  1,      0       )
776  *    $(SV  +$(INFIN),+$(INFIN))
777  *  )
778  */
779 real acosh(real x)
780 {
781     if (x > 1/real.epsilon)
782     return LN2 + log(x);
783     else
784     return log(x + sqrt(x*x - 1));
785 }
786 
787 unittest
788 {
789     test(isNaN(acosh(0.9)));
790     test(isNaN(acosh(real.nan)));
791     test(acosh(1)==0.0);
792     test(acosh(real.infinity) == real.infinity);
793     // NaN payloads
794     test(isIdentical(acosh(NaN(0xABC)), NaN(0xABC)));
795 }
796 
797 /***********************************
798  * Calculates the inverse hyperbolic sine of x.
799  *
800  *  Mathematically,
801  *  ---------------
802  *  asinh(x) =  log( x + sqrt( x*x + 1 )) // if x >= +0
803  *  asinh(x) = -log(-x + sqrt( x*x + 1 )) // if x <= -0
804  *  -------------
805  *
806  *    $(TABLE_SV
807  *    $(SVH x,                asinh(x)       )
808  *    $(SV  $(NAN),           $(NAN)         )
809  *    $(SV  $(PLUSMN)0,       $(PLUSMN)0      )
810  *    $(SV  $(PLUSMN)$(INFIN),$(PLUSMN)$(INFIN))
811  *    )
812  */
813 real asinh(real x)
814 {
815     if (ocean.math.IEEE.fabs(x) > 1 / real.epsilon) // beyond this point, x*x + 1 == x*x
816     return ocean.math.IEEE.copysign(LN2 + log(ocean.math.IEEE.fabs(x)), x);
817     else
818     {
819     // sqrt(x*x + 1) ==  1 + x * x / ( 1 + sqrt(x*x + 1) )
820     return ocean.math.IEEE.copysign(log1p(ocean.math.IEEE.fabs(x) + x*x / (1 + sqrt(x*x + 1)) ), x);
821     }
822 }
823 
824 unittest
825 {
826     test(isIdentical(0.0L,asinh(0.0)));
827     test(isIdentical(-0.0L,asinh(-0.0)));
828     test(asinh(real.infinity) == real.infinity);
829     test(asinh(-real.infinity) == -real.infinity);
830     test(isNaN(asinh(real.nan)));
831     // NaN payloads
832     test(isIdentical(asinh(NaN(0xABC)), NaN(0xABC)));
833 }
834 
835 /***********************************
836  * Calculates the inverse hyperbolic tangent of x,
837  * returning a value from ranging from -1 to 1.
838  *
839  * Mathematically, atanh(x) = log( (1+x)/(1-x) ) / 2
840  *
841  *
842  *    $(TABLE_SV
843  *    $(SVH  x,     acosh(x) )
844  *    $(SV  $(NAN), $(NAN) )
845  *    $(SV  $(PLUSMN)0, $(PLUSMN)0)
846  *    $(SV  -$(INFIN), -0)
847  *    )
848  */
849 real atanh(real x)
850 {
851     // log( (1+x)/(1-x) ) == log ( 1 + (2*x)/(1-x) )
852     return  0.5 * log1p( 2 * x / (1 - x) );
853 }
854 
855 unittest
856 {
857     test(isIdentical(0.0L, atanh(0.0)));
858     test(isIdentical(-0.0L,atanh(-0.0)));
859     test(isIdentical(atanh(-1),-real.infinity));
860     // Fails with -O because of DMD : https://issues.dlang.org/show_bug.cgi?id=13743
861     // Nothing can be done about it without feedback from DMD upstream
862     //assert(isIdentical(atanh(1),real.infinity));
863     test(isNaN(atanh(-real.infinity)));
864     // NaN payloads
865     test(isIdentical(atanh(NaN(0xABC)), NaN(0xABC)));
866 }
867 
868 /** ditto */
869 creal atanh(ireal y)
870 {
871     // Not optimised for accuracy or speed
872     return 0.5*(log(1+y) - log(1-y));
873 }
874 
875 /** ditto */
876 creal atanh(creal z)
877 {
878     // Not optimised for accuracy or speed
879     return 0.5 * (log(1 + z) - log(1-z));
880 }
881 
882 /*
883  * Powers and Roots
884  */
885 
886 /***************************************
887  * Compute square root of x.
888  *
889  *      $(TABLE_SV
890  *      $(TR $(TH x)         $(TH sqrt(x))   $(TH invalid?))
891  *      $(TR $(TD -0.0)      $(TD -0.0)      $(TD no))
892  *      $(TR $(TD $(LT)0.0)  $(TD $(NAN))    $(TD yes))
893  *      $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no))
894  *      )
895  */
896 float sqrt(float x) /* intrinsic */
897 {
898     version(D_InlineAsm_X86)
899     {
900         asm
901         {
902             fld x;
903             fsqrt;
904         }
905     }
906     else
907     {
908         return core.stdc.math.sqrtf(x);
909     }
910 }
911 
912 double sqrt(double x) /* intrinsic */ /// ditto
913 {
914     version(D_InlineAsm_X86)
915     {
916         asm
917         {
918             fld x;
919             fsqrt;
920         }
921     }
922     else
923     {
924         return core.stdc.math.sqrt(x);
925     }
926 }
927 
928 real sqrt(real x) /* intrinsic */ /// ditto
929 {
930     version(D_InlineAsm_X86)
931     {
932         asm
933         {
934             fld x;
935             fsqrt;
936         }
937     }
938     else
939     {
940         return core.stdc.math.sqrtl(x);
941     }
942 }
943 
944 /** ditto */
945 creal sqrt(creal z)
946 {
947 
948     if (z == 0.0) return z;
949     real x,y,w,r;
950     creal c;
951 
952     x = ocean.math.IEEE.fabs(z.re);
953     y = ocean.math.IEEE.fabs(z.im);
954     if (x >= y) {
955         r = y / x;
956         w = sqrt(x) * sqrt(0.5 * (1 + sqrt(1 + r * r)));
957     } else  {
958         r = x / y;
959         w = sqrt(y) * sqrt(0.5 * (r + sqrt(1 + r * r)));
960     }
961 
962     if (z.re >= 0) {
963         c = w + (z.im / (w + w)) * 1.0i;
964     } else {
965         if (z.im < 0)  w = -w;
966         c = z.im / (w + w) + w * 1.0i;
967     }
968     return c;
969 }
970 
971 unittest {
972     // NaN payloads
973     test(isIdentical(sqrt(NaN(0xABC)), NaN(0xABC)));
974     test(sqrt(-1+0i) == 1i);
975     test(isIdentical(sqrt(0-0i), 0-0i));
976     test(cfeqrel(sqrt(4+16i)*sqrt(4+16i), 4+16i)>=real.mant_dig-2);
977 }
978 
979 /***************
980  * Calculates the cube root of x.
981  *
982  *      $(TABLE_SV
983  *      $(TR $(TH $(I x))            $(TH cbrt(x))           $(TH invalid?))
984  *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(PLUSMN)0.0)      $(TD no) )
985  *      $(TR $(TD $(NAN))            $(TD $(NAN))            $(TD yes) )
986  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)$(INFIN)) $(TD no) )
987  *      )
988  */
989 real cbrt(real x)
990 {
991     return core.stdc.math.cbrtl(x);
992 }
993 
994 
995 unittest {
996     // NaN payloads
997     test(isIdentical(cbrt(NaN(0xABC)), NaN(0xABC)));
998 }
999 
1000 public:
1001 
1002 /**
1003  * Calculates e$(SUP x).
1004  *
1005  *  $(TABLE_SV
1006  *    $(TR $(TH x)                   $(TH e$(SUP x)) )
1007  *    $(TR $(TD +$(INFIN))           $(TD +$(INFIN)) )
1008  *    $(TR $(TD -$(INFIN))           $(TD +0.0)      )
1009  *    $(TR $(TD $(NAN))              $(TD $(NAN))    )
1010  *  )
1011  */
1012 real exp(real x) {
1013     version(Naked_D_InlineAsm_X86) {
1014    //  e^x = 2^(LOG2E*x)
1015    // (This is valid because the overflow & underflow limits for exp
1016    // and exp2 are so similar).
1017     return exp2(LOG2E*x);
1018     } else {
1019         return core.stdc.math.expl(x);
1020     }
1021 }
1022 
1023 /**
1024  * Calculates the value of the natural logarithm base (e)
1025  * raised to the power of x, minus 1.
1026  *
1027  * For very small x, expm1(x) is more accurate
1028  * than exp(x)-1.
1029  *
1030  *  $(TABLE_SV
1031  *    $(TR $(TH x)             $(TH e$(SUP x)-1)  )
1032  *    $(TR $(TD $(PLUSMN)0.0)  $(TD $(PLUSMN)0.0) )
1033  *    $(TR $(TD +$(INFIN))     $(TD +$(INFIN))    )
1034  *    $(TR $(TD -$(INFIN))     $(TD -1.0)         )
1035  *    $(TR $(TD $(NAN))        $(TD $(NAN))       )
1036  *  )
1037  */
1038 real expm1(real x)
1039 {
1040     version(Naked_D_InlineAsm_X86) {
1041       enum { PARAMSIZE = (real.sizeof+3)&(0xFFFF_FFFC) } // always a multiple of 4
1042       asm {
1043         /*  expm1() for x87 80-bit reals, IEEE754-2008 conformant.
1044          * Author: Don Clugston.
1045          *
1046          *    expm1(x) = 2^(rndint(y))* 2^(y-rndint(y)) - 1 where y = LN2*x.
1047          *    = 2rndy * 2ym1 + 2rndy - 1, where 2rndy = 2^(rndint(y))
1048          *     and 2ym1 = (2^(y-rndint(y))-1).
1049          *    If 2rndy  < 0.5*real.epsilon, result is -1.
1050          *    Implementation is otherwise the same as for exp2()
1051          */
1052         naked;
1053         fld real ptr [ESP+4] ; // x
1054         mov AX, [ESP+4+8]; // AX = exponent and sign
1055         sub ESP, 12+8; // Create scratch space on the stack
1056         // [ESP,ESP+2] = scratchint
1057         // [ESP+4..+6, +8..+10, +10] = scratchreal
1058         // set scratchreal mantissa = 1.0
1059         mov dword ptr [ESP+8], 0;
1060         mov dword ptr [ESP+8+4], 0x80000000;
1061         and AX, 0x7FFF; // drop sign bit
1062         cmp AX, 0x401D; // avoid InvalidException in fist
1063         jae L_extreme;
1064         fldl2e;
1065         fmul ; // y = x*log2(e)
1066         fist dword ptr [ESP]; // scratchint = rndint(y)
1067         fisub dword ptr [ESP]; // y - rndint(y)
1068         // and now set scratchreal exponent
1069         mov EAX, [ESP];
1070         add EAX, 0x3fff;
1071         jle short L_largenegative;
1072         cmp EAX,0x8000;
1073         jge short L_largepositive;
1074         mov [ESP+8+8],AX;
1075         f2xm1; // 2^(y-rndint(y)) -1
1076         fld real ptr [ESP+8] ; // 2^rndint(y)
1077         fmul ST(1), ST;
1078         fld1;
1079         fsubp ST(1), ST;
1080         fadd;
1081         add ESP,12+8;
1082         ret PARAMSIZE;
1083 
1084 L_extreme: // Extreme exponent. X is very large positive, very
1085         // large negative, infinity, or NaN.
1086         fxam;
1087         fstsw AX;
1088         test AX, 0x0400; // NaN_or_zero, but we already know x!=0
1089         jz L_was_nan;  // if x is NaN, returns x
1090         test AX, 0x0200;
1091         jnz L_largenegative;
1092 L_largepositive:
1093         // Set scratchreal = real.max.
1094         // squaring it will create infinity, and set overflow flag.
1095         mov word  ptr [ESP+8+8], 0x7FFE;
1096         fstp ST(0);
1097         fld real ptr [ESP+8];  // load scratchreal
1098         fmul ST(0), ST;        // square it, to create havoc!
1099 L_was_nan:
1100         add ESP,12+8;
1101         ret PARAMSIZE;
1102 L_largenegative:
1103         fstp ST(0);
1104         fld1;
1105         fchs; // return -1. Underflow flag is not set.
1106         add ESP,12+8;
1107         ret PARAMSIZE;
1108       }
1109     } else version(D_InlineAsm_X86_64) {
1110         asm
1111         {
1112             naked;
1113         }
1114         asm
1115         {
1116             fld real ptr [RSP+8]; // x
1117             mov AX,[RSP+8+8]; // AX = exponent and sign
1118         }
1119         asm
1120         {
1121 /* expm1() for x87 80-bit reals, IEEE754-2008 conformant.
1122 * Author: Don Clugston.
1123 *
1124 * expm1(x) = 2^(rndint(y))* 2^(y-rndint(y)) - 1 where y = LN2*x.
1125 * = 2rndy * 2ym1 + 2rndy - 1, where 2rndy = 2^(rndint(y))
1126 * and 2ym1 = (2^(y-rndint(y))-1).
1127 * If 2rndy < 0.5*real.epsilon, result is -1.
1128 * Implementation is otherwise the same as for exp2()
1129 */
1130             sub RSP, 24; // Create scratch space on the stack
1131             // [RSP,RSP+2] = scratchint
1132             // [RSP+4..+6, +8..+10, +10] = scratchreal
1133             // set scratchreal mantissa = 1.0
1134             mov dword ptr [RSP+8], 0;
1135             mov dword ptr [RSP+8+4], 0x80000000;
1136             and AX, 0x7FFF; // drop sign bit
1137             cmp AX, 0x401D; // avoid InvalidException in fist
1138             jae L_extreme;
1139             fldl2e;
1140             fmul ; // y = x*log2(e)
1141             fist dword ptr [RSP]; // scratchint = rndint(y)
1142             fisub dword ptr [RSP]; // y - rndint(y)
1143             // and now set scratchreal exponent
1144             mov EAX, [RSP];
1145             add EAX, 0x3fff;
1146             jle short L_largenegative;
1147             cmp EAX,0x8000;
1148             jge short L_largepositive;
1149             mov [RSP+8+8],AX;
1150             f2xm1; // 2^(y-rndint(y)) -1
1151             fld real ptr [RSP+8] ; // 2^rndint(y)
1152             fmul ST(1), ST;
1153             fld1;
1154             fsubp ST(1), ST;
1155             fadd;
1156             add RSP,24;
1157             ret;
1158 
1159 L_extreme: // Extreme exponent. X is very large positive, very
1160             // large negative, infinity, or NaN.
1161             fxam;
1162             fstsw AX;
1163             test AX, 0x0400; // NaN_or_zero, but we already know x!=0
1164             jz L_was_nan; // if x is NaN, returns x
1165             test AX, 0x0200;
1166             jnz L_largenegative;
1167 L_largepositive:
1168             // Set scratchreal = real.max.
1169             // squaring it will create infinity, and set overflow flag.
1170             mov word ptr [RSP+8+8], 0x7FFE;
1171             fstp ST(0);
1172             fld real ptr [RSP+8]; // load scratchreal
1173             fmul ST(0), ST; // square it, to create havoc!
1174 L_was_nan:
1175             add RSP,24;
1176             ret;
1177 
1178 L_largenegative:
1179             fstp ST(0);
1180             fld1;
1181             fchs; // return -1. Underflow flag is not set.
1182             add RSP,24;
1183             ret;
1184         }
1185     } else {
1186         return core.stdc.math.expm1l(x);
1187     }
1188 }
1189 
1190 /**
1191  * Calculates 2$(SUP x).
1192  *
1193  *  $(TABLE_SV
1194  *    $(TR $(TH x)             $(TH exp2(x))   )
1195  *    $(TR $(TD +$(INFIN))     $(TD +$(INFIN)) )
1196  *    $(TR $(TD -$(INFIN))     $(TD +0.0)      )
1197  *    $(TR $(TD $(NAN))        $(TD $(NAN))    )
1198  *  )
1199  */
1200 real exp2(real x)
1201 {
1202     version(Naked_D_InlineAsm_X86) {
1203       enum { PARAMSIZE = (real.sizeof+3)&(0xFFFF_FFFC) } // always a multiple of 4
1204       asm {
1205         /*  exp2() for x87 80-bit reals, IEEE754-2008 conformant.
1206          * Author: Don Clugston.
1207          *
1208          * exp2(x) = 2^(rndint(x))* 2^(y-rndint(x))
1209          * The trick for high performance is to avoid the fscale(28cycles on core2),
1210          * frndint(19 cycles), leaving f2xm1(19 cycles) as the only slow instruction.
1211          *
1212          * We can do frndint by using fist. BUT we can't use it for huge numbers,
1213          * because it will set the Invalid Operation flag is overflow or NaN occurs.
1214          * Fortunately, whenever this happens the result would be zero or infinity.
1215          *
1216          * We can perform fscale by directly poking into the exponent. BUT this doesn't
1217          * work for the (very rare) cases where the result is subnormal. So we fall back
1218          * to the slow method in that case.
1219          */
1220         naked;
1221         fld real ptr [ESP+4] ; // x
1222         mov AX, [ESP+4+8]; // AX = exponent and sign
1223         sub ESP, 12+8; // Create scratch space on the stack
1224         // [ESP,ESP+2] = scratchint
1225         // [ESP+4..+6, +8..+10, +10] = scratchreal
1226         // set scratchreal mantissa = 1.0
1227         mov dword ptr [ESP+8], 0;
1228         mov dword ptr [ESP+8+4], 0x80000000;
1229         and AX, 0x7FFF; // drop sign bit
1230         cmp AX, 0x401D; // avoid InvalidException in fist
1231         jae L_extreme;
1232         fist dword ptr [ESP]; // scratchint = rndint(x)
1233         fisub dword ptr [ESP]; // x - rndint(x)
1234         // and now set scratchreal exponent
1235         mov EAX, [ESP];
1236         add EAX, 0x3fff;
1237         jle short L_subnormal;
1238         cmp EAX,0x8000;
1239         jge short L_overflow;
1240         mov [ESP+8+8],AX;
1241 L_normal:
1242         f2xm1;
1243         fld1;
1244         fadd; // 2^(x-rndint(x))
1245         fld real ptr [ESP+8] ; // 2^rndint(x)
1246         add ESP,12+8;
1247         fmulp ST(1), ST;
1248         ret PARAMSIZE;
1249 
1250 L_subnormal:
1251         // Result will be subnormal.
1252         // In this rare case, the simple poking method doesn't work.
1253         // The speed doesn't matter, so use the slow fscale method.
1254         fild dword ptr [ESP];  // scratchint
1255         fld1;
1256         fscale;
1257         fstp real ptr [ESP+8]; // scratchreal = 2^scratchint
1258         fstp ST(0);         // drop scratchint
1259         jmp L_normal;
1260 
1261 L_extreme: // Extreme exponent. X is very large positive, very
1262         // large negative, infinity, or NaN.
1263         fxam;
1264         fstsw AX;
1265         test AX, 0x0400; // NaN_or_zero, but we already know x!=0
1266         jz L_was_nan;  // if x is NaN, returns x
1267         // set scratchreal = real.min
1268         // squaring it will return 0, setting underflow flag
1269         mov word  ptr [ESP+8+8], 1;
1270         test AX, 0x0200;
1271         jnz L_waslargenegative;
1272 L_overflow:
1273         // Set scratchreal = real.max.
1274         // squaring it will create infinity, and set overflow flag.
1275         mov word  ptr [ESP+8+8], 0x7FFE;
1276 L_waslargenegative:
1277         fstp ST(0);
1278         fld real ptr [ESP+8];  // load scratchreal
1279         fmul ST(0), ST;        // square it, to create havoc!
1280 L_was_nan:
1281         add ESP,12+8;
1282         ret PARAMSIZE;
1283       }
1284     } else version(D_InlineAsm_X86_64) {
1285         asm
1286         {
1287             naked;
1288         }
1289         asm
1290         {
1291             fld real ptr [RSP+8]; // x
1292             mov AX,[RSP+8+8]; // AX = exponent and sign
1293         }
1294         asm
1295         {
1296 /* exp2() for x87 80-bit reals, IEEE754-2008 conformant.
1297  * Author: Don Clugston.
1298  *
1299  * exp2(x) = 2^(rndint(x))* 2^(y-rndint(x))
1300  * The trick for high performance is to avoid the fscale(28cycles on core2),
1301  * frndint(19 cycles), leaving f2xm1(19 cycles) as the only slow instruction.
1302  *
1303  * We can do frndint by using fist. BUT we can't use it for huge numbers,
1304  * because it will set the Invalid Operation flag is overflow or NaN occurs.
1305  * Fortunately, whenever this happens the result would be zero or infinity.
1306  *
1307  * We can perform fscale by directly poking into the exponent. BUT this doesn't
1308  * work for the (very rare) cases where the result is subnormal. So we fall back
1309  * to the slow method in that case.
1310  */
1311             sub RSP, 24; // Create scratch space on the stack
1312             // [RSP,RSP+2] = scratchint
1313             // [RSP+4..+6, +8..+10, +10] = scratchreal
1314             // set scratchreal mantissa = 1.0
1315             mov dword ptr [RSP+8], 0;
1316             mov dword ptr [RSP+8+4], 0x80000000;
1317             and AX, 0x7FFF; // drop sign bit
1318             cmp AX, 0x401D; // avoid InvalidException in fist
1319             jae L_extreme;
1320             fist dword ptr [RSP]; // scratchint = rndint(x)
1321             fisub dword ptr [RSP]; // x - rndint(x)
1322             // and now set scratchreal exponent
1323             mov EAX, [RSP];
1324             add EAX, 0x3fff;
1325             jle short L_subnormal;
1326             cmp EAX,0x8000;
1327             jge short L_overflow;
1328             mov [RSP+8+8],AX;
1329 L_normal:
1330             f2xm1;
1331             fld1;
1332             fadd; // 2^(x-rndint(x))
1333             fld real ptr [RSP+8] ; // 2^rndint(x)
1334             add RSP,24;
1335             fmulp ST(1), ST;
1336             ret;
1337 
1338 L_subnormal:
1339             // Result will be subnormal.
1340             // In this rare case, the simple poking method doesn't work.
1341             // The speed doesn't matter, so use the slow fscale method.
1342             fild dword ptr [RSP]; // scratchint
1343             fld1;
1344             fscale;
1345             fstp real ptr [RSP+8]; // scratchreal = 2^scratchint
1346             fstp ST(0); // drop scratchint
1347             jmp L_normal;
1348 
1349 L_extreme: // Extreme exponent. X is very large positive, very
1350             // large negative, infinity, or NaN.
1351             fxam;
1352             fstsw AX;
1353             test AX, 0x0400; // NaN_or_zero, but we already know x!=0
1354             jz L_was_nan; // if x is NaN, returns x
1355             // set scratchreal = real.min
1356             // squaring it will return 0, setting underflow flag
1357             mov word ptr [RSP+8+8], 1;
1358             test AX, 0x0200;
1359             jnz L_waslargenegative;
1360 L_overflow:
1361             // Set scratchreal = real.max.
1362             // squaring it will create infinity, and set overflow flag.
1363             mov word ptr [RSP+8+8], 0x7FFE;
1364 L_waslargenegative:
1365             fstp ST(0);
1366             fld real ptr [RSP+8]; // load scratchreal
1367             fmul ST(0), ST; // square it, to create havoc!
1368 L_was_nan:
1369             add RSP,24;
1370             ret;
1371         }
1372     } else {
1373         return core.stdc.math.exp2l(x);
1374     }
1375 }
1376 
1377 unittest {
1378     // NaN payloads
1379     test(isIdentical(exp(NaN(0xABC)), NaN(0xABC)));
1380 }
1381 
1382 unittest {
1383     // NaN payloads
1384     test(isIdentical(expm1(NaN(0xABC)), NaN(0xABC)));
1385 }
1386 
1387 unittest {
1388     // NaN payloads
1389     test(isIdentical(exp2(NaN(0xABC)), NaN(0xABC)));
1390 }
1391 
1392 /*
1393  * Powers and Roots
1394  */
1395 
1396 /**************************************
1397  * Calculate the natural logarithm of x.
1398  *
1399  *    $(TABLE_SV
1400  *    $(TR $(TH x)            $(TH log(x))    $(TH divide by 0?) $(TH invalid?))
1401  *    $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes)          $(TD no))
1402  *    $(TR $(TD $(LT)0.0)     $(TD $(NAN))    $(TD no)           $(TD yes))
1403  *    $(TR $(TD +$(INFIN))    $(TD +$(INFIN)) $(TD no)           $(TD no))
1404  *    )
1405  */
1406 real log(real x)
1407 {
1408     return core.stdc.math.logl(x);
1409 }
1410 
1411 unittest {
1412     // NaN payloads
1413     test(isIdentical(log(NaN(0xABC)), NaN(0xABC)));
1414 }
1415 
1416 /******************************************
1417  *      Calculates the natural logarithm of 1 + x.
1418  *
1419  *      For very small x, log1p(x) will be more accurate than
1420  *      log(1 + x).
1421  *
1422  *  $(TABLE_SV
1423  *  $(TR $(TH x)            $(TH log1p(x))     $(TH divide by 0?) $(TH invalid?))
1424  *  $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no)           $(TD no))
1425  *  $(TR $(TD -1.0)         $(TD -$(INFIN))    $(TD yes)          $(TD no))
1426  *  $(TR $(TD $(LT)-1.0)    $(TD $(NAN))       $(TD no)           $(TD yes))
1427  *  $(TR $(TD +$(INFIN))    $(TD -$(INFIN))    $(TD no)           $(TD no))
1428  *  )
1429  */
1430 real log1p(real x)
1431 {
1432     return core.stdc.math.log1pl(x);
1433 }
1434 
1435 unittest {
1436     // NaN payloads
1437     test(isIdentical(log1p(NaN(0xABC)), NaN(0xABC)));
1438 }
1439 
1440 /***************************************
1441  * Calculates the base-2 logarithm of x:
1442  * $(SUB log, 2)x
1443  *
1444  *  $(TABLE_SV
1445  *  $(TR $(TH x)            $(TH log2(x))   $(TH divide by 0?) $(TH invalid?))
1446  *  $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes)          $(TD no) )
1447  *  $(TR $(TD $(LT)0.0)     $(TD $(NAN))    $(TD no)           $(TD yes) )
1448  *  $(TR $(TD +$(INFIN))    $(TD +$(INFIN)) $(TD no)           $(TD no) )
1449  *  )
1450  */
1451 real log2(real x)
1452 {
1453     return core.stdc.math.log2l(x);
1454 }
1455 
1456 unittest {
1457     // NaN payloads
1458     test(isIdentical(log2(NaN(0xABC)), NaN(0xABC)));
1459 }
1460 
1461 /**************************************
1462  * Calculate the base-10 logarithm of x.
1463  *
1464  *      $(TABLE_SV
1465  *      $(TR $(TH x)            $(TH log10(x))  $(TH divide by 0?) $(TH invalid?))
1466  *      $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes)          $(TD no))
1467  *      $(TR $(TD $(LT)0.0)     $(TD $(NAN))    $(TD no)           $(TD yes))
1468  *      $(TR $(TD +$(INFIN))    $(TD +$(INFIN)) $(TD no)           $(TD no))
1469  *      )
1470  */
1471 real log10(real x)
1472 {
1473     return core.stdc.math.log10l(x);
1474 }
1475 
1476 unittest {
1477     // NaN payloads
1478     test(isIdentical(log10(NaN(0xABC)), NaN(0xABC)));
1479 }
1480 
1481 /***********************************
1482  * Exponential, complex and imaginary
1483  *
1484  * For complex numbers, the exponential function is defined as
1485  *
1486  *  exp(z) = exp(z.re)cos(z.im) + exp(z.re)sin(z.im)i.
1487  *
1488  *  For a pure imaginary argument,
1489  *  exp(&theta;i)  = cos(&theta;) + sin(&theta;)i.
1490  *
1491  */
1492 creal exp(ireal y)
1493 {
1494    return expi(y.im);
1495 }
1496 
1497 /** ditto */
1498 creal exp(creal z)
1499 {
1500   return expi(z.im) * exp(z.re);
1501 }
1502 
1503 unittest {
1504     test(exp(1.3e5Li)==cos(1.3e5L)+sin(1.3e5L)*1i);
1505     test(exp(0.0Li)==1L+0.0Li);
1506     test(exp(7.2 + 0.0i) == exp(7.2L));
1507     creal c = exp(ireal.nan);
1508     test(isNaN(c.re) && isNaN(c.im));
1509     c = exp(ireal.infinity);
1510     test(isNaN(c.re) && isNaN(c.im));
1511 }
1512 
1513 /***********************************
1514  *  Natural logarithm, complex
1515  *
1516  * Returns complex logarithm to the base e (2.718...) of
1517  * the complex argument x.
1518  *
1519  * If z = x + iy, then
1520  *       log(z) = log(abs(z)) + i arctan(y/x).
1521  *
1522  * The arctangent ranges from -PI to +PI.
1523  * There are branch cuts along both the negative real and negative
1524  * imaginary axes. For pure imaginary arguments, use one of the
1525  * following forms, depending on which branch is required.
1526  * ------------
1527  *    log( 0.0 + yi) = log(-y) + PI_2i  // y<=-0.0
1528  *    log(-0.0 + yi) = log(-y) - PI_2i  // y<=-0.0
1529  * ------------
1530  */
1531 creal log(creal z)
1532 {
1533   return log(abs(z)) + atan2(z.im, z.re)*1i;
1534 }
1535 
1536 /*
1537  * feqrel for complex numbers. Returns the worst relative
1538  * equality of the two components.
1539  */
1540 private int cfeqrel(creal a, creal b)
1541 {
1542     int intmin(int a, int b) { return a<b? a: b; }
1543     return intmin(feqrel(a.re, b.re), feqrel(a.im, b.im));
1544 }
1545 
1546 unittest {
1547 
1548   test(log(3.0L +0i) == log(3.0L)+0i);
1549   test(cfeqrel(log(0.0L-2i),( log(2.0L)-PI_2*1i)) >= real.mant_dig-10);
1550   test(cfeqrel(log(0.0L+2i),( log(2.0L)+PI_2*1i)) >= real.mant_dig-10);
1551 }
1552 
1553 /**
1554  * Fast integral powers.
1555  */
1556 real pow(real x, uint n)
1557 {
1558     real p;
1559 
1560     switch (n)
1561     {
1562     case 0:
1563         p = 1.0;
1564         break;
1565 
1566     case 1:
1567         p = x;
1568         break;
1569 
1570     case 2:
1571         p = x * x;
1572         break;
1573 
1574     default:
1575         p = 1.0;
1576         while (1){
1577             if (n & 1)
1578                 p *= x;
1579             n >>= 1;
1580             if (!n)
1581                 break;
1582             x *= x;
1583         }
1584         break;
1585     }
1586     return p;
1587 }
1588 
1589 /** ditto */
1590 real pow(real x, int n)
1591 {
1592     if (n < 0) return pow(x, cast(real)n);
1593     else return pow(x, cast(uint)n);
1594 }
1595 
1596 /*********************************************
1597  * Calculates x$(SUP y).
1598  *
1599  * $(TABLE_SV
1600  * $(TR $(TH x) $(TH y) $(TH pow(x, y))
1601  *      $(TH div 0) $(TH invalid?))
1602  * $(TR $(TD anything)      $(TD $(PLUSMN)0.0)                $(TD 1.0)
1603  *      $(TD no)        $(TD no) )
1604  * $(TR $(TD |x| $(GT) 1)    $(TD +$(INFIN))                  $(TD +$(INFIN))
1605  *      $(TD no)        $(TD no) )
1606  * $(TR $(TD |x| $(LT) 1)    $(TD +$(INFIN))                  $(TD +0.0)
1607  *      $(TD no)        $(TD no) )
1608  * $(TR $(TD |x| $(GT) 1)    $(TD -$(INFIN))                  $(TD +0.0)
1609  *      $(TD no)        $(TD no) )
1610  * $(TR $(TD |x| $(LT) 1)    $(TD -$(INFIN))                  $(TD +$(INFIN))
1611  *      $(TD no)        $(TD no) )
1612  * $(TR $(TD +$(INFIN))      $(TD $(GT) 0.0)                  $(TD +$(INFIN))
1613  *      $(TD no)        $(TD no) )
1614  * $(TR $(TD +$(INFIN))      $(TD $(LT) 0.0)                  $(TD +0.0)
1615  *      $(TD no)        $(TD no) )
1616  * $(TR $(TD -$(INFIN))      $(TD odd integer $(GT) 0.0)      $(TD -$(INFIN))
1617  *      $(TD no)        $(TD no) )
1618  * $(TR $(TD -$(INFIN))      $(TD $(GT) 0.0, not odd integer) $(TD +$(INFIN))
1619  *      $(TD no)        $(TD no))
1620  * $(TR $(TD -$(INFIN))      $(TD odd integer $(LT) 0.0)      $(TD -0.0)
1621  *      $(TD no)        $(TD no) )
1622  * $(TR $(TD -$(INFIN))      $(TD $(LT) 0.0, not odd integer) $(TD +0.0)
1623  *      $(TD no)        $(TD no) )
1624  * $(TR $(TD $(PLUSMN)1.0)   $(TD $(PLUSMN)$(INFIN))          $(TD $(NAN))
1625  *      $(TD no)        $(TD yes) )
1626  * $(TR $(TD $(LT) 0.0)      $(TD finite, nonintegral)        $(TD $(NAN))
1627  *      $(TD no)        $(TD yes))
1628  * $(TR $(TD $(PLUSMN)0.0)   $(TD odd integer $(LT) 0.0)      $(TD $(PLUSMNINF))
1629  *      $(TD yes)       $(TD no) )
1630  * $(TR $(TD $(PLUSMN)0.0)   $(TD $(LT) 0.0, not odd integer) $(TD +$(INFIN))
1631  *      $(TD yes)       $(TD no))
1632  * $(TR $(TD $(PLUSMN)0.0)   $(TD odd integer $(GT) 0.0)      $(TD $(PLUSMN)0.0)
1633  *      $(TD no)        $(TD no) )
1634  * $(TR $(TD $(PLUSMN)0.0)   $(TD $(GT) 0.0, not odd integer) $(TD +0.0)
1635  *      $(TD no)        $(TD no) )
1636  * )
1637  */
1638 real pow(real x, real y)
1639 {
1640     if (isNaN(y))
1641         return y;
1642 
1643     if (y == 0)
1644         return 1;       // even if x is $(NAN)
1645     if (isNaN(x) && y != 0)
1646         return x;
1647     if (isInfinity(y))
1648     {
1649         if (ocean.math.IEEE.fabs(x) > 1)
1650         {
1651             if (signbit(y))
1652                 return +0.0;
1653             else
1654                 return real.infinity;
1655         }
1656         else if (ocean.math.IEEE.fabs(x) == 1)
1657         {
1658             return NaN(TANGO_NAN.POW_DOMAIN);
1659         }
1660         else // < 1
1661         {
1662             if (signbit(y))
1663                 return real.infinity;
1664             else
1665                 return +0.0;
1666         }
1667     }
1668     if (isInfinity(x))
1669     {
1670         if (signbit(x))
1671         {
1672             long i;
1673             i = cast(long)y;
1674             if (y > 0)
1675             {
1676                 if (i == y && i & 1)
1677                 return -real.infinity;
1678                 else
1679                 return real.infinity;
1680             }
1681             else if (y < 0)
1682             {
1683                 if (i == y && i & 1)
1684                 return -0.0;
1685                 else
1686                 return +0.0;
1687             }
1688         }
1689         else
1690         {
1691             if (y > 0)
1692                 return real.infinity;
1693             else if (y < 0)
1694                 return +0.0;
1695         }
1696     }
1697 
1698     if (x == 0.0)
1699     {
1700         if (signbit(x))
1701         {
1702             long i;
1703 
1704             i = cast(long)y;
1705             if (y > 0)
1706             {
1707                 if (i == y && i & 1)
1708                 return -0.0;
1709                 else
1710                 return +0.0;
1711             }
1712             else if (y < 0)
1713             {
1714                 if (i == y && i & 1)
1715                 return -real.infinity;
1716                 else
1717                 return real.infinity;
1718             }
1719         }
1720         else
1721         {
1722             if (y > 0)
1723                 return +0.0;
1724             else if (y < 0)
1725                 return real.infinity;
1726         }
1727     }
1728 
1729     return core.stdc.math.powl(x, y);
1730 }
1731 
1732 unittest
1733 {
1734     real x = 46;
1735 
1736     test(pow(x,0) == 1.0);
1737     test(pow(x,1) == x);
1738     test(pow(x,2) == x * x);
1739     test(pow(x,3) == x * x * x);
1740     test(pow(x,8) == (x * x) * (x * x) * (x * x) * (x * x));
1741     // NaN payloads
1742     test(isIdentical(pow(NaN(0xABC), 19), NaN(0xABC)));
1743 }
1744 
1745 /***********************************************************************
1746  * Calculates the length of the
1747  * hypotenuse of a right-angled triangle with sides of length x and y.
1748  * The hypotenuse is the value of the square root of
1749  * the sums of the squares of x and y:
1750  *
1751  *      sqrt($(POW x, 2) + $(POW y, 2))
1752  *
1753  * Note that hypot(x, y), hypot(y, x) and
1754  * hypot(x, -y) are equivalent.
1755  *
1756  *  $(TABLE_SV
1757  *  $(TR $(TH x)            $(TH y)            $(TH hypot(x, y)) $(TH invalid?))
1758  *  $(TR $(TD x)            $(TD $(PLUSMN)0.0) $(TD |x|)         $(TD no))
1759  *  $(TR $(TD $(PLUSMNINF)) $(TD y)            $(TD +$(INFIN))   $(TD no))
1760  *  $(TR $(TD $(PLUSMNINF)) $(TD $(NAN))       $(TD +$(INFIN))   $(TD no))
1761  *  )
1762  */
1763 real hypot(real x, real y)
1764 {
1765     // Scale x and y to avoid underflow and overflow.
1766     // If one is huge and the other tiny, return the larger.
1767     // If both are huge, avoid overflow by scaling by 1/sqrt(real.max/2).
1768     // If both are tiny, avoid underflow by scaling by sqrt(real.min_normal*real.epsilon).
1769 
1770     static immutable real SQRTMIN = 0x8.0p-8195L; // 0.5 * sqrt(real.min_normal); // This is a power of 2.
1771     static immutable real SQRTMAX = 1.0L / SQRTMIN; // 2^^((max_exp)/2) = nextUp(sqrt(real.max))
1772 
1773     static assert(2 * (SQRTMAX / 2) * (SQRTMAX / 2) <= real.max);
1774 
1775     // Proves that sqrt(real.max) ~~  0.5/sqrt(real.min_normal)
1776     static assert(real.min_normal * real.max > 2 && real.min_normal * real.max <= 4);
1777 
1778     real u = fabs(x);
1779     real v = fabs(y);
1780     if (!(u >= v))  // check for NaN as well.
1781     {
1782         v = u;
1783         u = fabs(y);
1784         if (u == real.infinity) return u; // hypot(inf, nan) == inf
1785         if (v == real.infinity) return v; // hypot(nan, inf) == inf
1786     }
1787 
1788     // Now u >= v, or else one is NaN.
1789     if (v >= SQRTMAX * 0.5)
1790     {
1791             // hypot(huge, huge) -- avoid overflow
1792         u *= SQRTMIN * 0.5;
1793         v *= SQRTMIN * 0.5;
1794         return sqrt(u * u + v * v) * SQRTMAX * 2.0;
1795     }
1796 
1797     if (u <= SQRTMIN)
1798     {
1799         // hypot (tiny, tiny) -- avoid underflow
1800         // This is only necessary to avoid setting the underflow
1801         // flag.
1802         u *= SQRTMAX / real.epsilon;
1803         v *= SQRTMAX / real.epsilon;
1804         return sqrt(u * u + v * v) * SQRTMIN * real.epsilon;
1805     }
1806 
1807     if (u * real.epsilon > v)
1808     {
1809         // hypot (huge, tiny) = huge
1810         return u;
1811     }
1812 
1813     // both are in the normal range
1814     return sqrt(u * u + v * v);
1815 }
1816 
1817 unittest
1818 {
1819     static real[3][] vals = // x,y,hypot
1820     [
1821         [ 0.0,   0.0,  0.0],
1822         [ 0.0,  -0.0,  0.0],
1823         [-0.0,  -0.0,  0.0],
1824         [ 3.0,   4.0,  5.0],
1825         [-300,  -400,  500],
1826         [ 0.0,   7.0,  7.0],
1827         [ 9.0,   9.0 * real.epsilon, 9.0],
1828         [ 0xb.0p+8188L /+88 / (64 * sqrt(real.min_normal))+/, 0xd.2p+8188L /+105 / (64 * sqrt(real.min_normal))+/, 0x8.9p+8189L /+137 / (64 * sqrt(real.min_normal))+/],
1829         [ 0xb.0p+8187L /+88 / (128 * sqrt(real.min_normal))+/, 0xd.2p+8187L /+105 / (128 * sqrt(real.min_normal))+/, 0x8.9p+8188L /+137 / (128 * sqrt(real.min_normal))+/],
1830         [ 3 * real.min_normal * real.epsilon, 4 * real.min_normal * real.epsilon, 5 * real.min_normal * real.epsilon],
1831         [ real.min_normal, real.min_normal,  0x1.6a09e667f3bcc908p-16382L /+sqrt(2.0L) * real,min_normal+/],
1832         [ real.max / 2.0, real.max / 2.0, 0x1.6a09e667f3bcc908p+16383L /+real.max / sqrt(2.0L)+/],
1833         [ 0x1.6a09e667f3bcc908p+16383L /+real.max / sqrt(2.0L)+/, 0x1.6a09e667f3bcc908p+16383L /+real.max / sqrt(2.0L)+/, real.max],
1834         [ real.max, 1.0, real.max],
1835         [ real.infinity, real.nan, real.infinity],
1836         [ real.nan, real.infinity, real.infinity],
1837         [ real.nan, real.nan, real.nan],
1838         [ real.nan, real.max, real.nan],
1839         [ real.max, real.nan, real.nan]
1840     ];
1841 
1842     for (int i = 0; i < vals.length; i++)
1843     {
1844         real x = vals[i][0];
1845         real y = vals[i][1];
1846         real z = vals[i][2];
1847         real h = hypot(x, y);
1848 
1849         test(isIdentical(z, h));
1850     }
1851     // NaN payloads
1852     test(isIdentical(hypot(NaN(0xABC), 3.14), NaN(0xABC)));
1853     test(isIdentical(hypot(7.6e39, NaN(0xABC)), NaN(0xABC)));
1854 }
1855 
1856 /***********************************
1857  * Evaluate polynomial A(x) = $(SUB a, 0) + $(SUB a, 1)x + $(SUB a, 2)$(POWER x,2)
1858  *                          + $(SUB a,3)$(POWER x,3); ...
1859  *
1860  * Uses Horner's rule A(x) = $(SUB a, 0) + x($(SUB a, 1) + x($(SUB a, 2)
1861  *                         + x($(SUB a, 3) + ...)))
1862  * Params:
1863  *      A =     array of coefficients $(SUB a, 0), $(SUB a, 1), etc.
1864  *      x =  point in which to evaluate polynomial
1865  */
1866 T poly(T)(T x, T[] A)
1867 {
1868     verify(A.length > 0);
1869 
1870   version (Naked_D_InlineAsm_X86) {
1871       static immutable bool Use_D_InlineAsm_X86 = true;
1872   } else static immutable bool Use_D_InlineAsm_X86 = false;
1873 
1874   // BUG (Inherited from Phobos): This code assumes a frame pointer in EBP.
1875   // This is not in the spec.
1876   static if (Use_D_InlineAsm_X86 && is(T==real) && T.sizeof == 10) {
1877     asm // assembler by W. Bright
1878     {
1879         // EDX = (A.length - 1) * real.sizeof
1880         mov     ECX,A[EBP]          ; // ECX = A.length
1881         dec     ECX                 ;
1882         lea     EDX,[ECX][ECX*8]    ;
1883         add     EDX,ECX             ;
1884         add     EDX,A+4[EBP]        ;
1885         fld     real ptr [EDX]      ; // ST0 = coeff[ECX]
1886         jecxz   return_ST           ;
1887         fld     x[EBP]              ; // ST0 = x
1888         fxch    ST(1)               ; // ST1 = x, ST0 = r
1889         align   4                   ;
1890     L2:  fmul    ST,ST(1)           ; // r *= x
1891         fld     real ptr -10[EDX]   ;
1892         sub     EDX,10              ; // deg--
1893         faddp   ST(1),ST            ;
1894         dec     ECX                 ;
1895         jne     L2                  ;
1896         fxch    ST(1)               ; // ST1 = r, ST0 = x
1897         fstp    ST(0)               ; // dump x
1898         align   4                   ;
1899     return_ST:                      ;
1900         ;
1901     }
1902   } else static if ( Use_D_InlineAsm_X86 && is(T==real) && T.sizeof==12){
1903     asm // assembler by W. Bright
1904     {
1905         // EDX = (A.length - 1) * real.sizeof
1906         mov     ECX,A[EBP]          ; // ECX = A.length
1907         dec     ECX                 ;
1908         lea     EDX,[ECX*8]         ;
1909         lea     EDX,[EDX][ECX*4]    ;
1910         add     EDX,A+4[EBP]        ;
1911         fld     real ptr [EDX]      ; // ST0 = coeff[ECX]
1912         jecxz   return_ST           ;
1913         fld     x                   ; // ST0 = x
1914         fxch    ST(1)               ; // ST1 = x, ST0 = r
1915         align   4                   ;
1916     L2: fmul    ST,ST(1)            ; // r *= x
1917         fld     real ptr -12[EDX]   ;
1918         sub     EDX,12              ; // deg--
1919         faddp   ST(1),ST            ;
1920         dec     ECX                 ;
1921         jne     L2                  ;
1922         fxch    ST(1)               ; // ST1 = r, ST0 = x
1923         fstp    ST(0)               ; // dump x
1924         align   4                   ;
1925     return_ST:                      ;
1926         ;
1927         }
1928   } else {
1929         ptrdiff_t i = A.length - 1;
1930         real r = A[i];
1931         while (--i >= 0)
1932         {
1933             r *= x;
1934             r += A[i];
1935         }
1936         return r;
1937   }
1938 }
1939 
1940 unittest
1941 {
1942     real x = 3.1;
1943     static immutable real[] pp = [56.1L, 32.7L, 6L];
1944 
1945     test( poly(x, pp) == (56.1L + (32.7L + 6L * x) * x) );
1946 
1947     test(isIdentical(poly(NaN(0xABC), pp), NaN(0xABC)));
1948 }
1949 
1950 package {
1951     T rationalPoly(T)(T x, T [] numerator, T [] denominator)
1952     {
1953         return poly(x, numerator)/poly(x, denominator);
1954     }
1955 }
1956 
1957 /*
1958  * Rounding (returning real)
1959  */
1960 
1961 /**
1962  * Returns the value of x rounded downward to the next integer
1963  * (toward negative infinity).
1964  */
1965 real floor(real x)
1966 {
1967     return core.stdc.math.floorl(x);
1968 }
1969 
1970 unittest {
1971     test(isIdentical(floor(NaN(0xABC)), NaN(0xABC)));
1972 }
1973 
1974 /**
1975  * Returns the value of x rounded upward to the next integer
1976  * (toward positive infinity).
1977  */
1978 real ceil(real x)
1979 {
1980     return core.stdc.math.ceill(x);
1981 }
1982 
1983 unittest {
1984     test(isIdentical(ceil(NaN(0xABC)), NaN(0xABC)));
1985 }
1986 
1987 /**
1988  * Return the value of x rounded to the nearest integer.
1989  * If the fractional part of x is exactly 0.5, the return value is rounded to
1990  * the even integer.
1991  */
1992 real round(real x)
1993 {
1994     return core.stdc.math.roundl(x);
1995 }
1996 
1997 unittest {
1998     test(isIdentical(round(NaN(0xABC)), NaN(0xABC)));
1999 }
2000 
2001 /**
2002  * Returns the integer portion of x, dropping the fractional portion.
2003  *
2004  * This is also known as "chop" rounding.
2005  */
2006 real trunc(real x)
2007 {
2008     return core.stdc.math.truncl(x);
2009 }
2010 
2011 unittest {
2012     test(isIdentical(trunc(NaN(0xABC)), NaN(0xABC)));
2013 }
2014 
2015 /**
2016 * Rounds x to the nearest int or long.
2017 *
2018 * This is generally the fastest method to convert a floating-point number
2019 * to an integer. Note that the results from this function
2020 * depend on the rounding mode, if the fractional part of x is exactly 0.5.
2021 * If using the default rounding mode (ties round to even integers)
2022 * rndint(4.5) == 4, rndint(5.5)==6.
2023 */
2024 int rndint(real x)
2025 {
2026     version(Naked_D_InlineAsm_X86)
2027     {
2028         int n;
2029         asm
2030         {
2031             fld x;
2032             fistp n;
2033         }
2034         return n;
2035     }
2036     else
2037     {
2038         return cast(int) core.stdc.math.lrintl(x);
2039     }
2040 }
2041 
2042 /** ditto */
2043 long rndlong(real x)
2044 {
2045     version(Naked_D_InlineAsm_X86)
2046     {
2047         long n;
2048         asm
2049         {
2050             fld x;
2051             fistp n;
2052         }
2053         return n;
2054     }
2055     else
2056     {
2057         return core.stdc.math.llrintl(x);
2058     }
2059 }
2060 
2061 version(D_InlineAsm_X86)
2062 {
2063     // Won't work for anything else yet
2064 
2065     unittest
2066     {
2067         int r = getIeeeRounding;
2068         test(r==RoundingMode.ROUNDTONEAREST);
2069         real b = 5.5;
2070         int cnear = ocean.math.Math.rndint(b);
2071         test(cnear == 6);
2072         auto oldrounding = setIeeeRounding(RoundingMode.ROUNDDOWN);
2073         scope (exit) setIeeeRounding(oldrounding);
2074 
2075         test(getIeeeRounding==RoundingMode.ROUNDDOWN);
2076 
2077         int cdown = ocean.math.Math.rndint(b);
2078         test(cdown==5);
2079     }
2080 
2081     unittest
2082     {
2083         // Check that the previous test correctly restored the rounding mode
2084         test(getIeeeRounding==RoundingMode.ROUNDTONEAREST);
2085     }
2086 }
2087 
2088 /***************************************************************************
2089 
2090     Integer pow function. Returns the power'th power of base
2091 
2092     Params:
2093         base  = base number
2094         power = power
2095 
2096     Returns:
2097         the power'th power of base
2098 
2099 ***************************************************************************/
2100 
2101 public ulong pow ( ulong base, ulong power )
2102 {
2103     ulong res = void;
2104 
2105     switch (power)
2106     {
2107         case 0:
2108             res = 1;
2109             break;
2110         case 1:
2111             res = base;
2112             break;
2113         case 2:
2114             res = base * base;
2115             break;
2116 
2117         default:
2118             res = 1;
2119 
2120             while (1)
2121             {
2122                 if (power & 1) res *= base;
2123                 power >>= 1;
2124                 if (!power) break;
2125                 base *= base;
2126             }
2127             break;
2128     }
2129 
2130     return res;
2131 }
2132 
2133 unittest
2134 {
2135     ulong x = 46;
2136 
2137     test(pow(x,0UL) == 1);
2138     test(pow(x,1UL) == x);
2139     test(pow(x,2UL) == x * x);
2140     test(pow(x,3UL) == x * x * x);
2141     test(pow(x,8UL) == (x * x) * (x * x) * (x * x) * (x * x));
2142 }
2143 
2144 /*******************************************************************************
2145 
2146     Does an integer division, rounding towards the nearest integer.
2147     Rounds to the even one if both integers are equal near.
2148 
2149     Params:
2150         a = number to divide
2151         b = number to divide by
2152 
2153     Returns:
2154         number divided according to given description
2155 
2156 *******************************************************************************/
2157 
2158 T divRoundEven(T)(T a, T b)
2159 {
2160     // both integers equal near?
2161     if (b % 2 == 0 && (a % b == b / 2 || a % b == -b / 2))
2162     {
2163         auto div_rounded_down = a / b;
2164 
2165         auto add = div_rounded_down < 0 ? -1 : 1;
2166 
2167         return div_rounded_down % 2 == 0 ?
2168             div_rounded_down : div_rounded_down + add;
2169     }
2170 
2171     if ( (a >= 0) != (b >= 0) )
2172     {
2173         return (a - b / 2) / b;
2174     }
2175     else
2176     {
2177         return (a + b / 2) / b;
2178     }
2179 }
2180 
2181 version (unittest)
2182 {
2183     import ocean.math.Math : rndlong;
2184 }
2185 
2186 unittest
2187 {
2188     long roundDivCheat ( long a, long b )
2189     {
2190         real x = cast(real)a / cast(real)b;
2191         return rndlong(x);
2192     }
2193 
2194     test(divRoundEven(-3, 2)  == -2);
2195     test(divRoundEven(3, 2)   == 2);
2196     test(divRoundEven(-3, -2) == 2);
2197     test(divRoundEven(3, -2)  == -2);
2198 
2199     test(divRoundEven(7, 11) == 1);
2200     test(divRoundEven(11, 11) == 1);
2201     test(divRoundEven(16, 11) == 1);
2202     test(divRoundEven(-17, 11) == -2);
2203     test(divRoundEven(-17, 11) == -2);
2204     test(divRoundEven(-16, 11) == -1);
2205 
2206     test(divRoundEven(17, -11) == -2);
2207     test(divRoundEven(16, -11) == -1);
2208     test(divRoundEven(-17, -11) == 2);
2209     test(divRoundEven(-16, -11) == 1);
2210 
2211     for (int i = -100; i <= 100; ++i) for (int j = -100; j <= 100; ++j)
2212     {
2213         if (j != 0)
2214         {
2215             test (divRoundEven(i,j) == roundDivCheat(i,j));
2216         }
2217     }
2218 }