← Back to branch summary

~ubuntu-branches/ubuntu/vivid/gcl/vivid

« back to all changes in this revision

Viewing changes to o/num_sfun.c

  • Committer: Bazaar Package Importer
  • Author(s): Camm Maguire
  • Date: 2002-03-04 14:29:59 UTC
  • Revision ID: james.westby@ubuntu.com-20020304142959-dey14w08kr7lldu3
Tags: upstream-2.5.0.cvs20020219
Import upstream version 2.5.0.cvs20020219

Show diffs side-by-side

added added

removed removed

Lines of Context:
o/num_sfun.c
 
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(sLnumber, 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(sLnumber, 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(sLnumber, 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(sLnumber, 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(sLnumber, 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(sLnumber, 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
}