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 }