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 $(MREF 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 package(std): // Not public yet
 172 /* Return the value that lies halfway between x and y on the IEEE number line.
 173  *
 174  * Formally, the result is the arithmetic mean of the binary significands of x
 175  * and y, multiplied by the geometric mean of the binary exponents of x and y.
 176  * x and y must have the same sign, and must not be NaN.
 177  * Note: this function is useful for ensuring O(log n) behaviour in algorithms
 178  * involving a 'binary chop'.
 179  *
 180  * Special cases:
 181  * If x and y are within a factor of 2, (ie, feqrel(x, y) > 0), the return value
 182  * is the arithmetic mean (x + y) / 2.
 183  * If x and y are even powers of 2, the return value is the geometric mean,
 184  *   ieeeMean(x, y) = sqrt(x * y).
 185  *
 186  */
 187 T ieeeMean(T)(const T x, const T y)  @trusted pure nothrow @nogc
 188 in
 189 {
 190     // both x and y must have the same sign, and must not be NaN.
 191     assert(signbit(x) == signbit(y));
 192     assert(x == x && y == y);
 193 }
 194 do
 195 {
 196     // Runtime behaviour for contract violation:
 197     // If signs are opposite, or one is a NaN, return 0.
 198     if (!((x >= 0 && y >= 0) || (x <= 0 && y <= 0))) return 0.0;
 199
 200     // The implementation is simple: cast x and y to integers,
 201     // average them (avoiding overflow), and cast the result back to a floating-point number.
 202
 203     alias F = floatTraits!(T);
 204     T u;
 205     static if (F.realFormat == RealFormat.ieeeExtended ||
 206                F.realFormat == RealFormat.ieeeExtended53)
 207     {
 208         // There's slight additional complexity because they are actually
 209         // 79-bit reals...
 210         ushort *ue = cast(ushort *)&u;
 211         ulong *ul = cast(ulong *)&u;
 212         ushort *xe = cast(ushort *)&x;
 213         ulong *xl = cast(ulong *)&x;
 214         ushort *ye = cast(ushort *)&y;
 215         ulong *yl = cast(ulong *)&y;
 216
 217         // Ignore the useless implicit bit. (Bonus: this prevents overflows)
 218         ulong m = ((*xl) & 0x7FFF_FFFF_FFFF_FFFFL) + ((*yl) & 0x7FFF_FFFF_FFFF_FFFFL);
 219
 220         // @@@ BUG? @@@
 221         // Cast shouldn't be here
 222         ushort e = cast(ushort) ((xe[F.EXPPOS_SHORT] & F.EXPMASK)
 223                                  + (ye[F.EXPPOS_SHORT] & F.EXPMASK));
 224         if (m & 0x8000_0000_0000_0000L)
 225         {
 226             ++e;
 227             m &= 0x7FFF_FFFF_FFFF_FFFFL;
 228         }
 229         // Now do a multi-byte right shift
 230         const uint c = e & 1; // carry
 231         e >>= 1;
 232         m >>>= 1;
 233         if (c)
 234             m |= 0x4000_0000_0000_0000L; // shift carry into significand
 235         if (e)
 236             *ul = m | 0x8000_0000_0000_0000L; // set implicit bit...
 237         else
 238             *ul = m; // ... unless exponent is 0 (subnormal or zero).
 239
 240         ue[4]= e | (xe[F.EXPPOS_SHORT]& 0x8000); // restore sign bit
 241     }
 242     else static if (F.realFormat == RealFormat.ieeeQuadruple)
 243     {
 244         // This would be trivial if 'ucent' were implemented...
 245         ulong *ul = cast(ulong *)&u;
 246         ulong *xl = cast(ulong *)&x;
 247         ulong *yl = cast(ulong *)&y;
 248
 249         // Multi-byte add, then multi-byte right shift.
 250         import core.checkedint : addu;
 251         bool carry;
 252         ulong ml = addu(xl[MANTISSA_LSB], yl[MANTISSA_LSB], carry);
 253
 254         ulong mh = carry + (xl[MANTISSA_MSB] & 0x7FFF_FFFF_FFFF_FFFFL) +
 255             (yl[MANTISSA_MSB] & 0x7FFF_FFFF_FFFF_FFFFL);
 256
 257         ul[MANTISSA_MSB] = (mh >>> 1) | (xl[MANTISSA_MSB] & 0x8000_0000_0000_0000);
 258         ul[MANTISSA_LSB] = (ml >>> 1) | (mh & 1) << 63;
 259     }
 260     else static if (F.realFormat == RealFormat.ieeeDouble)
 261     {
 262         ulong *ul = cast(ulong *)&u;
 263         ulong *xl = cast(ulong *)&x;
 264         ulong *yl = cast(ulong *)&y;
 265         ulong m = (((*xl) & 0x7FFF_FFFF_FFFF_FFFFL)
 266                    + ((*yl) & 0x7FFF_FFFF_FFFF_FFFFL)) >>> 1;
 267         m |= ((*xl) & 0x8000_0000_0000_0000L);
 268         *ul = m;
 269     }
 270     else static if (F.realFormat == RealFormat.ieeeSingle)
 271     {
 272         uint *ul = cast(uint *)&u;
 273         uint *xl = cast(uint *)&x;
 274         uint *yl = cast(uint *)&y;
 275         uint m = (((*xl) & 0x7FFF_FFFF) + ((*yl) & 0x7FFF_FFFF)) >>> 1;
 276         m |= ((*xl) & 0x8000_0000);
 277         *ul = m;
 278     }
 279     else
 280     {
 281         assert(0, "Not implemented");
 282     }
 283     return u;
 284 }
 285
 286 @safe pure nothrow @nogc unittest
 287 {
 288     assert(ieeeMean(-0.0,-1e-20)<0);
 289     assert(ieeeMean(0.0,1e-20)>0);
 290
 291     assert(ieeeMean(1.0L,4.0L)==2L);
 292     assert(ieeeMean(2.0*1.013,8.0*1.013)==4*1.013);
 293     assert(ieeeMean(-1.0L,-4.0L)==-2L);
 294     assert(ieeeMean(-1.0,-4.0)==-2);
 295     assert(ieeeMean(-1.0f,-4.0f)==-2f);
 296     assert(ieeeMean(-1.0,-2.0)==-1.5);
 297     assert(ieeeMean(-1*(1+8*real.epsilon),-2*(1+8*real.epsilon))
 298                  ==-1.5*(1+5*real.epsilon));
 299     assert(ieeeMean(0x1p60,0x1p-10)==0x1p25);
 300
 301     static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended)
 302     {
 303       assert(ieeeMean(1.0L,real.infinity)==0x1p8192L);
 304       assert(ieeeMean(0.0L,real.infinity)==1.5);
 305     }
 306     assert(ieeeMean(0.5*real.min_normal*(1-4*real.epsilon),0.5*real.min_normal)
 307            == 0.5*real.min_normal*(1-2*real.epsilon));
 308 }
 309
 310
 311 // The following IEEE 'real' formats are currently supported.
 312 version (LittleEndian)
 313 {
 314     static assert(real.mant_dig == 53 || real.mant_dig == 64
 315                || real.mant_dig == 113,
 316       "Only 64-bit, 80-bit, and 128-bit reals"~
 317       " are supported for LittleEndian CPUs");
 318 }
 319 else
 320 {
 321     static assert(real.mant_dig == 53 || real.mant_dig == 113,
 322     "Only 64-bit and 128-bit reals are supported for BigEndian CPUs.");
 323 }