src/mp-trigonometric.c

   1 /*
   2  * Copyright (C) 1987-2008 Sun Microsystems, Inc. All Rights Reserved.
   3  * Copyright (C) 2008-2011 Robert Ancell.
   4  *
   5  * This program is free software: you can redistribute it and/or modify it under
   6  * the terms of the GNU General Public License as published by the Free Software
   7  * Foundation, either version 2 of the License, or (at your option) any later
   8  * version. See http://www.gnu.org/copyleft/gpl.html the full text of the
   9  * license.
  10  */
  11
  12 #include <stdlib.h>
  13 #include <string.h>
  14 #include <math.h>
  15 #include <libintl.h>
  16
  17 #include "mp.h"
  18 #include "mp-private.h"
  19
  20 static MPNumber pi;
  21 static gboolean have_pi = FALSE;
  22
  23 static int
  24 mp_compare_mp_to_int(const MPNumber *x, int i)
  25 {
  26     MPNumber t;
  27     mp_set_from_integer(i, &t);
  28     return mp_compare_mp_to_mp(x, &t);
  29 }
  30
  31
  32 /* Convert x to radians */
  33 void
  34 convert_to_radians(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
  35 {
  36     MPNumber t1, t2;
  37
  38     switch(unit) {
  39     default:
  40     case MP_RADIANS:
  41         mp_set_from_mp(x, z);
  42         break;
  43
  44     case MP_DEGREES:
  45         mp_get_pi(&t1);
  46         mp_multiply(x, &t1, &t2);
  47         mp_divide_integer(&t2, 180, z);
  48         break;
  49
  50     case MP_GRADIANS:
  51         mp_get_pi(&t1);
  52         mp_multiply(x, &t1, &t2);
  53         mp_divide_integer(&t2, 200, z);
  54         break;
  55     }
  56 }
  57
  58
  59 void
  60 mp_get_pi(MPNumber *z)
  61 {
  62     if (mp_is_zero(&pi)) {
  63         mp_set_from_string("3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679", 10, &pi);
  64         have_pi = TRUE;
  65     }
  66     mp_set_from_mp(&pi, z);
  67 }
  68
  69
  70 void
  71 convert_from_radians(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
  72 {
  73     MPNumber t1, t2;
  74
  75     switch (unit) {
  76     default:
  77     case MP_RADIANS:
  78         mp_set_from_mp(x, z);
  79         break;
  80
  81     case MP_DEGREES:
  82         mp_multiply_integer(x, 180, &t2);
  83         mp_get_pi(&t1);
  84         mp_divide(&t2, &t1, z);
  85         break;
  86
  87     case MP_GRADIANS:
  88         mp_multiply_integer(x, 200, &t2);
  89         mp_get_pi(&t1);
  90         mp_divide(&t2, &t1, z);
  91         break;
  92     }
  93 }
  94
  95
  96 /* z = sin(x) -1 >= x >= 1, do_sin = 1
  97  * z = cos(x) -1 >= x >= 1, do_sin = 0
  98  */
  99 static void
 100 mpsin1(const MPNumber *x, MPNumber *z, int do_sin)
 101 {
 102     int i, b2;
 103     MPNumber t1, t2;
 104
 105     /* sin(0) = 0, cos(0) = 1 */
 106     if (mp_is_zero(x)) {
 107         if (do_sin == 0)
 108             mp_set_from_integer(1, z);
 109         else
 110             mp_set_from_integer(0, z);
 111         return;
 112     }
 113
 114     mp_multiply(x, x, &t2);
 115     if (mp_compare_mp_to_int(&t2, 1) > 0) {
 116         mperr("*** ABS(X) > 1 IN CALL TO MPSIN1 ***");
 117     }
 118
 119     if (do_sin == 0) {
 120         mp_set_from_integer(1, &t1);
 121         mp_set_from_integer(0, z);
 122         i = 1;
 123     } else {
 124         mp_set_from_mp(x, &t1);
 125         mp_set_from_mp(&t1, z);
 126         i = 2;
 127     }
 128
 129     /* Taylor series */
 130     /* POWER SERIES LOOP.  REDUCE T IF POSSIBLE */
 131     b2 = 2 * max(MP_BASE, 64);
 132     do {
 133         if (MP_T + t1.exponent <= 0)
 134             break;
 135
 136         /*  IF I*(I+1) IS NOT REPRESENTABLE AS AN INTEGER, THE FOLLOWING
 137          *  DIVISION BY I*(I+1) HAS TO BE SPLIT UP.
 138          */
 139         mp_multiply(&t2, &t1, &t1);
 140         if (i > b2) {
 141             mp_divide_integer(&t1, -i, &t1);
 142             mp_divide_integer(&t1, i + 1, &t1);
 143         } else {
 144             mp_divide_integer(&t1, -i * (i + 1), &t1);
 145         }
 146         mp_add(&t1, z, z);
 147
 148         i += 2;
 149     } while (t1.sign != 0);
 150
 151     if (do_sin == 0)
 152         mp_add_integer(z, 1, z);
 153 }
 154
 155
 156 static void
 157 mp_sin_real(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
 158 {
 159     int xs;
 160     MPNumber x_radians;
 161
 162     /* sin(0) = 0 */
 163     if (mp_is_zero(x)) {
 164         mp_set_from_integer(0, z);
 165         return;
 166     }
 167
 168     convert_to_radians(x, unit, &x_radians);
 169
 170     xs = x_radians.sign;
 171     mp_abs(&x_radians, &x_radians);
 172
 173     /* USE MPSIN1 IF ABS(X) <= 1 */
 174     if (mp_compare_mp_to_int(&x_radians, 1) <= 0) {
 175         mpsin1(&x_radians, z, 1);
 176     }
 177     /* FIND ABS(X) MODULO 2PI */
 178     else {
 179         mp_get_pi(z);
 180         mp_divide_integer(z, 4, z);
 181         mp_divide(&x_radians, z, &x_radians);
 182         mp_divide_integer(&x_radians, 8, &x_radians);
 183         mp_fractional_component(&x_radians, &x_radians);
 184
 185         /* SUBTRACT 1/2, SAVE SIGN AND TAKE ABS */
 186         mp_add_fraction(&x_radians, -1, 2, &x_radians);
 187         xs = -xs * x_radians.sign;
 188         if (xs == 0) {
 189             mp_set_from_integer(0, z);
 190             return;
 191         }
 192
 193         x_radians.sign = 1;
 194         mp_multiply_integer(&x_radians, 4, &x_radians);
 195
 196         /* IF NOT LESS THAN 1, SUBTRACT FROM 2 */
 197         if (x_radians.exponent > 0)
 198             mp_add_integer(&x_radians, -2, &x_radians);
 199
 200         if (mp_is_zero(&x_radians)) {
 201             mp_set_from_integer(0, z);
 202             return;
 203         }
 204
 205         x_radians.sign = 1;
 206         mp_multiply_integer(&x_radians, 2, &x_radians);
 207
 208         /*  NOW REDUCED TO FIRST QUADRANT, IF LESS THAN PI/4 USE
 209          *  POWER SERIES, ELSE COMPUTE COS OF COMPLEMENT
 210          */
 211         if (x_radians.exponent > 0) {
 212             mp_add_integer(&x_radians, -2, &x_radians);
 213             mp_multiply(&x_radians, z, &x_radians);
 214             mpsin1(&x_radians, z, 0);
 215         } else {
 216             mp_multiply(&x_radians, z, &x_radians);
 217             mpsin1(&x_radians, z, 1);
 218         }
 219     }
 220
 221     z->sign = xs;
 222 }
 223
 224
 225 static void
 226 mp_cos_real(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
 227 {
 228     /* cos(0) = 1 */
 229     if (mp_is_zero(x)) {
 230         mp_set_from_integer(1, z);
 231         return;
 232     }
 233
 234     convert_to_radians(x, unit, z);
 235
 236     /* Use power series if |x| <= 1 */
 237     mp_abs(z, z);
 238     if (mp_compare_mp_to_int(z, 1) <= 0) {
 239         mpsin1(z, z, 0);
 240     } else {
 241         MPNumber t;
 242
 243         /* cos(x) = sin(π/2 - |x|) */
 244         mp_get_pi(&t);
 245         mp_divide_integer(&t, 2, &t);
 246         mp_subtract(&t, z, z);
 247         mp_sin(z, MP_RADIANS, z);
 248     }
 249 }
 250
 251
 252 void
 253 mp_sin(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
 254 {
 255     if (mp_is_complex(x)) {
 256         MPNumber x_real, x_im, z_real, z_im, t;
 257
 258         mp_real_component(x, &x_real);
 259         mp_imaginary_component(x, &x_im);
 260
 261         mp_sin_real(&x_real, unit, &z_real);
 262         mp_cosh(&x_im, &t);
 263         mp_multiply(&z_real, &t, &z_real);
 264
 265         mp_cos_real(&x_real, unit, &z_im);
 266         mp_sinh(&x_im, &t);
 267         mp_multiply(&z_im, &t, &z_im);
 268
 269         mp_set_from_complex(&z_real, &z_im, z);
 270     }
 271     else
 272        mp_sin_real(x, unit, z);
 273 }
 274
 275
 276 void
 277 mp_cos(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
 278 {
 279     if (mp_is_complex(x)) {
 280         MPNumber x_real, x_im, z_real, z_im, t;
 281
 282         mp_real_component(x, &x_real);
 283         mp_imaginary_component(x, &x_im);
 284
 285         mp_cos_real(&x_real, unit, &z_real);
 286         mp_cosh(&x_im, &t);
 287         mp_multiply(&z_real, &t, &z_real);
 288
 289         mp_sin_real(&x_real, unit, &z_im);
 290         mp_sinh(&x_im, &t);
 291         mp_multiply(&z_im, &t, &z_im);
 292         mp_invert_sign(&z_im, &z_im);
 293
 294         mp_set_from_complex(&z_real, &z_im, z);
 295     }
 296     else
 297        mp_cos_real(x, unit, z);
 298 }
 299
 300
 301 void
 302 mp_tan(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
 303 {
 304     MPNumber cos_x, sin_x;
 305
 306     /* Check for undefined values */
 307     mp_cos(x, unit, &cos_x);
 308     if (mp_is_zero(&cos_x)) {
 309         /* Translators: Error displayed when tangent value is undefined */
 310         mperr(_("Tangent is undefined for angles that are multiples of π (180°) from π∕2 (90°)"));
 311         mp_set_from_integer(0, z);
 312         return;
 313     }
 314
 315     /* tan(x) = sin(x) / cos(x) */
 316     mp_sin(x, unit, &sin_x);
 317     mp_divide(&sin_x, &cos_x, z);
 318 }
 319
 320
 321 void
 322 mp_asin(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
 323 {
 324     MPNumber t1, t2;
 325
 326     /* asin⁻¹(0) = 0 */
 327     if (mp_is_zero(x)) {
 328         mp_set_from_integer(0, z);
 329         return;
 330     }
 331
 332     /* sin⁻¹(x) = tan⁻¹(x / √(1 - x²)), |x| < 1 */
 333     if (x->exponent <= 0) {
 334         mp_set_from_integer(1, &t1);
 335         mp_set_from_mp(&t1, &t2);
 336         mp_subtract(&t1, x, &t1);
 337         mp_add(&t2, x, &t2);
 338         mp_multiply(&t1, &t2, &t2);
 339         mp_root(&t2, -2, &t2);
 340         mp_multiply(x, &t2, z);
 341         mp_atan(z, unit, z);
 342         return;
 343     }
 344
 345     /* sin⁻¹(1) = π/2, sin⁻¹(-1) = -π/2 */
 346     mp_set_from_integer(x->sign, &t2);
 347     if (mp_is_equal(x, &t2)) {
 348         mp_get_pi(z);
 349         mp_divide_integer(z, 2 * t2.sign, z);
 350         convert_from_radians(z, unit, z);
 351         return;
 352     }
 353
 354     /* Translators: Error displayed when inverse sine value is undefined */
 355     mperr(_("Inverse sine is undefined for values outside [-1, 1]"));
 356     mp_set_from_integer(0, z);
 357 }
 358
 359
 360 void
 361 mp_acos(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
 362 {
 363     MPNumber t1, t2;
 364     MPNumber MPn1, pi, MPy;
 365
 366     mp_get_pi(&pi);
 367     mp_set_from_integer(1, &t1);
 368     mp_set_from_integer(-1, &MPn1);
 369
 370     if (mp_is_greater_than(x, &t1) || mp_is_less_than(x, &MPn1)) {
 371         /* Translators: Error displayed when inverse cosine value is undefined */
 372         mperr(_("Inverse cosine is undefined for values outside [-1, 1]"));
 373         mp_set_from_integer(0, z);
 374     } else if (mp_is_zero(x)) {
 375         mp_divide_integer(&pi, 2, z);
 376     } else if (mp_is_equal(x, &t1)) {
 377         mp_set_from_integer(0, z);
 378     } else if (mp_is_equal(x, &MPn1)) {
 379         mp_set_from_mp(&pi, z);
 380     } else {
 381         /* cos⁻¹(x) = tan⁻¹(√(1 - x²) / x) */
 382         mp_multiply(x, x, &t2);
 383         mp_subtract(&t1, &t2, &t2);
 384         mp_sqrt(&t2, &t2);
 385         mp_divide(&t2, x, &t2);
 386         mp_atan(&t2, MP_RADIANS, &MPy);
 387         if (x->sign > 0) {
 388             mp_set_from_mp(&MPy, z);
 389         } else {
 390             mp_add(&MPy, &pi, z);
 391         }
 392     }
 393
 394     convert_from_radians(z, unit, z);
 395 }
 396
 397
 398 void
 399 mp_atan(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
 400 {
 401     int i, q;
 402     float rx = 0.0;
 403     MPNumber t1, t2;
 404
 405     if (mp_is_zero(x)) {
 406         mp_set_from_integer(0, z);
 407         return;
 408     }
 409
 410     mp_set_from_mp(x, &t2);
 411     if (abs(x->exponent) <= 2)
 412         rx = mp_cast_to_float(x);
 413
 414     /* REDUCE ARGUMENT IF NECESSARY BEFORE USING SERIES */
 415     q = 1;
 416     while (t2.exponent >= 0)
 417     {
 418         if (t2.exponent == 0 && 2 * (t2.fraction[0] + 1) <= MP_BASE)
 419             break;
 420
 421         q *= 2;
 422
 423         /* t = t / (√(t² + 1) + 1) */
 424         mp_multiply(&t2, &t2, z);
 425         mp_add_integer(z, 1, z);
 426         mp_sqrt(z, z);
 427         mp_add_integer(z, 1, z);
 428         mp_divide(&t2, z, &t2);
 429     }
 430
 431     /* USE POWER SERIES NOW ARGUMENT IN (-0.5, 0.5) */
 432     mp_set_from_mp(&t2, z);
 433     mp_multiply(&t2, &t2, &t1);
 434
 435     /* SERIES LOOP.  REDUCE T IF POSSIBLE. */
 436     for (i = 1; ; i += 2) {
 437         if (MP_T + 2 + t2.exponent <= 1)
 438             break;
 439
 440         mp_multiply(&t2, &t1, &t2);
 441         mp_multiply_fraction(&t2, -i, i + 2, &t2);
 442
 443         mp_add(z, &t2, z);
 444         if (mp_is_zero(&t2))
 445             break;
 446     }
 447
 448     /* CORRECT FOR ARGUMENT REDUCTION */
 449     mp_multiply_integer(z, q, z);
 450
 451     /*  CHECK THAT RELATIVE ERROR LESS THAN 0.01 UNLESS EXPONENT
 452      *  OF X IS LARGE (WHEN ATAN MIGHT NOT WORK)
 453      */
 454     if (abs(x->exponent) <= 2) {
 455         float ry = mp_cast_to_float(z);
 456         /* THE FOLLOWING MESSAGE MAY INDICATE THAT B**(T-1) IS TOO SMALL. */
 457         if (fabs(ry - atan(rx)) >= fabs(ry) * 0.01)
 458             mperr("*** ERROR OCCURRED IN MP_ATAN, RESULT INCORRECT ***");
 459     }
 460
 461     convert_from_radians(z, unit, z);
 462 }
 463
 464
 465 void
 466 mp_sinh(const MPNumber *x, MPNumber *z)
 467 {
 468     MPNumber abs_x;
 469
 470     /* sinh(0) = 0 */
 471     if (mp_is_zero(x)) {
 472         mp_set_from_integer(0, z);
 473         return;
 474     }
 475
 476     /* WORK WITH ABS(X) */
 477     mp_abs(x, &abs_x);
 478
 479     /* If |x| < 1 USE MPEXP TO AVOID CANCELLATION, otherwise IF TOO LARGE MP_EPOWY GIVES ERROR MESSAGE */
 480     if (abs_x.exponent <= 0) {
 481         MPNumber exp_x, a, b;
 482
 483         /* ((e^|x| + 1) * (e^|x| - 1)) / e^|x| */
 484         // FIXME: Solves to e^|x| - e^-|x|, why not lower branch always? */
 485         mp_epowy(&abs_x, &exp_x);
 486         mp_add_integer(&exp_x, 1, &a);
 487         mp_add_integer(&exp_x, -1, &b);
 488         mp_multiply(&a, &b, z);
 489         mp_divide(z, &exp_x, z);
 490     }
 491     else {
 492         MPNumber exp_x;
 493
 494         /* e^|x| - e^-|x| */
 495         mp_epowy(&abs_x, &exp_x);
 496         mp_reciprocal(&exp_x, z);
 497         mp_subtract(&exp_x, z, z);
 498     }
 499
 500     /* DIVIDE BY TWO AND RESTORE SIGN */
 501     mp_divide_integer(z, 2, z);
 502     mp_multiply_integer(z, x->sign, z);
 503 }
 504
 505
 506 void
 507 mp_cosh(const MPNumber *x, MPNumber *z)
 508 {
 509     MPNumber t;
 510
 511     /* cosh(0) = 1 */
 512     if (mp_is_zero(x)) {
 513         mp_set_from_integer(1, z);
 514         return;
 515     }
 516
 517     /* cosh(x) = (e^x + e^-x) / 2 */
 518     mp_abs(x, &t);
 519     mp_epowy(&t, &t);
 520     mp_reciprocal(&t, z);
 521     mp_add(&t, z, z);
 522     mp_divide_integer(z, 2, z);
 523 }
 524
 525
 526 void
 527 mp_tanh(const MPNumber *x, MPNumber *z)
 528 {
 529     float r__1;
 530     MPNumber t;
 531
 532     /* tanh(0) = 0 */
 533     if (mp_is_zero(x)) {
 534         mp_set_from_integer(0, z);
 535         return;
 536     }
 537
 538     mp_abs(x, &t);
 539
 540     /* SEE IF ABS(X) SO LARGE THAT RESULT IS +-1 */
 541     r__1 = (float) MP_T * 0.5 * log((float) MP_BASE);
 542     mp_set_from_float(r__1, z);
 543     if (mp_compare_mp_to_mp(&t, z) > 0) {
 544         mp_set_from_integer(x->sign, z);
 545         return;
 546     }
 547
 548     /* If |x| >= 1/2 use ?, otherwise use ? to avoid cancellation */
 549     /* |tanh(x)| = (e^|2x| - 1) / (e^|2x| + 1) */
 550     mp_multiply_integer(&t, 2, &t);
 551     if (t.exponent > 0) {
 552         mp_epowy(&t, &t);
 553         mp_add_integer(&t, -1, z);
 554         mp_add_integer(&t, 1, &t);
 555         mp_divide(z, &t, z);
 556     } else {
 557         mp_epowy(&t, &t);
 558         mp_add_integer(&t, 1, z);
 559         mp_add_integer(&t, -1, &t);
 560         mp_divide(&t, z, z);
 561     }
 562
 563     /* Restore sign */
 564     z->sign = x->sign * z->sign;
 565 }
 566
 567
 568 void
 569 mp_asinh(const MPNumber *x, MPNumber *z)
 570 {
 571     MPNumber t;
 572
 573     /* sinh⁻¹(x) = ln(x + √(x² + 1)) */
 574     mp_multiply(x, x, &t);
 575     mp_add_integer(&t, 1, &t);
 576     mp_sqrt(&t, &t);
 577     mp_add(x, &t, &t);
 578     mp_ln(&t, z);
 579 }
 580
 581
 582 void
 583 mp_acosh(const MPNumber *x, MPNumber *z)
 584 {
 585     MPNumber t;
 586
 587     /* Check x >= 1 */
 588     mp_set_from_integer(1, &t);
 589     if (mp_is_less_than(x, &t)) {
 590         /* Translators: Error displayed when inverse hyperbolic cosine value is undefined */
 591         mperr(_("Inverse hyperbolic cosine is undefined for values less than one"));
 592         mp_set_from_integer(0, z);
 593         return;
 594     }
 595
 596     /* cosh⁻¹(x) = ln(x + √(x² - 1)) */
 597     mp_multiply(x, x, &t);
 598     mp_add_integer(&t, -1, &t);
 599     mp_sqrt(&t, &t);
 600     mp_add(x, &t, &t);
 601     mp_ln(&t, z);
 602 }
 603
 604
 605 void
 606 mp_atanh(const MPNumber *x, MPNumber *z)
 607 {
 608     MPNumber one, minus_one, n, d;
 609
 610     /* Check -1 <= x <= 1 */
 611     mp_set_from_integer(1, &one);
 612     mp_set_from_integer(-1, &minus_one);
 613     if (mp_is_greater_equal(x, &one) || mp_is_less_equal(x, &minus_one)) {
 614         /* Translators: Error displayed when inverse hyperbolic tangent value is undefined */
 615         mperr(_("Inverse hyperbolic tangent is undefined for values outside [-1, 1]"));
 616         mp_set_from_integer(0, z);
 617         return;
 618     }
 619
 620     /* atanh(x) = 0.5 * ln((1 + x) / (1 - x)) */
 621     mp_add_integer(x, 1, &n);
 622     mp_set_from_mp(x, &d);
 623     mp_invert_sign(&d, &d);
 624     mp_add_integer(&d, 1, &d);
 625     mp_divide(&n, &d, z);
 626     mp_ln(z, z);
 627     mp_divide_integer(z, 2, z);
 628 }