redoing dev style to be more test centric, from lessons learned with lisp-matrix.
[CommonLispStat.git] / lib / num_sfun.c
blob6e328a3eeb3652d43e1da5c06c7941bd6ebe433e
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 */
7 #include "include.h"
8 #include "num_include.h"
10 object imag_unit, minus_imag_unit, imag_two;
12 int
13 fixnum_expt(x, y)
15 int z;
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;
26 return(z);
29 object
30 number_exp(x)
31 object x;
33 double exp();
35 switch (type_of(x)) {
37 case t_fixnum:
38 case t_bignum:
39 case t_ratio:
40 return(make_longfloat(exp(number_to_double(x))));
42 case t_shortfloat:
43 return(make_shortfloat((shortfloat)exp((double)(sf(x)))));
45 case t_longfloat:
46 return(make_longfloat(exp(lf(x))));
48 case t_complex:
50 object y, y1;
51 object number_sin(), number_cos();
52 vs_mark;
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);
69 default:
70 FEwrong_type_argument(Snumber, x);
74 object
75 number_expt(x, y)
76 object x, y;
78 enum type tx, ty;
79 object z, number_nlog();
80 vs_mark;
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));
89 case t_shortfloat:
90 return(make_shortfloat((shortfloat)1.0));
92 case t_longfloat:
93 return(make_longfloat(1.0));
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);
102 default:
103 FEwrong_type_argument(Snumber, x);
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));
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);
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;
136 vs_reset;
137 return(z);
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);
148 object
149 number_nlog(x)
150 object x;
152 double log();
153 object r, i, a, p, number_sqrt(), number_atan2();
154 vs_mark;
156 if (type_of(x) == t_complex) {
157 r = x->cmp.cmp_real;
158 i = x->cmp.cmp_imag;
159 goto COMPLEX;
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;
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))));
174 case t_shortfloat:
175 return(make_shortfloat((shortfloat)log((double)(sf(x)))));
177 case t_longfloat:
178 return(make_longfloat(log(lf(x))));
180 default:
181 FEwrong_type_argument(Snumber, x);
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);
202 object
203 number_log(x, y)
204 object x, y;
206 object z;
207 vs_mark;
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);
222 object
223 number_sqrt(x)
224 object x;
226 object z;
227 double sqrt();
228 vs_mark;
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))));
240 case t_shortfloat:
241 return(make_shortfloat((shortfloat)sqrt((double)(sf(x)))));
243 case t_longfloat:
244 return(make_longfloat(sqrt(lf(x))));
246 default:
247 FEwrong_type_argument(Snumber, x);
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);
258 object
259 number_atan2(y, x)
260 object y, x;
262 object z;
263 double atan(), dy, dx, dz;
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);
292 object
293 number_atan(y)
294 object y;
296 object z, z1;
297 vs_mark;
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);
318 return(number_atan2(y, small_fixnum(1)));
321 object
322 number_sin(x)
323 object x;
325 double sin();
327 switch (type_of(x)) {
329 case t_fixnum:
330 case t_bignum:
331 case t_ratio:
332 return(make_longfloat(sin(number_to_double(x))));
334 case t_shortfloat:
335 return(make_shortfloat((shortfloat)sin((double)(sf(x)))));
337 case t_longfloat:
338 return(make_longfloat(sin(lf(x))));
340 case t_complex:
342 object r;
343 object x0, x1, x2;
344 vs_mark;
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);
358 vs_reset;
359 return(r);
362 default:
363 FEwrong_type_argument(Snumber, x);
368 object
369 number_cos(x)
370 object x;
372 double cos();
374 switch (type_of(x)) {
376 case t_fixnum:
377 case t_bignum:
378 case t_ratio:
379 return(make_longfloat(cos(number_to_double(x))));
381 case t_shortfloat:
382 return(make_shortfloat((shortfloat)cos((double)(sf(x)))));
384 case t_longfloat:
385 return(make_longfloat(cos(lf(x))));
387 case t_complex:
389 object r;
390 object x0, x1, x2;
391 vs_mark;
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));
405 vs_reset;
406 return(r);
409 default:
410 FEwrong_type_argument(Snumber, x);
415 object
416 number_tan(x)
417 object x;
419 object r, s, c;
420 vs_mark;
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);
433 object
434 number_asin(x)
435 object x;
437 double asin();
438 double dx;
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);
457 if (-1.0 <= dx && dx <= 1.0) return(make_longfloat(asin(dx)));
460 /* treat as complex argument, result */
462 object r;
463 object x0, x1;
464 vs_mark;
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);
484 object
485 number_acos(x)
486 object x;
488 double acos();
489 double dx;
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);
508 if (-1.0 <= dx && dx <= 1.0) return(make_longfloat(acos(dx)));
511 /* treat as complex argument, result */
513 object r;
514 object x0;
515 vs_mark;
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);
535 Lexp()
537 check_arg(1);
538 check_type_number(&vs_base[0]);
539 vs_base[0] = number_exp(vs_base[0]);
542 Lexpt()
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;
551 Llog()
553 int narg;
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();
570 Lsqrt()
572 check_arg(1);
573 check_type_number(&vs_base[0]);
574 vs_base[0] = number_sqrt(vs_base[0]);
577 Lsin()
579 check_arg(1);
580 check_type_number(&vs_base[0]);
581 vs_base[0] = number_sin(vs_base[0]);
584 Lcos()
586 check_arg(1);
587 check_type_number(&vs_base[0]);
588 vs_base[0] = number_cos(vs_base[0]);
591 Ltan()
593 check_arg(1);
594 check_type_number(&vs_base[0]);
595 vs_base[0] = number_tan(vs_base[0]);
598 Latan()
600 int narg;
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();
617 Lasin()
619 check_arg(1);
620 check_type_number(&vs_base[0]);
621 vs_base[0] = number_asin(vs_base[0]);
624 Lacos()
626 check_arg(1);
627 check_type_number(&vs_base[0]);
628 vs_base[0] = number_acos(vs_base[0]);
631 init_num_sfun()
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);
643 make_constant("PI", make_longfloat(PI));
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);