num_sfun.c

   1 /*
   2 (c) Copyright Taiichi Yuasa and Masami Hagiya, 1984.  All rights reserved.
   3 Copying of this file is authorized to users who have executed the true and
   4 proper "License Agreement for Kyoto Common LISP" with SIGLISP.
   5 */
   6
   7 #include "include.h"
   8 #include "num_include.h"
   9
  10 object imag_unit, minus_imag_unit, imag_two;
  11
  12 int
  13 fixnum_expt(x, y)
  14 {
  15         int z;
  16
  17         z = 1;
  18         while (y > 0)
  19                 if (y%2 == 0) {
  20                         x *= x;
  21                         y /= 2;
  22                 } else {
  23                         z *= x;
  24                         --y;
  25                 }
  26         return(z);
  27 }
  28
  29 object
  30 number_exp(x)
  31 object x;
  32 {
  33         double exp();
  34
  35         switch (type_of(x)) {
  36
  37         case t_fixnum:
  38         case t_bignum:
  39         case t_ratio:
  40                 return(make_longfloat(exp(number_to_double(x))));
  41
  42         case t_shortfloat:
  43                 return(make_shortfloat((shortfloat)exp((double)(sf(x)))));
  44
  45         case t_longfloat:
  46                 return(make_longfloat(exp(lf(x))));
  47
  48         case t_complex:
  49         {
  50                 object y, y1;
  51                 object number_sin(), number_cos();
  52                 vs_mark;
  53
  54                 y = x->cmp.cmp_imag;
  55                 x = x->cmp.cmp_real;
  56                 x = number_exp(x);
  57                 vs_push(x);
  58                 y1 = number_cos(y);
  59                 vs_push(y1);
  60                 y = number_sin(y);
  61                 vs_push(y);
  62                 y = make_complex(y1, y);
  63                 vs_push(y);
  64                 x = number_times(x, y);
  65                 vs_reset;
  66                 return(x);
  67         }
  68
  69         default:
  70                 FEwrong_type_argument(Snumber, x);
  71         }
  72 }
  73
  74 object
  75 number_expt(x, y)
  76 object x, y;
  77 {
  78         enum type tx, ty;
  79         object z, number_nlog();
  80         vs_mark;
  81
  82         tx = type_of(x);
  83         ty = type_of(y);
  84         if (ty == t_fixnum && fix(y) == 0)
  85                 switch (tx) {
  86                 case t_fixnum:  case t_bignum:  case t_ratio:
  87                         return(small_fixnum(1));
  88
  89                 case t_shortfloat:
  90                         return(make_shortfloat((shortfloat)1.0));
  91
  92                 case t_longfloat:
  93                         return(make_longfloat(1.0));
  94
  95                 case t_complex:
  96                         z = number_expt(x->cmp.cmp_real, y);
  97                         vs_push(z);
  98                         z = make_complex(z, small_fixnum(0));
  99                         vs_reset;
 100                         return(z);
 101
 102                 default:
 103                         FEwrong_type_argument(Snumber, x);
 104                 }
 105         if (number_zerop(x)) {
 106                 if (!number_plusp(ty==t_complex?y->cmp.cmp_real:y))
 107                         FEerror("Cannot raise zero to the power ~S.", 1, y);
 108                 return(number_times(x, y));
 109         }
 110         if (ty == t_fixnum || ty == t_bignum) {
 111                 if (number_minusp(y)) {
 112                         z = number_negate(y);
 113                         vs_push(z);
 114                         z = number_expt(x, z);
 115                         vs_push(z);
 116                         z = number_divide(small_fixnum(1), z);
 117                         vs_reset;
 118                         return(z);
 119                 }
 120                 z = small_fixnum(1);
 121                 vs_push(z);
 122                 vs_push(Cnil);
 123                 vs_push(Cnil);
 124                 while (number_plusp(y))
 125                         if (number_evenp(y)) {
 126                                 x = number_times(x, x);
 127                                 vs_top[-1] = x;
 128                                 y = integer_divide1(y, small_fixnum(2));
 129                                 vs_top[-2] = y;
 130                         } else {
 131                                 z = number_times(z, x);
 132                                 vs_top[-3] = z;
 133                                 y = number_minus(y, small_fixnum(1));
 134                                 vs_top[-2] = y;
 135                         }
 136                 vs_reset;
 137                 return(z);
 138         }
 139         z = number_nlog(x);
 140         vs_push(z);
 141         z = number_times(z, y);
 142         vs_push(z);
 143         z = number_exp(z);
 144         vs_reset;
 145         return(z);
 146 }
 147
 148 object
 149 number_nlog(x)
 150 object x;
 151 {
 152         double log();
 153         object r, i, a, p, number_sqrt(), number_atan2();
 154         vs_mark;
 155
 156         if (type_of(x) == t_complex) {
 157                 r = x->cmp.cmp_real;
 158                 i = x->cmp.cmp_imag;
 159                 goto COMPLEX;
 160         }
 161         if (number_zerop(x))
 162                 FEerror("Zero is the logarithmic singularity.", 0);
 163         if (number_minusp(x)) {
 164                 r = x;
 165                 i = small_fixnum(0);
 166                 goto COMPLEX;
 167         }
 168         switch (type_of(x)) {
 169         case t_fixnum:
 170         case t_bignum:
 171         case t_ratio:
 172                 return(make_longfloat(log(number_to_double(x))));
 173
 174         case t_shortfloat:
 175                 return(make_shortfloat((shortfloat)log((double)(sf(x)))));
 176
 177         case t_longfloat:
 178                 return(make_longfloat(log(lf(x))));
 179
 180         default:
 181                 FEwrong_type_argument(Snumber, x);
 182         }
 183
 184 COMPLEX:
 185         a = number_times(r, r);
 186         vs_push(a);
 187         p = number_times(i, i);
 188         vs_push(p);
 189         a = number_plus(a, p);
 190         vs_push(a);
 191         a = number_nlog(a);
 192         vs_push(a);
 193         a = number_divide(a, small_fixnum(2));
 194         vs_push(a);
 195         p = number_atan2(i, r);
 196         vs_push(p);
 197         x = make_complex(a, p);
 198         vs_reset;
 199         return(x);
 200 }
 201
 202 object
 203 number_log(x, y)
 204 object x, y;
 205 {
 206         object z;
 207         vs_mark;
 208
 209         if (number_zerop(y))
 210                 FEerror("Zero is the logarithmic singularity.", 0);
 211         if (number_zerop(x))
 212                 return(number_times(x, y));
 213         x = number_nlog(x);
 214         vs_push(x);
 215         y = number_nlog(y);
 216         vs_push(y);
 217         z = number_divide(y, x);
 218         vs_reset;
 219         return(z);
 220 }
 221
 222 object
 223 number_sqrt(x)
 224 object x;
 225 {
 226         object z;
 227         double sqrt();
 228         vs_mark;
 229
 230         if (type_of(x) == t_complex)
 231                 goto COMPLEX;
 232         if (number_minusp(x))
 233                 goto COMPLEX;
 234         switch (type_of(x)) {
 235         case t_fixnum:
 236         case t_bignum:
 237         case t_ratio:
 238                 return(make_longfloat(sqrt(number_to_double(x))));
 239
 240         case t_shortfloat:
 241                 return(make_shortfloat((shortfloat)sqrt((double)(sf(x)))));
 242
 243         case t_longfloat:
 244                 return(make_longfloat(sqrt(lf(x))));
 245
 246         default:
 247                 FEwrong_type_argument(Snumber, x);
 248         }
 249
 250 COMPLEX:
 251         z = make_ratio(small_fixnum(1), small_fixnum(2));
 252         vs_push(z);
 253         z = number_expt(x, z);
 254         vs_reset;
 255         return(z);
 256 }
 257
 258 object
 259 number_atan2(y, x)
 260 object y, x;
 261 {
 262         object z;
 263         double atan(), dy, dx, dz;
 264
 265         dy = number_to_double(y);
 266         dx = number_to_double(x);
 267         if (dx > 0.0)
 268                 if (dy > 0.0)
 269                         dz = atan(dy / dx);
 270                 else if (dy == 0.0)
 271                         dz = 0.0;
 272                 else
 273                         dz = -atan(-dy / dx);
 274         else if (dx == 0.0)
 275                 if (dy > 0.0)
 276                         dz = PI / 2.0;
 277                 else if (dy == 0.0)
 278                         FEerror("Logarithmic singularity.", 0);
 279                 else
 280                         dz = -PI / 2.0;
 281         else
 282                 if (dy > 0.0)
 283                         dz = PI - atan(dy / -dx);
 284                 else if (dy == 0.0)
 285                         dz = PI;
 286                 else
 287                         dz = -PI + atan(-dy / -dx);
 288         z = make_longfloat(dz);
 289         return(z);
 290 }
 291
 292 object
 293 number_atan(y)
 294 object y;
 295 {
 296         object z, z1;
 297         vs_mark;
 298
 299         if (type_of(y) == t_complex) {
 300                 z = number_times(imag_unit, y);
 301                 vs_push(z);
 302                 z = one_plus(z);
 303                 vs_push(z);
 304                 z1 = number_times(y, y);
 305                 vs_push(z1);
 306                 z1 = one_plus(z1);
 307                 vs_push(z1);
 308                 z1 = number_sqrt(z1);
 309                 vs_push(z1);
 310                 z = number_divide(z, z1);
 311                 vs_push(z);
 312                 z = number_nlog(z);
 313                 vs_push(z);
 314                 z = number_times(minus_imag_unit, z);
 315                 vs_reset;
 316                 return(z);
 317         }
 318         return(number_atan2(y, small_fixnum(1)));
 319 }
 320
 321 object
 322 number_sin(x)
 323 object x;
 324 {
 325         double sin();
 326
 327         switch (type_of(x)) {
 328
 329         case t_fixnum:
 330         case t_bignum:
 331         case t_ratio:
 332                 return(make_longfloat(sin(number_to_double(x))));
 333
 334         case t_shortfloat:
 335                 return(make_shortfloat((shortfloat)sin((double)(sf(x)))));
 336
 337         case t_longfloat:
 338                 return(make_longfloat(sin(lf(x))));
 339
 340         case t_complex:
 341         {
 342                 object  r;
 343                 object  x0, x1, x2;
 344                 vs_mark;
 345
 346                 x0 = number_times(imag_unit, x);
 347                 vs_push(x0);
 348                 x0 = number_exp(x0);
 349                 vs_push(x0);
 350                 x1 = number_times(minus_imag_unit, x);
 351                 vs_push(x1);
 352                 x1 = number_exp(x1);
 353                 vs_push(x1);
 354                 x2 = number_minus(x0, x1);
 355                 vs_push(x2);
 356                 r = number_divide(x2, imag_two);
 357
 358                 vs_reset;
 359                 return(r);
 360         }
 361
 362         default:
 363                 FEwrong_type_argument(Snumber, x);
 364
 365         }
 366 }
 367
 368 object
 369 number_cos(x)
 370 object x;
 371 {
 372         double cos();
 373
 374         switch (type_of(x)) {
 375
 376         case t_fixnum:
 377         case t_bignum:
 378         case t_ratio:
 379                 return(make_longfloat(cos(number_to_double(x))));
 380
 381         case t_shortfloat:
 382                 return(make_shortfloat((shortfloat)cos((double)(sf(x)))));
 383
 384         case t_longfloat:
 385                 return(make_longfloat(cos(lf(x))));
 386
 387         case t_complex:
 388         {
 389                 object r;
 390                 object x0, x1, x2;
 391                 vs_mark;
 392
 393                 x0 = number_times(imag_unit, x);
 394                 vs_push(x0);
 395                 x0 = number_exp(x0);
 396                 vs_push(x0);
 397                 x1 = number_times(minus_imag_unit, x);
 398                 vs_push(x1);
 399                 x1 = number_exp(x1);
 400                 vs_push(x1);
 401                 x2 = number_plus(x0, x1);
 402                 vs_push(x2);
 403                 r = number_divide(x2, small_fixnum(2));
 404
 405                 vs_reset;
 406                 return(r);
 407         }
 408
 409         default:
 410                 FEwrong_type_argument(Snumber, x);
 411
 412         }
 413 }
 414
 415 object
 416 number_tan(x)
 417 object x;
 418 {
 419         object r, s, c;
 420         vs_mark;
 421
 422         s = number_sin(x);
 423         vs_push(s);
 424         c = number_cos(x);
 425         vs_push(c);
 426         if (number_zerop(c) == TRUE)
 427                 FEerror("Cannot compute the tangent of ~S.", 1, x);
 428         r = number_divide(s, c);
 429         vs_reset;
 430         return(r);
 431 }
 432
 433 object
 434 number_asin(x)
 435      object x;
 436 {
 437   double asin();
 438   double dx;
 439
 440   /* check for a real argument in [-1,1] */
 441   if (type_of(x) != t_complex) {
 442     switch (type_of(x)) {
 443     case t_fixnum:
 444     case t_bignum:
 445     case t_ratio:
 446       dx = number_to_double(x);
 447       break;
 448     case t_shortfloat:
 449       dx = sf(x);
 450       break;
 451     case t_longfloat:
 452       dx = lf(x);
 453       break;
 454     default:
 455       FEwrong_type_argument(Snumber, x);
 456     }
 457     if (-1.0 <= dx && dx <= 1.0) return(make_longfloat(asin(dx)));
 458   }
 459
 460   /* treat as complex argument, result */
 461   {
 462     object r;
 463     object x0, x1;
 464     vs_mark;
 465
 466     x0 = number_times(x, x);
 467     vs_push(x0);
 468     x0 = number_minus(small_fixnum(1), x0);
 469     vs_push(x0);
 470     x0 = number_sqrt(x0);
 471     vs_push(x0);
 472     x1 = number_times(imag_unit, x);
 473     vs_push(x1);
 474     x0 = number_plus(x0, x1);
 475     vs_push(x0);
 476     x0 = number_nlog(x0);
 477     vs_push(x0);
 478     r = number_times(minus_imag_unit, x0);
 479     vs_reset;
 480     return(r);
 481   }
 482 }
 483
 484 object
 485 number_acos(x)
 486      object x;
 487 {
 488   double acos();
 489   double dx;
 490
 491   /* check for a real argument in [-1,1] */
 492   if (type_of(x) != t_complex) {
 493     switch (type_of(x)) {
 494     case t_fixnum:
 495     case t_bignum:
 496     case t_ratio:
 497       dx = number_to_double(x);
 498       break;
 499     case t_shortfloat:
 500       dx = sf(x);
 501       break;
 502     case t_longfloat:
 503       dx = lf(x);
 504       break;
 505     default:
 506       FEwrong_type_argument(Snumber, x);
 507     }
 508     if (-1.0 <= dx && dx <= 1.0) return(make_longfloat(acos(dx)));
 509   }
 510
 511   /* treat as complex argument, result */
 512   {
 513     object r;
 514     object x0;
 515     vs_mark;
 516
 517     x0 = number_times(x, x);
 518     vs_push(x0);
 519     x0 = number_minus(small_fixnum(1), x0);
 520     vs_push(x0);
 521     x0 = number_sqrt(x0);
 522     vs_push(x0);
 523     x0 = number_times(imag_unit, x0);
 524     vs_push(x0);
 525     x0 = number_plus(x0, x);
 526     vs_push(x0);
 527     x0 = number_nlog(x0);
 528     vs_push(x0);
 529     r = number_times(minus_imag_unit, x0);
 530     vs_reset;
 531     return(r);
 532   }
 533 }
 534
 535 Lexp()
 536 {
 537         check_arg(1);
 538         check_type_number(&vs_base[0]);
 539         vs_base[0] = number_exp(vs_base[0]);
 540 }
 541
 542 Lexpt()
 543 {
 544         check_arg(2);
 545         check_type_number(&vs_base[0]);
 546         check_type_number(&vs_base[1]);
 547         vs_base[0] = number_expt(vs_base[0], vs_base[1]);
 548         vs_pop;
 549 }
 550
 551 Llog()
 552 {
 553         int narg;
 554
 555         narg = vs_top - vs_base;
 556         if (narg < 1)
 557                 too_few_arguments();
 558         else if (narg == 1) {
 559                 check_type_number(&vs_base[0]);
 560                 vs_base[0] = number_nlog(vs_base[0]);
 561         } else if (narg == 2) {
 562                 check_type_number(&vs_base[0]);
 563                 check_type_number(&vs_base[1]);
 564                 vs_base[0] = number_log(vs_base[1], vs_base[0]);
 565                 vs_pop;
 566         } else
 567                 too_many_arguments();
 568 }
 569
 570 Lsqrt()
 571 {
 572         check_arg(1);
 573         check_type_number(&vs_base[0]);
 574         vs_base[0] = number_sqrt(vs_base[0]);
 575 }
 576
 577 Lsin()
 578 {
 579         check_arg(1);
 580         check_type_number(&vs_base[0]);
 581         vs_base[0] = number_sin(vs_base[0]);
 582 }
 583
 584 Lcos()
 585 {
 586         check_arg(1);
 587         check_type_number(&vs_base[0]);
 588         vs_base[0] = number_cos(vs_base[0]);
 589 }
 590
 591 Ltan()
 592 {
 593         check_arg(1);
 594         check_type_number(&vs_base[0]);
 595         vs_base[0] = number_tan(vs_base[0]);
 596 }
 597
 598 Latan()
 599 {
 600         int narg;
 601
 602         narg = vs_top - vs_base;
 603         if (narg < 1)
 604                 too_few_arguments();
 605         if (narg == 1) {
 606                 check_type_number(&vs_base[0]);
 607                 vs_base[0] = number_atan(vs_base[0]);
 608         } else if (narg == 2) {
 609                 check_type_or_rational_float(&vs_base[0]);
 610                 check_type_or_rational_float(&vs_base[1]);
 611                 vs_base[0] = number_atan2(vs_base[0], vs_base[1]);
 612                 vs_pop;
 613         } else
 614                 too_many_arguments();
 615 }
 616
 617 Lasin()
 618 {
 619   check_arg(1);
 620   check_type_number(&vs_base[0]);
 621   vs_base[0] = number_asin(vs_base[0]);
 622 }
 623
 624 Lacos()
 625 {
 626   check_arg(1);
 627   check_type_number(&vs_base[0]);
 628   vs_base[0] = number_acos(vs_base[0]);
 629 }
 630
 631 init_num_sfun()
 632 {
 633         imag_unit
 634         = make_complex(make_longfloat(0.0), make_longfloat(1.0));
 635         enter_mark_origin(&imag_unit);
 636         minus_imag_unit
 637         = make_complex(make_longfloat(0.0), make_longfloat(-1.0));
 638         enter_mark_origin(&minus_imag_unit);
 639         imag_two
 640         = make_complex(make_longfloat(0.0), make_longfloat(2.0));
 641         enter_mark_origin(&imag_two);
 642
 643         make_constant("PI", make_longfloat(PI));
 644
 645         make_function("EXP", Lexp);
 646         make_function("EXPT", Lexpt);
 647         make_function("LOG", Llog);
 648         make_function("SQRT", Lsqrt);
 649         make_function("SIN", Lsin);
 650         make_function("COS", Lcos);
 651         make_function("TAN", Ltan);
 652         make_function("ATAN", Latan);
 653         make_function("ASIN", Lasin);
 654         make_function("ACOS", Lacos);
 655 }