libphobos/libdruntime/core/math.d

   1 // Written in the D programming language.
   2
   3 /**
   4  * Builtin mathematical intrinsics
   5  *
   6  * Source: $(DRUNTIMESRC core/_math.d)
   7  * Macros:
   8  *      TABLE_SV = <table border="1" cellpadding="4" cellspacing="0">
   9  *              <caption>Special Values</caption>
  10  *              $0</table>
  11  *
  12  *      NAN = $(RED NAN)
  13  *      SUP = <span style="vertical-align:super;font-size:smaller">$0</span>
  14  *      POWER = $1<sup>$2</sup>
  15  *      PLUSMN = &plusmn;
  16  *      INFIN = &infin;
  17  *      PLUSMNINF = &plusmn;&infin;
  18  *      LT = &lt;
  19  *      GT = &gt;
  20  *
  21  * Copyright: Copyright Digital Mars 2000 - 2011.
  22  * License:   $(HTTP www.boost.org/LICENSE_1_0.txt, Boost License 1.0).
  23  * Authors:   $(HTTP digitalmars.com, Walter Bright),
  24  *                        Don Clugston
  25  */
  26 module core.math;
  27
  28 public:
  29 @nogc:
  30 nothrow:
  31 @safe:
  32
  33 /*****************************************
  34  * Returns x rounded to a long value using the FE_TONEAREST rounding mode.
  35  * If the integer value of x is
  36  * greater than long.max, the result is
  37  * indeterminate.
  38  */
  39 extern (C) real rndtonl(real x);
  40
  41 pure:
  42 /***********************************
  43  * Returns cosine of x. x is in radians.
  44  *
  45  *      $(TABLE_SV
  46  *      $(TR $(TH x)                 $(TH cos(x)) $(TH invalid?))
  47  *      $(TR $(TD $(NAN))            $(TD $(NAN)) $(TD yes)     )
  48  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(NAN)) $(TD yes)     )
  49  *      )
  50  * Bugs:
  51  *      Results are undefined if |x| >= $(POWER 2,64).
  52  */
  53
  54 float cos(float x);     /* intrinsic */
  55 double cos(double x);   /* intrinsic */ /// ditto
  56 real cos(real x);       /* intrinsic */ /// ditto
  57
  58 /***********************************
  59  * Returns sine of x. x is in radians.
  60  *
  61  *      $(TABLE_SV
  62  *      $(TR $(TH x)               $(TH sin(x))      $(TH invalid?))
  63  *      $(TR $(TD $(NAN))          $(TD $(NAN))      $(TD yes))
  64  *      $(TR $(TD $(PLUSMN)0.0)    $(TD $(PLUSMN)0.0) $(TD no))
  65  *      $(TR $(TD $(PLUSMNINF))    $(TD $(NAN))      $(TD yes))
  66  *      )
  67  * Bugs:
  68  *      Results are undefined if |x| >= $(POWER 2,64).
  69  */
  70
  71 float sin(float x);     /* intrinsic */
  72 double sin(double x);   /* intrinsic */ /// ditto
  73 real sin(real x);       /* intrinsic */ /// ditto
  74
  75 /*****************************************
  76  * Returns x rounded to a long value using the current rounding mode.
  77  * If the integer value of x is
  78  * greater than long.max, the result is
  79  * indeterminate.
  80  */
  81
  82 long rndtol(float x);   /* intrinsic */
  83 long rndtol(double x);  /* intrinsic */ /// ditto
  84 long rndtol(real x);    /* intrinsic */ /// ditto
  85
  86 /***************************************
  87  * Compute square root of x.
  88  *
  89  *      $(TABLE_SV
  90  *      $(TR $(TH x)         $(TH sqrt(x))   $(TH invalid?))
  91  *      $(TR $(TD -0.0)      $(TD -0.0)      $(TD no))
  92  *      $(TR $(TD $(LT)0.0)  $(TD $(NAN))    $(TD yes))
  93  *      $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no))
  94  *      )
  95  */
  96
  97 float sqrt(float x);    /* intrinsic */
  98 double sqrt(double x);  /* intrinsic */ /// ditto
  99 real sqrt(real x);      /* intrinsic */ /// ditto
 100
 101 /*******************************************
 102  * Compute n * 2$(SUPERSCRIPT exp)
 103  * References: frexp
 104  */
 105
 106 float ldexp(float n, int exp);   /* intrinsic */
 107 double ldexp(double n, int exp); /* intrinsic */ /// ditto
 108 real ldexp(real n, int exp);     /* intrinsic */ /// ditto
 109
 110 unittest {
 111     static if (real.mant_dig == 113)
 112     {
 113         assert(ldexp(1.0L, -16384) == 0x1p-16384L);
 114         assert(ldexp(1.0L, -16382) == 0x1p-16382L);
 115     }
 116     else static if (real.mant_dig == 106)
 117     {
 118         assert(ldexp(1.0L,  1023) == 0x1p1023L);
 119         assert(ldexp(1.0L, -1022) == 0x1p-1022L);
 120         assert(ldexp(1.0L, -1021) == 0x1p-1021L);
 121     }
 122     else static if (real.mant_dig == 64)
 123     {
 124         assert(ldexp(1.0L, -16384) == 0x1p-16384L);
 125         assert(ldexp(1.0L, -16382) == 0x1p-16382L);
 126     }
 127     else static if (real.mant_dig == 53)
 128     {
 129         assert(ldexp(1.0L,  1023) == 0x1p1023L);
 130         assert(ldexp(1.0L, -1022) == 0x1p-1022L);
 131         assert(ldexp(1.0L, -1021) == 0x1p-1021L);
 132     }
 133     else
 134         assert(false, "Only 128bit, 80bit and 64bit reals expected here");
 135 }
 136
 137 /*******************************
 138  * Compute the absolute value.
 139  *      $(TABLE_SV
 140  *      $(TR $(TH x)                 $(TH fabs(x)))
 141  *      $(TR $(TD $(PLUSMN)0.0)      $(TD +0.0) )
 142  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD +$(INFIN)) )
 143  *      )
 144  * It is implemented as a compiler intrinsic.
 145  * Params:
 146  *      x = floating point value
 147  * Returns: |x|
 148  * References: equivalent to `std.math.fabs`
 149  */
 150 @safe pure nothrow @nogc
 151 {
 152     float  fabs(float  x);
 153     double fabs(double x); /// ditto
 154     real   fabs(real   x); /// ditto
 155 }
 156
 157 /**********************************
 158  * Rounds x to the nearest integer value, using the current rounding
 159  * mode.
 160  * If the return value is not equal to x, the FE_INEXACT
 161  * exception is raised.
 162  * $(B nearbyint) performs
 163  * the same operation, but does not set the FE_INEXACT exception.
 164  */
 165 float rint(float x);    /* intrinsic */
 166 double rint(double x);  /* intrinsic */ /// ditto
 167 real rint(real x);      /* intrinsic */ /// ditto
 168
 169 /***********************************
 170  * Building block functions, they
 171  * translate to a single x87 instruction.
 172  */
 173 // y * log2(x)
 174 float yl2x(float x, float y);    /* intrinsic */
 175 double yl2x(double x, double y);  /* intrinsic */ /// ditto
 176 real yl2x(real x, real y);      /* intrinsic */ /// ditto
 177 // y * log2(x +1)
 178 float yl2xp1(float x, float y);    /* intrinsic */
 179 double yl2xp1(double x, double y);  /* intrinsic */ /// ditto
 180 real yl2xp1(real x, real y);      /* intrinsic */ /// ditto
 181
 182 unittest
 183 {
 184     version (INLINE_YL2X)
 185     {
 186         assert(yl2x(1024.0L, 1) == 10);
 187         assert(yl2xp1(1023.0L, 1) == 10);
 188     }
 189 }
 190
 191 /*************************************
 192  * Round argument to a specific precision.
 193  *
 194  * D language types specify only a minimum precision, not a maximum. The
 195  * `toPrec()` function forces rounding of the argument `f` to the precision
 196  * of the specified floating point type `T`.
 197  * The rounding mode used is inevitably target-dependent, but will be done in
 198  * a way to maximize accuracy. In most cases, the default is round-to-nearest.
 199  *
 200  * Params:
 201  *      T = precision type to round to
 202  *      f = value to convert
 203  * Returns:
 204  *      f in precision of type `T`
 205  */
 206 T toPrec(T:float)(float f) { pragma(inline, false); return f; }
 207 /// ditto
 208 T toPrec(T:float)(double f) { pragma(inline, false); return cast(T) f; }
 209 /// ditto
 210 T toPrec(T:float)(real f)  { pragma(inline, false); return cast(T) f; }
 211 /// ditto
 212 T toPrec(T:double)(float f) { pragma(inline, false); return f; }
 213 /// ditto
 214 T toPrec(T:double)(double f) { pragma(inline, false); return f; }
 215 /// ditto
 216 T toPrec(T:double)(real f)  { pragma(inline, false); return cast(T) f; }
 217 /// ditto
 218 T toPrec(T:real)(float f) { pragma(inline, false); return f; }
 219 /// ditto
 220 T toPrec(T:real)(double f) { pragma(inline, false); return f; }
 221 /// ditto
 222 T toPrec(T:real)(real f)  { pragma(inline, false); return f; }
 223
 224 @safe unittest
 225 {
 226     // Test all instantiations work with all combinations of float.
 227     float f = 1.1f;
 228     double d = 1.1;
 229     real r = 1.1L;
 230     f = toPrec!float(f + f);
 231     f = toPrec!float(d + d);
 232     f = toPrec!float(r + r);
 233     d = toPrec!double(f + f);
 234     d = toPrec!double(d + d);
 235     d = toPrec!double(r + r);
 236     r = toPrec!real(f + f);
 237     r = toPrec!real(d + d);
 238     r = toPrec!real(r + r);
 239
 240     // Comparison tests.
 241     bool approxEqual(T)(T lhs, T rhs)
 242     {
 243         return fabs((lhs - rhs) / rhs) <= 1e-2 || fabs(lhs - rhs) <= 1e-5;
 244     }
 245
 246     enum real PIR = 0xc.90fdaa22168c235p-2;
 247     enum double PID = 0x1.921fb54442d18p+1;
 248     enum float PIF = 0x1.921fb6p+1;
 249     static assert(approxEqual(toPrec!float(PIR), PIF));
 250     static assert(approxEqual(toPrec!double(PIR), PID));
 251     static assert(approxEqual(toPrec!real(PIR), PIR));
 252     static assert(approxEqual(toPrec!float(PID), PIF));
 253     static assert(approxEqual(toPrec!double(PID), PID));
 254     static assert(approxEqual(toPrec!real(PID), PID));
 255     static assert(approxEqual(toPrec!float(PIF), PIF));
 256     static assert(approxEqual(toPrec!double(PIF), PIF));
 257     static assert(approxEqual(toPrec!real(PIF), PIF));
 258
 259     assert(approxEqual(toPrec!float(PIR), PIF));
 260     assert(approxEqual(toPrec!double(PIR), PID));
 261     assert(approxEqual(toPrec!real(PIR), PIR));
 262     assert(approxEqual(toPrec!float(PID), PIF));
 263     assert(approxEqual(toPrec!double(PID), PID));
 264     assert(approxEqual(toPrec!real(PID), PID));
 265     assert(approxEqual(toPrec!float(PIF), PIF));
 266     assert(approxEqual(toPrec!double(PIF), PIF));
 267     assert(approxEqual(toPrec!real(PIF), PIF));
 268 }