newlib/private/math.ajla

   1 {*
   2  * Copyright (C) 2024 Mikulas Patocka
   3  *
   4  * This file is part of Ajla.
   5  *
   6  * Ajla is free software: you can redistribute it and/or modify it under the
   7  * terms of the GNU General Public License as published by the Free Software
   8  * Foundation, either version 3 of the License, or (at your option) any later
   9  * version.
  10  *
  11  * Ajla is distributed in the hope that it will be useful, but WITHOUT ANY
  12  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  13  * A PARTICULAR PURPOSE. See the GNU General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License along with
  16  * Ajla. If not, see <https://www.gnu.org/licenses/>.
  17  *}
  18
  19 private unit private.math;
  20
  21 fn math_modulo(t : type, implicit cls : class_real_number(t), x y : t) : t;
  22 fn math_power(t : type, implicit cls : class_real_number(t), x y : t) : t;
  23 fn math_ldexp(t : type, implicit cls : class_real_number(t), x y : t) : t;
  24 fn math_atan2(t : type, implicit cls : class_real_number(t), x y : t) : t;
  25 fn math_pi(t : type, implicit cls : class_real_number(t)) : t;
  26 fn math_sqrt(t : type, implicit cls : class_real_number(t), x : t) : t;
  27 fn math_cbrt(t : type, implicit cls : class_real_number(t), x : t) : t;
  28 fn math_sin(t : type, implicit cls : class_real_number(t), x : t) : t;
  29 fn math_cos(t : type, implicit cls : class_real_number(t), x : t) : t;
  30 fn math_tan(t : type, implicit cls : class_real_number(t), x : t) : t;
  31 fn math_asin(t : type, implicit cls : class_real_number(t), x : t) : t;
  32 fn math_acos(t : type, implicit cls : class_real_number(t), x : t) : t;
  33 fn math_atan(t : type, implicit cls : class_real_number(t), x : t) : t;
  34 fn math_sinh(t : type, implicit cls : class_real_number(t), x : t) : t;
  35 fn math_cosh(t : type, implicit cls : class_real_number(t), x : t) : t;
  36 fn math_tanh(t : type, implicit cls : class_real_number(t), x : t) : t;
  37 fn math_asinh(t : type, implicit cls : class_real_number(t), x : t) : t;
  38 fn math_acosh(t : type, implicit cls : class_real_number(t), x : t) : t;
  39 fn math_atanh(t : type, implicit cls : class_real_number(t), x : t) : t;
  40 fn math_exp2(t : type, implicit cls : class_real_number(t), x : t) : t;
  41 fn math_exp(t : type, implicit cls : class_real_number(t), x : t) : t;
  42 fn math_exp10(t : type, implicit cls : class_real_number(t), x : t) : t;
  43 fn math_log2(t : type, implicit cls : class_real_number(t), x : t) : t;
  44 fn math_log(t : type, implicit cls : class_real_number(t), x : t) : t;
  45 fn math_log10(t : type, implicit cls : class_real_number(t), x : t) : t;
  46 fn math_round(t : type, implicit cls : class_real_number(t), x : t) : t;
  47 fn math_ceil(t : type, implicit cls : class_real_number(t), x : t) : t;
  48 fn math_floor(t : type, implicit cls : class_real_number(t), x : t) : t;
  49 fn math_trunc(t : type, implicit cls : class_real_number(t), x : t) : t;
  50 fn math_fract(t : type, implicit cls : class_real_number(t), x : t) : t;
  51 fn math_mantissa(t : type, implicit cls : class_real_number(t), x : t) : t;
  52 fn math_exponent(t : type, implicit cls : class_real_number(t), x : t) : t;
  53
  54 implementation
  55
  56 uses exception;
  57
  58 fn math_modulo(t : type, implicit cls : class_real_number(t), x y : t) : t
  59 [
  60         if is_infinity x or y = 0 then
  61                 return 0./0.;
  62         if is_infinity y then
  63                 return x;
  64         //var t := math_trunc(x / y);
  65         //return x - y * t;
  66         if x = 0 then
  67                 return x;
  68         var neg := sgn(x);
  69         x := abs(x);
  70         y := abs(y);
  71         var bit := 0;
  72         var yy := y;
  73         while yy < x do [
  74                 yy *= 2;
  75                 bit += 1;
  76         ]
  77         while bit >= 0 do [
  78                 if x - yy >= 0 then
  79                         x -= yy;
  80                 bit -= 1;
  81                 if is_infinity yy then
  82                         yy := ldexp(y, bit);
  83                 else
  84                         yy /= 2;
  85         ]
  86         return x * neg;
  87 ]
  88
  89 fn math_power(t : type, implicit cls : class_real_number(t), x y : t) : t
  90 [
  91         if x = 1 then [
  92                 xeval y;
  93                 return 1;
  94         ]
  95         if y = 0 then
  96                 return 1;
  97         if x < 0 then [
  98                 if math_fract(y) <> 0 then [
  99                         if is_infinity x then
 100                                 return select(y < 0, 1.0/0.0, 0);
 101                         return 0./0.;
 102                 ]
 103                 if y < 0 then [
 104                         y := -y;
 105                         x := cls.recip(x);
 106                 ]
 107                 var yi : int := y;
 108                 if is_exception yi then [
 109                         var m : t;
 110                         if is_infinity y then
 111                                 m := 0;
 112                         else
 113                                 m := math_modulo(y, 2);
 114                         if m = 0 then
 115                                 x := abs(x);
 116                         if abs(x) < 1 then
 117                                 return select(is_negative x, 0.0, -0.0);
 118                         else
 119                                 return select(is_negative x, 1.0/0.0, -1.0/0.0);
 120                 ]
 121                 var res := cls.one;
 122                 while yi > 0 do [
 123                         if yi bt 0 then
 124                                 res *= x;
 125                         x *= x;
 126                         yi shr= 1;
 127                 ]
 128                 return res;
 129         ]
 130         return math_exp(math_log(x) * y);
 131 ]
 132
 133 // warning: this logic is duplicated in floating_ldexp
 134 fn math_ldexp(t : type, implicit cls : class_real_number(t), x y : t) : t
 135 [
 136         xeval x, y;
 137         var yi : int := y;
 138         if not is_exception yi, y = yi then [
 139                 if is_infinity x then
 140                         return x;
 141                 var q : t := 1 shl abs(yi);
 142                 if is_infinity q then [
 143                         var y1 := yi div 2;
 144                         var y2 := yi - y1;
 145                         x := math_ldexp(x, y1);
 146                         x := math_ldexp(x, y2);
 147                         return x;
 148                 ]
 149                 if yi >= 0 then
 150                         x *= q;
 151                 else
 152                         x /= q;
 153                 return x;
 154         ]
 155         return x * math_exp2(y);
 156 ]
 157
 158 fn math_atan2(t : type, implicit cls : class_real_number(t), x y : t) : t
 159 [
 160         if x = 0, y = 0 then [
 161                 y := select(is_negative y, 1, -1);
 162         ]
 163         if is_infinity x, is_infinity y then [
 164                 x := sgn(x);
 165                 y := sgn(y);
 166         ]
 167         if y = 0 then [
 168                 if is_negative x = is_negative y then
 169                         return math_pi(t, cls) / 2;
 170                 else
 171                         return -math_pi(t, cls) / 2;
 172         ]
 173         if y < 0 then [
 174                 if not is_negative x then
 175                         return math_pi(t, cls) + math_atan(x / y);
 176                 else
 177                         return -math_pi(t, cls) + math_atan(x / y);
 178         ]
 179         return math_atan(x / y);
 180 ]
 181
 182 fn math_pi(t : type, implicit cls : class_real_number(t)) : t
 183 [
 184         //return math_atan(cls.one) * 4;
 185         var k : int := 0;
 186         var sum := cls.zero;
 187         while true do [
 188                 var q : t := 0.25;
 189                 var den0 := ipower(16, k);
 190                 var den1 := q + 2 * k;
 191                 var den2 := 2 + 4 * k;
 192                 var den3 := 5 + 8 * k;
 193                 var den4 := 6 + 8 * k;
 194                 var old_sum := sum;
 195                 var diff := ((cls.one / den1) - (cls.one / den2) - (cls.one / den3) - (cls.one / den4)) / den0;
 196                 //eval debug("k=" + ntos(k) + ",sum:" + ntos(sum) + ",diff=" + ntos(diff));
 197                 sum += diff;
 198                 if sum = old_sum then
 199                         break;
 200                 k += 1;
 201         ]
 202         return sum;
 203 ]
 204
 205 fn math_sqrt(t : type, implicit cls : class_real_number(t), x : t) : t
 206 [
 207         if x <= 0 then [
 208                 if x < 0 then
 209                         return 0./0.;
 210                 return x;
 211         ]
 212         if is_infinity x then
 213                 return x;
 214         var guess := x;
 215         while true do [
 216                 var last := guess;
 217                 guess -= (guess - x / guess) / 2;
 218                 if guess = last then
 219                         break;
 220         ]
 221         return guess;
 222 ]
 223
 224 fn math_cbrt(t : type, implicit cls : class_real_number(t), x : t) : t
 225 [
 226         if x = 0 then
 227                 return x;
 228         if is_infinity x then
 229                 return x;
 230         var guess := x;
 231         while true do [
 232                 var last := guess;
 233                 guess -= (guess - x / (guess * guess)) * (cls.one / 3);
 234                 if guess = last then
 235                         break;
 236         ]
 237         return guess;
 238 ]
 239
 240 fn math_sin(t : type, implicit cls : class_real_number(t), x : t) : t
 241 [
 242         var m : t := 1;
 243         if x <= 0 then [
 244                 if x = 0 then
 245                         return x;
 246                 x := -x;
 247                 m := -1;
 248         ]
 249         var pi := math_pi(t, cls);
 250         x := math_modulo(x, pi * 2);
 251         if x > pi then [
 252                 x := x - pi;
 253                 m := -m;
 254         ]
 255         if x >= pi / 2 then [
 256                 x := pi - x;
 257         ]
 258         var k := 1;
 259         var xp := x;
 260         var x2 := x * x;
 261         var fact := 1;
 262         var sum := cls.zero;
 263         var s := cls.one;
 264         while true do [
 265                 var diff := xp / fact;
 266                 xp *= x2;
 267                 fact *= (k + 1) * (k + 2);
 268                 var old_sum := sum;
 269                 sum += diff * s;
 270                 if sum = old_sum then
 271                         break;
 272                 k += 2;
 273                 s := -s;
 274         ]
 275         return sum * m;
 276 ]
 277
 278 fn math_cos(t : type, implicit cls : class_real_number(t), x : t) : t
 279 [
 280         return math_sin(x + math_pi(t, cls) / 2);
 281 ]
 282
 283 fn math_tan(t : type, implicit cls : class_real_number(t), x : t) : t
 284 [
 285         return math_sin(x) / math_cos(x);
 286 ]
 287
 288 fn math_asin(t : type, implicit cls : class_real_number(t), x : t) : t
 289 [
 290         if x = -1 or x = 1 then
 291                 return math_pi(t, cls) / 2 * x;
 292         return math_atan(x / math_sqrt(cls.one - x * x));
 293 ]
 294
 295 fn math_acos(t : type, implicit cls : class_real_number(t), x : t) : t
 296 [
 297         return math_pi(t, cls) / 2 - math_asin(x);
 298 ]
 299
 300 fn math_atan(t : type, implicit cls : class_real_number(t), x : t) : t
 301 [
 302         var m := cls.one;
 303         if x <= 0 then [
 304                 if x = 0 then
 305                         return x;
 306                 x := -x;
 307                 m := -m;
 308         ]
 309         if is_infinity x then
 310                 return math_pi(t, cls) / 2 * m;
 311         var a := cls.zero;
 312         var f := 0;
 313         var limit : t := 0.2;
 314         if x > limit then [
 315                 a := math_atan(limit);
 316                 while x > limit do [
 317                         f += 1;
 318                         x := (x - limit) / (x * limit + 1);
 319                 ]
 320         ]
 321         var k := 1;
 322         var xp := x;
 323         var x2 := x * x;
 324         var sum := cls.zero;
 325         var s := cls.one;
 326         while true do [
 327                 var diff := xp / k;
 328                 xp *= x2;
 329                 var old_sum := sum;
 330                 sum += diff * s;
 331                 if old_sum = sum then
 332                         break;
 333                 k += 2;
 334                 s := -s;
 335         ]
 336         return (a * f + sum) * m;
 337 ]
 338
 339 fn math_sinh(t : type, implicit cls : class_real_number(t), x : t) : t
 340 [
 341         return (math_exp(x) - math_exp(-x)) / 2;
 342 ]
 343
 344 fn math_cosh(t : type, implicit cls : class_real_number(t), x : t) : t
 345 [
 346         return (math_exp(x) + math_exp(-x)) / 2;
 347 ]
 348
 349 fn math_tanh(t : type, implicit cls : class_real_number(t), x : t) : t
 350 [
 351         var val1 := math_exp(x) - math_exp(-x);
 352         var val2 := math_exp(x) + math_exp(-x);
 353         if is_infinity val1, is_infinity val2 then
 354                 return sgn(x);
 355         return val1 / val2;
 356 ]
 357
 358 fn math_asinh(t : type, implicit cls : class_real_number(t), x : t) : t
 359 [
 360         if is_infinity x then
 361                 return x;
 362         return math_log(x + math_sqrt(x * x + 1));
 363 ]
 364
 365 fn math_acosh(t : type, implicit cls : class_real_number(t), x : t) : t
 366 [
 367         return math_log(x + math_sqrt(x * x - 1));
 368 ]
 369
 370 fn math_atanh(t : type, implicit cls : class_real_number(t), x : t) : t
 371 [
 372         return math_log((cls.one + x) / (cls.one - x)) / 2;
 373 ]
 374
 375 fn math_exp2(t : type, implicit cls : class_real_number(t), x : t) : t
 376 [
 377         var c2 : t := 2;
 378         return math_exp(x * math_log(c2));
 379 ]
 380
 381 fn math_exp(t : type, implicit cls : class_real_number(t), x : t) : t
 382 [
 383         var neg := false;
 384         if is_infinity x then [
 385                 if x > 0 then
 386                         return x;
 387                 else
 388                         return 0;
 389         ]
 390         if x < 0 then [
 391                 neg := true;
 392                 x := -x;
 393         ]
 394         var f := 0;
 395         while x > 1 do [
 396                 x *= 0.5;
 397                 f += 1;
 398         ]
 399         var k := 1;
 400         var fact := 1;
 401         var xp := cls.one;
 402         var sum := cls.zero;
 403         while true do [
 404                 var diff := xp / fact;
 405                 xp *= x;
 406                 fact *= k;
 407                 var old_sum := sum;
 408                 sum += diff;
 409                 if old_sum = sum then
 410                         break;
 411                 k += 1;
 412         ]
 413         while f > 0 do [
 414                 f -= 1;
 415                 sum *= sum;
 416         ]
 417         if not neg then
 418                 return sum;
 419         else
 420                 return cls.recip(sum);
 421 ]
 422
 423 fn math_exp10(t : type, implicit cls : class_real_number(t), x : t) : t
 424 [
 425         var c2 : t := 10;
 426         return math_exp(x * math_log(c2));
 427 ]
 428
 429 fn math_log2(t : type, implicit cls : class_real_number(t), x : t) : t
 430 [
 431         var c2 : t := 2;
 432         return math_log(x) / math_log(c2);
 433 ]
 434
 435 fn math_log(t : type, implicit cls : class_real_number(t), x : t) : t
 436 [
 437         if x <= 0 then [
 438                 if x = 0 then
 439                         return -1./0.;
 440                 else
 441                         return 0./0.;
 442         ]
 443         if is_infinity x then
 444                 return x;
 445
 446         var f := 2;
 447         while x <= 0.5 or x >= 2 do [
 448                 x := math_sqrt(x);
 449                 f *= 2;
 450         ]
 451
 452         var k := 1;
 453         var sum := cls.zero;
 454         var a := (x - 1) / (x + 1);
 455         var ap := a;
 456         var a2 := a * a;
 457         while true do [
 458                 var diff := ap / k;
 459                 ap *= a2;
 460                 var old_sum := sum;
 461                 sum += diff;
 462                 if old_sum = sum then
 463                         break;
 464                 k += 2;
 465         ]
 466         return sum * f;
 467 ]
 468
 469 fn math_log10(t : type, implicit cls : class_real_number(t), x : t) : t
 470 [
 471         var c10 : t := 10;
 472         return math_log(x) / math_log(c10);
 473 ]
 474
 475 fn math_round(t : type, implicit cls : class_real_number(t), x : t) : t
 476 [
 477         var i : int := x;
 478         if is_exception i then
 479                 return x;
 480         if abs(math_fract(x)) = 0.5 then [
 481                 if i bt 0 = is_negative x then
 482                         return i - (i and 1);
 483         ]
 484         var xx := x + 0.5;
 485         if xx = x + 1 then
 486                 return x;
 487         var fl := math_floor(xx);
 488         if fl = 0, is_negative(x) then
 489                 return -0.0;
 490         return fl;
 491 ]
 492
 493 fn math_ceil(t : type, implicit cls : class_real_number(t), x : t) : t
 494 [
 495         var i : int := x;
 496         if is_exception i then
 497                 return x;
 498         if x > i then
 499                 i += 1;
 500         if i = 0, is_negative(x) then
 501                 return -0.0;
 502         return i;
 503 ]
 504
 505 fn math_floor(t : type, implicit cls : class_real_number(t), x : t) : t
 506 [
 507         var i : int := x;
 508         if is_exception i then
 509                 return x;
 510         if x < i then
 511                 i -= 1;
 512         if i = 0, is_negative(x) then
 513                 return -0.0;
 514         return i;
 515 ]
 516
 517 fn math_trunc(t : type, implicit cls : class_real_number(t), x : t) : t
 518 [
 519         var i : int := x;
 520         if is_exception i then
 521                 return x;
 522         if i = 0, is_negative(x) then
 523                 return -0.0;
 524         return i;
 525 ]
 526
 527 fn math_fract(t : type, implicit cls : class_real_number(t), x : t) : t
 528 [
 529         var i : int := x;
 530         if is_exception i then
 531                 return select(is_negative x, 0.0, -0.0);
 532         var r := x - i;
 533         if r = 0, is_negative(x) then
 534                 return -0.0;
 535         return r;
 536 ]
 537
 538 fn math_mantissa(t : type, implicit cls : class_real_number(t), x : t) : t
 539 [
 540         if is_infinity x then
 541                 return x;
 542         var m := 1;
 543         if x <= 0 then [
 544                 if x = 0 then
 545                         return 0;
 546                 x := -x;
 547                 m := -1;
 548         ]
 549         while x < 0.5 do [
 550                 x *= 2;
 551         ]
 552         while x >= 1 do [
 553                 x *= 0.5;
 554         ]
 555         return x * m;
 556 ]
 557
 558 fn math_exponent(t : type, implicit cls : class_real_number(t), x : t) : t
 559 [
 560         if is_infinity x then
 561                 return 0;
 562         if x <= 0 then [
 563                 if x = 0 then
 564                         return 0;
 565                 x := -x;
 566         ]
 567         var e : int := 0;
 568         while x < 0.5 do [
 569                 x *= 2;
 570                 e -= 1;
 571         ]
 572         while x >= 1 do [
 573                 x *= 0.5;
 574                 e += 1;
 575         ]
 576         return e;
 577 ]