libphobos/src/std/math/package.d

   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 $(SUBMODULE Constants, constants)) $(TD
  14     $(SUBREF constants, E)
  15     $(SUBREF constants, PI)
  16     $(SUBREF constants, PI_2)
  17     $(SUBREF constants, PI_4)
  18     $(SUBREF constants, M_1_PI)
  19     $(SUBREF constants, M_2_PI)
  20     $(SUBREF constants, M_2_SQRTPI)
  21     $(SUBREF constants, LN10)
  22     $(SUBREF constants, LN2)
  23     $(SUBREF constants, LOG2)
  24     $(SUBREF constants, LOG2E)
  25     $(SUBREF constants, LOG2T)
  26     $(SUBREF constants, LOG10E)
  27     $(SUBREF constants, SQRT2)
  28     $(SUBREF constants, SQRT1_2)
  29 ))
  30 $(TR $(TDNW $(SUBMODULE Algebraic, algebraic)) $(TD
  31     $(SUBREF algebraic, abs)
  32     $(SUBREF algebraic, fabs)
  33     $(SUBREF algebraic, sqrt)
  34     $(SUBREF algebraic, cbrt)
  35     $(SUBREF algebraic, hypot)
  36     $(SUBREF algebraic, poly)
  37     $(SUBREF algebraic, nextPow2)
  38     $(SUBREF algebraic, truncPow2)
  39 ))
  40 $(TR $(TDNW $(SUBMODULE Trigonometry, trigonometry)) $(TD
  41     $(SUBREF trigonometry, sin)
  42     $(SUBREF trigonometry, cos)
  43     $(SUBREF trigonometry, tan)
  44     $(SUBREF trigonometry, asin)
  45     $(SUBREF trigonometry, acos)
  46     $(SUBREF trigonometry, atan)
  47     $(SUBREF trigonometry, atan2)
  48     $(SUBREF trigonometry, sinh)
  49     $(SUBREF trigonometry, cosh)
  50     $(SUBREF trigonometry, tanh)
  51     $(SUBREF trigonometry, asinh)
  52     $(SUBREF trigonometry, acosh)
  53     $(SUBREF trigonometry, atanh)
  54 ))
  55 $(TR $(TDNW $(SUBMODULE Rounding, rounding)) $(TD
  56     $(SUBREF rounding, ceil)
  57     $(SUBREF rounding, floor)
  58     $(SUBREF rounding, round)
  59     $(SUBREF rounding, lround)
  60     $(SUBREF rounding, trunc)
  61     $(SUBREF rounding, rint)
  62     $(SUBREF rounding, lrint)
  63     $(SUBREF rounding, nearbyint)
  64     $(SUBREF rounding, rndtol)
  65     $(SUBREF rounding, quantize)
  66 ))
  67 $(TR $(TDNW $(SUBMODULE Exponentiation & Logarithms, exponential)) $(TD
  68     $(SUBREF exponential, pow)
  69     $(SUBREF exponential, powmod)
  70     $(SUBREF exponential, exp)
  71     $(SUBREF exponential, exp2)
  72     $(SUBREF exponential, expm1)
  73     $(SUBREF exponential, ldexp)
  74     $(SUBREF exponential, frexp)
  75     $(SUBREF exponential, log)
  76     $(SUBREF exponential, log2)
  77     $(SUBREF exponential, log10)
  78     $(SUBREF exponential, logb)
  79     $(SUBREF exponential, ilogb)
  80     $(SUBREF exponential, log1p)
  81     $(SUBREF exponential, scalbn)
  82 ))
  83 $(TR $(TDNW $(SUBMODULE Remainder, remainder)) $(TD
  84     $(SUBREF remainder, fmod)
  85     $(SUBREF remainder, modf)
  86     $(SUBREF remainder, remainder)
  87     $(SUBREF remainder, remquo)
  88 ))
  89 $(TR $(TDNW $(SUBMODULE Floating-point operations, operations)) $(TD
  90     $(SUBREF operations, approxEqual)
  91     $(SUBREF operations, feqrel)
  92     $(SUBREF operations, fdim)
  93     $(SUBREF operations, fmax)
  94     $(SUBREF operations, fmin)
  95     $(SUBREF operations, fma)
  96     $(SUBREF operations, isClose)
  97     $(SUBREF operations, nextDown)
  98     $(SUBREF operations, nextUp)
  99     $(SUBREF operations, nextafter)
 100     $(SUBREF operations, NaN)
 101     $(SUBREF operations, getNaNPayload)
 102     $(SUBREF operations, cmp)
 103 ))
 104 $(TR $(TDNW $(SUBMODULE Introspection, traits)) $(TD
 105     $(SUBREF traits, isFinite)
 106     $(SUBREF traits, isIdentical)
 107     $(SUBREF traits, isInfinity)
 108     $(SUBREF traits, isNaN)
 109     $(SUBREF traits, isNormal)
 110     $(SUBREF traits, isSubnormal)
 111     $(SUBREF traits, signbit)
 112     $(SUBREF traits, sgn)
 113     $(SUBREF traits, copysign)
 114     $(SUBREF traits, isPowerOf2)
 115 ))
 116 $(TR $(TDNW $(SUBMODULE Hardware Control, hardware)) $(TD
 117     $(SUBREF hardware, IeeeFlags)
 118     $(SUBREF hardware, ieeeFlags)
 119     $(SUBREF hardware, resetIeeeFlags)
 120     $(SUBREF hardware, FloatingPointControl)
 121 ))
 122 )
 123 )
 124
 125  * The functionality closely follows the IEEE754-2008 standard for
 126  * floating-point arithmetic, including the use of camelCase names rather
 127  * than C99-style lower case names. All of these functions behave correctly
 128  * when presented with an infinity or NaN.
 129  *
 130  * The following IEEE 'real' formats are currently supported:
 131  * $(UL
 132  * $(LI 64 bit Big-endian  'double' (eg PowerPC))
 133  * $(LI 128 bit Big-endian 'quadruple' (eg SPARC))
 134  * $(LI 64 bit Little-endian 'double' (eg x86-SSE2))
 135  * $(LI 80 bit Little-endian, with implied bit 'real80' (eg x87, Itanium))
 136  * $(LI 128 bit Little-endian 'quadruple' (not implemented on any known processor!))
 137  * $(LI Non-IEEE 128 bit Big-endian 'doubledouble' (eg PowerPC) has partial support)
 138  * )
 139  * Unlike C, there is no global 'errno' variable. Consequently, almost all of
 140  * these functions are pure nothrow.
 141  *
 142  * Macros:
 143  *      SUBMODULE = $(MREF_ALTTEXT $1, std, math, $2)
 144  *      SUBREF = $(REF_ALTTEXT $(TT $2), $2, std, math, $1)$(NBSP)
 145  *
 146  * Copyright: Copyright The D Language Foundation 2000 - 2011.
 147  *            D implementations of tan, atan, atan2, exp, expm1, exp2, log, log10, log1p,
 148  *            log2, floor, ceil and lrint functions are based on the CEPHES math library,
 149  *            which is Copyright (C) 2001 Stephen L. Moshier $(LT)steve@moshier.net$(GT)
 150  *            and are incorporated herein by permission of the author.  The author
 151  *            reserves the right to distribute this material elsewhere under different
 152  *            copying permissions.  These modifications are distributed here under
 153  *            the following terms:
 154  * License:   $(HTTP www.boost.org/LICENSE_1_0.txt, Boost License 1.0).
 155  * Authors:   $(HTTP digitalmars.com, Walter Bright), Don Clugston,
 156  *            Conversion of CEPHES math library to D by Iain Buclaw and David Nadlinger
 157  * Source: $(PHOBOSSRC std/math/package.d)
 158  */
 159 module std.math;
 160
 161 public import std.math.algebraic;
 162 public import std.math.constants;
 163 public import std.math.exponential;
 164 public import std.math.operations;
 165 public import std.math.hardware;
 166 public import std.math.remainder;
 167 public import std.math.rounding;
 168 public import std.math.traits;
 169 public import std.math.trigonometry;
 170
 171 // @@@DEPRECATED_2.102@@@
 172 // Note: Exposed accidentally, should be deprecated / removed
 173 deprecated("std.meta.AliasSeq was unintentionally available from std.math "
 174            ~ "and will be removed after 2.102. Please import std.meta instead")
 175 public import std.meta : AliasSeq;
 176
 177 package(std): // Not public yet
 178 /* Return the value that lies halfway between x and y on the IEEE number line.
 179  *
 180  * Formally, the result is the arithmetic mean of the binary significands of x
 181  * and y, multiplied by the geometric mean of the binary exponents of x and y.
 182  * x and y must have the same sign, and must not be NaN.
 183  * Note: this function is useful for ensuring O(log n) behaviour in algorithms
 184  * involving a 'binary chop'.
 185  *
 186  * Special cases:
 187  * If x and y are within a factor of 2, (ie, feqrel(x, y) > 0), the return value
 188  * is the arithmetic mean (x + y) / 2.
 189  * If x and y are even powers of 2, the return value is the geometric mean,
 190  *   ieeeMean(x, y) = sqrt(x * y).
 191  *
 192  */
 193 T ieeeMean(T)(const T x, const T y)  @trusted pure nothrow @nogc
 194 in
 195 {
 196     // both x and y must have the same sign, and must not be NaN.
 197     assert(signbit(x) == signbit(y));
 198     assert(x == x && y == y);
 199 }
 200 do
 201 {
 202     // Runtime behaviour for contract violation:
 203     // If signs are opposite, or one is a NaN, return 0.
 204     if (!((x >= 0 && y >= 0) || (x <= 0 && y <= 0))) return 0.0;
 205
 206     // The implementation is simple: cast x and y to integers,
 207     // average them (avoiding overflow), and cast the result back to a floating-point number.
 208
 209     alias F = floatTraits!(T);
 210     T u;
 211     static if (F.realFormat == RealFormat.ieeeExtended ||
 212                F.realFormat == RealFormat.ieeeExtended53)
 213     {
 214         // There's slight additional complexity because they are actually
 215         // 79-bit reals...
 216         ushort *ue = cast(ushort *)&u;
 217         ulong *ul = cast(ulong *)&u;
 218         ushort *xe = cast(ushort *)&x;
 219         ulong *xl = cast(ulong *)&x;
 220         ushort *ye = cast(ushort *)&y;
 221         ulong *yl = cast(ulong *)&y;
 222
 223         // Ignore the useless implicit bit. (Bonus: this prevents overflows)
 224         ulong m = ((*xl) & 0x7FFF_FFFF_FFFF_FFFFL) + ((*yl) & 0x7FFF_FFFF_FFFF_FFFFL);
 225
 226         // @@@ BUG? @@@
 227         // Cast shouldn't be here
 228         ushort e = cast(ushort) ((xe[F.EXPPOS_SHORT] & F.EXPMASK)
 229                                  + (ye[F.EXPPOS_SHORT] & F.EXPMASK));
 230         if (m & 0x8000_0000_0000_0000L)
 231         {
 232             ++e;
 233             m &= 0x7FFF_FFFF_FFFF_FFFFL;
 234         }
 235         // Now do a multi-byte right shift
 236         const uint c = e & 1; // carry
 237         e >>= 1;
 238         m >>>= 1;
 239         if (c)
 240             m |= 0x4000_0000_0000_0000L; // shift carry into significand
 241         if (e)
 242             *ul = m | 0x8000_0000_0000_0000L; // set implicit bit...
 243         else
 244             *ul = m; // ... unless exponent is 0 (subnormal or zero).
 245
 246         ue[4]= e | (xe[F.EXPPOS_SHORT]& 0x8000); // restore sign bit
 247     }
 248     else static if (F.realFormat == RealFormat.ieeeQuadruple)
 249     {
 250         // This would be trivial if 'ucent' were implemented...
 251         ulong *ul = cast(ulong *)&u;
 252         ulong *xl = cast(ulong *)&x;
 253         ulong *yl = cast(ulong *)&y;
 254
 255         // Multi-byte add, then multi-byte right shift.
 256         import core.checkedint : addu;
 257         bool carry;
 258         ulong ml = addu(xl[MANTISSA_LSB], yl[MANTISSA_LSB], carry);
 259
 260         ulong mh = carry + (xl[MANTISSA_MSB] & 0x7FFF_FFFF_FFFF_FFFFL) +
 261             (yl[MANTISSA_MSB] & 0x7FFF_FFFF_FFFF_FFFFL);
 262
 263         ul[MANTISSA_MSB] = (mh >>> 1) | (xl[MANTISSA_MSB] & 0x8000_0000_0000_0000);
 264         ul[MANTISSA_LSB] = (ml >>> 1) | (mh & 1) << 63;
 265     }
 266     else static if (F.realFormat == RealFormat.ieeeDouble)
 267     {
 268         ulong *ul = cast(ulong *)&u;
 269         ulong *xl = cast(ulong *)&x;
 270         ulong *yl = cast(ulong *)&y;
 271         ulong m = (((*xl) & 0x7FFF_FFFF_FFFF_FFFFL)
 272                    + ((*yl) & 0x7FFF_FFFF_FFFF_FFFFL)) >>> 1;
 273         m |= ((*xl) & 0x8000_0000_0000_0000L);
 274         *ul = m;
 275     }
 276     else static if (F.realFormat == RealFormat.ieeeSingle)
 277     {
 278         uint *ul = cast(uint *)&u;
 279         uint *xl = cast(uint *)&x;
 280         uint *yl = cast(uint *)&y;
 281         uint m = (((*xl) & 0x7FFF_FFFF) + ((*yl) & 0x7FFF_FFFF)) >>> 1;
 282         m |= ((*xl) & 0x8000_0000);
 283         *ul = m;
 284     }
 285     else
 286     {
 287         assert(0, "Not implemented");
 288     }
 289     return u;
 290 }
 291
 292 @safe pure nothrow @nogc unittest
 293 {
 294     assert(ieeeMean(-0.0,-1e-20)<0);
 295     assert(ieeeMean(0.0,1e-20)>0);
 296
 297     assert(ieeeMean(1.0L,4.0L)==2L);
 298     assert(ieeeMean(2.0*1.013,8.0*1.013)==4*1.013);
 299     assert(ieeeMean(-1.0L,-4.0L)==-2L);
 300     assert(ieeeMean(-1.0,-4.0)==-2);
 301     assert(ieeeMean(-1.0f,-4.0f)==-2f);
 302     assert(ieeeMean(-1.0,-2.0)==-1.5);
 303     assert(ieeeMean(-1*(1+8*real.epsilon),-2*(1+8*real.epsilon))
 304                  ==-1.5*(1+5*real.epsilon));
 305     assert(ieeeMean(0x1p60,0x1p-10)==0x1p25);
 306
 307     static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended)
 308     {
 309       assert(ieeeMean(1.0L,real.infinity)==0x1p8192L);
 310       assert(ieeeMean(0.0L,real.infinity)==1.5);
 311     }
 312     assert(ieeeMean(0.5*real.min_normal*(1-4*real.epsilon),0.5*real.min_normal)
 313            == 0.5*real.min_normal*(1-2*real.epsilon));
 314 }
 315
 316
 317 // The following IEEE 'real' formats are currently supported.
 318 version (LittleEndian)
 319 {
 320     static assert(real.mant_dig == 53 || real.mant_dig == 64
 321                || real.mant_dig == 113,
 322       "Only 64-bit, 80-bit, and 128-bit reals"~
 323       " are supported for LittleEndian CPUs");
 324 }
 325 else
 326 {
 327     static assert(real.mant_dig == 53 || real.mant_dig == 113,
 328     "Only 64-bit and 128-bit reals are supported for BigEndian CPUs.");
 329 }
 330
 331 // Underlying format exposed through floatTraits
 332 enum RealFormat
 333 {
 334     ieeeHalf,
 335     ieeeSingle,
 336     ieeeDouble,
 337     ieeeExtended,   // x87 80-bit real
 338     ieeeExtended53, // x87 real rounded to precision of double.
 339     ibmExtended,    // IBM 128-bit extended
 340     ieeeQuadruple,
 341 }
 342
 343 // Constants used for extracting the components of the representation.
 344 // They supplement the built-in floating point properties.
 345 template floatTraits(T)
 346 {
 347     import std.traits : Unqual;
 348
 349     // EXPMASK is a ushort mask to select the exponent portion (without sign)
 350     // EXPSHIFT is the number of bits the exponent is left-shifted by in its ushort
 351     // EXPBIAS is the exponent bias - 1 (exp == EXPBIAS yields ×2^-1).
 352     // EXPPOS_SHORT is the index of the exponent when represented as a ushort array.
 353     // SIGNPOS_BYTE is the index of the sign when represented as a ubyte array.
 354     // RECIP_EPSILON is the value such that (smallest_subnormal) * RECIP_EPSILON == T.min_normal
 355     enum Unqual!T RECIP_EPSILON = (1/T.epsilon);
 356     static if (T.mant_dig == 24)
 357     {
 358         // Single precision float
 359         enum ushort EXPMASK = 0x7F80;
 360         enum ushort EXPSHIFT = 7;
 361         enum ushort EXPBIAS = 0x3F00;
 362         enum uint EXPMASK_INT = 0x7F80_0000;
 363         enum uint MANTISSAMASK_INT = 0x007F_FFFF;
 364         enum realFormat = RealFormat.ieeeSingle;
 365         version (LittleEndian)
 366         {
 367             enum EXPPOS_SHORT = 1;
 368             enum SIGNPOS_BYTE = 3;
 369         }
 370         else
 371         {
 372             enum EXPPOS_SHORT = 0;
 373             enum SIGNPOS_BYTE = 0;
 374         }
 375     }
 376     else static if (T.mant_dig == 53)
 377     {
 378         static if (T.sizeof == 8)
 379         {
 380             // Double precision float, or real == double
 381             enum ushort EXPMASK = 0x7FF0;
 382             enum ushort EXPSHIFT = 4;
 383             enum ushort EXPBIAS = 0x3FE0;
 384             enum uint EXPMASK_INT = 0x7FF0_0000;
 385             enum uint MANTISSAMASK_INT = 0x000F_FFFF; // for the MSB only
 386             enum realFormat = RealFormat.ieeeDouble;
 387             version (LittleEndian)
 388             {
 389                 enum EXPPOS_SHORT = 3;
 390                 enum SIGNPOS_BYTE = 7;
 391             }
 392             else
 393             {
 394                 enum EXPPOS_SHORT = 0;
 395                 enum SIGNPOS_BYTE = 0;
 396             }
 397         }
 398         else static if (T.sizeof == 12)
 399         {
 400             // Intel extended real80 rounded to double
 401             enum ushort EXPMASK = 0x7FFF;
 402             enum ushort EXPSHIFT = 0;
 403             enum ushort EXPBIAS = 0x3FFE;
 404             enum realFormat = RealFormat.ieeeExtended53;
 405             version (LittleEndian)
 406             {
 407                 enum EXPPOS_SHORT = 4;
 408                 enum SIGNPOS_BYTE = 9;
 409             }
 410             else
 411             {
 412                 enum EXPPOS_SHORT = 0;
 413                 enum SIGNPOS_BYTE = 0;
 414             }
 415         }
 416         else
 417             static assert(false, "No traits support for " ~ T.stringof);
 418     }
 419     else static if (T.mant_dig == 64)
 420     {
 421         // Intel extended real80
 422         enum ushort EXPMASK = 0x7FFF;
 423         enum ushort EXPSHIFT = 0;
 424         enum ushort EXPBIAS = 0x3FFE;
 425         enum realFormat = RealFormat.ieeeExtended;
 426         version (LittleEndian)
 427         {
 428             enum EXPPOS_SHORT = 4;
 429             enum SIGNPOS_BYTE = 9;
 430         }
 431         else
 432         {
 433             enum EXPPOS_SHORT = 0;
 434             enum SIGNPOS_BYTE = 0;
 435         }
 436     }
 437     else static if (T.mant_dig == 113)
 438     {
 439         // Quadruple precision float
 440         enum ushort EXPMASK = 0x7FFF;
 441         enum ushort EXPSHIFT = 0;
 442         enum ushort EXPBIAS = 0x3FFE;
 443         enum realFormat = RealFormat.ieeeQuadruple;
 444         version (LittleEndian)
 445         {
 446             enum EXPPOS_SHORT = 7;
 447             enum SIGNPOS_BYTE = 15;
 448         }
 449         else
 450         {
 451             enum EXPPOS_SHORT = 0;
 452             enum SIGNPOS_BYTE = 0;
 453         }
 454     }
 455     else static if (T.mant_dig == 106)
 456     {
 457         // IBM Extended doubledouble
 458         enum ushort EXPMASK = 0x7FF0;
 459         enum ushort EXPSHIFT = 4;
 460         enum realFormat = RealFormat.ibmExtended;
 461
 462         // For IBM doubledouble the larger magnitude double comes first.
 463         // It's really a double[2] and arrays don't index differently
 464         // between little and big-endian targets.
 465         enum DOUBLEPAIR_MSB = 0;
 466         enum DOUBLEPAIR_LSB = 1;
 467
 468         // The exponent/sign byte is for most significant part.
 469         version (LittleEndian)
 470         {
 471             enum EXPPOS_SHORT = 3;
 472             enum SIGNPOS_BYTE = 7;
 473         }
 474         else
 475         {
 476             enum EXPPOS_SHORT = 0;
 477             enum SIGNPOS_BYTE = 0;
 478         }
 479     }
 480     else
 481         static assert(false, "No traits support for " ~ T.stringof);
 482 }
 483
 484 // These apply to all floating-point types
 485 version (LittleEndian)
 486 {
 487     enum MANTISSA_LSB = 0;
 488     enum MANTISSA_MSB = 1;
 489 }
 490 else
 491 {
 492     enum MANTISSA_LSB = 1;
 493     enum MANTISSA_MSB = 0;
 494 }