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 = ± 16 * INFIN = ∞ 17 * PLUSMNINF = ±∞ 18 * LT = < 19 * 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 version (LDC) 29 { 30 import ldc.intrinsics; 31 32 private enum isRealX87 = (real.mant_dig == 64); 33 } 34 35 public: 36 @nogc: 37 nothrow: 38 @safe: 39 40 /***************************************** 41 * Returns x rounded to a long value using the FE_TONEAREST rounding mode. 42 * If the integer value of x is 43 * greater than long.max, the result is 44 * indeterminate. 45 */ 46 deprecated("rndtonl is to be removed by 2.100. Please use round instead") 47 extern (C) real rndtonl(real x); 48 49 pure: 50 /*********************************** 51 * Returns cosine of x. x is in radians. 52 * 53 * $(TABLE_SV 54 * $(TR $(TH x) $(TH cos(x)) $(TH invalid?)) 55 * $(TR $(TD $(NAN)) $(TD $(NAN)) $(TD yes) ) 56 * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(NAN)) $(TD yes) ) 57 * ) 58 * Bugs: 59 * Results are undefined if |x| >= $(POWER 2,64). 60 */ 61 62 version (LDC) 63 { 64 alias cos = llvm_cos!float; 65 alias cos = llvm_cos!double; 66 alias cos = llvm_cos!real; 67 } 68 else 69 { 70 float cos(float x); /* intrinsic */ 71 double cos(double x); /* intrinsic */ /// ditto 72 real cos(real x); /* intrinsic */ /// ditto 73 } 74 75 /*********************************** 76 * Returns sine of x. x is in radians. 77 * 78 * $(TABLE_SV 79 * $(TR $(TH x) $(TH sin(x)) $(TH invalid?)) 80 * $(TR $(TD $(NAN)) $(TD $(NAN)) $(TD yes)) 81 * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no)) 82 * $(TR $(TD $(PLUSMNINF)) $(TD $(NAN)) $(TD yes)) 83 * ) 84 * Bugs: 85 * Results are undefined if |x| >= $(POWER 2,64). 86 */ 87 88 version (LDC) 89 { 90 alias sin = llvm_sin!float; 91 alias sin = llvm_sin!double; 92 alias sin = llvm_sin!real; 93 } 94 else 95 { 96 float sin(float x); /* intrinsic */ 97 double sin(double x); /* intrinsic */ /// ditto 98 real sin(real x); /* intrinsic */ /// ditto 99 } 100 101 /***************************************** 102 * Returns x rounded to a long value using the current rounding mode. 103 * If the integer value of x is 104 * greater than long.max, the result is 105 * indeterminate. 106 */ 107 108 version (LDC) 109 { 110 private extern(C) 111 { 112 long llroundf(float x); 113 long llround(double x); 114 long llroundl(real x); 115 } 116 117 alias rndtol = llroundf; 118 alias rndtol = llround; 119 alias rndtol = llroundl; 120 } 121 else 122 { 123 long rndtol(float x); /* intrinsic */ 124 long rndtol(double x); /* intrinsic */ /// ditto 125 long rndtol(real x); /* intrinsic */ /// ditto 126 } 127 128 /*************************************** 129 * Compute square root of x. 130 * 131 * $(TABLE_SV 132 * $(TR $(TH x) $(TH sqrt(x)) $(TH invalid?)) 133 * $(TR $(TD -0.0) $(TD -0.0) $(TD no)) 134 * $(TR $(TD $(LT)0.0) $(TD $(NAN)) $(TD yes)) 135 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no)) 136 * ) 137 */ 138 139 version (LDC) 140 { 141 pragma(inline, true): 142 143 // http://llvm.org/docs/LangRef.html#llvm-sqrt-intrinsic 144 // sqrt(x) when x is less than zero is undefined 145 float sqrt(float x) { return x < 0 ? float.nan : llvm_sqrt(x); } 146 double sqrt(double x) { return x < 0 ? double.nan : llvm_sqrt(x); } 147 real sqrt(real x) { return x < 0 ? real.nan : llvm_sqrt(x); } 148 } 149 else 150 { 151 float sqrt(float x); /* intrinsic */ 152 double sqrt(double x); /* intrinsic */ /// ditto 153 real sqrt(real x); /* intrinsic */ /// ditto 154 } 155 156 /******************************************* 157 * Compute n * 2$(SUPERSCRIPT exp) 158 * References: frexp 159 */ 160 161 version (LDC) 162 { 163 pragma(inline, true): 164 165 // Implementation from libmir: 166 // https://github.com/libmir/mir-core/blob/master/source/mir/math/ieee.d 167 private T ldexpImpl(T)(const T n, int exp) @trusted pure nothrow 168 { 169 enum RealFormat { ieeeSingle, ieeeDouble, ieeeExtended, ieeeQuadruple } 170 171 static if (T.mant_dig == 24) enum realFormat = RealFormat.ieeeSingle; 172 else static if (T.mant_dig == 53) enum realFormat = RealFormat.ieeeDouble; 173 else static if (T.mant_dig == 64) enum realFormat = RealFormat.ieeeExtended; 174 else static if (T.mant_dig == 113) enum realFormat = RealFormat.ieeeQuadruple; 175 else static assert(false, "Unsupported format for " ~ T.stringof); 176 177 version (LittleEndian) 178 { 179 enum MANTISSA_LSB = 0; 180 enum MANTISSA_MSB = 1; 181 } 182 else 183 { 184 enum MANTISSA_LSB = 1; 185 enum MANTISSA_MSB = 0; 186 } 187 188 static if (realFormat == RealFormat.ieeeExtended) 189 { 190 alias S = int; 191 alias U = ushort; 192 enum sig_mask = U(1) << (U.sizeof * 8 - 1); 193 enum exp_shft = 0; 194 enum man_mask = 0; 195 version (LittleEndian) 196 enum idx = 4; 197 else 198 enum idx = 0; 199 } 200 else 201 { 202 static if (realFormat == RealFormat.ieeeQuadruple || realFormat == RealFormat.ieeeDouble && double.sizeof == size_t.sizeof) 203 { 204 alias S = long; 205 alias U = ulong; 206 } 207 else 208 { 209 alias S = int; 210 alias U = uint; 211 } 212 static if (realFormat == RealFormat.ieeeQuadruple) 213 alias M = ulong; 214 else 215 alias M = U; 216 enum sig_mask = U(1) << (U.sizeof * 8 - 1); 217 enum uint exp_shft = T.mant_dig - 1 - (T.sizeof > U.sizeof ? U.sizeof * 8 : 0); 218 enum man_mask = (U(1) << exp_shft) - 1; 219 enum idx = T.sizeof > U.sizeof ? MANTISSA_MSB : 0; 220 } 221 enum exp_mask = (U.max >> (exp_shft + 1)) << exp_shft; 222 enum int exp_msh = exp_mask >> exp_shft; 223 enum intPartMask = man_mask + 1; 224 225 import core.checkedint : adds; 226 alias _expect = llvm_expect; 227 228 enum norm_factor = 1 / T.epsilon; 229 T vf = n; 230 231 auto u = (cast(U*)&vf)[idx]; 232 int e = (u & exp_mask) >> exp_shft; 233 if (_expect(e != exp_msh, true)) 234 { 235 if (_expect(e == 0, false)) // subnormals input 236 { 237 bool overflow; 238 vf *= norm_factor; 239 u = (cast(U*)&vf)[idx]; 240 e = int((u & exp_mask) >> exp_shft) - (T.mant_dig - 1); 241 } 242 bool overflow; 243 exp = adds(exp, e, overflow); 244 if (_expect(overflow || exp >= exp_msh, false)) // infs 245 { 246 static if (realFormat == RealFormat.ieeeExtended) 247 { 248 return vf * T.infinity; 249 } 250 else 251 { 252 u &= sig_mask; 253 u ^= exp_mask; 254 static if (realFormat == RealFormat.ieeeExtended) 255 { 256 version (LittleEndian) 257 auto mp = cast(ulong*)&vf; 258 else 259 auto mp = cast(ulong*)((cast(ushort*)&vf) + 1); 260 *mp = 0; 261 } 262 else 263 static if (T.sizeof > U.sizeof) 264 { 265 (cast(U*)&vf)[MANTISSA_LSB] = 0; 266 } 267 } 268 } 269 else 270 if (_expect(exp > 0, true)) // normal 271 { 272 u = cast(U)((u & ~exp_mask) ^ (cast(typeof(U.init + 0))exp << exp_shft)); 273 } 274 else // subnormal output 275 { 276 exp = 1 - exp; 277 static if (realFormat != RealFormat.ieeeExtended) 278 { 279 auto m = u & man_mask; 280 if (exp > T.mant_dig) 281 { 282 exp = T.mant_dig; 283 static if (T.sizeof > U.sizeof) 284 (cast(U*)&vf)[MANTISSA_LSB] = 0; 285 } 286 } 287 u &= sig_mask; 288 static if (realFormat == RealFormat.ieeeExtended) 289 { 290 version (LittleEndian) 291 auto mp = cast(ulong*)&vf; 292 else 293 auto mp = cast(ulong*)((cast(ushort*)&vf) + 1); 294 if (exp >= ulong.sizeof * 8) 295 *mp = 0; 296 else 297 *mp >>>= exp; 298 } 299 else 300 { 301 m ^= intPartMask; 302 static if (T.sizeof > U.sizeof) 303 { 304 int exp2 = exp - int(U.sizeof) * 8; 305 if (exp2 < 0) 306 { 307 (cast(U*)&vf)[MANTISSA_LSB] = ((cast(U*)&vf)[MANTISSA_LSB] >> exp) ^ (m << (U.sizeof * 8 - exp)); 308 m >>>= exp; 309 u ^= cast(U) m; 310 } 311 else 312 { 313 exp = exp2; 314 (cast(U*)&vf)[MANTISSA_LSB] = (exp < U.sizeof * 8) ? m >> exp : 0; 315 } 316 } 317 else 318 { 319 m >>>= exp; 320 u ^= cast(U) m; 321 } 322 } 323 } 324 (cast(U*)&vf)[idx] = u; 325 } 326 return vf; 327 } 328 329 float ldexp(float n, int exp) { return ldexpImpl(n, exp); } 330 double ldexp(double n, int exp) { return ldexpImpl(n, exp); } 331 static if (isRealX87) 332 { 333 // Roughly 20% faster than ldexpImpl() on an i5-3550 CPU. 334 real ldexp(real n, int exp) 335 { 336 real r = void; 337 asm @trusted pure nothrow @nogc 338 { 339 `fildl %1 # push exp 340 fxch %%st(1) # swap ST(0) and ST(1) 341 fscale # ST(0) := ST(0) * (2 ^^ ST(1)) 342 fstp %%st(1) # pop and keep ST(0) value on top` 343 : "=st" (r) 344 : "m" (exp), "st" (n) 345 : "flags"; // might clobber x87 flags 346 } 347 return r; 348 } 349 } 350 else 351 { 352 real ldexp(real n, int exp) { return ldexpImpl(n, exp); } 353 } 354 } 355 else 356 { 357 float ldexp(float n, int exp); /* intrinsic */ 358 double ldexp(double n, int exp); /* intrinsic */ /// ditto 359 real ldexp(real n, int exp); /* intrinsic */ /// ditto 360 } 361 362 unittest { 363 static if (real.mant_dig == 113) 364 { 365 assert(ldexp(1.0L, -16384) == 0x1p-16384L); 366 assert(ldexp(1.0L, -16382) == 0x1p-16382L); 367 } 368 else static if (real.mant_dig == 106) 369 { 370 assert(ldexp(1.0L, 1023) == 0x1p1023L); 371 assert(ldexp(1.0L, -1022) == 0x1p-1022L); 372 assert(ldexp(1.0L, -1021) == 0x1p-1021L); 373 } 374 else static if (real.mant_dig == 64) 375 { 376 assert(ldexp(1.0L, -16384) == 0x1p-16384L); 377 assert(ldexp(1.0L, -16382) == 0x1p-16382L); 378 } 379 else static if (real.mant_dig == 53) 380 { 381 assert(ldexp(1.0L, 1023) == 0x1p1023L); 382 assert(ldexp(1.0L, -1022) == 0x1p-1022L); 383 assert(ldexp(1.0L, -1021) == 0x1p-1021L); 384 } 385 else 386 assert(false, "Only 128bit, 80bit and 64bit reals expected here"); 387 } 388 389 /******************************* 390 * Compute the absolute value. 391 * $(TABLE_SV 392 * $(TR $(TH x) $(TH fabs(x))) 393 * $(TR $(TD $(PLUSMN)0.0) $(TD +0.0) ) 394 * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD +$(INFIN)) ) 395 * ) 396 * It is implemented as a compiler intrinsic. 397 * Params: 398 * x = floating point value 399 * Returns: |x| 400 * References: equivalent to `std.math.fabs` 401 */ 402 version (LDC) 403 { 404 alias fabs = llvm_fabs!float; 405 alias fabs = llvm_fabs!double; 406 alias fabs = llvm_fabs!real; 407 } 408 else @safe pure nothrow @nogc 409 { 410 float fabs(float x); 411 double fabs(double x); /// ditto 412 real fabs(real x); /// ditto 413 } 414 415 /********************************** 416 * Rounds x to the nearest integer value, using the current rounding 417 * mode. 418 * If the return value is not equal to x, the FE_INEXACT 419 * exception is raised. 420 * $(B nearbyint) performs 421 * the same operation, but does not set the FE_INEXACT exception. 422 */ 423 version (LDC) 424 { 425 alias rint = llvm_rint!float; 426 alias rint = llvm_rint!double; 427 alias rint = llvm_rint!real; 428 } 429 else 430 { 431 float rint(float x); /* intrinsic */ 432 double rint(double x); /* intrinsic */ /// ditto 433 real rint(real x); /* intrinsic */ /// ditto 434 } 435 436 /*********************************** 437 * Building block functions, they 438 * translate to a single x87 instruction. 439 */ 440 441 version (LDC) 442 { 443 static if (isRealX87) 444 { 445 pragma(inline, true): 446 447 // y * log2(x) 448 real yl2x(real x, real y) 449 { 450 real r = void; 451 asm @trusted pure nothrow @nogc { "fyl2x" : "=st" (r) : "st(1)" (y), "st" (x) : "st(1)", "flags"; } 452 return r; 453 } 454 455 // y * log2(x + 1) 456 real yl2xp1(real x, real y) 457 { 458 real r = void; 459 asm @trusted pure nothrow @nogc { "fyl2xp1" : "=st" (r) : "st(1)" (y), "st" (x) : "st(1)", "flags"; } 460 return r; 461 } 462 } 463 } 464 else 465 { 466 // y * log2(x) 467 float yl2x(float x, float y); /* intrinsic */ 468 double yl2x(double x, double y); /* intrinsic */ /// ditto 469 real yl2x(real x, real y); /* intrinsic */ /// ditto 470 // y * log2(x +1) 471 float yl2xp1(float x, float y); /* intrinsic */ 472 double yl2xp1(double x, double y); /* intrinsic */ /// ditto 473 real yl2xp1(real x, real y); /* intrinsic */ /// ditto 474 } 475 476 unittest 477 { 478 version (INLINE_YL2X) 479 { 480 assert(yl2x(1024.0L, 1) == 10); 481 assert(yl2xp1(1023.0L, 1) == 10); 482 } 483 } 484 485 /************************************* 486 * Round argument to a specific precision. 487 * 488 * D language types specify only a minimum precision, not a maximum. The 489 * `toPrec()` function forces rounding of the argument `f` to the precision 490 * of the specified floating point type `T`. 491 * The rounding mode used is inevitably target-dependent, but will be done in 492 * a way to maximize accuracy. In most cases, the default is round-to-nearest. 493 * 494 * Params: 495 * T = precision type to round to 496 * f = value to convert 497 * Returns: 498 * f in precision of type `T` 499 */ 500 T toPrec(T:float)(float f) { pragma(inline, false); return f; } 501 /// ditto 502 T toPrec(T:float)(double f) { pragma(inline, false); return cast(T) f; } 503 /// ditto 504 T toPrec(T:float)(real f) { pragma(inline, false); return cast(T) f; } 505 /// ditto 506 T toPrec(T:double)(float f) { pragma(inline, false); return f; } 507 /// ditto 508 T toPrec(T:double)(double f) { pragma(inline, false); return f; } 509 /// ditto 510 T toPrec(T:double)(real f) { pragma(inline, false); return cast(T) f; } 511 /// ditto 512 T toPrec(T:real)(float f) { pragma(inline, false); return f; } 513 /// ditto 514 T toPrec(T:real)(double f) { pragma(inline, false); return f; } 515 /// ditto 516 T toPrec(T:real)(real f) { pragma(inline, false); return f; } 517 518 @safe unittest 519 { 520 // Test all instantiations work with all combinations of float. 521 float f = 1.1f; 522 double d = 1.1; 523 real r = 1.1L; 524 f = toPrec!float(f + f); 525 f = toPrec!float(d + d); 526 f = toPrec!float(r + r); 527 d = toPrec!double(f + f); 528 d = toPrec!double(d + d); 529 d = toPrec!double(r + r); 530 r = toPrec!real(f + f); 531 r = toPrec!real(d + d); 532 r = toPrec!real(r + r); 533 534 // Comparison tests. 535 bool approxEqual(T)(T lhs, T rhs) 536 { 537 return fabs((lhs - rhs) / rhs) <= 1e-2 || fabs(lhs - rhs) <= 1e-5; 538 } 539 540 enum real PIR = 0xc.90fdaa22168c235p-2; 541 enum double PID = 0x1.921fb54442d18p+1; 542 enum float PIF = 0x1.921fb6p+1; 543 static assert(approxEqual(toPrec!float(PIR), PIF)); 544 static assert(approxEqual(toPrec!double(PIR), PID)); 545 static assert(approxEqual(toPrec!real(PIR), PIR)); 546 static assert(approxEqual(toPrec!float(PID), PIF)); 547 static assert(approxEqual(toPrec!double(PID), PID)); 548 static assert(approxEqual(toPrec!real(PID), PID)); 549 static assert(approxEqual(toPrec!float(PIF), PIF)); 550 static assert(approxEqual(toPrec!double(PIF), PIF)); 551 static assert(approxEqual(toPrec!real(PIF), PIF)); 552 553 assert(approxEqual(toPrec!float(PIR), PIF)); 554 assert(approxEqual(toPrec!double(PIR), PID)); 555 assert(approxEqual(toPrec!real(PIR), PIR)); 556 assert(approxEqual(toPrec!float(PID), PIF)); 557 assert(approxEqual(toPrec!double(PID), PID)); 558 assert(approxEqual(toPrec!real(PID), PID)); 559 assert(approxEqual(toPrec!float(PIF), PIF)); 560 assert(approxEqual(toPrec!double(PIF), PIF)); 561 assert(approxEqual(toPrec!real(PIF), PIF)); 562 }