std.math.exponential source code

1 // Written in the D programming language.
2 
3 /**
4 This is a submodule of $(MREF std, math).
5 
6 It contains several exponential and logarithm functions.
7 
8 Copyright: Copyright The D Language Foundation 2000 - 2011.
9            D implementations of exp, expm1, exp2, log, log10, log1p, and log2
10            functions are based on the CEPHES math library, which is Copyright
11            (C) 2001 Stephen L. Moshier $(LT)steve@moshier.net$(GT) and are
12            incorporated herein by permission of the author. The author reserves
13            the right to distribute this material elsewhere under different
14            copying permissions. These modifications are distributed here under
15            the following terms:
16 License:   $(HTTP www.boost.org/LICENSE_1_0.txt, Boost License 1.0).
17 Authors:   $(HTTP digitalmars.com, Walter Bright), Don Clugston,
18            Conversion of CEPHES math library to D by Iain Buclaw and David Nadlinger
19 Source: $(PHOBOSSRC std/math/exponential.d)
20 
21 Macros:
22     TABLE_SV = <table border="1" cellpadding="4" cellspacing="0">
23                <caption>Special Values</caption>
24                $0</table>
25     NAN = $(RED NAN)
26     PLUSMN = &plusmn;
27     INFIN = &infin;
28     PLUSMNINF = &plusmn;&infin;
29     LT = &lt;
30     GT = &gt;
31  */
32 
33 module std.math.exponential;
34 
35 import std.traits :  isFloatingPoint, isIntegral, isSigned, isUnsigned, Largest, Unqual;
36 
37 static import core.math;
38 static import core.stdc.math;
39 
40 version (LDC)
41 {
42     import ldc.intrinsics;
43 
44     version (CRuntime_Microsoft) version = LDC_MSVCRT;
45 
46     version (LDC_MSVCRT)   {}
47     else version (Android) {}
48     else
49     {
50         version (X86)    version = INLINE_YL2X;
51         version (X86_64) version = INLINE_YL2X;
52     }
53 }
54 
55 version (DigitalMars)
56 {
57     version (OSX) { }             // macOS 13 (M1) has issues emulating instruction
58     else version = INLINE_YL2X;   // x87 has opcodes for these
59 }
60 
61 version (D_InlineAsm_X86)    version = InlineAsm_X86_Any;
62 version (D_InlineAsm_X86_64) version = InlineAsm_X86_Any;
63 
64 version (LDC_MSVCRT)   {}
65 else version (Android) {}
66 else version (InlineAsm_X86_Any) version = InlineAsm_X87;
67 version (InlineAsm_X87)
68 {
69     static assert(real.mant_dig == 64);
70     version (CRuntime_Microsoft) version = InlineAsm_X87_MSVC;
71 }
72 
73 version (D_HardFloat)
74 {
75     // FloatingPointControl.clearExceptions() depends on version IeeeFlagsSupport
76     version (IeeeFlagsSupport) version = FloatingPointControlSupport;
77 }
78 
79 /**
80  * Compute the value of x $(SUPERSCRIPT n), where n is an integer
81  */
82 Unqual!F pow(F, G)(F x, G n) @nogc @trusted pure nothrow
83 if (isFloatingPoint!(F) && isIntegral!(G))
84 {
85   version (none)
86   {
87     // LDC: Leads to linking error on MSVC x64 as the intrinsic maps to
88     // MSVC++ function `pow(double/float, int)` (a C++ template for
89     // Visual Studio 2015).
90     // Most likely not worth the effort anyway (and hindering CTFE).
91     pragma(inline, true);
92     return llvm_powi!(Unqual!F)(x, cast(int) n);
93   }
94   else
95   {
96     import std.traits : Unsigned;
97 
98     real p = 1.0, v = void;
99     Unsigned!(Unqual!G) m = n;
100 
101     if (n < 0)
102     {
103         if (n == -1) return 1 / x;
104 
105         m = cast(typeof(m))(0 - n);
106         v = p / x;
107     }
108     else
109     {
110         switch (n)
111         {
112         case 0:
113             return 1.0;
114         case 1:
115             return x;
116         case 2:
117             return x * x;
118         default:
119         }
120 
121         v = x;
122     }
123 
124     while (1)
125     {
126         if (m & 1)
127             p *= v;
128         m >>= 1;
129         if (!m)
130             break;
131         v *= v;
132     }
133     return p;
134   } // !none
135 }
136 
137 ///
138 @safe pure nothrow @nogc unittest
139 {
140     import std.math.operations : feqrel;
141 
142     assert(pow(2.0, 5) == 32.0);
143     assert(pow(1.5, 9).feqrel(38.4433) > 16);
144     assert(pow(real.nan, 2) is real.nan);
145     assert(pow(real.infinity, 2) == real.infinity);
146 }
147 
148 @safe pure nothrow @nogc unittest
149 {
150     import std.math.operations : isClose, feqrel;
151 
152     // Make sure it instantiates and works properly on immutable values and
153     // with various integer and float types.
154     immutable real x = 46;
155     immutable float xf = x;
156     immutable double xd = x;
157     immutable uint one = 1;
158     immutable ushort two = 2;
159     immutable ubyte three = 3;
160     immutable ulong eight = 8;
161 
162     immutable int neg1 = -1;
163     immutable short neg2 = -2;
164     immutable byte neg3 = -3;
165     immutable long neg8 = -8;
166 
167 
168     assert(pow(x,0) == 1.0);
169     assert(pow(xd,one) == x);
170     assert(pow(xf,two) == x * x);
171     assert(pow(x,three) == x * x * x);
172     assert(pow(x,eight) == (x * x) * (x * x) * (x * x) * (x * x));
173 
174     assert(pow(x, neg1) == 1 / x);
175 
176     assert(isClose(pow(xd, neg2), cast(double) (1 / (x * x)), 1e-15));
177     assert(isClose(pow(xf, neg8), cast(float) (1 / ((x * x) * (x * x) * (x * x) * (x * x))), 1e-15));
178 
179     assert(feqrel(pow(x, neg3),  1 / (x * x * x)) >= real.mant_dig - 1);
180 }
181 
182 @safe @nogc nothrow unittest
183 {
184     import std.math.operations : isClose;
185 
186     assert(isClose(pow(2.0L, 10L), 1024, 1e-18));
187 }
188 
189 // https://issues.dlang.org/show_bug.cgi?id=21601
190 @safe @nogc nothrow pure unittest
191 {
192     // When reals are large enough the results of pow(b, e) can be
193     // calculated correctly, if b is of type float or double and e is
194     // not too large.
195     static if (real.mant_dig >= 64)
196     {
197         // expected result: 3.790e-42
198         assert(pow(-513645318757045764096.0f, -2) > 0.0);
199 
200         // expected result: 3.763915357831797e-309
201         assert(pow(-1.6299717435255677e+154, -2) > 0.0);
202     }
203 }
204 
205 @safe @nogc nothrow unittest
206 {
207     import std.math.operations : isClose;
208     import std.math.traits : isInfinity;
209 
210     static float f1 = 19100.0f;
211     static float f2 = 0.000012f;
212 
213     assert(isClose(pow(f1,9), 3.3829868e+38f));
214     assert(isInfinity(pow(f1,10)));
215     assert(pow(f2,9) > 0.0f);
216     assert(isClose(pow(f2,10), 0.0f, 0.0, float.min_normal));
217 
218     static double d1 = 21800.0;
219     static double d2 = 0.000012;
220 
221     assert(isClose(pow(d1,71), 1.0725339442974e+308));
222     assert(isInfinity(pow(d1,72)));
223     assert(pow(d2,65) > 0.0f);
224     assert(isClose(pow(d2,66), 0.0, 0.0, double.min_normal));
225 
226     static if (real.mant_dig == 64) // x87
227     {
228         static real r1 = 21950.0L;
229         static real r2 = 0.000011L;
230 
231         assert(isClose(pow(r1,1136), 7.4066175654969242752260330529e+4931L));
232         assert(isInfinity(pow(r1,1137)));
233         assert(pow(r2,998) > 0.0L);
234         assert(isClose(pow(r2,999), 0.0L, 0.0, real.min_normal));
235     }
236 }
237 
238 @safe @nogc nothrow pure unittest
239 {
240     import std.math.operations : isClose;
241 
242     enum f1 = 19100.0f;
243     enum f2 = 0.000012f;
244 
245     static assert(isClose(pow(f1,9), 3.3829868e+38f));
246     static assert(pow(f1,10) > float.max);
247     static assert(pow(f2,9) > 0.0f);
248     static assert(isClose(pow(f2,10), 0.0f, 0.0, float.min_normal));
249 
250     enum d1 = 21800.0;
251     enum d2 = 0.000012;
252 
253     static assert(isClose(pow(d1,71), 1.0725339442974e+308));
254     static assert(pow(d1,72) > double.max);
255     static assert(pow(d2,65) > 0.0f);
256     static assert(isClose(pow(d2,66), 0.0, 0.0, double.min_normal));
257 
258     static if (real.mant_dig == 64) // x87
259     {
260         enum r1 = 21950.0L;
261         enum r2 = 0.000011L;
262 
263         static assert(isClose(pow(r1,1136), 7.4066175654969242752260330529e+4931L));
264         static assert(pow(r1,1137) > real.max);
265         static assert(pow(r2,998) > 0.0L);
266         static assert(isClose(pow(r2,999), 0.0L, 0.0, real.min_normal));
267     }
268 }
269 
270 /**
271  * Compute the power of two integral numbers.
272  *
273  * Params:
274  *     x = base
275  *     n = exponent
276  *
277  * Returns:
278  *     x raised to the power of n. If n is negative the result is 1 / pow(x, -n),
279  *     which is calculated as integer division with remainder. This may result in
280  *     a division by zero error.
281  *
282  *     If both x and n are 0, the result is 1.
283  *
284  * Throws:
285  *     If x is 0 and n is negative, the result is the same as the result of a
286  *     division by zero.
287  */
288 typeof(Unqual!(F).init * Unqual!(G).init) pow(F, G)(F x, G n) @nogc @trusted pure nothrow
289 if (isIntegral!(F) && isIntegral!(G))
290 {
291     import std.traits : isSigned;
292 
293     typeof(return) p, v = void;
294     Unqual!G m = n;
295 
296     static if (isSigned!(F))
297     {
298         if (x == -1) return cast(typeof(return)) (m & 1 ? -1 : 1);
299     }
300     static if (isSigned!(G))
301     {
302         if (x == 0 && m <= -1) return x / 0;
303     }
304     if (x == 1) return 1;
305     static if (isSigned!(G))
306     {
307         if (m < 0) return 0;
308     }
309 
310     switch (m)
311     {
312     case 0:
313         p = 1;
314         break;
315 
316     case 1:
317         p = x;
318         break;
319 
320     case 2:
321         p = x * x;
322         break;
323 
324     default:
325         v = x;
326         p = 1;
327         while (1)
328         {
329             if (m & 1)
330                 p *= v;
331             m >>= 1;
332             if (!m)
333                 break;
334             v *= v;
335         }
336         break;
337     }
338     return p;
339 }
340 
341 ///
342 @safe pure nothrow @nogc unittest
343 {
344     assert(pow(2, 3) == 8);
345     assert(pow(3, 2) == 9);
346 
347     assert(pow(2, 10) == 1_024);
348     assert(pow(2, 20) == 1_048_576);
349     assert(pow(2, 30) == 1_073_741_824);
350 
351     assert(pow(0, 0) == 1);
352 
353     assert(pow(1, -5) == 1);
354     assert(pow(1, -6) == 1);
355     assert(pow(-1, -5) == -1);
356     assert(pow(-1, -6) == 1);
357 
358     assert(pow(-2, 5) == -32);
359     assert(pow(-2, -5) == 0);
360     assert(pow(cast(double) -2, -5) == -0.03125);
361 }
362 
363 @safe pure nothrow @nogc unittest
364 {
365     immutable int one = 1;
366     immutable byte two = 2;
367     immutable ubyte three = 3;
368     immutable short four = 4;
369     immutable long ten = 10;
370 
371     assert(pow(two, three) == 8);
372     assert(pow(two, ten) == 1024);
373     assert(pow(one, ten) == 1);
374     assert(pow(ten, four) == 10_000);
375     assert(pow(four, 10) == 1_048_576);
376     assert(pow(three, four) == 81);
377 }
378 
379 // https://issues.dlang.org/show_bug.cgi?id=7006
380 @safe pure nothrow @nogc unittest
381 {
382     assert(pow(5, -1) == 0);
383     assert(pow(-5, -1) == 0);
384     assert(pow(5, -2) == 0);
385     assert(pow(-5, -2) == 0);
386     assert(pow(-1, int.min) == 1);
387     assert(pow(-2, int.min) == 0);
388 
389     assert(pow(4294967290UL,2) == 18446744022169944100UL);
390     assert(pow(0,uint.max) == 0);
391 }
392 
393 /**Computes integer to floating point powers.*/
394 real pow(I, F)(I x, F y) @nogc @trusted pure nothrow
395 if (isIntegral!I && isFloatingPoint!F)
396 {
397     return pow(cast(real) x, cast(Unqual!F) y);
398 }
399 
400 ///
401 @safe pure nothrow @nogc unittest
402 {
403     assert(pow(2, 5.0) == 32.0);
404     assert(pow(7, 3.0) == 343.0);
405     assert(pow(2, real.nan) is real.nan);
406     assert(pow(2, real.infinity) == real.infinity);
407 }
408 
409 /**
410  * Calculates x$(SUPERSCRIPT y).
411  *
412  * $(TABLE_SV
413  * $(TR $(TH x) $(TH y) $(TH pow(x, y))
414  *      $(TH div 0) $(TH invalid?))
415  * $(TR $(TD anything)      $(TD $(PLUSMN)0.0)                $(TD 1.0)
416  *      $(TD no)        $(TD no) )
417  * $(TR $(TD |x| $(GT) 1)    $(TD +$(INFIN))                  $(TD +$(INFIN))
418  *      $(TD no)        $(TD no) )
419  * $(TR $(TD |x| $(LT) 1)    $(TD +$(INFIN))                  $(TD +0.0)
420  *      $(TD no)        $(TD no) )
421  * $(TR $(TD |x| $(GT) 1)    $(TD -$(INFIN))                  $(TD +0.0)
422  *      $(TD no)        $(TD no) )
423  * $(TR $(TD |x| $(LT) 1)    $(TD -$(INFIN))                  $(TD +$(INFIN))
424  *      $(TD no)        $(TD no) )
425  * $(TR $(TD +$(INFIN))      $(TD $(GT) 0.0)                  $(TD +$(INFIN))
426  *      $(TD no)        $(TD no) )
427  * $(TR $(TD +$(INFIN))      $(TD $(LT) 0.0)                  $(TD +0.0)
428  *      $(TD no)        $(TD no) )
429  * $(TR $(TD -$(INFIN))      $(TD odd integer $(GT) 0.0)      $(TD -$(INFIN))
430  *      $(TD no)        $(TD no) )
431  * $(TR $(TD -$(INFIN))      $(TD $(GT) 0.0, not odd integer) $(TD +$(INFIN))
432  *      $(TD no)        $(TD no))
433  * $(TR $(TD -$(INFIN))      $(TD odd integer $(LT) 0.0)      $(TD -0.0)
434  *      $(TD no)        $(TD no) )
435  * $(TR $(TD -$(INFIN))      $(TD $(LT) 0.0, not odd integer) $(TD +0.0)
436  *      $(TD no)        $(TD no) )
437  * $(TR $(TD $(PLUSMN)1.0)   $(TD $(PLUSMN)$(INFIN))          $(TD -$(NAN))
438  *      $(TD no)        $(TD yes) )
439  * $(TR $(TD $(LT) 0.0)      $(TD finite, nonintegral)        $(TD $(NAN))
440  *      $(TD no)        $(TD yes))
441  * $(TR $(TD $(PLUSMN)0.0)   $(TD odd integer $(LT) 0.0)      $(TD $(PLUSMNINF))
442  *      $(TD yes)       $(TD no) )
443  * $(TR $(TD $(PLUSMN)0.0)   $(TD $(LT) 0.0, not odd integer) $(TD +$(INFIN))
444  *      $(TD yes)       $(TD no))
445  * $(TR $(TD $(PLUSMN)0.0)   $(TD odd integer $(GT) 0.0)      $(TD $(PLUSMN)0.0)
446  *      $(TD no)        $(TD no) )
447  * $(TR $(TD $(PLUSMN)0.0)   $(TD $(GT) 0.0, not odd integer) $(TD +0.0)
448  *      $(TD no)        $(TD no) )
449  * )
450  */
451 Unqual!(Largest!(F, G)) pow(F, G)(F x, G y) @nogc @trusted pure nothrow
452 if (isFloatingPoint!(F) && isFloatingPoint!(G))
453 {
454     import core.math : fabs, sqrt;
455     import std.math.traits : isInfinity, isNaN, signbit;
456 
457     alias Float = typeof(return);
458 
459   version (none) // LDC FIXME: Use of this LLVM intrinsic causes a unit test failure
460   {
461     pragma(inline, true);
462     return llvm_pow!(Float)(x, y);
463   }
464   else
465   {
466     static real impl(real x, real y) @nogc pure nothrow
467     {
468         // Special cases.
469         if (isNaN(y))
470             return y;
471         if (isNaN(x) && y != 0.0)
472             return x;
473 
474         // Even if x is NaN.
475         if (y == 0.0)
476             return 1.0;
477         if (y == 1.0)
478             return x;
479 
480         if (isInfinity(y))
481         {
482             if (isInfinity(x))
483             {
484                 if (!signbit(y) && !signbit(x))
485                     return F.infinity;
486                 else
487                     return F.nan;
488             }
489             else if (fabs(x) > 1)
490             {
491                 if (signbit(y))
492                     return +0.0;
493                 else
494                     return F.infinity;
495             }
496             else if (fabs(x) == 1)
497             {
498                 return F.nan;
499             }
500             else // < 1
501             {
502                 if (signbit(y))
503                     return F.infinity;
504                 else
505                     return +0.0;
506             }
507         }
508         if (isInfinity(x))
509         {
510             if (signbit(x))
511             {
512                 long i = cast(long) y;
513                 if (y > 0.0)
514                 {
515                     if (i == y && i & 1)
516                         return -F.infinity;
517                     else if (i == y)
518                         return F.infinity;
519                     else
520                         return -F.nan;
521                 }
522                 else if (y < 0.0)
523                 {
524                     if (i == y && i & 1)
525                         return -0.0;
526                     else if (i == y)
527                         return +0.0;
528                     else
529                         return F.nan;
530                 }
531             }
532             else
533             {
534                 if (y > 0.0)
535                     return F.infinity;
536                 else if (y < 0.0)
537                     return +0.0;
538             }
539         }
540 
541         if (x == 0.0)
542         {
543             if (signbit(x))
544             {
545                 long i = cast(long) y;
546                 if (y > 0.0)
547                 {
548                     if (i == y && i & 1)
549                         return -0.0;
550                     else
551                         return +0.0;
552                 }
553                 else if (y < 0.0)
554                 {
555                     if (i == y && i & 1)
556                         return -F.infinity;
557                     else
558                         return F.infinity;
559                 }
560             }
561             else
562             {
563                 if (y > 0.0)
564                     return +0.0;
565                 else if (y < 0.0)
566                     return F.infinity;
567             }
568         }
569         if (x == 1.0)
570             return 1.0;
571 
572         if (y >= F.max)
573         {
574             if ((x > 0.0 && x < 1.0) || (x > -1.0 && x < 0.0))
575                 return 0.0;
576             if (x > 1.0 || x < -1.0)
577                 return F.infinity;
578         }
579         if (y <= -F.max)
580         {
581             if ((x > 0.0 && x < 1.0) || (x > -1.0 && x < 0.0))
582                 return F.infinity;
583             if (x > 1.0 || x < -1.0)
584                 return 0.0;
585         }
586 
587         if (x >= F.max)
588         {
589             if (y > 0.0)
590                 return F.infinity;
591             else
592                 return 0.0;
593         }
594         if (x <= -F.max)
595         {
596             long i = cast(long) y;
597             if (y > 0.0)
598             {
599                 if (i == y && i & 1)
600                     return -F.infinity;
601                 else
602                     return F.infinity;
603             }
604             else if (y < 0.0)
605             {
606                 if (i == y && i & 1)
607                     return -0.0;
608                 else
609                     return +0.0;
610             }
611         }
612 
613         // Integer power of x.
614         long iy = cast(long) y;
615         if (iy == y && fabs(y) < 32_768.0)
616             return pow(x, iy);
617 
618         real sign = 1.0;
619         if (x < 0)
620         {
621             // Result is real only if y is an integer
622             // Check for a non-zero fractional part
623             enum maxOdd = pow(2.0L, real.mant_dig) - 1.0L;
624             static if (maxOdd > ulong.max)
625             {
626                 // Generic method, for any FP type
627                 import std.math.rounding : floor;
628                 if (floor(y) != y)
629                     return sqrt(x); // Complex result -- create a NaN
630 
631                 const hy = 0.5 * y;
632                 if (floor(hy) != hy)
633                     sign = -1.0;
634             }
635             else
636             {
637                 // Much faster, if ulong has enough precision
638                 const absY = fabs(y);
639                 if (absY <= maxOdd)
640                 {
641                     const uy = cast(ulong) absY;
642                     if (uy != absY)
643                         return sqrt(x); // Complex result -- create a NaN
644 
645                     if (uy & 1)
646                         sign = -1.0;
647                 }
648             }
649             x = -x;
650         }
651         version (INLINE_YL2X)
652         {
653             // If x > 0, x ^^ y == 2 ^^ ( y * log2(x) )
654             // TODO: This is not accurate in practice. A fast and accurate
655             // (though complicated) method is described in:
656             // "An efficient rounding boundary test for pow(x, y)
657             // in double precision", C.Q. Lauter and V. Lefèvre, INRIA (2007).
658             return sign * exp2( core.math.yl2x(x, y) );
659         }
660         else
661         {
662             // If x > 0, x ^^ y == 2 ^^ ( y * log2(x) )
663             // TODO: This is not accurate in practice. A fast and accurate
664             // (though complicated) method is described in:
665             // "An efficient rounding boundary test for pow(x, y)
666             // in double precision", C.Q. Lauter and V. Lefèvre, INRIA (2007).
667             Float w = exp2(y * log2(x));
668             return sign * w;
669         }
670     }
671     return impl(x, y);
672   } // !none
673 }
674 
675 ///
676 @safe pure nothrow @nogc unittest
677 {
678     import std.math.operations : isClose;
679 
680     assert(isClose(pow(2.0, 3.0), 8.0));
681     assert(isClose(pow(1.5, 10.0), 57.6650390625));
682 
683     // square root of 9
684     assert(isClose(pow(9.0, 0.5), 3.0));
685     // 10th root of 1024
686     assert(isClose(pow(1024.0, 0.1), 2.0));
687 
688     assert(isClose(pow(-4.0, 3.0), -64.0));
689 
690     // reciprocal of 4 ^^ 2
691     assert(isClose(pow(4.0, -2.0), 0.0625));
692     // reciprocal of (-2) ^^ 3
693     assert(isClose(pow(-2.0, -3.0), -0.125));
694 
695     assert(isClose(pow(-2.5, 3.0), -15.625));
696     // reciprocal of 2.5 ^^ 3
697     assert(isClose(pow(2.5, -3.0), 0.064));
698     // reciprocal of (-2.5) ^^ 3
699     assert(isClose(pow(-2.5, -3.0), -0.064));
700 
701     // reciprocal of square root of 4
702     assert(isClose(pow(4.0, -0.5), 0.5));
703 
704     // per definition
705     assert(isClose(pow(0.0, 0.0), 1.0));
706 }
707 
708 ///
709 @safe pure nothrow @nogc unittest
710 {
711     import std.math.operations : isClose;
712 
713     // the result is a complex number
714     // which cannot be represented as floating point number
715     import std.math.traits : isNaN;
716     assert(isNaN(pow(-2.5, -1.5)));
717 
718     // use the ^^-operator of std.complex instead
719     import std.complex : complex;
720     auto c1 = complex(-2.5, 0.0);
721     auto c2 = complex(-1.5, 0.0);
722     auto result = c1 ^^ c2;
723     // exact result apparently depends on `real` precision => increased tolerance
724     assert(isClose(result.re, -4.64705438e-17, 2e-4));
725     assert(isClose(result.im, 2.52982e-1, 2e-4));
726 }
727 
728 @safe pure nothrow @nogc unittest
729 {
730     import std.math.traits : isNaN;
731 
732     assert(pow(1.5, real.infinity) == real.infinity);
733     assert(pow(0.5, real.infinity) == 0.0);
734     assert(pow(1.5, -real.infinity) == 0.0);
735     assert(pow(0.5, -real.infinity) == real.infinity);
736     assert(pow(real.infinity, 1.0) == real.infinity);
737     assert(pow(real.infinity, -1.0) == 0.0);
738     assert(pow(real.infinity, real.infinity) == real.infinity);
739     assert(pow(-real.infinity, 1.0) == -real.infinity);
740     assert(pow(-real.infinity, 2.0) == real.infinity);
741     assert(pow(-real.infinity, -1.0) == -0.0);
742     assert(pow(-real.infinity, -2.0) == 0.0);
743     assert(isNaN(pow(1.0, real.infinity)));
744     assert(pow(0.0, -1.0) == real.infinity);
745     assert(pow(real.nan, 0.0) == 1.0);
746     assert(isNaN(pow(real.nan, 3.0)));
747     assert(isNaN(pow(3.0, real.nan)));
748 }
749 
750 @safe @nogc nothrow unittest
751 {
752     import std.math.operations : isClose;
753 
754     assert(isClose(pow(2.0L, 10.0L), 1024, 1e-18));
755 }
756 
757 @safe pure nothrow @nogc unittest
758 {
759     import std.math.operations : isClose;
760     import std.math.traits : isIdentical, isNaN;
761     import std.math.constants : PI;
762 
763     // Test all the special values.  These unittests can be run on Windows
764     // by temporarily changing the version (linux) to version (all).
765     immutable float zero = 0;
766     immutable real one = 1;
767     immutable double two = 2;
768     immutable float three = 3;
769     immutable float fnan = float.nan;
770     immutable double dnan = double.nan;
771     immutable real rnan = real.nan;
772     immutable dinf = double.infinity;
773     immutable rninf = -real.infinity;
774 
775     assert(pow(fnan, zero) == 1);
776     assert(pow(dnan, zero) == 1);
777     assert(pow(rnan, zero) == 1);
778 
779     assert(pow(two, dinf) == double.infinity);
780     assert(isIdentical(pow(0.2f, dinf), +0.0));
781     assert(pow(0.99999999L, rninf) == real.infinity);
782     assert(isIdentical(pow(1.000000001, rninf), +0.0));
783     assert(pow(dinf, 0.001) == dinf);
784     assert(isIdentical(pow(dinf, -0.001), +0.0));
785     assert(pow(rninf, 3.0L) == rninf);
786     assert(pow(rninf, 2.0L) == real.infinity);
787     assert(isIdentical(pow(rninf, -3.0), -0.0));
788     assert(isIdentical(pow(rninf, -2.0), +0.0));
789 
790     // @@@BUG@@@ somewhere
791     version (OSX) {} else assert(isNaN(pow(one, dinf)));
792     version (OSX) {} else assert(isNaN(pow(-one, dinf)));
793     assert(isNaN(pow(-0.2, PI)));
794     // boundary cases. Note that epsilon == 2^^-n for some n,
795     // so 1/epsilon == 2^^n is always even.
796     assert(pow(-1.0L, 1/real.epsilon - 1.0L) == -1.0L);
797     static if (LLVM_version >= 1300) { /* LDC: on x86, yields -1 with enabled optimizations */ } else
798         assert(pow(-1.0L, 1/real.epsilon) == 1.0L);
799     assert(isNaN(pow(-1.0L, 1/real.epsilon-0.5L)));
800     assert(isNaN(pow(-1.0L, -1/real.epsilon+0.5L)));
801 
802     assert(pow(0.0, -3.0) == double.infinity);
803     assert(pow(-0.0, -3.0) == -double.infinity);
804     assert(pow(0.0, -PI) == double.infinity);
805     assert(pow(-0.0, -PI) == double.infinity);
806     assert(isIdentical(pow(0.0, 5.0), 0.0));
807     assert(isIdentical(pow(-0.0, 5.0), -0.0));
808     assert(isIdentical(pow(0.0, 6.0), 0.0));
809     assert(isIdentical(pow(-0.0, 6.0), 0.0));
810 
811     // https://issues.dlang.org/show_bug.cgi?id=14786 fixed
812     immutable real maxOdd = pow(2.0L, real.mant_dig) - 1.0L;
813     assert(pow(-1.0L,  maxOdd) == -1.0L);
814     assert(pow(-1.0L, -maxOdd) == -1.0L);
815     assert(pow(-1.0L, maxOdd + 1.0L) == 1.0L);
816     assert(pow(-1.0L, -maxOdd + 1.0L) == 1.0L);
817     assert(pow(-1.0L, maxOdd - 1.0L) == 1.0L);
818     assert(pow(-1.0L, -maxOdd - 1.0L) == 1.0L);
819 
820     // Now, actual numbers.
821     assert(isClose(pow(two, three), 8.0));
822     assert(isClose(pow(two, -2.5), 0.1767766953));
823 
824     // Test integer to float power.
825     immutable uint twoI = 2;
826     assert(isClose(pow(twoI, three), 8.0));
827 }
828 
829 // https://issues.dlang.org/show_bug.cgi?id=20508
830 @safe pure nothrow @nogc unittest
831 {
832     import std.math.traits : isNaN;
833 
834     assert(isNaN(pow(-double.infinity, 0.5)));
835 
836     assert(isNaN(pow(-real.infinity, real.infinity)));
837     assert(isNaN(pow(-real.infinity, -real.infinity)));
838     assert(isNaN(pow(-real.infinity, 1.234)));
839     assert(isNaN(pow(-real.infinity, -0.751)));
840     assert(pow(-real.infinity, 0.0) == 1.0);
841 }
842 
843 /** Computes the value of a positive integer `x`, raised to the power `n`, modulo `m`.
844  *
845  *  Params:
846  *      x = base
847  *      n = exponent
848  *      m = modulus
849  *
850  *  Returns:
851  *      `x` to the power `n`, modulo `m`.
852  *      The return type is the largest of `x`'s and `m`'s type.
853  *
854  * The function requires that all values have unsigned types.
855  */
856 Unqual!(Largest!(F, H)) powmod(F, G, H)(F x, G n, H m)
857 if (isUnsigned!F && isUnsigned!G && isUnsigned!H)
858 {
859     import std.meta : AliasSeq;
860 
861     alias T = Unqual!(Largest!(F, H));
862     static if (T.sizeof <= 4)
863     {
864         alias DoubleT = AliasSeq!(void, ushort, uint, void, ulong)[T.sizeof];
865     }
866 
867     static T mulmod(T a, T b, T c)
868     {
869         static if (T.sizeof == 8)
870         {
871             static T addmod(T a, T b, T c)
872             {
873                 b = c - b;
874                 if (a >= b)
875                     return a - b;
876                 else
877                     return c - b + a;
878             }
879 
880             T result = 0, tmp;
881 
882             b %= c;
883             while (a > 0)
884             {
885                 if (a & 1)
886                     result = addmod(result, b, c);
887 
888                 a >>= 1;
889                 b = addmod(b, b, c);
890             }
891 
892             return result;
893         }
894         else
895         {
896             DoubleT result = cast(DoubleT) (cast(DoubleT) a * cast(DoubleT) b);
897             return result % c;
898         }
899     }
900 
901     T base = x, result = 1, modulus = m;
902     Unqual!G exponent = n;
903 
904     while (exponent > 0)
905     {
906         if (exponent & 1)
907             result = mulmod(result, base, modulus);
908 
909         base = mulmod(base, base, modulus);
910         exponent >>= 1;
911     }
912 
913     return result;
914 }
915 
916 ///
917 @safe pure nothrow @nogc unittest
918 {
919     assert(powmod(1U, 10U, 3U) == 1);
920     assert(powmod(3U, 2U, 6U) == 3);
921     assert(powmod(5U, 5U, 15U) == 5);
922     assert(powmod(2U, 3U, 5U) == 3);
923     assert(powmod(2U, 4U, 5U) == 1);
924     assert(powmod(2U, 5U, 5U) == 2);
925 }
926 
927 @safe pure nothrow @nogc unittest
928 {
929     ulong a = 18446744073709551615u, b = 20u, c = 18446744073709551610u;
930     assert(powmod(a, b, c) == 95367431640625u);
931     a = 100; b = 7919; c = 18446744073709551557u;
932     assert(powmod(a, b, c) == 18223853583554725198u);
933     a = 117; b = 7919; c = 18446744073709551557u;
934     assert(powmod(a, b, c) == 11493139548346411394u);
935     a = 134; b = 7919; c = 18446744073709551557u;
936     assert(powmod(a, b, c) == 10979163786734356774u);
937     a = 151; b = 7919; c = 18446744073709551557u;
938     assert(powmod(a, b, c) == 7023018419737782840u);
939     a = 168; b = 7919; c = 18446744073709551557u;
940     assert(powmod(a, b, c) == 58082701842386811u);
941     a = 185; b = 7919; c = 18446744073709551557u;
942     assert(powmod(a, b, c) == 17423478386299876798u);
943     a = 202; b = 7919; c = 18446744073709551557u;
944     assert(powmod(a, b, c) == 5522733478579799075u);
945     a = 219; b = 7919; c = 18446744073709551557u;
946     assert(powmod(a, b, c) == 15230218982491623487u);
947     a = 236; b = 7919; c = 18446744073709551557u;
948     assert(powmod(a, b, c) == 5198328724976436000u);
949 
950     a = 0; b = 7919; c = 18446744073709551557u;
951     assert(powmod(a, b, c) == 0);
952     a = 123; b = 0; c = 18446744073709551557u;
953     assert(powmod(a, b, c) == 1);
954 
955     immutable ulong a1 = 253, b1 = 7919, c1 = 18446744073709551557u;
956     assert(powmod(a1, b1, c1) == 3883707345459248860u);
957 
958     uint x = 100 ,y = 7919, z = 1844674407u;
959     assert(powmod(x, y, z) == 1613100340u);
960     x = 134; y = 7919; z = 1844674407u;
961     assert(powmod(x, y, z) == 734956622u);
962     x = 151; y = 7919; z = 1844674407u;
963     assert(powmod(x, y, z) == 1738696945u);
964     x = 168; y = 7919; z = 1844674407u;
965     assert(powmod(x, y, z) == 1247580927u);
966     x = 185; y = 7919; z = 1844674407u;
967     assert(powmod(x, y, z) == 1293855176u);
968     x = 202; y = 7919; z = 1844674407u;
969     assert(powmod(x, y, z) == 1566963682u);
970     x = 219; y = 7919; z = 1844674407u;
971     assert(powmod(x, y, z) == 181227807u);
972     x = 236; y = 7919; z = 1844674407u;
973     assert(powmod(x, y, z) == 217988321u);
974     x = 253; y = 7919; z = 1844674407u;
975     assert(powmod(x, y, z) == 1588843243u);
976 
977     x = 0; y = 7919; z = 184467u;
978     assert(powmod(x, y, z) == 0);
979     x = 123; y = 0; z = 1844674u;
980     assert(powmod(x, y, z) == 1);
981 
982     immutable ubyte x1 = 117;
983     immutable uint y1 = 7919;
984     immutable uint z1 = 1844674407u;
985     auto res = powmod(x1, y1, z1);
986     assert(is(typeof(res) == uint));
987     assert(res == 9479781u);
988 
989     immutable ushort x2 = 123;
990     immutable uint y2 = 203;
991     immutable ubyte z2 = 113;
992     auto res2 = powmod(x2, y2, z2);
993     assert(is(typeof(res2) == ushort));
994     assert(res2 == 42u);
995 }
996 
997 /**
998  * Calculates e$(SUPERSCRIPT x).
999  *
1000  *  $(TABLE_SV
1001  *    $(TR $(TH x)             $(TH e$(SUPERSCRIPT x)) )
1002  *    $(TR $(TD +$(INFIN))     $(TD +$(INFIN)) )
1003  *    $(TR $(TD -$(INFIN))     $(TD +0.0)      )
1004  *    $(TR $(TD $(NAN))        $(TD $(NAN))    )
1005  *  )
1006  */
1007 version (none) // LDC FIXME: Use of this LLVM intrinsic causes a unit test failure
1008 {
1009     pragma(inline, true):
1010     real   exp(real   x) @safe pure nothrow @nogc { return llvm_exp(x); }
1011     ///ditto
1012     double exp(double x) @safe pure nothrow @nogc { return llvm_exp(x); }
1013     ///ditto
1014     float  exp(float  x) @safe pure nothrow @nogc { return llvm_exp(x); }
1015 }
1016 else
1017 {
1018 
1019 pragma(inline, true)
1020 real exp(real x) @trusted pure nothrow @nogc // TODO: @safe
1021 {
1022     version (InlineAsm_X87)
1023     {
1024         import std.math.constants : LOG2E;
1025 
1026         //  e^^x = 2^^(LOG2E*x)
1027         // (This is valid because the overflow & underflow limits for exp
1028         // and exp2 are so similar).
1029         if (!__ctfe)
1030             return exp2Asm(LOG2E*x);
1031     }
1032     return expImpl(x);
1033 }
1034 
1035 /// ditto
1036 pragma(inline, true)
1037 double exp(double x) @safe pure nothrow @nogc { return __ctfe ? cast(double) exp(cast(real) x) : expImpl(x); }
1038 
1039 /// ditto
1040 pragma(inline, true)
1041 float exp(float x) @safe pure nothrow @nogc { return __ctfe ? cast(float) exp(cast(real) x) : expImpl(x); }
1042 
1043 } // !none
1044 
1045 ///
1046 @safe unittest
1047 {
1048     import std.math.operations : feqrel;
1049     import std.math.constants : E;
1050 
1051     assert(exp(0.0) == 1.0);
1052     assert(exp(3.0).feqrel(E * E * E) > 16);
1053 }
1054 
1055 private T expImpl(T)(T x) @safe pure nothrow @nogc
1056 {
1057     import std.math : floatTraits, RealFormat;
1058     import std.math.traits : isNaN;
1059     import std.math.rounding : floor;
1060     import std.math.algebraic : poly;
1061     import std.math.constants : LOG2E;
1062 
1063     alias F = floatTraits!T;
1064     static if (F.realFormat == RealFormat.ieeeSingle)
1065     {
1066         static immutable T[6] P = [
1067             5.0000001201E-1,
1068             1.6666665459E-1,
1069             4.1665795894E-2,
1070             8.3334519073E-3,
1071             1.3981999507E-3,
1072             1.9875691500E-4,
1073         ];
1074 
1075         enum T C1 = 0.693359375;
1076         enum T C2 = -2.12194440e-4;
1077 
1078         // Overflow and Underflow limits.
1079         enum T OF = 88.72283905206835;
1080         enum T UF = -103.278929903431851103; // ln(2^-149)
1081     }
1082     else static if (F.realFormat == RealFormat.ieeeDouble)
1083     {
1084         // Coefficients for exp(x)
1085         static immutable T[3] P = [
1086             9.99999999999999999910E-1L,
1087             3.02994407707441961300E-2L,
1088             1.26177193074810590878E-4L,
1089         ];
1090         static immutable T[4] Q = [
1091             2.00000000000000000009E0L,
1092             2.27265548208155028766E-1L,
1093             2.52448340349684104192E-3L,
1094             3.00198505138664455042E-6L,
1095         ];
1096 
1097         // C1 + C2 = LN2.
1098         enum T C1 = 6.93145751953125E-1;
1099         enum T C2 = 1.42860682030941723212E-6;
1100 
1101         // Overflow and Underflow limits.
1102         enum T OF =  7.09782712893383996732E2;  // ln((1-2^-53) * 2^1024)
1103         enum T UF = -7.451332191019412076235E2; // ln(2^-1075)
1104     }
1105     else static if (F.realFormat == RealFormat.ieeeExtended ||
1106                     F.realFormat == RealFormat.ieeeExtended53)
1107     {
1108         // Coefficients for exp(x)
1109         static immutable T[3] P = [
1110             9.9999999999999999991025E-1L,
1111             3.0299440770744196129956E-2L,
1112             1.2617719307481059087798E-4L,
1113         ];
1114         static immutable T[4] Q = [
1115             2.0000000000000000000897E0L,
1116             2.2726554820815502876593E-1L,
1117             2.5244834034968410419224E-3L,
1118             3.0019850513866445504159E-6L,
1119         ];
1120 
1121         // C1 + C2 = LN2.
1122         enum T C1 = 6.9314575195312500000000E-1L;
1123         enum T C2 = 1.4286068203094172321215E-6L;
1124 
1125         // Overflow and Underflow limits.
1126         enum T OF =  1.1356523406294143949492E4L;  // ln((1-2^-64) * 2^16384)
1127         enum T UF = -1.13994985314888605586758E4L; // ln(2^-16446)
1128     }
1129     else static if (F.realFormat == RealFormat.ieeeQuadruple)
1130     {
1131         // Coefficients for exp(x) - 1
1132         static immutable T[5] P = [
1133             9.999999999999999999999999999999999998502E-1L,
1134             3.508710990737834361215404761139478627390E-2L,
1135             2.708775201978218837374512615596512792224E-4L,
1136             6.141506007208645008909088812338454698548E-7L,
1137             3.279723985560247033712687707263393506266E-10L
1138         ];
1139         static immutable T[6] Q = [
1140             2.000000000000000000000000000000000000150E0,
1141             2.368408864814233538909747618894558968880E-1L,
1142             3.611828913847589925056132680618007270344E-3L,
1143             1.504792651814944826817779302637284053660E-5L,
1144             1.771372078166251484503904874657985291164E-8L,
1145             2.980756652081995192255342779918052538681E-12L
1146         ];
1147 
1148         // C1 + C2 = LN2.
1149         enum T C1 = 6.93145751953125E-1L;
1150         enum T C2 = 1.428606820309417232121458176568075500134E-6L;
1151 
1152         // Overflow and Underflow limits.
1153         enum T OF =  1.135583025911358400418251384584930671458833e4L;
1154         enum T UF = -1.143276959615573793352782661133116431383730e4L;
1155     }
1156     else
1157         static assert(0, "Not implemented for this architecture");
1158 
1159     // Special cases.
1160     if (isNaN(x))
1161         return x;
1162     if (x > OF)
1163         return real.infinity;
1164     if (x < UF)
1165         return 0.0;
1166 
1167     // Express: e^^x = e^^g * 2^^n
1168     //   = e^^g * e^^(n * LOG2E)
1169     //   = e^^(g + n * LOG2E)
1170     T xx = floor((cast(T) LOG2E) * x + cast(T) 0.5);
1171     const int n = cast(int) xx;
1172     x -= xx * C1;
1173     x -= xx * C2;
1174 
1175     static if (F.realFormat == RealFormat.ieeeSingle)
1176     {
1177         xx = x * x;
1178         x = poly(x, P) * xx + x + 1.0f;
1179     }
1180     else
1181     {
1182         // Rational approximation for exponential of the fractional part:
1183         //  e^^x = 1 + 2x P(x^^2) / (Q(x^^2) - P(x^^2))
1184         xx = x * x;
1185         const T px = x * poly(xx, P);
1186         x = px / (poly(xx, Q) - px);
1187         x = (cast(T) 1.0) + (cast(T) 2.0) * x;
1188     }
1189 
1190     // Scale by power of 2.
1191     x = core.math.ldexp(x, n);
1192 
1193     return x;
1194 }
1195 
1196 @safe @nogc nothrow unittest
1197 {
1198     import std.math : floatTraits, RealFormat;
1199     import std.math.operations : NaN, feqrel, isClose;
1200     import std.math.constants : E;
1201     import std.math.traits : isIdentical;
1202     import std.math.algebraic : abs;
1203 
1204     version (IeeeFlagsSupport) import std.math.hardware : IeeeFlags, resetIeeeFlags, ieeeFlags;
1205     version (FloatingPointControlSupport)
1206     {
1207         import std.math.hardware : FloatingPointControl;
1208 
1209         FloatingPointControl ctrl;
1210         if (FloatingPointControl.hasExceptionTraps)
1211             ctrl.disableExceptions(FloatingPointControl.allExceptions);
1212         ctrl.rounding = FloatingPointControl.roundToNearest;
1213     }
1214 
1215     static void testExp(T)()
1216     {
1217         enum realFormat = floatTraits!T.realFormat;
1218         static if (realFormat == RealFormat.ieeeQuadruple)
1219         {
1220             static immutable T[2][] exptestpoints =
1221             [ //  x               exp(x)
1222                 [ 1.0L,           E                                        ],
1223                 [ 0.5L,           0x1.a61298e1e069bc972dfefab6df34p+0L     ],
1224                 [ 3.0L,           E*E*E                                    ],
1225                 [ 0x1.6p+13L,     0x1.6e509d45728655cdb4840542acb5p+16250L ], // near overflow
1226                 [ 0x1.7p+13L,     T.infinity                               ], // close overflow
1227                 [ 0x1p+80L,       T.infinity                               ], // far overflow
1228                 [ T.infinity,     T.infinity                               ],
1229                 [-0x1.18p+13L,    0x1.5e4bf54b4807034ea97fef0059a6p-12927L ], // near underflow
1230                 [-0x1.625p+13L,   0x1.a6bd68a39d11fec3a250cd97f524p-16358L ], // ditto
1231                 [-0x1.62dafp+13L, 0x0.cb629e9813b80ed4d639e875be6cp-16382L ], // near underflow - subnormal
1232                 [-0x1.6549p+13L,  0x0.0000000000000000000000000001p-16382L ], // ditto
1233                 [-0x1.655p+13L,   0                                        ], // close underflow
1234                 [-0x1p+30L,       0                                        ], // far underflow
1235             ];
1236         }
1237         else static if (realFormat == RealFormat.ieeeExtended ||
1238                         realFormat == RealFormat.ieeeExtended53)
1239         {
1240             static immutable T[2][] exptestpoints =
1241             [ //  x               exp(x)
1242                 [ 1.0L,           E                            ],
1243                 [ 0.5L,           0x1.a61298e1e069bc97p+0L     ],
1244                 [ 3.0L,           E*E*E                        ],
1245                 [ 0x1.1p+13L,     0x1.29aeffefc8ec645p+12557L  ], // near overflow
1246                 [ 0x1.7p+13L,     T.infinity                   ], // close overflow
1247                 [ 0x1p+80L,       T.infinity                   ], // far overflow
1248                 [ T.infinity,     T.infinity                   ],
1249                 [-0x1.18p+13L,    0x1.5e4bf54b4806db9p-12927L  ], // near underflow
1250                 [-0x1.625p+13L,   0x1.a6bd68a39d11f35cp-16358L ], // ditto
1251                 [-0x1.62dafp+13L, 0x1.96c53d30277021dp-16383L  ], // near underflow - subnormal
1252                 [-0x1.643p+13L,   0x1p-16444L                  ], // ditto
1253                 [-0x1.645p+13L,   0                            ], // close underflow
1254                 [-0x1p+30L,       0                            ], // far underflow
1255             ];
1256         }
1257         else static if (realFormat == RealFormat.ieeeDouble)
1258         {
1259             static immutable T[2][] exptestpoints =
1260             [ //  x,             exp(x)
1261                 [ 1.0L,          E                        ],
1262                 [ 0.5L,          0x1.a61298e1e069cp+0L    ],
1263                 [ 3.0L,          E*E*E                    ],
1264                 [ 0x1.6p+9L,     0x1.93bf4ec282efbp+1015L ], // near overflow
1265                 [ 0x1.7p+9L,     T.infinity               ], // close overflow
1266                 [ 0x1p+80L,      T.infinity               ], // far overflow
1267                 [ T.infinity,    T.infinity               ],
1268                 [-0x1.6p+9L,     0x1.44a3824e5285fp-1016L ], // near underflow
1269                 [-0x1.64p+9L,    0x0.06f84920bb2d4p-1022L ], // near underflow - subnormal
1270                 [-0x1.743p+9L,   0x0.0000000000001p-1022L ], // ditto
1271                 [-0x1.8p+9L,     0                        ], // close underflow
1272                 [-0x1p+30L,      0                        ], // far underflow
1273             ];
1274         }
1275         else static if (realFormat == RealFormat.ieeeSingle)
1276         {
1277             static immutable T[2][] exptestpoints =
1278             [ //  x,             exp(x)
1279                 [ 1.0L,          E                ],
1280                 [ 0.5L,          0x1.a61299p+0L   ],
1281                 [ 3.0L,          E*E*E            ],
1282                 [ 0x1.62p+6L,    0x1.99b988p+127L ], // near overflow
1283                 [ 0x1.7p+6L,     T.infinity       ], // close overflow
1284                 [ 0x1p+80L,      T.infinity       ], // far overflow
1285                 [ T.infinity,    T.infinity       ],
1286                 [-0x1.5cp+6L,    0x1.666d0ep-126L ], // near underflow
1287                 [-0x1.7p+6L,     0x0.026a42p-126L ], // near underflow - subnormal
1288                 [-0x1.9cp+6L,    0x0.000002p-126L ], // ditto
1289                 [-0x1.ap+6L,     0                ], // close underflow
1290                 [-0x1p+30L,      0                ], // far underflow
1291             ];
1292         }
1293         else
1294             static assert(0, "No exp() tests for real type!");
1295 
1296         const minEqualMantissaBits = T.mant_dig - 2;
1297         T x;
1298         version (IeeeFlagsSupport) IeeeFlags f;
1299         foreach (ref pair; exptestpoints)
1300         {
1301             version (IeeeFlagsSupport) resetIeeeFlags();
1302             x = exp(pair[0]);
1303             //printf("exp(%La) = %La, should be %La\n", cast(real) pair[0], cast(real) x, cast(real) pair[1]);
1304             assert(feqrel(x, pair[1]) >= minEqualMantissaBits);
1305         }
1306 
1307         // Ideally, exp(0) would not set the inexact flag.
1308         // Unfortunately, fldl2e sets it!
1309         // So it's not realistic to avoid setting it.
1310         assert(exp(cast(T) 0.0) == 1.0);
1311 
1312         // NaN propagation. Doesn't set flags, bcos was already NaN.
1313         version (IeeeFlagsSupport)
1314         {
1315             resetIeeeFlags();
1316             x = exp(T.nan);
1317             f = ieeeFlags;
1318             assert(isIdentical(abs(x), T.nan));
1319             assert(f.flags == 0);
1320 
1321             resetIeeeFlags();
1322             x = exp(-T.nan);
1323             f = ieeeFlags;
1324             assert(isIdentical(abs(x), T.nan));
1325             assert(f.flags == 0);
1326         }
1327         else
1328         {
1329             x = exp(T.nan);
1330             assert(isIdentical(abs(x), T.nan));
1331 
1332             x = exp(-T.nan);
1333             assert(isIdentical(abs(x), T.nan));
1334         }
1335 
1336         x = exp(NaN(0x123));
1337         assert(isIdentical(x, NaN(0x123)));
1338     }
1339 
1340     import std.meta : AliasSeq;
1341     foreach (T; AliasSeq!(real, double, float))
1342         testExp!T();
1343 
1344     // High resolution test (verified against GNU MPFR/Mathematica).
1345     assert(exp(0.5L) == 0x1.A612_98E1_E069_BC97_2DFE_FAB6_DF34p+0L);
1346 
1347     assert(isClose(exp(3.0L), E * E * E, real.sizeof > double.sizeof ? 1e-15 : 1e-14));
1348 }
1349 
1350 /**
1351  * Calculates the value of the natural logarithm base (e)
1352  * raised to the power of x, minus 1.
1353  *
1354  * For very small x, expm1(x) is more accurate
1355  * than exp(x)-1.
1356  *
1357  *  $(TABLE_SV
1358  *    $(TR $(TH x)             $(TH e$(SUPERSCRIPT x)-1)  )
1359  *    $(TR $(TD $(PLUSMN)0.0)  $(TD $(PLUSMN)0.0) )
1360  *    $(TR $(TD +$(INFIN))     $(TD +$(INFIN))    )
1361  *    $(TR $(TD -$(INFIN))     $(TD -1.0)         )
1362  *    $(TR $(TD $(NAN))        $(TD $(NAN))       )
1363  *  )
1364  */
1365 pragma(inline, true)
1366 real expm1(real x) @trusted pure nothrow @nogc // TODO: @safe
1367 {
1368     version (InlineAsm_X87)
1369     {
1370         if (!__ctfe)
1371             return expm1Asm(x);
1372     }
1373     return expm1Impl(x);
1374 }
1375 
1376 /// ditto
1377 pragma(inline, true)
1378 double expm1(double x) @safe pure nothrow @nogc
1379 {
1380     return __ctfe ? cast(double) expm1(cast(real) x) : expm1Impl(x);
1381 }
1382 
1383 /// ditto
1384 pragma(inline, true)
1385 float expm1(float x) @safe pure nothrow @nogc
1386 {
1387     // no single-precision version in Cephes => use double precision
1388     return __ctfe ? cast(float) expm1(cast(real) x) : cast(float) expm1Impl(cast(double) x);
1389 }
1390 
1391 ///
1392 @safe unittest
1393 {
1394     import std.math.traits : isIdentical;
1395     import std.math.operations : feqrel;
1396 
1397     assert(isIdentical(expm1(0.0), 0.0));
1398     assert(expm1(1.0).feqrel(1.71828) > 16);
1399     assert(expm1(2.0).feqrel(6.3890) > 16);
1400 }
1401 
1402 version (InlineAsm_X87)
1403 private real expm1Asm(real x) @trusted pure nothrow @nogc
1404 {
1405     version (X86)
1406     {
1407         enum PARAMSIZE = (real.sizeof+3)&(0xFFFF_FFFC); // always a multiple of 4
1408         asm pure nothrow @nogc
1409         {
1410             /*  expm1() for x87 80-bit reals, IEEE754-2008 conformant.
1411              * Author: Don Clugston.
1412              *
1413              *    expm1(x) = 2^^(rndint(y))* 2^^(y-rndint(y)) - 1 where y = LN2*x.
1414              *    = 2rndy * 2ym1 + 2rndy - 1, where 2rndy = 2^^(rndint(y))
1415              *     and 2ym1 = (2^^(y-rndint(y))-1).
1416              *    If 2rndy  < 0.5*real.epsilon, result is -1.
1417              *    Implementation is otherwise the same as for exp2()
1418              */
1419             naked;
1420             fld real ptr [ESP+4] ; // x
1421             mov AX, [ESP+4+8]; // AX = exponent and sign
1422             sub ESP, 12+8; // Create scratch space on the stack
1423             // [ESP,ESP+2] = scratchint
1424             // [ESP+4..+6, +8..+10, +10] = scratchreal
1425             // set scratchreal mantissa = 1.0
1426             mov dword ptr [ESP+8], 0;
1427             mov dword ptr [ESP+8+4], 0x80000000;
1428             and AX, 0x7FFF; // drop sign bit
1429             cmp AX, 0x401D; // avoid InvalidException in fist
1430             jae L_extreme;
1431             fldl2e;
1432             fmulp ST(1), ST; // y = x*log2(e)
1433             fist dword ptr [ESP]; // scratchint = rndint(y)
1434             fisub dword ptr [ESP]; // y - rndint(y)
1435             // and now set scratchreal exponent
1436             mov EAX, [ESP];
1437             add EAX, 0x3fff;
1438             jle short L_largenegative;
1439             cmp EAX,0x8000;
1440             jge short L_largepositive;
1441             mov [ESP+8+8],AX;
1442             f2xm1; // 2ym1 = 2^^(y-rndint(y)) -1
1443             fld real ptr [ESP+8] ; // 2rndy = 2^^rndint(y)
1444             fmul ST(1), ST;  // ST=2rndy, ST(1)=2rndy*2ym1
1445             fld1;
1446             fsubp ST(1), ST; // ST = 2rndy-1, ST(1) = 2rndy * 2ym1 - 1
1447             faddp ST(1), ST; // ST = 2rndy * 2ym1 + 2rndy - 1
1448             add ESP,12+8;
1449             ret PARAMSIZE;
1450 
1451 L_extreme:  // Extreme exponent. X is very large positive, very
1452             // large negative, infinity, or NaN.
1453             fxam;
1454             fstsw AX;
1455             test AX, 0x0400; // NaN_or_zero, but we already know x != 0
1456             jz L_was_nan;  // if x is NaN, returns x
1457             test AX, 0x0200;
1458             jnz L_largenegative;
1459 L_largepositive:
1460             // Set scratchreal = real.max.
1461             // squaring it will create infinity, and set overflow flag.
1462             mov word  ptr [ESP+8+8], 0x7FFE;
1463             fstp ST(0);
1464             fld real ptr [ESP+8];  // load scratchreal
1465             fmul ST(0), ST;        // square it, to create havoc!
1466 L_was_nan:
1467             add ESP,12+8;
1468             ret PARAMSIZE;
1469 L_largenegative:
1470             fstp ST(0);
1471             fld1;
1472             fchs; // return -1. Underflow flag is not set.
1473             add ESP,12+8;
1474             ret PARAMSIZE;
1475         }
1476     }
1477     else version (X86_64)
1478     {
1479         asm pure nothrow @nogc
1480         {
1481             naked;
1482         }
1483         version (Win64)
1484         {
1485             asm pure nothrow @nogc
1486             {
1487                 fld   real ptr [RCX];  // x
1488                 mov   AX,[RCX+8];      // AX = exponent and sign
1489             }
1490         }
1491         else
1492         {
1493             asm pure nothrow @nogc
1494             {
1495                 fld   real ptr [RSP+8];  // x
1496                 mov   AX,[RSP+8+8];      // AX = exponent and sign
1497             }
1498         }
1499         asm pure nothrow @nogc
1500         {
1501             /*  expm1() for x87 80-bit reals, IEEE754-2008 conformant.
1502              * Author: Don Clugston.
1503              *
1504              *    expm1(x) = 2^(rndint(y))* 2^(y-rndint(y)) - 1 where y = LN2*x.
1505              *    = 2rndy * 2ym1 + 2rndy - 1, where 2rndy = 2^(rndint(y))
1506              *     and 2ym1 = (2^(y-rndint(y))-1).
1507              *    If 2rndy  < 0.5*real.epsilon, result is -1.
1508              *    Implementation is otherwise the same as for exp2()
1509              */
1510             sub RSP, 24;       // Create scratch space on the stack
1511             // [RSP,RSP+2] = scratchint
1512             // [RSP+4..+6, +8..+10, +10] = scratchreal
1513             // set scratchreal mantissa = 1.0
1514             mov dword ptr [RSP+8], 0;
1515             mov dword ptr [RSP+8+4], 0x80000000;
1516             and AX, 0x7FFF; // drop sign bit
1517             cmp AX, 0x401D; // avoid InvalidException in fist
1518             jae L_extreme;
1519             fldl2e;
1520             fmul ; // y = x*log2(e)
1521             fist dword ptr [RSP]; // scratchint = rndint(y)
1522             fisub dword ptr [RSP]; // y - rndint(y)
1523             // and now set scratchreal exponent
1524             mov EAX, [RSP];
1525             add EAX, 0x3fff;
1526             jle short L_largenegative;
1527             cmp EAX,0x8000;
1528             jge short L_largepositive;
1529             mov [RSP+8+8],AX;
1530             f2xm1; // 2^(y-rndint(y)) -1
1531             fld real ptr [RSP+8] ; // 2^rndint(y)
1532             fmul ST(1), ST;
1533             fld1;
1534             fsubp ST(1), ST;
1535             fadd;
1536             add RSP,24;
1537             ret;
1538 
1539 L_extreme:  // Extreme exponent. X is very large positive, very
1540             // large negative, infinity, or NaN.
1541             fxam;
1542             fstsw AX;
1543             test AX, 0x0400; // NaN_or_zero, but we already know x != 0
1544             jz L_was_nan;  // if x is NaN, returns x
1545             test AX, 0x0200;
1546             jnz L_largenegative;
1547 L_largepositive:
1548             // Set scratchreal = real.max.
1549             // squaring it will create infinity, and set overflow flag.
1550             mov word  ptr [RSP+8+8], 0x7FFE;
1551             fstp ST(0);
1552             fld real ptr [RSP+8];  // load scratchreal
1553             fmul ST(0), ST;        // square it, to create havoc!
1554 L_was_nan:
1555             add RSP,24;
1556             ret;
1557 
1558 L_largenegative:
1559             fstp ST(0);
1560             fld1;
1561             fchs; // return -1. Underflow flag is not set.
1562             add RSP,24;
1563             ret;
1564         }
1565     }
1566     else
1567         static assert(0);
1568 }
1569 
1570 private T expm1Impl(T)(T x) @safe pure nothrow @nogc
1571 {
1572     import std.math : floatTraits, RealFormat;
1573     import std.math.rounding : floor;
1574     import std.math.algebraic : poly;
1575     import std.math.constants : LN2;
1576 
1577     // Coefficients for exp(x) - 1 and overflow/underflow limits.
1578     enum realFormat = floatTraits!T.realFormat;
1579     static if (realFormat == RealFormat.ieeeQuadruple)
1580     {
1581         static immutable T[8] P = [
1582             2.943520915569954073888921213330863757240E8L,
1583             -5.722847283900608941516165725053359168840E7L,
1584             8.944630806357575461578107295909719817253E6L,
1585             -7.212432713558031519943281748462837065308E5L,
1586             4.578962475841642634225390068461943438441E4L,
1587             -1.716772506388927649032068540558788106762E3L,
1588             4.401308817383362136048032038528753151144E1L,
1589             -4.888737542888633647784737721812546636240E-1L
1590         ];
1591 
1592         static immutable T[9] Q = [
1593             1.766112549341972444333352727998584753865E9L,
1594             -7.848989743695296475743081255027098295771E8L,
1595             1.615869009634292424463780387327037251069E8L,
1596             -2.019684072836541751428967854947019415698E7L,
1597             1.682912729190313538934190635536631941751E6L,
1598             -9.615511549171441430850103489315371768998E4L,
1599             3.697714952261803935521187272204485251835E3L,
1600             -8.802340681794263968892934703309274564037E1L,
1601             1.0
1602         ];
1603 
1604         enum T OF = 1.1356523406294143949491931077970764891253E4L;
1605         enum T UF = -1.143276959615573793352782661133116431383730e4L;
1606     }
1607     else static if (realFormat == RealFormat.ieeeExtended)
1608     {
1609         static immutable T[5] P = [
1610            -1.586135578666346600772998894928250240826E4L,
1611             2.642771505685952966904660652518429479531E3L,
1612            -3.423199068835684263987132888286791620673E2L,
1613             1.800826371455042224581246202420972737840E1L,
1614            -5.238523121205561042771939008061958820811E-1L,
1615         ];
1616         static immutable T[6] Q = [
1617            -9.516813471998079611319047060563358064497E4L,
1618             3.964866271411091674556850458227710004570E4L,
1619            -7.207678383830091850230366618190187434796E3L,
1620             7.206038318724600171970199625081491823079E2L,
1621            -4.002027679107076077238836622982900945173E1L,
1622             1.0
1623         ];
1624 
1625         enum T OF =  1.1356523406294143949492E4L;
1626         enum T UF = -4.5054566736396445112120088E1L;
1627     }
1628     else static if (realFormat == RealFormat.ieeeDouble)
1629     {
1630         static immutable T[3] P = [
1631             9.9999999999999999991025E-1,
1632             3.0299440770744196129956E-2,
1633             1.2617719307481059087798E-4,
1634         ];
1635         static immutable T[4] Q = [
1636             2.0000000000000000000897E0,
1637             2.2726554820815502876593E-1,
1638             2.5244834034968410419224E-3,
1639             3.0019850513866445504159E-6,
1640         ];
1641     }
1642     else
1643         static assert(0, "no coefficients for expm1()");
1644 
1645     static if (realFormat == RealFormat.ieeeDouble) // special case for double precision
1646     {
1647         if (x < -0.5 || x > 0.5)
1648             return exp(x) - 1.0;
1649         if (x == 0.0)
1650             return x;
1651 
1652         const T xx = x * x;
1653         x = x * poly(xx, P);
1654         x = x / (poly(xx, Q) - x);
1655         return x + x;
1656     }
1657     else
1658     {
1659         // C1 + C2 = LN2.
1660         enum T C1 = 6.9314575195312500000000E-1L;
1661         enum T C2 = 1.428606820309417232121458176568075500134E-6L;
1662 
1663         // Special cases.
1664         if (x > OF)
1665             return real.infinity;
1666         if (x == cast(T) 0.0)
1667             return x;
1668         if (x < UF)
1669             return -1.0;
1670 
1671         // Express x = LN2 (n + remainder), remainder not exceeding 1/2.
1672         int n = cast(int) floor((cast(T) 0.5) + x / cast(T) LN2);
1673         x -= n * C1;
1674         x -= n * C2;
1675 
1676         // Rational approximation:
1677         //  exp(x) - 1 = x + 0.5 x^^2 + x^^3 P(x) / Q(x)
1678         T px = x * poly(x, P);
1679         T qx = poly(x, Q);
1680         const T xx = x * x;
1681         qx = x + ((cast(T) 0.5) * xx + xx * px / qx);
1682 
1683         // We have qx = exp(remainder LN2) - 1, so:
1684         //  exp(x) - 1 = 2^^n (qx + 1) - 1 = 2^^n qx + 2^^n - 1.
1685         px = core.math.ldexp(cast(T) 1.0, n);
1686         x = px * qx + (px - cast(T) 1.0);
1687 
1688         return x;
1689     }
1690 }
1691 
1692 @safe @nogc nothrow unittest
1693 {
1694     import std.math.traits : isNaN;
1695     import std.math.operations : isClose, CommonDefaultFor;
1696 
1697     static void testExpm1(T)()
1698     {
1699         // NaN
1700         assert(isNaN(expm1(cast(T) T.nan)));
1701 
1702         static immutable T[] xs = [ -2, -0.75, -0.3, 0.0, 0.1, 0.2, 0.5, 1.0 ];
1703         foreach (x; xs)
1704         {
1705             const T e = expm1(x);
1706             const T r = exp(x) - 1;
1707 
1708             //printf("expm1(%Lg) = %Lg, should approximately be %Lg\n", cast(real) x, cast(real) e, cast(real) r);
1709             assert(isClose(r, e, CommonDefaultFor!(T,T), CommonDefaultFor!(T,T)));
1710         }
1711     }
1712 
1713     import std.meta : AliasSeq;
1714     foreach (T; AliasSeq!(real, double))
1715         testExpm1!T();
1716 }
1717 
1718 /**
1719  * Calculates 2$(SUPERSCRIPT x).
1720  *
1721  *  $(TABLE_SV
1722  *    $(TR $(TH x)             $(TH exp2(x))   )
1723  *    $(TR $(TD +$(INFIN))     $(TD +$(INFIN)) )
1724  *    $(TR $(TD -$(INFIN))     $(TD +0.0)      )
1725  *    $(TR $(TD $(NAN))        $(TD $(NAN))    )
1726  *  )
1727  */
1728 version (none) // LDC FIXME: Use of this LLVM intrinsic causes a unit test failure
1729 {
1730     pragma(inline, true):
1731     real   exp2(real   x) @safe pure nothrow @nogc { return llvm_exp2(x); }
1732     ///ditto
1733     double exp2(double x) @safe pure nothrow @nogc { return llvm_exp2(x); }
1734     ///ditto
1735     float  exp2(float  x) @safe pure nothrow @nogc { return llvm_exp2(x); }
1736 }
1737 else
1738 {
1739 
1740 pragma(inline, true)
1741 real exp2(real x) @nogc @trusted pure nothrow // TODO: @safe
1742 {
1743     version (InlineAsm_X87)
1744     {
1745         if (!__ctfe)
1746             return exp2Asm(x);
1747     }
1748     return exp2Impl(x);
1749 }
1750 
1751 /// ditto
1752 pragma(inline, true)
1753 double exp2(double x) @nogc @safe pure nothrow { return __ctfe ? cast(double) exp2(cast(real) x) : exp2Impl(x); }
1754 
1755 /// ditto
1756 pragma(inline, true)
1757 float exp2(float x) @nogc @safe pure nothrow { return __ctfe ? cast(float) exp2(cast(real) x) : exp2Impl(x); }
1758 
1759 } // !none
1760 
1761 ///
1762 @safe unittest
1763 {
1764     import std.math.traits : isIdentical;
1765     import std.math.operations : feqrel;
1766 
1767     assert(isIdentical(exp2(0.0), 1.0));
1768     assert(exp2(2.0).feqrel(4.0) > 16);
1769     assert(exp2(8.0).feqrel(256.0) > 16);
1770 }
1771 
1772 @safe unittest
1773 {
1774     version (CRuntime_Microsoft) {} else // aexp2/exp2f/exp2l not implemented
1775     {
1776         assert( core.stdc.math.exp2f(0.0f) == 1 );
1777         assert( core.stdc.math.exp2 (0.0)  == 1 );
1778         assert( core.stdc.math.exp2l(0.0L) == 1 );
1779     }
1780 }
1781 
1782 version (InlineAsm_X87)
1783 private real exp2Asm(real x) @nogc @trusted pure nothrow
1784 {
1785     version (X86)
1786     {
1787         enum PARAMSIZE = (real.sizeof+3)&(0xFFFF_FFFC); // always a multiple of 4
1788 
1789         asm pure nothrow @nogc
1790         {
1791             /*  exp2() for x87 80-bit reals, IEEE754-2008 conformant.
1792              * Author: Don Clugston.
1793              *
1794              * exp2(x) = 2^^(rndint(x))* 2^^(y-rndint(x))
1795              * The trick for high performance is to avoid the fscale(28cycles on core2),
1796              * frndint(19 cycles), leaving f2xm1(19 cycles) as the only slow instruction.
1797              *
1798              * We can do frndint by using fist. BUT we can't use it for huge numbers,
1799              * because it will set the Invalid Operation flag if overflow or NaN occurs.
1800              * Fortunately, whenever this happens the result would be zero or infinity.
1801              *
1802              * We can perform fscale by directly poking into the exponent. BUT this doesn't
1803              * work for the (very rare) cases where the result is subnormal. So we fall back
1804              * to the slow method in that case.
1805              */
1806             naked;
1807             fld real ptr [ESP+4] ; // x
1808             mov AX, [ESP+4+8]; // AX = exponent and sign
1809             sub ESP, 12+8; // Create scratch space on the stack
1810             // [ESP,ESP+2] = scratchint
1811             // [ESP+4..+6, +8..+10, +10] = scratchreal
1812             // set scratchreal mantissa = 1.0
1813             mov dword ptr [ESP+8], 0;
1814             mov dword ptr [ESP+8+4], 0x80000000;
1815             and AX, 0x7FFF; // drop sign bit
1816             cmp AX, 0x401D; // avoid InvalidException in fist
1817             jae L_extreme;
1818             fist dword ptr [ESP]; // scratchint = rndint(x)
1819             fisub dword ptr [ESP]; // x - rndint(x)
1820             // and now set scratchreal exponent
1821             mov EAX, [ESP];
1822             add EAX, 0x3fff;
1823             jle short L_subnormal;
1824             cmp EAX,0x8000;
1825             jge short L_overflow;
1826             mov [ESP+8+8],AX;
1827 L_normal:
1828             f2xm1;
1829             fld1;
1830             faddp ST(1), ST; // 2^^(x-rndint(x))
1831             fld real ptr [ESP+8] ; // 2^^rndint(x)
1832             add ESP,12+8;
1833             fmulp ST(1), ST;
1834             ret PARAMSIZE;
1835 
1836 L_subnormal:
1837             // Result will be subnormal.
1838             // In this rare case, the simple poking method doesn't work.
1839             // The speed doesn't matter, so use the slow fscale method.
1840             fild dword ptr [ESP];  // scratchint
1841             fld1;
1842             fscale;
1843             fstp real ptr [ESP+8]; // scratchreal = 2^^scratchint
1844             fstp ST(0);         // drop scratchint
1845             jmp L_normal;
1846 
1847 L_extreme:  // Extreme exponent. X is very large positive, very
1848             // large negative, infinity, or NaN.
1849             fxam;
1850             fstsw AX;
1851             test AX, 0x0400; // NaN_or_zero, but we already know x != 0
1852             jz L_was_nan;  // if x is NaN, returns x
1853             // set scratchreal = real.min_normal
1854             // squaring it will return 0, setting underflow flag
1855             mov word  ptr [ESP+8+8], 1;
1856             test AX, 0x0200;
1857             jnz L_waslargenegative;
1858 L_overflow:
1859             // Set scratchreal = real.max.
1860             // squaring it will create infinity, and set overflow flag.
1861             mov word  ptr [ESP+8+8], 0x7FFE;
1862 L_waslargenegative:
1863             fstp ST(0);
1864             fld real ptr [ESP+8];  // load scratchreal
1865             fmul ST(0), ST;        // square it, to create havoc!
1866 L_was_nan:
1867             add ESP,12+8;
1868             ret PARAMSIZE;
1869         }
1870     }
1871     else version (X86_64)
1872     {
1873         asm pure nothrow @nogc
1874         {
1875             naked;
1876         }
1877         version (Win64)
1878         {
1879             asm pure nothrow @nogc
1880             {
1881                 fld   real ptr [RCX];  // x
1882                 mov   AX,[RCX+8];      // AX = exponent and sign
1883             }
1884         }
1885         else
1886         {
1887             asm pure nothrow @nogc
1888             {
1889                 fld   real ptr [RSP+8];  // x
1890                 mov   AX,[RSP+8+8];      // AX = exponent and sign
1891             }
1892         }
1893         asm pure nothrow @nogc
1894         {
1895             /*  exp2() for x87 80-bit reals, IEEE754-2008 conformant.
1896              * Author: Don Clugston.
1897              *
1898              * exp2(x) = 2^(rndint(x))* 2^(y-rndint(x))
1899              * The trick for high performance is to avoid the fscale(28cycles on core2),
1900              * frndint(19 cycles), leaving f2xm1(19 cycles) as the only slow instruction.
1901              *
1902              * We can do frndint by using fist. BUT we can't use it for huge numbers,
1903              * because it will set the Invalid Operation flag is overflow or NaN occurs.
1904              * Fortunately, whenever this happens the result would be zero or infinity.
1905              *
1906              * We can perform fscale by directly poking into the exponent. BUT this doesn't
1907              * work for the (very rare) cases where the result is subnormal. So we fall back
1908              * to the slow method in that case.
1909              */
1910             sub RSP, 24; // Create scratch space on the stack
1911             // [RSP,RSP+2] = scratchint
1912             // [RSP+4..+6, +8..+10, +10] = scratchreal
1913             // set scratchreal mantissa = 1.0
1914             mov dword ptr [RSP+8], 0;
1915             mov dword ptr [RSP+8+4], 0x80000000;
1916             and AX, 0x7FFF; // drop sign bit
1917             cmp AX, 0x401D; // avoid InvalidException in fist
1918             jae L_extreme;
1919             fist dword ptr [RSP]; // scratchint = rndint(x)
1920             fisub dword ptr [RSP]; // x - rndint(x)
1921             // and now set scratchreal exponent
1922             mov EAX, [RSP];
1923             add EAX, 0x3fff;
1924             jle short L_subnormal;
1925             cmp EAX,0x8000;
1926             jge short L_overflow;
1927             mov [RSP+8+8],AX;
1928 L_normal:
1929             f2xm1;
1930             fld1;
1931             fadd; // 2^(x-rndint(x))
1932             fld real ptr [RSP+8] ; // 2^rndint(x)
1933             add RSP,24;
1934             fmulp ST(1), ST;
1935             ret;
1936 
1937 L_subnormal:
1938             // Result will be subnormal.
1939             // In this rare case, the simple poking method doesn't work.
1940             // The speed doesn't matter, so use the slow fscale method.
1941             fild dword ptr [RSP];  // scratchint
1942             fld1;
1943             fscale;
1944             fstp real ptr [RSP+8]; // scratchreal = 2^scratchint
1945             fstp ST(0);         // drop scratchint
1946             jmp L_normal;
1947 
1948 L_extreme:  // Extreme exponent. X is very large positive, very
1949             // large negative, infinity, or NaN.
1950             fxam;
1951             fstsw AX;
1952             test AX, 0x0400; // NaN_or_zero, but we already know x != 0
1953             jz L_was_nan;  // if x is NaN, returns x
1954             // set scratchreal = real.min
1955             // squaring it will return 0, setting underflow flag
1956             mov word  ptr [RSP+8+8], 1;
1957             test AX, 0x0200;
1958             jnz L_waslargenegative;
1959 L_overflow:
1960             // Set scratchreal = real.max.
1961             // squaring it will create infinity, and set overflow flag.
1962             mov word  ptr [RSP+8+8], 0x7FFE;
1963 L_waslargenegative:
1964             fstp ST(0);
1965             fld real ptr [RSP+8];  // load scratchreal
1966             fmul ST(0), ST;        // square it, to create havoc!
1967 L_was_nan:
1968             add RSP,24;
1969             ret;
1970         }
1971     }
1972     else
1973         static assert(0);
1974 }
1975 
1976 private T exp2Impl(T)(T x) @nogc @safe pure nothrow
1977 {
1978     import std.math : floatTraits, RealFormat;
1979     import std.math.traits : isNaN;
1980     import std.math.rounding : floor;
1981     import std.math.algebraic : poly;
1982 
1983     // Coefficients for exp2(x)
1984     enum realFormat = floatTraits!T.realFormat;
1985     static if (realFormat == RealFormat.ieeeQuadruple)
1986     {
1987         static immutable T[5] P = [
1988             9.079594442980146270952372234833529694788E12L,
1989             1.530625323728429161131811299626419117557E11L,
1990             5.677513871931844661829755443994214173883E8L,
1991             6.185032670011643762127954396427045467506E5L,
1992             1.587171580015525194694938306936721666031E2L
1993         ];
1994 
1995         static immutable T[6] Q = [
1996             2.619817175234089411411070339065679229869E13L,
1997             1.490560994263653042761789432690793026977E12L,
1998             1.092141473886177435056423606755843616331E10L,
1999             2.186249607051644894762167991800811827835E7L,
2000             1.236602014442099053716561665053645270207E4L,
2001             1.0
2002         ];
2003     }
2004     else static if (realFormat == RealFormat.ieeeExtended)
2005     {
2006         static immutable T[3] P = [
2007             2.0803843631901852422887E6L,
2008             3.0286971917562792508623E4L,
2009             6.0614853552242266094567E1L,
2010         ];
2011         static immutable T[4] Q = [
2012             6.0027204078348487957118E6L,
2013             3.2772515434906797273099E5L,
2014             1.7492876999891839021063E3L,
2015             1.0000000000000000000000E0L,
2016         ];
2017     }
2018     else static if (realFormat == RealFormat.ieeeDouble)
2019     {
2020         static immutable T[3] P = [
2021             1.51390680115615096133E3L,
2022             2.02020656693165307700E1L,
2023             2.30933477057345225087E-2L,
2024         ];
2025         static immutable T[3] Q = [
2026             4.36821166879210612817E3L,
2027             2.33184211722314911771E2L,
2028             1.00000000000000000000E0L,
2029         ];
2030     }
2031     else static if (realFormat == RealFormat.ieeeSingle)
2032     {
2033         static immutable T[6] P = [
2034             6.931472028550421E-001L,
2035             2.402264791363012E-001L,
2036             5.550332471162809E-002L,
2037             9.618437357674640E-003L,
2038             1.339887440266574E-003L,
2039             1.535336188319500E-004L,
2040         ];
2041     }
2042     else
2043         static assert(0, "no coefficients for exp2()");
2044 
2045     // Overflow and Underflow limits.
2046     enum T OF = T.max_exp;
2047     enum T UF = T.min_exp - 1;
2048 
2049     // Special cases.
2050     if (isNaN(x))
2051         return x;
2052     if (x > OF)
2053         return real.infinity;
2054     if (x < UF)
2055         return 0.0;
2056 
2057     static if (realFormat == RealFormat.ieeeSingle) // special case for single precision
2058     {
2059         // The following is necessary because range reduction blows up.
2060         if (x == 0.0f)
2061             return 1.0f;
2062 
2063         // Separate into integer and fractional parts.
2064         const T i = floor(x);
2065         int n = cast(int) i;
2066         x -= i;
2067         if (x > 0.5f)
2068         {
2069             n += 1;
2070             x -= 1.0f;
2071         }
2072 
2073         // Rational approximation:
2074         //  exp2(x) = 1.0 + x P(x)
2075         x = 1.0f + x * poly(x, P);
2076     }
2077     else
2078     {
2079         // Separate into integer and fractional parts.
2080         const T i = floor(x + cast(T) 0.5);
2081         int n = cast(int) i;
2082         x -= i;
2083 
2084         // Rational approximation:
2085         //  exp2(x) = 1.0 + 2x P(x^^2) / (Q(x^^2) - P(x^^2))
2086         const T xx = x * x;
2087         const T px = x * poly(xx, P);
2088         x = px / (poly(xx, Q) - px);
2089         x = (cast(T) 1.0) + (cast(T) 2.0) * x;
2090     }
2091 
2092     // Scale by power of 2.
2093     x = core.math.ldexp(x, n);
2094 
2095     return x;
2096 }
2097 
2098 @safe @nogc nothrow unittest
2099 {
2100     import std.math.operations : feqrel, NaN, isClose;
2101     import std.math.traits : isIdentical;
2102     import std.math.constants : SQRT2;
2103 
2104     assert(feqrel(exp2(0.5L), SQRT2) >= real.mant_dig -1);
2105     assert(exp2(8.0L) == 256.0);
2106     assert(exp2(-9.0L)== 1.0L/512.0);
2107 
2108     static void testExp2(T)()
2109     {
2110         // NaN
2111         const T specialNaN = NaN(0x0123L);
2112         assert(isIdentical(exp2(specialNaN), specialNaN));
2113 
2114         // over-/underflow
2115         enum T OF = T.max_exp;
2116         enum T UF = T.min_exp - T.mant_dig;
2117         assert(isIdentical(exp2(OF + 1), cast(T) T.infinity));
2118         assert(isIdentical(exp2(UF - 1), cast(T) 0.0));
2119 
2120         static immutable T[2][] vals =
2121         [
2122             // x, exp2(x)
2123             [  0.0, 1.0 ],
2124             [ -0.0, 1.0 ],
2125             [  0.5, SQRT2 ],
2126             [  8.0, 256.0 ],
2127             [ -9.0, 1.0 / 512 ],
2128         ];
2129 
2130         foreach (ref val; vals)
2131         {
2132             const T x = val[0];
2133             const T r = val[1];
2134             const T e = exp2(x);
2135 
2136             //printf("exp2(%Lg) = %Lg, should be %Lg\n", cast(real) x, cast(real) e, cast(real) r);
2137             assert(isClose(r, e));
2138         }
2139     }
2140 
2141     import std.meta : AliasSeq;
2142     foreach (T; AliasSeq!(real, double, float))
2143         testExp2!T();
2144 }
2145 
2146 /*********************************************************************
2147  * Separate floating point value into significand and exponent.
2148  *
2149  * Returns:
2150  *      Calculate and return $(I x) and $(I exp) such that
2151  *      value =$(I x)*2$(SUPERSCRIPT exp) and
2152  *      .5 $(LT)= |$(I x)| $(LT) 1.0
2153  *
2154  *      $(I x) has same sign as value.
2155  *
2156  *      $(TABLE_SV
2157  *      $(TR $(TH value)           $(TH returns)         $(TH exp))
2158  *      $(TR $(TD $(PLUSMN)0.0)    $(TD $(PLUSMN)0.0)    $(TD 0))
2159  *      $(TR $(TD +$(INFIN))       $(TD +$(INFIN))       $(TD int.max))
2160  *      $(TR $(TD -$(INFIN))       $(TD -$(INFIN))       $(TD int.min))
2161  *      $(TR $(TD $(PLUSMN)$(NAN)) $(TD $(PLUSMN)$(NAN)) $(TD int.min))
2162  *      )
2163  */
2164 T frexp(T)(const T value, out int exp) @trusted pure nothrow @nogc
2165 if (isFloatingPoint!T)
2166 {
2167     import std.math : floatTraits, RealFormat, MANTISSA_MSB, MANTISSA_LSB;
2168     import std.math.traits : isSubnormal;
2169 
2170     if (__ctfe)
2171     {
2172         // Handle special cases.
2173         if (value == 0) { exp = 0; return value; }
2174         if (value == T.infinity) { exp = int.max; return value; }
2175         if (value == -T.infinity || value != value) { exp = int.min; return value; }
2176         // Handle ordinary cases.
2177         // In CTFE there is no performance advantage for having separate
2178         // paths for different floating point types.
2179         T absValue = value < 0 ? -value : value;
2180         int expCount;
2181         static if (T.mant_dig > double.mant_dig)
2182         {
2183             for (; absValue >= 0x1.0p+1024L; absValue *= 0x1.0p-1024L)
2184                 expCount += 1024;
2185             for (; absValue < 0x1.0p-1021L; absValue *= 0x1.0p+1021L)
2186                 expCount -= 1021;
2187         }
2188         const double dval = cast(double) absValue;
2189         int dexp = cast(int) (((*cast(const long*) &dval) >>> 52) & 0x7FF) + double.min_exp - 2;
2190         dexp++;
2191         expCount += dexp;
2192         absValue *= 2.0 ^^ -dexp;
2193         // If the original value was subnormal or if it was a real
2194         // then absValue can still be outside the [0.5, 1.0) range.
2195         if (absValue < 0.5)
2196         {
2197             assert(T.mant_dig > double.mant_dig || isSubnormal(value));
2198             do
2199             {
2200                 absValue += absValue;
2201                 expCount--;
2202             } while (absValue < 0.5);
2203         }
2204         else
2205         {
2206             assert(absValue < 1 || T.mant_dig > double.mant_dig);
2207             for (; absValue >= 1; absValue *= T(0.5))
2208                 expCount++;
2209         }
2210         exp = expCount;
2211         return value < 0 ? -absValue : absValue;
2212     }
2213 
2214     Unqual!T vf = value;
2215     ushort* vu = cast(ushort*)&vf;
2216     static if (is(immutable T == immutable float))
2217         int* vi = cast(int*)&vf;
2218     else
2219         long* vl = cast(long*)&vf;
2220     int ex;
2221     alias F = floatTraits!T;
2222 
2223     ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
2224     static if (F.realFormat == RealFormat.ieeeExtended ||
2225                F.realFormat == RealFormat.ieeeExtended53)
2226     {
2227         if (ex)
2228         {   // If exponent is non-zero
2229             if (ex == F.EXPMASK) // infinity or NaN
2230             {
2231                 if (*vl &  0x7FFF_FFFF_FFFF_FFFF)  // NaN
2232                 {
2233                     *vl |= 0xC000_0000_0000_0000;  // convert NaNS to NaNQ
2234                     exp = int.min;
2235                 }
2236                 else if (vu[F.EXPPOS_SHORT] & 0x8000)   // negative infinity
2237                     exp = int.min;
2238                 else   // positive infinity
2239                     exp = int.max;
2240 
2241             }
2242             else
2243             {
2244                 exp = ex - F.EXPBIAS;
2245                 vu[F.EXPPOS_SHORT] = (0x8000 & vu[F.EXPPOS_SHORT]) | 0x3FFE;
2246             }
2247         }
2248         else if (!*vl)
2249         {
2250             // vf is +-0.0
2251             exp = 0;
2252         }
2253         else
2254         {
2255             // subnormal
2256 
2257             vf *= F.RECIP_EPSILON;
2258             ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
2259             exp = ex - F.EXPBIAS - T.mant_dig + 1;
2260             vu[F.EXPPOS_SHORT] = ((-1 - F.EXPMASK) & vu[F.EXPPOS_SHORT]) | 0x3FFE;
2261         }
2262         return vf;
2263     }
2264     else static if (F.realFormat == RealFormat.ieeeQuadruple)
2265     {
2266         if (ex)     // If exponent is non-zero
2267         {
2268             if (ex == F.EXPMASK)
2269             {
2270                 // infinity or NaN
2271                 if (vl[MANTISSA_LSB] |
2272                     (vl[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF))  // NaN
2273                 {
2274                     // convert NaNS to NaNQ
2275                     vl[MANTISSA_MSB] |= 0x0000_8000_0000_0000;
2276                     exp = int.min;
2277                 }
2278                 else if (vu[F.EXPPOS_SHORT] & 0x8000)   // negative infinity
2279                     exp = int.min;
2280                 else   // positive infinity
2281                     exp = int.max;
2282             }
2283             else
2284             {
2285                 exp = ex - F.EXPBIAS;
2286                 vu[F.EXPPOS_SHORT] = F.EXPBIAS | (0x8000 & vu[F.EXPPOS_SHORT]);
2287             }
2288         }
2289         else if ((vl[MANTISSA_LSB] |
2290             (vl[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) == 0)
2291         {
2292             // vf is +-0.0
2293             exp = 0;
2294         }
2295         else
2296         {
2297             // subnormal
2298             vf *= F.RECIP_EPSILON;
2299             ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
2300             exp = ex - F.EXPBIAS - T.mant_dig + 1;
2301             vu[F.EXPPOS_SHORT] = F.EXPBIAS | (0x8000 & vu[F.EXPPOS_SHORT]);
2302         }
2303         return vf;
2304     }
2305     else static if (F.realFormat == RealFormat.ieeeDouble)
2306     {
2307         if (ex) // If exponent is non-zero
2308         {
2309             if (ex == F.EXPMASK)   // infinity or NaN
2310             {
2311                 if (*vl == 0x7FF0_0000_0000_0000)  // positive infinity
2312                 {
2313                     exp = int.max;
2314                 }
2315                 else if (*vl == 0xFFF0_0000_0000_0000) // negative infinity
2316                     exp = int.min;
2317                 else
2318                 { // NaN
2319                     *vl |= 0x0008_0000_0000_0000;  // convert NaNS to NaNQ
2320                     exp = int.min;
2321                 }
2322             }
2323             else
2324             {
2325                 exp = (ex - F.EXPBIAS) >> 4;
2326                 vu[F.EXPPOS_SHORT] = cast(ushort)((0x800F & vu[F.EXPPOS_SHORT]) | 0x3FE0);
2327             }
2328         }
2329         else if (!(*vl & 0x7FFF_FFFF_FFFF_FFFF))
2330         {
2331             // vf is +-0.0
2332             exp = 0;
2333         }
2334         else
2335         {
2336             // subnormal
2337             vf *= F.RECIP_EPSILON;
2338             ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
2339             exp = ((ex - F.EXPBIAS) >> 4) - T.mant_dig + 1;
2340             vu[F.EXPPOS_SHORT] =
2341                 cast(ushort)(((-1 - F.EXPMASK) & vu[F.EXPPOS_SHORT]) | 0x3FE0);
2342         }
2343         return vf;
2344     }
2345     else static if (F.realFormat == RealFormat.ieeeSingle)
2346     {
2347         if (ex) // If exponent is non-zero
2348         {
2349             if (ex == F.EXPMASK)   // infinity or NaN
2350             {
2351                 if (*vi == 0x7F80_0000)  // positive infinity
2352                 {
2353                     exp = int.max;
2354                 }
2355                 else if (*vi == 0xFF80_0000) // negative infinity
2356                     exp = int.min;
2357                 else
2358                 { // NaN
2359                     *vi |= 0x0040_0000;  // convert NaNS to NaNQ
2360                     exp = int.min;
2361                 }
2362             }
2363             else
2364             {
2365                 exp = (ex - F.EXPBIAS) >> 7;
2366                 vu[F.EXPPOS_SHORT] = cast(ushort)((0x807F & vu[F.EXPPOS_SHORT]) | 0x3F00);
2367             }
2368         }
2369         else if (!(*vi & 0x7FFF_FFFF))
2370         {
2371             // vf is +-0.0
2372             exp = 0;
2373         }
2374         else
2375         {
2376             // subnormal
2377             vf *= F.RECIP_EPSILON;
2378             ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
2379             exp = ((ex - F.EXPBIAS) >> 7) - T.mant_dig + 1;
2380             vu[F.EXPPOS_SHORT] =
2381                 cast(ushort)(((-1 - F.EXPMASK) & vu[F.EXPPOS_SHORT]) | 0x3F00);
2382         }
2383         return vf;
2384     }
2385     else // static if (F.realFormat == RealFormat.ibmExtended)
2386     {
2387         assert(0, "frexp not implemented");
2388     }
2389 }
2390 
2391 ///
2392 @safe unittest
2393 {
2394     import std.math.operations : isClose;
2395 
2396     int exp;
2397     real mantissa = frexp(123.456L, exp);
2398 
2399     assert(isClose(mantissa * pow(2.0L, cast(real) exp), 123.456L));
2400 
2401     assert(frexp(-real.nan, exp) && exp == int.min);
2402     assert(frexp(real.nan, exp) && exp == int.min);
2403     assert(frexp(-real.infinity, exp) == -real.infinity && exp == int.min);
2404     assert(frexp(real.infinity, exp) == real.infinity && exp == int.max);
2405     assert(frexp(-0.0, exp) == -0.0 && exp == 0);
2406     assert(frexp(0.0, exp) == 0.0 && exp == 0);
2407 }
2408 
2409 @safe @nogc nothrow unittest
2410 {
2411     import std.math.operations : isClose;
2412 
2413     int exp;
2414     real mantissa = frexp(123.456L, exp);
2415 
2416     // check if values are equal to 19 decimal digits of precision
2417     assert(isClose(mantissa * pow(2.0L, cast(real) exp), 123.456L, 1e-18));
2418 }
2419 
2420 @safe unittest
2421 {
2422     import std.math : floatTraits, RealFormat;
2423     import std.math.traits : isIdentical;
2424     import std.meta : AliasSeq;
2425     import std.typecons : tuple, Tuple;
2426 
2427     static foreach (T; AliasSeq!(real, double, float))
2428     {{
2429         Tuple!(T, T, int)[] vals = [   // x,frexp,exp
2430             tuple(T(0.0),  T( 0.0 ), 0),
2431             tuple(T(-0.0), T( -0.0), 0),
2432             tuple(T(1.0),  T( .5  ), 1),
2433             tuple(T(-1.0), T( -.5 ), 1),
2434             tuple(T(2.0),  T( .5  ), 2),
2435             tuple(T(float.min_normal/2.0f), T(.5), -126),
2436             tuple(T.infinity, T.infinity, int.max),
2437             tuple(-T.infinity, -T.infinity, int.min),
2438             tuple(T.nan, T.nan, int.min),
2439             tuple(-T.nan, -T.nan, int.min),
2440 
2441             // https://issues.dlang.org/show_bug.cgi?id=16026:
2442             tuple(3 * (T.min_normal * T.epsilon), T( .75), (T.min_exp - T.mant_dig) + 2)
2443         ];
2444 
2445         foreach (elem; vals)
2446         {
2447             T x = elem[0];
2448             T e = elem[1];
2449             int exp = elem[2];
2450             int eptr;
2451             T v = frexp(x, eptr);
2452             assert(isIdentical(e, v));
2453             assert(exp == eptr);
2454         }
2455 
2456         static if (floatTraits!(T).realFormat == RealFormat.ieeeExtended)
2457         {
2458             static T[3][] extendedvals = [ // x,frexp,exp
2459                 [0x1.a5f1c2eb3fe4efp+73L,    0x1.A5F1C2EB3FE4EFp-1L,     74],    // normal
2460                 [0x1.fa01712e8f0471ap-1064L, 0x1.fa01712e8f0471ap-1L, -1063],
2461                 [T.min_normal,      .5, -16381],
2462                 [T.min_normal/2.0L, .5, -16382]    // subnormal
2463             ];
2464             foreach (elem; extendedvals)
2465             {
2466                 T x = elem[0];
2467                 T e = elem[1];
2468                 int exp = cast(int) elem[2];
2469                 int eptr;
2470                 T v = frexp(x, eptr);
2471                 assert(isIdentical(e, v));
2472                 assert(exp == eptr);
2473             }
2474         }
2475     }}
2476 
2477     // CTFE
2478     alias CtfeFrexpResult= Tuple!(real, int);
2479     static CtfeFrexpResult ctfeFrexp(T)(const T value)
2480     {
2481         int exp;
2482         auto significand = frexp(value, exp);
2483         return CtfeFrexpResult(significand, exp);
2484     }
2485     static foreach (T; AliasSeq!(real, double, float))
2486     {{
2487         enum Tuple!(T, T, int)[] vals = [   // x,frexp,exp
2488             tuple(T(0.0),  T( 0.0 ), 0),
2489             tuple(T(-0.0), T( -0.0), 0),
2490             tuple(T(1.0),  T( .5  ), 1),
2491             tuple(T(-1.0), T( -.5 ), 1),
2492             tuple(T(2.0),  T( .5  ), 2),
2493             tuple(T(float.min_normal/2.0f), T(.5), -126),
2494             tuple(T.infinity, T.infinity, int.max),
2495             tuple(-T.infinity, -T.infinity, int.min),
2496             tuple(T.nan, T.nan, int.min),
2497             tuple(-T.nan, -T.nan, int.min),
2498 
2499             // https://issues.dlang.org/show_bug.cgi?id=16026:
2500             tuple(3 * (T.min_normal * T.epsilon), T( .75), (T.min_exp - T.mant_dig) + 2)
2501         ];
2502 
2503         static foreach (elem; vals)
2504         {
2505             static assert(ctfeFrexp(elem[0]) is CtfeFrexpResult(elem[1], elem[2]));
2506         }
2507 
2508         static if (floatTraits!(T).realFormat == RealFormat.ieeeExtended)
2509         {
2510             enum T[3][] extendedvals = [ // x,frexp,exp
2511                 [0x1.a5f1c2eb3fe4efp+73L,    0x1.A5F1C2EB3FE4EFp-1L,     74],    // normal
2512                 [0x1.fa01712e8f0471ap-1064L, 0x1.fa01712e8f0471ap-1L, -1063],
2513                 [T.min_normal,      .5, -16381],
2514                 [T.min_normal/2.0L, .5, -16382]    // subnormal
2515             ];
2516             static foreach (elem; extendedvals)
2517             {
2518                 static assert(ctfeFrexp(elem[0]) is CtfeFrexpResult(elem[1], cast(int) elem[2]));
2519             }
2520         }
2521     }}
2522 }
2523 
2524 @safe unittest
2525 {
2526     import std.meta : AliasSeq;
2527     void foo() {
2528         static foreach (T; AliasSeq!(real, double, float))
2529         {{
2530             int exp;
2531             const T a = 1;
2532             immutable T b = 2;
2533             auto c = frexp(a, exp);
2534             auto d = frexp(b, exp);
2535         }}
2536     }
2537 }
2538 
2539 /******************************************
2540  * Extracts the exponent of x as a signed integral value.
2541  *
2542  * If x is not a special value, the result is the same as
2543  * $(D cast(int) logb(x)).
2544  *
2545  *      $(TABLE_SV
2546  *      $(TR $(TH x)                $(TH ilogb(x))     $(TH Range error?))
2547  *      $(TR $(TD 0)                 $(TD FP_ILOGB0)   $(TD yes))
2548  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD int.max)     $(TD no))
2549  *      $(TR $(TD $(NAN))            $(TD FP_ILOGBNAN) $(TD no))
2550  *      )
2551  */
2552 int ilogb(T)(const T x) @trusted pure nothrow @nogc
2553 if (isFloatingPoint!T)
2554 {
2555     import std.math : floatTraits, RealFormat, MANTISSA_MSB, MANTISSA_LSB;
2556 
2557     import core.bitop : bsr;
2558     alias F = floatTraits!T;
2559 
2560     union floatBits
2561     {
2562         T rv;
2563         ushort[T.sizeof/2] vu;
2564         uint[T.sizeof/4] vui;
2565         static if (T.sizeof >= 8)
2566             ulong[T.sizeof/8] vul;
2567     }
2568     floatBits y = void;
2569     y.rv = x;
2570 
2571     int ex = y.vu[F.EXPPOS_SHORT] & F.EXPMASK;
2572     static if (F.realFormat == RealFormat.ieeeExtended ||
2573                F.realFormat == RealFormat.ieeeExtended53)
2574     {
2575         if (ex)
2576         {
2577             // If exponent is non-zero
2578             if (ex == F.EXPMASK) // infinity or NaN
2579             {
2580                 if (y.vul[0] &  0x7FFF_FFFF_FFFF_FFFF)  // NaN
2581                     return FP_ILOGBNAN;
2582                 else // +-infinity
2583                     return int.max;
2584             }
2585             else
2586             {
2587                 return ex - F.EXPBIAS - 1;
2588             }
2589         }
2590         else if (!y.vul[0])
2591         {
2592             // vf is +-0.0
2593             return FP_ILOGB0;
2594         }
2595         else
2596         {
2597             // subnormal
2598             return ex - F.EXPBIAS - T.mant_dig + 1 + bsr(y.vul[0]);
2599         }
2600     }
2601     else static if (F.realFormat == RealFormat.ieeeQuadruple)
2602     {
2603         if (ex)    // If exponent is non-zero
2604         {
2605             if (ex == F.EXPMASK)
2606             {
2607                 // infinity or NaN
2608                 if (y.vul[MANTISSA_LSB] | ( y.vul[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF))  // NaN
2609                     return FP_ILOGBNAN;
2610                 else // +- infinity
2611                     return int.max;
2612             }
2613             else
2614             {
2615                 return ex - F.EXPBIAS - 1;
2616             }
2617         }
2618         else if ((y.vul[MANTISSA_LSB] | (y.vul[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) == 0)
2619         {
2620             // vf is +-0.0
2621             return FP_ILOGB0;
2622         }
2623         else
2624         {
2625             // subnormal
2626             const ulong msb = y.vul[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF;
2627             const ulong lsb = y.vul[MANTISSA_LSB];
2628             if (msb)
2629                 return ex - F.EXPBIAS - T.mant_dig + 1 + bsr(msb) + 64;
2630             else
2631                 return ex - F.EXPBIAS - T.mant_dig + 1 + bsr(lsb);
2632         }
2633     }
2634     else static if (F.realFormat == RealFormat.ieeeDouble)
2635     {
2636         if (ex) // If exponent is non-zero
2637         {
2638             if (ex == F.EXPMASK)   // infinity or NaN
2639             {
2640                 if ((y.vul[0] & 0x7FFF_FFFF_FFFF_FFFF) == 0x7FF0_0000_0000_0000)  // +- infinity
2641                     return int.max;
2642                 else // NaN
2643                     return FP_ILOGBNAN;
2644             }
2645             else
2646             {
2647                 return ((ex - F.EXPBIAS) >> 4) - 1;
2648             }
2649         }
2650         else if (!(y.vul[0] & 0x7FFF_FFFF_FFFF_FFFF))
2651         {
2652             // vf is +-0.0
2653             return FP_ILOGB0;
2654         }
2655         else
2656         {
2657             // subnormal
2658             enum MANTISSAMASK_64 = ((cast(ulong) F.MANTISSAMASK_INT) << 32) | 0xFFFF_FFFF;
2659             return ((ex - F.EXPBIAS) >> 4) - T.mant_dig + 1 + bsr(y.vul[0] & MANTISSAMASK_64);
2660         }
2661     }
2662     else static if (F.realFormat == RealFormat.ieeeSingle)
2663     {
2664         if (ex) // If exponent is non-zero
2665         {
2666             if (ex == F.EXPMASK)   // infinity or NaN
2667             {
2668                 if ((y.vui[0] & 0x7FFF_FFFF) == 0x7F80_0000)  // +- infinity
2669                     return int.max;
2670                 else // NaN
2671                     return FP_ILOGBNAN;
2672             }
2673             else
2674             {
2675                 return ((ex - F.EXPBIAS) >> 7) - 1;
2676             }
2677         }
2678         else if (!(y.vui[0] & 0x7FFF_FFFF))
2679         {
2680             // vf is +-0.0
2681             return FP_ILOGB0;
2682         }
2683         else
2684         {
2685             // subnormal
2686             const uint mantissa = y.vui[0] & F.MANTISSAMASK_INT;
2687             return ((ex - F.EXPBIAS) >> 7) - T.mant_dig + 1 + bsr(mantissa);
2688         }
2689     }
2690     else // static if (F.realFormat == RealFormat.ibmExtended)
2691     {
2692         assert(0, "ilogb not implemented");
2693     }
2694 }
2695 /// ditto
2696 int ilogb(T)(const T x) @safe pure nothrow @nogc
2697 if (isIntegral!T && isUnsigned!T)
2698 {
2699     import core.bitop : bsr;
2700     if (x == 0)
2701         return FP_ILOGB0;
2702     else
2703     {
2704         static assert(T.sizeof <= ulong.sizeof, "integer size too large for the current ilogb implementation");
2705         return bsr(x);
2706     }
2707 }
2708 /// ditto
2709 int ilogb(T)(const T x) @safe pure nothrow @nogc
2710 if (isIntegral!T && isSigned!T)
2711 {
2712     import std.traits : Unsigned;
2713     // Note: abs(x) can not be used because the return type is not Unsigned and
2714     //       the return value would be wrong for x == int.min
2715     Unsigned!T absx =  x >= 0 ? x : -x;
2716     return ilogb(absx);
2717 }
2718 
2719 ///
2720 @safe pure unittest
2721 {
2722     assert(ilogb(1) == 0);
2723     assert(ilogb(3) == 1);
2724     assert(ilogb(3.0) == 1);
2725     assert(ilogb(100_000_000) == 26);
2726 
2727     assert(ilogb(0) == FP_ILOGB0);
2728     assert(ilogb(0.0) == FP_ILOGB0);
2729     assert(ilogb(double.nan) == FP_ILOGBNAN);
2730     assert(ilogb(double.infinity) == int.max);
2731 }
2732 
2733 /**
2734 Special return values of $(LREF ilogb).
2735  */
2736 alias FP_ILOGB0   = core.stdc.math.FP_ILOGB0;
2737 /// ditto
2738 alias FP_ILOGBNAN = core.stdc.math.FP_ILOGBNAN;
2739 
2740 ///
2741 @safe pure unittest
2742 {
2743     assert(ilogb(0) == FP_ILOGB0);
2744     assert(ilogb(0.0) == FP_ILOGB0);
2745     assert(ilogb(double.nan) == FP_ILOGBNAN);
2746 }
2747 
2748 @safe nothrow @nogc unittest
2749 {
2750     import std.math : floatTraits, RealFormat;
2751     import std.math.operations : nextUp;
2752     import std.meta : AliasSeq;
2753     import std.typecons : Tuple;
2754     static foreach (F; AliasSeq!(float, double, real))
2755     {{
2756         alias T = Tuple!(F, int);
2757         T[13] vals =   // x, ilogb(x)
2758         [
2759             T(  F.nan     , FP_ILOGBNAN ),
2760             T( -F.nan     , FP_ILOGBNAN ),
2761             T(  F.infinity, int.max     ),
2762             T( -F.infinity, int.max     ),
2763             T(  0.0       , FP_ILOGB0   ),
2764             T( -0.0       , FP_ILOGB0   ),
2765             T(  2.0       , 1           ),
2766             T(  2.0001    , 1           ),
2767             T(  1.9999    , 0           ),
2768             T(  0.5       , -1          ),
2769             T(  123.123   , 6           ),
2770             T( -123.123   , 6           ),
2771             T(  0.123     , -4          ),
2772         ];
2773 
2774         foreach (elem; vals)
2775         {
2776             assert(ilogb(elem[0]) == elem[1]);
2777         }
2778     }}
2779 
2780     // min_normal and subnormals
2781     assert(ilogb(-float.min_normal) == -126);
2782     assert(ilogb(nextUp(-float.min_normal)) == -127);
2783     assert(ilogb(nextUp(-float(0.0))) == -149);
2784     assert(ilogb(-double.min_normal) == -1022);
2785     assert(ilogb(nextUp(-double.min_normal)) == -1023);
2786     assert(ilogb(nextUp(-double(0.0))) == -1074);
2787     static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended)
2788     {
2789         assert(ilogb(-real.min_normal) == -16382);
2790         assert(ilogb(nextUp(-real.min_normal)) == -16383);
2791         assert(ilogb(nextUp(-real(0.0))) == -16445);
2792     }
2793     else static if (floatTraits!(real).realFormat == RealFormat.ieeeDouble)
2794     {
2795         assert(ilogb(-real.min_normal) == -1022);
2796         assert(ilogb(nextUp(-real.min_normal)) == -1023);
2797         assert(ilogb(nextUp(-real(0.0))) == -1074);
2798     }
2799 
2800     // test integer types
2801     assert(ilogb(0) == FP_ILOGB0);
2802     assert(ilogb(int.max) == 30);
2803     assert(ilogb(int.min) == 31);
2804     assert(ilogb(uint.max) == 31);
2805     assert(ilogb(long.max) == 62);
2806     assert(ilogb(long.min) == 63);
2807     assert(ilogb(ulong.max) == 63);
2808 }
2809 
2810 /*******************************************
2811  * Compute n * 2$(SUPERSCRIPT exp)
2812  * References: frexp
2813  */
2814 
2815 pragma(inline, true)
2816 real ldexp(real n, int exp)     @safe pure nothrow @nogc { return core.math.ldexp(n, exp); }
2817 ///ditto
2818 pragma(inline, true)
2819 double ldexp(double n, int exp) @safe pure nothrow @nogc { return core.math.ldexp(n, exp); }
2820 ///ditto
2821 pragma(inline, true)
2822 float ldexp(float n, int exp)   @safe pure nothrow @nogc { return core.math.ldexp(n, exp); }
2823 
2824 ///
2825 @nogc @safe pure nothrow unittest
2826 {
2827     import std.meta : AliasSeq;
2828     static foreach (T; AliasSeq!(float, double, real))
2829     {{
2830         T r;
2831 
2832         r = ldexp(3.0L, 3);
2833         assert(r == 24);
2834 
2835         r = ldexp(cast(T) 3.0, cast(int) 3);
2836         assert(r == 24);
2837 
2838         T n = 3.0;
2839         int exp = 3;
2840         r = ldexp(n, exp);
2841         assert(r == 24);
2842     }}
2843 }
2844 
2845 @safe pure nothrow @nogc unittest
2846 {
2847     import std.math : floatTraits, RealFormat;
2848 
2849     static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended ||
2850                floatTraits!(real).realFormat == RealFormat.ieeeExtended53 ||
2851                floatTraits!(real).realFormat == RealFormat.ieeeQuadruple)
2852     {
2853         assert(ldexp(1.0L, -16384) == 0x1p-16384L);
2854         assert(ldexp(1.0L, -16382) == 0x1p-16382L);
2855         int x;
2856         real n = frexp(0x1p-16384L, x);
2857         assert(n == 0.5L);
2858         assert(x==-16383);
2859         assert(ldexp(n, x)==0x1p-16384L);
2860     }
2861     else static if (floatTraits!(real).realFormat == RealFormat.ieeeDouble)
2862     {
2863         assert(ldexp(1.0L, -1024) == 0x1p-1024L);
2864         assert(ldexp(1.0L, -1022) == 0x1p-1022L);
2865         int x;
2866         real n = frexp(0x1p-1024L, x);
2867         assert(n == 0.5L);
2868         assert(x==-1023);
2869         assert(ldexp(n, x)==0x1p-1024L);
2870     }
2871     else static assert(false, "Floating point type real not supported");
2872 }
2873 
2874 /* workaround https://issues.dlang.org/show_bug.cgi?id=14718
2875    float parsing depends on platform strtold
2876 @safe pure nothrow @nogc unittest
2877 {
2878     assert(ldexp(1.0, -1024) == 0x1p-1024);
2879     assert(ldexp(1.0, -1022) == 0x1p-1022);
2880     int x;
2881     double n = frexp(0x1p-1024, x);
2882     assert(n == 0.5);
2883     assert(x==-1023);
2884     assert(ldexp(n, x)==0x1p-1024);
2885 }
2886 
2887 @safe pure nothrow @nogc unittest
2888 {
2889     assert(ldexp(1.0f, -128) == 0x1p-128f);
2890     assert(ldexp(1.0f, -126) == 0x1p-126f);
2891     int x;
2892     float n = frexp(0x1p-128f, x);
2893     assert(n == 0.5f);
2894     assert(x==-127);
2895     assert(ldexp(n, x)==0x1p-128f);
2896 }
2897 */
2898 
2899 @safe @nogc nothrow unittest
2900 {
2901     import std.math.operations : isClose;
2902 
2903     static real[3][] vals =    // value,exp,ldexp
2904     [
2905     [    0,    0,    0],
2906     [    1,    0,    1],
2907     [    -1,    0,    -1],
2908     [    1,    1,    2],
2909     [    123,    10,    125952],
2910     [    real.max,    int.max,    real.infinity],
2911     [    real.max,    -int.max,    0],
2912     [    real.min_normal,    -int.max,    0],
2913     ];
2914     int i;
2915 
2916     for (i = 0; i < vals.length; i++)
2917     {
2918         real x = vals[i][0];
2919         int exp = cast(int) vals[i][1];
2920         real z = vals[i][2];
2921         real l = ldexp(x, exp);
2922 
2923         assert(isClose(z, l, 1e-6));
2924     }
2925 
2926     real function(real, int) pldexp = &ldexp;
2927     assert(pldexp != null);
2928 }
2929 
2930 private
2931 {
2932     // Coefficients shared across log(), log2(), log10(), log1p().
2933     template LogCoeffs(T)
2934     {
2935         import std.math : floatTraits, RealFormat;
2936 
2937         static if (floatTraits!T.realFormat == RealFormat.ieeeQuadruple)
2938         {
2939             // Coefficients for log(1 + x) = x - x^^2/2 + x^^3 P(x)/Q(x)
2940             // Theoretical peak relative error = 5.3e-37
2941             static immutable real[13] logP = [
2942                 1.313572404063446165910279910527789794488E4L,
2943                 7.771154681358524243729929227226708890930E4L,
2944                 2.014652742082537582487669938141683759923E5L,
2945                 3.007007295140399532324943111654767187848E5L,
2946                 2.854829159639697837788887080758954924001E5L,
2947                 1.797628303815655343403735250238293741397E5L,
2948                 7.594356839258970405033155585486712125861E4L,
2949                 2.128857716871515081352991964243375186031E4L,
2950                 3.824952356185897735160588078446136783779E3L,
2951                 4.114517881637811823002128927449878962058E2L,
2952                 2.321125933898420063925789532045674660756E1L,
2953                 4.998469661968096229986658302195402690910E-1L,
2954                 1.538612243596254322971797716843006400388E-6L
2955             ];
2956             static immutable real[13] logQ = [
2957                 3.940717212190338497730839731583397586124E4L,
2958                 2.626900195321832660448791748036714883242E5L,
2959                 7.777690340007566932935753241556479363645E5L,
2960                 1.347518538384329112529391120390701166528E6L,
2961                 1.514882452993549494932585972882995548426E6L,
2962                 1.158019977462989115839826904108208787040E6L,
2963                 6.132189329546557743179177159925690841200E5L,
2964                 2.248234257620569139969141618556349415120E5L,
2965                 5.605842085972455027590989944010492125825E4L,
2966                 9.147150349299596453976674231612674085381E3L,
2967                 9.104928120962988414618126155557301584078E2L,
2968                 4.839208193348159620282142911143429644326E1L,
2969                 1.0
2970             ];
2971 
2972             // log2 uses the same coefficients as log.
2973             alias log2P = logP;
2974             alias log2Q = logQ;
2975 
2976             // log10 uses the same coefficients as log.
2977             alias log10P = logP;
2978             alias log10Q = logQ;
2979 
2980             // Coefficients for log(x) = z + z^^3 R(z^^2)/S(z^^2)
2981             // where z = 2(x-1)/(x+1)
2982             // Theoretical peak relative error = 1.1e-35
2983             static immutable real[6] logR = [
2984                 1.418134209872192732479751274970992665513E5L,
2985                 -8.977257995689735303686582344659576526998E4L,
2986                 2.048819892795278657810231591630928516206E4L,
2987                 -2.024301798136027039250415126250455056397E3L,
2988                 8.057002716646055371965756206836056074715E1L,
2989                 -8.828896441624934385266096344596648080902E-1L
2990             ];
2991             static immutable real[7] logS = [
2992                 1.701761051846631278975701529965589676574E6L,
2993                 -1.332535117259762928288745111081235577029E6L,
2994                 4.001557694070773974936904547424676279307E5L,
2995                 -5.748542087379434595104154610899551484314E4L,
2996                 3.998526750980007367835804959888064681098E3L,
2997                 -1.186359407982897997337150403816839480438E2L,
2998                 1.0
2999             ];
3000         }
3001         else static if (floatTraits!T.realFormat == RealFormat.ieeeExtended ||
3002                         floatTraits!T.realFormat == RealFormat.ieeeExtended53)
3003         {
3004             // Coefficients for log(1 + x) = x - x^^2/2 + x^^3 P(x)/Q(x)
3005             // Theoretical peak relative error = 2.32e-20
3006             static immutable real[7] logP = [
3007                 2.0039553499201281259648E1L,
3008                 5.7112963590585538103336E1L,
3009                 6.0949667980987787057556E1L,
3010                 2.9911919328553073277375E1L,
3011                 6.5787325942061044846969E0L,
3012                 4.9854102823193375972212E-1L,
3013                 4.5270000862445199635215E-5L,
3014             ];
3015             static immutable real[7] logQ = [
3016                 6.0118660497603843919306E1L,
3017                 2.1642788614495947685003E2L,
3018                 3.0909872225312059774938E2L,
3019                 2.2176239823732856465394E2L,
3020                 8.3047565967967209469434E1L,
3021                 1.5062909083469192043167E1L,
3022                 1.0000000000000000000000E0L,
3023             ];
3024 
3025             // Coefficients for log(1 + x) = x - x^^2/2 + x^^3 P(x)/Q(x)
3026             // Theoretical peak relative error = 6.2e-22
3027             static immutable real[7] log2P = [
3028                 1.0747524399916215149070E2L,
3029                 3.4258224542413922935104E2L,
3030                 4.2401812743503691187826E2L,
3031                 2.5620629828144409632571E2L,
3032                 7.7671073698359539859595E1L,
3033                 1.0767376367209449010438E1L,
3034                 4.9962495940332550844739E-1L,
3035             ];
3036             static immutable real[8] log2Q = [
3037                 3.2242573199748645407652E2L,
3038                 1.2695660352705325274404E3L,
3039                 2.0307734695595183428202E3L,
3040                 1.6911722418503949084863E3L,
3041                 7.7952888181207260646090E2L,
3042                 1.9444210022760132894510E2L,
3043                 2.3479774160285863271658E1L,
3044                 1.0000000000000000000000E0,
3045             ];
3046 
3047             // log10 uses the same coefficients as log2.
3048             alias log10P = log2P;
3049             alias log10Q = log2Q;
3050 
3051             // Coefficients for log(x) = z + z^^3 R(z^^2)/S(z^^2)
3052             // where z = 2(x-1)/(x+1)
3053             // Theoretical peak relative error = 6.16e-22
3054             static immutable real[4] logR = [
3055                -3.5717684488096787370998E1L,
3056                 1.0777257190312272158094E1L,
3057                -7.1990767473014147232598E-1L,
3058                 1.9757429581415468984296E-3L,
3059             ];
3060             static immutable real[4] logS = [
3061                -4.2861221385716144629696E2L,
3062                 1.9361891836232102174846E2L,
3063                -2.6201045551331104417768E1L,
3064                 1.0000000000000000000000E0L,
3065             ];
3066         }
3067         else static if (floatTraits!T.realFormat == RealFormat.ieeeDouble)
3068         {
3069             // Coefficients for log(1 + x) = x - x^^2/2 + x^^3 P(x)/Q(x)
3070             static immutable double[6] logP = [
3071                 7.70838733755885391666E0,
3072                 1.79368678507819816313E1,
3073                 1.44989225341610930846E1,
3074                 4.70579119878881725854E0,
3075                 4.97494994976747001425E-1,
3076                 1.01875663804580931796E-4,
3077             ];
3078             static immutable double[6] logQ = [
3079                 2.31251620126765340583E1,
3080                 7.11544750618563894466E1,
3081                 8.29875266912776603211E1,
3082                 4.52279145837532221105E1,
3083                 1.12873587189167450590E1,
3084                 1.00000000000000000000E0,
3085             ];
3086 
3087             // log2 uses the same coefficients as log.
3088             alias log2P = logP;
3089             alias log2Q = logQ;
3090 
3091             // Coefficients for log(1 + x) = x - x^^2/2 + x^^3 P(x)/Q(x)
3092             static immutable double[7] logp1P = [
3093                 2.0039553499201281259648E1,
3094                 5.7112963590585538103336E1,
3095                 6.0949667980987787057556E1,
3096                 2.9911919328553073277375E1,
3097                 6.5787325942061044846969E0,
3098                 4.9854102823193375972212E-1,
3099                 4.5270000862445199635215E-5,
3100             ];
3101             static immutable double[7] logp1Q = [
3102                 6.0118660497603843919306E1,
3103                 2.1642788614495947685003E2,
3104                 3.0909872225312059774938E2,
3105                 2.2176239823732856465394E2,
3106                 8.3047565967967209469434E1,
3107                 1.5062909083469192043167E1,
3108                 1.0000000000000000000000E0,
3109             ];
3110 
3111             static immutable double[7] log10P = [
3112                 1.98892446572874072159E1,
3113                 5.67349287391754285487E1,
3114                 6.06127134467767258030E1,
3115                 2.97877425097986925891E1,
3116                 6.56312093769992875930E0,
3117                 4.98531067254050724270E-1,
3118                 4.58482948458143443514E-5,
3119             ];
3120             static immutable double[7] log10Q = [
3121                 5.96677339718622216300E1,
3122                 2.14955586696422947765E2,
3123                 3.07254189979530058263E2,
3124                 2.20664384982121929218E2,
3125                 8.27410449222435217021E1,
3126                 1.50314182634250003249E1,
3127                 1.00000000000000000000E0,
3128             ];
3129 
3130             // Coefficients for log(x) = z + z^^3 R(z)/S(z)
3131             // where z = 2(x-1)/(x+1)
3132             static immutable double[3] logR = [
3133                 -6.41409952958715622951E1,
3134                 1.63866645699558079767E1,
3135                 -7.89580278884799154124E-1,
3136             ];
3137             static immutable double[4] logS = [
3138                 -7.69691943550460008604E2,
3139                 3.12093766372244180303E2,
3140                 -3.56722798256324312549E1,
3141                 1.00000000000000000000E0,
3142             ];
3143         }
3144         else static if (floatTraits!T.realFormat == RealFormat.ieeeSingle)
3145         {
3146             // Coefficients for log(1 + x) = x - x^^2/2 + x^^3 P(x)
3147             static immutable float[9] logP = [
3148                  3.3333331174E-1,
3149                 -2.4999993993E-1,
3150                  2.0000714765E-1,
3151                 -1.6668057665E-1,
3152                  1.4249322787E-1,
3153                 -1.2420140846E-1,
3154                  1.1676998740E-1,
3155                 -1.1514610310E-1,
3156                  7.0376836292E-2,
3157             ];
3158 
3159             // Coefficients for log(1 + x) = x - x^^2/2 + x^^3 P(x)/Q(x)
3160             static immutable float[7] logp1P = [
3161                  2.0039553499E1,
3162                  5.7112963590E1,
3163                  6.0949667980E1,
3164                  2.9911919328E1,
3165                  6.5787325942E0,
3166                  4.9854102823E-1,
3167                  4.5270000862E-5,
3168             ];
3169             static immutable float[7] logp1Q = [
3170                 6.01186604976E1,
3171                 2.16427886144E2,
3172                 3.09098722253E2,
3173                 2.21762398237E2,
3174                 8.30475659679E1,
3175                 1.50629090834E1,
3176                 1.00000000000E0,
3177             ];
3178 
3179             // log2 and log10 uses the same coefficients as log.
3180             alias log2P = logP;
3181             alias log10P = logP;
3182         }
3183         else
3184             static assert(0, "no coefficients for log()");
3185     }
3186 }
3187 
3188 /**************************************
3189  * Calculate the natural logarithm of x.
3190  *
3191  *    $(TABLE_SV
3192  *    $(TR $(TH x)            $(TH log(x))    $(TH divide by 0?) $(TH invalid?))
3193  *    $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes)          $(TD no))
3194  *    $(TR $(TD $(LT)0.0)     $(TD $(NAN))    $(TD no)           $(TD yes))
3195  *    $(TR $(TD +$(INFIN))    $(TD +$(INFIN)) $(TD no)           $(TD no))
3196  *    )
3197  */
3198 pragma(inline, true)
3199 real log(real x) @safe pure nothrow @nogc
3200 {
3201     version (INLINE_YL2X)
3202     {
3203         import std.math.constants : LN2;
3204         return core.math.yl2x(x, LN2);
3205     }
3206     else
3207         return logImpl(x);
3208 }
3209 
3210 /// ditto
3211 pragma(inline, true)
3212 double log(double x) @safe pure nothrow @nogc { return __ctfe ? cast(double) log(cast(real) x) : logImpl(x); }
3213 
3214 /// ditto
3215 pragma(inline, true)
3216 float log(float x) @safe pure nothrow @nogc { return __ctfe ? cast(float) log(cast(real) x) : logImpl(x); }
3217 
3218 // @@@DEPRECATED_[2.112.0]@@@
3219 deprecated("`std.math.exponential.log` called with argument types `(int)` matches both "
3220            ~ "`log(real)`, `log(double)`, and `log(float)`. Cast argument to floating point type instead.")
3221 real log(int x) @safe pure nothrow @nogc { return log(cast(real) x); }
3222 // @@@DEPRECATED_[2.112.0]@@@
3223 deprecated("`std.math.exponential.log` called with argument types `(uint)` matches both "
3224            ~ "`log(real)`, `log(double)`, and `log(float)`. Cast argument to floating point type instead.")
3225 real log(uint x) @safe pure nothrow @nogc { return log(cast(real) x); }
3226 // @@@DEPRECATED_[2.112.0]@@@
3227 deprecated("`std.math.exponential.log` called with argument types `(long)` matches both "
3228            ~ "`log(real)`, `log(double)`, and `log(float)`. Cast argument to floating point type instead.")
3229 real log(long x) @safe pure nothrow @nogc { return log(cast(real) x); }
3230 // @@@DEPRECATED_[2.112.0]@@@
3231 deprecated("`std.math.exponential.log` called with argument types `(ulong)` matches both "
3232            ~ "`log(real)`, `log(double)`, and `log(float)`. Cast argument to floating point type instead.")
3233 real log(ulong x) @safe pure nothrow @nogc { return log(cast(real) x); }
3234 
3235 ///
3236 @safe pure nothrow @nogc unittest
3237 {
3238     import std.math.operations : feqrel;
3239     import std.math.constants : E;
3240 
3241     assert(feqrel(log(E), 1) >= real.mant_dig - 1);
3242 }
3243 
3244 private T logImpl(T, bool LOG1P = false)(T x) @safe pure nothrow @nogc
3245 {
3246     import std.math.constants : SQRT1_2;
3247     import std.math.algebraic : poly;
3248     import std.math.traits : isInfinity, isNaN, signbit;
3249     import std.math : floatTraits, RealFormat;
3250 
3251     alias coeffs = LogCoeffs!T;
3252     alias F = floatTraits!T;
3253 
3254     static if (LOG1P)
3255     {
3256         const T xm1 = x;
3257         x = x + 1.0;
3258     }
3259 
3260     static if (F.realFormat == RealFormat.ieeeExtended ||
3261                F.realFormat == RealFormat.ieeeExtended53 ||
3262                F.realFormat == RealFormat.ieeeQuadruple)
3263     {
3264         // C1 + C2 = LN2.
3265         enum T C1 = 6.93145751953125E-1L;
3266         enum T C2 = 1.428606820309417232121458176568075500134E-6L;
3267     }
3268     else static if (F.realFormat == RealFormat.ieeeDouble)
3269     {
3270         enum T C1 = 0.693359375;
3271         enum T C2 = -2.121944400546905827679e-4;
3272     }
3273     else static if (F.realFormat == RealFormat.ieeeSingle)
3274     {
3275         enum T C1 = 0.693359375;
3276         enum T C2 = -2.12194440e-4;
3277     }
3278     else
3279         static assert(0, "Not implemented for this architecture");
3280 
3281     // Special cases.
3282     if (isNaN(x))
3283         return x;
3284     if (isInfinity(x) && !signbit(x))
3285         return x;
3286     if (x == 0.0)
3287         return -T.infinity;
3288     if (x < 0.0)
3289         return T.nan;
3290 
3291     // Separate mantissa from exponent.
3292     // Note, frexp is used so that denormal numbers will be handled properly.
3293     T y, z;
3294     int exp;
3295 
3296     x = frexp(x, exp);
3297 
3298     static if (F.realFormat == RealFormat.ieeeDouble ||
3299                F.realFormat == RealFormat.ieeeExtended ||
3300                F.realFormat == RealFormat.ieeeExtended53 ||
3301                F.realFormat == RealFormat.ieeeQuadruple)
3302     {
3303         // Logarithm using log(x) = z + z^^3 R(z) / S(z),
3304         // where z = 2(x - 1)/(x + 1)
3305         if ((exp > 2) || (exp < -2))
3306         {
3307             if (x < SQRT1_2)
3308             {   // 2(2x - 1)/(2x + 1)
3309                 exp -= 1;
3310                 z = x - 0.5;
3311                 y = 0.5 * z + 0.5;
3312             }
3313             else
3314             {   // 2(x - 1)/(x + 1)
3315                 z = x - 0.5;
3316                 z -= 0.5;
3317                 y = 0.5 * x  + 0.5;
3318             }
3319             x = z / y;
3320             z = x * x;
3321             z = x * (z * poly(z, coeffs.logR) / poly(z, coeffs.logS));
3322             z += exp * C2;
3323             z += x;
3324             z += exp * C1;
3325 
3326             return z;
3327         }
3328     }
3329 
3330     // Logarithm using log(1 + x) = x - .5x^^2 + x^^3 P(x) / Q(x)
3331     if (x < SQRT1_2)
3332     {
3333         exp -= 1;
3334         static if (LOG1P)
3335         {
3336             if (exp != 0)
3337                 x = 2.0 * x - 1.0;
3338             else
3339                 x = xm1;
3340         }
3341         else
3342             x = 2.0 * x - 1.0;
3343 
3344     }
3345     else
3346     {
3347         static if (LOG1P)
3348         {
3349             if (exp != 0)
3350                 x = x - 1.0;
3351             else
3352                 x = xm1;
3353         }
3354         else
3355             x = x - 1.0;
3356     }
3357     z = x * x;
3358     static if (F.realFormat == RealFormat.ieeeSingle)
3359         y = x * (z * poly(x, coeffs.logP));
3360     else
3361         y = x * (z * poly(x, coeffs.logP) / poly(x, coeffs.logQ));
3362     y += exp * C2;
3363     z = y - 0.5 * z;
3364 
3365     // Note, the sum of above terms does not exceed x/4,
3366     // so it contributes at most about 1/4 lsb to the error.
3367     z += x;
3368     z += exp * C1;
3369 
3370     return z;
3371 }
3372 
3373 @safe @nogc nothrow unittest
3374 {
3375     import std.math : floatTraits, RealFormat;
3376     import std.meta : AliasSeq;
3377 
3378     static void testLog(T)(T[2][] vals)
3379     {
3380         import std.math.operations : isClose;
3381         import std.math.traits : isNaN;
3382         foreach (ref pair; vals)
3383         {
3384             if (isNaN(pair[1]))
3385                 assert(isNaN(log(pair[0])));
3386             else
3387                 assert(isClose(log(pair[0]), pair[1]));
3388         }
3389     }
3390     static foreach (F; AliasSeq!(float, double, real))
3391     {{
3392         scope F[2][] vals = [
3393             [F(1), F(0x0p+0)], [F(2), F(0x1.62e42fefa39ef358p-1)],
3394             [F(4), F(0x1.62e42fefa39ef358p+0)], [F(8), F(0x1.0a2b23f3bab73682p+1)],
3395             [F(16), F(0x1.62e42fefa39ef358p+1)], [F(32), F(0x1.bb9d3beb8c86b02ep+1)],
3396             [F(64), F(0x1.0a2b23f3bab73682p+2)], [F(128), F(0x1.3687a9f1af2b14ecp+2)],
3397             [F(256), F(0x1.62e42fefa39ef358p+2)], [F(512), F(0x1.8f40b5ed9812d1c2p+2)],
3398             [F(1024), F(0x1.bb9d3beb8c86b02ep+2)], [F(2048), F(0x1.e7f9c1e980fa8e98p+2)],
3399             [F(3), F(0x1.193ea7aad030a976p+0)], [F(5), F(0x1.9c041f7ed8d336bp+0)],
3400             [F(7), F(0x1.f2272ae325a57546p+0)], [F(15), F(0x1.5aa16394d481f014p+1)],
3401             [F(17), F(0x1.6aa6bc1fa7f79cfp+1)], [F(31), F(0x1.b78ce48912b59f12p+1)],
3402             [F(33), F(0x1.bf8d8f4d5b8d1038p+1)], [F(63), F(0x1.09291e8e3181b20ep+2)],
3403             [F(65), F(0x1.0b292939429755ap+2)], [F(-0), -F.infinity], [F(0), -F.infinity],
3404             [F(0.1), F(-0x1.26bb1bbb5551582ep+1)], [F(0.25), F(-0x1.62e42fefa39ef358p+0)],
3405             [F(0.75), F(-0x1.269621134db92784p-2)], [F(0.875), F(-0x1.1178e8227e47bde4p-3)],
3406             [F(10), F(0x1.26bb1bbb5551582ep+1)], [F(100), F(0x1.26bb1bbb5551582ep+2)],
3407             [F(10000), F(0x1.26bb1bbb5551582ep+3)],
3408         ];
3409         testLog(vals);
3410     }}
3411     {
3412         float[2][16] vals = [
3413             [float.nan, float.nan],[-float.nan, float.nan],
3414             [float.infinity, float.infinity], [-float.infinity, float.nan],
3415             [float.min_normal, -0x1.5d58ap+6f], [-float.min_normal, float.nan],
3416             [float.max, 0x1.62e43p+6f], [-float.max, float.nan],
3417             [float.min_normal / 2, -0x1.601e68p+6f], [-float.min_normal / 2, float.nan],
3418             [float.max / 2, 0x1.601e68p+6f], [-float.max / 2, float.nan],
3419             [float.min_normal / 3, -0x1.61bd9ap+6f], [-float.min_normal / 3, float.nan],
3420             [float.max / 3, 0x1.5e7f36p+6f], [-float.max / 3, float.nan],
3421         ];
3422         testLog(vals);
3423     }
3424     {
3425         double[2][16] vals = [
3426             [double.nan, double.nan],[-double.nan, double.nan],
3427             [double.infinity, double.infinity], [-double.infinity, double.nan],
3428             [double.min_normal, -0x1.6232bdd7abcd2p+9], [-double.min_normal, double.nan],
3429             [double.max, 0x1.62e42fefa39efp+9], [-double.max, double.nan],
3430             [double.min_normal / 2, -0x1.628b76e3a7b61p+9], [-double.min_normal / 2, double.nan],
3431             [double.max / 2, 0x1.628b76e3a7b61p+9], [-double.max / 2, double.nan],
3432             [double.min_normal / 3, -0x1.62bf5d2b81354p+9], [-double.min_normal / 3, double.nan],
3433             [double.max / 3, 0x1.6257909bce36ep+9], [-double.max / 3, double.nan],
3434         ];
3435         testLog(vals);
3436     }
3437     alias F = floatTraits!real;
3438     static if (F.realFormat == RealFormat.ieeeExtended || F.realFormat == RealFormat.ieeeQuadruple)
3439     {{
3440         real[2][16] vals = [
3441             [real.nan, real.nan],[-real.nan, real.nan],
3442             [real.infinity, real.infinity], [-real.infinity, real.nan],
3443             [real.min_normal, -0x1.62d918ce2421d66p+13L], [-real.min_normal, real.nan],
3444             [real.max, 0x1.62e42fefa39ef358p+13L], [-real.max, real.nan],
3445             [real.min_normal / 2, -0x1.62dea45ee3e064dcp+13L], [-real.min_normal / 2, real.nan],
3446             [real.max / 2, 0x1.62dea45ee3e064dcp+13L], [-real.max / 2, real.nan],
3447             [real.min_normal / 3, -0x1.62e1e2c3617857e6p+13L], [-real.min_normal / 3, real.nan],
3448             [real.max / 3, 0x1.62db65fa664871d2p+13L], [-real.max / 3, real.nan],
3449         ];
3450         testLog(vals);
3451     }}
3452 }
3453 
3454 /**************************************
3455  * Calculate the base-10 logarithm of x.
3456  *
3457  *      $(TABLE_SV
3458  *      $(TR $(TH x)            $(TH log10(x))  $(TH divide by 0?) $(TH invalid?))
3459  *      $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes)          $(TD no))
3460  *      $(TR $(TD $(LT)0.0)     $(TD $(NAN))    $(TD no)           $(TD yes))
3461  *      $(TR $(TD +$(INFIN))    $(TD +$(INFIN)) $(TD no)           $(TD no))
3462  *      )
3463  */
3464 pragma(inline, true)
3465 real log10(real x) @safe pure nothrow @nogc
3466 {
3467     version (INLINE_YL2X)
3468     {
3469         import std.math.constants : LOG2;
3470         return core.math.yl2x(x, LOG2);
3471     }
3472     else
3473         return log10Impl(x);
3474 }
3475 
3476 /// ditto
3477 pragma(inline, true)
3478 double log10(double x) @safe pure nothrow @nogc { return __ctfe ? cast(double) log10(cast(real) x) : log10Impl(x); }
3479 
3480 /// ditto
3481 pragma(inline, true)
3482 float log10(float x) @safe pure nothrow @nogc { return __ctfe ? cast(float) log10(cast(real) x) : log10Impl(x); }
3483 
3484 // @@@DEPRECATED_[2.112.0]@@@
3485 deprecated("`std.math.exponential.log10` called with argument types `(int)` matches both "
3486            ~ "`log10(real)`, `log10(double)`, and `log10(float)`. Cast argument to floating point type instead.")
3487 real log10(int x) @safe pure nothrow @nogc { return log10(cast(real) x); }
3488 // @@@DEPRECATED_[2.112.0]@@@
3489 deprecated("`std.math.exponential.log10` called with argument types `(uint)` matches both "
3490            ~ "`log10(real)`, `log10(double)`, and `log10(float)`. Cast argument to floating point type instead.")
3491 real log10(uint x) @safe pure nothrow @nogc { return log10(cast(real) x); }
3492 // @@@DEPRECATED_[2.112.0]@@@
3493 deprecated("`std.math.exponential.log10` called with argument types `(long)` matches both "
3494            ~ "`log10(real)`, `log10(double)`, and `log10(float)`. Cast argument to floating point type instead.")
3495 real log10(long x) @safe pure nothrow @nogc { return log10(cast(real) x); }
3496 // @@@DEPRECATED_[2.112.0]@@@
3497 deprecated("`std.math.exponential.log10` called with argument types `(ulong)` matches both "
3498            ~ "`log10(real)`, `log10(double)`, and `log10(float)`. Cast argument to floating point type instead.")
3499 real log10(ulong x) @safe pure nothrow @nogc { return log10(cast(real) x); }
3500 
3501 ///
3502 @safe pure nothrow @nogc unittest
3503 {
3504     import std.math.algebraic : fabs;
3505 
3506     assert(fabs(log10(1000.0L) - 3) < .000001);
3507 }
3508 
3509 @safe pure nothrow @nogc unittest
3510 {
3511     import std.math.algebraic : fabs;
3512 
3513     assert(fabs(log10(1000.0) - 3) < .000001);
3514     assert(fabs(log10(1000.0f) - 3) < .000001);
3515 }
3516 
3517 private T log10Impl(T)(T x) @safe pure nothrow @nogc
3518 {
3519     import std.math.constants : SQRT1_2;
3520     import std.math.algebraic : poly;
3521     import std.math.traits : isNaN, isInfinity, signbit;
3522     import std.math : floatTraits, RealFormat;
3523 
3524     alias coeffs = LogCoeffs!T;
3525     alias F = floatTraits!T;
3526 
3527     static if (F.realFormat == RealFormat.ieeeExtended ||
3528                F.realFormat == RealFormat.ieeeExtended53 ||
3529                F.realFormat == RealFormat.ieeeQuadruple)
3530     {
3531         // log10(2) split into two parts.
3532         enum T L102A =  0.3125L;
3533         enum T L102B = -1.14700043360188047862611052755069732318101185E-2L;
3534 
3535         // log10(e) split into two parts.
3536         enum T L10EA =  0.5L;
3537         enum T L10EB = -6.570551809674817234887108108339491770560299E-2L;
3538     }
3539     else static if (F.realFormat == RealFormat.ieeeDouble ||
3540                     F.realFormat == RealFormat.ieeeSingle)
3541     {
3542         enum T L102A =  3.0078125E-1;
3543         enum T L102B = 2.48745663981195213739E-4;
3544 
3545         enum T L10EA =  4.3359375E-1;
3546         enum T L10EB = 7.00731903251827651129E-4;
3547     }
3548     else
3549         static assert(0, "Not implemented for this architecture");
3550 
3551     // Special cases are the same as for log.
3552     if (isNaN(x))
3553         return x;
3554     if (isInfinity(x) && !signbit(x))
3555         return x;
3556     if (x == 0.0)
3557         return -T.infinity;
3558     if (x < 0.0)
3559         return T.nan;
3560 
3561     // Separate mantissa from exponent.
3562     // Note, frexp is used so that denormal numbers will be handled properly.
3563     T y, z;
3564     int exp;
3565 
3566     x = frexp(x, exp);
3567 
3568     static if (F.realFormat == RealFormat.ieeeExtended ||
3569                F.realFormat == RealFormat.ieeeExtended53 ||
3570                F.realFormat == RealFormat.ieeeQuadruple)
3571     {
3572         // Logarithm using log(x) = z + z^^3 R(z) / S(z),
3573         // where z = 2(x - 1)/(x + 1)
3574         if ((exp > 2) || (exp < -2))
3575         {
3576             if (x < SQRT1_2)
3577             {   // 2(2x - 1)/(2x + 1)
3578                 exp -= 1;
3579                 z = x - 0.5;
3580                 y = 0.5 * z + 0.5;
3581             }
3582             else
3583             {   // 2(x - 1)/(x + 1)
3584                 z = x - 0.5;
3585                 z -= 0.5;
3586                 y = 0.5 * x  + 0.5;
3587             }
3588             x = z / y;
3589             z = x * x;
3590             y = x * (z * poly(z, coeffs.logR) / poly(z, coeffs.logS));
3591             goto Ldone;
3592         }
3593     }
3594 
3595     // Logarithm using log(1 + x) = x - .5x^^2 + x^^3 P(x) / Q(x)
3596     if (x < SQRT1_2)
3597     {
3598         exp -= 1;
3599         x = 2.0 * x - 1.0;
3600     }
3601     else
3602         x = x - 1.0;
3603 
3604     z = x * x;
3605     static if (F.realFormat == RealFormat.ieeeSingle)
3606         y = x * (z * poly(x, coeffs.log10P));
3607     else
3608         y = x * (z * poly(x, coeffs.log10P) / poly(x, coeffs.log10Q));
3609     y = y - 0.5 * z;
3610 
3611     // Multiply log of fraction by log10(e) and base 2 exponent by log10(2).
3612     // This sequence of operations is critical and it may be horribly
3613     // defeated by some compiler optimizers.
3614 Ldone:
3615     z = y * L10EB;
3616     z += x * L10EB;
3617     z += exp * L102B;
3618     z += y * L10EA;
3619     z += x * L10EA;
3620     z += exp * L102A;
3621 
3622     return z;
3623 }
3624 
3625 @safe @nogc nothrow unittest
3626 {
3627     import std.math : floatTraits, RealFormat;
3628     import std.meta : AliasSeq;
3629 
3630     static void testLog10(T)(T[2][] vals)
3631     {
3632         import std.math.operations : isClose;
3633         import std.math.traits : isNaN;
3634         foreach (ref pair; vals)
3635         {
3636             if (isNaN(pair[1]))
3637                 assert(isNaN(log10(pair[0])));
3638             else
3639                 assert(isClose(log10(pair[0]), pair[1]));
3640         }
3641     }
3642     static foreach (F; AliasSeq!(float, double, real))
3643     {{
3644         scope F[2][] vals = [
3645             [F(1), F(0x0p+0)], [F(2), F(0x1.34413509f79fef32p-2)],
3646             [F(4), F(0x1.34413509f79fef32p-1)], [F(8), F(0x1.ce61cf8ef36fe6cap-1)],
3647             [F(16), F(0x1.34413509f79fef32p+0)], [F(32), F(0x1.8151824c7587eafep+0)],
3648             [F(64), F(0x1.ce61cf8ef36fe6cap+0)], [F(128), F(0x1.0db90e68b8abf14cp+1)],
3649             [F(256), F(0x1.34413509f79fef32p+1)], [F(512), F(0x1.5ac95bab3693ed18p+1)],
3650             [F(1024), F(0x1.8151824c7587eafep+1)], [F(2048), F(0x1.a7d9a8edb47be8e4p+1)],
3651             [F(3), F(0x1.e8927964fd5fd08cp-2)], [F(5), F(0x1.65df657b04300868p-1)],
3652             [F(7), F(0x1.b0b0b0b78cc3f296p-1)], [F(15), F(0x1.2d145116c16ff856p+0)],
3653             [F(17), F(0x1.3afeb354b7d9731ap+0)], [F(31), F(0x1.7dc9e145867e62eap+0)],
3654             [F(33), F(0x1.84bd545e4baeddp+0)], [F(63), F(0x1.cca1950e4511e192p+0)],
3655             [F(65), F(0x1.d01b16f9433cf7b8p+0)], [F(-0), -F.infinity], [F(0), -F.infinity],
3656             [F(0.1), F(-0x1p+0)], [F(0.25), F(-0x1.34413509f79fef32p-1)],
3657             [F(0.75), F(-0x1.ffbfc2bbc7803758p-4)], [F(0.875), F(-0x1.db11ed766abf432ep-5)],
3658             [F(10), F(0x1p+0)], [F(100), F(0x1p+1)], [F(10000), F(0x1p+2)],
3659         ];
3660         testLog10(vals);
3661     }}
3662     {
3663         float[2][16] vals = [
3664             [float.nan, float.nan],[-float.nan, float.nan],
3665             [float.infinity, float.infinity], [-float.infinity, float.nan],
3666             [float.min_normal, -0x1.2f703p+5f], [-float.min_normal, float.nan],
3667             [float.max, 0x1.344136p+5f], [-float.max, float.nan],
3668             [float.min_normal / 2, -0x1.31d8b2p+5f], [-float.min_normal / 2, float.nan],
3669             [float.max / 2, 0x1.31d8b2p+5f], [-float.max / 2, float.nan],
3670             [float.min_normal / 3, -0x1.334156p+5f], [-float.min_normal / 3, float.nan],
3671             [float.max / 3, 0x1.30701p+5f], [-float.max / 3, float.nan],
3672         ];
3673         testLog10(vals);
3674     }
3675     {
3676         double[2][16] vals = [
3677             [double.nan, double.nan],[-double.nan, double.nan],
3678             [double.infinity, double.infinity], [-double.infinity, double.nan],
3679             [double.min_normal, -0x1.33a7146f72a42p+8], [-double.min_normal, double.nan],
3680             [double.max, 0x1.34413509f79ffp+8], [-double.max, double.nan],
3681             [double.min_normal / 2, -0x1.33f424bcb522p+8], [-double.min_normal / 2, double.nan],
3682             [double.max / 2, 0x1.33f424bcb522p+8], [-double.max / 2, double.nan],
3683             [double.min_normal / 3, -0x1.3421390dcbe37p+8], [-double.min_normal / 3, double.nan],
3684             [double.max / 3, 0x1.33c7106b9e609p+8], [-double.max / 3, double.nan],
3685         ];
3686         testLog10(vals);
3687     }
3688     alias F = floatTraits!real;
3689     static if (F.realFormat == RealFormat.ieeeExtended || F.realFormat == RealFormat.ieeeQuadruple)
3690     {{
3691         real[2][16] vals = [
3692             [real.nan, real.nan],[-real.nan, real.nan],
3693             [real.infinity, real.infinity], [-real.infinity, real.nan],
3694             [real.min_normal, -0x1.343793004f503232p+12L], [-real.min_normal, real.nan],
3695             [real.max, 0x1.34413509f79fef32p+12L], [-real.max, real.nan],
3696             [real.min_normal / 2, -0x1.343c6405237810b2p+12L], [-real.min_normal / 2, real.nan],
3697             [real.max / 2, 0x1.343c6405237810b2p+12L], [-real.max / 2, real.nan],
3698             [real.min_normal / 3, -0x1.343f354a34e427bp+12L], [-real.min_normal / 3, real.nan],
3699             [real.max / 3, 0x1.343992c0120bf9b2p+12L], [-real.max / 3, real.nan],
3700         ];
3701         testLog10(vals);
3702     }}
3703 }
3704 
3705 /**
3706  * Calculates the natural logarithm of 1 + x.
3707  *
3708  * For very small x, log1p(x) will be more accurate than
3709  * log(1 + x).
3710  *
3711  *  $(TABLE_SV
3712  *  $(TR $(TH x)            $(TH log1p(x))     $(TH divide by 0?) $(TH invalid?))
3713  *  $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no)           $(TD no))
3714  *  $(TR $(TD -1.0)         $(TD -$(INFIN))    $(TD yes)          $(TD no))
3715  *  $(TR $(TD $(LT)-1.0)    $(TD -$(NAN))      $(TD no)           $(TD yes))
3716  *  $(TR $(TD +$(INFIN))    $(TD +$(INFIN))    $(TD no)           $(TD no))
3717  *  )
3718  */
3719 pragma(inline, true)
3720 real log1p(real x) @safe pure nothrow @nogc
3721 {
3722     version (INLINE_YL2X)
3723     {
3724         // On x87, yl2xp1 is valid if and only if -0.5 <= lg(x) <= 0.5,
3725         //    ie if -0.29 <= x <= 0.414
3726         import std.math.constants : LN2;
3727         return (core.math.fabs(x) <= 0.25)  ? core.math.yl2xp1(x, LN2) : core.math.yl2x(x+1, LN2);
3728     }
3729     else
3730         return log1pImpl(x);
3731 }
3732 
3733 /// ditto
3734 pragma(inline, true)
3735 double log1p(double x) @safe pure nothrow @nogc { return __ctfe ? cast(double) log1p(cast(real) x) : log1pImpl(x); }
3736 
3737 /// ditto
3738 pragma(inline, true)
3739 float log1p(float x) @safe pure nothrow @nogc { return __ctfe ? cast(float) log1p(cast(real) x) : log1pImpl(x); }
3740 
3741 // @@@DEPRECATED_[2.112.0]@@@
3742 deprecated("`std.math.exponential.log1p` called with argument types `(int)` matches both "
3743            ~ "`log1p(real)`, `log1p(double)`, and `log1p(float)`. Cast argument to floating point type instead.")
3744 real log1p(int x) @safe pure nothrow @nogc { return log1p(cast(real) x); }
3745 // @@@DEPRECATED_[2.112.0]@@@
3746 deprecated("`std.math.exponential.log1p` called with argument types `(uint)` matches both "
3747            ~ "`log1p(real)`, `log1p(double)`, and `log1p(float)`. Cast argument to floating point type instead.")
3748 real log1p(uint x) @safe pure nothrow @nogc { return log1p(cast(real) x); }
3749 // @@@DEPRECATED_[2.112.0]@@@
3750 deprecated("`std.math.exponential.log1p` called with argument types `(long)` matches both "
3751            ~ "`log1p(real)`, `log1p(double)`, and `log1p(float)`. Cast argument to floating point type instead.")
3752 real log1p(long x) @safe pure nothrow @nogc { return log1p(cast(real) x); }
3753 // @@@DEPRECATED_[2.112.0]@@@
3754 deprecated("`std.math.exponential.log1p` called with argument types `(ulong)` matches both "
3755            ~ "`log1p(real)`, `log1p(double)`, and `log1p(float)`. Cast argument to floating point type instead.")
3756 real log1p(ulong x) @safe pure nothrow @nogc { return log1p(cast(real) x); }
3757 
3758 ///
3759 @safe pure unittest
3760 {
3761     import std.math.traits : isIdentical, isNaN;
3762     import std.math.operations : feqrel;
3763 
3764     assert(isIdentical(log1p(0.0), 0.0));
3765     assert(log1p(1.0).feqrel(0.69314) > 16);
3766 
3767     assert(log1p(-1.0) == -real.infinity);
3768     assert(isNaN(log1p(-2.0)));
3769     assert(log1p(real.nan) is real.nan);
3770     assert(log1p(-real.nan) is -real.nan);
3771     assert(log1p(real.infinity) == real.infinity);
3772 }
3773 
3774 private T log1pImpl(T)(T x) @safe pure nothrow @nogc
3775 {
3776     import std.math.traits : isNaN, isInfinity, signbit;
3777     import std.math.algebraic : poly;
3778     import std.math.constants : SQRT1_2, SQRT2;
3779     import std.math : floatTraits, RealFormat;
3780 
3781     // Special cases.
3782     if (isNaN(x) || x == 0.0)
3783         return x;
3784     if (isInfinity(x) && !signbit(x))
3785         return x;
3786     if (x == -1.0)
3787         return -T.infinity;
3788     if (x < -1.0)
3789         return T.nan;
3790 
3791     alias F = floatTraits!T;
3792     static if (F.realFormat == RealFormat.ieeeSingle ||
3793                F.realFormat == RealFormat.ieeeDouble)
3794     {
3795         // When the input is within the range 1/sqrt(2) <= x+1 <= sqrt(2), compute
3796         // log1p inline. Forwarding to log() would otherwise result in inaccuracies.
3797         const T xp1 = x + 1.0;
3798         if (xp1 >= SQRT1_2 && xp1 <= SQRT2)
3799         {
3800             alias coeffs = LogCoeffs!T;
3801 
3802             T px = poly(x, coeffs.logp1P);
3803             T qx = poly(x, coeffs.logp1Q);
3804             const T xx = x * x;
3805             qx = x + ((cast(T) -0.5) * xx + x * (xx * px / qx));
3806             return qx;
3807         }
3808     }
3809 
3810     return logImpl!(T, true)(x);
3811 }
3812 
3813 @safe @nogc nothrow unittest
3814 {
3815     import std.math : floatTraits, RealFormat;
3816     import std.meta : AliasSeq;
3817 
3818     static void testLog1p(T)(T[2][] vals)
3819     {
3820         import std.math.operations : isClose;
3821         import std.math.traits : isNaN;
3822         foreach (ref pair; vals)
3823         {
3824             if (isNaN(pair[1]))
3825                 assert(isNaN(log1p(pair[0])));
3826             else
3827                 assert(isClose(log1p(pair[0]), pair[1]));
3828         }
3829     }
3830     static foreach (F; AliasSeq!(float, double, real))
3831     {{
3832         scope F[2][] vals = [
3833             [F(1), F(0x1.62e42fefa39ef358p-1)], [F(2), F(0x1.193ea7aad030a976p+0)],
3834             [F(4), F(0x1.9c041f7ed8d336bp+0)], [F(8), F(0x1.193ea7aad030a976p+1)],
3835             [F(16), F(0x1.6aa6bc1fa7f79cfp+1)], [F(32), F(0x1.bf8d8f4d5b8d1038p+1)],
3836             [F(64), F(0x1.0b292939429755ap+2)], [F(128), F(0x1.37072a9b5b6cb31p+2)],
3837             [F(256), F(0x1.63241004e9010ad8p+2)], [F(512), F(0x1.8f60adf041bde2a8p+2)],
3838             [F(1024), F(0x1.bbad39ebe1cc08b6p+2)], [F(2048), F(0x1.e801c1698ba4395cp+2)],
3839             [F(3), F(0x1.62e42fefa39ef358p+0)], [F(5), F(0x1.cab0bfa2a2002322p+0)],
3840             [F(7), F(0x1.0a2b23f3bab73682p+1)], [F(15), F(0x1.62e42fefa39ef358p+1)],
3841             [F(17), F(0x1.71f7b3a6b918664cp+1)], [F(31), F(0x1.bb9d3beb8c86b02ep+1)],
3842             [F(33), F(0x1.c35fc81b90df59c6p+1)], [F(63), F(0x1.0a2b23f3bab73682p+2)],
3843             [F(65), F(0x1.0c234da4a23a6686p+2)], [F(-0), F(-0x0p+0)], [F(0), F(0x0p+0)],
3844             [F(0.1), F(0x1.8663f793c46c69cp-4)], [F(0.25), F(0x1.c8ff7c79a9a21ac4p-3)],
3845             [F(0.75), F(0x1.1e85f5e7040d03ep-1)], [F(0.875), F(0x1.41d8fe84672ae646p-1)],
3846             [F(10), F(0x1.32ee3b77f374bb7cp+1)], [F(100), F(0x1.275e2271bba30be4p+2)],
3847             [F(10000), F(0x1.26bbed6fbd84182ep+3)],
3848         ];
3849         testLog1p(vals);
3850     }}
3851     {
3852         float[2][16] vals = [
3853             [float.nan, float.nan],[-float.nan, float.nan],
3854             [float.infinity, float.infinity], [-float.infinity, float.nan],
3855             [float.min_normal, 0x1p-126f], [-float.min_normal, -0x1p-126f],
3856             [float.max, 0x1.62e43p+6f], [-float.max, float.nan],
3857             [float.min_normal / 2, 0x0.8p-126f], [-float.min_normal / 2, -0x0.8p-126f],
3858             [float.max / 2, 0x1.601e68p+6f], [-float.max / 2, float.nan],
3859             [float.min_normal / 3, 0x0.555556p-126f], [-float.min_normal / 3, -0x0.555556p-126f],
3860             [float.max / 3, 0x1.5e7f36p+6f], [-float.max / 3, float.nan],
3861         ];
3862         testLog1p(vals);
3863     }
3864     {
3865         double[2][16] vals = [
3866             [double.nan, double.nan],[-double.nan, double.nan],
3867             [double.infinity, double.infinity], [-double.infinity, double.nan],
3868             [double.min_normal, 0x1p-1022], [-double.min_normal, -0x1p-1022],
3869             [double.max, 0x1.62e42fefa39efp+9], [-double.max, double.nan],
3870             [double.min_normal / 2, 0x0.8p-1022], [-double.min_normal / 2, -0x0.8p-1022],
3871             [double.max / 2, 0x1.628b76e3a7b61p+9], [-double.max / 2, double.nan],
3872             [double.min_normal / 3, 0x0.5555555555555p-1022], [-double.min_normal / 3, -0x0.5555555555555p-1022],
3873             [double.max / 3, 0x1.6257909bce36ep+9], [-double.max / 3, double.nan],
3874         ];
3875         testLog1p(vals);
3876     }
3877     alias F = floatTraits!real;
3878     static if (F.realFormat == RealFormat.ieeeExtended || F.realFormat == RealFormat.ieeeQuadruple)
3879     {{
3880         real[2][16] vals = [
3881             [real.nan, real.nan],[-real.nan, real.nan],
3882             [real.infinity, real.infinity], [-real.infinity, real.nan],
3883             [real.min_normal, 0x1p-16382L], [-real.min_normal, -0x1p-16382L],
3884             [real.max, 0x1.62e42fefa39ef358p+13L], [-real.max, real.nan],
3885             [real.min_normal / 2, 0x0.8p-16382L], [-real.min_normal / 2, -0x0.8p-16382L],
3886             [real.max / 2, 0x1.62dea45ee3e064dcp+13L], [-real.max / 2, real.nan],
3887             [real.min_normal / 3, 0x0.5555555555555556p-16382L], [-real.min_normal / 3, -0x0.5555555555555556p-16382L],
3888             [real.max / 3, 0x1.62db65fa664871d2p+13L], [-real.max / 3, real.nan],
3889         ];
3890         testLog1p(vals);
3891     }}
3892 }
3893 
3894 /***************************************
3895  * Calculates the base-2 logarithm of x:
3896  * $(SUB log, 2)x
3897  *
3898  *  $(TABLE_SV
3899  *  $(TR $(TH x)            $(TH log2(x))   $(TH divide by 0?) $(TH invalid?))
3900  *  $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes)          $(TD no) )
3901  *  $(TR $(TD $(LT)0.0)     $(TD $(NAN))    $(TD no)           $(TD yes) )
3902  *  $(TR $(TD +$(INFIN))    $(TD +$(INFIN)) $(TD no)           $(TD no) )
3903  *  )
3904  */
3905 pragma(inline, true)
3906 real log2(real x) @safe pure nothrow @nogc
3907 {
3908     version (INLINE_YL2X)
3909         return core.math.yl2x(x, 1.0L);
3910     else
3911         return log2Impl(x);
3912 }
3913 
3914 /// ditto
3915 pragma(inline, true)
3916 double log2(double x) @safe pure nothrow @nogc { return __ctfe ? cast(double) log2(cast(real) x) : log2Impl(x); }
3917 
3918 /// ditto
3919 pragma(inline, true)
3920 float log2(float x) @safe pure nothrow @nogc { return __ctfe ? cast(float) log2(cast(real) x) : log2Impl(x); }
3921 
3922 // @@@DEPRECATED_[2.112.0]@@@
3923 deprecated("`std.math.exponential.log2` called with argument types `(int)` matches both "
3924            ~ "`log2(real)`, `log2(double)`, and `log2(float)`. Cast argument to floating point type instead.")
3925 real log2(int x) @safe pure nothrow @nogc { return log2(cast(real) x); }
3926 // @@@DEPRECATED_[2.112.0]@@@
3927 deprecated("`std.math.exponential.log2` called with argument types `(uint)` matches both "
3928            ~ "`log2(real)`, `log2(double)`, and `log2(float)`. Cast argument to floating point type instead.")
3929 real log2(uint x) @safe pure nothrow @nogc { return log2(cast(real) x); }
3930 // @@@DEPRECATED_[2.112.0]@@@
3931 deprecated("`std.math.exponential.log2` called with argument types `(long)` matches both "
3932            ~ "`log2(real)`, `log2(double)`, and `log2(float)`. Cast argument to floating point type instead.")
3933 real log2(long x) @safe pure nothrow @nogc { return log2(cast(real) x); }
3934 // @@@DEPRECATED_[2.112.0]@@@
3935 deprecated("`std.math.exponential.log2` called with argument types `(ulong)` matches both "
3936            ~ "`log2(real)`, `log2(double)`, and `log2(float)`. Cast argument to floating point type instead.")
3937 real log2(ulong x) @safe pure nothrow @nogc { return log2(cast(real) x); }
3938 
3939 ///
3940 @safe unittest
3941 {
3942     import std.math.operations : isClose;
3943 
3944     assert(isClose(log2(1024.0L), 10));
3945 }
3946 
3947 @safe @nogc nothrow unittest
3948 {
3949     import std.math.operations : isClose;
3950 
3951     // check if values are equal to 19 decimal digits of precision
3952     assert(isClose(log2(1024.0L), 10, 1e-18));
3953 }
3954 
3955 private T log2Impl(T)(T x) @safe pure nothrow @nogc
3956 {
3957     import std.math.traits : isNaN, isInfinity, signbit;
3958     import std.math.constants : SQRT1_2, LOG2E;
3959     import std.math.algebraic : poly;
3960     import std.math : floatTraits, RealFormat;
3961 
3962     alias coeffs = LogCoeffs!T;
3963     alias F = floatTraits!T;
3964 
3965     // Special cases are the same as for log.
3966     if (isNaN(x))
3967         return x;
3968     if (isInfinity(x) && !signbit(x))
3969         return x;
3970     if (x == 0.0)
3971         return -T.infinity;
3972     if (x < 0.0)
3973         return T.nan;
3974 
3975     // Separate mantissa from exponent.
3976     // Note, frexp is used so that denormal numbers will be handled properly.
3977     T y, z;
3978     int exp;
3979 
3980     x = frexp(x, exp);
3981 
3982     static if (F.realFormat == RealFormat.ieeeDouble ||
3983                F.realFormat == RealFormat.ieeeExtended ||
3984                F.realFormat == RealFormat.ieeeExtended53 ||
3985                F.realFormat == RealFormat.ieeeQuadruple)
3986     {
3987         // Logarithm using log(x) = z + z^^3 R(z) / S(z),
3988         // where z = 2(x - 1)/(x + 1)
3989         if ((exp > 2) || (exp < -2))
3990         {
3991             if (x < SQRT1_2)
3992             {   // 2(2x - 1)/(2x + 1)
3993                 exp -= 1;
3994                 z = x - 0.5;
3995                 y = 0.5 * z + 0.5;
3996             }
3997             else
3998             {   // 2(x - 1)/(x + 1)
3999                 z = x - 0.5;
4000                 z -= 0.5;
4001                 y = 0.5 * x  + 0.5;
4002             }
4003             x = z / y;
4004             z = x * x;
4005             y = x * (z * poly(z, coeffs.logR) / poly(z, coeffs.logS));
4006             goto Ldone;
4007         }
4008     }
4009 
4010     // Logarithm using log(1 + x) = x - .5x^^2 + x^^3 P(x) / Q(x)
4011     if (x < SQRT1_2)
4012     {
4013         exp -= 1;
4014         x = 2.0 * x - 1.0;
4015     }
4016     else
4017         x = x - 1.0;
4018 
4019     z = x * x;
4020     static if (F.realFormat == RealFormat.ieeeSingle)
4021         y = x * (z * poly(x, coeffs.log2P));
4022     else
4023         y = x * (z * poly(x, coeffs.log2P) / poly(x, coeffs.log2Q));
4024     y = y - 0.5 * z;
4025 
4026     // Multiply log of fraction by log10(e) and base 2 exponent by log10(2).
4027     // This sequence of operations is critical and it may be horribly
4028     // defeated by some compiler optimizers.
4029 Ldone:
4030     z = y * (LOG2E - 1.0);
4031     z += x * (LOG2E - 1.0);
4032     z += y;
4033     z += x;
4034     z += exp;
4035 
4036     return z;
4037 }
4038 
4039 @safe @nogc nothrow unittest
4040 {
4041     import std.math : floatTraits, RealFormat;
4042     import std.meta : AliasSeq;
4043 
4044     static void testLog2(T)(T[2][] vals)
4045     {
4046         import std.math.operations : isClose;
4047         import std.math.traits : isNaN;
4048         foreach (ref pair; vals)
4049         {
4050             if (isNaN(pair[1]))
4051                 assert(isNaN(log2(pair[0])));
4052             else
4053                 assert(isClose(log2(pair[0]), pair[1]));
4054         }
4055     }
4056     static foreach (F; AliasSeq!(float, double, real))
4057     {{
4058         scope F[2][] vals = [
4059             [F(1), F(0x0p+0)], [F(2), F(0x1p+0)],
4060             [F(4), F(0x1p+1)], [F(8), F(0x1.8p+1)],
4061             [F(16), F(0x1p+2)], [F(32), F(0x1.4p+2)],
4062             [F(64), F(0x1.8p+2)], [F(128), F(0x1.cp+2)],
4063             [F(256), F(0x1p+3)], [F(512), F(0x1.2p+3)],
4064             [F(1024), F(0x1.4p+3)], [F(2048), F(0x1.6p+3)],
4065             [F(3), F(0x1.95c01a39fbd687ap+0)], [F(5), F(0x1.2934f0979a3715fcp+1)],
4066             [F(7), F(0x1.675767f54042cd9ap+1)], [F(15), F(0x1.f414fdb4982259ccp+1)],
4067             [F(17), F(0x1.0598fdbeb244c5ap+2)], [F(31), F(0x1.3d118d66c4d4e554p+2)],
4068             [F(33), F(0x1.42d75a6eb1dfb0e6p+2)], [F(63), F(0x1.7e8bc1179e0caa9cp+2)],
4069             [F(65), F(0x1.816e79685c2d2298p+2)], [F(-0), -F.infinity], [F(0), -F.infinity],
4070             [F(0.1), F(-0x1.a934f0979a3715fcp+1)], [F(0.25), F(-0x1p+1)],
4071             [F(0.75), F(-0x1.a8ff971810a5e182p-2)], [F(0.875), F(-0x1.8a8980abfbd32668p-3)],
4072             [F(10), F(0x1.a934f0979a3715fcp+1)], [F(100), F(0x1.a934f0979a3715fcp+2)],
4073             [F(10000), F(0x1.a934f0979a3715fcp+3)],
4074         ];
4075         testLog2(vals);
4076     }}
4077     {
4078         float[2][16] vals = [
4079             [float.nan, float.nan],[-float.nan, float.nan],
4080             [float.infinity, float.infinity], [-float.infinity, float.nan],
4081             [float.min_normal, -0x1.f8p+6f], [-float.min_normal, float.nan],
4082             [float.max, 0x1p+7f], [-float.max, float.nan],
4083             [float.min_normal / 2, -0x1.fcp+6f], [-float.min_normal / 2, float.nan],
4084             [float.max / 2, 0x1.fcp+6f], [-float.max / 2, float.nan],
4085             [float.min_normal / 3, -0x1.fe57p+6f], [-float.min_normal / 3, float.nan],
4086             [float.max / 3, 0x1.f9a9p+6f], [-float.max / 3, float.nan],
4087         ];
4088         testLog2(vals);
4089     }
4090     {
4091         double[2][16] vals = [
4092             [double.nan, double.nan],[-double.nan, double.nan],
4093             [double.infinity, double.infinity], [-double.infinity, double.nan],
4094             [double.min_normal, -0x1.ffp+9], [-double.min_normal, double.nan],
4095             [double.max, 0x1p+10], [-double.max, double.nan],
4096             [double.min_normal / 2, -0x1.ff8p+9], [-double.min_normal / 2, double.nan],
4097             [double.max / 2, 0x1.ff8p+9], [-double.max / 2, double.nan],
4098             [double.min_normal / 3, -0x1.ffcae00d1cfdfp+9], [-double.min_normal / 3, double.nan],
4099             [double.max / 3, 0x1.ff351ff2e3021p+9], [-double.max / 3, double.nan],
4100         ];
4101         testLog2(vals);
4102     }
4103     alias F = floatTraits!real;
4104     static if (F.realFormat == RealFormat.ieeeExtended || F.realFormat == RealFormat.ieeeQuadruple)
4105     {{
4106         real[2][16] vals = [
4107             [real.nan, real.nan],[-real.nan, real.nan],
4108             [real.infinity, real.infinity], [-real.infinity, real.nan],
4109             [real.min_normal, -0x1.fffp+13L], [-real.min_normal, real.nan],
4110             [real.max, 0x1p+14L], [-real.max, real.nan],
4111             [real.min_normal / 2, -0x1.fff8p+13L], [-real.min_normal / 2, real.nan],
4112             [real.max / 2, 0x1.fff8p+13L], [-real.max / 2, real.nan],
4113             [real.min_normal / 3, -0x1.fffcae00d1cfdeb4p+13L], [-real.min_normal / 3, real.nan],
4114             [real.max / 3, 0x1.fff351ff2e30214cp+13L], [-real.max / 3, real.nan],
4115         ];
4116         testLog2(vals);
4117     }}
4118 }
4119 
4120 /*****************************************
4121  * Extracts the exponent of x as a signed integral value.
4122  *
4123  * If x is subnormal, it is treated as if it were normalized.
4124  * For a positive, finite x:
4125  *
4126  * 1 $(LT)= $(I x) * FLT_RADIX$(SUPERSCRIPT -logb(x)) $(LT) FLT_RADIX
4127  *
4128  *      $(TABLE_SV
4129  *      $(TR $(TH x)                 $(TH logb(x))   $(TH divide by 0?) )
4130  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD +$(INFIN)) $(TD no))
4131  *      $(TR $(TD $(PLUSMN)0.0)      $(TD -$(INFIN)) $(TD yes) )
4132  *      )
4133  */
4134 pragma(inline, true)
4135 real logb(real x) @trusted pure nothrow @nogc
4136 {
4137     version (InlineAsm_X87_MSVC)
4138         return logbAsm(x);
4139     else
4140         return logbImpl(x);
4141 }
4142 
4143 /// ditto
4144 pragma(inline, true)
4145 double logb(double x) @trusted pure nothrow @nogc { return logbImpl(x); }
4146 
4147 /// ditto
4148 pragma(inline, true)
4149 float logb(float x) @trusted pure nothrow @nogc { return logbImpl(x); }
4150 
4151 ///
4152 @safe @nogc nothrow unittest
4153 {
4154     assert(logb(1.0) == 0);
4155     assert(logb(100.0) == 6);
4156 
4157     assert(logb(0.0) == -real.infinity);
4158     assert(logb(real.infinity) == real.infinity);
4159     assert(logb(-real.infinity) == real.infinity);
4160 }
4161 
4162 @safe @nogc nothrow unittest
4163 {
4164     import std.meta : AliasSeq;
4165     import std.typecons : Tuple;
4166     import std.math.traits : isNaN;
4167     static foreach (F; AliasSeq!(float, double, real))
4168     {{
4169         alias T = Tuple!(F, F);
4170         T[17] vals =   // x, logb(x)
4171         [
4172             T(1.0          , 0          ),
4173             T(100.0        , 6          ),
4174             T(0.0          , -F.infinity),
4175             T(-0.0         , -F.infinity),
4176             T(1024         , 10         ),
4177             T(-2000        , 10         ),
4178             T(0x0.1p-127   , -131       ),
4179             T(0x0.01p-127  , -135       ),
4180             T(0x0.011p-127 , -135       ),
4181             T(F.nan        , F.nan      ),
4182             T(-F.nan       , F.nan      ),
4183             T(F.infinity   , F.infinity ),
4184             T(-F.infinity  , F.infinity ),
4185             T(F.min_normal , F.min_exp-1),
4186             T(-F.min_normal, F.min_exp-1),
4187             T(F.max        , F.max_exp-1),
4188             T(-F.max       , F.max_exp-1),
4189         ];
4190 
4191         foreach (elem; vals)
4192         {
4193             if (isNaN(elem[1]))
4194                 assert(isNaN(logb(elem[1])));
4195             else
4196                 assert(logb(elem[0]) == elem[1]);
4197         }
4198     }}
4199 }
4200 
4201 version (InlineAsm_X87_MSVC)
4202 private T logbAsm(T)(T x) @trusted pure nothrow @nogc
4203 {
4204     version (X86_64)
4205     {
4206         asm pure nothrow @nogc
4207         {
4208             naked                       ;
4209             fld     real ptr [RCX]      ;
4210             fxtract                     ;
4211             fstp    ST(0)               ;
4212             ret                         ;
4213         }
4214     }
4215     else
4216     {
4217         asm pure nothrow @nogc
4218         {
4219             fld     x                   ;
4220             fxtract                     ;
4221             fstp    ST(0)               ;
4222         }
4223     }
4224 }
4225 
4226 private T logbImpl(T)(T x) @trusted pure nothrow @nogc
4227 {
4228     import std.math.traits : isFinite;
4229 
4230     // Handle special cases.
4231     if (!isFinite(x))
4232         return x * x;
4233     if (x == 0)
4234         return -1 / (x * x);
4235 
4236     return ilogb(x);
4237 }
4238 
4239 @safe @nogc nothrow unittest
4240 {
4241     import std.math : floatTraits, RealFormat;
4242     import std.meta : AliasSeq;
4243 
4244     static void testLogb(T)(T[2][] vals)
4245     {
4246         import std.math.operations : isClose;
4247         import std.math.traits : isNaN;
4248         foreach (ref pair; vals)
4249         {
4250             if (isNaN(pair[1]))
4251                 assert(isNaN(logb(pair[0])));
4252             else
4253                 assert(isClose(logb(pair[0]), pair[1]));
4254         }
4255     }
4256     static foreach (F; AliasSeq!(float, double, real))
4257     {{
4258         scope F[2][] vals = [
4259             [F(1), F(0x0p+0)], [F(2), F(0x1p+0)],
4260             [F(4), F(0x1p+1)], [F(8), F(0x1.8p+1)],
4261             [F(16), F(0x1p+2)], [F(32), F(0x1.4p+2)],
4262             [F(64), F(0x1.8p+2)], [F(128), F(0x1.cp+2)],
4263             [F(256), F(0x1p+3)], [F(512), F(0x1.2p+3)],
4264             [F(1024), F(0x1.4p+3)], [F(2048), F(0x1.6p+3)],
4265             [F(3), F(0x1p+0)], [F(5), F(0x1p+1)],
4266             [F(7), F(0x1p+1)], [F(15), F(0x1.8p+1)],
4267             [F(17), F(0x1p+2)], [F(31), F(0x1p+2)],
4268             [F(33), F(0x1.4p+2)], [F(63), F(0x1.4p+2)],
4269             [F(65), F(0x1.8p+2)], [F(-0), -F.infinity], [F(0), -F.infinity],
4270             [F(0.1), F(-0x1p+2)], [F(0.25), F(-0x1p+1)],
4271             [F(0.75), F(-0x1p+0)], [F(0.875), F(-0x1p+0)],
4272             [F(10), F(0x1.8p+1)], [F(100), F(0x1.8p+2)],
4273             [F(10000), F(0x1.ap+3)],
4274         ];
4275         testLogb(vals);
4276     }}
4277     {
4278         float[2][16] vals = [
4279             [float.nan, float.nan],[-float.nan, float.nan],
4280             [float.infinity, float.infinity], [-float.infinity, float.infinity],
4281             [float.min_normal, -0x1.f8p+6f], [-float.min_normal, -0x1.f8p+6f],
4282             [float.max, 0x1.fcp+6f], [-float.max, 0x1.fcp+6f],
4283             [float.min_normal / 2, -0x1.fcp+6f], [-float.min_normal / 2, -0x1.fcp+6f],
4284             [float.max / 2, 0x1.f8p+6f], [-float.max / 2, 0x1.f8p+6f],
4285             [float.min_normal / 3, -0x1p+7f], [-float.min_normal / 3, -0x1p+7f],
4286             [float.max / 3, 0x1.f8p+6f], [-float.max / 3, 0x1.f8p+6f],
4287         ];
4288         testLogb(vals);
4289     }
4290     {
4291         double[2][16] vals = [
4292             [double.nan, double.nan],[-double.nan, double.nan],
4293             [double.infinity, double.infinity], [-double.infinity, double.infinity],
4294             [double.min_normal, -0x1.ffp+9], [-double.min_normal, -0x1.ffp+9],
4295             [double.max, 0x1.ff8p+9], [-double.max, 0x1.ff8p+9],
4296             [double.min_normal / 2, -0x1.ff8p+9], [-double.min_normal / 2, -0x1.ff8p+9],
4297             [double.max / 2, 0x1.ffp+9], [-double.max / 2, 0x1.ffp+9],
4298             [double.min_normal / 3, -0x1p+10], [-double.min_normal / 3, -0x1p+10],
4299             [double.max / 3, 0x1.ffp+9], [-double.max / 3, 0x1.ffp+9],
4300         ];
4301         testLogb(vals);
4302     }
4303     alias F = floatTraits!real;
4304     static if (F.realFormat == RealFormat.ieeeExtended || F.realFormat == RealFormat.ieeeQuadruple)
4305     {{
4306         real[2][16] vals = [
4307             [real.nan, real.nan],[-real.nan, real.nan],
4308             [real.infinity, real.infinity], [-real.infinity, real.infinity],
4309             [real.min_normal, -0x1.fffp+13L], [-real.min_normal, -0x1.fffp+13L],
4310             [real.max, 0x1.fff8p+13L], [-real.max, 0x1.fff8p+13L],
4311             [real.min_normal / 2, -0x1.fff8p+13L], [-real.min_normal / 2, -0x1.fff8p+13L],
4312             [real.max / 2, 0x1.fffp+13L], [-real.max / 2, 0x1.fffp+13L],
4313             [real.min_normal / 3, -0x1p+14L], [-real.min_normal / 3, -0x1p+14L],
4314             [real.max / 3, 0x1.fffp+13L], [-real.max / 3, 0x1.fffp+13L],
4315         ];
4316         testLogb(vals);
4317     }}
4318 }
4319 
4320 /*************************************
4321  * Efficiently calculates x * 2$(SUPERSCRIPT n).
4322  *
4323  * scalbn handles underflow and overflow in
4324  * the same fashion as the basic arithmetic operators.
4325  *
4326  *      $(TABLE_SV
4327  *      $(TR $(TH x)                 $(TH scalb(x)))
4328  *      $(TR $(TD $(PLUSMNINF))      $(TD $(PLUSMNINF)) )
4329  *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(PLUSMN)0.0) )
4330  *      )
4331  */
4332 pragma(inline, true)
4333 real scalbn(real x, int n) @safe pure nothrow @nogc { return _scalbn(x,n); }
4334 
4335 /// ditto
4336 pragma(inline, true)
4337 double scalbn(double x, int n) @safe pure nothrow @nogc { return _scalbn(x,n); }
4338 
4339 /// ditto
4340 pragma(inline, true)
4341 float scalbn(float x, int n) @safe pure nothrow @nogc { return _scalbn(x,n); }
4342 
4343 ///
4344 @safe pure nothrow @nogc unittest
4345 {
4346     assert(scalbn(0x1.2345678abcdefp0L, 999) == 0x1.2345678abcdefp999L);
4347     assert(scalbn(-real.infinity, 5) == -real.infinity);
4348     assert(scalbn(2.0,10) == 2048.0);
4349     assert(scalbn(2048.0f,-10) == 2.0f);
4350 }
4351 
4352 pragma(inline, true)
4353 private F _scalbn(F)(F x, int n)
4354 {
4355     import std.math.traits : isInfinity;
4356 
4357     if (__ctfe)
4358     {
4359         // Handle special cases.
4360         if (x == F(0.0) || isInfinity(x))
4361             return x;
4362     }
4363     return core.math.ldexp(x, n);
4364 }
4365 
4366 @safe pure nothrow @nogc unittest
4367 {
4368     // CTFE-able test
4369     static assert(scalbn(0x1.2345678abcdefp0L, 999) == 0x1.2345678abcdefp999L);
4370     static assert(scalbn(-real.infinity, 5) == -real.infinity);
4371     // Test with large exponent delta n where the result is in bounds but 2.0L ^^ n is not.
4372     enum initialExponent = real.min_exp + 2, resultExponent = real.max_exp - 2;
4373     enum int n = resultExponent - initialExponent;
4374     enum real x = 0x1.2345678abcdefp0L * (2.0L ^^ initialExponent);
4375     enum staticResult = scalbn(x, n);
4376     static assert(staticResult == 0x1.2345678abcdefp0L * (2.0L ^^ resultExponent));
4377     assert(scalbn(x, n) == staticResult);
4378 }
4379
The OpenD Programming Language
