1 module ctstdmath; 2 3 import std.math : LN2, LOG2E, PI, PI_2, abs, fabs, sqrt, hypot, poly; 4 unittest 5 { 6 enum a = real(1).abs, b = real(1).fabs, c = real(1).sqrt, d = real(1).hypot(real(1)), e = real(1).poly([real(0)]); 7 } 8 unittest 9 { 10 enum a = real(1).cbrt, b = real(1).exp2, c = real(1).exp, d = real(1).lg, e = real(1).log, f = real(1).expm1; 11 } 12 13 /// 14 real cbrt(in real a) 15 { 16 if (a < 0) 17 return -(-a).cbrt; 18 if (a.isNaN || a.isInfinity) 19 return a; 20 if (a == 0) 21 return 0; 22 // safe initial guess. 23 real x = a.sqrt, y = real.infinity; 24 while (x < y) 25 { 26 y = x; 27 x = (x * 2 + (a / (x * x))) / 3; 28 } 29 return x; 30 } 31 /// 32 unittest 33 { 34 static assert(8.cbrt == 2); 35 static assert(0.cbrt == 0); 36 static assert((-8).cbrt == -2); 37 static assert(real.nan.cbrt.isNaN); 38 static assert(real.infinity.cbrt.isInfinity); 39 } 40 41 /// exponential function. 42 real exp(real x) 43 { 44 return (x * LOG2E).exp2; 45 } 46 47 /// power of 2. 48 real exp2(real x) 49 { 50 if (x.isNaN) return x; 51 if (x < -(14).exp2small) return 0; 52 if (14.exp2small < x) return real.infinity; 53 real ret = 1; 54 if (0 < x) 55 { 56 foreach_reverse (byte i; 0..15) 57 { 58 if (x <= i.exp2small) 59 continue; 60 x -= i.exp2small; 61 ret *= i.exp2exp2; 62 } 63 } 64 else 65 { 66 foreach_reverse (byte i; 0..15) 67 { 68 if (-(i.exp2small) <= x) 69 continue; 70 x += i.exp2small; 71 ret /= i.exp2exp2; 72 } 73 } 74 if (ret == 0 || ret == real.infinity) 75 return ret; 76 if (+0.5 < x) 77 { 78 x -= 1; 79 ret *= 2; 80 } 81 if (x < -0.5) 82 { 83 x += 1; 84 ret /= 2; 85 } 86 return (x * LN2).expsmall * ret; 87 } 88 /// 89 unittest 90 { 91 static assert(real.nan.exp2.isNaN); 92 static assert(real.infinity.exp2.isInfinity); 93 static assert((-real.infinity).exp2 == 0); 94 static assert(0.exp2 == 1); 95 static assert(2.exp2 == 4); 96 static assert(8.exp2 == 256); 97 static assert((-1).exp2 == 0.5); 98 } 99 100 private auto expsmall(in real x) 101 in 102 { 103 assert (-1 < x); 104 assert (x < +1); 105 } 106 body 107 { 108 real ret = 0, nextTerm = 1, old = real.nan; 109 size_t i; 110 while (ret != old) 111 { 112 old = ret; 113 ret += nextTerm; 114 i += 1; 115 nextTerm *= x / i; 116 } 117 return ret; 118 } 119 120 private auto exp2small(byte n) 121 { 122 real ret = 1; 123 if (0 <= n) 124 foreach (i; 0..n) 125 ret *= 2; 126 else 127 foreach (i; 0..-n) 128 ret *= real(0.5); 129 return ret; 130 } 131 private auto doubling(real x, in ubyte n) 132 { 133 foreach (i; 0..n) 134 x *= x; 135 return x; 136 } 137 private auto exp2exp2(byte n) 138 { 139 return (0 <= n) ? 2.doubling(n) : 0.5.doubling(-n); 140 } 141 unittest 142 { 143 static assert(2.exp2small == 4); 144 static assert(2.exp2exp2 == 16); 145 static assert((-2).exp2small * 4 == 1); 146 static assert((-2).exp2exp2 * 16 == 1); 147 } 148 149 /// equivalent to exp(x) - 1, but more accurate when |x| is small. 150 real expm1(real x) 151 { 152 if (!(-1 < x && x < 1)) 153 return x.exp - 1; 154 real ret = 0, nextTerm = 1, old = real.nan; 155 size_t i; 156 while (ret != old) 157 { 158 i += 1; 159 nextTerm *= x / i; 160 old = ret; 161 ret += nextTerm; 162 } 163 return ret; 164 } 165 166 /// Natural logarithm. 167 real log(real x) 168 { 169 return x.lg * LN2; 170 } 171 alias ln = log; /// ditto 172 173 /// Base-2 logarithm. 174 real lg(real x) 175 { 176 real ret = 0; 177 if (x < 0) return real.nan; 178 if (x == 0) return -real.infinity; 179 if (x.isNaN) return x; 180 if (x.isInfinity) return real.infinity; // x = +inf 181 enum left = real(2) / 3, right = real(4) / 3; 182 if (right < x) 183 foreach_reverse (byte i; 0..15) 184 { 185 if (x / i.exp2exp2 < left) 186 continue; 187 x /= i.exp2exp2; 188 ret += i.exp2small; 189 } 190 else if (x < left) 191 foreach_reverse (byte i; 0..15) 192 { 193 if (right < i.exp2exp2 * x) 194 continue; 195 x *= i.exp2exp2; 196 ret -= i.exp2small; 197 } 198 return x.logsmall / LN2 + ret; 199 } 200 /// ditto 201 alias log2 = lg; 202 unittest 203 { 204 static assert((-1).lg.isNaN); 205 static assert(real.infinity.lg.isInfinity); 206 static assert(real.nan.lg.isNaN); 207 static assert(0.lg < 0 && 0.lg.isInfinity); 208 static assert(0.5.lg == -1); 209 static assert(1.lg == 0); 210 static assert(2.lg == 1); 211 } 212 213 private real logsmall(real x) 214 in 215 { 216 assert (0.5 < x); 217 assert (x < 1.5); 218 } 219 body 220 { 221 return log1p(x - 1); 222 } 223 224 /// equivalent to log(x + 1), but more accurate when |x| is small. 225 real log1p(real x) 226 { 227 real ret = 0, power = -1, old = real.nan; 228 size_t i; 229 while (ret != old) 230 { 231 i += 1; 232 power *= -x; 233 old = ret; 234 ret += power / i; 235 } 236 return ret; 237 } 238 239 /// floor function. 240 real floor(real x) 241 { 242 if (x.isNaN) 243 return x; 244 if (x <= -0x1p+63 || +0x1p+63 <= x) 245 return x; 246 if (x < 0) 247 return -(-x).ceil; 248 ulong left = 0, right = 1UL << 63; 249 // binary search [) 250 foreach (i; 0..63) 251 { 252 auto mid = (left & right) + ((left ^ right) >> 1); 253 if (x < mid) 254 right = mid; 255 else 256 left = mid; 257 } 258 return left; 259 } 260 /// 261 unittest 262 { 263 static assert((-0.01).floor == -1); 264 static assert(real(0).floor == 0); 265 static assert(0.01.floor == 0); 266 } 267 268 /// ceiling function. 269 real ceil(real x) 270 { 271 if (x.isNaN) 272 return x; 273 if (x <= -0x1p+63 || +0x1p+63 <= x) 274 return x; 275 if (x <= 0) 276 return -(-x).floor; 277 ulong left = 0, right = 1UL << 63; 278 // binary search (] 279 foreach (i; 0..63) 280 { 281 auto mid = (left & right) + ((left ^ right) >> 1); 282 if (x <= mid) 283 right = mid; 284 else 285 left = mid; 286 } 287 return right; 288 } 289 /// ditto 290 alias ceiling = ceil; 291 /// 292 unittest 293 { 294 static assert((-0.01).ceil == 0); 295 static assert(real(0).ceil == 0); 296 static assert(0.01.ceil == 1); 297 } 298 299 /// reduce x into [-y .. y] modulo 2y. 300 real modHalf(real x, real y) 301 { 302 if (x.isNaN || y.isNaN || y == 0) 303 return real.nan; 304 if (y.isInfinity) 305 return x; 306 if (y < 0) 307 return x.modHalf(-y); 308 if (x < 0) 309 return -(-x).modHalf(y); 310 void dfs(real mod) 311 { 312 if (mod * 2 <= x) 313 return dfs(mod * 2); 314 if (mod <= x) 315 x -= mod; 316 } 317 dfs(y * 2); 318 if (y < x) 319 return x - 2 * y; 320 return x; 321 } 322 /// 323 unittest 324 { 325 static assert(0.modHalf(1) == 0); 326 static assert(0.5.modHalf(1) == 0.5); 327 static assert(1.modHalf(1) == 1); 328 static assert(1.5.modHalf(1) == -0.5); 329 static assert(2.modHalf(1) == 0); 330 static assert(2.5.modHalf(1) == 0.5); 331 } 332 333 /// cosine. 334 real cos(real x) 335 { 336 return (x.modHalf(PI) + PI_2).sin; 337 } 338 /// 339 unittest 340 { 341 static assert(PI_2.cos == 0); 342 } 343 344 /// sine. 345 real sin(real x) 346 { 347 x = x.modHalf(PI); 348 if (x < 0) 349 return -(-x).sin; 350 if (PI < x * 2) 351 return (PI - x).sin; 352 real ret = 0, nextTerm = x, old = real.nan; 353 size_t i; 354 x = -(x * x); 355 while (ret != old) 356 { 357 old = ret; 358 ret += nextTerm; 359 i += 2; 360 nextTerm *= x / (i * (i + 1)); 361 } 362 return ret; 363 } 364 /// 365 unittest 366 { 367 static assert(0.sin == 0); 368 } 369 370 /// tangent. 371 real tan(real x) 372 { 373 return x.sin / x.cos; 374 } 375 /// cotangent. 376 real cot(real x) 377 { 378 return x.cos / x.sin; 379 } 380 /// secant. 381 real sec(real x) 382 { 383 return 1 / x.cos; 384 } 385 /// cosecant. 386 real csc(real x) 387 { 388 return 1 / x.sin; 389 } 390 /// ditto 391 alias cosec = csc; 392 393 /// arctangent. 394 real atan(real x) 395 { 396 if (x < 0) 397 return -(-x).atan; 398 if (x.isInfinity) 399 return PI_2; 400 if (x == 1) 401 return PI_2 / 2; 402 if (1 < x) 403 return PI_2 - (1 / x).atan; 404 if (1 < 3 * x * x) // without this, very slow when |x| ~ 1. 405 return PI_2 / 3 + ((x - 1 / real(3).sqrt) / (1 + x / real(3).sqrt)).atan; 406 real ret = 0, numeratorX = x, old = real.nan; 407 x *= x; 408 real i = 2; 409 while (ret != old) 410 { 411 old = ret; 412 ret += numeratorX * (i + 1 - (i - 1) * x) / (i * i - 1); 413 i += 4; 414 numeratorX *= x * x; 415 } 416 import std.experimental.logger; 417 i.trace; 418 return ret; 419 } 420 421 /// 422 bool isInfinity(real x) 423 { 424 return x < -real.max || real.max < x; 425 } 426 /// 427 unittest 428 { 429 static assert(!isInfinity(real.init)); 430 static assert(!isInfinity(-real.init)); 431 static assert(!isInfinity(real.nan)); 432 static assert(!isInfinity(-real.nan)); 433 static assert(isInfinity(real.infinity)); 434 static assert(isInfinity(-real.infinity)); 435 static assert(isInfinity(-1.0L / 0.0L)); 436 static assert(!isInfinity(real(0))); 437 } 438 439 /// 440 bool isNaN(real x) 441 { 442 return !(0 <= x) && !(x <= 0); 443 } 444 /// 445 unittest 446 { 447 static assert(isNaN(real.nan)); 448 static assert(isNaN(-real.nan)); 449 static assert(!isNaN(real.infinity)); 450 static assert(!isNaN(-real.infinity)); 451 static assert(!isNaN(-1.0L / 0.0L)); 452 static assert(!isNaN(real(0))); 453 }