Mercurial > projects > ldc
annotate tango/tango/math/Math.d @ 164:a64becf2a702 trunk
[svn r180] Fixed complex negation, and tango.math.Math now compiles.
| author | lindquist |
|---|---|
| date | Mon, 05 May 2008 20:28:59 +0200 |
| parents | 1700239cab2e |
| children |
| rev | line source |
|---|---|
| 132 | 1 /** |
| 2 * Elementary Mathematical Functions | |
| 3 * | |
| 4 * Copyright: Portions Copyright (C) 2001-2005 Digital Mars. | |
| 5 * License: BSD style: $(LICENSE), Digital Mars. | |
| 6 * Authors: Walter Bright, Don Clugston, Sean Kelly | |
| 7 */ | |
| 8 /* Portions of this code were taken from Phobos std.math, which has the following | |
| 9 * copyright notice: | |
| 10 * | |
| 11 * Author: | |
| 12 * Walter Bright | |
| 13 * Copyright: | |
| 14 * Copyright (c) 2001-2005 by Digital Mars, | |
| 15 * All Rights Reserved, | |
| 16 * www.digitalmars.com | |
| 17 * License: | |
| 18 * This software is provided 'as-is', without any express or implied | |
| 19 * warranty. In no event will the authors be held liable for any damages | |
| 20 * arising from the use of this software. | |
| 21 * | |
| 22 * Permission is granted to anyone to use this software for any purpose, | |
| 23 * including commercial applications, and to alter it and redistribute it | |
| 24 * freely, subject to the following restrictions: | |
| 25 * | |
| 26 * <ul> | |
| 27 * <li> The origin of this software must not be misrepresented; you must not | |
| 28 * claim that you wrote the original software. If you use this software | |
| 29 * in a product, an acknowledgment in the product documentation would be | |
| 30 * appreciated but is not required. | |
| 31 * </li> | |
| 32 * <li> Altered source versions must be plainly marked as such, and must not | |
| 33 * be misrepresented as being the original software. | |
| 34 * </li> | |
| 35 * <li> This notice may not be removed or altered from any source | |
| 36 * distribution. | |
| 37 * </li> | |
| 38 * </ul> | |
| 39 */ | |
| 40 | |
| 41 /** | |
| 42 * Macros: | |
| 43 * NAN = $(RED NAN) | |
| 44 * TEXTNAN = $(RED NAN:$1 ) | |
| 45 * SUP = <span style="vertical-align:super;font-size:smaller">$0</span> | |
| 46 * GAMMA = Γ | |
| 47 * INTEGRAL = ∫ | |
| 48 * INTEGRATE = $(BIG ∫<sub>$(SMALL $1)</sub><sup>$2</sup>) | |
| 49 * POWER = $1<sup>$2</sup> | |
| 50 * BIGSUM = $(BIG Σ <sup>$2</sup><sub>$(SMALL $1)</sub>) | |
| 51 * CHOOSE = $(BIG () <sup>$(SMALL $1)</sup><sub>$(SMALL $2)</sub> $(BIG )) | |
| 52 * TABLE_SV = <table border=1 cellpadding=4 cellspacing=0> | |
| 53 * <caption>Special Values</caption> | |
| 54 * $0</table> | |
| 55 * SVH = $(TR $(TH $1) $(TH $2)) | |
| 56 * SV = $(TR $(TD $1) $(TD $2)) | |
| 57 * TABLE_DOMRG = <table border=1 cellpadding=4 cellspacing=0>$0</table> | |
| 58 * DOMAIN = $(TR $(TD Domain) $(TD $0)) | |
| 59 * RANGE = $(TR $(TD Range) $(TD $0)) | |
| 60 */ | |
| 61 | |
| 62 module tango.math.Math; | |
| 63 | |
| 64 static import tango.stdc.math; | |
| 65 private import tango.math.IEEE; | |
| 66 | |
| 67 version(DigitalMars) | |
| 68 { | |
| 69 version(D_InlineAsm_X86) | |
| 70 { | |
| 71 version = DigitalMars_D_InlineAsm_X86; | |
| 72 } | |
| 73 } | |
| 74 | |
| 75 /* | |
| 76 * Constants | |
| 77 */ | |
| 78 | |
| 79 const real E = 2.7182818284590452354L; /** e */ | |
| 80 const real LOG2T = 0x1.a934f0979a3715fcp+1; /** log<sub>2</sub>10 */ // 3.32193 fldl2t | |
| 81 const real LOG2E = 0x1.71547652b82fe178p+0; /** log<sub>2</sub>e */ // 1.4427 fldl2e | |
| 82 const real LOG2 = 0x1.34413509f79fef32p-2; /** log<sub>10</sub>2 */ // 0.30103 fldlg2 | |
| 83 const real LOG10E = 0.43429448190325182765L; /** log<sub>10</sub>e */ | |
| 84 const real LN2 = 0x1.62e42fefa39ef358p-1; /** ln 2 */ // 0.693147 fldln2 | |
| 85 const real LN10 = 2.30258509299404568402L; /** ln 10 */ | |
| 86 const real PI = 0x1.921fb54442d1846ap+1; /** π */ // 3.14159 fldpi | |
| 87 const real PI_2 = 1.57079632679489661923L; /** π / 2 */ | |
| 88 const real PI_4 = 0.78539816339744830962L; /** π / 4 */ | |
| 89 const real M_1_PI = 0.31830988618379067154L; /** 1 / π */ | |
| 90 const real M_2_PI = 0.63661977236758134308L; /** 2 / π */ | |
| 91 const real M_2_SQRTPI = 1.12837916709551257390L; /** 2 / √π */ | |
| 92 const real SQRT2 = 1.41421356237309504880L; /** √2 */ | |
| 93 const real SQRT1_2 = 0.70710678118654752440L; /** √½ */ | |
| 94 | |
| 95 //const real SQRTPI = 1.77245385090551602729816748334114518279754945612238L; /** √π */ | |
| 96 //const real SQRT2PI = 2.50662827463100050242E0L; /** √(2 π) */ | |
| 97 | |
| 98 const real MAXLOG = 0x1.62e42fefa39ef358p+13; /** log(real.max) */ | |
| 99 const real MINLOG = -0x1.6436716d5406e6d8p+13; /** log(real.min*real.epsilon) */ | |
| 100 const real EULERGAMMA = 0.57721_56649_01532_86060_65120_90082_40243_10421_59335_93992; /** Euler-Mascheroni constant 0.57721566.. */ | |
| 101 | |
| 102 /* | |
| 103 * Primitives | |
| 104 */ | |
| 105 | |
| 106 /** | |
| 107 * Calculates the absolute value | |
| 108 * | |
| 109 * For complex numbers, abs(z) = sqrt( $(POWER z.re, 2) + $(POWER z.im, 2) ) | |
| 110 * = hypot(z.re, z.im). | |
| 111 */ | |
| 112 real abs(real x) | |
| 113 { | |
| 114 return tango.math.IEEE.fabs(x); | |
| 115 } | |
| 116 | |
| 117 /** ditto */ | |
| 118 long abs(long x) | |
| 119 { | |
| 120 return x>=0 ? x : -x; | |
| 121 } | |
| 122 | |
| 123 /** ditto */ | |
| 124 int abs(int x) | |
| 125 { | |
| 126 return x>=0 ? x : -x; | |
| 127 } | |
| 128 | |
| 129 /** ditto */ | |
| 130 real abs(creal z) | |
| 131 { | |
| 132 return hypot(z.re, z.im); | |
| 133 } | |
| 134 | |
| 135 /** ditto */ | |
| 136 real abs(ireal y) | |
| 137 { | |
| 138 return tango.math.IEEE.fabs(y.im); | |
| 139 } | |
| 140 | |
| 141 debug(UnitTest) { | |
| 142 unittest | |
| 143 { | |
| 144 assert(isIdentical(0.0L,abs(-0.0L))); | |
| 145 assert(isNaN(abs(real.nan))); | |
| 146 assert(abs(-real.infinity) == real.infinity); | |
| 147 assert(abs(-3.2Li) == 3.2L); | |
| 148 assert(abs(71.6Li) == 71.6L); | |
| 149 assert(abs(-56) == 56); | |
| 150 assert(abs(2321312L) == 2321312L); | |
| 151 assert(abs(-1.0L+1.0Li) == sqrt(2.0L)); | |
| 152 } | |
| 153 } | |
| 154 | |
| 155 /** | |
| 156 * Complex conjugate | |
| 157 * | |
| 158 * conj(x + iy) = x - iy | |
| 159 * | |
| 160 * Note that z * conj(z) = $(POWER z.re, 2) + $(POWER z.im, 2) | |
| 161 * is always a real number | |
| 162 */ | |
| 163 creal conj(creal z) | |
| 164 { | |
| 165 return z.re - z.im*1i; | |
| 166 } | |
| 167 | |
| 168 /** ditto */ | |
| 169 ireal conj(ireal y) | |
| 170 { | |
| 171 return -y; | |
| 172 } | |
| 173 | |
| 174 debug(UnitTest) { | |
| 175 unittest | |
| 176 { | |
| 177 assert(conj(7 + 3i) == 7-3i); | |
| 178 ireal z = -3.2Li; | |
| 179 assert(conj(z) == -z); | |
| 180 } | |
| 181 } | |
| 182 | |
| 183 private { | |
| 184 // Return the type which would be returned by a max or min operation | |
| 185 template minmaxtype(T...){ | |
| 186 static if(T.length == 1) alias typeof(T[0]) minmaxtype; | |
| 187 else static if(T.length > 2) | |
| 188 alias minmaxtype!(minmaxtype!(T[0..2]), T[2..$]) minmaxtype; | |
| 189 else alias typeof (T[1] > T[0] ? T[1] : T[0]) minmaxtype; | |
| 190 } | |
| 191 } | |
| 192 | |
| 193 /** Return the minimum of the supplied arguments. | |
| 194 * | |
| 195 * Note: If the arguments are floating-point numbers, and at least one is a NaN, | |
| 196 * the result is undefined. | |
| 197 */ | |
| 198 minmaxtype!(T) min(T...)(T arg){ | |
| 199 static if(arg.length == 1) return arg[0]; | |
| 200 else static if(arg.length == 2) return arg[1] < arg[0] ? arg[1] : arg[0]; | |
| 201 static if(arg.length > 2) return min(arg[1] < arg[0] ? arg[1] : arg[0], arg[2..$]); | |
| 202 } | |
| 203 | |
| 204 /** Return the maximum of the supplied arguments. | |
| 205 * | |
| 206 * Note: If the arguments are floating-point numbers, and at least one is a NaN, | |
| 207 * the result is undefined. | |
| 208 */ | |
| 209 minmaxtype!(T) max(T...)(T arg){ | |
| 210 static if(arg.length == 1) return arg[0]; | |
| 211 else static if(arg.length == 2) return arg[1] > arg[0] ? arg[1] : arg[0]; | |
| 212 static if(arg.length > 2) return max(arg[1] > arg[0] ? arg[1] : arg[0], arg[2..$]); | |
| 213 } | |
| 214 debug(UnitTest) { | |
| 215 unittest | |
| 216 { | |
| 217 assert(max('e', 'f')=='f'); | |
| 218 assert(min(3.5, 3.8)==3.5); | |
| 219 // check implicit conversion to integer. | |
| 220 assert(min(3.5, 18)==3.5); | |
| 221 | |
| 222 } | |
| 223 } | |
| 224 | |
| 225 /** Returns the minimum number of x and y, favouring numbers over NaNs. | |
| 226 * | |
| 227 * If both x and y are numbers, the minimum is returned. | |
| 228 * If both parameters are NaN, either will be returned. | |
| 229 * If one parameter is a NaN and the other is a number, the number is | |
| 230 * returned (this behaviour is mandated by IEEE 754R, and is useful | |
| 231 * for determining the range of a function). | |
| 232 */ | |
| 233 real minNum(real x, real y) { | |
| 234 if (x<=y || isNaN(y)) return x; else return y; | |
| 235 } | |
| 236 | |
| 237 /** Returns the maximum number of x and y, favouring numbers over NaNs. | |
| 238 * | |
| 239 * If both x and y are numbers, the maximum is returned. | |
| 240 * If both parameters are NaN, either will be returned. | |
| 241 * If one parameter is a NaN and the other is a number, the number is | |
| 242 * returned (this behaviour is mandated by IEEE 754R, and is useful | |
| 243 * for determining the range of a function). | |
| 244 */ | |
| 245 real maxNum(real x, real y) { | |
| 246 if (x>=y || isNaN(y)) return x; else return y; | |
| 247 } | |
| 248 | |
| 249 /** Returns the minimum of x and y, favouring NaNs over numbers | |
| 250 * | |
| 251 * If both x and y are numbers, the minimum is returned. | |
| 252 * If both parameters are NaN, either will be returned. | |
| 253 * If one parameter is a NaN and the other is a number, the NaN is returned. | |
| 254 */ | |
| 255 real minNaN(real x, real y) { | |
| 256 return (x<=y || isNaN(x))? x : y; | |
| 257 } | |
| 258 | |
| 259 /** Returns the maximum of x and y, favouring NaNs over numbers | |
| 260 * | |
| 261 * If both x and y are numbers, the maximum is returned. | |
| 262 * If both parameters are NaN, either will be returned. | |
| 263 * If one parameter is a NaN and the other is a number, the NaN is returned. | |
| 264 */ | |
| 265 real maxNaN(real x, real y) { | |
| 266 return (x>=y || isNaN(x))? x : y; | |
| 267 } | |
| 268 | |
| 269 debug(UnitTest) { | |
| 270 unittest | |
| 271 { | |
| 272 assert(maxNum(NaN(0xABC), 56.1L)== 56.1L); | |
| 273 assert(isIdentical(maxNaN(NaN(1389), 56.1L), NaN(1389))); | |
| 274 assert(maxNum(28.0, NaN(0xABC))== 28.0); | |
| 275 assert(minNum(1e12, NaN(0xABC))== 1e12); | |
| 276 assert(isIdentical(minNaN(1e12, NaN(23454)), NaN(23454))); | |
| 277 assert(isIdentical(minNum(NaN(489), NaN(23)), NaN(489))); | |
| 278 } | |
| 279 } | |
| 280 | |
| 281 /* | |
| 282 * Trig Functions | |
| 283 */ | |
| 284 | |
| 285 /** | |
| 286 * Returns cosine of x. x is in radians. | |
| 287 * | |
| 288 * $(TABLE_SV | |
| 289 * $(TR $(TH x) $(TH cos(x)) $(TH invalid?) ) | |
| 290 * $(TR $(TD $(NAN)) $(TD $(NAN)) $(TD yes) ) | |
| 291 * $(TR $(TD ±∞) $(TD $(NAN)) $(TD yes) ) | |
| 292 * ) | |
| 293 * Bugs: | |
| 294 * Results are undefined if |x| >= $(POWER 2,64). | |
| 295 */ | |
| 296 real cos(real x) /* intrinsic */ | |
| 297 { | |
| 298 version(D_InlineAsm_X86) | |
| 299 { | |
| 300 asm | |
| 301 { | |
| 302 fld x; | |
| 303 fcos; | |
| 304 } | |
| 305 } | |
| 306 else | |
| 307 { | |
| 308 return tango.stdc.math.cosl(x); | |
| 309 } | |
| 310 } | |
| 311 | |
| 312 debug(UnitTest) { | |
| 313 unittest { | |
| 314 // NaN payloads | |
| 315 assert(isIdentical(cos(NaN(314)), NaN(314))); | |
| 316 } | |
| 317 } | |
| 318 | |
| 319 /** | |
| 320 * Returns sine of x. x is in radians. | |
| 321 * | |
| 322 * $(TABLE_SV | |
| 323 * <tr> <th> x <th> sin(x) <th>invalid? | |
| 324 * <tr> <td> $(NAN) <td> $(NAN) <td> yes | |
| 325 * <tr> <td> ±0.0 <td> ±0.0 <td> no | |
| 326 * <tr> <td> ±∞ <td> $(NAN) <td> yes | |
| 327 * ) | |
| 328 * Bugs: | |
| 329 * Results are undefined if |x| >= $(POWER 2,64). | |
| 330 */ | |
| 331 real sin(real x) /* intrinsic */ | |
| 332 { | |
| 333 version(D_InlineAsm_X86) | |
| 334 { | |
| 335 asm | |
| 336 { | |
| 337 fld x; | |
| 338 fsin; | |
| 339 } | |
| 340 } | |
| 341 else | |
| 342 { | |
| 343 return tango.stdc.math.sinl(x); | |
| 344 } | |
| 345 } | |
| 346 | |
| 347 debug(UnitTest) { | |
| 348 unittest { | |
| 349 // NaN payloads | |
| 350 assert(isIdentical(sin(NaN(314)), NaN(314))); | |
| 351 } | |
| 352 } | |
| 353 | |
| 354 version (GNU) { | |
| 355 extern (C) real tanl(real); | |
| 356 } | |
| 357 | |
| 358 /** | |
| 359 * Returns tangent of x. x is in radians. | |
| 360 * | |
| 361 * $(TABLE_SV | |
| 362 * <tr> <th> x <th> tan(x) <th> invalid? | |
| 363 * <tr> <td> $(NAN) <td> $(NAN) <td> yes | |
| 364 * <tr> <td> ±0.0 <td> ±0.0 <td> no | |
| 365 * <tr> <td> ±∞ <td> $(NAN) <td> yes | |
| 366 * ) | |
| 367 */ | |
| 368 real tan(real x) | |
| 369 { | |
| 370 version (GNU) { | |
| 371 return tanl(x); | |
|
164
a64becf2a702
[svn r180] Fixed complex negation, and tango.math.Math now compiles.
lindquist
parents:
132
diff
changeset
|
372 } else version (LLVMDC) { |
|
a64becf2a702
[svn r180] Fixed complex negation, and tango.math.Math now compiles.
lindquist
parents:
132
diff
changeset
|
373 return tango.stdc.math.tanl(x); |
| 132 | 374 } else { |
| 375 asm | |
| 376 { | |
| 377 fld x[EBP] ; // load theta | |
| 378 fxam ; // test for oddball values | |
| 379 fstsw AX ; | |
| 380 sahf ; | |
| 381 jc trigerr ; // x is NAN, infinity, or empty | |
| 382 // 387's can handle denormals | |
| 383 SC18: fptan ; | |
| 384 fstp ST(0) ; // dump X, which is always 1 | |
| 385 fstsw AX ; | |
| 386 sahf ; | |
| 387 jnp Lret ; // C2 = 1 (x is out of range) | |
| 388 | |
| 389 // Do argument reduction to bring x into range | |
| 390 fldpi ; | |
| 391 fxch ; | |
| 392 SC17: fprem1 ; | |
| 393 fstsw AX ; | |
| 394 sahf ; | |
| 395 jp SC17 ; | |
| 396 fstp ST(1) ; // remove pi from stack | |
| 397 jmp SC18 ; | |
| 398 | |
| 399 trigerr: | |
| 400 jnp Lret ; // if x is NaN, return x. | |
| 401 fstp ST(0) ; // dump x, which will be infinity | |
| 402 } | |
| 403 return NaN(TANGO_NAN.TAN_DOMAIN); | |
| 404 Lret: | |
| 405 ; | |
| 406 } | |
| 407 } | |
| 408 | |
| 409 debug(UnitTest) { | |
| 410 unittest | |
| 411 { | |
| 412 // Returns true if equal to precision, false if not | |
| 413 // (Used only in unit test for tan()) | |
| 414 bool mfeq(real x, real y, real precision) | |
| 415 { | |
| 416 if (x == y) | |
| 417 return true; | |
| 418 if (isNaN(x) || isNaN(y)) | |
| 419 return false; | |
| 420 return fabs(x - y) <= precision; | |
| 421 } | |
| 422 | |
| 423 | |
| 424 static real vals[][2] = // angle,tan | |
| 425 [ | |
| 426 [ 0, 0], | |
| 427 [ .5, .5463024898], | |
| 428 [ 1, 1.557407725], | |
| 429 [ 1.5, 14.10141995], | |
| 430 [ 2, -2.185039863], | |
| 431 [ 2.5,-.7470222972], | |
| 432 [ 3, -.1425465431], | |
| 433 [ 3.5, .3745856402], | |
| 434 [ 4, 1.157821282], | |
| 435 [ 4.5, 4.637332055], | |
| 436 [ 5, -3.380515006], | |
| 437 [ 5.5,-.9955840522], | |
| 438 [ 6, -.2910061914], | |
| 439 [ 6.5, .2202772003], | |
| 440 [ 10, .6483608275], | |
| 441 | |
| 442 // special angles | |
| 443 [ PI_4, 1], | |
| 444 //[ PI_2, real.infinity], | |
| 445 [ 3*PI_4, -1], | |
| 446 [ PI, 0], | |
| 447 [ 5*PI_4, 1], | |
| 448 //[ 3*PI_2, -real.infinity], | |
| 449 [ 7*PI_4, -1], | |
| 450 [ 2*PI, 0], | |
| 451 | |
| 452 // overflow | |
| 453 [ real.infinity, real.nan], | |
| 454 [ real.nan, real.nan], | |
| 455 ]; | |
| 456 int i; | |
| 457 | |
| 458 for (i = 0; i < vals.length; i++) | |
| 459 { | |
| 460 real x = vals[i][0]; | |
| 461 real r = vals[i][1]; | |
| 462 real t = tan(x); | |
| 463 | |
| 464 //printf("tan(%Lg) = %Lg, should be %Lg\n", x, t, r); | |
| 465 if (isNaN(r)) assert(isNaN(t)); | |
| 466 else assert(mfeq(r, t, .0000001)); | |
| 467 | |
| 468 x = -x; | |
| 469 r = -r; | |
| 470 t = tan(x); | |
| 471 //printf("tan(%Lg) = %Lg, should be %Lg\n", x, t, r); | |
| 472 if (isNaN(r)) assert(isNaN(t)); | |
| 473 else assert(mfeq(r, t, .0000001)); | |
| 474 } | |
| 475 assert(isIdentical(tan(NaN(735)), NaN(735))); | |
| 476 assert(isNaN(tan(real.infinity))); | |
| 477 } | |
| 478 } | |
| 479 | |
| 480 /***************************************** | |
| 481 * Sine, cosine, and arctangent of multiple of π | |
| 482 * | |
| 483 * Accuracy is preserved for large values of x. | |
| 484 */ | |
| 485 real cosPi(real x) | |
| 486 { | |
| 487 return cos((x%2.0)*PI); | |
| 488 } | |
| 489 | |
| 490 /** ditto */ | |
| 491 real sinPi(real x) | |
| 492 { | |
| 493 return sin((x%2.0)*PI); | |
| 494 } | |
| 495 | |
| 496 /** ditto */ | |
| 497 real atanPi(real x) | |
| 498 { | |
| 499 return PI * atan(x); // BUG: Fix this. | |
| 500 } | |
| 501 | |
| 502 debug(UnitTest) { | |
| 503 unittest { | |
| 504 assert(isIdentical(sinPi(0.0), 0.0)); | |
| 505 assert(isIdentical(sinPi(-0.0), -0.0)); | |
| 506 assert(isIdentical(atanPi(0.0), 0.0)); | |
| 507 assert(isIdentical(atanPi(-0.0), -0.0)); | |
| 508 } | |
| 509 } | |
| 510 | |
| 511 /*********************************** | |
| 512 * sine, complex and imaginary | |
| 513 * | |
| 514 * sin(z) = sin(z.re)*cosh(z.im) + cos(z.re)*sinh(z.im)i | |
| 515 * | |
| 516 * If both sin(θ) and cos(θ) are required, | |
| 517 * it is most efficient to use expi(&theta). | |
| 518 */ | |
| 519 creal sin(creal z) | |
| 520 { | |
| 521 creal cs = expi(z.re); | |
| 522 return cs.im * cosh(z.im) + cs.re * sinh(z.im) * 1i; | |
| 523 } | |
| 524 | |
| 525 /** ditto */ | |
| 526 ireal sin(ireal y) | |
| 527 { | |
| 528 return cosh(y.im)*1i; | |
| 529 } | |
| 530 | |
| 531 debug(UnitTest) { | |
| 532 unittest { | |
| 533 assert(sin(0.0+0.0i) == 0.0); | |
| 534 assert(sin(2.0+0.0i) == sin(2.0L) ); | |
| 535 } | |
| 536 } | |
| 537 | |
| 538 /*********************************** | |
| 539 * cosine, complex and imaginary | |
| 540 * | |
| 541 * cos(z) = cos(z.re)*cosh(z.im) + sin(z.re)*sinh(z.im)i | |
| 542 */ | |
| 543 creal cos(creal z) | |
| 544 { | |
| 545 creal cs = expi(z.re); | |
| 546 return cs.re * cosh(z.im) + cs.im * sinh(z.im) * 1i; | |
| 547 } | |
| 548 | |
| 549 /** ditto */ | |
| 550 real cos(ireal y) | |
| 551 { | |
| 552 return cosh(y.im); | |
| 553 } | |
| 554 | |
| 555 debug(UnitTest) { | |
| 556 unittest{ | |
| 557 assert(cos(0.0+0.0i)==1.0); | |
| 558 assert(cos(1.3L+0.0i)==cos(1.3L)); | |
| 559 assert(cos(5.2Li)== cosh(5.2L)); | |
| 560 } | |
| 561 } | |
| 562 | |
| 563 /** | |
| 564 * Calculates the arc cosine of x, | |
| 565 * returning a value ranging from -π/2 to π/2. | |
| 566 * | |
| 567 * $(TABLE_SV | |
| 568 * <tr> <th> x <th> acos(x) <th> invalid? | |
| 569 * <tr> <td> >1.0 <td> $(NAN) <td> yes | |
| 570 * <tr> <td> <-1.0 <td> $(NAN) <td> yes | |
| 571 * <tr> <td> $(NAN) <td> $(NAN) <td> yes | |
| 572 * ) | |
| 573 */ | |
| 574 real acos(real x) | |
| 575 { | |
| 576 return tango.stdc.math.acosl(x); | |
| 577 } | |
| 578 | |
| 579 debug(UnitTest) { | |
| 580 unittest { | |
| 581 // NaN payloads | |
| 582 assert(isIdentical(acos(NaN(254)), NaN(254))); | |
| 583 } | |
| 584 } | |
| 585 | |
| 586 /** | |
| 587 * Calculates the arc sine of x, | |
| 588 * returning a value ranging from -π/2 to π/2. | |
| 589 * | |
| 590 * $(TABLE_SV | |
| 591 * <tr> <th> x <th> asin(x) <th> invalid? | |
| 592 * <tr> <td> ±0.0 <td> ±0.0 <td> no | |
| 593 * <tr> <td> >1.0 <td> $(NAN) <td> yes | |
| 594 * <tr> <td> <-1.0 <td> $(NAN) <td> yes | |
| 595 * ) | |
| 596 */ | |
| 597 real asin(real x) | |
| 598 { | |
| 599 return tango.stdc.math.asinl(x); | |
| 600 } | |
| 601 | |
| 602 debug(UnitTest) { | |
| 603 unittest { | |
| 604 // NaN payloads | |
| 605 assert(isIdentical(asin(NaN(7249)), NaN(7249))); | |
| 606 } | |
| 607 } | |
| 608 | |
| 609 /** | |
| 610 * Calculates the arc tangent of x, | |
| 611 * returning a value ranging from -π/2 to π/2. | |
| 612 * | |
| 613 * $(TABLE_SV | |
| 614 * <tr> <th> x <th> atan(x) <th> invalid? | |
| 615 * <tr> <td> ±0.0 <td> ±0.0 <td> no | |
| 616 * <tr> <td> ±∞ <td> $(NAN) <td> yes | |
| 617 * ) | |
| 618 */ | |
| 619 real atan(real x) | |
| 620 { | |
| 621 return tango.stdc.math.atanl(x); | |
| 622 } | |
| 623 | |
| 624 debug(UnitTest) { | |
| 625 unittest { | |
| 626 // NaN payloads | |
| 627 assert(isIdentical(atan(NaN(9876)), NaN(9876))); | |
| 628 } | |
| 629 } | |
| 630 | |
| 631 /** | |
| 632 * Calculates the arc tangent of y / x, | |
| 633 * returning a value ranging from -π/2 to π/2. | |
| 634 * | |
| 635 * Remarks: | |
| 636 * The Complex Argument of a complex number z is given by | |
| 637 * Arg(z) = atan2(z.re, z.im) | |
| 638 * | |
| 639 * $(TABLE_SV | |
| 640 * <tr> <th> y <th> x <th> atan(y, x) | |
| 641 * <tr> <td> $(NAN) <td> anything <td> $(NAN) | |
| 642 * <tr> <td> anything <td> $(NAN) <td> $(NAN) | |
| 643 * <tr> <td> ±0.0 <td> > 0.0 <td> ±0.0 | |
| 644 * <tr> <td> ±0.0 <td> ±0.0 <td> ±0.0 | |
| 645 * <tr> <td> ±0.0 <td> < 0.0 <td> ±π | |
| 646 * <tr> <td> ±0.0 <td> -0.0 <td> ±π | |
| 647 * <tr> <td> > 0.0 <td> ±0.0 <td> π/2 | |
| 648 * <tr> <td> < 0.0 <td> ±0.0 <td> π/2 | |
| 649 * <tr> <td> > 0.0 <td> ∞ <td> ±0.0 | |
| 650 * <tr> <td> ±∞ <td> anything <td> ±π/2 | |
| 651 * <tr> <td> > 0.0 <td> -∞ <td> ±π | |
| 652 * <tr> <td> ±∞ <td> ∞ <td> ±π/4 | |
| 653 * <tr> <td> ±∞ <td> -∞ <td> ±3π/4 | |
| 654 * ) | |
| 655 */ | |
| 656 real atan2(real y, real x) | |
| 657 { | |
| 658 return tango.stdc.math.atan2l(y,x); | |
| 659 } | |
| 660 | |
| 661 debug(UnitTest) { | |
| 662 unittest { | |
| 663 // NaN payloads | |
| 664 assert(isIdentical(atan2(5.3, NaN(9876)), NaN(9876))); | |
| 665 assert(isIdentical(atan2(NaN(9876), 2.18), NaN(9876))); | |
| 666 } | |
| 667 } | |
| 668 | |
| 669 /*********************************** | |
| 670 * Complex inverse sine | |
| 671 * | |
| 672 * asin(z) = -i log( sqrt(1-$(POWER z, 2)) + iz) | |
| 673 * where both log and sqrt are complex. | |
| 674 */ | |
| 675 creal asin(creal z) | |
| 676 { | |
| 677 return -log(sqrt(1-z*z) + z*1i)*1i; | |
| 678 } | |
| 679 | |
| 680 debug(UnitTest) { | |
| 681 unittest { | |
| 682 assert(asin(sin(0+0i)) == 0 + 0i); | |
| 683 } | |
| 684 } | |
| 685 | |
| 686 /*********************************** | |
| 687 * Complex inverse cosine | |
| 688 * | |
| 689 * acos(z) = &pi/2 - asin(z) | |
| 690 */ | |
| 691 creal acos(creal z) | |
| 692 { | |
| 693 return PI_2 - asin(z); | |
| 694 } | |
| 695 | |
| 696 | |
| 697 /** | |
| 698 * Calculates the hyperbolic cosine of x. | |
| 699 * | |
| 700 * $(TABLE_SV | |
| 701 * <tr> <th> x <th> cosh(x) <th> invalid? | |
| 702 * <tr> <td> ±∞ <td> ±0.0 <td> no | |
| 703 * ) | |
| 704 */ | |
| 705 real cosh(real x) | |
| 706 { | |
| 707 return tango.stdc.math.coshl(x); | |
| 708 } | |
| 709 | |
| 710 debug(UnitTest) { | |
| 711 unittest { | |
| 712 // NaN payloads | |
| 713 assert(isIdentical(cosh(NaN(432)), NaN(432))); | |
| 714 } | |
| 715 } | |
| 716 | |
| 717 /** | |
| 718 * Calculates the hyperbolic sine of x. | |
| 719 * | |
| 720 * $(TABLE_SV | |
| 721 * <tr> <th> x <th> sinh(x) <th> invalid? | |
| 722 * <tr> <td> ±0.0 <td> ±0.0 <td> no | |
| 723 * <tr> <td> ±∞ <td> ±∞ <td> no | |
| 724 * ) | |
| 725 */ | |
| 726 real sinh(real x) | |
| 727 { | |
| 728 return tango.stdc.math.sinhl(x); | |
| 729 } | |
| 730 | |
| 731 debug(UnitTest) { | |
| 732 unittest { | |
| 733 // NaN payloads | |
| 734 assert(isIdentical(sinh(NaN(0xABC)), NaN(0xABC))); | |
| 735 } | |
| 736 } | |
| 737 | |
| 738 /** | |
| 739 * Calculates the hyperbolic tangent of x. | |
| 740 * | |
| 741 * $(TABLE_SV | |
| 742 * <tr> <th> x <th> tanh(x) <th> invalid? | |
| 743 * <tr> <td> ±0.0 <td> ±0.0 <td> no | |
| 744 * <tr> <td> ±∞ <td> ±1.0 <td> no | |
| 745 * ) | |
| 746 */ | |
| 747 real tanh(real x) | |
| 748 { | |
| 749 return tango.stdc.math.tanhl(x); | |
| 750 } | |
| 751 | |
| 752 debug(UnitTest) { | |
| 753 unittest { | |
| 754 // NaN payloads | |
| 755 assert(isIdentical(tanh(NaN(0xABC)), NaN(0xABC))); | |
| 756 } | |
| 757 } | |
| 758 | |
| 759 /*********************************** | |
| 760 * hyperbolic sine, complex and imaginary | |
| 761 * | |
| 762 * sinh(z) = cos(z.im)*sinh(z.re) + sin(z.im)*cosh(z.re)i | |
| 763 */ | |
| 764 creal sinh(creal z) | |
| 765 { | |
| 766 creal cs = expi(z.im); | |
| 767 return cs.re * sinh(z.re) + cs.im * cosh(z.re) * 1i; | |
| 768 } | |
| 769 | |
| 770 /** ditto */ | |
| 771 ireal sinh(ireal y) | |
| 772 { | |
| 773 return sin(y.im)*1i; | |
| 774 } | |
| 775 | |
| 776 debug(UnitTest) { | |
| 777 unittest { | |
| 778 assert(sinh(4.2L + 0i)==sinh(4.2L)); | |
| 779 } | |
| 780 } | |
| 781 | |
| 782 /*********************************** | |
| 783 * hyperbolic cosine, complex and imaginary | |
| 784 * | |
| 785 * cosh(z) = cos(z.im)*cosh(z.re) + sin(z.im)*sinh(z.re)i | |
| 786 */ | |
| 787 creal cosh(creal z) | |
| 788 { | |
| 789 creal cs = expi(z.im); | |
| 790 return cs.re * cosh(z.re) + cs.im * sinh(z.re) * 1i; | |
| 791 } | |
| 792 | |
| 793 /** ditto */ | |
| 794 real cosh(ireal y) | |
| 795 { | |
| 796 return cos(y.im); | |
| 797 } | |
| 798 | |
| 799 debug(UnitTest) { | |
| 800 unittest { | |
| 801 assert(cosh(8.3L + 0i)==cosh(8.3L)); | |
| 802 } | |
| 803 } | |
| 804 | |
| 805 | |
| 806 /** | |
| 807 * Calculates the inverse hyperbolic cosine of x. | |
| 808 * | |
| 809 * Mathematically, acosh(x) = log(x + sqrt( x*x - 1)) | |
| 810 * | |
| 811 * $(TABLE_DOMRG | |
| 812 * $(DOMAIN 1..∞) | |
| 813 * $(RANGE 1..log(real.max), ∞ ) | |
| 814 * ) | |
| 815 * $(TABLE_SV | |
| 816 * $(SVH x, acosh(x) ) | |
| 817 * $(SV $(NAN), $(NAN) ) | |
| 818 * $(SV <1, $(NAN) ) | |
| 819 * $(SV 1, 0 ) | |
| 820 * $(SV +∞,+∞) | |
| 821 * ) | |
| 822 */ | |
| 823 real acosh(real x) | |
| 824 { | |
| 825 if (x > 1/real.epsilon) | |
| 826 return LN2 + log(x); | |
| 827 else | |
| 828 return log(x + sqrt(x*x - 1)); | |
| 829 } | |
| 830 | |
| 831 debug(UnitTest) { | |
| 832 unittest | |
| 833 { | |
| 834 assert(isNaN(acosh(0.9))); | |
| 835 assert(isNaN(acosh(real.nan))); | |
| 836 assert(acosh(1)==0.0); | |
| 837 assert(acosh(real.infinity) == real.infinity); | |
| 838 // NaN payloads | |
| 839 assert(isIdentical(acosh(NaN(0xABC)), NaN(0xABC))); | |
| 840 } | |
| 841 } | |
| 842 | |
| 843 /** | |
| 844 * Calculates the inverse hyperbolic sine of x. | |
| 845 * | |
| 846 * Mathematically, | |
| 847 * --------------- | |
| 848 * asinh(x) = log( x + sqrt( x*x + 1 )) // if x >= +0 | |
| 849 * asinh(x) = -log(-x + sqrt( x*x + 1 )) // if x <= -0 | |
| 850 * ------------- | |
| 851 * | |
| 852 * $(TABLE_SV | |
| 853 * $(SVH x, asinh(x) ) | |
| 854 * $(SV $(NAN), $(NAN) ) | |
| 855 * $(SV ±0, ±0 ) | |
| 856 * $(SV ±∞,±∞) | |
| 857 * ) | |
| 858 */ | |
| 859 real asinh(real x) | |
| 860 { | |
| 861 if (tango.math.IEEE.fabs(x) > 1 / real.epsilon) // beyond this point, x*x + 1 == x*x | |
| 862 return tango.math.IEEE.copysign(LN2 + log(tango.math.IEEE.fabs(x)), x); | |
| 863 else | |
| 864 { | |
| 865 // sqrt(x*x + 1) == 1 + x * x / ( 1 + sqrt(x*x + 1) ) | |
| 866 return tango.math.IEEE.copysign(log1p(tango.math.IEEE.fabs(x) + x*x / (1 + sqrt(x*x + 1)) ), x); | |
| 867 } | |
| 868 } | |
| 869 | |
| 870 debug(UnitTest) { | |
| 871 unittest | |
| 872 { | |
| 873 assert(isIdentical(0.0L,asinh(0.0))); | |
| 874 assert(isIdentical(-0.0L,asinh(-0.0))); | |
| 875 assert(asinh(real.infinity) == real.infinity); | |
| 876 assert(asinh(-real.infinity) == -real.infinity); | |
| 877 assert(isNaN(asinh(real.nan))); | |
| 878 // NaN payloads | |
| 879 assert(isIdentical(asinh(NaN(0xABC)), NaN(0xABC))); | |
| 880 } | |
| 881 } | |
| 882 | |
| 883 /** | |
| 884 * Calculates the inverse hyperbolic tangent of x, | |
| 885 * returning a value from ranging from -1 to 1. | |
| 886 * | |
| 887 * Mathematically, atanh(x) = log( (1+x)/(1-x) ) / 2 | |
| 888 * | |
| 889 * | |
| 890 * $(TABLE_DOMRG | |
| 891 * $(DOMAIN -∞..∞) | |
| 892 * $(RANGE -1..1) ) | |
| 893 * $(TABLE_SV | |
| 894 * $(SVH x, atanh(x) ) | |
| 895 * $(SV $(NAN), $(NAN) ) | |
| 896 * $(SV ±0, ±0) | |
| 897 * $(SV ±1, ±∞) | |
| 898 * ) | |
| 899 */ | |
| 900 real atanh(real x) | |
| 901 { | |
| 902 // log( (1+x)/(1-x) ) == log ( 1 + (2*x)/(1-x) ) | |
| 903 return 0.5 * log1p( 2 * x / (1 - x) ); | |
| 904 } | |
| 905 | |
| 906 debug(UnitTest) { | |
| 907 unittest | |
| 908 { | |
| 909 assert(isIdentical(0.0L, atanh(0.0))); | |
| 910 assert(isIdentical(-0.0L,atanh(-0.0))); | |
| 911 assert(isIdentical(atanh(-1),-real.infinity)); | |
| 912 assert(isIdentical(atanh(1),real.infinity)); | |
| 913 assert(isNaN(atanh(-real.infinity))); | |
| 914 // NaN payloads | |
| 915 assert(isIdentical(atanh(NaN(0xABC)), NaN(0xABC))); | |
| 916 } | |
| 917 } | |
| 918 | |
| 919 /** ditto */ | |
| 920 creal atanh(ireal y) | |
| 921 { | |
| 922 // Not optimised for accuracy or speed | |
| 923 return 0.5*(log(1+y) - log(1-y)); | |
| 924 } | |
| 925 | |
| 926 /** ditto */ | |
| 927 creal atanh(creal z) | |
| 928 { | |
| 929 // Not optimised for accuracy or speed | |
| 930 return 0.5 * (log(1 + z) - log(1-z)); | |
| 931 } | |
| 932 | |
| 933 /* | |
| 934 * Powers and Roots | |
| 935 */ | |
| 936 | |
| 937 /** | |
| 938 * Compute square root of x. | |
| 939 * | |
| 940 * $(TABLE_SV | |
| 941 * <tr> <th> x <th> sqrt(x) <th> invalid? | |
| 942 * <tr> <td> -0.0 <td> -0.0 <td> no | |
| 943 * <tr> <td> <0.0 <td> $(NAN) <td> yes | |
| 944 * <tr> <td> +∞ <td> +∞ <td> no | |
| 945 * ) | |
| 946 */ | |
| 947 float sqrt(float x) /* intrinsic */ | |
| 948 { | |
| 949 version(D_InlineAsm_X86) | |
| 950 { | |
| 951 asm | |
| 952 { | |
| 953 fld x; | |
| 954 fsqrt; | |
| 955 } | |
| 956 } | |
| 957 else | |
| 958 { | |
| 959 return tango.stdc.math.sqrtf(x); | |
| 960 } | |
| 961 } | |
| 962 | |
| 963 double sqrt(double x) /* intrinsic */ /// ditto | |
| 964 { | |
| 965 version(D_InlineAsm_X86) | |
| 966 { | |
| 967 asm | |
| 968 { | |
| 969 fld x; | |
| 970 fsqrt; | |
| 971 } | |
| 972 } | |
| 973 else | |
| 974 { | |
| 975 return tango.stdc.math.sqrt(x); | |
| 976 } | |
| 977 } | |
| 978 | |
| 979 real sqrt(real x) /* intrinsic */ /// ditto | |
| 980 { | |
| 981 version(D_InlineAsm_X86) | |
| 982 { | |
| 983 asm | |
| 984 { | |
| 985 fld x; | |
| 986 fsqrt; | |
| 987 } | |
| 988 } | |
| 989 else | |
| 990 { | |
| 991 return tango.stdc.math.sqrtl(x); | |
| 992 } | |
| 993 } | |
| 994 | |
| 995 /** ditto */ | |
| 996 creal sqrt(creal z) | |
| 997 { | |
| 998 | |
| 999 if (z == 0.0) return z; | |
| 1000 real x,y,w,r; | |
| 1001 creal c; | |
| 1002 | |
| 1003 x = tango.math.IEEE.fabs(z.re); | |
| 1004 y = tango.math.IEEE.fabs(z.im); | |
| 1005 if (x >= y) { | |
| 1006 r = y / x; | |
| 1007 w = sqrt(x) * sqrt(0.5 * (1 + sqrt(1 + r * r))); | |
| 1008 } else { | |
| 1009 r = x / y; | |
| 1010 w = sqrt(y) * sqrt(0.5 * (r + sqrt(1 + r * r))); | |
| 1011 } | |
| 1012 | |
| 1013 if (z.re >= 0) { | |
| 1014 c = w + (z.im / (w + w)) * 1.0i; | |
| 1015 } else { | |
| 1016 if (z.im < 0) w = -w; | |
| 1017 c = z.im / (w + w) + w * 1.0i; | |
| 1018 } | |
| 1019 return c; | |
| 1020 } | |
| 1021 | |
| 1022 import tango.stdc.stdio; | |
| 1023 debug(UnitTest) { | |
| 1024 unittest { | |
| 1025 // NaN payloads | |
| 1026 assert(isIdentical(sqrt(NaN(0xABC)), NaN(0xABC))); | |
| 1027 assert(sqrt(-1+0i) == 1i); | |
| 1028 assert(isIdentical(sqrt(0-0i), 0-0i)); | |
| 1029 assert(cfeqrel(sqrt(4+16i)*sqrt(4+16i), 4+16i)>=real.mant_dig-2); | |
| 1030 } | |
| 1031 } | |
| 1032 | |
| 1033 /** | |
| 1034 * Calculates the cube root of x. | |
| 1035 * | |
| 1036 * $(TABLE_SV | |
| 1037 * <tr> <th> <i>x</i> <th> cbrt(x) <th> invalid? | |
| 1038 * <tr> <td> ±0.0 <td> ±0.0 <td> no | |
| 1039 * <tr> <td> $(NAN) <td> $(NAN) <td> yes | |
| 1040 * <tr> <td> ±∞ <td> ±∞ <td> no | |
| 1041 * ) | |
| 1042 */ | |
| 1043 real cbrt(real x) | |
| 1044 { | |
| 1045 return tango.stdc.math.cbrtl(x); | |
| 1046 } | |
| 1047 | |
| 1048 | |
| 1049 debug(UnitTest) { | |
| 1050 unittest { | |
| 1051 // NaN payloads | |
| 1052 assert(isIdentical(cbrt(NaN(0xABC)), NaN(0xABC))); | |
| 1053 } | |
| 1054 } | |
| 1055 | |
| 1056 /** | |
| 1057 * Calculates e$(SUP x). | |
| 1058 * | |
| 1059 * $(TABLE_SV | |
| 1060 * <tr> <th> x <th> exp(x) | |
| 1061 * <tr> <td> +∞ <td> +∞ | |
| 1062 * <tr> <td> -∞ <td> +0.0 | |
| 1063 * ) | |
| 1064 */ | |
| 1065 real exp(real x) | |
| 1066 { | |
| 1067 return tango.stdc.math.expl(x); | |
| 1068 } | |
| 1069 | |
| 1070 debug(UnitTest) { | |
| 1071 unittest { | |
| 1072 // NaN payloads | |
| 1073 assert(isIdentical(exp(NaN(0xABC)), NaN(0xABC))); | |
| 1074 } | |
| 1075 } | |
| 1076 | |
| 1077 /** | |
| 1078 * Calculates the value of the natural logarithm base (e) | |
| 1079 * raised to the power of x, minus 1. | |
| 1080 * | |
| 1081 * For very small x, expm1(x) is more accurate | |
| 1082 * than exp(x)-1. | |
| 1083 * | |
| 1084 * $(TABLE_SV | |
| 1085 * <tr> <th> x <th> e$(SUP x)-1 | |
| 1086 * <tr> <td> ±0.0 <td> ±0.0 | |
| 1087 * <tr> <td> +∞ <td> +∞ | |
| 1088 * <tr> <td> -∞ <td> -1.0 | |
| 1089 * ) | |
| 1090 */ | |
| 1091 real expm1(real x) | |
| 1092 { | |
| 1093 return tango.stdc.math.expm1l(x); | |
| 1094 } | |
| 1095 | |
| 1096 debug(UnitTest) { | |
| 1097 unittest { | |
| 1098 // NaN payloads | |
| 1099 assert(isIdentical(expm1(NaN(0xABC)), NaN(0xABC))); | |
| 1100 } | |
| 1101 } | |
| 1102 | |
| 1103 | |
| 1104 /** | |
| 1105 * Calculates 2$(SUP x). | |
| 1106 * | |
| 1107 * $(TABLE_SV | |
| 1108 * <tr> <th> x <th> exp2(x) | |
| 1109 * <tr> <td> +∞ <td> +∞ | |
| 1110 * <tr> <td> -∞ <td> +0.0 | |
| 1111 * ) | |
| 1112 */ | |
| 1113 real exp2(real x) | |
| 1114 { | |
| 1115 return tango.stdc.math.exp2l(x); | |
| 1116 } | |
| 1117 | |
| 1118 debug(UnitTest) { | |
| 1119 unittest { | |
| 1120 // NaN payloads | |
| 1121 assert(isIdentical(exp2(NaN(0xABC)), NaN(0xABC))); | |
| 1122 } | |
| 1123 } | |
| 1124 | |
| 1125 /* | |
| 1126 * Powers and Roots | |
| 1127 */ | |
| 1128 | |
| 1129 /** | |
| 1130 * Calculate the natural logarithm of x. | |
| 1131 * | |
| 1132 * $(TABLE_SV | |
| 1133 * <tr> <th> x <th> log(x) <th> divide by 0? <th> invalid? | |
| 1134 * <tr> <td> ±0.0 <td> -∞ <td> yes <td> no | |
| 1135 * <tr> <td> < 0.0 <td> $(NAN) <td> no <td> yes | |
| 1136 * <tr> <td> +∞ <td> +∞ <td> no <td> no | |
| 1137 * ) | |
| 1138 */ | |
| 1139 real log(real x) | |
| 1140 { | |
| 1141 return tango.stdc.math.logl(x); | |
| 1142 } | |
| 1143 | |
| 1144 debug(UnitTest) { | |
| 1145 unittest { | |
| 1146 // NaN payloads | |
| 1147 assert(isIdentical(log(NaN(0xABC)), NaN(0xABC))); | |
| 1148 } | |
| 1149 } | |
| 1150 | |
| 1151 /** | |
| 1152 * Calculates the natural logarithm of 1 + x. | |
| 1153 * | |
| 1154 * For very small x, log1p(x) will be more accurate than | |
| 1155 * log(1 + x). | |
| 1156 * | |
| 1157 * $(TABLE_SV | |
| 1158 * <tr> <th> x <th> log1p(x) <th> divide by 0? <th> invalid? | |
| 1159 * <tr> <td> ±0.0 <td> ±0.0 <td> no <td> no | |
| 1160 * <tr> <td> -1.0 <td> -∞ <td> yes <td> no | |
| 1161 * <tr> <td> <-1.0 <td> $(NAN) <td> no <td> yes | |
| 1162 * <tr> <td> +∞ <td> -∞ <td> no <td> no | |
| 1163 * ) | |
| 1164 */ | |
| 1165 real log1p(real x) | |
| 1166 { | |
| 1167 return tango.stdc.math.log1pl(x); | |
| 1168 } | |
| 1169 | |
| 1170 debug(UnitTest) { | |
| 1171 unittest { | |
| 1172 // NaN payloads | |
| 1173 assert(isIdentical(log1p(NaN(0xABC)), NaN(0xABC))); | |
| 1174 } | |
| 1175 } | |
| 1176 | |
| 1177 /** | |
| 1178 * Calculates the base-2 logarithm of x: | |
| 1179 * log<sub>2</sub>x | |
| 1180 * | |
| 1181 * $(TABLE_SV | |
| 1182 * <tr> <th> x <th> log2(x) <th> divide by 0? <th> invalid? | |
| 1183 * <tr> <td> ±0.0 <td> -∞ <td> yes <td> no | |
| 1184 * <tr> <td> < 0.0 <td> $(NAN) <td> no <td> yes | |
| 1185 * <tr> <td> +∞ <td> +∞ <td> no <td> no | |
| 1186 * ) | |
| 1187 */ | |
| 1188 real log2(real x) | |
| 1189 { | |
| 1190 return tango.stdc.math.log2l(x); | |
| 1191 } | |
| 1192 | |
| 1193 debug(UnitTest) { | |
| 1194 unittest { | |
| 1195 // NaN payloads | |
| 1196 assert(isIdentical(log2(NaN(0xABC)), NaN(0xABC))); | |
| 1197 } | |
| 1198 } | |
| 1199 | |
| 1200 /** | |
| 1201 * Calculate the base-10 logarithm of x. | |
| 1202 * | |
| 1203 * $(TABLE_SV | |
| 1204 * <tr> <th> x <th> log10(x) <th> divide by 0? <th> invalid? | |
| 1205 * <tr> <td> ±0.0 <td> -∞ <td> yes <td> no | |
| 1206 * <tr> <td> < 0.0 <td> $(NAN) <td> no <td> yes | |
| 1207 * <tr> <td> +∞ <td> +∞ <td> no <td> no | |
| 1208 * ) | |
| 1209 */ | |
| 1210 real log10(real x) | |
| 1211 { | |
| 1212 return tango.stdc.math.log10l(x); | |
| 1213 } | |
| 1214 | |
| 1215 debug(UnitTest) { | |
| 1216 unittest { | |
| 1217 // NaN payloads | |
| 1218 assert(isIdentical(log10(NaN(0xABC)), NaN(0xABC))); | |
| 1219 } | |
| 1220 } | |
| 1221 | |
| 1222 /*********************************** | |
| 1223 * Exponential, complex and imaginary | |
| 1224 * | |
| 1225 * For complex numbers, the exponential function is defined as | |
| 1226 * | |
| 1227 * exp(z) = exp(z.re)cos(z.im) + exp(z.re)sin(z.im)i. | |
| 1228 * | |
| 1229 * For a pure imaginary argument, | |
| 1230 * exp(θi) = cos(θ) + sin(θ)i. | |
| 1231 * | |
| 1232 */ | |
| 1233 creal exp(ireal y) | |
| 1234 { | |
| 1235 return expi(y.im); | |
| 1236 } | |
| 1237 | |
| 1238 /** ditto */ | |
| 1239 creal exp(creal z) | |
| 1240 { | |
| 1241 return expi(z.im) * exp(z.re); | |
| 1242 } | |
| 1243 | |
| 1244 debug(UnitTest) { | |
| 1245 unittest { | |
| 1246 assert(exp(1.3e5Li)==cos(1.3e5L)+sin(1.3e5L)*1i); | |
| 1247 assert(exp(0.0Li)==1L+0.0Li); | |
| 1248 assert(exp(7.2 + 0.0i) == exp(7.2L)); | |
| 1249 creal c = exp(ireal.nan); | |
| 1250 assert(isNaN(c.re) && isNaN(c.im)); | |
| 1251 c = exp(ireal.infinity); | |
| 1252 assert(isNaN(c.re) && isNaN(c.im)); | |
| 1253 } | |
| 1254 } | |
| 1255 | |
| 1256 /*********************************** | |
| 1257 * Natural logarithm, complex | |
| 1258 * | |
| 1259 * Returns complex logarithm to the base e (2.718...) of | |
| 1260 * the complex argument x. | |
| 1261 * | |
| 1262 * If z = x + iy, then | |
| 1263 * log(z) = log(abs(z)) + i arctan(y/x). | |
| 1264 * | |
| 1265 * The arctangent ranges from -PI to +PI. | |
| 1266 * There are branch cuts along both the negative real and negative | |
| 1267 * imaginary axes. For pure imaginary arguments, use one of the | |
| 1268 * following forms, depending on which branch is required. | |
| 1269 * ------------ | |
| 1270 * log( 0.0 + yi) = log(-y) + PI_2i // y<=-0.0 | |
| 1271 * log(-0.0 + yi) = log(-y) - PI_2i // y<=-0.0 | |
| 1272 * ------------ | |
| 1273 */ | |
| 1274 creal log(creal z) | |
| 1275 { | |
| 1276 return log(abs(z)) + atan2(z.im, z.re)*1i; | |
| 1277 } | |
| 1278 | |
| 1279 debug(UnitTest) { | |
| 1280 private { | |
| 1281 /* | |
| 1282 * feqrel for complex numbers. Returns the worst relative | |
| 1283 * equality of the two components. | |
| 1284 */ | |
| 1285 int cfeqrel(creal a, creal b) | |
| 1286 { | |
| 1287 int intmin(int a, int b) { return a<b? a: b; } | |
| 1288 return intmin(feqrel(a.re, b.re), feqrel(a.im, b.im)); | |
| 1289 } | |
| 1290 } | |
| 1291 unittest { | |
| 1292 | |
| 1293 assert(log(3.0L +0i) == log(3.0L)+0i); | |
| 1294 assert(cfeqrel(log(0.0L-2i),( log(2.0L)-PI_2*1i)) >= real.mant_dig-10); | |
| 1295 assert(cfeqrel(log(0.0L+2i),( log(2.0L)+PI_2*1i)) >= real.mant_dig-10); | |
| 1296 } | |
| 1297 } | |
| 1298 | |
| 1299 /** | |
| 1300 * Fast integral powers. | |
| 1301 */ | |
| 1302 real pow(real x, uint n) | |
| 1303 { | |
| 1304 real p; | |
| 1305 | |
| 1306 switch (n) | |
| 1307 { | |
| 1308 case 0: | |
| 1309 p = 1.0; | |
| 1310 break; | |
| 1311 | |
| 1312 case 1: | |
| 1313 p = x; | |
| 1314 break; | |
| 1315 | |
| 1316 case 2: | |
| 1317 p = x * x; | |
| 1318 break; | |
| 1319 | |
| 1320 default: | |
| 1321 p = 1.0; | |
| 1322 while (1){ | |
| 1323 if (n & 1) | |
| 1324 p *= x; | |
| 1325 n >>= 1; | |
| 1326 if (!n) | |
| 1327 break; | |
| 1328 x *= x; | |
| 1329 } | |
| 1330 break; | |
| 1331 } | |
| 1332 return p; | |
| 1333 } | |
| 1334 | |
| 1335 /** ditto */ | |
| 1336 real pow(real x, int n) | |
| 1337 { | |
| 1338 if (n < 0) return pow(x, cast(real)n); | |
| 1339 else return pow(x, cast(uint)n); | |
| 1340 } | |
| 1341 | |
| 1342 /** | |
| 1343 * Calculates x$(SUP y). | |
| 1344 * | |
| 1345 * $(TABLE_SV | |
| 1346 * <tr> | |
| 1347 * <th> x <th> y <th> pow(x, y) <th> div 0 <th> invalid? | |
| 1348 * <tr> | |
| 1349 * <td> anything <td> ±0.0 <td> 1.0 <td> no <td> no | |
| 1350 * <tr> | |
| 1351 * <td> |x| > 1 <td> +∞ <td> +∞ <td> no <td> no | |
| 1352 * <tr> | |
| 1353 * <td> |x| < 1 <td> +∞ <td> +0.0 <td> no <td> no | |
| 1354 * <tr> | |
| 1355 * <td> |x| > 1 <td> -∞ <td> +0.0 <td> no <td> no | |
| 1356 * <tr> | |
| 1357 * <td> |x| < 1 <td> -∞ <td> +∞ <td> no <td> no | |
| 1358 * <tr> | |
| 1359 * <td> +∞ <td> > 0.0 <td> +∞ <td> no <td> no | |
| 1360 * <tr> | |
| 1361 * <td> +∞ <td> < 0.0 <td> +0.0 <td> no <td> no | |
| 1362 * <tr> | |
| 1363 * <td> -∞ <td> odd integer > 0.0 <td> -∞ <td> no <td> no | |
| 1364 * <tr> | |
| 1365 * <td> -∞ <td> > 0.0, not odd integer <td> +∞ <td> no <td> no | |
| 1366 * <tr> | |
| 1367 * <td> -∞ <td> odd integer < 0.0 <td> -0.0 <td> no <td> no | |
| 1368 * <tr> | |
| 1369 * <td> -∞ <td> < 0.0, not odd integer <td> +0.0 <td> no <td> no | |
| 1370 * <tr> | |
| 1371 * <td> ±1.0 <td> ±∞ <td> $(NAN) <td> no <td> yes | |
| 1372 * <tr> | |
| 1373 * <td> < 0.0 <td> finite, nonintegral <td> $(NAN) <td> no <td> yes | |
| 1374 * <tr> | |
| 1375 * <td> ±0.0 <td> odd integer < 0.0 <td> ±∞ <td> yes <td> no | |
| 1376 * <tr> | |
| 1377 * <td> ±0.0 <td> < 0.0, not odd integer <td> +∞ <td> yes <td> no | |
| 1378 * <tr> | |
| 1379 * <td> ±0.0 <td> odd integer > 0.0 <td> ±0.0 <td> no <td> no | |
| 1380 * <tr> | |
| 1381 * <td> ±0.0 <td> > 0.0, not odd integer <td> +0.0 <td> no <td> no | |
| 1382 * ) | |
| 1383 */ | |
| 1384 real pow(real x, real y) | |
| 1385 { | |
| 1386 version (linux) // C pow() often does not handle special values correctly | |
| 1387 { | |
| 1388 if (isNaN(y)) | |
| 1389 return y; | |
| 1390 | |
| 1391 if (y == 0) | |
| 1392 return 1; // even if x is $(NAN) | |
| 1393 if (isNaN(x) && y != 0) | |
| 1394 return x; | |
| 1395 if (isInfinity(y)) | |
| 1396 { | |
| 1397 if (tango.math.IEEE.fabs(x) > 1) | |
| 1398 { | |
| 1399 if (signbit(y)) | |
| 1400 return +0.0; | |
| 1401 else | |
| 1402 return real.infinity; | |
| 1403 } | |
| 1404 else if (tango.math.IEEE.fabs(x) == 1) | |
| 1405 { | |
| 1406 return NaN(TANGO_NAN.POW_DOMAIN); | |
| 1407 } | |
| 1408 else // < 1 | |
| 1409 { | |
| 1410 if (signbit(y)) | |
| 1411 return real.infinity; | |
| 1412 else | |
| 1413 return +0.0; | |
| 1414 } | |
| 1415 } | |
| 1416 if (isInfinity(x)) | |
| 1417 { | |
| 1418 if (signbit(x)) | |
| 1419 { | |
| 1420 long i; | |
| 1421 i = cast(long)y; | |
| 1422 if (y > 0) | |
| 1423 { | |
| 1424 if (i == y && i & 1) | |
| 1425 return -real.infinity; | |
| 1426 else | |
| 1427 return real.infinity; | |
| 1428 } | |
| 1429 else if (y < 0) | |
| 1430 { | |
| 1431 if (i == y && i & 1) | |
| 1432 return -0.0; | |
| 1433 else | |
| 1434 return +0.0; | |
| 1435 } | |
| 1436 } | |
| 1437 else | |
| 1438 { | |
| 1439 if (y > 0) | |
| 1440 return real.infinity; | |
| 1441 else if (y < 0) | |
| 1442 return +0.0; | |
| 1443 } | |
| 1444 } | |
| 1445 | |
| 1446 if (x == 0.0) | |
| 1447 { | |
| 1448 if (signbit(x)) | |
| 1449 { | |
| 1450 long i; | |
| 1451 | |
| 1452 i = cast(long)y; | |
| 1453 if (y > 0) | |
| 1454 { | |
| 1455 if (i == y && i & 1) | |
| 1456 return -0.0; | |
| 1457 else | |
| 1458 return +0.0; | |
| 1459 } | |
| 1460 else if (y < 0) | |
| 1461 { | |
| 1462 if (i == y && i & 1) | |
| 1463 return -real.infinity; | |
| 1464 else | |
| 1465 return real.infinity; | |
| 1466 } | |
| 1467 } | |
| 1468 else | |
| 1469 { | |
| 1470 if (y > 0) | |
| 1471 return +0.0; | |
| 1472 else if (y < 0) | |
| 1473 return real.infinity; | |
| 1474 } | |
| 1475 } | |
| 1476 } | |
| 1477 return tango.stdc.math.powl(x, y); | |
| 1478 } | |
| 1479 | |
| 1480 debug(UnitTest) { | |
| 1481 unittest | |
| 1482 { | |
| 1483 real x = 46; | |
| 1484 | |
| 1485 assert(pow(x,0) == 1.0); | |
| 1486 assert(pow(x,1) == x); | |
| 1487 assert(pow(x,2) == x * x); | |
| 1488 assert(pow(x,3) == x * x * x); | |
| 1489 assert(pow(x,8) == (x * x) * (x * x) * (x * x) * (x * x)); | |
| 1490 // NaN payloads | |
| 1491 assert(isIdentical(pow(NaN(0xABC), 19), NaN(0xABC))); | |
| 1492 } | |
| 1493 } | |
| 1494 | |
| 1495 /** | |
| 1496 * Calculates the length of the | |
| 1497 * hypotenuse of a right-angled triangle with sides of length x and y. | |
| 1498 * The hypotenuse is the value of the square root of | |
| 1499 * the sums of the squares of x and y: | |
| 1500 * | |
| 1501 * sqrt(x² + y²) | |
| 1502 * | |
| 1503 * Note that hypot(x, y), hypot(y, x) and | |
| 1504 * hypot(x, -y) are equivalent. | |
| 1505 * | |
| 1506 * $(TABLE_SV | |
| 1507 * <tr> <th> x <th> y <th> hypot(x, y) <th> invalid? | |
| 1508 * <tr> <td> x <td> ±0.0 <td> |x| <td> no | |
| 1509 * <tr> <td> ±∞ <td> y <td> +∞ <td> no | |
| 1510 * <tr> <td> ±∞ <td> $(NAN) <td> +∞ <td> no | |
| 1511 * ) | |
| 1512 */ | |
| 1513 real hypot(real x, real y) | |
| 1514 { | |
| 1515 /* | |
| 1516 * This is based on code from: | |
| 1517 * Cephes Math Library Release 2.1: January, 1989 | |
| 1518 * Copyright 1984, 1987, 1989 by Stephen L. Moshier | |
| 1519 * Direct inquiries to 30 Frost Street, Cambridge, MA 02140 | |
| 1520 */ | |
| 1521 | |
| 1522 const int PRECL = real.mant_dig/2; // = 32 | |
| 1523 | |
| 1524 real xx, yy, b, re, im; | |
| 1525 int ex, ey, e; | |
| 1526 | |
| 1527 // Note, hypot(INFINITY, NAN) = INFINITY. | |
| 1528 if (tango.math.IEEE.isInfinity(x) || tango.math.IEEE.isInfinity(y)) | |
| 1529 return real.infinity; | |
| 1530 | |
| 1531 if (tango.math.IEEE.isNaN(x)) | |
| 1532 return x; | |
| 1533 if (tango.math.IEEE.isNaN(y)) | |
| 1534 return y; | |
| 1535 | |
| 1536 re = tango.math.IEEE.fabs(x); | |
| 1537 im = tango.math.IEEE.fabs(y); | |
| 1538 | |
| 1539 if (re == 0.0) | |
| 1540 return im; | |
| 1541 if (im == 0.0) | |
| 1542 return re; | |
| 1543 | |
| 1544 // Get the exponents of the numbers | |
| 1545 xx = tango.math.IEEE.frexp(re, ex); | |
| 1546 yy = tango.math.IEEE.frexp(im, ey); | |
| 1547 | |
| 1548 // Check if one number is tiny compared to the other | |
| 1549 e = ex - ey; | |
| 1550 if (e > PRECL) | |
| 1551 return re; | |
| 1552 if (e < -PRECL) | |
| 1553 return im; | |
| 1554 | |
| 1555 // Find approximate exponent e of the geometric mean. | |
| 1556 e = (ex + ey) >> 1; | |
| 1557 | |
| 1558 // Rescale so mean is about 1 | |
| 1559 xx = tango.math.IEEE.ldexp(re, -e); | |
| 1560 yy = tango.math.IEEE.ldexp(im, -e); | |
| 1561 | |
| 1562 // Hypotenuse of the right triangle | |
| 1563 b = sqrt(xx * xx + yy * yy); | |
| 1564 | |
| 1565 // Compute the exponent of the answer. | |
| 1566 yy = tango.math.IEEE.frexp(b, ey); | |
| 1567 ey = e + ey; | |
| 1568 | |
| 1569 // Check it for overflow and underflow. | |
| 1570 if (ey > real.max_exp + 2) { | |
| 1571 return real.infinity; | |
| 1572 } | |
| 1573 if (ey < real.min_exp - 2) | |
| 1574 return 0.0; | |
| 1575 | |
| 1576 // Undo the scaling | |
| 1577 b = tango.math.IEEE.ldexp(b, e); | |
| 1578 return b; | |
| 1579 } | |
| 1580 | |
| 1581 debug(UnitTest) { | |
| 1582 unittest | |
| 1583 { | |
| 1584 static real vals[][3] = // x,y,hypot | |
| 1585 [ | |
| 1586 [ 0, 0, 0], | |
| 1587 [ 0, -0, 0], | |
| 1588 [ 3, 4, 5], | |
| 1589 [ -300, -400, 500], | |
| 1590 [ real.min, real.min, 0x1.6a09e667f3bcc908p-16382L], | |
| 1591 [ real.max/2, real.max/2, 0x1.6a09e667f3bcc908p+16383L /*8.41267e+4931L*/], | |
| 1592 [ real.max, 1, real.max], | |
| 1593 [ real.infinity, real.nan, real.infinity], | |
| 1594 [ real.nan, real.nan, real.nan], | |
| 1595 ]; | |
| 1596 | |
| 1597 for (int i = 0; i < vals.length; i++) | |
| 1598 { | |
| 1599 real x = vals[i][0]; | |
| 1600 real y = vals[i][1]; | |
| 1601 real z = vals[i][2]; | |
| 1602 real h = hypot(x, y); | |
| 1603 | |
| 1604 // printf("hypot(%La, %La) = %La, should be %La\n", x, y, h, z); | |
| 1605 assert(isIdentical(z, h)); | |
| 1606 } | |
| 1607 // NaN payloads | |
| 1608 assert(isIdentical(hypot(NaN(0xABC), 3.14), NaN(0xABC))); | |
| 1609 assert(isIdentical(hypot(7.6e39, NaN(0xABC)), NaN(0xABC))); | |
| 1610 } | |
| 1611 } | |
| 1612 | |
| 1613 /** | |
| 1614 * Evaluate polynomial A(x) = a<sub>0</sub> + a<sub>1</sub>x + a<sub>2</sub>x² + a<sub>3</sub>x³ ... | |
| 1615 * | |
| 1616 * Uses Horner's rule A(x) = a<sub>0</sub> + x(a<sub>1</sub> + x(a<sub>2</sub> + x(a<sub>3</sub> + ...))) | |
| 1617 * Params: | |
| 1618 * A = array of coefficients a<sub>0</sub>, a<sub>1</sub>, etc. | |
| 1619 */ | |
| 1620 T poly(T)(T x, T[] A) | |
| 1621 in | |
| 1622 { | |
| 1623 assert(A.length > 0); | |
| 1624 } | |
| 1625 body | |
| 1626 { | |
| 1627 version (DigitalMars_D_InlineAsm_X86) { | |
| 1628 const bool Use_D_InlineAsm_X86 = true; | |
| 1629 } else const bool Use_D_InlineAsm_X86 = false; | |
| 1630 | |
| 1631 // BUG (Inherited from Phobos): This code assumes a frame pointer in EBP. | |
| 1632 // This is not in the spec. | |
| 1633 static if (Use_D_InlineAsm_X86 && is(T==real) && T.sizeof == 10) { | |
| 1634 asm // assembler by W. Bright | |
| 1635 { | |
| 1636 // EDX = (A.length - 1) * real.sizeof | |
| 1637 mov ECX,A[EBP] ; // ECX = A.length | |
| 1638 dec ECX ; | |
| 1639 lea EDX,[ECX][ECX*8] ; | |
| 1640 add EDX,ECX ; | |
| 1641 add EDX,A+4[EBP] ; | |
| 1642 fld real ptr [EDX] ; // ST0 = coeff[ECX] | |
| 1643 jecxz return_ST ; | |
| 1644 fld x[EBP] ; // ST0 = x | |
| 1645 fxch ST(1) ; // ST1 = x, ST0 = r | |
| 1646 align 4 ; | |
| 1647 L2: fmul ST,ST(1) ; // r *= x | |
| 1648 fld real ptr -10[EDX] ; | |
| 1649 sub EDX,10 ; // deg-- | |
| 1650 faddp ST(1),ST ; | |
| 1651 dec ECX ; | |
| 1652 jne L2 ; | |
| 1653 fxch ST(1) ; // ST1 = r, ST0 = x | |
| 1654 fstp ST(0) ; // dump x | |
| 1655 align 4 ; | |
| 1656 return_ST: ; | |
| 1657 ; | |
| 1658 } | |
| 1659 } else static if ( Use_D_InlineAsm_X86 && is(T==real) && T.sizeof==12){ | |
| 1660 asm // assembler by W. Bright | |
| 1661 { | |
| 1662 // EDX = (A.length - 1) * real.sizeof | |
| 1663 mov ECX,A[EBP] ; // ECX = A.length | |
| 1664 dec ECX ; | |
| 1665 lea EDX,[ECX*8] ; | |
| 1666 lea EDX,[EDX][ECX*4] ; | |
| 1667 add EDX,A+4[EBP] ; | |
| 1668 fld real ptr [EDX] ; // ST0 = coeff[ECX] | |
| 1669 jecxz return_ST ; | |
| 1670 fld x ; // ST0 = x | |
| 1671 fxch ST(1) ; // ST1 = x, ST0 = r | |
| 1672 align 4 ; | |
| 1673 L2: fmul ST,ST(1) ; // r *= x | |
| 1674 fld real ptr -12[EDX] ; | |
| 1675 sub EDX,12 ; // deg-- | |
| 1676 faddp ST(1),ST ; | |
| 1677 dec ECX ; | |
| 1678 jne L2 ; | |
| 1679 fxch ST(1) ; // ST1 = r, ST0 = x | |
| 1680 fstp ST(0) ; // dump x | |
| 1681 align 4 ; | |
| 1682 return_ST: ; | |
| 1683 ; | |
| 1684 } | |
| 1685 } else { | |
| 1686 ptrdiff_t i = A.length - 1; | |
| 1687 real r = A[i]; | |
| 1688 while (--i >= 0) | |
| 1689 { | |
| 1690 r *= x; | |
| 1691 r += A[i]; | |
| 1692 } | |
| 1693 return r; | |
| 1694 } | |
| 1695 } | |
| 1696 | |
| 1697 debug(UnitTest) { | |
| 1698 unittest | |
| 1699 { | |
| 1700 debug (math) printf("math.poly.unittest\n"); | |
| 1701 real x = 3.1; | |
| 1702 const real pp[] = [56.1L, 32.7L, 6L]; | |
| 1703 | |
| 1704 assert( poly(x, pp) == (56.1L + (32.7L + 6L * x) * x) ); | |
| 1705 | |
| 1706 assert(isIdentical(poly(NaN(0xABC), pp), NaN(0xABC))); | |
| 1707 } | |
| 1708 } | |
| 1709 | |
| 1710 package { | |
| 1711 T rationalPoly(T)(T x, T [] numerator, T [] denominator) | |
| 1712 { | |
| 1713 return poly(x, numerator)/poly(x, denominator); | |
| 1714 } | |
| 1715 } | |
| 1716 | |
| 1717 private enum : int { MANTDIG_2 = real.mant_dig/2 } // Compiler workaround | |
| 1718 | |
| 1719 /** Floating point "approximate equality". | |
| 1720 * | |
| 1721 * Return true if x is equal to y, to within the specified precision | |
| 1722 * If roundoffbits is not specified, a reasonable default is used. | |
| 1723 */ | |
| 1724 bool feq(int precision = MANTDIG_2, XReal=real, YReal=real)(XReal x, YReal y) | |
| 1725 { | |
| 1726 static assert(is( XReal: real) && is(YReal : real)); | |
| 1727 return tango.math.IEEE.feqrel(x, y) >= precision; | |
| 1728 } | |
| 1729 | |
| 1730 unittest{ | |
| 1731 assert(!feq(1.0,2.0)); | |
| 1732 real y = 58.0000000001; | |
| 1733 assert(feq!(20)(58, y)); | |
| 1734 } | |
| 1735 | |
| 1736 /* | |
| 1737 * Rounding (returning real) | |
| 1738 */ | |
| 1739 | |
| 1740 /** | |
| 1741 * Returns the value of x rounded downward to the next integer | |
| 1742 * (toward negative infinity). | |
| 1743 */ | |
| 1744 real floor(real x) | |
| 1745 { | |
| 1746 return tango.stdc.math.floorl(x); | |
| 1747 } | |
| 1748 | |
| 1749 debug(UnitTest) { | |
| 1750 unittest { | |
| 1751 assert(isIdentical(floor(NaN(0xABC)), NaN(0xABC))); | |
| 1752 } | |
| 1753 } | |
| 1754 | |
| 1755 /** | |
| 1756 * Returns the value of x rounded upward to the next integer | |
| 1757 * (toward positive infinity). | |
| 1758 */ | |
| 1759 real ceil(real x) | |
| 1760 { | |
| 1761 return tango.stdc.math.ceill(x); | |
| 1762 } | |
| 1763 | |
| 1764 unittest { | |
| 1765 assert(isIdentical(ceil(NaN(0xABC)), NaN(0xABC))); | |
| 1766 } | |
| 1767 | |
| 1768 /** | |
| 1769 * Return the value of x rounded to the nearest integer. | |
| 1770 * If the fractional part of x is exactly 0.5, the return value is rounded to | |
| 1771 * the even integer. | |
| 1772 */ | |
| 1773 real round(real x) | |
| 1774 { | |
| 1775 return tango.stdc.math.roundl(x); | |
| 1776 } | |
| 1777 | |
| 1778 debug(UnitTest) { | |
| 1779 unittest { | |
| 1780 assert(isIdentical(round(NaN(0xABC)), NaN(0xABC))); | |
| 1781 } | |
| 1782 } | |
| 1783 | |
| 1784 /** | |
| 1785 * Returns the integer portion of x, dropping the fractional portion. | |
| 1786 * | |
| 1787 * This is also known as "chop" rounding. | |
| 1788 */ | |
| 1789 real trunc(real x) | |
| 1790 { | |
| 1791 return tango.stdc.math.truncl(x); | |
| 1792 } | |
| 1793 | |
| 1794 debug(UnitTest) { | |
| 1795 unittest { | |
| 1796 assert(isIdentical(trunc(NaN(0xABC)), NaN(0xABC))); | |
| 1797 } | |
| 1798 } | |
| 1799 | |
| 1800 /** | |
| 1801 * Rounds x to the nearest int or long. | |
| 1802 * | |
| 1803 * This is generally the fastest method to convert a floating-point number | |
| 1804 * to an integer. Note that the results from this function | |
| 1805 * depend on the rounding mode, if the fractional part of x is exactly 0.5. | |
| 1806 * If using the default rounding mode (ties round to even integers) | |
| 1807 * rndint(4.5) == 4, rndint(5.5)==6. | |
| 1808 */ | |
| 1809 int rndint(real x) | |
| 1810 { | |
| 1811 version(DigitalMars_D_InlineAsm_X86) | |
| 1812 { | |
| 1813 int n; | |
| 1814 asm | |
| 1815 { | |
| 1816 fld x; | |
| 1817 fistp n; | |
| 1818 } | |
| 1819 return n; | |
| 1820 } | |
| 1821 else | |
| 1822 { | |
| 1823 return tango.stdc.math.lrintl(x); | |
| 1824 } | |
| 1825 } | |
| 1826 | |
| 1827 /** ditto */ | |
| 1828 long rndlong(real x) | |
| 1829 { | |
| 1830 version(DigitalMars_D_InlineAsm_X86) | |
| 1831 { | |
| 1832 long n; | |
| 1833 asm | |
| 1834 { | |
| 1835 fld x; | |
| 1836 fistp n; | |
| 1837 } | |
| 1838 return n; | |
| 1839 } | |
| 1840 else | |
| 1841 { | |
| 1842 return tango.stdc.math.llrintl(x); | |
| 1843 } | |
| 1844 } | |
| 1845 | |
| 1846 debug(UnitTest) { | |
| 1847 version(D_InlineAsm_X86) { // Won't work for anything else yet | |
| 1848 | |
| 1849 unittest { | |
| 1850 | |
| 1851 int r = getIeeeRounding; | |
| 1852 assert(r==RoundingMode.ROUNDTONEAREST); | |
| 1853 real b = 5.5; | |
| 1854 int cnear = tango.math.Math.rndint(b); | |
| 1855 assert(cnear == 6); | |
| 1856 auto oldrounding = setIeeeRounding(RoundingMode.ROUNDDOWN); | |
| 1857 scope (exit) setIeeeRounding(oldrounding); | |
| 1858 | |
| 1859 assert(getIeeeRounding==RoundingMode.ROUNDDOWN); | |
| 1860 | |
| 1861 int cdown = tango.math.Math.rndint(b); | |
| 1862 assert(cdown==5); | |
| 1863 } | |
| 1864 | |
| 1865 unittest { | |
| 1866 // Check that the previous test correctly restored the rounding mode | |
| 1867 assert(getIeeeRounding==RoundingMode.ROUNDTONEAREST); | |
| 1868 } | |
| 1869 } | |
| 1870 } |