Logo Search packages:      
Sourcecode: ldc version File versions

math.d

// Written in the D programming language

/**
 * Macros:
 *    WIKI = Phobos/StdMath
 *
 *    TABLE_SV = <table border=1 cellpadding=4 cellspacing=0>
 *          <caption>Special Values</caption>
 *          $0</table>
 *    SVH = $(TR $(TH $1) $(TH $2))
 *    SV  = $(TR $(TD $1) $(TD $2))
 *
 *    NAN = $(RED NAN)
 *    SUP = <span style="vertical-align:super;font-size:smaller">$0</span>
 *    GAMMA =  &#915;
 *    INTEGRAL = &#8747;
 *    INTEGRATE = $(BIG &#8747;<sub>$(SMALL $1)</sub><sup>$2</sup>)
 *    POWER = $1<sup>$2</sup>
 *    BIGSUM = $(BIG &Sigma; <sup>$2</sup><sub>$(SMALL $1)</sub>)
 *    CHOOSE = $(BIG &#40;) <sup>$(SMALL $1)</sup><sub>$(SMALL $2)</sub> $(BIG &#41;)
 *    PLUSMN = &plusmn;
 *    INFIN = &infin;
 *    PI = &pi;
 *    LT = &lt;
 *    GT = &gt;
 */

/*
 * Author:
 *    Walter Bright
 * Copyright:
 *    Copyright (c) 2001-2005 by Digital Mars,
 *    All Rights Reserved,
 *    www.digitalmars.com
 *  Copyright (c) 2007 by Tomas Lindquist Olsen
 * License:
 *  This software is provided 'as-is', without any express or implied
 *  warranty. In no event will the authors be held liable for any damages
 *  arising from the use of this software.
 *
 *  Permission is granted to anyone to use this software for any purpose,
 *  including commercial applications, and to alter it and redistribute it
 *  freely, subject to the following restrictions:
 *
 *  <ul>
 *  <li> The origin of this software must not be misrepresented; you must not
 *       claim that you wrote the original software. If you use this software
 *       in a product, an acknowledgment in the product documentation would be
 *       appreciated but is not required.
 *  </li>
 *  <li> Altered source versions must be plainly marked as such, and must not
 *       be misrepresented as being the original software.
 *  </li>
 *  <li> This notice may not be removed or altered from any source
 *       distribution.
 *  </li>
 *  </ul>
 */


00061 module std.math;

//debug=math;           // uncomment to turn on debugging printf's

private import std.c.stdio;
private import std.c.math;

class NotImplemented : Error
{
    this(string msg)
    {
      super(msg ~ "not implemented");
    }
}

const real E =          2.7182818284590452354L;  /** e */
00077 const real LOG2T =      0x1.a934f0979a3715fcp+1; /** log<sub>2</sub>10 */ // 3.32193 fldl2t
00078 const real LOG2E =      0x1.71547652b82fe178p+0; /** log<sub>2</sub>e */ // 1.4427 fldl2e
00079 const real LOG2 = 0x1.34413509f79fef32p-2; /** log<sub>10</sub>2 */ // 0.30103 fldlg2
00080 const real LOG10E =     0.43429448190325182765;  /** log<sub>10</sub>e */
00081 const real LN2 =  0x1.62e42fefa39ef358p-1; /** ln 2 */      // 0.693147 fldln2
00082 const real LN10 = 2.30258509299404568402;  /** ln 10 */
00083 const real PI =         0x1.921fb54442d1846ap+1; /** $(PI) */ // 3.14159 fldpi
00084 const real PI_2 = 1.57079632679489661923;  /** $(PI) / 2 */
00085 const real PI_4 = 0.78539816339744830962;  /** $(PI) / 4 */
00086 const real M_1_PI =     0.31830988618379067154;  /** 1 / $(PI) */
00087 const real M_2_PI =     0.63661977236758134308;  /** 2 / $(PI) */
00088 const real M_2_SQRTPI = 1.12837916709551257390;  /** 2 / &radic;$(PI) */
00089 const real SQRT2 =      1.41421356237309504880;  /** &radic;2 */
00090 const real SQRT1_2 =    0.70710678118654752440;  /** &radic;&frac12; */

/*
      Octal versions:
      PI/64800    0.00001 45530 36176 77347 02143 15351 61441 26767
      PI/180            0.01073 72152 11224 72344 25603 54276 63351 22056
      PI/8        0.31103 75524 21026 43021 51423 06305 05600 67016
      SQRT(1/PI)  0.44067 27240 41233 33210 65616 51051 77327 77303
      2/PI        0.50574 60333 44710 40522 47741 16537 21752 32335
      PI/4        0.62207 73250 42055 06043 23046 14612 13401 56034
      SQRT(2/PI)  0.63041 05147 52066 24106 41762 63612 00272 56161

      PI          3.11037 55242 10264 30215 14230 63050 56006 70163
      LOG2        0.23210 11520 47674 77674 61076 11263 26013 37111
 */


/***********************************
 * Calculates the absolute value
 *
 * For complex numbers, abs(z) = sqrt( $(POWER z.re, 2) + $(POWER z.im, 2) )
 * = hypot(z.re, z.im).
 */
00113 real abs(real x)
{
    return fabs(x);
}

/** ditto */
00119 long abs(long x)
{
    return x>=0 ? x : -x;
}

/** ditto */
00125 int abs(int x)
{
    return x>=0 ? x : -x;
}

/** ditto */
00131 real abs(creal z)
{
    return hypot(z.re, z.im);
}

/** ditto */
00137 real abs(ireal y)
{
    return fabs(y.im);
}


unittest
{
    assert(isPosZero(abs(-0.0L)));
    assert(isnan(abs(real.nan)));
    assert(abs(-real.infinity) == real.infinity);
    assert(abs(-3.2Li) == 3.2L);
    assert(abs(71.6Li) == 71.6L);
    assert(abs(-56) == 56);
    assert(abs(2321312L)  == 2321312L);
    assert(abs(-1+1i) == sqrt(2.0));
}

/***********************************
 * Complex conjugate
 *
 *  conj(x + iy) = x - iy
 *
 * Note that z * conj(z) = $(POWER z.re, 2) - $(POWER z.im, 2)
 * is always a real number
 */
creal conj(creal z)
{
    return z.re - z.im*1i;
}

/** ditto */
ireal conj(ireal y)
{
    return -y;
}

unittest
{
    assert(conj(7 + 3i) == 7-3i);
    ireal z = -3.2Li;
    assert(conj(z) == -z);
}

/***********************************
 * Returns cosine of x. x is in radians.
 *
 *    $(TABLE_SV
 *    $(TR $(TH x)                 $(TH cos(x)) $(TH invalid?))
 *    $(TR $(TD $(NAN))            $(TD $(NAN)) $(TD yes)   )
 *    $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(NAN)) $(TD yes)   )
 *    )
 * Bugs:
 *    Results are undefined if |x| >= $(POWER 2,64).
 */

pragma(intrinsic, "llvm.cos.f32")
float cos(float x);

pragma(intrinsic, "llvm.cos.f64")
double cos(double x); // ditto

pragma(intrinsic, "llvm.cos.f80")
real cos(real x); /// ditto


/***********************************
 * Returns sine of x. x is in radians.
 *
 *    $(TABLE_SV
 *    $(TR $(TH x)               $(TH sin(x))      $(TH invalid?))
 *    $(TR $(TD $(NAN))          $(TD $(NAN))      $(TD yes))
 *    $(TR $(TD $(PLUSMN)0.0)    $(TD $(PLUSMN)0.0) $(TD no))
 *    $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(NAN))      $(TD yes))
 *    )
 * Bugs:
 *    Results are undefined if |x| >= $(POWER 2,64).
 */

pragma(intrinsic, "llvm.sin.f32")
float sin(float x);

pragma(intrinsic, "llvm.sin.f64")
double sin(double x); // ditto

pragma(intrinsic, "llvm.sin.f80")
real sin(real x); /// ditto


/****************************************************************************
 * Returns tangent of x. x is in radians.
 *
 *    $(TABLE_SV
 *    $(TR $(TH x)               $(TH tan(x))       $(TH invalid?))
 *    $(TR $(TD $(NAN))          $(TD $(NAN))       $(TD yes))
 *    $(TR $(TD $(PLUSMN)0.0)    $(TD $(PLUSMN)0.0) $(TD no))
 *    $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(NAN))     $(TD yes))
 *    )
 */

00237 version(D_InlineAsm_X86)
real tan(real x)
{

    asm
    {
      fld   x[EBP]                  ; // load theta
      fxam                    ; // test for oddball values
      fstsw AX                ;
      sahf                    ;
      jc    trigerr                 ; // x is NAN, infinity, or empty
                                // 387's can handle denormals
SC18: fptan                   ;
      fstp  ST(0)             ; // dump X, which is always 1
      fstsw AX                ;
      sahf                    ;
      jnp   Lret              ; // C2 = 1 (x is out of range)

      // Do argument reduction to bring x into range
      fldpi                   ;
      fxch                    ;
SC17: fprem1                        ;
      fstsw AX                ;
      sahf                    ;
      jp    SC17              ;
      fstp  ST(1)             ; // remove pi from stack
      jmp   SC18              ;

trigerr:
      jnp   Lret              ; // if theta is NAN, return theta
      fstp  ST(0)             ; // dump theta
    }
    return real.nan;

Lret:
    ;
}
else
{
real tan(real x) { return std.c.math.atan(x); }
}


unittest
{
    static real vals[][2] =   // angle,tan
    [
          [   0,   0],
          [   .5,  .5463024898],
          [   1,   1.557407725],
          [   1.5, 14.10141995],
          [   2,  -2.185039863],
          [   2.5,-.7470222972],
          [   3,  -.1425465431],
          [   3.5, .3745856402],
          [   4,   1.157821282],
          [   4.5, 4.637332055],
          [   5,  -3.380515006],
          [   5.5,-.9955840522],
          [   6,  -.2910061914],
          [   6.5, .2202772003],
          [   10,  .6483608275],

          // special angles
          [   PI_4,     1],
          //[     PI_2, real.infinity],
          [   3*PI_4,   -1],
          [   PI, 0],
          [   5*PI_4,   1],
          //[     3*PI_2,     -real.infinity],
          [   7*PI_4,   -1],
          [   2*PI,     0],

          // overflow
          [   real.infinity,  real.nan],
          [   real.nan, real.nan],
          //[   1e+100, real.nan],
    ];
    int i;

    for (i = 0; i < vals.length; i++)
    {
      real x = vals[i][0];
      real r = vals[i][1];
      real t = tan(x);

      //printf("tan(%Lg) = %Lg, should be %Lg\n", x, t, r);
      assert(mfeq(r, t, .0000001));

      x = -x;
      r = -r;
      t = tan(x);
      //printf("tan(%Lg) = %Lg, should be %Lg\n", x, t, r);
      assert(mfeq(r, t, .0000001));
    }
}

/***************
 * Calculates the arc cosine of x,
 * returning a value ranging from -$(PI)/2 to $(PI)/2.
 *
 *    $(TABLE_SV
 *      $(TR $(TH x)         $(TH acos(x)) $(TH invalid?))
 *      $(TR $(TD $(GT)1.0)  $(TD $(NAN))  $(TD yes))
 *      $(TR $(TD $(LT)-1.0) $(TD $(NAN))  $(TD yes))
 *      $(TR $(TD $(NAN))    $(TD $(NAN))  $(TD yes))
 *      )
 */
real acos(real x)       { return std.c.math.acos(x); }

/***************
 * Calculates the arc sine of x,
 * returning a value ranging from -$(PI)/2 to $(PI)/2.
 *
 *    $(TABLE_SV
 *    $(TR $(TH x)            $(TH asin(x))      $(TH invalid?))
 *    $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no))
 *    $(TR $(TD $(GT)1.0)     $(TD $(NAN))       $(TD yes))
 *    $(TR $(TD $(LT)-1.0)    $(TD $(NAN))       $(TD yes))
 *       )
 */
real asin(real x)       { return std.c.math.asin(x); }

/***************
 * Calculates the arc tangent of x,
 * returning a value ranging from -$(PI)/2 to $(PI)/2.
 *
 *    $(TABLE_SV
 *    $(TR $(TH x)                 $(TH atan(x))      $(TH invalid?))
 *    $(TR $(TD $(PLUSMN)0.0)      $(TD $(PLUSMN)0.0) $(TD no))
 *    $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(NAN))       $(TD yes))
 *       )
 */
real atan(real x)       { return std.c.math.atan(x); }

/***************
 * Calculates the arc tangent of y / x,
 * returning a value ranging from -$(PI)/2 to $(PI)/2.
 *
 *      $(TABLE_SV
 *      $(TR $(TH y)                 $(TH x)            $(TH atan(y, x)))
 *      $(TR $(TD $(NAN))            $(TD anything)     $(TD $(NAN)) )
 *      $(TR $(TD anything)          $(TD $(NAN))       $(TD $(NAN)) )
 *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(GT)0.0)     $(TD $(PLUSMN)0.0) )
 *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) )
 *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(LT)0.0)     $(TD $(PLUSMN)$(PI)))
 *      $(TR $(TD $(PLUSMN)0.0)      $(TD -0.0)         $(TD $(PLUSMN)$(PI)))
 *      $(TR $(TD $(GT)0.0)          $(TD $(PLUSMN)0.0) $(TD $(PI)/2) )
 *      $(TR $(TD $(LT)0.0)          $(TD $(PLUSMN)0.0) $(TD $(PI)/2))
 *      $(TR $(TD $(GT)0.0)          $(TD $(INFIN))     $(TD $(PLUSMN)0.0) )
 *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD anything)     $(TD $(PLUSMN)$(PI)/2))
 *      $(TR $(TD $(GT)0.0)          $(TD -$(INFIN))    $(TD $(PLUSMN)$(PI)) )
 *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(INFIN))     $(TD $(PLUSMN)$(PI)/4))
 *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD -$(INFIN))    $(TD $(PLUSMN)3$(PI)/4))
 *      )
 */
real atan2(real y, real x)      { return std.c.math.atan2(y,x); }

/***********************************
 * Calculates the hyperbolic cosine of x.
 *
 *    $(TABLE_SV
 *    $(TR $(TH x)                 $(TH cosh(x))      $(TH invalid?))
 *    $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)0.0) $(TD no) )
 *      )
 */
real cosh(real x)       { return std.c.math.cosh(x); }

/***********************************
 * Calculates the hyperbolic sine of x.
 *
 *    $(TABLE_SV
 *    $(TR $(TH x)                 $(TH sinh(x))           $(TH invalid?))
 *    $(TR $(TD $(PLUSMN)0.0)      $(TD $(PLUSMN)0.0)      $(TD no))
 *    $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)$(INFIN)) $(TD no))
 *      )
 */
real sinh(real x)       { return std.c.math.sinh(x); }

/***********************************
 * Calculates the hyperbolic tangent of x.
 *
 *    $(TABLE_SV
 *    $(TR $(TH x)                 $(TH tanh(x))      $(TH invalid?))
 *    $(TR $(TD $(PLUSMN)0.0)      $(TD $(PLUSMN)0.0) $(TD no) )
 *    $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)1.0) $(TD no))
 *      )
 */
real tanh(real x)       { return std.c.math.tanh(x); }

//real acosh(real x)          { return std.c.math.acoshl(x); }
//real asinh(real x)          { return std.c.math.asinhl(x); }
//real atanh(real x)          { return std.c.math.atanhl(x); }

/***********************************
 * Calculates the inverse hyperbolic cosine of x.
 *
 *  Mathematically, acosh(x) = log(x + sqrt( x*x - 1))
 *
 * $(TABLE_DOMRG
 *  $(DOMAIN 1..$(INFIN))
 *  $(RANGE  1..log(real.max), $(INFIN)) )
 *    $(TABLE_SV
 *    $(SVH  x,     acosh(x) )
 *    $(SV  $(NAN), $(NAN) )
 *    $(SV  <1,     $(NAN) )
 *    $(SV  1,      0       )
 *    $(SV  +$(INFIN),+$(INFIN))
 *  )
 */   
real acosh(real x)
{
    if (x > 1/real.epsilon)
      return LN2 + log(x);
    else
      return log(x + sqrt(x*x - 1));
}

unittest
{
    assert(isnan(acosh(0.9)));
    assert(isnan(acosh(real.nan)));
    assert(acosh(1)==0.0);
    assert(acosh(real.infinity) == real.infinity);
}

/***********************************
 * Calculates the inverse hyperbolic sine of x.
 *
 *  Mathematically,
 *  ---------------
 *  asinh(x) =  log( x + sqrt( x*x + 1 )) // if x >= +0
 *  asinh(x) = -log(-x + sqrt( x*x + 1 )) // if x <= -0
 *  -------------
 *
 *    $(TABLE_SV
 *    $(SVH x,                asinh(x)       )
 *    $(SV  $(NAN),           $(NAN)         )
 *    $(SV  $(PLUSMN)0,       $(PLUSMN)0      )
 *    $(SV  $(PLUSMN)$(INFIN),$(PLUSMN)$(INFIN))
 *  )
 */
real asinh(real x)
{   
    if (fabs(x) > 1 / real.epsilon)   // beyond this point, x*x + 1 == x*x
      return copysign(LN2 + log(fabs(x)), x);
    else
    {
      // sqrt(x*x + 1) ==  1 + x * x / ( 1 + sqrt(x*x + 1) )
      return copysign(log1p(fabs(x) + x*x / (1 + sqrt(x*x + 1)) ), x);
    }
}

unittest
{
    assert(isPosZero(asinh(0.0)));
    assert(isNegZero(asinh(-0.0)));
    assert(asinh(real.infinity) == real.infinity);
    assert(asinh(-real.infinity) == -real.infinity);
    assert(isnan(asinh(real.nan)));
}

/***********************************
 * Calculates the inverse hyperbolic tangent of x,
 * returning a value from ranging from -1 to 1.
 *  
 * Mathematically, atanh(x) = log( (1+x)/(1-x) ) / 2
 *  
 *
 * $(TABLE_DOMRG
 *  $(DOMAIN -$(INFIN)..$(INFIN))
 *  $(RANGE  -1..1) )
 *    $(TABLE_SV
 *    $(SVH  x,     acosh(x) )
 *    $(SV  $(NAN), $(NAN) )
 *    $(SV  $(PLUSMN)0, $(PLUSMN)0)
 *    $(SV  -$(INFIN), -0)
 *  )
 */   
real atanh(real x)
{
    // log( (1+x)/(1-x) ) == log ( 1 + (2*x)/(1-x) )
    return  0.5 * log1p( 2 * x / (1 - x) );
}

unittest
{
    assert(isPosZero(atanh(0.0)));
    assert(isNegZero(atanh(-0.0)));
    assert(isnan(atanh(real.nan)));
    assert(isnan(atanh(-real.infinity))); 
}

/*****************************************
 * Returns x rounded to a long value using the current rounding mode.
 * If the integer value of x is
 * greater than long.max, the result is
 * indeterminate.
 */
long rndtol(real x);    /* intrinsic */


/*****************************************
 * Returns x rounded to a long value using the FE_TONEAREST rounding mode.
 * If the integer value of x is
 * greater than long.max, the result is
 * indeterminate.
 */
extern (C) real rndtonl(real x);

/***************************************
 * Compute square root of x.
 *
 *    $(TABLE_SV
 *    $(TR $(TH x)         $(TH sqrt(x))   $(TH invalid?))
 *    $(TR $(TD -0.0)      $(TD -0.0)      $(TD no))
 *    $(TR $(TD $(LT)0.0)  $(TD $(NAN))    $(TD yes))
 *    $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no))
 *    )
 */

pragma(intrinsic, "llvm.sqrt.f32")
float sqrt(float x);    /* intrinsic */

pragma(intrinsic, "llvm.sqrt.f64")
double sqrt(double x);  /* intrinsic */   /// ditto

pragma(intrinsic, "llvm.sqrt.f80")
real sqrt(real x);      /* intrinsic */ /// ditto

00567 creal sqrt(creal z)
{
    creal c;
    real x,y,w,r;

    if (z == 0)
    {
      c = 0 + 0i;
    }
    else
    { real z_re = z.re;
      real z_im = z.im;

      x = fabs(z_re);
      y = fabs(z_im);
      if (x >= y)
      {
          r = y / x;
          w = sqrt(x) * sqrt(0.5 * (1 + sqrt(1 + r * r)));
      }
      else
      {
          r = x / y;
          w = sqrt(y) * sqrt(0.5 * (r + sqrt(1 + r * r)));
      }

      if (z_re >= 0)
      {
          c = w + (z_im / (w + w)) * 1.0i;
      }
      else
      {
          if (z_im < 0)
            w = -w;
          c = z_im / (w + w) + w * 1.0i;
      }
    }
    return c;
}

/**********************
 * Calculates e$(SUP x).
 *
 *    $(TABLE_SV
 *    $(TR $(TH x)         $(TH exp(x)))
 *    $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) )
 *    $(TR $(TD -$(INFIN)) $(TD +0.0) )
 *    )
 */
real exp(real x)        { return std.c.math.exp(x); }

/**********************
 * Calculates 2$(SUP x).
 *
 *    $(TABLE_SV
 *    $(TR $(TH x)         $(TH exp2(x)))
 *    $(TR $(TD +$(INFIN)) $(TD +$(INFIN)))
 *    $(TR $(TD -$(INFIN)) $(TD +0.0))
 *    )
 */
real exp2(real x)       { return std.c.math.exp2(x); }

/******************************************
 * Calculates the value of the natural logarithm base (e)
 * raised to the power of x, minus 1.
 *
 * For very small x, expm1(x) is more accurate 
 * than exp(x)-1. 
 *
 *    $(TABLE_SV
 *    $(TR $(TH x)            $(TH e$(SUP x)-1))
 *    $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0))
 *    $(TR $(TD +$(INFIN))    $(TD +$(INFIN)))
 *    $(TR $(TD -$(INFIN))    $(TD -1.0))
 *    )
 */

real expm1(real x)            { return std.c.math.expm1(x); }


/*********************************************************************
 * Separate floating point value into significand and exponent.
 *
 * Returns:
 *    Calculate and return <i>x</i> and exp such that
 *    value =<i>x</i>*2$(SUP exp) and
 *    .5 $(LT)= |<i>x</i>| $(LT) 1.0<br>
 *    <i>x</i> has same sign as value.
 *
 *    $(TABLE_SV
 *    $(TR $(TH value)           $(TH returns)         $(TH exp))
 *    $(TR $(TD $(PLUSMN)0.0)    $(TD $(PLUSMN)0.0)    $(TD 0))
 *    $(TR $(TD +$(INFIN))       $(TD +$(INFIN))       $(TD int.max))
 *    $(TR $(TD -$(INFIN))       $(TD -$(INFIN))       $(TD int.min))
 *    $(TR $(TD $(PLUSMN)$(NAN)) $(TD $(PLUSMN)$(NAN)) $(TD int.min))
 *    )
 */


real frexp(real value, out int exp)
{
    ushort* vu = cast(ushort*)&value;
    long* vl = cast(long*)&value;
    uint ex;

    // If exponent is non-zero
    ex = vu[4] & 0x7FFF;
    if (ex)
    {
      if (ex == 0x7FFF)
      {   // infinity or NaN
          if (*vl &  0x7FFFFFFFFFFFFFFF)  // if NaN
          { *vl |= 0xC000000000000000;    // convert $(NAN)S to $(NAN)Q
            exp = int.min;
          }
          else if (vu[4] & 0x8000)
          { // negative infinity
            exp = int.min;
          }
          else
          { // positive infinity
            exp = int.max;
          }
      }
      else
      {
          exp = ex - 0x3FFE;
          vu[4] = cast(ushort)((0x8000 & vu[4]) | 0x3FFE);
      }
    }
    else if (!*vl)
    {
      // value is +-0.0
      exp = 0;
    }
    else
    { // denormal
      int i = -0x3FFD;

      do
      {
          i--;
          *vl <<= 1;
      } while (*vl > 0);
      exp = i;
        vu[4] = cast(ushort)((0x8000 & vu[4]) | 0x3FFE);
    }
    return value;
}


unittest
{
    static real vals[][3] =   // x,frexp,exp
    [
      [0.0, 0.0,  0],
      [-0.0,      -0.0, 0],
      [1.0, .5,   1],
      [-1.0,      -.5,  1],
      [2.0, .5,   2],
      [155.67e20, 0x1.A5F1C2EB3FE4Fp-1,   74],  // normal
      [1.0e-320,  0.98829225,       -1063],
      [real.min,  .5,         -16381],
      [real.min/2.0L,   .5,         -16382],    // denormal

      [real.infinity,real.infinity,int.max],
      [-real.infinity,-real.infinity,int.min],
      [real.nan,real.nan,int.min],
      [-real.nan,-real.nan,int.min],

      // Don't really support signalling nan's in D
      //[real.nans,real.nan,int.min],
      //[-real.nans,-real.nan,int.min],
    ];
    int i;

    for (i = 0; i < vals.length; i++)
    {
      real x = vals[i][0];
      real e = vals[i][1];
      int exp = cast(int)vals[i][2];
      int eptr;
      real v = frexp(x, eptr);

      //printf("frexp(%Lg) = %.8Lg, should be %.8Lg, eptr = %d, should be %d\n", x, v, e, eptr, exp);
      assert(mfeq(e, v, .0000001));
      assert(exp == eptr);
    }
}


/******************************************
 * Extracts the exponent of x as a signed integral value.
 *
 * If x is not a special value, the result is the same as
 * <tt>cast(int)logb(x)</tt>.
 *
 *    $(TABLE_SV
 *    $(TR $(TH x)                $(TH ilogb(x))     $(TH Range error?))
 *    $(TR $(TD 0)                 $(TD FP_ILOGB0)   $(TD yes))
 *    $(TR $(TD $(PLUSMN)$(INFIN)) $(TD int.max)     $(TD no))
 *    $(TR $(TD $(NAN))            $(TD FP_ILOGBNAN) $(TD no))
 *    )
 */
int ilogb(real x)       { return std.c.math.ilogb(x); }

alias std.c.math.FP_ILOGB0   FP_ILOGB0;
alias std.c.math.FP_ILOGBNAN FP_ILOGBNAN;


/*******************************************
 * Compute n * 2$(SUP exp)
 * References: frexp
 */

real ldexp(real n, int exp)   { return std.c.math.ldexp(n, exp); }

/**************************************
 * Calculate the natural logarithm of x.
 *
 *    $(TABLE_SV
 *    $(TR $(TH x)            $(TH log(x))    $(TH divide by 0?) $(TH invalid?))
 *    $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes)          $(TD no))
 *    $(TR $(TD $(LT)0.0)     $(TD $(NAN))    $(TD no)           $(TD yes))
 *    $(TR $(TD +$(INFIN))    $(TD +$(INFIN)) $(TD no)           $(TD no))
 *    )
 */

real log(real x)        { return std.c.math.log(x); }

/**************************************
 * Calculate the base-10 logarithm of x.
 *
 *    $(TABLE_SV
 *    $(TR $(TH x)            $(TH log10(x))  $(TH divide by 0?) $(TH invalid?))
 *    $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes)          $(TD no))
 *    $(TR $(TD $(LT)0.0)     $(TD $(NAN))    $(TD no)           $(TD yes))
 *    $(TR $(TD +$(INFIN))    $(TD +$(INFIN)) $(TD no)           $(TD no))
 *    )
 */

real log10(real x)            { return std.c.math.log10(x); }

/******************************************
 *    Calculates the natural logarithm of 1 + x.
 *
 *    For very small x, log1p(x) will be more accurate than 
 *    log(1 + x). 
 *
 *    $(TABLE_SV
 *    $(TR $(TH x)            $(TH log1p(x))     $(TH divide by 0?) $(TH invalid?))
 *    $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no)           $(TD no))
 *    $(TR $(TD -1.0)         $(TD -$(INFIN))    $(TD yes)          $(TD no))
 *    $(TR $(TD $(LT)-1.0)    $(TD $(NAN))       $(TD no)           $(TD yes))
 *    $(TR $(TD +$(INFIN))    $(TD -$(INFIN))    $(TD no)           $(TD no))
 *    )
 */

real log1p(real x)            { return std.c.math.log1p(x); }

/***************************************
 * Calculates the base-2 logarithm of x:
 * log<sub>2</sub>x
 *
 *    $(TABLE_SV
 *    $(TR $(TH x)              $(TH log2(x))   $(TH divide by 0?) $(TH invalid?))
 *    $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes)          $(TD no) )
 *    $(TR $(TD $(LT)0.0)     $(TD $(NAN))    $(TD no)           $(TD yes) )
 *    $(TR $(TD +$(INFIN))    $(TD +$(INFIN)) $(TD no)           $(TD no) )
 *    )
 */
real log2(real x)       { return std.c.math.log2(x); }

/*****************************************
 * Extracts the exponent of x as a signed integral value.
 *
 * If x is subnormal, it is treated as if it were normalized.
 * For a positive, finite x: 
 *
 * -----
 * 1 <= $(I x) * FLT_RADIX$(SUP -logb(x)) < FLT_RADIX 
 * -----
 *
 *    $(TABLE_SV
 *    $(TR $(TH x)                 $(TH logb(x))   $(TH divide by 0?) )
 *    $(TR $(TD $(PLUSMN)$(INFIN)) $(TD +$(INFIN)) $(TD no))
 *    $(TR $(TD $(PLUSMN)0.0)      $(TD -$(INFIN)) $(TD yes) )
 *    )
 */
real logb(real x)       { return std.c.math.logb(x); }

/************************************
 * Calculates the remainder from the calculation x/y.
 * Returns:
 * The value of x - i * y, where i is the number of times that y can 
 * be completely subtracted from x. The result has the same sign as x. 
 *
 *    $(TABLE_SV
 *    $(TR $(TH x)                   $(TH y)                 $(TH modf(x, y))   $(TH invalid?))
 *    $(TR $(TD $(PLUSMN)0.0)        $(TD no)t 0.0           $(TD $(PLUSMN)0.0) $(TD no))
 *    $(TR $(TD $(PLUSMN)$(INFIN))   $(TD anything)          $(TD $(NAN))       $(TD yes))
 *    $(TR $(TD anything)            $(TD $(PLUSMN)0.0)      $(TD $(NAN))       $(TD yes))
 *    $(TR $(TD !=$(PLUSMN)$(INFIN)) $(TD $(PLUSMN)$(INFIN)) $(TD x)            $(TD no))
 *    )
 */
real modf(real x, inout real y)
{
double Y = y;
auto tmp = std.c.math.modf(x,&Y);
y = Y;
return tmp;
}

/*************************************
 * Efficiently calculates x * 2$(SUP n).
 *
 * scalbn handles underflow and overflow in 
 * the same fashion as the basic arithmetic operators. 
 *
 *    $(TABLE_SV
 *    $(TR $(TH x)                 $(TH scalb(x)))
 *    $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)$(INFIN)) )
 *    $(TR $(TD $(PLUSMN)0.0)      $(TD $(PLUSMN)0.0) )
 *    )
 */
real scalbn(real x, int n)
{
    version (linux)
      return std.c.math.scalbn(x, n);
    else
      throw new NotImplemented("scalbn");
}

/***************
 * Calculates the cube root x.
 *
 *    $(TABLE_SV
 *    $(TR $(TH $(I x))      $(TH cbrt(x))           $(TH invalid?))
 *    $(TR $(TD $(PLUSMN)0.0)      $(TD $(PLUSMN)0.0)      $(TD no) )
 *    $(TR $(TD $(NAN))      $(TD $(NAN))            $(TD yes) )
 *    $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)$(INFIN)) $(TD no) )
 *    )
 */
real cbrt(real x)       { return std.c.math.cbrt(x); }


/*******************************
 * Returns |x|
 *
 *    $(TABLE_SV
 *    $(TR $(TH x)                 $(TH fabs(x)))
 *    $(TR $(TD $(PLUSMN)0.0)      $(TD +0.0) )
 *    $(TR $(TD $(PLUSMN)$(INFIN)) $(TD +$(INFIN)) )
 *    )
 */
real fabs(real x) { return std.c.math.fabs(x); }


/***********************************************************************
 * Calculates the length of the 
 * hypotenuse of a right-angled triangle with sides of length x and y. 
 * The hypotenuse is the value of the square root of 
 * the sums of the squares of x and y:
 *
 *    sqrt(x&sup2; + y&sup2;)
 *
 * Note that hypot(x, y), hypot(y, x) and
 * hypot(x, -y) are equivalent.
 *
 *    $(TABLE_SV
 *    $(TR $(TH x)                 $(TH y)            $(TH hypot(x, y)) $(TH invalid?))
 *    $(TR $(TD x)                 $(TD $(PLUSMN)0.0) $(TD |x|)         $(TD no))
 *    $(TR $(TD $(PLUSMN)$(INFIN)) $(TD y)            $(TD +$(INFIN))   $(TD no))
 *    $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(NAN))       $(TD +$(INFIN))   $(TD no))
 *    )
 */

real hypot(real x, real y)
{
    /*
     * This is based on code from:
     * Cephes Math Library Release 2.1:  January, 1989
     * Copyright 1984, 1987, 1989 by Stephen L. Moshier
     * Direct inquiries to 30 Frost Street, Cambridge, MA 02140
     */

    const int PRECL = 32;
    const int MAXEXPL = real.max_exp; //16384;
    const int MINEXPL = real.min_exp; //-16384;

    real xx, yy, b, re, im;
    int ex, ey, e;

    // Note, hypot(INFINITY, NAN) = INFINITY.
    if (isinf(x) || isinf(y))
      return real.infinity;

    if (isnan(x))
      return x;
    if (isnan(y))
      return y;

    re = fabs(x);
    im = fabs(y);

    if (re == 0.0)
      return im;
    if (im == 0.0)
      return re;

    // Get the exponents of the numbers
    xx = frexp(re, ex);
    yy = frexp(im, ey);

    // Check if one number is tiny compared to the other
    e = ex - ey;
    if (e > PRECL)
      return re;
    if (e < -PRECL)
      return im;

    // Find approximate exponent e of the geometric mean.
    e = (ex + ey) >> 1;

    // Rescale so mean is about 1
    xx = ldexp(re, -e);
    yy = ldexp(im, -e);

    // Hypotenuse of the right triangle
    b = sqrt(xx * xx  +  yy * yy);

    // Compute the exponent of the answer.
    yy = frexp(b, ey);
    ey = e + ey;

    // Check it for overflow and underflow.
    if (ey > MAXEXPL + 2)
    {
      //return __matherr(_OVERFLOW, INFINITY, x, y, "hypotl");
      return real.infinity;
    }
    if (ey < MINEXPL - 2)
      return 0.0;

    // Undo the scaling
    b = ldexp(b, e);
    return b;
}

unittest
{
    static real vals[][3] =   // x,y,hypot
    [
      [     0,    0,    0],
      [     0,    -0,   0],
      [     3,    4,    5],
      [     -300, -400, 500],
      [     real.min, real.min, 4.75473e-4932L],
      [     real.max/2, real.max/2, 0x1.6a09e667f3bcc908p+16383L /*8.41267e+4931L*/],
      [     real.infinity, real.nan, real.infinity],
      [     real.nan, real.nan, real.nan],
    ];
    int i;

    for (i = 0; i < vals.length; i++)
    {
      real x = vals[i][0];
      real y = vals[i][1];
      real z = vals[i][2];
      real h = hypot(x, y);

      //printf("hypot(%Lg, %Lg) = %Lg, should be %Lg\n", x, y, h, z);
      //if (!mfeq(z, h, .0000001))
          //printf("%La\n", h);
      assert(mfeq(z, h, .0000001));
    }
}

/**********************************
 * Returns the error function of x.
 *
 * <img src="erf.gif" alt="error function">
 */
real erf(real x)        { return std.c.math.erf(x); }

/**********************************
 * Returns the complementary error function of x, which is 1 - erf(x).
 *
 * <img src="erfc.gif" alt="complementary error function">
 */
real erfc(real x)       { return std.c.math.erfc(x); }

/***********************************
 * Natural logarithm of gamma function.
 *
 * Returns the base e (2.718...) logarithm of the absolute
 * value of the gamma function of the argument.
 *
 * For reals, lgamma is equivalent to log(fabs(gamma(x))).
 *
 *    $(TABLE_SV
 *    $(TR $(TH x)                 $(TH lgamma(x)) $(TH invalid?))
 *    $(TR $(TD $(NAN))            $(TD $(NAN))    $(TD yes))
 *    $(TR $(TD integer <= 0)      $(TD +$(INFIN)) $(TD yes))
 *    $(TR $(TD $(PLUSMN)$(INFIN)) $(TD +$(INFIN)) $(TD no))
 *    )
 */
/* Documentation prepared by Don Clugston */
real lgamma(real x)
{
    return std.c.math.lgamma(x);

    // Use etc.gamma.lgamma for those C systems that are missing it
}

/***********************************
 *  The Gamma function, $(GAMMA)(x)
 *
 *  $(GAMMA)(x) is a generalisation of the factorial function
 *  to real and complex numbers.
 *  Like x!, $(GAMMA)(x+1) = x*$(GAMMA)(x).
 *
 *  Mathematically, if z.re > 0 then
 *   $(GAMMA)(z) =<big>$(INTEGRAL)<sub><small>0</small></sub><sup>$(INFIN)</sup></big>t<sup>z-1</sup>e<sup>-t</sup>dt
 *
 *    $(TABLE_SV
 *    $(TR $(TH x)              $(TH $(GAMMA)(x))       $(TH invalid?))
 *    $(TR $(TD $(NAN))         $(TD $(NAN))            $(TD yes))
 *    $(TR $(TD $(PLUSMN)0.0)   $(TD $(PLUSMN)$(INFIN)) $(TD yes))
 *    $(TR $(TD integer $(GT)0) $(TD (x-1)!)            $(TD no))
 *    $(TR $(TD integer $(LT)0) $(TD $(NAN))            $(TD yes))
 *    $(TR $(TD +$(INFIN))      $(TD +$(INFIN))         $(TD no))
 *    $(TR $(TD -$(INFIN))      $(TD $(NAN))            $(TD yes))
 *    )
 *
 *  References:
 *    $(LINK http://en.wikipedia.org/wiki/Gamma_function),
 *    $(LINK http://www.netlib.org/cephes/ldoubdoc.html#gamma)
 */
/* Documentation prepared by Don Clugston */
real tgamma(real x)
{
    return std.c.math.tgamma(x);

    // Use etc.gamma.tgamma for those C systems that are missing it
}

/**************************************
 * Returns the value of x rounded upward to the next integer
 * (toward positive infinity).
 */
real ceil(real x)       { return std.c.math.ceil(x); }

/**************************************
 * Returns the value of x rounded downward to the next integer
 * (toward negative infinity).
 */
real floor(real x)            { return std.c.math.floor(x); }

/******************************************
 * Rounds x to the nearest integer value, using the current rounding 
 * mode.
 *
 * Unlike the rint functions, nearbyint does not raise the 
 * FE_INEXACT exception. 
 */
real nearbyint(real x) { return std.c.math.nearbyint(x); }

/**********************************
 * Rounds x to the nearest integer value, using the current rounding
 * mode.
 * If the return value is not equal to x, the FE_INEXACT
 * exception is raised.
 * <b>nearbyint</b> performs
 * the same operation, but does not set the FE_INEXACT exception.
 */
real rint(real x) { return std.c.math.rint(x); }

/***************************************
 * Rounds x to the nearest integer value, using the current rounding
 * mode.
 */
long lrint(real x)
{
    version (linux)
      return std.c.math.llrint(x);
    else
      throw new NotImplemented("lrint");
}

/*******************************************
 * Return the value of x rounded to the nearest integer.
 * If the fractional part of x is exactly 0.5, the return value is rounded to
 * the even integer. 
 */
real round(real x) { return std.c.math.round(x); }

/**********************************************
 * Return the value of x rounded to the nearest integer.
 *
 * If the fractional part of x is exactly 0.5, the return value is rounded
 * away from zero.
 */
long lround(real x)
{
    version (linux)
      return std.c.math.llround(x);
    else
      throw new NotImplemented("lround");
}

/****************************************************
 * Returns the integer portion of x, dropping the fractional portion. 
 *
 * This is also known as "chop" rounding. 
 */
real trunc(real x) { return std.c.math.trunc(x); }

/****************************************************
 * Calculate the remainder x REM y, following IEC 60559.
 *
 * REM is the value of x - y * n, where n is the integer nearest the exact 
 * value of x / y.
 * If |n - x / y| == 0.5, n is even.
 * If the result is zero, it has the same sign as x.
 * Otherwise, the sign of the result is the sign of x / y.
 * Precision mode has no effect on the remainder functions.
 *
 * remquo returns n in the parameter n.
 *
 *    $(TABLE_SV
 *    $(TR $(TH x)                    $(TH y)                 $(TH remainder(x, y)) $(TH n)   $(TH invalid?))
 *    $(TR $(TD $(PLUSMN)0.0)         $(TD no)t 0.0           $(TD $(PLUSMN)0.0)    $(TD 0.0) $(TD no))
 *    $(TR $(TD $(PLUSMN)$(INFIN))    $(TD anything)          $(TD $(NAN))          $(TD ?)   $(TD yes))
 *    $(TR $(TD anything)             $(TD $(PLUSMN)0.0)      $(TD $(NAN))          $(TD ?)   $(TD yes))
 *    $(TR $(TD != $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)$(INFIN)) $(TD x)               $(TD ?)   $(TD no))
 *    )
 */
real remainder(real x, real y) { return std.c.math.remainder(x, y); }

real remquo(real x, real y, out int n)    /// ditto
{
    version (linux)
      return std.c.math.remquo(x, y, &n);
    else
      throw new NotImplemented("remquo");
}

/*********************************
 * Returns !=0 if e is a NaN.
 */

int isnan(real e)
{
    ushort* pe = cast(ushort *)&e;
    ulong*  ps = cast(ulong *)&e;

    return (pe[4] & 0x7FFF) == 0x7FFF &&
          *ps & 0x7FFFFFFFFFFFFFFF;
}

unittest
{
    assert(isnan(float.nan));
    assert(isnan(-double.nan));
    assert(isnan(real.nan));

    assert(!isnan(53.6));
    assert(!isnan(float.infinity));
}

/*********************************
 * Returns !=0 if e is finite.
 */

int isfinite(real e)
{
    ushort* pe = cast(ushort *)&e;

    return (pe[4] & 0x7FFF) != 0x7FFF;
}

unittest
{
    assert(isfinite(1.23));
    assert(!isfinite(double.infinity));
    assert(!isfinite(float.nan));
}


/*********************************
 * Returns !=0 if x is normalized.
 */

/* Need one for each format because subnormal floats might
 * be converted to normal reals.
 */

int isnormal(float x)
{
    uint *p = cast(uint *)&x;
    uint e;

    e = *p & 0x7F800000;
    //printf("e = x%x, *p = x%x\n", e, *p);
    return e && e != 0x7F800000;
}

/// ditto

int isnormal(double d)
{
    uint *p = cast(uint *)&d;
    uint e;

    e = p[1] & 0x7FF00000;
    return e && e != 0x7FF00000;
}

/// ditto

int isnormal(real e)
{
    ushort* pe = cast(ushort *)&e;
    long*   ps = cast(long *)&e;

    return (pe[4] & 0x7FFF) != 0x7FFF && *ps < 0;
}

unittest
{
    float f = 3;
    double d = 500;
    real e = 10e+48;

    assert(isnormal(f));
    assert(isnormal(d));
    assert(isnormal(e));
}

/*********************************
 * Is number subnormal? (Also called "denormal".)
 * Subnormals have a 0 exponent and a 0 most significant mantissa bit.
 */

/* Need one for each format because subnormal floats might
 * be converted to normal reals.
 */

int issubnormal(float f)
{
    uint *p = cast(uint *)&f;

    //printf("*p = x%x\n", *p);
    return (*p & 0x7F800000) == 0 && *p & 0x007FFFFF;
}

unittest
{
    float f = 3.0;

    for (f = 1.0; !issubnormal(f); f /= 2)
      assert(f != 0);
}

/// ditto

int issubnormal(double d)
{
    uint *p = cast(uint *)&d;

    return (p[1] & 0x7FF00000) == 0 && (p[0] || p[1] & 0x000FFFFF);
}

unittest
{
    double f;

    for (f = 1; !issubnormal(f); f /= 2)
      assert(f != 0);
}

/// ditto

int issubnormal(real e)
{
    ushort* pe = cast(ushort *)&e;
    long*   ps = cast(long *)&e;

    return (pe[4] & 0x7FFF) == 0 && *ps > 0;
}

unittest
{
    real f;

    for (f = 1; !issubnormal(f); f /= 2)
      assert(f != 0);
}

/*********************************
 * Return !=0 if e is $(PLUSMN)$(INFIN).
 */

int isinf(real e)
{
    ushort* pe = cast(ushort *)&e;
    ulong*  ps = cast(ulong *)&e;

    return (pe[4] & 0x7FFF) == 0x7FFF &&
          *ps == 0x8000000000000000;
}

unittest
{
    assert(isinf(float.infinity));
    assert(!isinf(float.nan));
    assert(isinf(double.infinity));
    assert(isinf(-real.infinity));

    assert(isinf(-1.0 / 0.0));
}

/*********************************
 * Return 1 if sign bit of e is set, 0 if not.
 */

int signbit(real e)
{
    ubyte* pe = cast(ubyte *)&e;

//printf("e = %Lg\n", e);
    return (pe[9] & 0x80) != 0;
}

unittest
{
    debug (math) printf("math.signbit.unittest\n");
    assert(!signbit(float.nan));
    assert(signbit(-float.nan));
    assert(!signbit(168.1234));
    assert(signbit(-168.1234));
    assert(!signbit(0.0));
    assert(signbit(-0.0));
}

/*********************************
 * Return a value composed of to with from's sign bit.
 */

real copysign(real to, real from)
{
    ubyte* pto   = cast(ubyte *)&to;
    ubyte* pfrom = cast(ubyte *)&from;

    pto[9] &= 0x7F;
    pto[9] |= pfrom[9] & 0x80;

    return to;
}

unittest
{
    real e;

    e = copysign(21, 23.8);
    assert(e == 21);

    e = copysign(-21, 23.8);
    assert(e == 21);

    e = copysign(21, -23.8);
    assert(e == -21);

    e = copysign(-21, -23.8);
    assert(e == -21);

    e = copysign(real.nan, -23.8);
    assert(isnan(e) && signbit(e));
}

/******************************************
 * Creates a quiet NAN with the information from tagp[] embedded in it.
 */
real nan(char[] tagp) { return std.c.math.nan((tagp~\0).ptr); }

/******************************************
 * Calculates the next representable value after x in the direction of y. 
 *
 * If y $(GT) x, the result will be the next largest floating-point value;
 * if y $(LT) x, the result will be the next smallest value.
 * If x == y, the result is y.
 * The FE_INEXACT and FE_OVERFLOW exceptions will be raised if x is finite and
 * the function result is infinite. The FE_INEXACT and FE_UNDERFLOW 
 * exceptions will be raised if the function value is subnormal, and x is 
 * not equal to y. 
 */
real nextafter(real x, real y)
{
    version (linux)
      return std.c.math.nextafterl(x, y);
    else
      throw new NotImplemented("nextafter");
}

//real nexttoward(real x, real y) { return std.c.math.nexttowardl(x, y); }

/*******************************************
 * Returns the positive difference between x and y.
 * Returns:
 *    $(TABLE_SV
 *    $(TR $(TH x, y)       $(TH fdim(x, y)))
 *    $(TR $(TD x $(GT) y)  $(TD x - y))
 *    $(TR $(TD x $(LT)= y) $(TD +0.0))
 *    )
 */
real fdim(real x, real y) { return (x > y) ? x - y : +0.0; }

/****************************************
 * Returns the larger of x and y.
 */
real fmax(real x, real y) { return x > y ? x : y; }

/****************************************
 * Returns the smaller of x and y.
 */
real fmin(real x, real y) { return x < y ? x : y; }

/**************************************
 * Returns (x * y) + z, rounding only once according to the
 * current rounding mode.
 */
real fma(real x, real y, real z) { return (x * y) + z; }

/*******************************************************************
 * Fast integral powers.
 */

pragma(intrinsic, "llvm.powi.f32")
{
float pow(float x, uint n);
/// ditto
float pow(float x, int n);
}

pragma(intrinsic, "llvm.powi.f64")
{
/// ditto
double pow(double x, uint n);
/// ditto
double pow(double x, int n);
}

pragma(intrinsic, "llvm.powi.f80")
{
/// ditto
real pow(real x, uint n);
/// ditto
real pow(real x, int n);
}

/+
real pow(real x, uint n);
{
    real p;

    switch (n)
    {
      case 0:
          p = 1.0;
          break;

      case 1:
          p = x;
          break;

      case 2:
          p = x * x;
          break;

      default:
          p = 1.0;
          while (1)
          {
            if (n & 1)
                p *= x;
            n >>= 1;
            if (!n)
                break;
            x *= x;
          }
          break;
    }
    return p;
}

/// ditto
real pow(real x, int n);
{
    if (n < 0)
      return pow(x, cast(real)n);
    else
      return pow(x, cast(uint)n);
}
+/

/*********************************************
 * Calculates x$(SUP y).
 *
 * $(TABLE_SV
 * $(TR
 * $(TH x) $(TH y) $(TH pow(x, y)) $(TH div 0) $(TH invalid?))
 * $(TR
 * $(TD anything)       $(TD $(PLUSMN)0.0)      $(TD 1.0)   $(TD no)    $(TD no) )
 * $(TR
 * $(TD |x| $(GT) 1)    $(TD +$(INFIN))   $(TD +$(INFIN))   $(TD no)    $(TD no) )
 * $(TR
 * $(TD |x| $(LT) 1)    $(TD +$(INFIN))   $(TD +0.0)  $(TD no)    $(TD no) )
 * $(TR
 * $(TD |x| $(GT) 1)    $(TD -$(INFIN))   $(TD +0.0)  $(TD no)    $(TD no) )
 * $(TR
 * $(TD |x| $(LT) 1)    $(TD -$(INFIN))   $(TD +$(INFIN))   $(TD no)    $(TD no) )
 * $(TR
 * $(TD +$(INFIN))      $(TD $(GT) 0.0)   $(TD +$(INFIN))   $(TD no)    $(TD no) )
 * $(TR
 * $(TD +$(INFIN))      $(TD $(LT) 0.0)   $(TD +0.0)  $(TD no)    $(TD no) )
 * $(TR
 * $(TD -$(INFIN))      $(TD odd integer $(GT) 0.0)   $(TD -$(INFIN))   $(TD no)    $(TD no) )
 * $(TR
 * $(TD -$(INFIN))      $(TD $(GT) 0.0, not odd integer) $(TD +$(INFIN))      $(TD no)    $(TD no))
 * $(TR
 * $(TD -$(INFIN))      $(TD odd integer $(LT) 0.0)   $(TD -0.0)  $(TD no)    $(TD no) )
 * $(TR
 * $(TD -$(INFIN))      $(TD $(LT) 0.0, not odd integer) $(TD +0.0)     $(TD no)    $(TD no) )
 * $(TR
 * $(TD $(PLUSMN)1.0)   $(TD $(PLUSMN)$(INFIN))       $(TD $(NAN))      $(TD no)    $(TD yes) )
 * $(TR
 * $(TD $(LT) 0.0)      $(TD finite, nonintegral)     $(TD $(NAN))      $(TD no)    $(TD yes))
 * $(TR
 * $(TD $(PLUSMN)0.0)   $(TD odd integer $(LT) 0.0)   $(TD $(PLUSMN)$(INFIN)) $(TD yes)   $(TD no) )
 * $(TR
 * $(TD $(PLUSMN)0.0)   $(TD $(LT) 0.0, not odd integer) $(TD +$(INFIN))      $(TD yes)   $(TD no))
 * $(TR
 * $(TD $(PLUSMN)0.0)   $(TD odd integer $(GT) 0.0)   $(TD $(PLUSMN)0.0) $(TD no)   $(TD no) )
 * $(TR
 * $(TD $(PLUSMN)0.0)   $(TD $(GT) 0.0, not odd integer) $(TD +0.0)     $(TD no)    $(TD no) )
 * )
 */

pragma(intrinsic, "llvm.pow.f32")
float pow(float x, float y);

pragma(intrinsic, "llvm.pow.f64")
double pow(double x, double y);

pragma(intrinsic, "llvm.pow.f80")
real pow(real x, real y);

/+
real pow(real x, real y);
{
    version (linux) // C pow() often does not handle special values correctly
    {
      if (isnan(y))
          return y;

      if (y == 0)
          return 1;           // even if x is $(NAN)
      if (isnan(x) && y != 0)
          return x;
      if (isinf(y))
      {
          if (fabs(x) > 1)
          {
            if (signbit(y))
                return +0.0;
            else
                return real.infinity;
          }
          else if (fabs(x) == 1)
          {
            return real.nan;
          }
          else // < 1
          {
            if (signbit(y))
                return real.infinity;
            else
                return +0.0;
          }
      }
      if (isinf(x))
      {
          if (signbit(x))
          {   long i;

            i = cast(long)y;
            if (y > 0)
            {
                if (i == y && i & 1)
                  return -real.infinity;
                else
                  return real.infinity;
            }
            else if (y < 0)
            {
                if (i == y && i & 1)
                  return -0.0;
                else
                  return +0.0;
            }
          }
          else
          {
            if (y > 0)
                return real.infinity;
            else if (y < 0)
                return +0.0;
          }
      }

      if (x == 0.0)
      {
          if (signbit(x))
          {   long i;

            i = cast(long)y;
            if (y > 0)
            {
                if (i == y && i & 1)
                  return -0.0;
                else
                  return +0.0;
            }
            else if (y < 0)
            {
                if (i == y && i & 1)
                  return -real.infinity;
                else
                  return real.infinity;
            }
          }
          else
          {
            if (y > 0)
                return +0.0;
            else if (y < 0)
                return real.infinity;
          }
      }
    }
    return std.c.math.powl(x, y);
}
+/

unittest
{
    real x = 46;

    assert(pow(x,0) == 1.0);
    assert(pow(x,1) == x);
    assert(pow(x,2) == x * x);
    assert(pow(x,3) == x * x * x);
    assert(pow(x,8) == (x * x) * (x * x) * (x * x) * (x * x));
}

/****************************************
 * Simple function to compare two floating point values
 * to a specified precision.
 * Returns:
 *    1     match
 *    0     nomatch
 */

private int mfeq(real x, real y, real precision)
{
    if (x == y)
      return 1;
    if (isnan(x))
      return isnan(y);
    if (isnan(y))
      return 0;
    return fabs(x - y) <= precision;
}

// Returns true if x is +0.0 (This function is used in unit tests)
bool isPosZero(real x)
{
    return (x == 0) && (signbit(x) == 0);
}

// Returns true if x is -0.0 (This function is used in unit tests)
bool isNegZero(real x)
{
    return (x == 0) && signbit(x);
}

/**************************************
 * To what precision is x equal to y?
 *
 * Returns: the number of mantissa bits which are equal in x and y.
 * eg, 0x1.F8p+60 and 0x1.F1p+60 are equal to 5 bits of precision.
 *
 *    $(TABLE_SV
 *    $(TR $(TH x)      $(TH y)          $(TH feqrel(x, y)))
 *    $(TR $(TD x)      $(TD x)          $(TD real.mant_dig))
 *    $(TR $(TD x)      $(TD $(GT)= 2*x) $(TD 0))
 *    $(TR $(TD x)      $(TD $(LT)= x/2) $(TD 0))
 *    $(TR $(TD $(NAN)) $(TD any)        $(TD 0))
 *    $(TR $(TD any)    $(TD $(NAN))     $(TD 0))
 *    )
 */

int feqrel(real x, real y)
{
    /* Public Domain. Author: Don Clugston, 18 Aug 2005.
     */

    if (x == y)
      return real.mant_dig; // ensure diff!=0, cope with INF.

    real diff = fabs(x - y);

    ushort *pa = cast(ushort *)(&x);
    ushort *pb = cast(ushort *)(&y);
    ushort *pd = cast(ushort *)(&diff);

    // The difference in abs(exponent) between x or y and abs(x-y)
    // is equal to the number of mantissa bits of x which are
    // equal to y. If negative, x and y have different exponents.
    // If positive, x and y are equal to 'bitsdiff' bits.
    // AND with 0x7FFF to form the absolute value.
    // To avoid out-by-1 errors, we subtract 1 so it rounds down
    // if the exponents were different. This means 'bitsdiff' is
    // always 1 lower than we want, except that if bitsdiff==0,
    // they could have 0 or 1 bits in common.
    int bitsdiff = ( ((pa[4]&0x7FFF) + (pb[4]&0x7FFF)-1)>>1) - pd[4];

    if (pd[4] == 0)
    { // Difference is denormal
      // For denormals, we need to add the number of zeros that
      // lie at the start of diff's mantissa.
      // We do this by multiplying by 2^real.mant_dig
      diff *= 0x1p+63;
      return bitsdiff + real.mant_dig - pd[4];
    }

    if (bitsdiff > 0)
      return bitsdiff + 1; // add the 1 we subtracted before

    // Avoid out-by-1 errors when factor is almost 2.
    return (bitsdiff == 0) ? (pa[4] == pb[4]) : 0;
}

unittest
{
   // Exact equality
   assert(feqrel(real.max,real.max)==real.mant_dig);
   assert(feqrel(0,0)==real.mant_dig);
   assert(feqrel(7.1824,7.1824)==real.mant_dig);
   assert(feqrel(real.infinity,real.infinity)==real.mant_dig);

   // a few bits away from exact equality
   real w=1;
   for (int i=1; i<real.mant_dig-1; ++i) {
      assert(feqrel(1+w*real.epsilon,1)==real.mant_dig-i);
      assert(feqrel(1-w*real.epsilon,1)==real.mant_dig-i);
      assert(feqrel(1,1+(w-1)*real.epsilon)==real.mant_dig-i+1);
      w*=2;
   }
   assert(feqrel(1.5+real.epsilon,1.5)==real.mant_dig-1);
   assert(feqrel(1.5-real.epsilon,1.5)==real.mant_dig-1);
   assert(feqrel(1.5-real.epsilon,1.5+real.epsilon)==real.mant_dig-2);

   // Numbers that are close
   assert(feqrel(0x1.Bp+84, 0x1.B8p+84)==5);
   assert(feqrel(0x1.8p+10, 0x1.Cp+10)==2);
   assert(feqrel(1.5*(1-real.epsilon), 1)==2);
   assert(feqrel(1.5, 1)==1);
   assert(feqrel(2*(1-real.epsilon), 1)==1);

   // Factors of 2
   assert(feqrel(real.max,real.infinity)==0);
   assert(feqrel(2*(1-real.epsilon), 1)==1);
   assert(feqrel(1, 2)==0);
   assert(feqrel(4, 1)==0);

   // Extreme inequality
   assert(feqrel(real.nan,real.nan)==0);
   assert(feqrel(0,-real.nan)==0);
   assert(feqrel(real.nan,real.infinity)==0);
   assert(feqrel(real.infinity,-real.infinity)==0);
   assert(feqrel(-real.max,real.infinity)==0);
   assert(feqrel(real.max,-real.max)==0);
}


/***********************************
 * Evaluate polynomial A(x) = a<sub>0</sub> + a<sub>1</sub>x + a<sub>2</sub>x&sup2; + a<sub>3</sub>x&sup3; ...
 *
 * Uses Horner's rule A(x) = a<sub>0</sub> + x(a<sub>1</sub> + x(a<sub>2</sub> + x(a<sub>3</sub> + ...)))
 * Params:
 *    A =   array of coefficients a<sub>0</sub>, a<sub>1</sub>, etc.
 */ 
real poly(real x, real[] A)
in
{
    assert(A.length > 0);
}
body
{
    version (D_InlineAsm_X86)
    {
      version (Windows)
      {
          asm     // assembler by W. Bright
          {
            // EDX = (A.length - 1) * real.sizeof
            mov     ECX,A[EBP]            ; // ECX = A.length
            dec     ECX             ;
            lea     EDX,[ECX][ECX*8]      ;
            add     EDX,ECX               ;
            add     EDX,A+4[EBP]          ;
            fld     real ptr [EDX]        ; // ST0 = coeff[ECX]
            jecxz   return_ST       ;
            fld     x[EBP]                ; // ST0 = x
            fxch    ST(1)                 ; // ST1 = x, ST0 = r
            align   4               ;
      L2:     fmul    ST,ST(1)            ; // r *= x
            fld     real ptr -10[EDX]     ;
            sub     EDX,10                ; // deg--
            faddp   ST(1),ST        ;
            dec     ECX             ;
            jne     L2              ;
            fxch    ST(1)                 ; // ST1 = r, ST0 = x
            fstp    ST(0)                 ; // dump x
            align   4               ;
      return_ST:                    ;
            ;
          }
      }
      else
      {
          asm     // assembler by W. Bright
          {
            // EDX = (A.length - 1) * real.sizeof
            mov     ECX,A[EBP]            ; // ECX = A.length
            dec     ECX             ;
            lea     EDX,[ECX*8]           ;
            lea   EDX,[EDX][ECX*4]  ;
            add     EDX,A+4[EBP]          ;
            fld     real ptr [EDX]        ; // ST0 = coeff[ECX]
            jecxz   return_ST       ;
            fld     x[EBP]                ; // ST0 = x
            fxch    ST(1)                 ; // ST1 = x, ST0 = r
            align   4               ;
      L2:     fmul    ST,ST(1)            ; // r *= x
            fld     real ptr -12[EDX]     ;
            sub     EDX,12                ; // deg--
            faddp   ST(1),ST        ;
            dec     ECX             ;
            jne     L2              ;
            fxch    ST(1)                 ; // ST1 = r, ST0 = x
            fstp    ST(0)                 ; // dump x
            align   4               ;
      return_ST:                    ;
            ;
          }
      }
    }
    else
    {
      int i = A.length - 1;
      real r = A[i];
      while (--i >= 0)
      {
          r *= x;
          r += A[i];
      }
      return r;
    }
}

unittest
{
    debug (math) printf("math.poly.unittest\n");
    real x = 3.1;
    static real pp[] = [56.1, 32.7, 6];

    assert( poly(x, pp) == (56.1L + (32.7L + 6L * x) * x) );
}



Generated by  Doxygen 1.6.0   Back to index