"Fossies" - the Fresh Open Source Software Archive

Member "dmd2/src/phobos/std/math.d" (20 Nov 2020, 303605 Bytes) of package /linux/misc/dmd.2.094.2.linux.tar.xz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) D source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file.

    1 // Written in the D programming language.
    2 
    3 /**
    4  * Contains the elementary mathematical functions (powers, roots,
    5  * and trigonometric functions), and low-level floating-point operations.
    6  * Mathematical special functions are available in `std.mathspecial`.
    7  *
    8 $(SCRIPT inhibitQuickIndex = 1;)
    9 
   10 $(DIVC quickindex,
   11 $(BOOKTABLE ,
   12 $(TR $(TH Category) $(TH Members) )
   13 $(TR $(TDNW Constants) $(TD
   14     $(MYREF E) $(MYREF PI) $(MYREF PI_2) $(MYREF PI_4) $(MYREF M_1_PI)
   15     $(MYREF M_2_PI) $(MYREF M_2_SQRTPI) $(MYREF LN10) $(MYREF LN2)
   16     $(MYREF LOG2) $(MYREF LOG2E) $(MYREF LOG2T) $(MYREF LOG10E)
   17     $(MYREF SQRT2) $(MYREF SQRT1_2)
   18 ))
   19 $(TR $(TDNW Classics) $(TD
   20     $(MYREF abs) $(MYREF fabs) $(MYREF sqrt) $(MYREF cbrt) $(MYREF hypot)
   21     $(MYREF poly) $(MYREF nextPow2) $(MYREF truncPow2)
   22 ))
   23 $(TR $(TDNW Trigonometry) $(TD
   24     $(MYREF sin) $(MYREF cos) $(MYREF tan) $(MYREF asin) $(MYREF acos)
   25     $(MYREF atan) $(MYREF atan2) $(MYREF sinh) $(MYREF cosh) $(MYREF tanh)
   26     $(MYREF asinh) $(MYREF acosh) $(MYREF atanh)
   27 ))
   28 $(TR $(TDNW Rounding) $(TD
   29     $(MYREF ceil) $(MYREF floor) $(MYREF round) $(MYREF lround)
   30     $(MYREF trunc) $(MYREF rint) $(MYREF lrint) $(MYREF nearbyint)
   31     $(MYREF rndtol) $(MYREF quantize)
   32 ))
   33 $(TR $(TDNW Exponentiation & Logarithms) $(TD
   34     $(MYREF pow) $(MYREF exp) $(MYREF exp2) $(MYREF expm1) $(MYREF ldexp)
   35     $(MYREF frexp) $(MYREF log) $(MYREF log2) $(MYREF log10) $(MYREF logb)
   36     $(MYREF ilogb) $(MYREF log1p) $(MYREF scalbn)
   37 ))
   38 $(TR $(TDNW Modulus) $(TD
   39     $(MYREF fmod) $(MYREF modf) $(MYREF remainder)
   40 ))
   41 $(TR $(TDNW Floating-point operations) $(TD
   42     $(MYREF approxEqual) $(MYREF feqrel) $(MYREF fdim) $(MYREF fmax)
   43     $(MYREF fmin) $(MYREF fma) $(MYREF isClose) $(MYREF nextDown) $(MYREF nextUp)
   44     $(MYREF nextafter) $(MYREF NaN) $(MYREF getNaNPayload)
   45     $(MYREF cmp)
   46 ))
   47 $(TR $(TDNW Introspection) $(TD
   48     $(MYREF isFinite) $(MYREF isIdentical) $(MYREF isInfinity) $(MYREF isNaN)
   49     $(MYREF isNormal) $(MYREF isSubnormal) $(MYREF signbit) $(MYREF sgn)
   50     $(MYREF copysign) $(MYREF isPowerOf2)
   51 ))
   52 $(TR $(TDNW Hardware Control) $(TD
   53     $(MYREF IeeeFlags) $(MYREF FloatingPointControl)
   54 ))
   55 )
   56 )
   57 
   58  * The functionality closely follows the IEEE754-2008 standard for
   59  * floating-point arithmetic, including the use of camelCase names rather
   60  * than C99-style lower case names. All of these functions behave correctly
   61  * when presented with an infinity or NaN.
   62  *
   63  * The following IEEE 'real' formats are currently supported:
   64  * $(UL
   65  * $(LI 64 bit Big-endian  'double' (eg PowerPC))
   66  * $(LI 128 bit Big-endian 'quadruple' (eg SPARC))
   67  * $(LI 64 bit Little-endian 'double' (eg x86-SSE2))
   68  * $(LI 80 bit Little-endian, with implied bit 'real80' (eg x87, Itanium))
   69  * $(LI 128 bit Little-endian 'quadruple' (not implemented on any known processor!))
   70  * $(LI Non-IEEE 128 bit Big-endian 'doubledouble' (eg PowerPC) has partial support)
   71  * )
   72  * Unlike C, there is no global 'errno' variable. Consequently, almost all of
   73  * these functions are pure nothrow.
   74  *
   75  * Macros:
   76  *      TABLE_SV = <table border="1" cellpadding="4" cellspacing="0">
   77  *              <caption>Special Values</caption>
   78  *              $0</table>
   79  *      SVH = $(TR $(TH $1) $(TH $2))
   80  *      SV  = $(TR $(TD $1) $(TD $2))
   81  *      TH3 = $(TR $(TH $1) $(TH $2) $(TH $3))
   82  *      TD3 = $(TR $(TD $1) $(TD $2) $(TD $3))
   83  *      TABLE_DOMRG = <table border="1" cellpadding="4" cellspacing="0">
   84  *              $(SVH Domain X, Range Y)
   85                 $(SV $1, $2)
   86  *              </table>
   87  *      DOMAIN=$1
   88  *      RANGE=$1
   89 
   90  *      NAN = $(RED NAN)
   91  *      SUP = <span style="vertical-align:super;font-size:smaller">$0</span>
   92  *      GAMMA = &#915;
   93  *      THETA = &theta;
   94  *      INTEGRAL = &#8747;
   95  *      INTEGRATE = $(BIG &#8747;<sub>$(SMALL $1)</sub><sup>$2</sup>)
   96  *      POWER = $1<sup>$2</sup>
   97  *      SUB = $1<sub>$2</sub>
   98  *      BIGSUM = $(BIG &Sigma; <sup>$2</sup><sub>$(SMALL $1)</sub>)
   99  *      CHOOSE = $(BIG &#40;) <sup>$(SMALL $1)</sup><sub>$(SMALL $2)</sub> $(BIG &#41;)
  100  *      PLUSMN = &plusmn;
  101  *      INFIN = &infin;
  102  *      PLUSMNINF = &plusmn;&infin;
  103  *      PI = &pi;
  104  *      LT = &lt;
  105  *      GT = &gt;
  106  *      SQRT = &radic;
  107  *      HALF = &frac12;
  108  *
  109  * Copyright: Copyright The D Language Foundation 2000 - 2011.
  110  *            D implementations of tan, atan, atan2, exp, expm1, exp2, log, log10, log1p,
  111  *            log2, floor, ceil and lrint functions are based on the CEPHES math library,
  112  *            which is Copyright (C) 2001 Stephen L. Moshier $(LT)steve@moshier.net$(GT)
  113  *            and are incorporated herein by permission of the author.  The author
  114  *            reserves the right to distribute this material elsewhere under different
  115  *            copying permissions.  These modifications are distributed here under
  116  *            the following terms:
  117  * License:   $(HTTP www.boost.org/LICENSE_1_0.txt, Boost License 1.0).
  118  * Authors:   $(HTTP digitalmars.com, Walter Bright), Don Clugston,
  119  *            Conversion of CEPHES math library to D by Iain Buclaw and David Nadlinger
  120  * Source: $(PHOBOSSRC std/math.d)
  121  */
  122 module std.math;
  123 
  124 version (Win64)
  125 {
  126     version (D_InlineAsm_X86_64)
  127         version = Win64_DMD_InlineAsm;
  128 }
  129 
  130 static import core.math;
  131 static import core.stdc.math;
  132 static import core.stdc.fenv;
  133 import std.traits :  CommonType, isFloatingPoint, isIntegral, isNumeric,
  134     isSigned, isUnsigned, Largest, Unqual;
  135 
  136 // @@@DEPRECATED_2.102@@@
  137 // Note: Exposed accidentally, should be deprecated / removed
  138 deprecated("std.meta.AliasSeq was unintentionally available from std.math "
  139            ~ "and will be removed after 2.102. Please import std.meta instead")
  140 public import std.meta : AliasSeq;
  141 
  142 version (DigitalMars)
  143 {
  144     version = INLINE_YL2X;        // x87 has opcodes for these
  145 }
  146 
  147 version (X86)       version = X86_Any;
  148 version (X86_64)    version = X86_Any;
  149 version (PPC)       version = PPC_Any;
  150 version (PPC64)     version = PPC_Any;
  151 version (MIPS32)    version = MIPS_Any;
  152 version (MIPS64)    version = MIPS_Any;
  153 version (AArch64)   version = ARM_Any;
  154 version (ARM)       version = ARM_Any;
  155 version (S390)      version = IBMZ_Any;
  156 version (SPARC)     version = SPARC_Any;
  157 version (SPARC64)   version = SPARC_Any;
  158 version (SystemZ)   version = IBMZ_Any;
  159 version (RISCV32)   version = RISCV_Any;
  160 version (RISCV64)   version = RISCV_Any;
  161 
  162 version (D_InlineAsm_X86)
  163 {
  164     version = InlineAsm_X86_Any;
  165 }
  166 else version (D_InlineAsm_X86_64)
  167 {
  168     version = InlineAsm_X86_Any;
  169 }
  170 
  171 version (CRuntime_Microsoft)
  172 {
  173     version (InlineAsm_X86_Any)
  174         version = MSVC_InlineAsm;
  175 }
  176 
  177 version (X86_64) version = StaticallyHaveSSE;
  178 version (X86) version (OSX) version = StaticallyHaveSSE;
  179 
  180 version (StaticallyHaveSSE)
  181 {
  182     private enum bool haveSSE = true;
  183 }
  184 else version (X86)
  185 {
  186     static import core.cpuid;
  187     private alias haveSSE = core.cpuid.sse;
  188 }
  189 
  190 version (D_SoftFloat)
  191 {
  192     // Some soft float implementations may support IEEE floating flags.
  193     // The implementation here supports hardware flags only and is so currently
  194     // only available for supported targets.
  195 }
  196 else version (X86_Any)   version = IeeeFlagsSupport;
  197 else version (PPC_Any)   version = IeeeFlagsSupport;
  198 else version (RISCV_Any) version = IeeeFlagsSupport;
  199 else version (MIPS_Any)  version = IeeeFlagsSupport;
  200 else version (ARM_Any)   version = IeeeFlagsSupport;
  201 
  202 // Struct FloatingPointControl is only available if hardware FP units are available.
  203 version (D_HardFloat)
  204 {
  205     // FloatingPointControl.clearExceptions() depends on version IeeeFlagsSupport
  206     version (IeeeFlagsSupport) version = FloatingPointControlSupport;
  207 }
  208 
  209 version (StdUnittest) private
  210 {
  211     static if (real.sizeof > double.sizeof)
  212         enum uint useDigits = 16;
  213     else
  214         enum uint useDigits = 15;
  215 
  216     /******************************************
  217      * Compare floating point numbers to n decimal digits of precision.
  218      * Returns:
  219      *  1       match
  220      *  0       nomatch
  221      */
  222 
  223     private bool equalsDigit(real x, real y, uint ndigits) @safe nothrow @nogc
  224     {
  225         import core.stdc.stdio : sprintf;
  226 
  227         if (signbit(x) != signbit(y))
  228             return 0;
  229 
  230         if (isInfinity(x) && isInfinity(y))
  231             return 1;
  232         if (isInfinity(x) || isInfinity(y))
  233             return 0;
  234 
  235         if (isNaN(x) && isNaN(y))
  236             return 1;
  237         if (isNaN(x) || isNaN(y))
  238             return 0;
  239 
  240         char[30] bufx;
  241         char[30] bufy;
  242         assert(ndigits < bufx.length);
  243 
  244         int ix;
  245         int iy;
  246         version (CRuntime_Microsoft)
  247             alias real_t = double;
  248         else
  249             alias real_t = real;
  250 
  251         () @trusted {
  252             ix = sprintf(bufx.ptr, is(real_t == real) ? "%.*Lg" : "%.*g", ndigits, cast(real_t) x);
  253             iy = sprintf(bufy.ptr, is(real_t == real) ? "%.*Lg" : "%.*g", ndigits, cast(real_t) y);
  254         } ();
  255 
  256         assert(ix < bufx.length && ix > 0);
  257         assert(ix < bufy.length && ix > 0);
  258 
  259         return bufx[0 .. ix] == bufy[0 .. iy];
  260     }
  261 }
  262 
  263 
  264 // Values obtained from Wolfram Alpha. 116 bits ought to be enough for anybody.
  265 // Wolfram Alpha LLC. 2011. Wolfram|Alpha. http://www.wolframalpha.com/input/?i=e+in+base+16 (access July 6, 2011).
  266 enum real E =          0x1.5bf0a8b1457695355fb8ac404e7a8p+1L; /** e = 2.718281... */
  267 enum real LOG2T =      0x1.a934f0979a3715fc9257edfe9b5fbp+1L; /** $(SUB log, 2)10 = 3.321928... */
  268 enum real LOG2E =      0x1.71547652b82fe1777d0ffda0d23a8p+0L; /** $(SUB log, 2)e = 1.442695... */
  269 enum real LOG2 =       0x1.34413509f79fef311f12b35816f92p-2L; /** $(SUB log, 10)2 = 0.301029... */
  270 enum real LOG10E =     0x1.bcb7b1526e50e32a6ab7555f5a67cp-2L; /** $(SUB log, 10)e = 0.434294... */
  271 enum real LN2 =        0x1.62e42fefa39ef35793c7673007e5fp-1L; /** ln 2  = 0.693147... */
  272 enum real LN10 =       0x1.26bb1bbb5551582dd4adac5705a61p+1L; /** ln 10 = 2.302585... */
  273 enum real PI =         0x1.921fb54442d18469898cc51701b84p+1L; /** &pi; = 3.141592... */
  274 enum real PI_2 =       PI/2;                                  /** $(PI) / 2 = 1.570796... */
  275 enum real PI_4 =       PI/4;                                  /** $(PI) / 4 = 0.785398... */
  276 enum real M_1_PI =     0x1.45f306dc9c882a53f84eafa3ea69cp-2L; /** 1 / $(PI) = 0.318309... */
  277 enum real M_2_PI =     2*M_1_PI;                              /** 2 / $(PI) = 0.636619... */
  278 enum real M_2_SQRTPI = 0x1.20dd750429b6d11ae3a914fed7fd8p+0L; /** 2 / $(SQRT)$(PI) = 1.128379... */
  279 enum real SQRT2 =      0x1.6a09e667f3bcc908b2fb1366ea958p+0L; /** $(SQRT)2 = 1.414213... */
  280 enum real SQRT1_2 =    SQRT2/2;                               /** $(SQRT)$(HALF) = 0.707106... */
  281 // Note: Make sure the magic numbers in compiler backend for x87 match these.
  282 
  283 /***********************************
  284  * Calculates the absolute value of a number.
  285  *
  286  * Params:
  287  *     Num = (template parameter) type of number
  288  *       x = real number value
  289  *
  290  * Returns:
  291  *     The absolute value of the number. If floating-point or integral,
  292  *     the return type will be the same as the input.
  293  *
  294  * Limitations:
  295  *     Does not work correctly for signed intergal types and value `Num`.min.
  296  */
  297 auto abs(Num)(Num x) @nogc pure nothrow
  298 if ((is(immutable Num == immutable short) || is(immutable Num == immutable byte)) ||
  299     (is(typeof(Num.init >= 0)) && is(typeof(-Num.init))))
  300 {
  301     static if (isFloatingPoint!(Num))
  302         return fabs(x);
  303     else
  304     {
  305         static if (is(immutable Num == immutable short) || is(immutable Num == immutable byte))
  306             return x >= 0 ? x : cast(Num) -int(x);
  307         else
  308             return x >= 0 ? x : -x;
  309     }
  310 }
  311 
  312 /// ditto
  313 @safe pure nothrow @nogc unittest
  314 {
  315     assert(isIdentical(abs(-0.0L), 0.0L));
  316     assert(isNaN(abs(real.nan)));
  317     assert(abs(-real.infinity) == real.infinity);
  318     assert(abs(-56) == 56);
  319     assert(abs(2321312L)  == 2321312L);
  320 }
  321 
  322 @safe pure nothrow @nogc unittest
  323 {
  324     short s = -8;
  325     byte b = -8;
  326     assert(abs(s) == 8);
  327     assert(abs(b) == 8);
  328     immutable(byte) c = -8;
  329     assert(abs(c) == 8);
  330 }
  331 
  332 @safe pure nothrow @nogc unittest
  333 {
  334     import std.meta : AliasSeq;
  335     static foreach (T; AliasSeq!(float, double, real))
  336     {{
  337         T f = 3;
  338         assert(abs(f) == f);
  339         assert(abs(-f) == f);
  340     }}
  341 }
  342 
  343 // see https://issues.dlang.org/show_bug.cgi?id=20205
  344 // to avoid falling into the trap again
  345 @safe pure nothrow @nogc unittest
  346 {
  347     assert(50 - abs(-100) == -50);
  348 }
  349 
  350 // https://issues.dlang.org/show_bug.cgi?id=19162
  351 @safe unittest
  352 {
  353     struct Vector(T, int size)
  354     {
  355         T x, y, z;
  356     }
  357 
  358     static auto abs(T, int size)(auto ref const Vector!(T, size) v)
  359     {
  360         return v;
  361     }
  362     Vector!(int, 3) v;
  363     assert(abs(v) == v);
  364 }
  365 
  366 /***********************************
  367  * Returns cosine of x. x is in radians.
  368  *
  369  *      $(TABLE_SV
  370  *      $(TR $(TH x)                 $(TH cos(x)) $(TH invalid?))
  371  *      $(TR $(TD $(NAN))            $(TD $(NAN)) $(TD yes)     )
  372  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(NAN)) $(TD yes)     )
  373  *      )
  374  * Bugs:
  375  *      Results are undefined if |x| >= $(POWER 2,64).
  376  */
  377 
  378 real cos(real x) @safe pure nothrow @nogc { pragma(inline, true); return core.math.cos(x); }
  379 //FIXME
  380 ///ditto
  381 double cos(double x) @safe pure nothrow @nogc { return cos(cast(real) x); }
  382 //FIXME
  383 ///ditto
  384 float cos(float x) @safe pure nothrow @nogc { return cos(cast(real) x); }
  385 
  386 ///
  387 @safe unittest
  388 {
  389     assert(cos(0.0) == 1.0);
  390     assert(cos(1.0).approxEqual(0.540));
  391     assert(cos(3.0).approxEqual(-0.989));
  392 }
  393 
  394 @safe unittest
  395 {
  396     real function(real) pcos = &cos;
  397     assert(pcos != null);
  398 }
  399 
  400 /***********************************
  401  * Returns $(HTTP en.wikipedia.org/wiki/Sine, sine) of x. x is in $(HTTP en.wikipedia.org/wiki/Radian, radians).
  402  *
  403  *      $(TABLE_SV
  404  *      $(TH3 x           ,  sin(x)      ,  invalid?)
  405  *      $(TD3 $(NAN)      ,  $(NAN)      ,  yes     )
  406  *      $(TD3 $(PLUSMN)0.0,  $(PLUSMN)0.0,  no      )
  407  *      $(TD3 $(PLUSMNINF),  $(NAN)      ,  yes     )
  408  *      )
  409  *
  410  * Params:
  411  *      x = angle in radians (not degrees)
  412  * Returns:
  413  *      sine of x
  414  * See_Also:
  415  *      $(MYREF cos), $(MYREF tan), $(MYREF asin)
  416  * Bugs:
  417  *      Results are undefined if |x| >= $(POWER 2,64).
  418  */
  419 
  420 real sin(real x) @safe pure nothrow @nogc { pragma(inline, true); return core.math.sin(x); }
  421 //FIXME
  422 ///ditto
  423 double sin(double x) @safe pure nothrow @nogc { return sin(cast(real) x); }
  424 //FIXME
  425 ///ditto
  426 float sin(float x) @safe pure nothrow @nogc { return sin(cast(real) x); }
  427 
  428 ///
  429 @safe unittest
  430 {
  431     import std.math : sin, PI;
  432     import std.stdio : writefln;
  433 
  434     void someFunc()
  435     {
  436       real x = 30.0;
  437       auto result = sin(x * (PI / 180)); // convert degrees to radians
  438       writefln("The sine of %s degrees is %s", x, result);
  439     }
  440 }
  441 
  442 @safe unittest
  443 {
  444     real function(real) psin = &sin;
  445     assert(psin != null);
  446 }
  447 
  448 /****************************************************************************
  449  * Returns tangent of x. x is in radians.
  450  *
  451  *      $(TABLE_SV
  452  *      $(TR $(TH x)             $(TH tan(x))       $(TH invalid?))
  453  *      $(TR $(TD $(NAN))        $(TD $(NAN))       $(TD yes))
  454  *      $(TR $(TD $(PLUSMN)0.0)  $(TD $(PLUSMN)0.0) $(TD no))
  455  *      $(TR $(TD $(PLUSMNINF))  $(TD $(NAN))       $(TD yes))
  456  *      )
  457  */
  458 real tan(real x) @trusted pure nothrow @nogc // TODO: @safe
  459 {
  460     version (InlineAsm_X86_Any)
  461     {
  462         if (!__ctfe)
  463             return tanAsm(x);
  464     }
  465     return tanImpl(x);
  466 }
  467 
  468 /// ditto
  469 double tan(double x) @safe pure nothrow @nogc { return __ctfe ? cast(double) tan(cast(real) x) : tanImpl(x); }
  470 
  471 /// ditto
  472 float tan(float x) @safe pure nothrow @nogc { return __ctfe ? cast(float) tan(cast(real) x) : tanImpl(x); }
  473 
  474 ///
  475 @safe unittest
  476 {
  477     assert(isIdentical(tan(0.0), 0.0));
  478     assert(tan(PI).approxEqual(0));
  479     assert(tan(PI / 3).approxEqual(sqrt(3.0)));
  480 }
  481 
  482 version (InlineAsm_X86_Any)
  483 private real tanAsm(real x) @trusted pure nothrow @nogc
  484 {
  485     version (D_InlineAsm_X86)
  486     {
  487     asm pure nothrow @nogc
  488     {
  489         fld     x[EBP]                  ; // load theta
  490         fxam                            ; // test for oddball values
  491         fstsw   AX                      ;
  492         sahf                            ;
  493         jc      trigerr                 ; // x is NAN, infinity, or empty
  494                                           // 387's can handle subnormals
  495 SC18:   fptan                           ;
  496         fstsw   AX                      ;
  497         sahf                            ;
  498         jnp     Clear1                  ; // C2 = 1 (x is out of range)
  499 
  500         // Do argument reduction to bring x into range
  501         fldpi                           ;
  502         fxch                            ;
  503 SC17:   fprem1                          ;
  504         fstsw   AX                      ;
  505         sahf                            ;
  506         jp      SC17                    ;
  507         fstp    ST(1)                   ; // remove pi from stack
  508         jmp     SC18                    ;
  509 
  510 trigerr:
  511         jnp     Lret                    ; // if theta is NAN, return theta
  512         fstp    ST(0)                   ; // dump theta
  513     }
  514     return real.nan;
  515 
  516 Clear1: asm pure nothrow @nogc{
  517         fstp    ST(0)                   ; // dump X, which is always 1
  518     }
  519 
  520 Lret: {}
  521     }
  522     else version (D_InlineAsm_X86_64)
  523     {
  524         version (Win64)
  525         {
  526             asm pure nothrow @nogc
  527             {
  528                 fld     real ptr [RCX]  ; // load theta
  529             }
  530         }
  531         else
  532         {
  533             asm pure nothrow @nogc
  534             {
  535                 fld     x[RBP]          ; // load theta
  536             }
  537         }
  538     asm pure nothrow @nogc
  539     {
  540         fxam                            ; // test for oddball values
  541         fstsw   AX                      ;
  542         test    AH,1                    ;
  543         jnz     trigerr                 ; // x is NAN, infinity, or empty
  544                                           // 387's can handle subnormals
  545 SC18:   fptan                           ;
  546         fstsw   AX                      ;
  547         test    AH,4                    ;
  548         jz      Clear1                  ; // C2 = 1 (x is out of range)
  549 
  550         // Do argument reduction to bring x into range
  551         fldpi                           ;
  552         fxch                            ;
  553 SC17:   fprem1                          ;
  554         fstsw   AX                      ;
  555         test    AH,4                    ;
  556         jnz     SC17                    ;
  557         fstp    ST(1)                   ; // remove pi from stack
  558         jmp     SC18                    ;
  559 
  560 trigerr:
  561         test    AH,4                    ;
  562         jz      Lret                    ; // if theta is NAN, return theta
  563         fstp    ST(0)                   ; // dump theta
  564     }
  565     return real.nan;
  566 
  567 Clear1: asm pure nothrow @nogc{
  568         fstp    ST(0)                   ; // dump X, which is always 1
  569     }
  570 
  571 Lret: {}
  572     }
  573     else
  574         static assert(0);
  575 }
  576 
  577 private T tanImpl(T)(T x) @safe pure nothrow @nogc
  578 {
  579     // Coefficients for tan(x) and PI/4 split into three parts.
  580     enum realFormat = floatTraits!T.realFormat;
  581     static if (realFormat == RealFormat.ieeeQuadruple)
  582     {
  583         static immutable T[6] P = [
  584             2.883414728874239697964612246732416606301E10L,
  585             -2.307030822693734879744223131873392503321E9L,
  586             5.160188250214037865511600561074819366815E7L,
  587             -4.249691853501233575668486667664718192660E5L,
  588             1.272297782199996882828849455156962260810E3L,
  589             -9.889929415807650724957118893791829849557E-1L
  590         ];
  591         static immutable T[7] Q = [
  592             8.650244186622719093893836740197250197602E10L,
  593             -4.152206921457208101480801635640958361612E10L,
  594             2.758476078803232151774723646710890525496E9L,
  595             -5.733709132766856723608447733926138506824E7L,
  596             4.529422062441341616231663543669583527923E5L,
  597             -1.317243702830553658702531997959756728291E3L,
  598             1.0
  599         ];
  600 
  601         enum T P1 =
  602             7.853981633974483067550664827649598009884357452392578125E-1L;
  603         enum T P2 =
  604             2.8605943630549158983813312792950660807511260829685741796657E-18L;
  605         enum T P3 =
  606             2.1679525325309452561992610065108379921905808E-35L;
  607     }
  608     else static if (realFormat == RealFormat.ieeeExtended ||
  609                     realFormat == RealFormat.ieeeDouble)
  610     {
  611         static immutable T[3] P = [
  612            -1.7956525197648487798769E7L,
  613             1.1535166483858741613983E6L,
  614            -1.3093693918138377764608E4L,
  615         ];
  616         static immutable T[5] Q = [
  617            -5.3869575592945462988123E7L,
  618             2.5008380182335791583922E7L,
  619            -1.3208923444021096744731E6L,
  620             1.3681296347069295467845E4L,
  621             1.0000000000000000000000E0L,
  622         ];
  623 
  624         enum T P1 = 7.853981554508209228515625E-1L;
  625         enum T P2 = 7.946627356147928367136046290398E-9L;
  626         enum T P3 = 3.061616997868382943065164830688E-17L;
  627     }
  628     else static if (realFormat == RealFormat.ieeeSingle)
  629     {
  630         static immutable T[6] P = [
  631             3.33331568548E-1,
  632             1.33387994085E-1,
  633             5.34112807005E-2,
  634             2.44301354525E-2,
  635             3.11992232697E-3,
  636             9.38540185543E-3,
  637         ];
  638 
  639         enum T P1 = 0.78515625;
  640         enum T P2 = 2.4187564849853515625E-4;
  641         enum T P3 = 3.77489497744594108E-8;
  642     }
  643     else
  644         static assert(0, "no coefficients for tan()");
  645 
  646     // Special cases.
  647     if (x == cast(T) 0.0 || isNaN(x))
  648         return x;
  649     if (isInfinity(x))
  650         return T.nan;
  651 
  652     // Make argument positive but save the sign.
  653     bool sign = false;
  654     if (signbit(x))
  655     {
  656         sign = true;
  657         x = -x;
  658     }
  659 
  660     // Compute x mod PI/4.
  661     static if (realFormat == RealFormat.ieeeSingle)
  662     {
  663         enum T FOPI = 4 / PI;
  664         int j = cast(int) (FOPI * x);
  665         T y = j;
  666         T z;
  667     }
  668     else
  669     {
  670         T y = floor(x / cast(T) PI_4);
  671         // Strip high bits of integer part.
  672         enum T highBitsFactor = (realFormat == RealFormat.ieeeDouble ? 0x1p3 : 0x1p4);
  673         enum T highBitsInv = 1.0 / highBitsFactor;
  674         T z = y * highBitsInv;
  675         // Compute y - 2^numHighBits * (y / 2^numHighBits).
  676         z = y - highBitsFactor * floor(z);
  677 
  678         // Integer and fraction part modulo one octant.
  679         int j = cast(int)(z);
  680     }
  681 
  682     // Map zeros and singularities to origin.
  683     if (j & 1)
  684     {
  685         j += 1;
  686         y += cast(T) 1.0;
  687     }
  688 
  689     z = ((x - y * P1) - y * P2) - y * P3;
  690     const T zz = z * z;
  691 
  692     enum T zzThreshold = (realFormat == RealFormat.ieeeSingle ? 1.0e-4L :
  693                           realFormat == RealFormat.ieeeDouble ? 1.0e-14L : 1.0e-20L);
  694     if (zz > zzThreshold)
  695     {
  696         static if (realFormat == RealFormat.ieeeSingle)
  697             y = z + z * (zz * poly(zz, P));
  698         else
  699             y = z + z * (zz * poly(zz, P) / poly(zz, Q));
  700     }
  701     else
  702         y = z;
  703 
  704     if (j & 2)
  705         y = (cast(T) -1.0) / y;
  706 
  707     return (sign) ? -y : y;
  708 }
  709 
  710 @safe @nogc nothrow unittest
  711 {
  712     static void testTan(T)()
  713     {
  714         // ±0
  715         const T zero = 0.0;
  716         assert(isIdentical(tan(zero), zero));
  717         assert(isIdentical(tan(-zero), -zero));
  718         // ±∞
  719         const T inf = T.infinity;
  720         assert(isNaN(tan(inf)));
  721         assert(isNaN(tan(-inf)));
  722         // NaN
  723         const T specialNaN = NaN(0x0123L);
  724         assert(isIdentical(tan(specialNaN), specialNaN));
  725 
  726         static immutable T[2][] vals =
  727         [
  728             // angle, tan
  729             [   .5,  .5463024898],
  730             [   1,   1.557407725],
  731             [   1.5, 14.10141995],
  732             [   2,  -2.185039863],
  733             [   2.5,-.7470222972],
  734             [   3,  -.1425465431],
  735             [   3.5, .3745856402],
  736             [   4,   1.157821282],
  737             [   4.5, 4.637332055],
  738             [   5,  -3.380515006],
  739             [   5.5,-.9955840522],
  740             [   6,  -.2910061914],
  741             [   6.5, .2202772003],
  742             [   10,  .6483608275],
  743 
  744             // special angles
  745             [   PI_4,   1],
  746             //[   PI_2,   T.infinity], // PI_2 is not _exactly_ pi/2.
  747             [   3*PI_4, -1],
  748             [   PI,     0],
  749             [   5*PI_4, 1],
  750             //[   3*PI_2, -T.infinity],
  751             [   7*PI_4, -1],
  752             [   2*PI,   0],
  753          ];
  754 
  755         foreach (ref val; vals)
  756         {
  757             T x = val[0];
  758             T r = val[1];
  759             T t = tan(x);
  760 
  761             //printf("tan(%Lg) = %Lg, should be %Lg\n", cast(real) x, cast(real) t, cast(real) r);
  762             assert(approxEqual(r, t));
  763 
  764             x = -x;
  765             r = -r;
  766             t = tan(x);
  767             //printf("tan(%Lg) = %Lg, should be %Lg\n", cast(real) x, cast(real) t, cast(real) r);
  768             assert(approxEqual(r, t));
  769         }
  770     }
  771 
  772     import std.meta : AliasSeq;
  773     foreach (T; AliasSeq!(real, double, float))
  774         testTan!T();
  775 
  776     assert(equalsDigit(tan(PI / 3), std.math.sqrt(3.0L), useDigits));
  777 }
  778 
  779 /***************
  780  * Calculates the arc cosine of x,
  781  * returning a value ranging from 0 to $(PI).
  782  *
  783  *      $(TABLE_SV
  784  *      $(TR $(TH x)         $(TH acos(x)) $(TH invalid?))
  785  *      $(TR $(TD $(GT)1.0)  $(TD $(NAN))  $(TD yes))
  786  *      $(TR $(TD $(LT)-1.0) $(TD $(NAN))  $(TD yes))
  787  *      $(TR $(TD $(NAN))    $(TD $(NAN))  $(TD yes))
  788  *  )
  789  */
  790 real acos(real x) @safe pure nothrow @nogc
  791 {
  792     return atan2(sqrt(1-x*x), x);
  793 }
  794 
  795 /// ditto
  796 double acos(double x) @safe pure nothrow @nogc { return acos(cast(real) x); }
  797 
  798 /// ditto
  799 float acos(float x) @safe pure nothrow @nogc  { return acos(cast(real) x); }
  800 
  801 ///
  802 @safe unittest
  803 {
  804     assert(acos(0.0).approxEqual(1.570));
  805     assert(acos(0.5).approxEqual(std.math.PI / 3));
  806     assert(acos(PI).isNaN);
  807 }
  808 
  809 @safe @nogc nothrow unittest
  810 {
  811     assert(equalsDigit(acos(0.5), std.math.PI / 3, useDigits));
  812 }
  813 
  814 /***************
  815  * Calculates the arc sine of x,
  816  * returning a value ranging from -$(PI)/2 to $(PI)/2.
  817  *
  818  *      $(TABLE_SV
  819  *      $(TR $(TH x)            $(TH asin(x))      $(TH invalid?))
  820  *      $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no))
  821  *      $(TR $(TD $(GT)1.0)     $(TD $(NAN))       $(TD yes))
  822  *      $(TR $(TD $(LT)-1.0)    $(TD $(NAN))       $(TD yes))
  823  *  )
  824  */
  825 real asin(real x) @safe pure nothrow @nogc
  826 {
  827     return atan2(x, sqrt(1-x*x));
  828 }
  829 
  830 /// ditto
  831 double asin(double x) @safe pure nothrow @nogc { return asin(cast(real) x); }
  832 
  833 /// ditto
  834 float asin(float x) @safe pure nothrow @nogc  { return asin(cast(real) x); }
  835 
  836 ///
  837 @safe unittest
  838 {
  839     assert(isIdentical(asin(0.0), 0.0));
  840     assert(asin(0.5).approxEqual(PI / 6));
  841     assert(asin(PI).isNaN);
  842 }
  843 
  844 @safe @nogc nothrow unittest
  845 {
  846     assert(equalsDigit(asin(0.5), PI / 6, useDigits));
  847 }
  848 
  849 /***************
  850  * Calculates the arc tangent of x,
  851  * returning a value ranging from -$(PI)/2 to $(PI)/2.
  852  *
  853  *      $(TABLE_SV
  854  *      $(TR $(TH x)                 $(TH atan(x))      $(TH invalid?))
  855  *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(PLUSMN)0.0) $(TD no))
  856  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(NAN))       $(TD yes))
  857  *  )
  858  */
  859 real atan(real x) @safe pure nothrow @nogc
  860 {
  861     version (InlineAsm_X86_Any)
  862     {
  863         if (!__ctfe)
  864             return atan2Asm(x, 1.0L);
  865     }
  866     return atanImpl(x);
  867 }
  868 
  869 /// ditto
  870 double atan(double x) @safe pure nothrow @nogc { return __ctfe ? cast(double) atan(cast(real) x) : atanImpl(x); }
  871 
  872 /// ditto
  873 float atan(float x) @safe pure nothrow @nogc { return __ctfe ? cast(float) atan(cast(real) x) : atanImpl(x); }
  874 
  875 ///
  876 @safe unittest
  877 {
  878     assert(isIdentical(atan(0.0), 0.0));
  879     assert(atan(sqrt(3.0)).approxEqual(PI / 3));
  880 }
  881 
  882 private T atanImpl(T)(T x) @safe pure nothrow @nogc
  883 {
  884     // Coefficients for atan(x)
  885     enum realFormat = floatTraits!T.realFormat;
  886     static if (realFormat == RealFormat.ieeeQuadruple)
  887     {
  888         static immutable T[9] P = [
  889             -6.880597774405940432145577545328795037141E2L,
  890             -2.514829758941713674909996882101723647996E3L,
  891             -3.696264445691821235400930243493001671932E3L,
  892             -2.792272753241044941703278827346430350236E3L,
  893             -1.148164399808514330375280133523543970854E3L,
  894             -2.497759878476618348858065206895055957104E2L,
  895             -2.548067867495502632615671450650071218995E1L,
  896             -8.768423468036849091777415076702113400070E-1L,
  897             -6.635810778635296712545011270011752799963E-4L
  898         ];
  899         static immutable T[9] Q = [
  900             2.064179332321782129643673263598686441900E3L,
  901             8.782996876218210302516194604424986107121E3L,
  902             1.547394317752562611786521896296215170819E4L,
  903             1.458510242529987155225086911411015961174E4L,
  904             7.928572347062145288093560392463784743935E3L,
  905             2.494680540950601626662048893678584497900E3L,
  906             4.308348370818927353321556740027020068897E2L,
  907             3.566239794444800849656497338030115886153E1L,
  908             1.0
  909         ];
  910     }
  911     else static if (realFormat == RealFormat.ieeeExtended)
  912     {
  913         static immutable T[5] P = [
  914            -5.0894116899623603312185E1L,
  915            -9.9988763777265819915721E1L,
  916            -6.3976888655834347413154E1L,
  917            -1.4683508633175792446076E1L,
  918            -8.6863818178092187535440E-1L,
  919         ];
  920         static immutable T[6] Q = [
  921             1.5268235069887081006606E2L,
  922             3.9157570175111990631099E2L,
  923             3.6144079386152023162701E2L,
  924             1.4399096122250781605352E2L,
  925             2.2981886733594175366172E1L,
  926             1.0000000000000000000000E0L,
  927         ];
  928     }
  929     else static if (realFormat == RealFormat.ieeeDouble)
  930     {
  931         static immutable T[5] P = [
  932            -6.485021904942025371773E1L,
  933            -1.228866684490136173410E2L,
  934            -7.500855792314704667340E1L,
  935            -1.615753718733365076637E1L,
  936            -8.750608600031904122785E-1L,
  937         ];
  938         static immutable T[6] Q = [
  939             1.945506571482613964425E2L,
  940             4.853903996359136964868E2L,
  941             4.328810604912902668951E2L,
  942             1.650270098316988542046E2L,
  943             2.485846490142306297962E1L,
  944             1.000000000000000000000E0L,
  945         ];
  946 
  947         enum T MOREBITS = 6.123233995736765886130E-17L;
  948     }
  949     else static if (realFormat == RealFormat.ieeeSingle)
  950     {
  951         static immutable T[4] P = [
  952            -3.33329491539E-1,
  953             1.99777106478E-1,
  954            -1.38776856032E-1,
  955             8.05374449538E-2,
  956         ];
  957     }
  958     else
  959         static assert(0, "no coefficients for atan()");
  960 
  961     // tan(PI/8)
  962     enum T TAN_PI_8 = 0.414213562373095048801688724209698078569672L;
  963     // tan(3 * PI/8)
  964     enum T TAN3_PI_8 = 2.414213562373095048801688724209698078569672L;
  965 
  966     // Special cases.
  967     if (x == cast(T) 0.0)
  968         return x;
  969     if (isInfinity(x))
  970         return copysign(cast(T) PI_2, x);
  971 
  972     // Make argument positive but save the sign.
  973     bool sign = false;
  974     if (signbit(x))
  975     {
  976         sign = true;
  977         x = -x;
  978     }
  979 
  980     static if (realFormat == RealFormat.ieeeDouble) // special case for double precision
  981     {
  982         short flag = 0;
  983         T y;
  984         if (x > TAN3_PI_8)
  985         {
  986             y = PI_2;
  987             flag = 1;
  988             x = -(1.0 / x);
  989         }
  990         else if (x <= 0.66)
  991         {
  992             y = 0.0;
  993         }
  994         else
  995         {
  996             y = PI_4;
  997             flag = 2;
  998             x = (x - 1.0)/(x + 1.0);
  999         }
 1000 
 1001         T z = x * x;
 1002         z = z * poly(z, P) / poly(z, Q);
 1003         z = x * z + x;
 1004         if (flag == 2)
 1005             z += 0.5 * MOREBITS;
 1006         else if (flag == 1)
 1007             z += MOREBITS;
 1008         y = y + z;
 1009     }
 1010     else
 1011     {
 1012         // Range reduction.
 1013         T y;
 1014         if (x > TAN3_PI_8)
 1015         {
 1016             y = PI_2;
 1017             x = -((cast(T) 1.0) / x);
 1018         }
 1019         else if (x > TAN_PI_8)
 1020         {
 1021             y = PI_4;
 1022             x = (x - cast(T) 1.0)/(x + cast(T) 1.0);
 1023         }
 1024         else
 1025             y = 0.0;
 1026 
 1027         // Rational form in x^^2.
 1028         const T z = x * x;
 1029         static if (realFormat == RealFormat.ieeeSingle)
 1030             y += poly(z, P) * z * x + x;
 1031         else
 1032             y = y + (poly(z, P) / poly(z, Q)) * z * x + x;
 1033     }
 1034 
 1035     return (sign) ? -y : y;
 1036 }
 1037 
 1038 @safe @nogc nothrow unittest
 1039 {
 1040     static void testAtan(T)()
 1041     {
 1042         // ±0
 1043         const T zero = 0.0;
 1044         assert(isIdentical(atan(zero), zero));
 1045         assert(isIdentical(atan(-zero), -zero));
 1046         // ±∞
 1047         const T inf = T.infinity;
 1048         assert(approxEqual(atan(inf), cast(T) PI_2));
 1049         assert(approxEqual(atan(-inf), cast(T) -PI_2));
 1050         // NaN
 1051         const T specialNaN = NaN(0x0123L);
 1052         assert(isIdentical(atan(specialNaN), specialNaN));
 1053 
 1054         static immutable T[2][] vals =
 1055         [
 1056             // x, atan(x)
 1057             [ 0.25, 0.2449786631 ],
 1058             [ 0.5,  0.4636476090 ],
 1059             [ 1,    PI_4         ],
 1060             [ 1.5,  0.9827937232 ],
 1061             [ 10,   1.4711276743 ],
 1062         ];
 1063 
 1064         foreach (ref val; vals)
 1065         {
 1066             T x = val[0];
 1067             T r = val[1];
 1068             T a = atan(x);
 1069 
 1070             //printf("atan(%Lg) = %Lg, should be %Lg\n", cast(real) x, cast(real) a, cast(real) r);
 1071             assert(approxEqual(r, a));
 1072 
 1073             x = -x;
 1074             r = -r;
 1075             a = atan(x);
 1076             //printf("atan(%Lg) = %Lg, should be %Lg\n", cast(real) x, cast(real) a, cast(real) r);
 1077             assert(approxEqual(r, a));
 1078         }
 1079     }
 1080 
 1081     import std.meta : AliasSeq;
 1082     foreach (T; AliasSeq!(real, double, float))
 1083         testAtan!T();
 1084 
 1085     assert(equalsDigit(atan(std.math.sqrt(3.0L)), PI / 3, useDigits));
 1086 }
 1087 
 1088 /***************
 1089  * Calculates the arc tangent of y / x,
 1090  * returning a value ranging from -$(PI) to $(PI).
 1091  *
 1092  *      $(TABLE_SV
 1093  *      $(TR $(TH y)                 $(TH x)            $(TH atan(y, x)))
 1094  *      $(TR $(TD $(NAN))            $(TD anything)     $(TD $(NAN)) )
 1095  *      $(TR $(TD anything)          $(TD $(NAN))       $(TD $(NAN)) )
 1096  *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(GT)0.0)     $(TD $(PLUSMN)0.0) )
 1097  *      $(TR $(TD $(PLUSMN)0.0)      $(TD +0.0)         $(TD $(PLUSMN)0.0) )
 1098  *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(LT)0.0)     $(TD $(PLUSMN)$(PI)))
 1099  *      $(TR $(TD $(PLUSMN)0.0)      $(TD -0.0)         $(TD $(PLUSMN)$(PI)))
 1100  *      $(TR $(TD $(GT)0.0)          $(TD $(PLUSMN)0.0) $(TD $(PI)/2) )
 1101  *      $(TR $(TD $(LT)0.0)          $(TD $(PLUSMN)0.0) $(TD -$(PI)/2) )
 1102  *      $(TR $(TD $(GT)0.0)          $(TD $(INFIN))     $(TD $(PLUSMN)0.0) )
 1103  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD anything)     $(TD $(PLUSMN)$(PI)/2))
 1104  *      $(TR $(TD $(GT)0.0)          $(TD -$(INFIN))    $(TD $(PLUSMN)$(PI)) )
 1105  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(INFIN))     $(TD $(PLUSMN)$(PI)/4))
 1106  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD -$(INFIN))    $(TD $(PLUSMN)3$(PI)/4))
 1107  *      )
 1108  */
 1109 real atan2(real y, real x) @trusted pure nothrow @nogc // TODO: @safe
 1110 {
 1111     version (InlineAsm_X86_Any)
 1112     {
 1113         if (!__ctfe)
 1114             return atan2Asm(y, x);
 1115     }
 1116     return atan2Impl(y, x);
 1117 }
 1118 
 1119 /// ditto
 1120 double atan2(double y, double x) @safe pure nothrow @nogc
 1121 {
 1122     return __ctfe ? cast(double) atan2(cast(real) y, cast(real) x) : atan2Impl(y, x);
 1123 }
 1124 
 1125 /// ditto
 1126 float atan2(float y, float x) @safe pure nothrow @nogc
 1127 {
 1128     return __ctfe ? cast(float) atan2(cast(real) y, cast(real) x) : atan2Impl(y, x);
 1129 }
 1130 
 1131 ///
 1132 @safe unittest
 1133 {
 1134     assert(atan2(1.0, sqrt(3.0)).approxEqual(PI / 6));
 1135 }
 1136 
 1137 version (InlineAsm_X86_Any)
 1138 private real atan2Asm(real y, real x) @trusted pure nothrow @nogc
 1139 {
 1140     version (Win64)
 1141     {
 1142         asm pure nothrow @nogc {
 1143             naked;
 1144             fld real ptr [RDX]; // y
 1145             fld real ptr [RCX]; // x
 1146             fpatan;
 1147             ret;
 1148         }
 1149     }
 1150     else
 1151     {
 1152         asm pure nothrow @nogc {
 1153             fld y;
 1154             fld x;
 1155             fpatan;
 1156         }
 1157     }
 1158 }
 1159 
 1160 private T atan2Impl(T)(T y, T x) @safe pure nothrow @nogc
 1161 {
 1162     // Special cases.
 1163     if (isNaN(x) || isNaN(y))
 1164         return T.nan;
 1165     if (y == cast(T) 0.0)
 1166     {
 1167         if (x >= 0 && !signbit(x))
 1168             return copysign(0, y);
 1169         else
 1170             return copysign(cast(T) PI, y);
 1171     }
 1172     if (x == cast(T) 0.0)
 1173         return copysign(cast(T) PI_2, y);
 1174     if (isInfinity(x))
 1175     {
 1176         if (signbit(x))
 1177         {
 1178             if (isInfinity(y))
 1179                 return copysign(3 * cast(T) PI_4, y);
 1180             else
 1181                 return copysign(cast(T) PI, y);
 1182         }
 1183         else
 1184         {
 1185             if (isInfinity(y))
 1186                 return copysign(cast(T) PI_4, y);
 1187             else
 1188                 return copysign(cast(T) 0.0, y);
 1189         }
 1190     }
 1191     if (isInfinity(y))
 1192         return copysign(cast(T) PI_2, y);
 1193 
 1194     // Call atan and determine the quadrant.
 1195     T z = atan(y / x);
 1196 
 1197     if (signbit(x))
 1198     {
 1199         if (signbit(y))
 1200             z = z - cast(T) PI;
 1201         else
 1202             z = z + cast(T) PI;
 1203     }
 1204 
 1205     if (z == cast(T) 0.0)
 1206         return copysign(z, y);
 1207 
 1208     return z;
 1209 }
 1210 
 1211 @safe @nogc nothrow unittest
 1212 {
 1213     static void testAtan2(T)()
 1214     {
 1215         // NaN
 1216         const T nan = T.nan;
 1217         assert(isNaN(atan2(nan, cast(T) 1)));
 1218         assert(isNaN(atan2(cast(T) 1, nan)));
 1219 
 1220         const T inf = T.infinity;
 1221         static immutable T[3][] vals =
 1222         [
 1223             // y, x, atan2(y, x)
 1224 
 1225             // ±0
 1226             [  0.0,  1.0,  0.0 ],
 1227             [ -0.0,  1.0, -0.0 ],
 1228             [  0.0,  0.0,  0.0 ],
 1229             [ -0.0,  0.0, -0.0 ],
 1230             [  0.0, -1.0,  PI ],
 1231             [ -0.0, -1.0, -PI ],
 1232             [  0.0, -0.0,  PI ],
 1233             [ -0.0, -0.0, -PI ],
 1234             [  1.0,  0.0,  PI_2 ],
 1235             [  1.0, -0.0,  PI_2 ],
 1236             [ -1.0,  0.0, -PI_2 ],
 1237             [ -1.0, -0.0, -PI_2 ],
 1238 
 1239             // ±∞
 1240             [  1.0,  inf,  0.0 ],
 1241             [ -1.0,  inf, -0.0 ],
 1242             [  1.0, -inf,  PI ],
 1243             [ -1.0, -inf, -PI ],
 1244             [  inf,  1.0,  PI_2 ],
 1245             [  inf, -1.0,  PI_2 ],
 1246             [ -inf,  1.0, -PI_2 ],
 1247             [ -inf, -1.0, -PI_2 ],
 1248             [  inf,  inf,  PI_4 ],
 1249             [ -inf,  inf, -PI_4 ],
 1250             [  inf, -inf,  3 * PI_4 ],
 1251             [ -inf, -inf, -3 * PI_4 ],
 1252 
 1253             [  1.0,  1.0,  PI_4 ],
 1254             [ -2.0,  2.0, -PI_4 ],
 1255             [  3.0, -3.0,  3 * PI_4 ],
 1256             [ -4.0, -4.0, -3 * PI_4 ],
 1257 
 1258             [  0.75,  0.25,   1.249045772398 ],
 1259             [ -0.5,   0.375, -0.927295218002 ],
 1260             [  0.5,  -0.125,  1.815774989922 ],
 1261             [ -0.75, -0.5,   -2.158798930342 ],
 1262         ];
 1263 
 1264         foreach (ref val; vals)
 1265         {
 1266             const T y = val[0];
 1267             const T x = val[1];
 1268             const T r = val[2];
 1269             const T a = atan2(y, x);
 1270 
 1271             //printf("atan2(%Lg, %Lg) = %Lg, should be %Lg\n", cast(real) y, cast(real) x, cast(real) a, cast(real) r);
 1272             if (r == 0)
 1273                 assert(isIdentical(r, a)); // check sign
 1274             else
 1275                 assert(approxEqual(r, a));
 1276         }
 1277     }
 1278 
 1279     import std.meta : AliasSeq;
 1280     foreach (T; AliasSeq!(real, double, float))
 1281         testAtan2!T();
 1282 
 1283     assert(equalsDigit(atan2(1.0L, std.math.sqrt(3.0L)), PI / 6, useDigits));
 1284 }
 1285 
 1286 /***********************************
 1287  * Calculates the hyperbolic cosine of x.
 1288  *
 1289  *      $(TABLE_SV
 1290  *      $(TR $(TH x)                 $(TH cosh(x))      $(TH invalid?))
 1291  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)0.0) $(TD no) )
 1292  *      )
 1293  */
 1294 real cosh(real x) @safe pure nothrow @nogc
 1295 {
 1296     //  cosh = (exp(x)+exp(-x))/2.
 1297     // The naive implementation works correctly.
 1298     const real y = exp(x);
 1299     return (y + 1.0/y) * 0.5;
 1300 }
 1301 
 1302 /// ditto
 1303 double cosh(double x) @safe pure nothrow @nogc { return cosh(cast(real) x); }
 1304 
 1305 /// ditto
 1306 float cosh(float x) @safe pure nothrow @nogc  { return cosh(cast(real) x); }
 1307 
 1308 ///
 1309 @safe unittest
 1310 {
 1311     assert(cosh(0.0) == 1.0);
 1312     assert(cosh(1.0).approxEqual((E + 1.0 / E) / 2));
 1313 }
 1314 
 1315 @safe @nogc nothrow unittest
 1316 {
 1317     assert(equalsDigit(cosh(1.0), (E + 1.0 / E) / 2, useDigits));
 1318 }
 1319 
 1320 /***********************************
 1321  * Calculates the hyperbolic sine of x.
 1322  *
 1323  *      $(TABLE_SV
 1324  *      $(TR $(TH x)                 $(TH sinh(x))           $(TH invalid?))
 1325  *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(PLUSMN)0.0)      $(TD no))
 1326  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)$(INFIN)) $(TD no))
 1327  *      )
 1328  */
 1329 real sinh(real x) @safe pure nothrow @nogc
 1330 {
 1331     //  sinh(x) =  (exp(x)-exp(-x))/2;
 1332     // Very large arguments could cause an overflow, but
 1333     // the maximum value of x for which exp(x) + exp(-x)) != exp(x)
 1334     // is x = 0.5 * (real.mant_dig) * LN2. // = 22.1807 for real80.
 1335     if (fabs(x) > real.mant_dig * LN2)
 1336     {
 1337         return copysign(0.5 * exp(fabs(x)), x);
 1338     }
 1339 
 1340     const real y = expm1(x);
 1341     return 0.5 * y / (y+1) * (y+2);
 1342 }
 1343 
 1344 /// ditto
 1345 double sinh(double x) @safe pure nothrow @nogc { return sinh(cast(real) x); }
 1346 
 1347 /// ditto
 1348 float sinh(float x) @safe pure nothrow @nogc  { return sinh(cast(real) x); }
 1349 
 1350 ///
 1351 @safe unittest
 1352 {
 1353     assert(isIdentical(sinh(0.0), 0.0));
 1354     assert(sinh(1.0).approxEqual((E - 1.0 / E) / 2));
 1355 }
 1356 
 1357 @safe @nogc nothrow unittest
 1358 {
 1359     assert(equalsDigit(sinh(1.0), (E - 1.0 / E) / 2, useDigits));
 1360 }
 1361 
 1362 /***********************************
 1363  * Calculates the hyperbolic tangent of x.
 1364  *
 1365  *      $(TABLE_SV
 1366  *      $(TR $(TH x)                 $(TH tanh(x))      $(TH invalid?))
 1367  *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(PLUSMN)0.0) $(TD no) )
 1368  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)1.0) $(TD no))
 1369  *      )
 1370  */
 1371 real tanh(real x) @safe pure nothrow @nogc
 1372 {
 1373     //  tanh(x) = (exp(x) - exp(-x))/(exp(x)+exp(-x))
 1374     if (fabs(x) > real.mant_dig * LN2)
 1375     {
 1376         return copysign(1, x);
 1377     }
 1378 
 1379     const real y = expm1(2*x);
 1380     return y / (y + 2);
 1381 }
 1382 
 1383 /// ditto
 1384 double tanh(double x) @safe pure nothrow @nogc { return tanh(cast(real) x); }
 1385 
 1386 /// ditto
 1387 float tanh(float x) @safe pure nothrow @nogc { return tanh(cast(real) x); }
 1388 
 1389 ///
 1390 @safe unittest
 1391 {
 1392     assert(isIdentical(tanh(0.0), 0.0));
 1393     assert(tanh(1.0).approxEqual(sinh(1.0) / cosh(1.0)));
 1394 }
 1395 
 1396 @safe @nogc nothrow unittest
 1397 {
 1398     assert(equalsDigit(tanh(1.0), sinh(1.0) / cosh(1.0), 15));
 1399 }
 1400 
 1401 /***********************************
 1402  * Calculates the inverse hyperbolic cosine of x.
 1403  *
 1404  *  Mathematically, acosh(x) = log(x + sqrt( x*x - 1))
 1405  *
 1406  * $(TABLE_DOMRG
 1407  *    $(DOMAIN 1..$(INFIN)),
 1408  *    $(RANGE  0..$(INFIN))
 1409  * )
 1410  *
 1411  *  $(TABLE_SV
 1412  *    $(SVH  x,     acosh(x) )
 1413  *    $(SV  $(NAN), $(NAN) )
 1414  *    $(SV  $(LT)1,     $(NAN) )
 1415  *    $(SV  1,      0       )
 1416  *    $(SV  +$(INFIN),+$(INFIN))
 1417  *  )
 1418  */
 1419 real acosh(real x) @safe pure nothrow @nogc
 1420 {
 1421     if (x > 1/real.epsilon)
 1422         return LN2 + log(x);
 1423     else
 1424         return log(x + sqrt(x*x - 1));
 1425 }
 1426 
 1427 /// ditto
 1428 double acosh(double x) @safe pure nothrow @nogc { return acosh(cast(real) x); }
 1429 
 1430 /// ditto
 1431 float acosh(float x) @safe pure nothrow @nogc  { return acosh(cast(real) x); }
 1432 
 1433 ///
 1434 @safe @nogc nothrow unittest
 1435 {
 1436     assert(isNaN(acosh(0.9)));
 1437     assert(isNaN(acosh(real.nan)));
 1438     assert(isIdentical(acosh(1.0), 0.0));
 1439     assert(acosh(real.infinity) == real.infinity);
 1440     assert(isNaN(acosh(0.5)));
 1441 }
 1442 
 1443 @safe @nogc nothrow unittest
 1444 {
 1445     assert(equalsDigit(acosh(cosh(3.0)), 3, useDigits));
 1446 }
 1447 
 1448 /***********************************
 1449  * Calculates the inverse hyperbolic sine of x.
 1450  *
 1451  *  Mathematically,
 1452  *  ---------------
 1453  *  asinh(x) =  log( x + sqrt( x*x + 1 )) // if x >= +0
 1454  *  asinh(x) = -log(-x + sqrt( x*x + 1 )) // if x <= -0
 1455  *  -------------
 1456  *
 1457  *    $(TABLE_SV
 1458  *    $(SVH x,                asinh(x)       )
 1459  *    $(SV  $(NAN),           $(NAN)         )
 1460  *    $(SV  $(PLUSMN)0,       $(PLUSMN)0      )
 1461  *    $(SV  $(PLUSMN)$(INFIN),$(PLUSMN)$(INFIN))
 1462  *    )
 1463  */
 1464 real asinh(real x) @safe pure nothrow @nogc
 1465 {
 1466     return (fabs(x) > 1 / real.epsilon)
 1467        // beyond this point, x*x + 1 == x*x
 1468        ?  copysign(LN2 + log(fabs(x)), x)
 1469        // sqrt(x*x + 1) ==  1 + x * x / ( 1 + sqrt(x*x + 1) )
 1470        : copysign(log1p(fabs(x) + x*x / (1 + sqrt(x*x + 1)) ), x);
 1471 }
 1472 
 1473 /// ditto
 1474 double asinh(double x) @safe pure nothrow @nogc { return asinh(cast(real) x); }
 1475 
 1476 /// ditto
 1477 float asinh(float x) @safe pure nothrow @nogc { return asinh(cast(real) x); }
 1478 
 1479 ///
 1480 @safe @nogc nothrow unittest
 1481 {
 1482     assert(isIdentical(asinh(0.0), 0.0));
 1483     assert(isIdentical(asinh(-0.0), -0.0));
 1484     assert(asinh(real.infinity) == real.infinity);
 1485     assert(asinh(-real.infinity) == -real.infinity);
 1486     assert(isNaN(asinh(real.nan)));
 1487 }
 1488 
 1489 @safe unittest
 1490 {
 1491     assert(equalsDigit(asinh(sinh(3.0)), 3, useDigits));
 1492 }
 1493 
 1494 /***********************************
 1495  * Calculates the inverse hyperbolic tangent of x,
 1496  * returning a value from ranging from -1 to 1.
 1497  *
 1498  * Mathematically, atanh(x) = log( (1+x)/(1-x) ) / 2
 1499  *
 1500  * $(TABLE_DOMRG
 1501  *    $(DOMAIN -$(INFIN)..$(INFIN)),
 1502  *    $(RANGE  -1 .. 1)
 1503  * )
 1504  * $(BR)
 1505  * $(TABLE_SV
 1506  *    $(SVH  x,     acosh(x) )
 1507  *    $(SV  $(NAN), $(NAN) )
 1508  *    $(SV  $(PLUSMN)0, $(PLUSMN)0)
 1509  *    $(SV  -$(INFIN), -0)
 1510  * )
 1511  */
 1512 real atanh(real x) @safe pure nothrow @nogc
 1513 {
 1514     // log( (1+x)/(1-x) ) == log ( 1 + (2*x)/(1-x) )
 1515     return  0.5 * log1p( 2 * x / (1 - x) );
 1516 }
 1517 
 1518 /// ditto
 1519 double atanh(double x) @safe pure nothrow @nogc { return atanh(cast(real) x); }
 1520 
 1521 /// ditto
 1522 float atanh(float x) @safe pure nothrow @nogc { return atanh(cast(real) x); }
 1523 
 1524 ///
 1525 @safe @nogc nothrow unittest
 1526 {
 1527     assert(isIdentical(atanh(0.0), 0.0));
 1528     assert(isIdentical(atanh(-0.0),-0.0));
 1529     assert(isNaN(atanh(real.nan)));
 1530     assert(isNaN(atanh(-real.infinity)));
 1531     assert(atanh(0.0) == 0);
 1532 }
 1533 
 1534 @safe unittest
 1535 {
 1536     assert(equalsDigit(atanh(tanh(0.5L)), 0.5, useDigits));
 1537 }
 1538 
 1539 /*****************************************
 1540  * Returns x rounded to a long value using the current rounding mode.
 1541  * If the integer value of x is
 1542  * greater than long.max, the result is
 1543  * indeterminate.
 1544  */
 1545 long rndtol(real x) @nogc @safe pure nothrow { pragma(inline, true); return core.math.rndtol(x); }
 1546 //FIXME
 1547 ///ditto
 1548 long rndtol(double x) @safe pure nothrow @nogc { return rndtol(cast(real) x); }
 1549 //FIXME
 1550 ///ditto
 1551 long rndtol(float x) @safe pure nothrow @nogc { return rndtol(cast(real) x); }
 1552 
 1553 ///
 1554 @safe unittest
 1555 {
 1556     assert(rndtol(1.0) == 1L);
 1557     assert(rndtol(1.2) == 1L);
 1558     assert(rndtol(1.7) == 2L);
 1559     assert(rndtol(1.0001) == 1L);
 1560 }
 1561 
 1562 @safe unittest
 1563 {
 1564     long function(real) prndtol = &rndtol;
 1565     assert(prndtol != null);
 1566 }
 1567 
 1568 /**
 1569 $(RED Deprecated. Please use $(LREF round) instead.)
 1570 
 1571 Returns `x` rounded to a `long` value using the `FE_TONEAREST` rounding mode.
 1572 If the integer value of `x` is greater than `long.max`, the result is
 1573 indeterminate.
 1574 
 1575 Only works with the Digital Mars C Runtime.
 1576 
 1577 Params:
 1578     x = the number to round
 1579 Returns:
 1580     `x` rounded to an integer value
 1581  */
 1582 deprecated("rndtonl is to be removed by 2.089. Please use round instead")
 1583 extern (C) real rndtonl(real x);
 1584 
 1585 ///
 1586 deprecated @system unittest
 1587 {
 1588     version (CRuntime_DigitalMars)
 1589     {
 1590         assert(rndtonl(1.0) is -real.nan);
 1591         assert(rndtonl(1.2) is -real.nan);
 1592         assert(rndtonl(1.7) is -real.nan);
 1593         assert(rndtonl(1.0001) is -real.nan);
 1594     }
 1595 }
 1596 
 1597 /***************************************
 1598  * Compute square root of x.
 1599  *
 1600  *      $(TABLE_SV
 1601  *      $(TR $(TH x)         $(TH sqrt(x))   $(TH invalid?))
 1602  *      $(TR $(TD -0.0)      $(TD -0.0)      $(TD no))
 1603  *      $(TR $(TD $(LT)0.0)  $(TD $(NAN))    $(TD yes))
 1604  *      $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no))
 1605  *      )
 1606  */
 1607 float sqrt(float x) @nogc @safe pure nothrow { pragma(inline, true); return core.math.sqrt(x); }
 1608 
 1609 /// ditto
 1610 double sqrt(double x) @nogc @safe pure nothrow { pragma(inline, true); return core.math.sqrt(x); }
 1611 
 1612 /// ditto
 1613 real sqrt(real x) @nogc @safe pure nothrow { pragma(inline, true); return core.math.sqrt(x); }
 1614 
 1615 ///
 1616 @safe pure nothrow @nogc unittest
 1617 {
 1618     assert(sqrt(2.0).feqrel(1.4142) > 16);
 1619     assert(sqrt(9.0).feqrel(3.0) > 16);
 1620 
 1621     assert(isNaN(sqrt(-1.0f)));
 1622     assert(isNaN(sqrt(-1.0)));
 1623     assert(isNaN(sqrt(-1.0L)));
 1624 }
 1625 
 1626 @safe unittest
 1627 {
 1628     float function(float) psqrtf = &sqrt;
 1629     assert(psqrtf != null);
 1630     double function(double) psqrtd = &sqrt;
 1631     assert(psqrtd != null);
 1632     real function(real) psqrtr = &sqrt;
 1633     assert(psqrtr != null);
 1634 
 1635     //ctfe
 1636     enum ZX80 = sqrt(7.0f);
 1637     enum ZX81 = sqrt(7.0);
 1638     enum ZX82 = sqrt(7.0L);
 1639 }
 1640 
 1641 /**
 1642  * Calculates e$(SUPERSCRIPT x).
 1643  *
 1644  *  $(TABLE_SV
 1645  *    $(TR $(TH x)             $(TH e$(SUPERSCRIPT x)) )
 1646  *    $(TR $(TD +$(INFIN))     $(TD +$(INFIN)) )
 1647  *    $(TR $(TD -$(INFIN))     $(TD +0.0)      )
 1648  *    $(TR $(TD $(NAN))        $(TD $(NAN))    )
 1649  *  )
 1650  */
 1651 real exp(real x) @trusted pure nothrow @nogc // TODO: @safe
 1652 {
 1653     version (D_InlineAsm_X86)
 1654     {
 1655         //  e^^x = 2^^(LOG2E*x)
 1656         // (This is valid because the overflow & underflow limits for exp
 1657         // and exp2 are so similar).
 1658         if (!__ctfe)
 1659             return exp2Asm(LOG2E*x);
 1660     }
 1661     else version (D_InlineAsm_X86_64)
 1662     {
 1663         //  e^^x = 2^^(LOG2E*x)
 1664         // (This is valid because the overflow & underflow limits for exp
 1665         // and exp2 are so similar).
 1666         if (!__ctfe)
 1667             return exp2Asm(LOG2E*x);
 1668     }
 1669     return expImpl(x);
 1670 }
 1671 
 1672 /// ditto
 1673 double exp(double x) @safe pure nothrow @nogc { return __ctfe ? cast(double) exp(cast(real) x) : expImpl(x); }
 1674 
 1675 /// ditto
 1676 float exp(float x) @safe pure nothrow @nogc { return __ctfe ? cast(float) exp(cast(real) x) : expImpl(x); }
 1677 
 1678 ///
 1679 @safe unittest
 1680 {
 1681     assert(exp(0.0) == 1.0);
 1682     assert(exp(3.0).feqrel(E * E * E) > 16);
 1683 }
 1684 
 1685 private T expImpl(T)(T x) @safe pure nothrow @nogc
 1686 {
 1687     alias F = floatTraits!T;
 1688     static if (F.realFormat == RealFormat.ieeeSingle)
 1689     {
 1690         static immutable T[6] P = [
 1691             5.0000001201E-1,
 1692             1.6666665459E-1,
 1693             4.1665795894E-2,
 1694             8.3334519073E-3,
 1695             1.3981999507E-3,
 1696             1.9875691500E-4,
 1697         ];
 1698 
 1699         enum T C1 = 0.693359375;
 1700         enum T C2 = -2.12194440e-4;
 1701 
 1702         // Overflow and Underflow limits.
 1703         enum T OF = 88.72283905206835;
 1704         enum T UF = -103.278929903431851103; // ln(2^-149)
 1705     }
 1706     else static if (F.realFormat == RealFormat.ieeeDouble)
 1707     {
 1708         // Coefficients for exp(x)
 1709         static immutable T[3] P = [
 1710             9.99999999999999999910E-1L,
 1711             3.02994407707441961300E-2L,
 1712             1.26177193074810590878E-4L,
 1713         ];
 1714         static immutable T[4] Q = [
 1715             2.00000000000000000009E0L,
 1716             2.27265548208155028766E-1L,
 1717             2.52448340349684104192E-3L,
 1718             3.00198505138664455042E-6L,
 1719         ];
 1720 
 1721         // C1 + C2 = LN2.
 1722         enum T C1 = 6.93145751953125E-1;
 1723         enum T C2 = 1.42860682030941723212E-6;
 1724 
 1725         // Overflow and Underflow limits.
 1726         enum T OF =  7.09782712893383996732E2;  // ln((1-2^-53) * 2^1024)
 1727         enum T UF = -7.451332191019412076235E2; // ln(2^-1075)
 1728     }
 1729     else static if (F.realFormat == RealFormat.ieeeExtended)
 1730     {
 1731         // Coefficients for exp(x)
 1732         static immutable T[3] P = [
 1733             9.9999999999999999991025E-1L,
 1734             3.0299440770744196129956E-2L,
 1735             1.2617719307481059087798E-4L,
 1736         ];
 1737         static immutable T[4] Q = [
 1738             2.0000000000000000000897E0L,
 1739             2.2726554820815502876593E-1L,
 1740             2.5244834034968410419224E-3L,
 1741             3.0019850513866445504159E-6L,
 1742         ];
 1743 
 1744         // C1 + C2 = LN2.
 1745         enum T C1 = 6.9314575195312500000000E-1L;
 1746         enum T C2 = 1.4286068203094172321215E-6L;
 1747 
 1748         // Overflow and Underflow limits.
 1749         enum T OF =  1.1356523406294143949492E4L;  // ln((1-2^-64) * 2^16384)
 1750         enum T UF = -1.13994985314888605586758E4L; // ln(2^-16446)
 1751     }
 1752     else static if (F.realFormat == RealFormat.ieeeQuadruple)
 1753     {
 1754         // Coefficients for exp(x) - 1
 1755         static immutable T[5] P = [
 1756             9.999999999999999999999999999999999998502E-1L,
 1757             3.508710990737834361215404761139478627390E-2L,
 1758             2.708775201978218837374512615596512792224E-4L,
 1759             6.141506007208645008909088812338454698548E-7L,
 1760             3.279723985560247033712687707263393506266E-10L
 1761         ];
 1762         static immutable T[6] Q = [
 1763             2.000000000000000000000000000000000000150E0,
 1764             2.368408864814233538909747618894558968880E-1L,
 1765             3.611828913847589925056132680618007270344E-3L,
 1766             1.504792651814944826817779302637284053660E-5L,
 1767             1.771372078166251484503904874657985291164E-8L,
 1768             2.980756652081995192255342779918052538681E-12L
 1769         ];
 1770 
 1771         // C1 + C2 = LN2.
 1772         enum T C1 = 6.93145751953125E-1L;
 1773         enum T C2 = 1.428606820309417232121458176568075500134E-6L;
 1774 
 1775         // Overflow and Underflow limits.
 1776         enum T OF =  1.135583025911358400418251384584930671458833e4L;
 1777         enum T UF = -1.143276959615573793352782661133116431383730e4L;
 1778     }
 1779     else
 1780         static assert(0, "Not implemented for this architecture");
 1781 
 1782     // Special cases.
 1783     if (isNaN(x))
 1784         return x;
 1785     if (x > OF)
 1786         return real.infinity;
 1787     if (x < UF)
 1788         return 0.0;
 1789 
 1790     // Express: e^^x = e^^g * 2^^n
 1791     //   = e^^g * e^^(n * LOG2E)
 1792     //   = e^^(g + n * LOG2E)
 1793     T xx = floor((cast(T) LOG2E) * x + cast(T) 0.5);
 1794     const int n = cast(int) xx;
 1795     x -= xx * C1;
 1796     x -= xx * C2;
 1797 
 1798     static if (F.realFormat == RealFormat.ieeeSingle)
 1799     {
 1800         xx = x * x;
 1801         x = poly(x, P) * xx + x + 1.0f;
 1802     }
 1803     else
 1804     {
 1805         // Rational approximation for exponential of the fractional part:
 1806         //  e^^x = 1 + 2x P(x^^2) / (Q(x^^2) - P(x^^2))
 1807         xx = x * x;
 1808         const T px = x * poly(xx, P);
 1809         x = px / (poly(xx, Q) - px);
 1810         x = (cast(T) 1.0) + (cast(T) 2.0) * x;
 1811     }
 1812 
 1813     // Scale by power of 2.
 1814     x = ldexp(x, n);
 1815 
 1816     return x;
 1817 }
 1818 
 1819 @safe @nogc nothrow unittest
 1820 {
 1821     version (FloatingPointControlSupport)
 1822     {
 1823         FloatingPointControl ctrl;
 1824         if (FloatingPointControl.hasExceptionTraps)
 1825             ctrl.disableExceptions(FloatingPointControl.allExceptions);
 1826         ctrl.rounding = FloatingPointControl.roundToNearest;
 1827     }
 1828 
 1829     static void testExp(T)()
 1830     {
 1831         enum realFormat = floatTraits!T.realFormat;
 1832         static if (realFormat == RealFormat.ieeeQuadruple)
 1833         {
 1834             static immutable T[2][] exptestpoints =
 1835             [ //  x               exp(x)
 1836                 [ 1.0L,           E                                        ],
 1837                 [ 0.5L,           0x1.a61298e1e069bc972dfefab6df34p+0L     ],
 1838                 [ 3.0L,           E*E*E                                    ],
 1839                 [ 0x1.6p+13L,     0x1.6e509d45728655cdb4840542acb5p+16250L ], // near overflow
 1840                 [ 0x1.7p+13L,     T.infinity                               ], // close overflow
 1841                 [ 0x1p+80L,       T.infinity                               ], // far overflow
 1842                 [ T.infinity,     T.infinity                               ],
 1843                 [-0x1.18p+13L,    0x1.5e4bf54b4807034ea97fef0059a6p-12927L ], // near underflow
 1844                 [-0x1.625p+13L,   0x1.a6bd68a39d11fec3a250cd97f524p-16358L ], // ditto
 1845                 [-0x1.62dafp+13L, 0x0.cb629e9813b80ed4d639e875be6cp-16382L ], // near underflow - subnormal
 1846                 [-0x1.6549p+13L,  0x0.0000000000000000000000000001p-16382L ], // ditto
 1847                 [-0x1.655p+13L,   0                                        ], // close underflow
 1848                 [-0x1p+30L,       0                                        ], // far underflow
 1849             ];
 1850         }
 1851         else static if (realFormat == RealFormat.ieeeExtended)
 1852         {
 1853             static immutable T[2][] exptestpoints =
 1854             [ //  x               exp(x)
 1855                 [ 1.0L,           E                            ],
 1856                 [ 0.5L,           0x1.a61298e1e069bc97p+0L     ],
 1857                 [ 3.0L,           E*E*E                        ],
 1858                 [ 0x1.1p+13L,     0x1.29aeffefc8ec645p+12557L  ], // near overflow
 1859                 [ 0x1.7p+13L,     T.infinity                   ], // close overflow
 1860                 [ 0x1p+80L,       T.infinity                   ], // far overflow
 1861                 [ T.infinity,     T.infinity                   ],
 1862                 [-0x1.18p+13L,    0x1.5e4bf54b4806db9p-12927L  ], // near underflow
 1863                 [-0x1.625p+13L,   0x1.a6bd68a39d11f35cp-16358L ], // ditto
 1864                 [-0x1.62dafp+13L, 0x1.96c53d30277021dp-16383L  ], // near underflow - subnormal
 1865                 [-0x1.643p+13L,   0x1p-16444L                  ], // ditto
 1866                 [-0x1.645p+13L,   0                            ], // close underflow
 1867                 [-0x1p+30L,       0                            ], // far underflow
 1868             ];
 1869         }
 1870         else static if (realFormat == RealFormat.ieeeDouble)
 1871         {
 1872             static immutable T[2][] exptestpoints =
 1873             [ //  x,             exp(x)
 1874                 [ 1.0L,          E                        ],
 1875                 [ 0.5L,          0x1.a61298e1e069cp+0L    ],
 1876                 [ 3.0L,          E*E*E                    ],
 1877                 [ 0x1.6p+9L,     0x1.93bf4ec282efbp+1015L ], // near overflow
 1878                 [ 0x1.7p+9L,     T.infinity               ], // close overflow
 1879                 [ 0x1p+80L,      T.infinity               ], // far overflow
 1880                 [ T.infinity,    T.infinity               ],
 1881                 [-0x1.6p+9L,     0x1.44a3824e5285fp-1016L ], // near underflow
 1882                 [-0x1.64p+9L,    0x0.06f84920bb2d4p-1022L ], // near underflow - subnormal
 1883                 [-0x1.743p+9L,   0x0.0000000000001p-1022L ], // ditto
 1884                 [-0x1.8p+9L,     0                        ], // close underflow
 1885                 [-0x1p+30L,      0                        ], // far underflow
 1886             ];
 1887         }
 1888         else static if (realFormat == RealFormat.ieeeSingle)
 1889         {
 1890             static immutable T[2][] exptestpoints =
 1891             [ //  x,             exp(x)
 1892                 [ 1.0L,          E                ],
 1893                 [ 0.5L,          0x1.a61299p+0L   ],
 1894                 [ 3.0L,          E*E*E            ],
 1895                 [ 0x1.62p+6L,    0x1.99b988p+127L ], // near overflow
 1896                 [ 0x1.7p+6L,     T.infinity       ], // close overflow
 1897                 [ 0x1p+80L,      T.infinity       ], // far overflow
 1898                 [ T.infinity,    T.infinity       ],
 1899                 [-0x1.5cp+6L,    0x1.666d0ep-126L ], // near underflow
 1900                 [-0x1.7p+6L,     0x0.026a42p-126L ], // near underflow - subnormal
 1901                 [-0x1.9cp+6L,    0x0.000002p-126L ], // ditto
 1902                 [-0x1.ap+6L,     0                ], // close underflow
 1903                 [-0x1p+30L,      0                ], // far underflow
 1904             ];
 1905         }
 1906         else
 1907             static assert(0, "No exp() tests for real type!");
 1908 
 1909         const minEqualMantissaBits = T.mant_dig - 2;
 1910         T x;
 1911         version (IeeeFlagsSupport) IeeeFlags f;
 1912         foreach (ref pair; exptestpoints)
 1913         {
 1914             version (IeeeFlagsSupport) resetIeeeFlags();
 1915             x = exp(pair[0]);
 1916             //printf("exp(%La) = %La, should be %La\n", cast(real) pair[0], cast(real) x, cast(real) pair[1]);
 1917             assert(feqrel(x, pair[1]) >= minEqualMantissaBits);
 1918         }
 1919 
 1920         // Ideally, exp(0) would not set the inexact flag.
 1921         // Unfortunately, fldl2e sets it!
 1922         // So it's not realistic to avoid setting it.
 1923         assert(exp(cast(T) 0.0) == 1.0);
 1924 
 1925         // NaN propagation. Doesn't set flags, bcos was already NaN.
 1926         version (IeeeFlagsSupport)
 1927         {
 1928             resetIeeeFlags();
 1929             x = exp(T.nan);
 1930             f = ieeeFlags;
 1931             assert(isIdentical(abs(x), T.nan));
 1932             assert(f.flags == 0);
 1933 
 1934             resetIeeeFlags();
 1935             x = exp(-T.nan);
 1936             f = ieeeFlags;
 1937             assert(isIdentical(abs(x), T.nan));
 1938             assert(f.flags == 0);
 1939         }
 1940         else
 1941         {
 1942             x = exp(T.nan);
 1943             assert(isIdentical(abs(x), T.nan));
 1944 
 1945             x = exp(-T.nan);
 1946             assert(isIdentical(abs(x), T.nan));
 1947         }
 1948 
 1949         x = exp(NaN(0x123));
 1950         assert(isIdentical(x, NaN(0x123)));
 1951     }
 1952 
 1953     import std.meta : AliasSeq;
 1954     foreach (T; AliasSeq!(real, double, float))
 1955         testExp!T();
 1956 
 1957     // High resolution test (verified against GNU MPFR/Mathematica).
 1958     assert(exp(0.5L) == 0x1.A612_98E1_E069_BC97_2DFE_FAB6_DF34p+0L);
 1959 
 1960     assert(equalsDigit(exp(3.0L), E * E * E, useDigits));
 1961 }
 1962 
 1963 /**
 1964  * Calculates the value of the natural logarithm base (e)
 1965  * raised to the power of x, minus 1.
 1966  *
 1967  * For very small x, expm1(x) is more accurate
 1968  * than exp(x)-1.
 1969  *
 1970  *  $(TABLE_SV
 1971  *    $(TR $(TH x)             $(TH e$(SUPERSCRIPT x)-1)  )
 1972  *    $(TR $(TD $(PLUSMN)0.0)  $(TD $(PLUSMN)0.0) )
 1973  *    $(TR $(TD +$(INFIN))     $(TD +$(INFIN))    )
 1974  *    $(TR $(TD -$(INFIN))     $(TD -1.0)         )
 1975  *    $(TR $(TD $(NAN))        $(TD $(NAN))       )
 1976  *  )
 1977  */
 1978 real expm1(real x) @trusted pure nothrow @nogc // TODO: @safe
 1979 {
 1980     version (InlineAsm_X86_Any)
 1981     {
 1982         if (!__ctfe)
 1983             return expm1Asm(x);
 1984     }
 1985     return expm1Impl(x);
 1986 }
 1987 
 1988 /// ditto
 1989 double expm1(double x) @safe pure nothrow @nogc
 1990 {
 1991     return __ctfe ? cast(double) expm1(cast(real) x) : expm1Impl(x);
 1992 }
 1993 
 1994 /// ditto
 1995 float expm1(float x) @safe pure nothrow @nogc
 1996 {
 1997     // no single-precision version in Cephes => use double precision
 1998     return __ctfe ? cast(float) expm1(cast(real) x) : cast(float) expm1Impl(cast(double) x);
 1999 }
 2000 
 2001 ///
 2002 @safe unittest
 2003 {
 2004     assert(isIdentical(expm1(0.0), 0.0));
 2005     assert(expm1(1.0).feqrel(1.71828) > 16);
 2006     assert(expm1(2.0).feqrel(6.3890) > 16);
 2007 }
 2008 
 2009 version (InlineAsm_X86_Any)
 2010 private real expm1Asm(real x) @trusted pure nothrow @nogc
 2011 {
 2012     version (D_InlineAsm_X86)
 2013     {
 2014         enum PARAMSIZE = (real.sizeof+3)&(0xFFFF_FFFC); // always a multiple of 4
 2015         asm pure nothrow @nogc
 2016         {
 2017             /*  expm1() for x87 80-bit reals, IEEE754-2008 conformant.
 2018              * Author: Don Clugston.
 2019              *
 2020              *    expm1(x) = 2^^(rndint(y))* 2^^(y-rndint(y)) - 1 where y = LN2*x.
 2021              *    = 2rndy * 2ym1 + 2rndy - 1, where 2rndy = 2^^(rndint(y))
 2022              *     and 2ym1 = (2^^(y-rndint(y))-1).
 2023              *    If 2rndy  < 0.5*real.epsilon, result is -1.
 2024              *    Implementation is otherwise the same as for exp2()
 2025              */
 2026             naked;
 2027             fld real ptr [ESP+4] ; // x
 2028             mov AX, [ESP+4+8]; // AX = exponent and sign
 2029             sub ESP, 12+8; // Create scratch space on the stack
 2030             // [ESP,ESP+2] = scratchint
 2031             // [ESP+4..+6, +8..+10, +10] = scratchreal
 2032             // set scratchreal mantissa = 1.0
 2033             mov dword ptr [ESP+8], 0;
 2034             mov dword ptr [ESP+8+4], 0x80000000;
 2035             and AX, 0x7FFF; // drop sign bit
 2036             cmp AX, 0x401D; // avoid InvalidException in fist
 2037             jae L_extreme;
 2038             fldl2e;
 2039             fmulp ST(1), ST; // y = x*log2(e)
 2040             fist dword ptr [ESP]; // scratchint = rndint(y)
 2041             fisub dword ptr [ESP]; // y - rndint(y)
 2042             // and now set scratchreal exponent
 2043             mov EAX, [ESP];
 2044             add EAX, 0x3fff;
 2045             jle short L_largenegative;
 2046             cmp EAX,0x8000;
 2047             jge short L_largepositive;
 2048             mov [ESP+8+8],AX;
 2049             f2xm1; // 2ym1 = 2^^(y-rndint(y)) -1
 2050             fld real ptr [ESP+8] ; // 2rndy = 2^^rndint(y)
 2051             fmul ST(1), ST;  // ST=2rndy, ST(1)=2rndy*2ym1
 2052             fld1;
 2053             fsubp ST(1), ST; // ST = 2rndy-1, ST(1) = 2rndy * 2ym1 - 1
 2054             faddp ST(1), ST; // ST = 2rndy * 2ym1 + 2rndy - 1
 2055             add ESP,12+8;
 2056             ret PARAMSIZE;
 2057 
 2058 L_extreme:  // Extreme exponent. X is very large positive, very
 2059             // large negative, infinity, or NaN.
 2060             fxam;
 2061             fstsw AX;
 2062             test AX, 0x0400; // NaN_or_zero, but we already know x != 0
 2063             jz L_was_nan;  // if x is NaN, returns x
 2064             test AX, 0x0200;
 2065             jnz L_largenegative;
 2066 L_largepositive:
 2067             // Set scratchreal = real.max.
 2068             // squaring it will create infinity, and set overflow flag.
 2069             mov word  ptr [ESP+8+8], 0x7FFE;
 2070             fstp ST(0);
 2071             fld real ptr [ESP+8];  // load scratchreal
 2072             fmul ST(0), ST;        // square it, to create havoc!
 2073 L_was_nan:
 2074             add ESP,12+8;
 2075             ret PARAMSIZE;
 2076 L_largenegative:
 2077             fstp ST(0);
 2078             fld1;
 2079             fchs; // return -1. Underflow flag is not set.
 2080             add ESP,12+8;
 2081             ret PARAMSIZE;
 2082         }
 2083     }
 2084     else version (D_InlineAsm_X86_64)
 2085     {
 2086         asm pure nothrow @nogc
 2087         {
 2088             naked;
 2089         }
 2090         version (Win64)
 2091         {
 2092             asm pure nothrow @nogc
 2093             {
 2094                 fld   real ptr [RCX];  // x
 2095                 mov   AX,[RCX+8];      // AX = exponent and sign
 2096             }
 2097         }
 2098         else
 2099         {
 2100             asm pure nothrow @nogc
 2101             {
 2102                 fld   real ptr [RSP+8];  // x
 2103                 mov   AX,[RSP+8+8];      // AX = exponent and sign
 2104             }
 2105         }
 2106         asm pure nothrow @nogc
 2107         {
 2108             /*  expm1() for x87 80-bit reals, IEEE754-2008 conformant.
 2109              * Author: Don Clugston.
 2110              *
 2111              *    expm1(x) = 2^(rndint(y))* 2^(y-rndint(y)) - 1 where y = LN2*x.
 2112              *    = 2rndy * 2ym1 + 2rndy - 1, where 2rndy = 2^(rndint(y))
 2113              *     and 2ym1 = (2^(y-rndint(y))-1).
 2114              *    If 2rndy  < 0.5*real.epsilon, result is -1.
 2115              *    Implementation is otherwise the same as for exp2()
 2116              */
 2117             sub RSP, 24;       // Create scratch space on the stack
 2118             // [RSP,RSP+2] = scratchint
 2119             // [RSP+4..+6, +8..+10, +10] = scratchreal
 2120             // set scratchreal mantissa = 1.0
 2121             mov dword ptr [RSP+8], 0;
 2122             mov dword ptr [RSP+8+4], 0x80000000;
 2123             and AX, 0x7FFF; // drop sign bit
 2124             cmp AX, 0x401D; // avoid InvalidException in fist
 2125             jae L_extreme;
 2126             fldl2e;
 2127             fmul ; // y = x*log2(e)
 2128             fist dword ptr [RSP]; // scratchint = rndint(y)
 2129             fisub dword ptr [RSP]; // y - rndint(y)
 2130             // and now set scratchreal exponent
 2131             mov EAX, [RSP];
 2132             add EAX, 0x3fff;
 2133             jle short L_largenegative;
 2134             cmp EAX,0x8000;
 2135             jge short L_largepositive;
 2136             mov [RSP+8+8],AX;
 2137             f2xm1; // 2^(y-rndint(y)) -1
 2138             fld real ptr [RSP+8] ; // 2^rndint(y)
 2139             fmul ST(1), ST;
 2140             fld1;
 2141             fsubp ST(1), ST;
 2142             fadd;
 2143             add RSP,24;
 2144             ret;
 2145 
 2146 L_extreme:  // Extreme exponent. X is very large positive, very
 2147             // large negative, infinity, or NaN.
 2148             fxam;
 2149             fstsw AX;
 2150             test AX, 0x0400; // NaN_or_zero, but we already know x != 0
 2151             jz L_was_nan;  // if x is NaN, returns x
 2152             test AX, 0x0200;
 2153             jnz L_largenegative;
 2154 L_largepositive:
 2155             // Set scratchreal = real.max.
 2156             // squaring it will create infinity, and set overflow flag.
 2157             mov word  ptr [RSP+8+8], 0x7FFE;
 2158             fstp ST(0);
 2159             fld real ptr [RSP+8];  // load scratchreal
 2160             fmul ST(0), ST;        // square it, to create havoc!
 2161 L_was_nan:
 2162             add RSP,24;
 2163             ret;
 2164 
 2165 L_largenegative:
 2166             fstp ST(0);
 2167             fld1;
 2168             fchs; // return -1. Underflow flag is not set.
 2169             add RSP,24;
 2170             ret;
 2171         }
 2172     }
 2173     else
 2174         static assert(0);
 2175 }
 2176 
 2177 private T expm1Impl(T)(T x) @safe pure nothrow @nogc
 2178 {
 2179     // Coefficients for exp(x) - 1 and overflow/underflow limits.
 2180     enum realFormat = floatTraits!T.realFormat;
 2181     static if (realFormat == RealFormat.ieeeQuadruple)
 2182     {
 2183         static immutable T[8] P = [
 2184             2.943520915569954073888921213330863757240E8L,
 2185             -5.722847283900608941516165725053359168840E7L,
 2186             8.944630806357575461578107295909719817253E6L,
 2187             -7.212432713558031519943281748462837065308E5L,
 2188             4.578962475841642634225390068461943438441E4L,
 2189             -1.716772506388927649032068540558788106762E3L,
 2190             4.401308817383362136048032038528753151144E1L,
 2191             -4.888737542888633647784737721812546636240E-1L
 2192         ];
 2193 
 2194         static immutable T[9] Q = [
 2195             1.766112549341972444333352727998584753865E9L,
 2196             -7.848989743695296475743081255027098295771E8L,
 2197             1.615869009634292424463780387327037251069E8L,
 2198             -2.019684072836541751428967854947019415698E7L,
 2199             1.682912729190313538934190635536631941751E6L,
 2200             -9.615511549171441430850103489315371768998E4L,
 2201             3.697714952261803935521187272204485251835E3L,
 2202             -8.802340681794263968892934703309274564037E1L,
 2203             1.0
 2204         ];
 2205 
 2206         enum T OF = 1.1356523406294143949491931077970764891253E4L;
 2207         enum T UF = -1.143276959615573793352782661133116431383730e4L;
 2208     }
 2209     else static if (realFormat == RealFormat.ieeeExtended)
 2210     {
 2211         static immutable T[5] P = [
 2212            -1.586135578666346600772998894928250240826E4L,
 2213             2.642771505685952966904660652518429479531E3L,
 2214            -3.423199068835684263987132888286791620673E2L,
 2215             1.800826371455042224581246202420972737840E1L,
 2216            -5.238523121205561042771939008061958820811E-1L,
 2217         ];
 2218         static immutable T[6] Q = [
 2219            -9.516813471998079611319047060563358064497E4L,
 2220             3.964866271411091674556850458227710004570E4L,
 2221            -7.207678383830091850230366618190187434796E3L,
 2222             7.206038318724600171970199625081491823079E2L,
 2223            -4.002027679107076077238836622982900945173E1L,
 2224             1.0
 2225         ];
 2226 
 2227         enum T OF =  1.1356523406294143949492E4L;
 2228         enum T UF = -4.5054566736396445112120088E1L;
 2229     }
 2230     else static if (realFormat == RealFormat.ieeeDouble)
 2231     {
 2232         static immutable T[3] P = [
 2233             9.9999999999999999991025E-1,
 2234             3.0299440770744196129956E-2,
 2235             1.2617719307481059087798E-4,
 2236         ];
 2237         static immutable T[4] Q = [
 2238             2.0000000000000000000897E0,
 2239             2.2726554820815502876593E-1,
 2240             2.5244834034968410419224E-3,
 2241             3.0019850513866445504159E-6,
 2242         ];
 2243     }
 2244     else
 2245         static assert(0, "no coefficients for expm1()");
 2246 
 2247     static if (realFormat == RealFormat.ieeeDouble) // special case for double precision
 2248     {
 2249         if (x < -0.5 || x > 0.5)
 2250             return exp(x) - 1.0;
 2251         if (x == 0.0)
 2252             return x;
 2253 
 2254         const T xx = x * x;
 2255         x = x * poly(xx, P);
 2256         x = x / (poly(xx, Q) - x);
 2257         return x + x;
 2258     }
 2259     else
 2260     {
 2261         // C1 + C2 = LN2.
 2262         enum T C1 = 6.9314575195312500000000E-1L;
 2263         enum T C2 = 1.428606820309417232121458176568075500134E-6L;
 2264 
 2265         // Special cases.
 2266         if (x > OF)
 2267             return real.infinity;
 2268         if (x == cast(T) 0.0)
 2269             return x;
 2270         if (x < UF)
 2271             return -1.0;
 2272 
 2273         // Express x = LN2 (n + remainder), remainder not exceeding 1/2.
 2274         int n = cast(int) floor((cast(T) 0.5) + x / cast(T) LN2);
 2275         x -= n * C1;
 2276         x -= n * C2;
 2277 
 2278         // Rational approximation:
 2279         //  exp(x) - 1 = x + 0.5 x^^2 + x^^3 P(x) / Q(x)
 2280         T px = x * poly(x, P);
 2281         T qx = poly(x, Q);
 2282         const T xx = x * x;
 2283         qx = x + ((cast(T) 0.5) * xx + xx * px / qx);
 2284 
 2285         // We have qx = exp(remainder LN2) - 1, so:
 2286         //  exp(x) - 1 = 2^^n (qx + 1) - 1 = 2^^n qx + 2^^n - 1.
 2287         px = ldexp(cast(T) 1.0, n);
 2288         x = px * qx + (px - cast(T) 1.0);
 2289 
 2290         return x;
 2291     }
 2292 }
 2293 
 2294 @safe @nogc nothrow unittest
 2295 {
 2296     static void testExpm1(T)()
 2297     {
 2298         // NaN
 2299         assert(isNaN(expm1(cast(T) T.nan)));
 2300 
 2301         static immutable T[] xs = [ -2, -0.75, -0.3, 0.0, 0.1, 0.2, 0.5, 1.0 ];
 2302         foreach (x; xs)
 2303         {
 2304             const T e = expm1(x);
 2305             const T r = exp(x) - 1;
 2306 
 2307             //printf("expm1(%Lg) = %Lg, should approximately be %Lg\n", cast(real) x, cast(real) e, cast(real) r);
 2308             assert(approxEqual(r, e));
 2309         }
 2310     }
 2311 
 2312     import std.meta : AliasSeq;
 2313     foreach (T; AliasSeq!(real, double))
 2314         testExpm1!T();
 2315 }
 2316 
 2317 /**
 2318  * Calculates 2$(SUPERSCRIPT x).
 2319  *
 2320  *  $(TABLE_SV
 2321  *    $(TR $(TH x)             $(TH exp2(x))   )
 2322  *    $(TR $(TD +$(INFIN))     $(TD +$(INFIN)) )
 2323  *    $(TR $(TD -$(INFIN))     $(TD +0.0)      )
 2324  *    $(TR $(TD $(NAN))        $(TD $(NAN))    )
 2325  *  )
 2326  */
 2327 real exp2(real x) @nogc @trusted pure nothrow // TODO: @safe
 2328 {
 2329     version (InlineAsm_X86_Any)
 2330     {
 2331         if (!__ctfe)
 2332             return exp2Asm(x);
 2333     }
 2334     return exp2Impl(x);
 2335 }
 2336 
 2337 /// ditto
 2338 double exp2(double x) @nogc @safe pure nothrow { return __ctfe ? cast(double) exp2(cast(real) x) : exp2Impl(x); }
 2339 
 2340 /// ditto
 2341 float exp2(float x) @nogc @safe pure nothrow { return __ctfe ? cast(float) exp2(cast(real) x) : exp2Impl(x); }
 2342 
 2343 ///
 2344 @safe unittest
 2345 {
 2346     assert(isIdentical(exp2(0.0), 1.0));
 2347     assert(exp2(2.0).feqrel(4.0) > 16);
 2348     assert(exp2(8.0).feqrel(256.0) > 16);
 2349 }
 2350 
 2351 @safe unittest
 2352 {
 2353     version (CRuntime_Microsoft) {} else // aexp2/exp2f/exp2l not implemented
 2354     {
 2355         assert( core.stdc.math.exp2f(0.0f) == 1 );
 2356         assert( core.stdc.math.exp2 (0.0)  == 1 );
 2357         assert( core.stdc.math.exp2l(0.0L) == 1 );
 2358     }
 2359 }
 2360 
 2361 version (InlineAsm_X86_Any)
 2362 private real exp2Asm(real x) @nogc @trusted pure nothrow
 2363 {
 2364     version (D_InlineAsm_X86)
 2365     {
 2366         enum PARAMSIZE = (real.sizeof+3)&(0xFFFF_FFFC); // always a multiple of 4
 2367 
 2368         asm pure nothrow @nogc
 2369         {
 2370             /*  exp2() for x87 80-bit reals, IEEE754-2008 conformant.
 2371              * Author: Don Clugston.
 2372              *
 2373              * exp2(x) = 2^^(rndint(x))* 2^^(y-rndint(x))
 2374              * The trick for high performance is to avoid the fscale(28cycles on core2),
 2375              * frndint(19 cycles), leaving f2xm1(19 cycles) as the only slow instruction.
 2376              *
 2377              * We can do frndint by using fist. BUT we can't use it for huge numbers,
 2378              * because it will set the Invalid Operation flag if overflow or NaN occurs.
 2379              * Fortunately, whenever this happens the result would be zero or infinity.
 2380              *
 2381              * We can perform fscale by directly poking into the exponent. BUT this doesn't
 2382              * work for the (very rare) cases where the result is subnormal. So we fall back
 2383              * to the slow method in that case.
 2384              */
 2385             naked;
 2386             fld real ptr [ESP+4] ; // x
 2387             mov AX, [ESP+4+8]; // AX = exponent and sign
 2388             sub ESP, 12+8; // Create scratch space on the stack
 2389             // [ESP,ESP+2] = scratchint
 2390             // [ESP+4..+6, +8..+10, +10] = scratchreal
 2391             // set scratchreal mantissa = 1.0
 2392             mov dword ptr [ESP+8], 0;
 2393             mov dword ptr [ESP+8+4], 0x80000000;
 2394             and AX, 0x7FFF; // drop sign bit
 2395             cmp AX, 0x401D; // avoid InvalidException in fist
 2396             jae L_extreme;
 2397             fist dword ptr [ESP]; // scratchint = rndint(x)
 2398             fisub dword ptr [ESP]; // x - rndint(x)
 2399             // and now set scratchreal exponent
 2400             mov EAX, [ESP];
 2401             add EAX, 0x3fff;
 2402             jle short L_subnormal;
 2403             cmp EAX,0x8000;
 2404             jge short L_overflow;
 2405             mov [ESP+8+8],AX;
 2406 L_normal:
 2407             f2xm1;
 2408             fld1;
 2409             faddp ST(1), ST; // 2^^(x-rndint(x))
 2410             fld real ptr [ESP+8] ; // 2^^rndint(x)
 2411             add ESP,12+8;
 2412             fmulp ST(1), ST;
 2413             ret PARAMSIZE;
 2414 
 2415 L_subnormal:
 2416             // Result will be subnormal.
 2417             // In this rare case, the simple poking method doesn't work.
 2418             // The speed doesn't matter, so use the slow fscale method.
 2419             fild dword ptr [ESP];  // scratchint
 2420             fld1;
 2421             fscale;
 2422             fstp real ptr [ESP+8]; // scratchreal = 2^^scratchint
 2423             fstp ST(0);         // drop scratchint
 2424             jmp L_normal;
 2425 
 2426 L_extreme:  // Extreme exponent. X is very large positive, very
 2427             // large negative, infinity, or NaN.
 2428             fxam;
 2429             fstsw AX;
 2430             test AX, 0x0400; // NaN_or_zero, but we already know x != 0
 2431             jz L_was_nan;  // if x is NaN, returns x
 2432             // set scratchreal = real.min_normal
 2433             // squaring it will return 0, setting underflow flag
 2434             mov word  ptr [ESP+8+8], 1;
 2435             test AX, 0x0200;
 2436             jnz L_waslargenegative;
 2437 L_overflow:
 2438             // Set scratchreal = real.max.
 2439             // squaring it will create infinity, and set overflow flag.
 2440             mov word  ptr [ESP+8+8], 0x7FFE;
 2441 L_waslargenegative:
 2442             fstp ST(0);
 2443             fld real ptr [ESP+8];  // load scratchreal
 2444             fmul ST(0), ST;        // square it, to create havoc!
 2445 L_was_nan:
 2446             add ESP,12+8;
 2447             ret PARAMSIZE;
 2448         }
 2449     }
 2450     else version (D_InlineAsm_X86_64)
 2451     {
 2452         asm pure nothrow @nogc
 2453         {
 2454             naked;
 2455         }
 2456         version (Win64)
 2457         {
 2458             asm pure nothrow @nogc
 2459             {
 2460                 fld   real ptr [RCX];  // x
 2461                 mov   AX,[RCX+8];      // AX = exponent and sign
 2462             }
 2463         }
 2464         else
 2465         {
 2466             asm pure nothrow @nogc
 2467             {
 2468                 fld   real ptr [RSP+8];  // x
 2469                 mov   AX,[RSP+8+8];      // AX = exponent and sign
 2470             }
 2471         }
 2472         asm pure nothrow @nogc
 2473         {
 2474             /*  exp2() for x87 80-bit reals, IEEE754-2008 conformant.
 2475              * Author: Don Clugston.
 2476              *
 2477              * exp2(x) = 2^(rndint(x))* 2^(y-rndint(x))
 2478              * The trick for high performance is to avoid the fscale(28cycles on core2),
 2479              * frndint(19 cycles), leaving f2xm1(19 cycles) as the only slow instruction.
 2480              *
 2481              * We can do frndint by using fist. BUT we can't use it for huge numbers,
 2482              * because it will set the Invalid Operation flag is overflow or NaN occurs.
 2483              * Fortunately, whenever this happens the result would be zero or infinity.
 2484              *
 2485              * We can perform fscale by directly poking into the exponent. BUT this doesn't
 2486              * work for the (very rare) cases where the result is subnormal. So we fall back
 2487              * to the slow method in that case.
 2488              */
 2489             sub RSP, 24; // Create scratch space on the stack
 2490             // [RSP,RSP+2] = scratchint
 2491             // [RSP+4..+6, +8..+10, +10] = scratchreal
 2492             // set scratchreal mantissa = 1.0
 2493             mov dword ptr [RSP+8], 0;
 2494             mov dword ptr [RSP+8+4], 0x80000000;
 2495             and AX, 0x7FFF; // drop sign bit
 2496             cmp AX, 0x401D; // avoid InvalidException in fist
 2497             jae L_extreme;
 2498             fist dword ptr [RSP]; // scratchint = rndint(x)
 2499             fisub dword ptr [RSP]; // x - rndint(x)
 2500             // and now set scratchreal exponent
 2501             mov EAX, [RSP];
 2502             add EAX, 0x3fff;
 2503             jle short L_subnormal;
 2504             cmp EAX,0x8000;
 2505             jge short L_overflow;
 2506             mov [RSP+8+8],AX;
 2507 L_normal:
 2508             f2xm1;
 2509             fld1;
 2510             fadd; // 2^(x-rndint(x))
 2511             fld real ptr [RSP+8] ; // 2^rndint(x)
 2512             add RSP,24;
 2513             fmulp ST(1), ST;
 2514             ret;
 2515 
 2516 L_subnormal:
 2517             // Result will be subnormal.
 2518             // In this rare case, the simple poking method doesn't work.
 2519             // The speed doesn't matter, so use the slow fscale method.
 2520             fild dword ptr [RSP];  // scratchint
 2521             fld1;
 2522             fscale;
 2523             fstp real ptr [RSP+8]; // scratchreal = 2^scratchint
 2524             fstp ST(0);         // drop scratchint
 2525             jmp L_normal;
 2526 
 2527 L_extreme:  // Extreme exponent. X is very large positive, very
 2528             // large negative, infinity, or NaN.
 2529             fxam;
 2530             fstsw AX;
 2531             test AX, 0x0400; // NaN_or_zero, but we already know x != 0
 2532             jz L_was_nan;  // if x is NaN, returns x
 2533             // set scratchreal = real.min
 2534             // squaring it will return 0, setting underflow flag
 2535             mov word  ptr [RSP+8+8], 1;
 2536             test AX, 0x0200;
 2537             jnz L_waslargenegative;
 2538 L_overflow:
 2539             // Set scratchreal = real.max.
 2540             // squaring it will create infinity, and set overflow flag.
 2541             mov word  ptr [RSP+8+8], 0x7FFE;
 2542 L_waslargenegative:
 2543             fstp ST(0);
 2544             fld real ptr [RSP+8];  // load scratchreal
 2545             fmul ST(0), ST;        // square it, to create havoc!
 2546 L_was_nan:
 2547             add RSP,24;
 2548             ret;
 2549         }
 2550     }
 2551     else
 2552         static assert(0);
 2553 }
 2554 
 2555 private T exp2Impl(T)(T x) @nogc @safe pure nothrow
 2556 {
 2557     // Coefficients for exp2(x)
 2558     enum realFormat = floatTraits!T.realFormat;
 2559     static if (realFormat == RealFormat.ieeeQuadruple)
 2560     {
 2561         static immutable T[5] P = [
 2562             9.079594442980146270952372234833529694788E12L,
 2563             1.530625323728429161131811299626419117557E11L,
 2564             5.677513871931844661829755443994214173883E8L,
 2565             6.185032670011643762127954396427045467506E5L,
 2566             1.587171580015525194694938306936721666031E2L
 2567         ];
 2568 
 2569         static immutable T[6] Q = [
 2570             2.619817175234089411411070339065679229869E13L,
 2571             1.490560994263653042761789432690793026977E12L,
 2572             1.092141473886177435056423606755843616331E10L,
 2573             2.186249607051644894762167991800811827835E7L,
 2574             1.236602014442099053716561665053645270207E4L,
 2575             1.0
 2576         ];
 2577     }
 2578     else static if (realFormat == RealFormat.ieeeExtended)
 2579     {
 2580         static immutable T[3] P = [
 2581             2.0803843631901852422887E6L,
 2582             3.0286971917562792508623E4L,
 2583             6.0614853552242266094567E1L,
 2584         ];
 2585         static immutable T[4] Q = [
 2586             6.0027204078348487957118E6L,
 2587             3.2772515434906797273099E5L,
 2588             1.7492876999891839021063E3L,
 2589             1.0000000000000000000000E0L,
 2590         ];
 2591     }
 2592     else static if (realFormat == RealFormat.ieeeDouble)
 2593     {
 2594         static immutable T[3] P = [
 2595             1.51390680115615096133E3L,
 2596             2.02020656693165307700E1L,
 2597             2.30933477057345225087E-2L,
 2598         ];
 2599         static immutable T[3] Q = [
 2600             4.36821166879210612817E3L,
 2601             2.33184211722314911771E2L,
 2602             1.00000000000000000000E0L,
 2603         ];
 2604     }
 2605     else static if (realFormat == RealFormat.ieeeSingle)
 2606     {
 2607         static immutable T[6] P = [
 2608             6.931472028550421E-001L,
 2609             2.402264791363012E-001L,
 2610             5.550332471162809E-002L,
 2611             9.618437357674640E-003L,
 2612             1.339887440266574E-003L,
 2613             1.535336188319500E-004L,
 2614         ];
 2615     }
 2616     else
 2617         static assert(0, "no coefficients for exp2()");
 2618 
 2619     // Overflow and Underflow limits.
 2620     enum T OF = T.max_exp;
 2621     enum T UF = T.min_exp - 1;
 2622 
 2623     // Special cases.
 2624     if (isNaN(x))
 2625         return x;
 2626     if (x > OF)
 2627         return real.infinity;
 2628     if (x < UF)
 2629         return 0.0;
 2630 
 2631     static if (realFormat == RealFormat.ieeeSingle) // special case for single precision
 2632     {
 2633         // The following is necessary because range reduction blows up.
 2634         if (x == 0.0f)
 2635             return 1.0f;
 2636 
 2637         // Separate into integer and fractional parts.
 2638         const T i = floor(x);
 2639         int n = cast(int) i;
 2640         x -= i;
 2641         if (x > 0.5f)
 2642         {
 2643             n += 1;
 2644             x -= 1.0f;
 2645         }
 2646 
 2647         // Rational approximation:
 2648         //  exp2(x) = 1.0 + x P(x)
 2649         x = 1.0f + x * poly(x, P);
 2650     }
 2651     else
 2652     {
 2653         // Separate into integer and fractional parts.
 2654         const T i = floor(x + cast(T) 0.5);
 2655         int n = cast(int) i;
 2656         x -= i;
 2657 
 2658         // Rational approximation:
 2659         //  exp2(x) = 1.0 + 2x P(x^^2) / (Q(x^^2) - P(x^^2))
 2660         const T xx = x * x;
 2661         const T px = x * poly(xx, P);
 2662         x = px / (poly(xx, Q) - px);
 2663         x = (cast(T) 1.0) + (cast(T) 2.0) * x;
 2664     }
 2665 
 2666     // Scale by power of 2.
 2667     x = ldexp(x, n);
 2668 
 2669     return x;
 2670 }
 2671 
 2672 @safe @nogc nothrow unittest
 2673 {
 2674     assert(feqrel(exp2(0.5L), SQRT2) >= real.mant_dig -1);
 2675     assert(exp2(8.0L) == 256.0);
 2676     assert(exp2(-9.0L)== 1.0L/512.0);
 2677 
 2678     static void testExp2(T)()
 2679     {
 2680         // NaN
 2681         const T specialNaN = NaN(0x0123L);
 2682         assert(isIdentical(exp2(specialNaN), specialNaN));
 2683 
 2684         // over-/underflow
 2685         enum T OF = T.max_exp;
 2686         enum T UF = T.min_exp - T.mant_dig;
 2687         assert(isIdentical(exp2(OF + 1), cast(T) T.infinity));
 2688         assert(isIdentical(exp2(UF - 1), cast(T) 0.0));
 2689 
 2690         static immutable T[2][] vals =
 2691         [
 2692             // x, exp2(x)
 2693             [  0.0, 1.0 ],
 2694             [ -0.0, 1.0 ],
 2695             [  0.5, SQRT2 ],
 2696             [  8.0, 256.0 ],
 2697             [ -9.0, 1.0 / 512 ],
 2698         ];
 2699 
 2700         foreach (ref val; vals)
 2701         {
 2702             const T x = val[0];
 2703             const T r = val[1];
 2704             const T e = exp2(x);
 2705 
 2706             //printf("exp2(%Lg) = %Lg, should be %Lg\n", cast(real) x, cast(real) e, cast(real) r);
 2707             assert(approxEqual(r, e));
 2708         }
 2709     }
 2710 
 2711     import std.meta : AliasSeq;
 2712     foreach (T; AliasSeq!(real, double, float))
 2713         testExp2!T();
 2714 }
 2715 
 2716 /*********************************************************************
 2717  * Separate floating point value into significand and exponent.
 2718  *
 2719  * Returns:
 2720  *      Calculate and return $(I x) and $(I exp) such that
 2721  *      value =$(I x)*2$(SUPERSCRIPT exp) and
 2722  *      .5 $(LT)= |$(I x)| $(LT) 1.0
 2723  *
 2724  *      $(I x) has same sign as value.
 2725  *
 2726  *      $(TABLE_SV
 2727  *      $(TR $(TH value)           $(TH returns)         $(TH exp))
 2728  *      $(TR $(TD $(PLUSMN)0.0)    $(TD $(PLUSMN)0.0)    $(TD 0))
 2729  *      $(TR $(TD +$(INFIN))       $(TD +$(INFIN))       $(TD int.max))
 2730  *      $(TR $(TD -$(INFIN))       $(TD -$(INFIN))       $(TD int.min))
 2731  *      $(TR $(TD $(PLUSMN)$(NAN)) $(TD $(PLUSMN)$(NAN)) $(TD int.min))
 2732  *      )
 2733  */
 2734 T frexp(T)(const T value, out int exp) @trusted pure nothrow @nogc
 2735 if (isFloatingPoint!T)
 2736 {
 2737     if (__ctfe)
 2738     {
 2739         // Handle special cases.
 2740         if (value == 0) { exp = 0; return value; }
 2741         if (value == T.infinity) { exp = int.max; return value; }
 2742         if (value == -T.infinity || value != value) { exp = int.min; return value; }
 2743         // Handle ordinary cases.
 2744         // In CTFE there is no performance advantage for having separate
 2745         // paths for different floating point types.
 2746         T absValue = value < 0 ? -value : value;
 2747         int expCount;
 2748         static if (T.mant_dig > double.mant_dig)
 2749         {
 2750             for (; absValue >= 0x1.0p+1024L; absValue *= 0x1.0p-1024L)
 2751                 expCount += 1024;
 2752             for (; absValue < 0x1.0p-1021L; absValue *= 0x1.0p+1021L)
 2753                 expCount -= 1021;
 2754         }
 2755         const double dval = cast(double) absValue;
 2756         int dexp = cast(int) (((*cast(const long*) &dval) >>> 52) & 0x7FF) + double.min_exp - 2;
 2757         dexp++;
 2758         expCount += dexp;
 2759         absValue *= 2.0 ^^ -dexp;
 2760         // If the original value was subnormal or if it was a real
 2761         // then absValue can still be outside the [0.5, 1.0) range.
 2762         if (absValue < 0.5)
 2763         {
 2764             assert(T.mant_dig > double.mant_dig || isSubnormal(value));
 2765             do
 2766             {
 2767                 absValue += absValue;
 2768                 expCount--;
 2769             } while (absValue < 0.5);
 2770         }
 2771         else
 2772         {
 2773             assert(absValue < 1 || T.mant_dig > double.mant_dig);
 2774             for (; absValue >= 1; absValue *= T(0.5))
 2775                 expCount++;
 2776         }
 2777         exp = expCount;
 2778         return value < 0 ? -absValue : absValue;
 2779     }
 2780 
 2781     Unqual!T vf = value;
 2782     ushort* vu = cast(ushort*)&vf;
 2783     static if (is(immutable T == immutable float))
 2784         int* vi = cast(int*)&vf;
 2785     else
 2786         long* vl = cast(long*)&vf;
 2787     int ex;
 2788     alias F = floatTraits!T;
 2789 
 2790     ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
 2791     static if (F.realFormat == RealFormat.ieeeExtended)
 2792     {
 2793         if (ex)
 2794         {   // If exponent is non-zero
 2795             if (ex == F.EXPMASK) // infinity or NaN
 2796             {
 2797                 if (*vl &  0x7FFF_FFFF_FFFF_FFFF)  // NaN
 2798                 {
 2799                     *vl |= 0xC000_0000_0000_0000;  // convert NaNS to NaNQ
 2800                     exp = int.min;
 2801                 }
 2802                 else if (vu[F.EXPPOS_SHORT] & 0x8000)   // negative infinity
 2803                     exp = int.min;
 2804                 else   // positive infinity
 2805                     exp = int.max;
 2806 
 2807             }
 2808             else
 2809             {
 2810                 exp = ex - F.EXPBIAS;
 2811                 vu[F.EXPPOS_SHORT] = (0x8000 & vu[F.EXPPOS_SHORT]) | 0x3FFE;
 2812             }
 2813         }
 2814         else if (!*vl)
 2815         {
 2816             // vf is +-0.0
 2817             exp = 0;
 2818         }
 2819         else
 2820         {
 2821             // subnormal
 2822 
 2823             vf *= F.RECIP_EPSILON;
 2824             ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
 2825             exp = ex - F.EXPBIAS - T.mant_dig + 1;
 2826             vu[F.EXPPOS_SHORT] = ((-1 - F.EXPMASK) & vu[F.EXPPOS_SHORT]) | 0x3FFE;
 2827         }
 2828         return vf;
 2829     }
 2830     else static if (F.realFormat == RealFormat.ieeeQuadruple)
 2831     {
 2832         if (ex)     // If exponent is non-zero
 2833         {
 2834             if (ex == F.EXPMASK)
 2835             {
 2836                 // infinity or NaN
 2837                 if (vl[MANTISSA_LSB] |
 2838                     (vl[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF))  // NaN
 2839                 {
 2840                     // convert NaNS to NaNQ
 2841                     vl[MANTISSA_MSB] |= 0x0000_8000_0000_0000;
 2842                     exp = int.min;
 2843                 }
 2844                 else if (vu[F.EXPPOS_SHORT] & 0x8000)   // negative infinity
 2845                     exp = int.min;
 2846                 else   // positive infinity
 2847                     exp = int.max;
 2848             }
 2849             else
 2850             {
 2851                 exp = ex - F.EXPBIAS;
 2852                 vu[F.EXPPOS_SHORT] = F.EXPBIAS | (0x8000 & vu[F.EXPPOS_SHORT]);
 2853             }
 2854         }
 2855         else if ((vl[MANTISSA_LSB] |
 2856             (vl[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) == 0)
 2857         {
 2858             // vf is +-0.0
 2859             exp = 0;
 2860         }
 2861         else
 2862         {
 2863             // subnormal
 2864             vf *= F.RECIP_EPSILON;
 2865             ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
 2866             exp = ex - F.EXPBIAS - T.mant_dig + 1;
 2867             vu[F.EXPPOS_SHORT] = F.EXPBIAS | (0x8000 & vu[F.EXPPOS_SHORT]);
 2868         }
 2869         return vf;
 2870     }
 2871     else static if (F.realFormat == RealFormat.ieeeDouble)
 2872     {
 2873         if (ex) // If exponent is non-zero
 2874         {
 2875             if (ex == F.EXPMASK)   // infinity or NaN
 2876             {
 2877                 if (*vl == 0x7FF0_0000_0000_0000)  // positive infinity
 2878                 {
 2879                     exp = int.max;
 2880                 }
 2881                 else if (*vl == 0xFFF0_0000_0000_0000) // negative infinity
 2882                     exp = int.min;
 2883                 else
 2884                 { // NaN
 2885                     *vl |= 0x0008_0000_0000_0000;  // convert NaNS to NaNQ
 2886                     exp = int.min;
 2887                 }
 2888             }
 2889             else
 2890             {
 2891                 exp = (ex - F.EXPBIAS) >> 4;
 2892                 vu[F.EXPPOS_SHORT] = cast(ushort)((0x800F & vu[F.EXPPOS_SHORT]) | 0x3FE0);
 2893             }
 2894         }
 2895         else if (!(*vl & 0x7FFF_FFFF_FFFF_FFFF))
 2896         {
 2897             // vf is +-0.0
 2898             exp = 0;
 2899         }
 2900         else
 2901         {
 2902             // subnormal
 2903             vf *= F.RECIP_EPSILON;
 2904             ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
 2905             exp = ((ex - F.EXPBIAS) >> 4) - T.mant_dig + 1;
 2906             vu[F.EXPPOS_SHORT] =
 2907                 cast(ushort)(((-1 - F.EXPMASK) & vu[F.EXPPOS_SHORT]) | 0x3FE0);
 2908         }
 2909         return vf;
 2910     }
 2911     else static if (F.realFormat == RealFormat.ieeeSingle)
 2912     {
 2913         if (ex) // If exponent is non-zero
 2914         {
 2915             if (ex == F.EXPMASK)   // infinity or NaN
 2916             {
 2917                 if (*vi == 0x7F80_0000)  // positive infinity
 2918                 {
 2919                     exp = int.max;
 2920                 }
 2921                 else if (*vi == 0xFF80_0000) // negative infinity
 2922                     exp = int.min;
 2923                 else
 2924                 { // NaN
 2925                     *vi |= 0x0040_0000;  // convert NaNS to NaNQ
 2926                     exp = int.min;
 2927                 }
 2928             }
 2929             else
 2930             {
 2931                 exp = (ex - F.EXPBIAS) >> 7;
 2932                 vu[F.EXPPOS_SHORT] = cast(ushort)((0x807F & vu[F.EXPPOS_SHORT]) | 0x3F00);
 2933             }
 2934         }
 2935         else if (!(*vi & 0x7FFF_FFFF))
 2936         {
 2937             // vf is +-0.0
 2938             exp = 0;
 2939         }
 2940         else
 2941         {
 2942             // subnormal
 2943             vf *= F.RECIP_EPSILON;
 2944             ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
 2945             exp = ((ex - F.EXPBIAS) >> 7) - T.mant_dig + 1;
 2946             vu[F.EXPPOS_SHORT] =
 2947                 cast(ushort)(((-1 - F.EXPMASK) & vu[F.EXPPOS_SHORT]) | 0x3F00);
 2948         }
 2949         return vf;
 2950     }
 2951     else // static if (F.realFormat == RealFormat.ibmExtended)
 2952     {
 2953         assert(0, "frexp not implemented");
 2954     }
 2955 }
 2956 
 2957 ///
 2958 @safe unittest
 2959 {
 2960     int exp;
 2961     real mantissa = frexp(123.456L, exp);
 2962 
 2963     assert(approxEqual(mantissa * pow(2.0L, cast(real) exp), 123.456L));
 2964 
 2965     assert(frexp(-real.nan, exp) && exp == int.min);
 2966     assert(frexp(real.nan, exp) && exp == int.min);
 2967     assert(frexp(-real.infinity, exp) == -real.infinity && exp == int.min);
 2968     assert(frexp(real.infinity, exp) == real.infinity && exp == int.max);
 2969     assert(frexp(-0.0, exp) == -0.0 && exp == 0);
 2970     assert(frexp(0.0, exp) == 0.0 && exp == 0);
 2971 }
 2972 
 2973 @safe @nogc nothrow unittest
 2974 {
 2975     int exp;
 2976     real mantissa = frexp(123.456L, exp);
 2977 
 2978     // check if values are equal to 19 decimal digits of precision
 2979     assert(equalsDigit(mantissa * pow(2.0L, cast(real) exp), 123.456L, 19));
 2980 }
 2981 
 2982 @safe unittest
 2983 {
 2984     import std.meta : AliasSeq;
 2985     import std.typecons : tuple, Tuple;
 2986 
 2987     static foreach (T; AliasSeq!(real, double, float))
 2988     {{
 2989         Tuple!(T, T, int)[] vals =     // x,frexp,exp
 2990             [
 2991              tuple(T(0.0),  T( 0.0 ), 0),
 2992              tuple(T(-0.0), T( -0.0), 0),
 2993              tuple(T(1.0),  T( .5  ), 1),
 2994              tuple(T(-1.0), T( -.5 ), 1),
 2995              tuple(T(2.0),  T( .5  ), 2),
 2996              tuple(T(float.min_normal/2.0f), T(.5), -126),
 2997              tuple(T.infinity, T.infinity, int.max),
 2998              tuple(-T.infinity, -T.infinity, int.min),
 2999              tuple(T.nan, T.nan, int.min),
 3000              tuple(-T.nan, -T.nan, int.min),
 3001 
 3002              // https://issues.dlang.org/show_bug.cgi?id=16026:
 3003              tuple(3 * (T.min_normal * T.epsilon), T( .75), (T.min_exp - T.mant_dig) + 2)
 3004              ];
 3005 
 3006         foreach (elem; vals)
 3007         {
 3008             T x = elem[0];
 3009             T e = elem[1];
 3010             int exp = elem[2];
 3011             int eptr;
 3012             T v = frexp(x, eptr);
 3013             assert(isIdentical(e, v));
 3014             assert(exp == eptr);
 3015 
 3016         }
 3017 
 3018         static if (floatTraits!(T).realFormat == RealFormat.ieeeExtended)
 3019         {
 3020             static T[3][] extendedvals = [ // x,frexp,exp
 3021                 [0x1.a5f1c2eb3fe4efp+73L,    0x1.A5F1C2EB3FE4EFp-1L,     74],    // normal
 3022                 [0x1.fa01712e8f0471ap-1064L, 0x1.fa01712e8f0471ap-1L, -1063],
 3023                 [T.min_normal,      .5, -16381],
 3024                 [T.min_normal/2.0L, .5, -16382]    // subnormal
 3025             ];
 3026             foreach (elem; extendedvals)
 3027             {
 3028                 T x = elem[0];
 3029                 T e = elem[1];
 3030                 int exp = cast(int) elem[2];
 3031                 int eptr;
 3032                 T v = frexp(x, eptr);
 3033                 assert(isIdentical(e, v));
 3034                 assert(exp == eptr);
 3035 
 3036             }
 3037         }
 3038     }}
 3039 
 3040     // CTFE
 3041     alias CtfeFrexpResult= Tuple!(real, int);
 3042     static CtfeFrexpResult ctfeFrexp(T)(const T value)
 3043     {
 3044         int exp;
 3045         auto significand = frexp(value, exp);
 3046         return CtfeFrexpResult(significand, exp);
 3047     }
 3048     static foreach (T; AliasSeq!(real, double, float))
 3049     {{
 3050         enum Tuple!(T, T, int)[] vals =     // x,frexp,exp
 3051             [
 3052              tuple(T(0.0),  T( 0.0 ), 0),
 3053              tuple(T(-0.0), T( -0.0), 0),
 3054              tuple(T(1.0),  T( .5  ), 1),
 3055              tuple(T(-1.0), T( -.5 ), 1),
 3056              tuple(T(2.0),  T( .5  ), 2),
 3057              tuple(T(float.min_normal/2.0f), T(.5), -126),
 3058              tuple(T.infinity, T.infinity, int.max),
 3059              tuple(-T.infinity, -T.infinity, int.min),
 3060              tuple(T.nan, T.nan, int.min),
 3061              tuple(-T.nan, -T.nan, int.min),
 3062 
 3063              // https://issues.dlang.org/show_bug.cgi?id=16026:
 3064              tuple(3 * (T.min_normal * T.epsilon), T( .75), (T.min_exp - T.mant_dig) + 2)
 3065              ];
 3066 
 3067         static foreach (elem; vals)
 3068         {
 3069             static assert(ctfeFrexp(elem[0]) is CtfeFrexpResult(elem[1], elem[2]));
 3070         }
 3071 
 3072         static if (floatTraits!(T).realFormat == RealFormat.ieeeExtended)
 3073         {
 3074             enum T[3][] extendedvals = [ // x,frexp,exp
 3075                 [0x1.a5f1c2eb3fe4efp+73L,    0x1.A5F1C2EB3FE4EFp-1L,     74],    // normal
 3076                 [0x1.fa01712e8f0471ap-1064L, 0x1.fa01712e8f0471ap-1L, -1063],
 3077                 [T.min_normal,      .5, -16381],
 3078                 [T.min_normal/2.0L, .5, -16382]    // subnormal
 3079             ];
 3080             static foreach (elem; extendedvals)
 3081             {
 3082                 static assert(ctfeFrexp(elem[0]) is CtfeFrexpResult(elem[1], cast(int) elem[2]));
 3083             }
 3084         }
 3085     }}
 3086 }
 3087 
 3088 @safe unittest
 3089 {
 3090     import std.meta : AliasSeq;
 3091     void foo() {
 3092         static foreach (T; AliasSeq!(real, double, float))
 3093         {{
 3094             int exp;
 3095             const T a = 1;
 3096             immutable T b = 2;
 3097             auto c = frexp(a, exp);
 3098             auto d = frexp(b, exp);
 3099         }}
 3100     }
 3101 }
 3102 
 3103 /******************************************
 3104  * Extracts the exponent of x as a signed integral value.
 3105  *
 3106  * If x is not a special value, the result is the same as
 3107  * $(D cast(int) logb(x)).
 3108  *
 3109  *      $(TABLE_SV
 3110  *      $(TR $(TH x)                $(TH ilogb(x))     $(TH Range error?))
 3111  *      $(TR $(TD 0)                 $(TD FP_ILOGB0)   $(TD yes))
 3112  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD int.max)     $(TD no))
 3113  *      $(TR $(TD $(NAN))            $(TD FP_ILOGBNAN) $(TD no))
 3114  *      )
 3115  */
 3116 int ilogb(T)(const T x) @trusted pure nothrow @nogc
 3117 if (isFloatingPoint!T)
 3118 {
 3119     import core.bitop : bsr;
 3120     alias F = floatTraits!T;
 3121 
 3122     union floatBits
 3123     {
 3124         T rv;
 3125         ushort[T.sizeof/2] vu;
 3126         uint[T.sizeof/4] vui;
 3127         static if (T.sizeof >= 8)
 3128             ulong[T.sizeof/8] vul;
 3129     }
 3130     floatBits y = void;
 3131     y.rv = x;
 3132 
 3133     int ex = y.vu[F.EXPPOS_SHORT] & F.EXPMASK;
 3134     static if (F.realFormat == RealFormat.ieeeExtended)
 3135     {
 3136         if (ex)
 3137         {
 3138             // If exponent is non-zero
 3139             if (ex == F.EXPMASK) // infinity or NaN
 3140             {
 3141                 if (y.vul[0] &  0x7FFF_FFFF_FFFF_FFFF)  // NaN
 3142                     return FP_ILOGBNAN;
 3143                 else // +-infinity
 3144                     return int.max;
 3145             }
 3146             else
 3147             {
 3148                 return ex - F.EXPBIAS - 1;
 3149             }
 3150         }
 3151         else if (!y.vul[0])
 3152         {
 3153             // vf is +-0.0
 3154             return FP_ILOGB0;
 3155         }
 3156         else
 3157         {
 3158             // subnormal
 3159             return ex - F.EXPBIAS - T.mant_dig + 1 + bsr(y.vul[0]);
 3160         }
 3161     }
 3162     else static if (F.realFormat == RealFormat.ieeeQuadruple)
 3163     {
 3164         if (ex)    // If exponent is non-zero
 3165         {
 3166             if (ex == F.EXPMASK)
 3167             {
 3168                 // infinity or NaN
 3169                 if (y.vul[MANTISSA_LSB] | ( y.vul[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF))  // NaN
 3170                     return FP_ILOGBNAN;
 3171                 else // +- infinity
 3172                     return int.max;
 3173             }
 3174             else
 3175             {
 3176                 return ex - F.EXPBIAS - 1;
 3177             }
 3178         }
 3179         else if ((y.vul[MANTISSA_LSB] | (y.vul[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) == 0)
 3180         {
 3181             // vf is +-0.0
 3182             return FP_ILOGB0;
 3183         }
 3184         else
 3185         {
 3186             // subnormal
 3187             const ulong msb = y.vul[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF;
 3188             const ulong lsb = y.vul[MANTISSA_LSB];
 3189             if (msb)
 3190                 return ex - F.EXPBIAS - T.mant_dig + 1 + bsr(msb) + 64;
 3191             else
 3192                 return ex - F.EXPBIAS - T.mant_dig + 1 + bsr(lsb);
 3193         }
 3194     }
 3195     else static if (F.realFormat == RealFormat.ieeeDouble)
 3196     {
 3197         if (ex) // If exponent is non-zero
 3198         {
 3199             if (ex == F.EXPMASK)   // infinity or NaN
 3200             {
 3201                 if ((y.vul[0] & 0x7FFF_FFFF_FFFF_FFFF) == 0x7FF0_0000_0000_0000)  // +- infinity
 3202                     return int.max;
 3203                 else // NaN
 3204                     return FP_ILOGBNAN;
 3205             }
 3206             else
 3207             {
 3208                 return ((ex - F.EXPBIAS) >> 4) - 1;
 3209             }
 3210         }
 3211         else if (!(y.vul[0] & 0x7FFF_FFFF_FFFF_FFFF))
 3212         {
 3213             // vf is +-0.0
 3214             return FP_ILOGB0;
 3215         }
 3216         else
 3217         {
 3218             // subnormal
 3219             enum MANTISSAMASK_64 = ((cast(ulong) F.MANTISSAMASK_INT) << 32) | 0xFFFF_FFFF;
 3220             return ((ex - F.EXPBIAS) >> 4) - T.mant_dig + 1 + bsr(y.vul[0] & MANTISSAMASK_64);
 3221         }
 3222     }
 3223     else static if (F.realFormat == RealFormat.ieeeSingle)
 3224     {
 3225         if (ex) // If exponent is non-zero
 3226         {
 3227             if (ex == F.EXPMASK)   // infinity or NaN
 3228             {
 3229                 if ((y.vui[0] & 0x7FFF_FFFF) == 0x7F80_0000)  // +- infinity
 3230                     return int.max;
 3231                 else // NaN
 3232                     return FP_ILOGBNAN;
 3233             }
 3234             else
 3235             {
 3236                 return ((ex - F.EXPBIAS) >> 7) - 1;
 3237             }
 3238         }
 3239         else if (!(y.vui[0] & 0x7FFF_FFFF))
 3240         {
 3241             // vf is +-0.0
 3242             return FP_ILOGB0;
 3243         }
 3244         else
 3245         {
 3246             // subnormal
 3247             const uint mantissa = y.vui[0] & F.MANTISSAMASK_INT;
 3248             return ((ex - F.EXPBIAS) >> 7) - T.mant_dig + 1 + bsr(mantissa);
 3249         }
 3250     }
 3251     else // static if (F.realFormat == RealFormat.ibmExtended)
 3252     {
 3253         assert(0, "ilogb not implemented");
 3254     }
 3255 }
 3256 /// ditto
 3257 int ilogb(T)(const T x) @safe pure nothrow @nogc
 3258 if (isIntegral!T && isUnsigned!T)
 3259 {
 3260     import core.bitop : bsr;
 3261     if (x == 0)
 3262         return FP_ILOGB0;
 3263     else
 3264     {
 3265         static assert(T.sizeof <= ulong.sizeof, "integer size too large for the current ilogb implementation");
 3266         return bsr(x);
 3267     }
 3268 }
 3269 /// ditto
 3270 int ilogb(T)(const T x) @safe pure nothrow @nogc
 3271 if (isIntegral!T && isSigned!T)
 3272 {
 3273     import std.traits : Unsigned;
 3274     // Note: abs(x) can not be used because the return type is not Unsigned and
 3275     //       the return value would be wrong for x == int.min
 3276     Unsigned!T absx =  x >= 0 ? x : -x;
 3277     return ilogb(absx);
 3278 }
 3279 
 3280 ///
 3281 @safe pure unittest
 3282 {
 3283     assert(ilogb(1) == 0);
 3284     assert(ilogb(3) == 1);
 3285     assert(ilogb(3.0) == 1);
 3286     assert(ilogb(100_000_000) == 26);
 3287 
 3288     assert(ilogb(0) == FP_ILOGB0);
 3289     assert(ilogb(0.0) == FP_ILOGB0);
 3290     assert(ilogb(double.nan) == FP_ILOGBNAN);
 3291     assert(ilogb(double.infinity) == int.max);
 3292 }
 3293 
 3294 /**
 3295 Special return values of $(LREF ilogb).
 3296  */
 3297 alias FP_ILOGB0   = core.stdc.math.FP_ILOGB0;
 3298 /// ditto
 3299 alias FP_ILOGBNAN = core.stdc.math.FP_ILOGBNAN;
 3300 
 3301 ///
 3302 @safe pure unittest
 3303 {
 3304     assert(ilogb(0) == FP_ILOGB0);
 3305     assert(ilogb(0.0) == FP_ILOGB0);
 3306     assert(ilogb(double.nan) == FP_ILOGBNAN);
 3307 }
 3308 
 3309 @safe nothrow @nogc unittest
 3310 {
 3311     import std.meta : AliasSeq;
 3312     import std.typecons : Tuple;
 3313     static foreach (F; AliasSeq!(float, double, real))
 3314     {{
 3315         alias T = Tuple!(F, int);
 3316         T[13] vals =   // x, ilogb(x)
 3317         [
 3318             T(  F.nan     , FP_ILOGBNAN ),
 3319             T( -F.nan     , FP_ILOGBNAN ),
 3320             T(  F.infinity, int.max     ),
 3321             T( -F.infinity, int.max     ),
 3322             T(  0.0       , FP_ILOGB0   ),
 3323             T( -0.0       , FP_ILOGB0   ),
 3324             T(  2.0       , 1           ),
 3325             T(  2.0001    , 1           ),
 3326             T(  1.9999    , 0           ),
 3327             T(  0.5       , -1          ),
 3328             T(  123.123   , 6           ),
 3329             T( -123.123   , 6           ),
 3330             T(  0.123     , -4          ),
 3331         ];
 3332 
 3333         foreach (elem; vals)
 3334         {
 3335             assert(ilogb(elem[0]) == elem[1]);
 3336         }
 3337     }}
 3338 
 3339     // min_normal and subnormals
 3340     assert(ilogb(-float.min_normal) == -126);
 3341     assert(ilogb(nextUp(-float.min_normal)) == -127);
 3342     assert(ilogb(nextUp(-float(0.0))) == -149);
 3343     assert(ilogb(-double.min_normal) == -1022);
 3344     assert(ilogb(nextUp(-double.min_normal)) == -1023);
 3345     assert(ilogb(nextUp(-double(0.0))) == -1074);
 3346     static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended)
 3347     {
 3348         assert(ilogb(-real.min_normal) == -16382);
 3349         assert(ilogb(nextUp(-real.min_normal)) == -16383);
 3350         assert(ilogb(nextUp(-real(0.0))) == -16445);
 3351     }
 3352     else static if (floatTraits!(real).realFormat == RealFormat.ieeeDouble)
 3353     {
 3354         assert(ilogb(-real.min_normal) == -1022);
 3355         assert(ilogb(nextUp(-real.min_normal)) == -1023);
 3356         assert(ilogb(nextUp(-real(0.0))) == -1074);
 3357     }
 3358 
 3359     // test integer types
 3360     assert(ilogb(0) == FP_ILOGB0);
 3361     assert(ilogb(int.max) == 30);
 3362     assert(ilogb(int.min) == 31);
 3363     assert(ilogb(uint.max) == 31);
 3364     assert(ilogb(long.max) == 62);
 3365     assert(ilogb(long.min) == 63);
 3366     assert(ilogb(ulong.max) == 63);
 3367 }
 3368 
 3369 /*******************************************
 3370  * Compute n * 2$(SUPERSCRIPT exp)
 3371  * References: frexp
 3372  */
 3373 
 3374 pragma(inline, true)
 3375 real ldexp(real n, int exp)     @safe pure nothrow @nogc { return core.math.ldexp(n, exp); }
 3376 ///ditto
 3377 pragma(inline, true)
 3378 double ldexp(double n, int exp) @safe pure nothrow @nogc { return core.math.ldexp(n, exp); }
 3379 ///ditto
 3380 pragma(inline, true)
 3381 float ldexp(float n, int exp)   @safe pure nothrow @nogc { return core.math.ldexp(n, exp); }
 3382 
 3383 ///
 3384 @nogc @safe pure nothrow unittest
 3385 {
 3386     import std.meta : AliasSeq;
 3387     static foreach (T; AliasSeq!(float, double, real))
 3388     {{
 3389         T r;
 3390 
 3391         r = ldexp(3.0L, 3);
 3392         assert(r == 24);
 3393 
 3394         r = ldexp(cast(T) 3.0, cast(int) 3);
 3395         assert(r == 24);
 3396 
 3397         T n = 3.0;
 3398         int exp = 3;
 3399         r = ldexp(n, exp);
 3400         assert(r == 24);
 3401     }}
 3402 }
 3403 
 3404 @safe pure nothrow @nogc unittest
 3405 {
 3406     static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended ||
 3407                floatTraits!(real).realFormat == RealFormat.ieeeQuadruple)
 3408     {
 3409         assert(ldexp(1.0L, -16384) == 0x1p-16384L);
 3410         assert(ldexp(1.0L, -16382) == 0x1p-16382L);
 3411         int x;
 3412         real n = frexp(0x1p-16384L, x);
 3413         assert(n == 0.5L);
 3414         assert(x==-16383);
 3415         assert(ldexp(n, x)==0x1p-16384L);
 3416     }
 3417     else static if (floatTraits!(real).realFormat == RealFormat.ieeeDouble)
 3418     {
 3419         assert(ldexp(1.0L, -1024) == 0x1p-1024L);
 3420         assert(ldexp(1.0L, -1022) == 0x1p-1022L);
 3421         int x;
 3422         real n = frexp(0x1p-1024L, x);
 3423         assert(n == 0.5L);
 3424         assert(x==-1023);
 3425         assert(ldexp(n, x)==0x1p-1024L);
 3426     }
 3427     else static assert(false, "Floating point type real not supported");
 3428 }
 3429 
 3430 /* workaround https://issues.dlang.org/show_bug.cgi?id=14718
 3431    float parsing depends on platform strtold
 3432 @safe pure nothrow @nogc unittest
 3433 {
 3434     assert(ldexp(1.0, -1024) == 0x1p-1024);
 3435     assert(ldexp(1.0, -1022) == 0x1p-1022);
 3436     int x;
 3437     double n = frexp(0x1p-1024, x);
 3438     assert(n == 0.5);
 3439     assert(x==-1023);
 3440     assert(ldexp(n, x)==0x1p-1024);
 3441 }
 3442 
 3443 @safe pure nothrow @nogc unittest
 3444 {
 3445     assert(ldexp(1.0f, -128) == 0x1p-128f);
 3446     assert(ldexp(1.0f, -126) == 0x1p-126f);
 3447     int x;
 3448     float n = frexp(0x1p-128f, x);
 3449     assert(n == 0.5f);
 3450     assert(x==-127);
 3451     assert(ldexp(n, x)==0x1p-128f);
 3452 }
 3453 */
 3454 
 3455 @safe @nogc nothrow unittest
 3456 {
 3457     static real[3][] vals =    // value,exp,ldexp
 3458     [
 3459     [    0,    0,    0],
 3460     [    1,    0,    1],
 3461     [    -1,    0,    -1],
 3462     [    1,    1,    2],
 3463     [    123,    10,    125952],
 3464     [    real.max,    int.max,    real.infinity],
 3465     [    real.max,    -int.max,    0],
 3466     [    real.min_normal,    -int.max,    0],
 3467     ];
 3468     int i;
 3469 
 3470     for (i = 0; i < vals.length; i++)
 3471     {
 3472         real x = vals[i][0];
 3473         int exp = cast(int) vals[i][1];
 3474         real z = vals[i][2];
 3475         real l = ldexp(x, exp);
 3476 
 3477         assert(equalsDigit(z, l, 7));
 3478     }
 3479 
 3480     real function(real, int) pldexp = &ldexp;
 3481     assert(pldexp != null);
 3482 }
 3483 
 3484 private
 3485 {
 3486     version (INLINE_YL2X) {} else
 3487     {
 3488         static if (floatTraits!real.realFormat == RealFormat.ieeeQuadruple)
 3489         {
 3490             // Coefficients for log(1 + x) = x - x**2/2 + x**3 P(x)/Q(x)
 3491             static immutable real[13] logCoeffsP = [
 3492                 1.313572404063446165910279910527789794488E4L,
 3493                 7.771154681358524243729929227226708890930E4L,
 3494                 2.014652742082537582487669938141683759923E5L,
 3495                 3.007007295140399532324943111654767187848E5L,
 3496                 2.854829159639697837788887080758954924001E5L,
 3497                 1.797628303815655343403735250238293741397E5L,
 3498                 7.594356839258970405033155585486712125861E4L,
 3499                 2.128857716871515081352991964243375186031E4L,
 3500                 3.824952356185897735160588078446136783779E3L,
 3501                 4.114517881637811823002128927449878962058E2L,
 3502                 2.321125933898420063925789532045674660756E1L,
 3503                 4.998469661968096229986658302195402690910E-1L,
 3504                 1.538612243596254322971797716843006400388E-6L
 3505             ];
 3506             static immutable real[13] logCoeffsQ = [
 3507                 3.940717212190338497730839731583397586124E4L,
 3508                 2.626900195321832660448791748036714883242E5L,
 3509                 7.777690340007566932935753241556479363645E5L,
 3510                 1.347518538384329112529391120390701166528E6L,
 3511                 1.514882452993549494932585972882995548426E6L,
 3512                 1.158019977462989115839826904108208787040E6L,
 3513                 6.132189329546557743179177159925690841200E5L,
 3514                 2.248234257620569139969141618556349415120E5L,
 3515                 5.605842085972455027590989944010492125825E4L,
 3516                 9.147150349299596453976674231612674085381E3L,
 3517                 9.104928120962988414618126155557301584078E2L,
 3518                 4.839208193348159620282142911143429644326E1L,
 3519                 1.0
 3520             ];
 3521 
 3522             // Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2)
 3523             // where z = 2(x-1)/(x+1)
 3524             static immutable real[6] logCoeffsR = [
 3525                 -8.828896441624934385266096344596648080902E-1L,
 3526                 8.057002716646055371965756206836056074715E1L,
 3527                 -2.024301798136027039250415126250455056397E3L,
 3528                 2.048819892795278657810231591630928516206E4L,
 3529                 -8.977257995689735303686582344659576526998E4L,
 3530                 1.418134209872192732479751274970992665513E5L
 3531             ];
 3532             static immutable real[6] logCoeffsS = [
 3533                 1.701761051846631278975701529965589676574E6L
 3534                 -1.332535117259762928288745111081235577029E6L,
 3535                 4.001557694070773974936904547424676279307E5L,
 3536                 -5.748542087379434595104154610899551484314E4L,
 3537                 3.998526750980007367835804959888064681098E3L,
 3538                 -1.186359407982897997337150403816839480438E2L,
 3539                 1.0
 3540             ];
 3541         }
 3542         else
 3543         {
 3544             // Coefficients for log(1 + x) = x - x**2/2 + x**3 P(x)/Q(x)
 3545             static immutable real[7] logCoeffsP = [
 3546                 2.0039553499201281259648E1L,
 3547                 5.7112963590585538103336E1L,
 3548                 6.0949667980987787057556E1L,
 3549                 2.9911919328553073277375E1L,
 3550                 6.5787325942061044846969E0L,
 3551                 4.9854102823193375972212E-1L,
 3552                 4.5270000862445199635215E-5L,
 3553             ];
 3554             static immutable real[7] logCoeffsQ = [
 3555                 6.0118660497603843919306E1L,
 3556                 2.1642788614495947685003E2L,
 3557                 3.0909872225312059774938E2L,
 3558                 2.2176239823732856465394E2L,
 3559                 8.3047565967967209469434E1L,
 3560                 1.5062909083469192043167E1L,
 3561                 1.0000000000000000000000E0L,
 3562             ];
 3563 
 3564             // Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2)
 3565             // where z = 2(x-1)/(x+1)
 3566             static immutable real[4] logCoeffsR = [
 3567                -3.5717684488096787370998E1L,
 3568                 1.0777257190312272158094E1L,
 3569                -7.1990767473014147232598E-1L,
 3570                 1.9757429581415468984296E-3L,
 3571             ];
 3572             static immutable real[4] logCoeffsS = [
 3573                -4.2861221385716144629696E2L,
 3574                 1.9361891836232102174846E2L,
 3575                -2.6201045551331104417768E1L,
 3576                 1.0000000000000000000000E0L,
 3577             ];
 3578         }
 3579     }
 3580 }
 3581 
 3582 /**************************************
 3583  * Calculate the natural logarithm of x.
 3584  *
 3585  *    $(TABLE_SV
 3586  *    $(TR $(TH x)            $(TH log(x))    $(TH divide by 0?) $(TH invalid?))
 3587  *    $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes)          $(TD no))
 3588  *    $(TR $(TD $(LT)0.0)     $(TD $(NAN))    $(TD no)           $(TD yes))
 3589  *    $(TR $(TD +$(INFIN))    $(TD +$(INFIN)) $(TD no)           $(TD no))
 3590  *    )
 3591  */
 3592 real log(real x) @safe pure nothrow @nogc
 3593 {
 3594     version (INLINE_YL2X)
 3595         return core.math.yl2x(x, LN2);
 3596     else
 3597     {
 3598         // C1 + C2 = LN2.
 3599         enum real C1 = 6.93145751953125E-1L;
 3600         enum real C2 = 1.428606820309417232121458176568075500134E-6L;
 3601 
 3602         // Special cases.
 3603         if (isNaN(x))
 3604             return x;
 3605         if (isInfinity(x) && !signbit(x))
 3606             return x;
 3607         if (x == 0.0)
 3608             return -real.infinity;
 3609         if (x < 0.0)
 3610             return real.nan;
 3611 
 3612         // Separate mantissa from exponent.
 3613         // Note, frexp is used so that denormal numbers will be handled properly.
 3614         real y, z;
 3615         int exp;
 3616 
 3617         x = frexp(x, exp);
 3618 
 3619         // Logarithm using log(x) = z + z^^3 R(z) / S(z),
 3620         // where z = 2(x - 1)/(x + 1)
 3621         if ((exp > 2) || (exp < -2))
 3622         {
 3623             if (x < SQRT1_2)
 3624             {   // 2(2x - 1)/(2x + 1)
 3625                 exp -= 1;
 3626                 z = x - 0.5;
 3627                 y = 0.5 * z + 0.5;
 3628             }
 3629             else
 3630             {   // 2(x - 1)/(x + 1)
 3631                 z = x - 0.5;
 3632                 z -= 0.5;
 3633                 y = 0.5 * x  + 0.5;
 3634             }
 3635             x = z / y;
 3636             z = x * x;
 3637             z = x * (z * poly(z, logCoeffsR) / poly(z, logCoeffsS));
 3638             z += exp * C2;
 3639             z += x;
 3640             z += exp * C1;
 3641 
 3642             return z;
 3643         }
 3644 
 3645         // Logarithm using log(1 + x) = x - .5x^^2 + x^^3 P(x) / Q(x)
 3646         if (x < SQRT1_2)
 3647         {
 3648             exp -= 1;
 3649             x = 2.0 * x - 1.0;
 3650         }
 3651         else
 3652         {
 3653             x = x - 1.0;
 3654         }
 3655         z = x * x;
 3656         y = x * (z * poly(x, logCoeffsP) / poly(x, logCoeffsQ));
 3657         y += exp * C2;
 3658         z = y - 0.5 * z;
 3659 
 3660         // Note, the sum of above terms does not exceed x/4,
 3661         // so it contributes at most about 1/4 lsb to the error.
 3662         z += x;
 3663         z += exp * C1;
 3664 
 3665         return z;
 3666     }
 3667 }
 3668 
 3669 ///
 3670 @safe pure nothrow @nogc unittest
 3671 {
 3672     assert(feqrel(log(E), 1) >= real.mant_dig - 1);
 3673 }
 3674 
 3675 /**************************************
 3676  * Calculate the base-10 logarithm of x.
 3677  *
 3678  *      $(TABLE_SV
 3679  *      $(TR $(TH x)            $(TH log10(x))  $(TH divide by 0?) $(TH invalid?))
 3680  *      $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes)          $(TD no))
 3681  *      $(TR $(TD $(LT)0.0)     $(TD $(NAN))    $(TD no)           $(TD yes))
 3682  *      $(TR $(TD +$(INFIN))    $(TD +$(INFIN)) $(TD no)           $(TD no))
 3683  *      )
 3684  */
 3685 real log10(real x) @safe pure nothrow @nogc
 3686 {
 3687     version (INLINE_YL2X)
 3688         return core.math.yl2x(x, LOG2);
 3689     else
 3690     {
 3691         // log10(2) split into two parts.
 3692         enum real L102A =  0.3125L;
 3693         enum real L102B = -1.14700043360188047862611052755069732318101185E-2L;
 3694 
 3695         // log10(e) split into two parts.
 3696         enum real L10EA =  0.5L;
 3697         enum real L10EB = -6.570551809674817234887108108339491770560299E-2L;
 3698 
 3699         // Special cases are the same as for log.
 3700         if (isNaN(x))
 3701             return x;
 3702         if (isInfinity(x) && !signbit(x))
 3703             return x;
 3704         if (x == 0.0)
 3705             return -real.infinity;
 3706         if (x < 0.0)
 3707             return real.nan;
 3708 
 3709         // Separate mantissa from exponent.
 3710         // Note, frexp is used so that denormal numbers will be handled properly.
 3711         real y, z;
 3712         int exp;
 3713 
 3714         x = frexp(x, exp);
 3715 
 3716         // Logarithm using log(x) = z + z^^3 R(z) / S(z),
 3717         // where z = 2(x - 1)/(x + 1)
 3718         if ((exp > 2) || (exp < -2))
 3719         {
 3720             if (x < SQRT1_2)
 3721             {   // 2(2x - 1)/(2x + 1)
 3722                 exp -= 1;
 3723                 z = x - 0.5;
 3724                 y = 0.5 * z + 0.5;
 3725             }
 3726             else
 3727             {   // 2(x - 1)/(x + 1)
 3728                 z = x - 0.5;
 3729                 z -= 0.5;
 3730                 y = 0.5 * x  + 0.5;
 3731             }
 3732             x = z / y;
 3733             z = x * x;
 3734             y = x * (z * poly(z, logCoeffsR) / poly(z, logCoeffsS));
 3735             goto Ldone;
 3736         }
 3737 
 3738         // Logarithm using log(1 + x) = x - .5x^^2 + x^^3 P(x) / Q(x)
 3739         if (x < SQRT1_2)
 3740         {
 3741             exp -= 1;
 3742             x = 2.0 * x - 1.0;
 3743         }
 3744         else
 3745             x = x - 1.0;
 3746 
 3747         z = x * x;
 3748         y = x * (z * poly(x, logCoeffsP) / poly(x, logCoeffsQ));
 3749         y = y - 0.5 * z;
 3750 
 3751         // Multiply log of fraction by log10(e) and base 2 exponent by log10(2).
 3752         // This sequence of operations is critical and it may be horribly
 3753         // defeated by some compiler optimizers.
 3754     Ldone:
 3755         z = y * L10EB;
 3756         z += x * L10EB;
 3757         z += exp * L102B;
 3758         z += y * L10EA;
 3759         z += x * L10EA;
 3760         z += exp * L102A;
 3761 
 3762         return z;
 3763     }
 3764 }
 3765 
 3766 ///
 3767 @safe pure nothrow @nogc unittest
 3768 {
 3769     assert(fabs(log10(1000) - 3) < .000001);
 3770 }
 3771 
 3772 /**
 3773  * Calculates the natural logarithm of 1 + x.
 3774  *
 3775  * For very small x, log1p(x) will be more accurate than
 3776  * log(1 + x).
 3777  *
 3778  *  $(TABLE_SV
 3779  *  $(TR $(TH x)            $(TH log1p(x))     $(TH divide by 0?) $(TH invalid?))
 3780  *  $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no)           $(TD no))
 3781  *  $(TR $(TD -1.0)         $(TD -$(INFIN))    $(TD yes)          $(TD no))
 3782  *  $(TR $(TD $(LT)-1.0)    $(TD -$(NAN))      $(TD no)           $(TD yes))
 3783  *  $(TR $(TD +$(INFIN))    $(TD +$(INFIN))    $(TD no)           $(TD no))
 3784  *  )
 3785  */
 3786 real log1p(real x) @safe pure nothrow @nogc
 3787 {
 3788     version (INLINE_YL2X)
 3789     {
 3790         // On x87, yl2xp1 is valid if and only if -0.5 <= lg(x) <= 0.5,
 3791         //    ie if -0.29 <= x <= 0.414
 3792         return (fabs(x) <= 0.25)  ? core.math.yl2xp1(x, LN2) : core.math.yl2x(x+1, LN2);
 3793     }
 3794     else
 3795     {
 3796         // Special cases.
 3797         if (isNaN(x) || x == 0.0)
 3798             return x;
 3799         if (isInfinity(x) && !signbit(x))
 3800             return x;
 3801         if (x == -1.0)
 3802             return -real.infinity;
 3803         if (x < -1.0)
 3804             return real.nan;
 3805 
 3806         return log(x + 1.0);
 3807     }
 3808 }
 3809 
 3810 ///
 3811 @safe pure unittest
 3812 {
 3813     assert(isIdentical(log1p(0.0), 0.0));
 3814     assert(log1p(1.0).feqrel(0.69314) > 16);
 3815 
 3816     assert(log1p(-1.0) == -real.infinity);
 3817     assert(isNaN(log1p(-2.0)));
 3818     assert(log1p(real.nan) is real.nan);
 3819     assert(log1p(-real.nan) is -real.nan);
 3820     assert(log1p(real.infinity) == real.infinity);
 3821 }
 3822 
 3823 /***************************************
 3824  * Calculates the base-2 logarithm of x:
 3825  * $(SUB log, 2)x
 3826  *
 3827  *  $(TABLE_SV
 3828  *  $(TR $(TH x)            $(TH log2(x))   $(TH divide by 0?) $(TH invalid?))
 3829  *  $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes)          $(TD no) )
 3830  *  $(TR $(TD $(LT)0.0)     $(TD $(NAN))    $(TD no)           $(TD yes) )
 3831  *  $(TR $(TD +$(INFIN))    $(TD +$(INFIN)) $(TD no)           $(TD no) )
 3832  *  )
 3833  */
 3834 real log2(real x) @safe pure nothrow @nogc
 3835 {
 3836     version (INLINE_YL2X)
 3837         return core.math.yl2x(x, 1.0L);
 3838     else
 3839     {
 3840         // Special cases are the same as for log.
 3841         if (isNaN(x))
 3842             return x;
 3843         if (isInfinity(x) && !signbit(x))
 3844             return x;
 3845         if (x == 0.0)
 3846             return -real.infinity;
 3847         if (x < 0.0)
 3848             return real.nan;
 3849 
 3850         // Separate mantissa from exponent.
 3851         // Note, frexp is used so that denormal numbers will be handled properly.
 3852         real y, z;
 3853         int exp;
 3854 
 3855         x = frexp(x, exp);
 3856 
 3857         // Logarithm using log(x) = z + z^^3 R(z) / S(z),
 3858         // where z = 2(x - 1)/(x + 1)
 3859         if ((exp > 2) || (exp < -2))
 3860         {
 3861             if (x < SQRT1_2)
 3862             {   // 2(2x - 1)/(2x + 1)
 3863                 exp -= 1;
 3864                 z = x - 0.5;
 3865                 y = 0.5 * z + 0.5;
 3866             }
 3867             else
 3868             {   // 2(x - 1)/(x + 1)
 3869                 z = x - 0.5;
 3870                 z -= 0.5;
 3871                 y = 0.5 * x  + 0.5;
 3872             }
 3873             x = z / y;
 3874             z = x * x;
 3875             y = x * (z * poly(z, logCoeffsR) / poly(z, logCoeffsS));
 3876             goto Ldone;
 3877         }
 3878 
 3879         // Logarithm using log(1 + x) = x - .5x^^2 + x^^3 P(x) / Q(x)
 3880         if (x < SQRT1_2)
 3881         {
 3882             exp -= 1;
 3883             x = 2.0 * x - 1.0;
 3884         }
 3885         else
 3886             x = x - 1.0;
 3887 
 3888         z = x * x;
 3889         y = x * (z * poly(x, logCoeffsP) / poly(x, logCoeffsQ));
 3890         y = y - 0.5 * z;
 3891 
 3892         // Multiply log of fraction by log10(e) and base 2 exponent by log10(2).
 3893         // This sequence of operations is critical and it may be horribly
 3894         // defeated by some compiler optimizers.
 3895     Ldone:
 3896         z = y * (LOG2E - 1.0);
 3897         z += x * (LOG2E - 1.0);
 3898         z += y;
 3899         z += x;
 3900         z += exp;
 3901 
 3902         return z;
 3903     }
 3904 }
 3905 
 3906 ///
 3907 @safe unittest
 3908 {
 3909     assert(approxEqual(log2(1024.0L), 10));
 3910 }
 3911 
 3912 @safe @nogc nothrow unittest
 3913 {
 3914     // check if values are equal to 19 decimal digits of precision
 3915     assert(equalsDigit(log2(1024.0L), 10, 19));
 3916 }
 3917 
 3918 /*****************************************
 3919  * Extracts the exponent of x as a signed integral value.
 3920  *
 3921  * If x is subnormal, it is treated as if it were normalized.
 3922  * For a positive, finite x:
 3923  *
 3924  * 1 $(LT)= $(I x) * FLT_RADIX$(SUPERSCRIPT -logb(x)) $(LT) FLT_RADIX
 3925  *
 3926  *      $(TABLE_SV
 3927  *      $(TR $(TH x)                 $(TH logb(x))   $(TH divide by 0?) )
 3928  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD +$(INFIN)) $(TD no))
 3929  *      $(TR $(TD $(PLUSMN)0.0)      $(TD -$(INFIN)) $(TD yes) )
 3930  *      )
 3931  */
 3932 real logb(real x) @trusted nothrow @nogc
 3933 {
 3934     version (Win64_DMD_InlineAsm)
 3935     {
 3936         asm pure nothrow @nogc
 3937         {
 3938             naked                       ;
 3939             fld     real ptr [RCX]      ;
 3940             fxtract                     ;
 3941             fstp    ST(0)               ;
 3942             ret                         ;
 3943         }
 3944     }
 3945     else version (MSVC_InlineAsm)
 3946     {
 3947         asm pure nothrow @nogc
 3948         {
 3949             fld     x                   ;
 3950             fxtract                     ;
 3951             fstp    ST(0)               ;
 3952         }
 3953     }
 3954     else
 3955         return core.stdc.math.logbl(x);
 3956 }
 3957 
 3958 ///
 3959 @safe @nogc nothrow unittest
 3960 {
 3961     assert(logb(1.0) == 0);
 3962     assert(logb(100.0) == 6);
 3963 
 3964     assert(logb(0.0) == -real.infinity);
 3965     assert(logb(real.infinity) == real.infinity);
 3966     assert(logb(-real.infinity) == real.infinity);
 3967 }
 3968 
 3969 /************************************
 3970  * Calculates the remainder from the calculation x/y.
 3971  * Returns:
 3972  * The value of x - i * y, where i is the number of times that y can
 3973  * be completely subtracted from x. The result has the same sign as x.
 3974  *
 3975  * $(TABLE_SV
 3976  *  $(TR $(TH x)              $(TH y)             $(TH fmod(x, y))   $(TH invalid?))
 3977  *  $(TR $(TD $(PLUSMN)0.0)   $(TD not 0.0)       $(TD $(PLUSMN)0.0) $(TD no))
 3978  *  $(TR $(TD $(PLUSMNINF))   $(TD anything)      $(TD $(NAN))       $(TD yes))
 3979  *  $(TR $(TD anything)       $(TD $(PLUSMN)0.0)  $(TD $(NAN))       $(TD yes))
 3980  *  $(TR $(TD !=$(PLUSMNINF)) $(TD $(PLUSMNINF))  $(TD x)            $(TD no))
 3981  * )
 3982  */
 3983 real fmod(real x, real y) @trusted nothrow @nogc
 3984 {
 3985     version (CRuntime_Microsoft)
 3986     {
 3987         return x % y;
 3988     }
 3989     else
 3990         return core.stdc.math.fmodl(x, y);
 3991 }
 3992 
 3993 ///
 3994 @safe unittest
 3995 {
 3996     assert(isIdentical(fmod(0.0, 1.0), 0.0));
 3997     assert(fmod(5.0, 3.0).feqrel(2.0) > 16);
 3998     assert(isNaN(fmod(5.0, 0.0)));
 3999 }
 4000 
 4001 /************************************
 4002  * Breaks x into an integral part and a fractional part, each of which has
 4003  * the same sign as x. The integral part is stored in i.
 4004  * Returns:
 4005  * The fractional part of x.
 4006  *
 4007  * $(TABLE_SV
 4008  *  $(TR $(TH x)              $(TH i (on input))  $(TH modf(x, i))   $(TH i (on return)))
 4009  *  $(TR $(TD $(PLUSMNINF))   $(TD anything)      $(TD $(PLUSMN)0.0) $(TD $(PLUSMNINF)))
 4010  * )
 4011  */
 4012 real modf(real x, ref real i) @trusted nothrow @nogc
 4013 {
 4014     version (CRuntime_Microsoft)
 4015     {
 4016         i = trunc(x);
 4017         return copysign(isInfinity(x) ? 0.0 : x - i, x);
 4018     }
 4019     else
 4020         return core.stdc.math.modfl(x,&i);
 4021 }
 4022 
 4023 ///
 4024 @safe unittest
 4025 {
 4026     real frac;
 4027     real intpart;
 4028 
 4029     frac = modf(3.14159, intpart);
 4030     assert(intpart.feqrel(3.0) > 16);
 4031     assert(frac.feqrel(0.14159) > 16);
 4032 }
 4033 
 4034 /*************************************
 4035  * Efficiently calculates x * 2$(SUPERSCRIPT n).
 4036  *
 4037  * scalbn handles underflow and overflow in
 4038  * the same fashion as the basic arithmetic operators.
 4039  *
 4040  *      $(TABLE_SV
 4041  *      $(TR $(TH x)                 $(TH scalb(x)))
 4042  *      $(TR $(TD $(PLUSMNINF))      $(TD $(PLUSMNINF)) )
 4043  *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(PLUSMN)0.0) )
 4044  *      )
 4045  */
 4046 real scalbn(real x, int n) @safe pure nothrow @nogc
 4047 {
 4048     if (__ctfe)
 4049     {
 4050         // Handle special cases.
 4051         if (x == 0.0 || isInfinity(x))
 4052             return x;
 4053     }
 4054     return core.math.ldexp(x, n);
 4055 }
 4056 
 4057 ///
 4058 @safe pure nothrow @nogc unittest
 4059 {
 4060     assert(scalbn(0x1.2345678abcdefp0, 999) == 0x1.2345678abcdefp999);
 4061     assert(scalbn(-real.infinity, 5) == -real.infinity);
 4062 }
 4063 
 4064 @safe pure nothrow @nogc unittest
 4065 {
 4066     // CTFE-able test
 4067     static assert(scalbn(0x1.2345678abcdefp0, 999) == 0x1.2345678abcdefp999);
 4068     static assert(scalbn(-real.infinity, 5) == -real.infinity);
 4069     // Test with large exponent delta n where the result is in bounds but 2.0L ^^ n is not.
 4070     enum initialExponent = real.min_exp + 2, resultExponent = real.max_exp - 2;
 4071     enum int n = resultExponent - initialExponent;
 4072     enum real x = 0x1.2345678abcdefp0 * (2.0L ^^ initialExponent);
 4073     enum staticResult = scalbn(x, n);
 4074     static assert(staticResult == 0x1.2345678abcdefp0 * (2.0L ^^ resultExponent));
 4075     assert(scalbn(x, n) == staticResult);
 4076 }
 4077 
 4078 /***************
 4079  * Calculates the cube root of x.
 4080  *
 4081  *      $(TABLE_SV
 4082  *      $(TR $(TH $(I x))            $(TH cbrt(x))           $(TH invalid?))
 4083  *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(PLUSMN)0.0)      $(TD no) )
 4084  *      $(TR $(TD $(NAN))            $(TD $(NAN))            $(TD yes) )
 4085  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)$(INFIN)) $(TD no) )
 4086  *      )
 4087  */
 4088 real cbrt(real x) @trusted nothrow @nogc
 4089 {
 4090     version (CRuntime_Microsoft)
 4091     {
 4092         version (INLINE_YL2X)
 4093             return copysign(exp2(core.math.yl2x(fabs(x), 1.0L/3.0L)), x);
 4094         else
 4095             return core.stdc.math.cbrtl(x);
 4096     }
 4097     else
 4098         return core.stdc.math.cbrtl(x);
 4099 }
 4100 
 4101 ///
 4102 @safe unittest
 4103 {
 4104     assert(cbrt(1.0).feqrel(1.0) > 16);
 4105     assert(cbrt(27.0).feqrel(3.0) > 16);
 4106     assert(cbrt(15.625).feqrel(2.5) > 16);
 4107 }
 4108 
 4109 /*******************************
 4110  * Returns |x|
 4111  *
 4112  *      $(TABLE_SV
 4113  *      $(TR $(TH x)                 $(TH fabs(x)))
 4114  *      $(TR $(TD $(PLUSMN)0.0)      $(TD +0.0) )
 4115  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD +$(INFIN)) )
 4116  *      )
 4117  */
 4118 real fabs(real x) @safe pure nothrow @nogc { pragma(inline, true); return core.math.fabs(x); }
 4119 
 4120 ///ditto
 4121 pragma(inline, true)
 4122 double fabs(double d) @trusted pure nothrow @nogc
 4123 {
 4124     ulong tmp = *cast(ulong*)&d & 0x7FFF_FFFF_FFFF_FFFF;
 4125     return *cast(double*)&tmp;
 4126 }
 4127 
 4128 ///ditto
 4129 pragma(inline, true)
 4130 float fabs(float f) @trusted pure nothrow @nogc
 4131 {
 4132     uint tmp = *cast(uint*)&f & 0x7FFF_FFFF;
 4133     return *cast(float*)&tmp;
 4134 }
 4135 
 4136 ///
 4137 @safe unittest
 4138 {
 4139 
 4140     assert(isIdentical(fabs(0.0f), 0.0f));
 4141     assert(isIdentical(fabs(-0.0f), 0.0f));
 4142     assert(fabs(-10.0f) == 10.0f);
 4143 
 4144     assert(isIdentical(fabs(0.0), 0.0));
 4145     assert(isIdentical(fabs(-0.0), 0.0));
 4146     assert(fabs(-10.0) == 10.0);
 4147 
 4148     assert(isIdentical(fabs(0.0L), 0.0L));
 4149     assert(isIdentical(fabs(-0.0L), 0.0L));
 4150     assert(fabs(-10.0L) == 10.0L);
 4151 }
 4152 
 4153 @safe unittest
 4154 {
 4155     real function(real) pfabs = &fabs;
 4156     assert(pfabs != null);
 4157 }
 4158 
 4159 /***********************************************************************
 4160  * Calculates the length of the
 4161  * hypotenuse of a right-angled triangle with sides of length x and y.
 4162  * The hypotenuse is the value of the square root of
 4163  * the sums of the squares of x and y:
 4164  *
 4165  *      sqrt($(POWER x, 2) + $(POWER y, 2))
 4166  *
 4167  * Note that hypot(x, y), hypot(y, x) and
 4168  * hypot(x, -y) are equivalent.
 4169  *
 4170  *  $(TABLE_SV
 4171  *  $(TR $(TH x)            $(TH y)            $(TH hypot(x, y)) $(TH invalid?))
 4172  *  $(TR $(TD x)            $(TD $(PLUSMN)0.0) $(TD |x|)         $(TD no))
 4173  *  $(TR $(TD $(PLUSMNINF)) $(TD y)            $(TD +$(INFIN))   $(TD no))
 4174  *  $(TR $(TD $(PLUSMNINF)) $(TD $(NAN))       $(TD +$(INFIN))   $(TD no))
 4175  *  )
 4176  */
 4177 
 4178 real hypot(real x, real y) @safe pure nothrow @nogc
 4179 {
 4180     // Scale x and y to avoid underflow and overflow.
 4181     // If one is huge and the other tiny, return the larger.
 4182     // If both are huge, avoid overflow by scaling by 1/sqrt(real.max/2).
 4183     // If both are tiny, avoid underflow by scaling by sqrt(real.min_normal*real.epsilon).
 4184 
 4185     enum real SQRTMIN = 0.5 * sqrt(real.min_normal); // This is a power of 2.
 4186     enum real SQRTMAX = 1.0L / SQRTMIN; // 2^^((max_exp)/2) = nextUp(sqrt(real.max))
 4187 
 4188     static assert(2*(SQRTMAX/2)*(SQRTMAX/2) <= real.max);
 4189 
 4190     // Proves that sqrt(real.max) ~~  0.5/sqrt(real.min_normal)
 4191     static assert(real.min_normal*real.max > 2 && real.min_normal*real.max <= 4);
 4192 
 4193     real u = fabs(x);
 4194     real v = fabs(y);
 4195     if (!(u >= v))  // check for NaN as well.
 4196     {
 4197         v = u;
 4198         u = fabs(y);
 4199         if (u == real.infinity) return u; // hypot(inf, nan) == inf
 4200         if (v == real.infinity) return v; // hypot(nan, inf) == inf
 4201     }
 4202 
 4203     // Now u >= v, or else one is NaN.
 4204     if (v >= SQRTMAX*0.5)
 4205     {
 4206             // hypot(huge, huge) -- avoid overflow
 4207         u *= SQRTMIN*0.5;
 4208         v *= SQRTMIN*0.5;
 4209         return sqrt(u*u + v*v) * SQRTMAX * 2.0;
 4210     }
 4211 
 4212     if (u <= SQRTMIN)
 4213     {
 4214         // hypot (tiny, tiny) -- avoid underflow
 4215         // This is only necessary to avoid setting the underflow
 4216         // flag.
 4217         u *= SQRTMAX / real.epsilon;
 4218         v *= SQRTMAX / real.epsilon;
 4219         return sqrt(u*u + v*v) * SQRTMIN * real.epsilon;
 4220     }
 4221 
 4222     if (u * real.epsilon > v)
 4223     {
 4224         // hypot (huge, tiny) = huge
 4225         return u;
 4226     }
 4227 
 4228     // both are in the normal range
 4229     return sqrt(u*u + v*v);
 4230 }
 4231 
 4232 ///
 4233 @safe unittest
 4234 {
 4235     assert(hypot(1.0, 1.0).feqrel(1.4142) > 16);
 4236     assert(hypot(3.0, 4.0).feqrel(5.0) > 16);
 4237     assert(hypot(real.infinity, 1.0) == real.infinity);
 4238     assert(hypot(real.infinity, real.nan) == real.infinity);
 4239 }
 4240 
 4241 @safe unittest
 4242 {
 4243     static real[3][] vals =     // x,y,hypot
 4244         [
 4245             [ 0.0,     0.0,   0.0],
 4246             [ 0.0,    -0.0,   0.0],
 4247             [ -0.0,   -0.0,   0.0],
 4248             [ 3.0,     4.0,   5.0],
 4249             [ -300,   -400,   500],
 4250             [0.0,      7.0,   7.0],
 4251             [9.0,   9*real.epsilon,   9.0],
 4252             [88/(64*sqrt(real.min_normal)), 105/(64*sqrt(real.min_normal)), 137/(64*sqrt(real.min_normal))],
 4253             [88/(128*sqrt(real.min_normal)), 105/(128*sqrt(real.min_normal)), 137/(128*sqrt(real.min_normal))],
 4254             [3*real.min_normal*real.epsilon, 4*real.min_normal*real.epsilon, 5*real.min_normal*real.epsilon],
 4255             [ real.min_normal, real.min_normal, sqrt(2.0L)*real.min_normal],
 4256             [ real.max/sqrt(2.0L), real.max/sqrt(2.0L), real.max],
 4257             [ real.infinity, real.nan, real.infinity],
 4258             [ real.nan, real.infinity, real.infinity],
 4259             [ real.nan, real.nan, real.nan],
 4260             [ real.nan, real.max, real.nan],
 4261             [ real.max, real.nan, real.nan],
 4262         ];
 4263         for (int i = 0; i < vals.length; i++)
 4264         {
 4265             real x = vals[i][0];
 4266             real y = vals[i][1];
 4267             real z = vals[i][2];
 4268             real h = hypot(x, y);
 4269             assert(isIdentical(z,h) || feqrel(z, h) >= real.mant_dig - 1);
 4270         }
 4271 }
 4272 
 4273 /**************************************
 4274  * Returns the value of x rounded upward to the next integer
 4275  * (toward positive infinity).
 4276  */
 4277 real ceil(real x) @trusted pure nothrow @nogc
 4278 {
 4279     version (Win64_DMD_InlineAsm)
 4280     {
 4281         asm pure nothrow @nogc
 4282         {
 4283             naked                       ;
 4284             fld     real ptr [RCX]      ;
 4285             fstcw   8[RSP]              ;
 4286             mov     AL,9[RSP]           ;
 4287             mov     DL,AL               ;
 4288             and     AL,0xC3             ;
 4289             or      AL,0x08             ; // round to +infinity
 4290             mov     9[RSP],AL           ;
 4291             fldcw   8[RSP]              ;
 4292             frndint                     ;
 4293             mov     9[RSP],DL           ;
 4294             fldcw   8[RSP]              ;
 4295             ret                         ;
 4296         }
 4297     }
 4298     else version (MSVC_InlineAsm)
 4299     {
 4300         short cw;
 4301         asm pure nothrow @nogc
 4302         {
 4303             fld     x                   ;
 4304             fstcw   cw                  ;
 4305             mov     AL,byte ptr cw+1    ;
 4306             mov     DL,AL               ;
 4307             and     AL,0xC3             ;
 4308             or      AL,0x08             ; // round to +infinity
 4309             mov     byte ptr cw+1,AL    ;
 4310             fldcw   cw                  ;
 4311             frndint                     ;
 4312             mov     byte ptr cw+1,DL    ;
 4313             fldcw   cw                  ;
 4314         }
 4315     }
 4316     else
 4317     {
 4318         // Special cases.
 4319         if (isNaN(x) || isInfinity(x))
 4320             return x;
 4321 
 4322         real y = floorImpl(x);
 4323         if (y < x)
 4324             y += 1.0;
 4325 
 4326         return y;
 4327     }
 4328 }
 4329 
 4330 ///
 4331 @safe pure nothrow @nogc unittest
 4332 {
 4333     assert(ceil(+123.456L) == +124);
 4334     assert(ceil(-123.456L) == -123);
 4335     assert(ceil(-1.234L) == -1);
 4336     assert(ceil(-0.123L) == 0);
 4337     assert(ceil(0.0L) == 0);
 4338     assert(ceil(+0.123L) == 1);
 4339     assert(ceil(+1.234L) == 2);
 4340     assert(ceil(real.infinity) == real.infinity);
 4341     assert(isNaN(ceil(real.nan)));
 4342     assert(isNaN(ceil(real.init)));
 4343 }
 4344 
 4345 /// ditto
 4346 double ceil(double x) @trusted pure nothrow @nogc
 4347 {
 4348     // Special cases.
 4349     if (isNaN(x) || isInfinity(x))
 4350         return x;
 4351 
 4352     double y = floorImpl(x);
 4353     if (y < x)
 4354         y += 1.0;
 4355 
 4356     return y;
 4357 }
 4358 
 4359 @safe pure nothrow @nogc unittest
 4360 {
 4361     assert(ceil(+123.456) == +124);
 4362     assert(ceil(-123.456) == -123);
 4363     assert(ceil(-1.234) == -1);
 4364     assert(ceil(-0.123) == 0);
 4365     assert(ceil(0.0) == 0);
 4366     assert(ceil(+0.123) == 1);
 4367     assert(ceil(+1.234) == 2);
 4368     assert(ceil(double.infinity) == double.infinity);
 4369     assert(isNaN(ceil(double.nan)));
 4370     assert(isNaN(ceil(double.init)));
 4371 }
 4372 
 4373 /// ditto
 4374 float ceil(float x) @trusted pure nothrow @nogc
 4375 {
 4376     // Special cases.
 4377     if (isNaN(x) || isInfinity(x))
 4378         return x;
 4379 
 4380     float y = floorImpl(x);
 4381     if (y < x)
 4382         y += 1.0;
 4383 
 4384     return y;
 4385 }
 4386 
 4387 @safe pure nothrow @nogc unittest
 4388 {
 4389     assert(ceil(+123.456f) == +124);
 4390     assert(ceil(-123.456f) == -123);
 4391     assert(ceil(-1.234f) == -1);
 4392     assert(ceil(-0.123f) == 0);
 4393     assert(ceil(0.0f) == 0);
 4394     assert(ceil(+0.123f) == 1);
 4395     assert(ceil(+1.234f) == 2);
 4396     assert(ceil(float.infinity) == float.infinity);
 4397     assert(isNaN(ceil(float.nan)));
 4398     assert(isNaN(ceil(float.init)));
 4399 }
 4400 
 4401 /**************************************
 4402  * Returns the value of x rounded downward to the next integer
 4403  * (toward negative infinity).
 4404  */
 4405 real floor(real x) @trusted pure nothrow @nogc
 4406 {
 4407     version (Win64_DMD_InlineAsm)
 4408     {
 4409         asm pure nothrow @nogc
 4410         {
 4411             naked                       ;
 4412             fld     real ptr [RCX]      ;
 4413             fstcw   8[RSP]              ;
 4414             mov     AL,9[RSP]           ;
 4415             mov     DL,AL               ;
 4416             and     AL,0xC3             ;
 4417             or      AL,0x04             ; // round to -infinity
 4418             mov     9[RSP],AL           ;
 4419             fldcw   8[RSP]              ;
 4420             frndint                     ;
 4421             mov     9[RSP],DL           ;
 4422             fldcw   8[RSP]              ;
 4423             ret                         ;
 4424         }
 4425     }
 4426     else version (MSVC_InlineAsm)
 4427     {
 4428         short cw;
 4429         asm pure nothrow @nogc
 4430         {
 4431             fld     x                   ;
 4432             fstcw   cw                  ;
 4433             mov     AL,byte ptr cw+1    ;
 4434             mov     DL,AL               ;
 4435             and     AL,0xC3             ;
 4436             or      AL,0x04             ; // round to -infinity
 4437             mov     byte ptr cw+1,AL    ;
 4438             fldcw   cw                  ;
 4439             frndint                     ;
 4440             mov     byte ptr cw+1,DL    ;
 4441             fldcw   cw                  ;
 4442         }
 4443     }
 4444     else
 4445     {
 4446         // Special cases.
 4447         if (isNaN(x) || isInfinity(x) || x == 0.0)
 4448             return x;
 4449 
 4450         return floorImpl(x);
 4451     }
 4452 }
 4453 
 4454 ///
 4455 @safe pure nothrow @nogc unittest
 4456 {
 4457     assert(floor(+123.456L) == +123);
 4458     assert(floor(-123.456L) == -124);
 4459     assert(floor(+123.0L) == +123);
 4460     assert(floor(-124.0L) == -124);
 4461     assert(floor(-1.234L) == -2);
 4462     assert(floor(-0.123L) == -1);
 4463     assert(floor(0.0L) == 0);
 4464     assert(floor(+0.123L) == 0);
 4465     assert(floor(+1.234L) == 1);
 4466     assert(floor(real.infinity) == real.infinity);
 4467     assert(isNaN(floor(real.nan)));
 4468     assert(isNaN(floor(real.init)));
 4469 }
 4470 
 4471 /// ditto
 4472 double floor(double x) @trusted pure nothrow @nogc
 4473 {
 4474     // Special cases.
 4475     if (isNaN(x) || isInfinity(x) || x == 0.0)
 4476         return x;
 4477 
 4478     return floorImpl(x);
 4479 }
 4480 
 4481 @safe pure nothrow @nogc unittest
 4482 {
 4483     assert(floor(+123.456) == +123);
 4484     assert(floor(-123.456) == -124);
 4485     assert(floor(+123.0) == +123);
 4486     assert(floor(-124.0) == -124);
 4487     assert(floor(-1.234) == -2);
 4488     assert(floor(-0.123) == -1);
 4489     assert(floor(0.0) == 0);
 4490     assert(floor(+0.123) == 0);
 4491     assert(floor(+1.234) == 1);
 4492     assert(floor(double.infinity) == double.infinity);
 4493     assert(isNaN(floor(double.nan)));
 4494     assert(isNaN(floor(double.init)));
 4495 }
 4496 
 4497 /// ditto
 4498 float floor(float x) @trusted pure nothrow @nogc
 4499 {
 4500     // Special cases.
 4501     if (isNaN(x) || isInfinity(x) || x == 0.0)
 4502         return x;
 4503 
 4504     return floorImpl(x);
 4505 }
 4506 
 4507 @safe pure nothrow @nogc unittest
 4508 {
 4509     assert(floor(+123.456f) == +123);
 4510     assert(floor(-123.456f) == -124);
 4511     assert(floor(+123.0f) == +123);
 4512     assert(floor(-124.0f) == -124);
 4513     assert(floor(-1.234f) == -2);
 4514     assert(floor(-0.123f) == -1);
 4515     assert(floor(0.0f) == 0);
 4516     assert(floor(+0.123f) == 0);
 4517     assert(floor(+1.234f) == 1);
 4518     assert(floor(float.infinity) == float.infinity);
 4519     assert(isNaN(floor(float.nan)));
 4520     assert(isNaN(floor(float.init)));
 4521 }
 4522 
 4523 /**
 4524  * Round `val` to a multiple of `unit`. `rfunc` specifies the rounding
 4525  * function to use; by default this is `rint`, which uses the current
 4526  * rounding mode.
 4527  */
 4528 Unqual!F quantize(alias rfunc = rint, F)(const F val, const F unit)
 4529 if (is(typeof(rfunc(F.init)) : F) && isFloatingPoint!F)
 4530 {
 4531     typeof(return) ret = val;
 4532     if (unit != 0)
 4533     {
 4534         const scaled = val / unit;
 4535         if (!scaled.isInfinity)
 4536             ret = rfunc(scaled) * unit;
 4537     }
 4538     return ret;
 4539 }
 4540 
 4541 ///
 4542 @safe pure nothrow @nogc unittest
 4543 {
 4544     assert(12345.6789L.quantize(0.01L) == 12345.68L);
 4545     assert(12345.6789L.quantize!floor(0.01L) == 12345.67L);
 4546     assert(12345.6789L.quantize(22.0L) == 12342.0L);
 4547 }
 4548 
 4549 ///
 4550 @safe pure nothrow @nogc unittest
 4551 {
 4552     assert(12345.6789L.quantize(0) == 12345.6789L);
 4553     assert(12345.6789L.quantize(real.infinity).isNaN);
 4554     assert(12345.6789L.quantize(real.nan).isNaN);
 4555     assert(real.infinity.quantize(0.01L) == real.infinity);
 4556     assert(real.infinity.quantize(real.nan).isNaN);
 4557     assert(real.nan.quantize(0.01L).isNaN);
 4558     assert(real.nan.quantize(real.infinity).isNaN);
 4559     assert(real.nan.quantize(real.nan).isNaN);
 4560 }
 4561 
 4562 /**
 4563  * Round `val` to a multiple of `pow(base, exp)`. `rfunc` specifies the
 4564  * rounding function to use; by default this is `rint`, which uses the
 4565  * current rounding mode.
 4566  */
 4567 Unqual!F quantize(real base, alias rfunc = rint, F, E)(const F val, const E exp)
 4568 if (is(typeof(rfunc(F.init)) : F) && isFloatingPoint!F && isIntegral!E)
 4569 {
 4570     // TODO: Compile-time optimization for power-of-two bases?
 4571     return quantize!rfunc(val, pow(cast(F) base, exp));
 4572 }
 4573 
 4574 /// ditto
 4575 Unqual!F quantize(real base, long exp = 1, alias rfunc = rint, F)(const F val)
 4576 if (is(typeof(rfunc(F.init)) : F) && isFloatingPoint!F)
 4577 {
 4578     enum unit = cast(F) pow(base, exp);
 4579     return quantize!rfunc(val, unit);
 4580 }
 4581 
 4582 ///
 4583 @safe pure nothrow @nogc unittest
 4584 {
 4585     assert(12345.6789L.quantize!10(-2) == 12345.68L);
 4586     assert(12345.6789L.quantize!(10, -2) == 12345.68L);
 4587     assert(12345.6789L.quantize!(10, floor)(-2) == 12345.67L);
 4588     assert(12345.6789L.quantize!(10, -2, floor) == 12345.67L);
 4589 
 4590     assert(12345.6789L.quantize!22(1) == 12342.0L);
 4591     assert(12345.6789L.quantize!22 == 12342.0L);
 4592 }
 4593 
 4594 @safe pure nothrow @nogc unittest
 4595 {
 4596     import std.meta : AliasSeq;
 4597 
 4598     static foreach (F; AliasSeq!(real, double, float))
 4599     {{
 4600         const maxL10 = cast(int) F.max.log10.floor;
 4601         const maxR10 = pow(cast(F) 10, maxL10);
 4602         assert((cast(F) 0.9L * maxR10).quantize!10(maxL10) ==  maxR10);
 4603         assert((cast(F)-0.9L * maxR10).quantize!10(maxL10) == -maxR10);
 4604 
 4605         assert(F.max.quantize(F.min_normal) == F.max);
 4606         assert((-F.max).quantize(F.min_normal) == -F.max);
 4607         assert(F.min_normal.quantize(F.max) == 0);
 4608         assert((-F.min_normal).quantize(F.max) == 0);
 4609         assert(F.min_normal.quantize(F.min_normal) == F.min_normal);
 4610         assert((-F.min_normal).quantize(F.min_normal) == -F.min_normal);
 4611     }}
 4612 }
 4613 
 4614 /******************************************
 4615  * Rounds x to the nearest integer value, using the current rounding
 4616  * mode.
 4617  *
 4618  * Unlike the rint functions, nearbyint does not raise the
 4619  * FE_INEXACT exception.
 4620  */
 4621 real nearbyint(real x) @safe pure nothrow @nogc
 4622 {
 4623     return core.stdc.math.nearbyintl(x);
 4624 }
 4625 
 4626 ///
 4627 @safe pure unittest
 4628 {
 4629     assert(nearbyint(0.4) == 0);
 4630     assert(nearbyint(0.5) == 0);
 4631     assert(nearbyint(0.6) == 1);
 4632     assert(nearbyint(100.0) == 100);
 4633 
 4634     assert(isNaN(nearbyint(real.nan)));
 4635     assert(nearbyint(real.infinity) == real.infinity);
 4636     assert(nearbyint(-real.infinity) == -real.infinity);
 4637 }
 4638 
 4639 /**********************************
 4640  * Rounds x to the nearest integer value, using the current rounding
 4641  * mode.
 4642  *
 4643  * If the return value is not equal to x, the FE_INEXACT
 4644  * exception is raised.
 4645  *
 4646  * $(LREF nearbyint) performs the same operation, but does
 4647  * not set the FE_INEXACT exception.
 4648  */
 4649 real rint(real x) @safe pure nothrow @nogc { pragma(inline, true); return core.math.rint(x); }
 4650 //FIXME
 4651 ///ditto
 4652 double rint(double x) @safe pure nothrow @nogc { return rint(cast(real) x); }
 4653 //FIXME
 4654 ///ditto
 4655 float rint(float x) @safe pure nothrow @nogc { return rint(cast(real) x); }
 4656 
 4657 ///
 4658 @safe unittest
 4659 {
 4660     version (IeeeFlagsSupport) resetIeeeFlags();
 4661     assert(rint(0.4) == 0);
 4662     version (IeeeFlagsSupport) assert(ieeeFlags.inexact);
 4663 
 4664     assert(rint(0.5) == 0);
 4665     assert(rint(0.6) == 1);
 4666     assert(rint(100.0) == 100);
 4667 
 4668     assert(isNaN(rint(real.nan)));
 4669     assert(rint(real.infinity) == real.infinity);
 4670     assert(rint(-real.infinity) == -real.infinity);
 4671 }
 4672 
 4673 @safe unittest
 4674 {
 4675     real function(real) print = &rint;
 4676     assert(print != null);
 4677 }
 4678 
 4679 /***************************************
 4680  * Rounds x to the nearest integer value, using the current rounding
 4681  * mode.
 4682  *
 4683  * This is generally the fastest method to convert a floating-point number
 4684  * to an integer. Note that the results from this function
 4685  * depend on the rounding mode, if the fractional part of x is exactly 0.5.
 4686  * If using the default rounding mode (ties round to even integers)
 4687  * lrint(4.5) == 4, lrint(5.5)==6.
 4688  */
 4689 long lrint(real x) @trusted pure nothrow @nogc
 4690 {
 4691     version (InlineAsm_X86_Any)
 4692     {
 4693         version (Win64)
 4694         {
 4695             asm pure nothrow @nogc
 4696             {
 4697                 naked;
 4698                 fld     real ptr [RCX];
 4699                 fistp   qword ptr 8[RSP];
 4700                 mov     RAX,8[RSP];
 4701                 ret;
 4702             }
 4703         }
 4704         else
 4705         {
 4706             long n;
 4707             asm pure nothrow @nogc
 4708             {
 4709                 fld x;
 4710                 fistp n;
 4711             }
 4712             return n;
 4713         }
 4714     }
 4715     else
 4716     {
 4717         alias F = floatTraits!(real);
 4718         static if (F.realFormat == RealFormat.ieeeDouble)
 4719         {
 4720             long result;
 4721 
 4722             // Rounding limit when casting from real(double) to ulong.
 4723             enum real OF = 4.50359962737049600000E15L;
 4724 
 4725             uint* vi = cast(uint*)(&x);
 4726 
 4727             // Find the exponent and sign
 4728             uint msb = vi[MANTISSA_MSB];
 4729             uint lsb = vi[MANTISSA_LSB];
 4730             int exp = ((msb >> 20) & 0x7ff) - 0x3ff;
 4731             const int sign = msb >> 31;
 4732             msb &= 0xfffff;
 4733             msb |= 0x100000;
 4734 
 4735             if (exp < 63)
 4736             {
 4737                 if (exp >= 52)
 4738                     result = (cast(long) msb << (exp - 20)) | (lsb << (exp - 52));
 4739                 else
 4740                 {
 4741                     // Adjust x and check result.
 4742                     const real j = sign ? -OF : OF;
 4743                     x = (j + x) - j;
 4744                     msb = vi[MANTISSA_MSB];
 4745                     lsb = vi[MANTISSA_LSB];
 4746                     exp = ((msb >> 20) & 0x7ff) - 0x3ff;
 4747                     msb &= 0xfffff;
 4748                     msb |= 0x100000;
 4749 
 4750                     if (exp < 0)
 4751                         result = 0;
 4752                     else if (exp < 20)
 4753                         result = cast(long) msb >> (20 - exp);
 4754                     else if (exp == 20)
 4755                         result = cast(long) msb;
 4756                     else
 4757                         result = (cast(long) msb << (exp - 20)) | (lsb >> (52 - exp));
 4758                 }
 4759             }
 4760             else
 4761             {
 4762                 // It is left implementation defined when the number is too large.
 4763                 return cast(long) x;
 4764             }
 4765 
 4766             return sign ? -result : result;
 4767         }
 4768         else static if (F.realFormat == RealFormat.ieeeExtended)
 4769         {
 4770             long result;
 4771 
 4772             // Rounding limit when casting from real(80-bit) to ulong.
 4773             enum real OF = 9.22337203685477580800E18L;
 4774 
 4775             ushort* vu = cast(ushort*)(&x);
 4776             uint* vi = cast(uint*)(&x);
 4777 
 4778             // Find the exponent and sign
 4779             int exp = (vu[F.EXPPOS_SHORT] & 0x7fff) - 0x3fff;
 4780             const int sign = (vu[F.EXPPOS_SHORT] >> 15) & 1;
 4781 
 4782             if (exp < 63)
 4783             {
 4784                 // Adjust x and check result.
 4785                 const real j = sign ? -OF : OF;
 4786                 x = (j + x) - j;
 4787                 exp = (vu[F.EXPPOS_SHORT] & 0x7fff) - 0x3fff;
 4788 
 4789                 version (LittleEndian)
 4790                 {
 4791                     if (exp < 0)
 4792                         result = 0;
 4793                     else if (exp <= 31)
 4794                         result = vi[1] >> (31 - exp);
 4795                     else
 4796                         result = (cast(long) vi[1] << (exp - 31)) | (vi[0] >> (63 - exp));
 4797                 }
 4798                 else
 4799                 {
 4800                     if (exp < 0)
 4801                         result = 0;
 4802                     else if (exp <= 31)
 4803                         result = vi[1] >> (31 - exp);
 4804                     else
 4805                         result = (cast(long) vi[1] << (exp - 31)) | (vi[2] >> (63 - exp));
 4806                 }
 4807             }
 4808             else
 4809             {
 4810                 // It is left implementation defined when the number is too large
 4811                 // to fit in a 64bit long.
 4812                 return cast(long) x;
 4813             }
 4814 
 4815             return sign ? -result : result;
 4816         }
 4817         else static if (F.realFormat == RealFormat.ieeeQuadruple)
 4818         {
 4819             const vu = cast(ushort*)(&x);
 4820 
 4821             // Find the exponent and sign
 4822             const sign = (vu[F.EXPPOS_SHORT] >> 15) & 1;
 4823             if ((vu[F.EXPPOS_SHORT] & F.EXPMASK) - (F.EXPBIAS + 1) > 63)
 4824             {
 4825                 // The result is left implementation defined when the number is
 4826                 // too large to fit in a 64 bit long.
 4827                 return cast(long) x;
 4828             }
 4829 
 4830             // Force rounding of lower bits according to current rounding
 4831             // mode by adding ±2^-112 and subtracting it again.
 4832             enum OF = 5.19229685853482762853049632922009600E33L;
 4833             const j = sign ? -OF : OF;
 4834             x = (j + x) - j;
 4835 
 4836             const exp = (vu[F.EXPPOS_SHORT] & F.EXPMASK) - (F.EXPBIAS + 1);
 4837             const implicitOne = 1UL << 48;
 4838             auto vl = cast(ulong*)(&x);
 4839             vl[MANTISSA_MSB] &= implicitOne - 1;
 4840             vl[MANTISSA_MSB] |= implicitOne;
 4841 
 4842             long result;
 4843 
 4844             if (exp < 0)
 4845                 result = 0;
 4846             else if (exp <= 48)
 4847                 result = vl[MANTISSA_MSB] >> (48 - exp);
 4848             else
 4849                 result = (vl[MANTISSA_MSB] << (exp - 48)) | (vl[MANTISSA_LSB] >> (112 - exp));
 4850 
 4851             return sign ? -result : result;
 4852         }
 4853         else
 4854         {
 4855             static assert(false, "real type not supported by lrint()");
 4856         }
 4857     }
 4858 }
 4859 
 4860 ///
 4861 @safe pure nothrow @nogc unittest
 4862 {
 4863     assert(lrint(4.5) == 4);
 4864     assert(lrint(5.5) == 6);
 4865     assert(lrint(-4.5) == -4);
 4866     assert(lrint(-5.5) == -6);
 4867 
 4868     assert(lrint(int.max - 0.5) == 2147483646L);
 4869     assert(lrint(int.max + 0.5) == 2147483648L);
 4870     assert(lrint(int.min - 0.5) == -2147483648L);
 4871     assert(lrint(int.min + 0.5) == -2147483648L);
 4872 }
 4873 
 4874 static if (real.mant_dig >= long.sizeof * 8)
 4875 {
 4876     @safe pure nothrow @nogc unittest
 4877     {
 4878         assert(lrint(long.max - 1.5L) == long.max - 1);
 4879         assert(lrint(long.max - 0.5L) == long.max - 1);
 4880         assert(lrint(long.min + 0.5L) == long.min);
 4881         assert(lrint(long.min + 1.5L) == long.min + 2);
 4882     }
 4883 }
 4884 
 4885 /*******************************************
 4886  * Return the value of x rounded to the nearest integer.
 4887  * If the fractional part of x is exactly 0.5, the return value is
 4888  * rounded away from zero.
 4889  *
 4890  * Returns:
 4891  *     A `real`.
 4892  */
 4893 auto round(real x) @trusted nothrow @nogc
 4894 {
 4895     version (CRuntime_Microsoft)
 4896     {
 4897         auto old = FloatingPointControl.getControlState();
 4898         FloatingPointControl.setControlState(
 4899             (old & (-1 - FloatingPointControl.roundingMask)) | FloatingPointControl.roundToZero
 4900         );
 4901         x = rint((x >= 0) ? x + 0.5 : x - 0.5);
 4902         FloatingPointControl.setControlState(old);
 4903         return x;
 4904     }
 4905     else
 4906         return core.stdc.math.roundl(x);
 4907 }
 4908 
 4909 ///
 4910 @safe nothrow @nogc unittest
 4911 {
 4912     assert(round(4.5) == 5);
 4913     assert(round(5.4) == 5);
 4914     assert(round(-4.5) == -5);
 4915     assert(round(-5.1) == -5);
 4916 }
 4917 
 4918 // assure purity on Posix
 4919 version (Posix)
 4920 {
 4921     @safe pure nothrow @nogc unittest
 4922     {
 4923         assert(round(4.5) == 5);
 4924     }
 4925 }
 4926 
 4927 /**********************************************
 4928  * Return the value of x rounded to the nearest integer.
 4929  *
 4930  * If the fractional part of x is exactly 0.5, the return value is rounded
 4931  * away from zero.
 4932  *
 4933  * $(BLUE This function is not implemented for Digital Mars C runtime.)
 4934  */
 4935 long lround(real x) @trusted nothrow @nogc
 4936 {
 4937     version (CRuntime_DigitalMars)
 4938         assert(0, "lround not implemented");
 4939     else
 4940         return core.stdc.math.llroundl(x);
 4941 }
 4942 
 4943 ///
 4944 @safe nothrow @nogc unittest
 4945 {
 4946     version (CRuntime_DigitalMars) {}
 4947     else
 4948     {
 4949         assert(lround(0.49) == 0);
 4950         assert(lround(0.5) == 1);
 4951         assert(lround(1.5) == 2);
 4952     }
 4953 }
 4954 
 4955 /**
 4956  Returns the integer portion of x, dropping the fractional portion.
 4957  This is also known as "chop" rounding.
 4958  `pure` on all platforms.
 4959  */
 4960 real trunc(real x) @trusted nothrow @nogc pure
 4961 {
 4962     version (Win64_DMD_InlineAsm)
 4963     {
 4964         asm pure nothrow @nogc
 4965         {
 4966             naked                       ;
 4967             fld     real ptr [RCX]      ;
 4968             fstcw   8[RSP]              ;
 4969             mov     AL,9[RSP]           ;
 4970             mov     DL,AL               ;
 4971             and     AL,0xC3             ;
 4972             or      AL,0x0C             ; // round to 0
 4973             mov     9[RSP],AL           ;
 4974             fldcw   8[RSP]              ;
 4975             frndint                     ;
 4976             mov     9[RSP],DL           ;
 4977             fldcw   8[RSP]              ;
 4978             ret                         ;
 4979         }
 4980     }
 4981     else version (MSVC_InlineAsm)
 4982     {
 4983         short cw;
 4984         asm pure nothrow @nogc
 4985         {
 4986             fld     x                   ;
 4987             fstcw   cw                  ;
 4988             mov     AL,byte ptr cw+1    ;
 4989             mov     DL,AL               ;
 4990             and     AL,0xC3             ;
 4991             or      AL,0x0C             ; // round to 0
 4992             mov     byte ptr cw+1,AL    ;
 4993             fldcw   cw                  ;
 4994             frndint                     ;
 4995             mov     byte ptr cw+1,DL    ;
 4996             fldcw   cw                  ;
 4997         }
 4998     }
 4999     else
 5000         return core.stdc.math.truncl(x);
 5001 }
 5002 
 5003 ///
 5004 @safe pure unittest
 5005 {
 5006     assert(trunc(0.01) == 0);
 5007     assert(trunc(0.49) == 0);
 5008     assert(trunc(0.5) == 0);
 5009     assert(trunc(1.5) == 1);
 5010 }
 5011 
 5012 /****************************************************
 5013  * Calculate the remainder x REM y, following IEC 60559.
 5014  *
 5015  * REM is the value of x - y * n, where n is the integer nearest the exact
 5016  * value of x / y.
 5017  * If |n - x / y| == 0.5, n is even.
 5018  * If the result is zero, it has the same sign as x.
 5019  * Otherwise, the sign of the result is the sign of x / y.
 5020  * Precision mode has no effect on the remainder functions.
 5021  *
 5022  * remquo returns `n` in the parameter `n`.
 5023  *
 5024  * $(TABLE_SV
 5025  *  $(TR $(TH x)               $(TH y)            $(TH remainder(x, y)) $(TH n)   $(TH invalid?))
 5026  *  $(TR $(TD $(PLUSMN)0.0)    $(TD not 0.0)      $(TD $(PLUSMN)0.0)    $(TD 0.0) $(TD no))
 5027  *  $(TR $(TD $(PLUSMNINF))    $(TD anything)     $(TD -$(NAN))         $(TD ?)   $(TD yes))
 5028  *  $(TR $(TD anything)        $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)$(NAN)) $(TD ?)   $(TD yes))
 5029  *  $(TR $(TD != $(PLUSMNINF)) $(TD $(PLUSMNINF)) $(TD x)               $(TD ?)   $(TD no))
 5030  * )
 5031  */
 5032 real remainder(real x, real y) @trusted nothrow @nogc
 5033 {
 5034     return core.stdc.math.remainderl(x, y);
 5035 }
 5036 
 5037 /// ditto
 5038 real remquo(real x, real y, out int n) @trusted nothrow @nogc  /// ditto
 5039 {
 5040     return core.stdc.math.remquol(x, y, &n);
 5041 }
 5042 
 5043 ///
 5044 @safe @nogc nothrow unittest
 5045 {
 5046     assert(remainder(5.1, 3.0).feqrel(-0.9) > 16);
 5047     assert(remainder(-5.1, 3.0).feqrel(0.9) > 16);
 5048     assert(remainder(0.0, 3.0) == 0.0);
 5049 
 5050     assert(isNaN(remainder(1.0, 0.0)));
 5051     assert(isNaN(remainder(-1.0, 0.0)));
 5052 }
 5053 
 5054 ///
 5055 @safe @nogc nothrow unittest
 5056 {
 5057     int n;
 5058 
 5059     assert(remquo(5.1, 3.0, n).feqrel(-0.9) > 16 && n == 2);
 5060     assert(remquo(-5.1, 3.0, n).feqrel(0.9) > 16 && n == -2);
 5061     assert(remquo(0.0, 3.0, n) == 0.0 && n == 0);
 5062 }
 5063 
 5064 
 5065 version (IeeeFlagsSupport)
 5066 {
 5067 
 5068 /** IEEE exception status flags ('sticky bits')
 5069 
 5070  These flags indicate that an exceptional floating-point condition has occurred.
 5071  They indicate that a NaN or an infinity has been generated, that a result
 5072  is inexact, or that a signalling NaN has been encountered. If floating-point
 5073  exceptions are enabled (unmasked), a hardware exception will be generated
 5074  instead of setting these flags.
 5075  */
 5076 struct IeeeFlags
 5077 {
 5078 nothrow @nogc:
 5079 
 5080 private:
 5081     // The x87 FPU status register is 16 bits.
 5082     // The Pentium SSE2 status register is 32 bits.
 5083     // The ARM and PowerPC FPSCR is a 32-bit register.
 5084     // The SPARC FSR is a 32bit register (64 bits for SPARC 7 & 8, but high bits are uninteresting).
 5085     // The RISC-V (32 & 64 bit) fcsr is 32-bit register.
 5086     uint flags;
 5087 
 5088     version (CRuntime_Microsoft)
 5089     {
 5090         // Microsoft uses hardware-incompatible custom constants in fenv.h (core.stdc.fenv).
 5091         // Applies to both x87 status word (16 bits) and SSE2 status word(32 bits).
 5092         enum : int
 5093         {
 5094             INEXACT_MASK   = 0x20,
 5095             UNDERFLOW_MASK = 0x10,
 5096             OVERFLOW_MASK  = 0x08,
 5097             DIVBYZERO_MASK = 0x04,
 5098             INVALID_MASK   = 0x01,
 5099 
 5100             EXCEPTIONS_MASK = 0b11_1111
 5101         }
 5102         // Don't bother about subnormals, they are not supported on most CPUs.
 5103         //  SUBNORMAL_MASK = 0x02;
 5104     }
 5105     else
 5106     {
 5107         enum : int
 5108         {
 5109             INEXACT_MASK    = core.stdc.fenv.FE_INEXACT,
 5110             UNDERFLOW_MASK  = core.stdc.fenv.FE_UNDERFLOW,
 5111             OVERFLOW_MASK   = core.stdc.fenv.FE_OVERFLOW,
 5112             DIVBYZERO_MASK  = core.stdc.fenv.FE_DIVBYZERO,
 5113             INVALID_MASK    = core.stdc.fenv.FE_INVALID,
 5114             EXCEPTIONS_MASK = core.stdc.fenv.FE_ALL_EXCEPT,
 5115         }
 5116     }
 5117 
 5118     static uint getIeeeFlags() @trusted pure
 5119     {
 5120         version (InlineAsm_X86_Any)
 5121         {
 5122             ushort sw;
 5123             asm pure nothrow @nogc { fstsw sw; }
 5124 
 5125             // OR the result with the SSE2 status register (MXCSR).
 5126             if (haveSSE)
 5127             {
 5128                 uint m