math.odin 18 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821
  1. package math
  2. import "intrinsics"
  3. _ :: intrinsics;
  4. Float_Class :: enum {
  5. Normal, // an ordinary nonzero floating point value
  6. Subnormal, // a subnormal floating point value
  7. Zero, // zero
  8. Neg_Zero, // the negative zero
  9. NaN, // Not-A-Number (NaN)
  10. Inf, // positive infinity
  11. Neg_Inf // negative infinity
  12. };
  13. TAU :: 6.28318530717958647692528676655900576;
  14. PI :: 3.14159265358979323846264338327950288;
  15. E :: 2.71828182845904523536;
  16. τ :: TAU;
  17. π :: PI;
  18. e :: E;
  19. SQRT_TWO :: 1.41421356237309504880168872420969808;
  20. SQRT_THREE :: 1.73205080756887729352744634150587236;
  21. SQRT_FIVE :: 2.23606797749978969640917366873127623;
  22. LN2 :: 0.693147180559945309417232121458176568;
  23. LN10 :: 2.30258509299404568401799145468436421;
  24. MAX_F64_PRECISION :: 16; // Maximum number of meaningful digits after the decimal point for 'f64'
  25. MAX_F32_PRECISION :: 8; // Maximum number of meaningful digits after the decimal point for 'f32'
  26. RAD_PER_DEG :: TAU/360.0;
  27. DEG_PER_RAD :: 360.0/TAU;
  28. @(default_calling_convention="none")
  29. foreign _ {
  30. @(link_name="llvm.sqrt.f32")
  31. sqrt_f32 :: proc(x: f32) -> f32 ---;
  32. @(link_name="llvm.sqrt.f64")
  33. sqrt_f64 :: proc(x: f64) -> f64 ---;
  34. @(link_name="llvm.sin.f32")
  35. sin_f32 :: proc(θ: f32) -> f32 ---;
  36. @(link_name="llvm.sin.f64")
  37. sin_f64 :: proc(θ: f64) -> f64 ---;
  38. @(link_name="llvm.cos.f32")
  39. cos_f32 :: proc(θ: f32) -> f32 ---;
  40. @(link_name="llvm.cos.f64")
  41. cos_f64 :: proc(θ: f64) -> f64 ---;
  42. @(link_name="llvm.pow.f32")
  43. pow_f32 :: proc(x, power: f32) -> f32 ---;
  44. @(link_name="llvm.pow.f64")
  45. pow_f64 :: proc(x, power: f64) -> f64 ---;
  46. @(link_name="llvm.fmuladd.f32")
  47. fmuladd_f32 :: proc(a, b, c: f32) -> f32 ---;
  48. @(link_name="llvm.fmuladd.f64")
  49. fmuladd_f64 :: proc(a, b, c: f64) -> f64 ---;
  50. @(link_name="llvm.log.f32")
  51. ln_f32 :: proc(x: f32) -> f32 ---;
  52. @(link_name="llvm.log.f64")
  53. ln_f64 :: proc(x: f64) -> f64 ---;
  54. @(link_name="llvm.exp.f32")
  55. exp_f32 :: proc(x: f32) -> f32 ---;
  56. @(link_name="llvm.exp.f64")
  57. exp_f64 :: proc(x: f64) -> f64 ---;
  58. @(link_name="llvm.ldexp.f32")
  59. ldexp_f32 :: proc(val: f32, exp: i32) -> f32 ---;
  60. @(link_name="llvm.ldexp.f64")
  61. ldexp_f64 :: proc(val: f64, exp: i32) -> f64 ---;
  62. }
  63. sqrt :: proc{sqrt_f32, sqrt_f64};
  64. sin :: proc{sin_f32, sin_f64};
  65. cos :: proc{cos_f32, cos_f64};
  66. pow :: proc{pow_f32, pow_f64};
  67. fmuladd :: proc{fmuladd_f32, fmuladd_f64};
  68. ln :: proc{ln_f32, ln_f64};
  69. exp :: proc{exp_f32, exp_f64};
  70. ldexp :: proc{ldexp_f32, ldexp_f64};
  71. log_f32 :: proc(x, base: f32) -> f32 { return ln(x) / ln(base); }
  72. log_f64 :: proc(x, base: f64) -> f64 { return ln(x) / ln(base); }
  73. log :: proc{log_f32, log_f64};
  74. log2_f32 :: proc(x: f32) -> f32 { return ln(x)/LN2; }
  75. log2_f64 :: proc(x: f64) -> f64 { return ln(x)/LN2; }
  76. log2 :: proc{log2_f32, log2_f64};
  77. log10_f32 :: proc(x: f32) -> f32 { return ln(x)/LN10; }
  78. log10_f64 :: proc(x: f64) -> f64 { return ln(x)/LN10; }
  79. log10 :: proc{log10_f32, log10_f64};
  80. tan_f32 :: proc "c" (θ: f32) -> f32 { return sin(θ)/cos(θ); }
  81. tan_f64 :: proc "c" (θ: f64) -> f64 { return sin(θ)/cos(θ); }
  82. tan :: proc{tan_f32, tan_f64};
  83. lerp :: proc(a, b: $T, t: $E) -> (x: T) { return a*(1-t) + b*t; }
  84. unlerp_f32 :: proc(a, b, x: f32) -> (t: f32) { return (x-a)/(b-a); }
  85. unlerp_f64 :: proc(a, b, x: f64) -> (t: f64) { return (x-a)/(b-a); }
  86. unlerp :: proc{unlerp_f32, unlerp_f64};
  87. wrap :: proc(x, y: $T) -> T where intrinsics.type_is_numeric(T), !intrinsics.type_is_array(T) {
  88. tmp := mod(x, y);
  89. return y + tmp if tmp < 0 else tmp;
  90. }
  91. angle_diff :: proc(a, b: $T) -> T where intrinsics.type_is_numeric(T), !intrinsics.type_is_array(T) {
  92. dist := wrap(b - a, TAU);
  93. return wrap(dist*2, TAU) - dist;
  94. }
  95. angle_lerp :: proc(a, b, t: $T) -> T where intrinsics.type_is_numeric(T), !intrinsics.type_is_array(T) {
  96. return a + angle_diff(a, b) * t;
  97. }
  98. step :: proc(edge, x: $T) -> T where intrinsics.type_is_numeric(T), !intrinsics.type_is_array(T) {
  99. return 0 if x < edge else 1;
  100. }
  101. smoothstep :: proc(edge0, edge1, x: $T) -> T where intrinsics.type_is_numeric(T), !intrinsics.type_is_array(T) {
  102. t := clamp((x - edge0) / (edge1 - edge0), 0, 1);
  103. return t * t * (3 - 2*t);
  104. }
  105. bias :: proc(t, b: $T) -> T where intrinsics.type_is_numeric(T) {
  106. return t / (((1/b) - 2) * (1 - t) + 1);
  107. }
  108. gain :: proc(t, g: $T) -> T where intrinsics.type_is_numeric(T) {
  109. if t < 0.5 {
  110. return bias(t*2, g)*0.5;
  111. }
  112. return bias(t*2 - 1, 1 - g)*0.5 + 0.5;
  113. }
  114. sign_f32 :: proc(x: f32) -> f32 { return f32(int(0 < x) - int(x < 0)); }
  115. sign_f64 :: proc(x: f64) -> f64 { return f64(int(0 < x) - int(x < 0)); }
  116. sign :: proc{sign_f32, sign_f64};
  117. sign_bit_f32 :: proc(x: f32) -> bool {
  118. return (transmute(u32)x) & (1<<31) != 0;
  119. }
  120. sign_bit_f64 :: proc(x: f64) -> bool {
  121. return (transmute(u64)x) & (1<<63) != 0;
  122. }
  123. sign_bit :: proc{sign_bit_f32, sign_bit_f64};
  124. copy_sign_f32 :: proc(x, y: f32) -> f32 {
  125. ix := transmute(u32)x;
  126. iy := transmute(u32)y;
  127. ix &= 0x7fff_ffff;
  128. ix |= iy & 0x8000_0000;
  129. return transmute(f32)ix;
  130. }
  131. copy_sign_f64 :: proc(x, y: f64) -> f64 {
  132. ix := transmute(u64)x;
  133. iy := transmute(u64)y;
  134. ix &= 0x7fff_ffff_ffff_ffff;
  135. ix |= iy & 0x8000_0000_0000_0000;
  136. return transmute(f64)ix;
  137. }
  138. copy_sign :: proc{copy_sign_f32, copy_sign_f64};
  139. to_radians_f32 :: proc(degrees: f32) -> f32 { return degrees * RAD_PER_DEG; }
  140. to_radians_f64 :: proc(degrees: f64) -> f64 { return degrees * RAD_PER_DEG; }
  141. to_degrees_f32 :: proc(radians: f32) -> f32 { return radians * DEG_PER_RAD; }
  142. to_degrees_f64 :: proc(radians: f64) -> f64 { return radians * DEG_PER_RAD; }
  143. to_radians :: proc{to_radians_f32, to_radians_f64};
  144. to_degrees :: proc{to_degrees_f32, to_degrees_f64};
  145. trunc_f32 :: proc(x: f32) -> f32 {
  146. trunc_internal :: proc(f: f32) -> f32 {
  147. mask :: 0xff;
  148. shift :: 32 - 9;
  149. bias :: 0x7f;
  150. if f < 1 {
  151. switch {
  152. case f < 0: return -trunc_internal(-f);
  153. case f == 0: return f;
  154. case: return 0;
  155. }
  156. }
  157. x := transmute(u32)f;
  158. e := (x >> shift) & mask - bias;
  159. if e < shift {
  160. x &= ~(1 << (shift-e)) - 1;
  161. }
  162. return transmute(f32)x;
  163. }
  164. switch classify(x) {
  165. case .Zero, .Neg_Zero, .NaN, .Inf, .Neg_Inf:
  166. return x;
  167. case .Normal, .Subnormal: // carry on
  168. }
  169. return trunc_internal(x);
  170. }
  171. trunc_f64 :: proc(x: f64) -> f64 {
  172. trunc_internal :: proc(f: f64) -> f64 {
  173. mask :: 0x7ff;
  174. shift :: 64 - 12;
  175. bias :: 0x3ff;
  176. if f < 1 {
  177. switch {
  178. case f < 0: return -trunc_internal(-f);
  179. case f == 0: return f;
  180. case: return 0;
  181. }
  182. }
  183. x := transmute(u64)f;
  184. e := (x >> shift) & mask - bias;
  185. if e < shift {
  186. x &= ~(1 << (shift-e)) - 1;
  187. }
  188. return transmute(f64)x;
  189. }
  190. switch classify(x) {
  191. case .Zero, .Neg_Zero, .NaN, .Inf, .Neg_Inf:
  192. return x;
  193. case .Normal, .Subnormal: // carry on
  194. }
  195. return trunc_internal(x);
  196. }
  197. trunc :: proc{trunc_f32, trunc_f64};
  198. round_f32 :: proc(x: f32) -> f32 {
  199. return ceil(x - 0.5) if x < 0 else floor(x + 0.5);
  200. }
  201. round_f64 :: proc(x: f64) -> f64 {
  202. return ceil(x - 0.5) if x < 0 else floor(x + 0.5);
  203. }
  204. round :: proc{round_f32, round_f64};
  205. ceil_f32 :: proc(x: f32) -> f32 { return -floor(-x); }
  206. ceil_f64 :: proc(x: f64) -> f64 { return -floor(-x); }
  207. ceil :: proc{ceil_f32, ceil_f64};
  208. floor_f32 :: proc(x: f32) -> f32 {
  209. if x == 0 || is_nan(x) || is_inf(x) {
  210. return x;
  211. }
  212. if x < 0 {
  213. d, fract := modf(-x);
  214. if fract != 0.0 {
  215. d = d + 1;
  216. }
  217. return -d;
  218. }
  219. d, _ := modf(x);
  220. return d;
  221. }
  222. floor_f64 :: proc(x: f64) -> f64 {
  223. if x == 0 || is_nan(x) || is_inf(x) {
  224. return x;
  225. }
  226. if x < 0 {
  227. d, fract := modf(-x);
  228. if fract != 0.0 {
  229. d = d + 1;
  230. }
  231. return -d;
  232. }
  233. d, _ := modf(x);
  234. return d;
  235. }
  236. floor :: proc{floor_f32, floor_f64};
  237. floor_div :: proc(x, y: $T) -> T
  238. where intrinsics.type_is_integer(T) {
  239. a := x / y;
  240. r := x % y;
  241. if (r > 0 && y < 0) || (r < 0 && y > 0) {
  242. a -= 1;
  243. }
  244. return a;
  245. }
  246. floor_mod :: proc(x, y: $T) -> T
  247. where intrinsics.type_is_integer(T) {
  248. r := x % y;
  249. if (r > 0 && y < 0) || (r < 0 && y > 0) {
  250. r += y;
  251. }
  252. return r;
  253. }
  254. modf_f32 :: proc(x: f32) -> (int: f32, frac: f32) {
  255. shift :: 32 - 8 - 1;
  256. mask :: 0xff;
  257. bias :: 127;
  258. if x < 1 {
  259. switch {
  260. case x < 0:
  261. int, frac = modf(-x);
  262. return -int, -frac;
  263. case x == 0:
  264. return x, x;
  265. }
  266. return 0, x;
  267. }
  268. i := transmute(u32)x;
  269. e := uint(i>>shift)&mask - bias;
  270. if e < shift {
  271. i &~= 1<<(shift-e) - 1;
  272. }
  273. int = transmute(f32)i;
  274. frac = x - int;
  275. return;
  276. }
  277. modf_f64 :: proc(x: f64) -> (int: f64, frac: f64) {
  278. shift :: 64 - 11 - 1;
  279. mask :: 0x7ff;
  280. bias :: 1023;
  281. if x < 1 {
  282. switch {
  283. case x < 0:
  284. int, frac = modf(-x);
  285. return -int, -frac;
  286. case x == 0:
  287. return x, x;
  288. }
  289. return 0, x;
  290. }
  291. i := transmute(u64)x;
  292. e := uint(i>>shift)&mask - bias;
  293. if e < shift {
  294. i &~= 1<<(shift-e) - 1;
  295. }
  296. int = transmute(f64)i;
  297. frac = x - int;
  298. return;
  299. }
  300. modf :: proc{modf_f32, modf_f64};
  301. split_decimal :: modf;
  302. mod_f32 :: proc(x, y: f32) -> (n: f32) {
  303. z := abs(y);
  304. n = remainder(abs(x), z);
  305. if sign(n) < 0 {
  306. n += z;
  307. }
  308. return copy_sign(n, x);
  309. }
  310. mod_f64 :: proc(x, y: f64) -> (n: f64) {
  311. z := abs(y);
  312. n = remainder(abs(x), z);
  313. if sign(n) < 0 {
  314. n += z;
  315. }
  316. return copy_sign(n, x);
  317. }
  318. mod :: proc{mod_f32, mod_f64};
  319. remainder_f32 :: proc(x, y: f32) -> f32 { return x - round(x/y) * y; }
  320. remainder_f64 :: proc(x, y: f64) -> f64 { return x - round(x/y) * y; }
  321. remainder :: proc{remainder_f32, remainder_f64};
  322. gcd :: proc(x, y: $T) -> T
  323. where intrinsics.type_is_ordered_numeric(T) {
  324. x, y := x, y;
  325. for y != 0 {
  326. x %= y;
  327. x, y = y, x;
  328. }
  329. return abs(x);
  330. }
  331. lcm :: proc(x, y: $T) -> T
  332. where intrinsics.type_is_ordered_numeric(T) {
  333. return x / gcd(x, y) * y;
  334. }
  335. frexp_f32 :: proc(x: f32) -> (significand: f32, exponent: int) {
  336. switch {
  337. case x == 0:
  338. return 0, 0;
  339. case x < 0:
  340. significand, exponent = frexp(-x);
  341. return -significand, exponent;
  342. }
  343. ex := trunc(log2(x));
  344. exponent = int(ex);
  345. significand = x / pow(2.0, ex);
  346. if abs(significand) >= 1 {
  347. exponent += 1;
  348. significand /= 2;
  349. }
  350. if exponent == 1024 && significand == 0 {
  351. significand = 0.99999999999999988898;
  352. }
  353. return;
  354. }
  355. frexp_f64 :: proc(x: f64) -> (significand: f64, exponent: int) {
  356. switch {
  357. case x == 0:
  358. return 0, 0;
  359. case x < 0:
  360. significand, exponent = frexp(-x);
  361. return -significand, exponent;
  362. }
  363. ex := trunc(log2(x));
  364. exponent = int(ex);
  365. significand = x / pow(2.0, ex);
  366. if abs(significand) >= 1 {
  367. exponent += 1;
  368. significand /= 2;
  369. }
  370. if exponent == 1024 && significand == 0 {
  371. significand = 0.99999999999999988898;
  372. }
  373. return;
  374. }
  375. frexp :: proc{frexp_f32, frexp_f64};
  376. binomial :: proc(n, k: int) -> int {
  377. switch {
  378. case k <= 0: return 1;
  379. case 2*k > n: return binomial(n, n-k);
  380. }
  381. b := n;
  382. for i in 2..<k {
  383. b = (b * (n+1-i))/i;
  384. }
  385. return b;
  386. }
  387. factorial :: proc(n: int) -> int {
  388. when size_of(int) == size_of(i64) {
  389. @static table := [21]int{
  390. 1,
  391. 1,
  392. 2,
  393. 6,
  394. 24,
  395. 120,
  396. 720,
  397. 5_040,
  398. 40_320,
  399. 362_880,
  400. 3_628_800,
  401. 39_916_800,
  402. 479_001_600,
  403. 6_227_020_800,
  404. 87_178_291_200,
  405. 1_307_674_368_000,
  406. 20_922_789_888_000,
  407. 355_687_428_096_000,
  408. 6_402_373_705_728_000,
  409. 121_645_100_408_832_000,
  410. 2_432_902_008_176_640_000,
  411. };
  412. } else {
  413. @static table := [13]int{
  414. 1,
  415. 1,
  416. 2,
  417. 6,
  418. 24,
  419. 120,
  420. 720,
  421. 5_040,
  422. 40_320,
  423. 362_880,
  424. 3_628_800,
  425. 39_916_800,
  426. 479_001_600,
  427. };
  428. }
  429. assert(n >= 0, "parameter must not be negative");
  430. assert(n < len(table), "parameter is too large to lookup in the table");
  431. return 0;
  432. }
  433. classify_f32 :: proc(x: f32) -> Float_Class {
  434. switch {
  435. case x == 0:
  436. i := transmute(i32)x;
  437. if i < 0 {
  438. return .Neg_Zero;
  439. }
  440. return .Zero;
  441. case x*0.5 == x:
  442. if x < 0 {
  443. return .Neg_Inf;
  444. }
  445. return .Inf;
  446. case !(x == x):
  447. return .NaN;
  448. }
  449. u := transmute(u32)x;
  450. exp := int(u>>23) & (1<<8 - 1);
  451. if exp == 0 {
  452. return .Subnormal;
  453. }
  454. return .Normal;
  455. }
  456. classify_f64 :: proc(x: f64) -> Float_Class {
  457. switch {
  458. case x == 0:
  459. i := transmute(i64)x;
  460. if i < 0 {
  461. return .Neg_Zero;
  462. }
  463. return .Zero;
  464. case x*0.5 == x:
  465. if x < 0 {
  466. return .Neg_Inf;
  467. }
  468. return .Inf;
  469. case !(x == x):
  470. return .NaN;
  471. }
  472. u := transmute(u64)x;
  473. exp := int(u>>52) & (1<<11 - 1);
  474. if exp == 0 {
  475. return .Subnormal;
  476. }
  477. return .Normal;
  478. }
  479. classify :: proc{classify_f32, classify_f64};
  480. is_nan_f32 :: proc(x: f32) -> bool { return classify(x) == .NaN; }
  481. is_nan_f64 :: proc(x: f64) -> bool { return classify(x) == .NaN; }
  482. is_nan :: proc{is_nan_f32, is_nan_f64};
  483. // is_inf reports whether f is an infinity, according to sign.
  484. // If sign > 0, is_inf reports whether f is positive infinity.
  485. // If sign < 0, is_inf reports whether f is negative infinity.
  486. // If sign == 0, is_inf reports whether f is either infinity.
  487. is_inf_f32 :: proc(x: f32, sign: int = 0) -> bool {
  488. class := classify(abs(x));
  489. switch {
  490. case sign > 0:
  491. return class == .Inf;
  492. case sign < 0:
  493. return class == .Neg_Inf;
  494. }
  495. return class == .Inf || class == .Neg_Inf;
  496. }
  497. is_inf_f64 :: proc(x: f64, sign: int = 0) -> bool {
  498. class := classify(abs(x));
  499. switch {
  500. case sign > 0:
  501. return class == .Inf;
  502. case sign < 0:
  503. return class == .Neg_Inf;
  504. }
  505. return class == .Inf || class == .Neg_Inf;
  506. }
  507. is_inf :: proc{is_inf_f32, is_inf_f64};
  508. is_power_of_two :: proc(x: int) -> bool {
  509. return x > 0 && (x & (x-1)) == 0;
  510. }
  511. next_power_of_two :: proc(x: int) -> int {
  512. k := x -1;
  513. when size_of(int) == 8 {
  514. k = k | (k >> 32);
  515. }
  516. k = k | (k >> 16);
  517. k = k | (k >> 8);
  518. k = k | (k >> 4);
  519. k = k | (k >> 2);
  520. k = k | (k >> 1);
  521. k += 1 + int(x <= 0);
  522. return k;
  523. }
  524. sum :: proc(x: $T/[]$E) -> (res: E)
  525. where intrinsics.type_is_numeric(E) {
  526. for i in x {
  527. res += i;
  528. }
  529. return;
  530. }
  531. prod :: proc(x: $T/[]$E) -> (res: E)
  532. where intrinsics.type_is_numeric(E) {
  533. for i in x {
  534. res *= i;
  535. }
  536. return;
  537. }
  538. cumsum_inplace :: proc(x: $T/[]$E) -> T
  539. where intrinsics.type_is_numeric(E) {
  540. for i in 1..<len(x) {
  541. x[i] = x[i-1] + x[i];
  542. }
  543. }
  544. cumsum :: proc(dst, src: $T/[]$E) -> T
  545. where intrinsics.type_is_numeric(E) {
  546. N := min(len(dst), len(src));
  547. if N > 0 {
  548. dst[0] = src[0];
  549. for i in 1..<N {
  550. dst[i] = dst[i-1] + src[i];
  551. }
  552. }
  553. return dst[:N];
  554. }
  555. atan2_f32 :: proc(y, x: f32) -> f32 {
  556. // TODO(bill): Better atan2_f32
  557. return f32(atan2_f64(f64(y), f64(x)));
  558. }
  559. atan2_f64 :: proc(y, x: f64) -> f64 {
  560. // TODO(bill): Faster atan2_f64 if possible
  561. // The original C code:
  562. // Stephen L. Moshier
  563. // [email protected]
  564. NAN :: 0h7fff_ffff_ffff_ffff;
  565. INF :: 0h7FF0_0000_0000_0000;
  566. PI :: 0h4009_21fb_5444_2d18;
  567. atan :: proc(x: f64) -> f64 {
  568. if x == 0 {
  569. return x;
  570. }
  571. if x > 0 {
  572. return s_atan(x);
  573. }
  574. return -s_atan(-x);
  575. }
  576. // s_atan reduces its argument (known to be positive) to the range [0, 0.66] and calls x_atan.
  577. s_atan :: proc(x: f64) -> f64 {
  578. MORE_BITS :: 6.123233995736765886130e-17; // pi/2 = PIO2 + MORE_BITS
  579. TAN3PI08 :: 2.41421356237309504880; // tan(3*pi/8)
  580. if x <= 0.66 {
  581. return x_atan(x);
  582. }
  583. if x > TAN3PI08 {
  584. return PI/2 - x_atan(1/x) + MORE_BITS;
  585. }
  586. return PI/4 + x_atan((x-1)/(x+1)) + 0.5*MORE_BITS;
  587. }
  588. // x_atan evaluates a series valid in the range [0, 0.66].
  589. x_atan :: proc(x: f64) -> f64 {
  590. P0 :: -8.750608600031904122785e-01;
  591. P1 :: -1.615753718733365076637e+01;
  592. P2 :: -7.500855792314704667340e+01;
  593. P3 :: -1.228866684490136173410e+02;
  594. P4 :: -6.485021904942025371773e+01;
  595. Q0 :: +2.485846490142306297962e+01;
  596. Q1 :: +1.650270098316988542046e+02;
  597. Q2 :: +4.328810604912902668951e+02;
  598. Q3 :: +4.853903996359136964868e+02;
  599. Q4 :: +1.945506571482613964425e+02;
  600. z := x * x;
  601. z = z * ((((P0*z+P1)*z+P2)*z+P3)*z + P4) / (((((z+Q0)*z+Q1)*z+Q2)*z+Q3)*z + Q4);
  602. z = x*z + x;
  603. return z;
  604. }
  605. switch {
  606. case is_nan(y) || is_nan(x):
  607. return NAN;
  608. case y == 0:
  609. if x >= 0 && !sign_bit(x) {
  610. return copy_sign(0.0, y);
  611. }
  612. return copy_sign(PI, y);
  613. case x == 0:
  614. return copy_sign(PI*0.5, y);
  615. case is_inf(x, 0):
  616. if is_inf(x, 1) {
  617. if is_inf(y, 0) {
  618. return copy_sign(PI*0.25, y);
  619. }
  620. return copy_sign(0, y);
  621. }
  622. if is_inf(y, 0) {
  623. return copy_sign(PI*0.75, y);
  624. }
  625. return copy_sign(PI, y);
  626. case is_inf(y, 0):
  627. return copy_sign(PI*0.5, y);
  628. }
  629. q := atan(y / x);
  630. if x < 0 {
  631. if q <= 0 {
  632. return q + PI;
  633. }
  634. return q - PI;
  635. }
  636. return q;
  637. }
  638. atan2 :: proc{atan2_f32, atan2_f64};
  639. atan_f32 :: proc(x: f32) -> f32 {
  640. return atan2_f32(x, 1);
  641. }
  642. atan_f64 :: proc(x: f64) -> f64 {
  643. return atan2_f64(x, 1);
  644. }
  645. atan :: proc{atan_f32, atan_f64};
  646. asin_f32 :: proc(x: f32) -> f32 {
  647. return atan2_f32(x, 1 + sqrt_f32(1 - x*x));
  648. }
  649. asin_f64 :: proc(x: f64) -> f64 {
  650. return atan2_f64(x, 1 + sqrt_f64(1 - x*x));
  651. }
  652. asin :: proc{asin_f32, asin_f64};
  653. acos_f32 :: proc(x: f32) -> f32 {
  654. return 2 * atan2_f32(sqrt_f32(1 - x), sqrt_f32(1 + x));
  655. }
  656. acos_f64 :: proc(x: f64) -> f64 {
  657. return 2 * atan2_f64(sqrt_f64(1 - x), sqrt_f64(1 + x));
  658. }
  659. acos :: proc{acos_f32, acos_f64};
  660. sinh_f32 :: proc(x: f32) -> f32 {
  661. return (exp(x) - exp(-x))*0.5;
  662. }
  663. sinh_f64 :: proc(x: f64) -> f64 {
  664. return (exp(x) - exp(-x))*0.5;
  665. }
  666. sinh :: proc{sinh_f32, sinh_f64};
  667. cosh_f32 :: proc(x: f32) -> f32 {
  668. return (exp(x) + exp(-x))*0.5;
  669. }
  670. cosh_f64 :: proc(x: f64) -> f64 {
  671. return (exp(x) + exp(-x))*0.5;
  672. }
  673. cosh :: proc{cosh_f32, cosh_f64};
  674. tanh_f32 :: proc(x: f32) -> f32 {
  675. t := exp(2*x);
  676. return (t - 1) / (t + 1);
  677. }
  678. tanh_f64 :: proc(x: f64) -> f64 {
  679. t := exp(2*x);
  680. return (t - 1) / (t + 1);
  681. }
  682. tanh :: proc{tanh_f32, tanh_f64};
  683. F32_DIG :: 6;
  684. F32_EPSILON :: 1.192092896e-07;
  685. F32_GUARD :: 0;
  686. F32_MANT_DIG :: 24;
  687. F32_MAX :: 3.402823466e+38;
  688. F32_MAX_10_EXP :: 38;
  689. F32_MAX_EXP :: 128;
  690. F32_MIN :: 1.175494351e-38;
  691. F32_MIN_10_EXP :: -37;
  692. F32_MIN_EXP :: -125;
  693. F32_NORMALIZE :: 0;
  694. F32_RADIX :: 2;
  695. F32_ROUNDS :: 1;
  696. F64_DIG :: 15; // # of decimal digits of precision
  697. F64_EPSILON :: 2.2204460492503131e-016; // smallest such that 1.0+F64_EPSILON != 1.0
  698. F64_MANT_DIG :: 53; // # of bits in mantissa
  699. F64_MAX :: 1.7976931348623158e+308; // max value
  700. F64_MAX_10_EXP :: 308; // max decimal exponent
  701. F64_MAX_EXP :: 1024; // max binary exponent
  702. F64_MIN :: 2.2250738585072014e-308; // min positive value
  703. F64_MIN_10_EXP :: -307; // min decimal exponent
  704. F64_MIN_EXP :: -1021; // min binary exponent
  705. F64_RADIX :: 2; // exponent radix
  706. F64_ROUNDS :: 1; // addition rounding: near