home *** CD-ROM | disk | FTP | other *** search
/ Linux Cubed Series 3: Developer Tools / Linux Cubed Series 3 - Developer Tools.iso / devel / lang / lisp / gcl-1.000 / gcl-1 / gcl-1.0 / c / num_sfun.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-05-07  |  9.9 KB  |  560 lines
  1. /*
  2.  Copyright (C) 1994 M. Hagiya, W. Schelter, T. Yuasa
  3.  
  4. This file is part of GNU Common Lisp, herein referred to as GCL
  5.  
  6. GCL is free software; you can redistribute it and/or modify it under
  7. the terms of the GNU LIBRARY GENERAL PUBLIC LICENSE as published by
  8. the Free Software Foundation; either version 2, or (at your option)
  9. any later version.
  10.  
  11. GCL is distributed in the hope that it will be useful, but WITHOUT
  12. ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  13. FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public 
  14. License for more details.
  15.  
  16. You should have received a copy of the GNU Library General Public License 
  17. along with GCL; see the file COPYING.  If not, write to the Free Software
  18. Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
  19.  
  20. */
  21.  
  22. #include "include.h"
  23. #include "num_include.h"
  24.  
  25. object imag_unit, minus_imag_unit, imag_two;
  26.  
  27. int
  28. fixnum_expt(x, y)
  29. {
  30.     int z;
  31.  
  32.     z = 1;
  33.     while (y > 0)
  34.         if (y%2 == 0) {
  35.             x *= x;
  36.             y /= 2;
  37.         } else {
  38.             z *= x;
  39.             --y;
  40.         }
  41.     return(z);
  42. }
  43.  
  44. object
  45. number_exp(x)
  46. object x;
  47. {
  48.     double exp();
  49.  
  50.     switch (type_of(x)) {
  51.  
  52.     case t_fixnum:
  53.     case t_bignum:
  54.     case t_ratio:
  55.         return(make_longfloat((longfloat)exp(number_to_double(x))));
  56.  
  57.     case t_shortfloat:
  58.         return(make_shortfloat((shortfloat)exp((double)(sf(x)))));
  59.  
  60.     case t_longfloat:
  61.         return(make_longfloat(exp(lf(x))));
  62.  
  63.     case t_complex:
  64.     {
  65.         object y, y1;
  66.         object number_sin(), number_cos();
  67.             vs_mark;
  68.     
  69.         y = x->cmp.cmp_imag;
  70.         x = x->cmp.cmp_real;
  71.         x = number_exp(x);
  72.         vs_push(x);
  73.         y1 = number_cos(y);
  74.         vs_push(y1);
  75.         y = number_sin(y);
  76.         vs_push(y);
  77.         y = make_complex(y1, y);
  78.         vs_push(y);
  79.         x = number_times(x, y);
  80.         vs_reset;
  81.         return(x);
  82.     }
  83.  
  84.     default:
  85.         FEwrong_type_argument(Snumber, x);
  86.     }
  87. }
  88.  
  89. object
  90. number_expt(x, y)
  91. object x, y;
  92. {
  93.     enum type tx, ty;
  94.     object z, number_nlog();
  95.     vs_mark;
  96.  
  97.     tx = type_of(x);
  98.     ty = type_of(y);
  99.     if (ty == t_fixnum && fix(y) == 0)
  100.         switch (tx) {
  101.         case t_fixnum:  case t_bignum:  case t_ratio:
  102.             return(small_fixnum(1));
  103.  
  104.         case t_shortfloat:
  105.             return(make_shortfloat((shortfloat)1.0));
  106.  
  107.         case t_longfloat:
  108.             return(make_longfloat(1.0));
  109.  
  110.         case t_complex:
  111.             z = number_expt(x->cmp.cmp_real, y);
  112.             vs_push(z);
  113.             z = make_complex(z, small_fixnum(0));
  114.             vs_reset;
  115.             return(z);
  116.  
  117.         default:
  118.             FEwrong_type_argument(Snumber, x);
  119.         }
  120.     if (number_zerop(x)) {
  121.         if (!number_plusp(ty==t_complex?y->cmp.cmp_real:y))
  122.             FEerror("Cannot raise zero to the power ~S.", 1, y);
  123.         return(number_times(x, y));
  124.     }
  125.     if (ty == t_fixnum || ty == t_bignum) {
  126.         if (number_minusp(y)) {
  127.             z = number_negate(y);
  128.             vs_push(z);
  129.             z = number_expt(x, z);
  130.             vs_push(z);
  131.             z = number_divide(small_fixnum(1), z);
  132.             vs_reset;
  133.             return(z);
  134.         }
  135.         z = small_fixnum(1);
  136.         vs_push(z);
  137.         vs_push(Cnil);
  138.         vs_push(Cnil);
  139.         while (number_plusp(y))
  140.             if (number_evenp(y)) {
  141.                 x = number_times(x, x);
  142.                 vs_top[-1] = x;
  143.                 y = integer_divide1(y, small_fixnum(2));
  144.                 vs_top[-2] = y;
  145.             } else {
  146.                 z = number_times(z, x);
  147.                 vs_top[-3] = z;
  148.                 y = number_minus(y, small_fixnum(1));
  149.                 vs_top[-2] = y;
  150.             }
  151.         vs_reset;
  152.         return(z);
  153.     }
  154.     z = number_nlog(x);
  155.     vs_push(z);
  156.     z = number_times(z, y);
  157.     vs_push(z);
  158.     z = number_exp(z);
  159.     vs_reset;
  160.     return(z);
  161. }
  162.  
  163. object
  164. number_nlog(x)
  165. object x;
  166. {
  167.     double log();
  168.     object r, i, a, p, number_sqrt(), number_atan2();
  169.     vs_mark;
  170.  
  171.     if (type_of(x) == t_complex) {
  172.         r = x->cmp.cmp_real;
  173.         i = x->cmp.cmp_imag;
  174.         goto COMPLEX;
  175.     }
  176.     if (number_zerop(x))
  177.         FEerror("Zero is the logarithmic singularity.", 0);
  178.     if (number_minusp(x)) {
  179.         r = x;
  180.         i = small_fixnum(0);
  181.         goto COMPLEX;
  182.     }
  183.     switch (type_of(x)) {
  184.     case t_fixnum:
  185.     case t_bignum:
  186.     case t_ratio:
  187.         return(make_longfloat(log(number_to_double(x))));
  188.  
  189.     case t_shortfloat:
  190.         return(make_shortfloat((shortfloat)log((double)(sf(x)))));
  191.  
  192.     case t_longfloat:
  193.         return(make_longfloat(log(lf(x))));
  194.  
  195.     default:
  196.         FEwrong_type_argument(Snumber, x);
  197.     }
  198.  
  199. COMPLEX:
  200.     a = number_times(r, r);
  201.     vs_push(a);
  202.     p = number_times(i, i);
  203.     vs_push(p);
  204.     a = number_plus(a, p);
  205.     vs_push(a);
  206.     a = number_nlog(a);
  207.     vs_push(a);
  208.     a = number_divide(a, small_fixnum(2));
  209.     vs_push(a);
  210.     p = number_atan2(i, r);
  211.     vs_push(p);
  212.     x = make_complex(a, p);
  213.     vs_reset;
  214.     return(x);
  215. }
  216.  
  217. object
  218. number_log(x, y)
  219. object x, y;
  220. {
  221.     object z;
  222.     vs_mark;
  223.  
  224.     if (number_zerop(y))
  225.         FEerror("Zero is the logarithmic singularity.", 0);
  226.     if (number_zerop(x))
  227.         return(number_times(x, y));
  228.     x = number_nlog(x);
  229.     vs_push(x);
  230.     y = number_nlog(y);
  231.     vs_push(y);
  232.     z = number_divide(y, x);
  233.     vs_reset;
  234.     return(z);
  235. }
  236.  
  237. object
  238. number_sqrt(x)
  239. object x;
  240. {
  241.     object z;
  242.     double sqrt();
  243.     vs_mark;
  244.  
  245.     if (type_of(x) == t_complex)
  246.         goto COMPLEX;
  247.     if (number_minusp(x))
  248.         goto COMPLEX;
  249.     switch (type_of(x)) {
  250.     case t_fixnum:
  251.     case t_bignum:
  252.     case t_ratio:
  253.         return(make_longfloat(
  254.             (longfloat)sqrt(number_to_double(x))));
  255.  
  256.     case t_shortfloat:
  257.         return(make_shortfloat((shortfloat)sqrt((double)(sf(x)))));
  258.  
  259.     case t_longfloat:
  260.         return(make_longfloat(sqrt(lf(x))));
  261.  
  262.     default:
  263.         FEwrong_type_argument(Snumber, x);
  264.     }
  265.  
  266. COMPLEX:
  267.     z = make_ratio(small_fixnum(1), small_fixnum(2));
  268.     vs_push(z);
  269.     z = number_expt(x, z);
  270.     vs_reset;
  271.     return(z);
  272. }
  273.  
  274. object
  275. number_atan2(y, x)
  276. object y, x;
  277. {
  278.     object z;
  279.     double atan(), dy, dx, dz;
  280.  
  281.     dy = number_to_double(y);
  282.     dx = number_to_double(x);
  283.     if (dx > 0.0)
  284.         if (dy > 0.0)
  285.             dz = atan(dy / dx);
  286.         else if (dy == 0.0)
  287.             dz = 0.0;
  288.         else
  289.             dz = -atan(-dy / dx);
  290.     else if (dx == 0.0)
  291.         if (dy > 0.0)
  292.             dz = PI / 2.0;
  293.         else if (dy == 0.0)
  294.             FEerror("Logarithmic singularity.", 0);
  295.         else
  296.             dz = -PI / 2.0;
  297.     else
  298.         if (dy > 0.0)
  299.             dz = PI - atan(dy / -dx);
  300.         else if (dy == 0.0)
  301.             dz = PI;
  302.         else
  303.             dz = -PI + atan(-dy / -dx);
  304.     if (type_of(x) == t_shortfloat)
  305.       z = make_shortfloat((shortfloat)dz);
  306.     else
  307.       z = make_longfloat(dz);
  308.     return(z);
  309. }
  310.  
  311. object
  312. number_atan(y)
  313. object y;
  314. {
  315.     object z, z1;
  316.         vs_mark;
  317.  
  318.     if (type_of(y) == t_complex) {
  319.         z = number_times(imag_unit, y);
  320.         vs_push(z);
  321.         z = one_plus(z);
  322.         vs_push(z);
  323.         z1 = number_times(y, y);
  324.         vs_push(z1);
  325.         z1 = one_plus(z1);
  326.         vs_push(z1);
  327.         z1 = number_sqrt(z1);
  328.         vs_push(z1);
  329.         z = number_divide(z, z1);
  330.         vs_push(z);
  331.         z = number_nlog(z);
  332.         vs_push(z);
  333.         z = number_times(minus_imag_unit, z);
  334.         vs_reset;
  335.         return(z);
  336.     }
  337.     return(number_atan2(y, small_fixnum(1)));
  338. }
  339.  
  340. object
  341. number_sin(x)
  342. object x;
  343. {
  344.     double sin();
  345.  
  346.     switch (type_of(x)) {
  347.  
  348.     case t_fixnum:
  349.     case t_bignum:
  350.     case t_ratio:
  351.         return(make_longfloat((longfloat)sin(number_to_double(x))));
  352.  
  353.     case t_shortfloat:
  354.         return(make_shortfloat((shortfloat)sin((double)(sf(x)))));
  355.  
  356.     case t_longfloat:
  357.         return(make_longfloat(sin(lf(x))));
  358.  
  359.     case t_complex:
  360.     {
  361.         object    r;
  362.         object    x0, x1, x2;
  363.         vs_mark;
  364.  
  365.         x0 = number_times(imag_unit, x);
  366.         vs_push(x0);
  367.         x0 = number_exp(x0);
  368.         vs_push(x0);
  369.         x1 = number_times(minus_imag_unit, x);
  370.         vs_push(x1);
  371.         x1 = number_exp(x1);
  372.         vs_push(x1);
  373.         x2 = number_minus(x0, x1);
  374.         vs_push(x2);
  375.         r = number_divide(x2, imag_two);
  376.  
  377.         vs_reset;
  378.         return(r);
  379.     }
  380.  
  381.     default:
  382.         FEwrong_type_argument(Snumber, x);
  383.  
  384.     }
  385. }
  386.  
  387. object
  388. number_cos(x)
  389. object x;
  390. {
  391.     double cos();
  392.  
  393.     switch (type_of(x)) {
  394.  
  395.     case t_fixnum:
  396.     case t_bignum:
  397.     case t_ratio:
  398.         return(make_longfloat((longfloat)cos(number_to_double(x))));
  399.  
  400.     case t_shortfloat:
  401.         return(make_shortfloat((shortfloat)cos((double)(sf(x)))));
  402.  
  403.     case t_longfloat:
  404.         return(make_longfloat(cos(lf(x))));
  405.  
  406.     case t_complex:
  407.     {
  408.         object r;
  409.         object x0, x1, x2;
  410.         vs_mark;
  411.  
  412.         x0 = number_times(imag_unit, x);
  413.         vs_push(x0);
  414.         x0 = number_exp(x0);
  415.         vs_push(x0);
  416.         x1 = number_times(minus_imag_unit, x);
  417.         vs_push(x1);
  418.         x1 = number_exp(x1);
  419.         vs_push(x1);
  420.         x2 = number_plus(x0, x1);
  421.         vs_push(x2);
  422.         r = number_divide(x2, small_fixnum(2));
  423.  
  424.         vs_reset;
  425.         return(r);
  426.     }
  427.  
  428.     default:
  429.         FEwrong_type_argument(Snumber, x);
  430.  
  431.     }
  432. }
  433.  
  434. object
  435. number_tan(x)
  436. object x;
  437. {
  438.     object r, s, c;
  439.     vs_mark;
  440.  
  441.     s = number_sin(x);
  442.     vs_push(s);
  443.     c = number_cos(x);
  444.     vs_push(c);
  445.     if (number_zerop(c) == TRUE)
  446.         FEerror("Cannot compute the tangent of ~S.", 1, x);
  447.     r = number_divide(s, c);
  448.     vs_reset;
  449.     return(r);
  450. }
  451.  
  452. Lexp()
  453. {
  454.     check_arg(1);
  455.     check_type_number(&vs_base[0]);
  456.     vs_base[0] = number_exp(vs_base[0]);
  457. }
  458.  
  459. Lexpt()
  460. {
  461.     check_arg(2);
  462.     check_type_number(&vs_base[0]);
  463.     check_type_number(&vs_base[1]);
  464.     vs_base[0] = number_expt(vs_base[0], vs_base[1]);
  465.     vs_pop;
  466. }
  467.  
  468. Llog()
  469. {
  470.     int narg;
  471.     
  472.     narg = vs_top - vs_base;
  473.     if (narg < 1)
  474.         too_few_arguments();
  475.     else if (narg == 1) {
  476.         check_type_number(&vs_base[0]);
  477.         vs_base[0] = number_nlog(vs_base[0]);
  478.     } else if (narg == 2) {
  479.         check_type_number(&vs_base[0]);
  480.         check_type_number(&vs_base[1]);
  481.         vs_base[0] = number_log(vs_base[1], vs_base[0]);
  482.         vs_pop;
  483.     } else
  484.         too_many_arguments();
  485. }
  486.  
  487. Lsqrt()
  488. {
  489.     check_arg(1);
  490.     check_type_number(&vs_base[0]);
  491.     vs_base[0] = number_sqrt(vs_base[0]);
  492. }
  493.  
  494. Lsin()
  495. {
  496.     check_arg(1);
  497.     check_type_number(&vs_base[0]);
  498.     vs_base[0] = number_sin(vs_base[0]);
  499. }
  500.  
  501. Lcos()
  502. {
  503.     check_arg(1);
  504.     check_type_number(&vs_base[0]);
  505.     vs_base[0] = number_cos(vs_base[0]);
  506. }
  507.  
  508. Ltan()
  509. {
  510.     check_arg(1);
  511.     check_type_number(&vs_base[0]);
  512.     vs_base[0] = number_tan(vs_base[0]);
  513. }
  514.  
  515. Latan()
  516. {
  517.     int narg;
  518.  
  519.     narg = vs_top - vs_base;
  520.     if (narg < 1)
  521.         too_few_arguments();
  522.     if (narg == 1) {
  523.         check_type_number(&vs_base[0]);
  524.         vs_base[0] = number_atan(vs_base[0]);
  525.     } else if (narg == 2) {
  526.         check_type_or_rational_float(&vs_base[0]);
  527.         check_type_or_rational_float(&vs_base[1]);
  528.         vs_base[0] = number_atan2(vs_base[0], vs_base[1]);
  529.         vs_pop;
  530.     } else
  531.         too_many_arguments();
  532. }
  533.  
  534. init_num_sfun()
  535. {
  536.     imag_unit
  537.     = make_complex(make_longfloat((longfloat)0.0),
  538.                make_longfloat((longfloat)1.0));
  539.     enter_mark_origin(&imag_unit);
  540.     minus_imag_unit
  541.     = make_complex(make_longfloat((longfloat)0.0),
  542.                make_longfloat((longfloat)-1.0));
  543.     enter_mark_origin(&minus_imag_unit);
  544.     imag_two
  545.     = make_complex(make_longfloat((longfloat)0.0),
  546.                make_longfloat((longfloat)2.0));
  547.     enter_mark_origin(&imag_two);
  548.  
  549.     make_constant("PI", make_longfloat(PI));
  550.  
  551.     make_function("EXP", Lexp);
  552.     make_function("EXPT", Lexpt);
  553.     make_function("LOG", Llog);
  554.     make_function("SQRT", Lsqrt);
  555.     make_function("SIN", Lsin);
  556.     make_function("COS", Lcos);
  557.     make_function("TAN", Ltan);
  558.     make_function("ATAN", Latan);
  559. }
  560.