source: branches/0.26.x/abcl/src/org/armedbear/lisp/MathFunctions.java

Last change on this file was 12940, checked in by ehuelsmann, 15 years ago

Fix loss of precision in (expt <non-double> <complex-double>),
fixes last Maxima failure.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 29.0 KB
Line 
1/*
2 * MathFunctions.java
3 *
4 * Copyright (C) 2004-2006 Peter Graves
5 * $Id: MathFunctions.java 12940 2010-10-02 21:39:52Z ehuelsmann $
6 *
7 * This program is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU General Public License
9 * as published by the Free Software Foundation; either version 2
10 * of the License, or (at your option) any later version.
11 *
12 * This program is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with this program; if not, write to the Free Software
19 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
20 *
21 * As a special exception, the copyright holders of this library give you
22 * permission to link this library with independent modules to produce an
23 * executable, regardless of the license terms of these independent
24 * modules, and to copy and distribute the resulting executable under
25 * terms of your choice, provided that you also meet, for each linked
26 * independent module, the terms and conditions of the license of that
27 * module.  An independent module is a module which is not derived from
28 * or based on this library.  If you modify this library, you may extend
29 * this exception to your version of the library, but you are not
30 * obligated to do so.  If you do not wish to do so, delete this
31 * exception statement from your version.
32 */
33
34package org.armedbear.lisp;
35
36import static org.armedbear.lisp.Lisp.*;
37
38public final class MathFunctions
39{
40
41    // Implementation of section 12.1.5.3, which says:
42    // "If the result of any computation would be a complex number whose
43    //  real part is of type rational and whose imaginary part is zero,
44    //  the result is converted to the rational which is the real part."
45    private static final LispObject complexToRealFixup(LispObject result,
46                                                       LispObject arg)
47    {
48        if (result instanceof Complex
49            && ! (arg instanceof Complex)) {
50            Complex c = (Complex)result;
51            LispObject im = c.getImaginaryPart();
52            if (im.zerop())
53                return c.getRealPart();
54        }
55        return result;
56    }
57
58    // ### sin
59    private static final Primitive SIN = new Primitive("sin", "radians")
60    {
61        @Override
62        public LispObject execute(LispObject arg)
63        {
64            return sin(arg);
65        }
66    };
67
68    static LispObject sin(LispObject arg)
69    {
70        if (arg instanceof DoubleFloat)
71            return new DoubleFloat(Math.sin(((DoubleFloat)arg).value));
72        if (arg.realp())
73            return new SingleFloat((float)Math.sin(SingleFloat.coerceToFloat(arg).value));
74        if (arg instanceof Complex) {
75            LispObject n = arg.multiplyBy(Complex.getInstance(Fixnum.ZERO,
76                                                              Fixnum.ONE));
77            LispObject result = exp(n);
78            result = result.subtract(exp(n.multiplyBy(Fixnum.MINUS_ONE)));
79            return result.divideBy(Fixnum.TWO.multiplyBy(Complex.getInstance(Fixnum.ZERO,
80                                                                             Fixnum.ONE)));
81        }
82        return type_error(arg, Symbol.NUMBER);
83    }
84
85    // ### cos
86    private static final Primitive COS = new Primitive("cos", "radians")
87    {
88        @Override
89        public LispObject execute(LispObject arg)
90        {
91            return cos(arg);
92        }
93    };
94
95    static LispObject cos(LispObject arg)
96    {
97        if (arg instanceof DoubleFloat)
98            return new DoubleFloat(Math.cos(((DoubleFloat)arg).value));
99        if (arg.realp())
100            return new SingleFloat((float)Math.cos(SingleFloat.coerceToFloat(arg).value));
101        if (arg instanceof Complex) {
102            LispObject n = arg.multiplyBy(Complex.getInstance(Fixnum.ZERO,
103                                                              Fixnum.ONE));
104            LispObject result = exp(n);
105            result = result.add(exp(n.multiplyBy(Fixnum.MINUS_ONE)));
106            return result.divideBy(Fixnum.TWO);
107        }
108        return type_error(arg, Symbol.NUMBER);
109    }
110
111    // ### tan
112    private static final Primitive TAN = new Primitive("tan", "radians")
113    {
114        @Override
115        public LispObject execute(LispObject arg)
116        {
117            if (arg instanceof DoubleFloat)
118                return new DoubleFloat(Math.tan(((DoubleFloat)arg).value));
119            if (arg.realp())
120                return new SingleFloat((float)Math.tan(SingleFloat.coerceToFloat(arg).value));
121            return sin(arg).divideBy(cos(arg));
122        }
123    };
124
125    // ### asin
126    private static final Primitive ASIN = new Primitive("asin", "number")
127    {
128        @Override
129        public LispObject execute(LispObject arg)
130        {
131            return asin(arg);
132        }
133    };
134
135    static LispObject asin(LispObject arg)
136    {
137        if (arg instanceof SingleFloat) {
138            float f = ((SingleFloat)arg).value;
139            if (Math.abs(f) <= 1)
140                return new SingleFloat((float)Math.asin(f));
141        }
142        if (arg instanceof DoubleFloat) {
143            double d = ((DoubleFloat)arg).value;
144            if (Math.abs(d) <= 1)
145                return new DoubleFloat(Math.asin(d));
146        }
147        LispObject result = arg.multiplyBy(arg);
148        result = Fixnum.ONE.subtract(result);
149        result = sqrt(result);
150        LispObject n = Complex.getInstance(Fixnum.ZERO, Fixnum.ONE);
151        n = n.multiplyBy(arg);
152        result = n.add(result);
153        result = log(result);
154        result = result.multiplyBy(Complex.getInstance(Fixnum.ZERO,
155                                                       Fixnum.MINUS_ONE));
156
157        return complexToRealFixup(result, arg);
158    }
159
160    // ### acos
161    private static final Primitive ACOS = new Primitive("acos", "number")
162    {
163        @Override
164        public LispObject execute(LispObject arg)
165        {
166            return acos(arg);
167        }
168    };
169
170    static LispObject acos(LispObject arg)
171    {
172        if (arg instanceof DoubleFloat) {
173            double d = ((DoubleFloat)arg).value;
174            if (Math.abs(d) <= 1)
175                return new DoubleFloat(Math.acos(d));
176        }
177        if (arg instanceof SingleFloat) {
178            float f = ((SingleFloat)arg).value;
179            if (Math.abs(f) <= 1)
180                return new SingleFloat((float)Math.acos(f));
181        }
182        LispObject result = new DoubleFloat(Math.PI/2);
183        if (!(arg instanceof DoubleFloat)) {
184            if (arg instanceof Complex &&
185                    ((Complex)arg).getRealPart() instanceof DoubleFloat) {
186                    // do nothing; we want to keep the double float value
187            }
188            else
189                result = new SingleFloat((float)((DoubleFloat)result).value);
190        }
191        result = result.subtract(asin(arg));
192
193        return complexToRealFixup(result, arg);
194    }
195
196    // ### atan
197    private static final Primitive ATAN =
198        new Primitive("atan", "number1 &optional number2")
199    {
200        @Override
201        public LispObject execute(LispObject arg)
202        {
203            if (arg.numberp())
204                return atan(arg);
205            return type_error(arg, Symbol.NUMBER);
206        }
207
208        // "If both number1 and number2 are supplied for atan, the result is
209        // the arc tangent of number1/number2."
210
211        // y = +0     x = +0       +0
212        // y = -0     x = +0       -0
213        // y = +0     x = -0       +<PI>
214        // y = -0     x = -0       -<PI>
215        @Override
216        public LispObject execute(LispObject y, LispObject x)
217
218        {
219            if (!y.realp())
220                return type_error(y, Symbol.REAL);
221            if (!x.realp())
222                return type_error(x, Symbol.REAL);
223            double d1, d2;
224            d1 = DoubleFloat.coerceToFloat(y).value;
225            d2 = DoubleFloat.coerceToFloat(x).value;
226            double result = Math.atan2(d1, d2);
227            if (y instanceof DoubleFloat || x instanceof DoubleFloat)
228                return new DoubleFloat(result);
229            else
230                return new SingleFloat((float)result);
231        }
232    };
233
234    static LispObject atan(LispObject arg)
235    {
236        if (arg instanceof Complex) {
237            LispObject im = ((Complex)arg).imagpart;
238            if (im.zerop())
239                return Complex.getInstance(atan(((Complex)arg).realpart),
240                                           im);
241            LispObject result = arg.multiplyBy(arg);
242            result = result.add(Fixnum.ONE);
243            result = Fixnum.ONE.divideBy(result);
244            result = sqrt(result);
245            LispObject n = Complex.getInstance(Fixnum.ZERO, Fixnum.ONE);
246            n = n.multiplyBy(arg);
247            n = n.add(Fixnum.ONE);
248            result = n.multiplyBy(result);
249            result = log(result);
250            result = result.multiplyBy(Complex.getInstance(Fixnum.ZERO, Fixnum.MINUS_ONE));
251            return result;
252        }
253        if (arg instanceof DoubleFloat)
254            return new DoubleFloat(Math.atan(((DoubleFloat)arg).value));
255        return new SingleFloat((float)Math.atan(SingleFloat.coerceToFloat(arg).value));
256    }
257
258    // ### sinh
259    private static final Primitive SINH = new Primitive("sinh", "number")
260    {
261        @Override
262        public LispObject execute(LispObject arg)
263        {
264            return sinh(arg);
265        }
266    };
267
268    static LispObject sinh(LispObject arg)
269    {
270        if (arg instanceof Complex) {
271            LispObject im = ((Complex)arg).getImaginaryPart();
272            if (im.zerop())
273                return Complex.getInstance(sinh(((Complex)arg).getRealPart()),
274                                           im);
275        }
276        if (arg instanceof SingleFloat) {
277            double d = Math.sinh(((SingleFloat)arg).value);
278            return new SingleFloat((float)d);
279        } else if (arg instanceof DoubleFloat) {
280            double d = Math.sinh(((DoubleFloat)arg).value);
281            return new DoubleFloat(d);
282        }
283        LispObject result = exp(arg);
284        result = result.subtract(exp(arg.multiplyBy(Fixnum.MINUS_ONE)));
285        result = result.divideBy(Fixnum.TWO);
286
287        return complexToRealFixup(result, arg);
288    }
289
290    // ### cosh
291    private static final Primitive COSH = new Primitive("cosh", "number")
292    {
293        @Override
294        public LispObject execute(LispObject arg)
295        {
296            return cosh(arg);
297        }
298    };
299
300    static LispObject cosh(LispObject arg)
301    {
302        if (arg instanceof Complex) {
303            LispObject im = ((Complex)arg).getImaginaryPart();
304            if (im.zerop())
305                return Complex.getInstance(cosh(((Complex)arg).getRealPart()),
306                                           im);
307        }
308        if (arg instanceof SingleFloat) {
309            double d = Math.cosh(((SingleFloat)arg).value);
310            return new SingleFloat((float)d);
311        } else if (arg instanceof DoubleFloat) {
312            double d = Math.cosh(((DoubleFloat)arg).value);
313            return new DoubleFloat(d);
314        }
315        LispObject result = exp(arg);
316        result = result.add(exp(arg.multiplyBy(Fixnum.MINUS_ONE)));
317        result = result.divideBy(Fixnum.TWO);
318
319        return complexToRealFixup(result, arg);
320    }
321
322    // ### tanh
323    private static final Primitive TANH = new Primitive("tanh", "number")
324    {
325        @Override
326        public LispObject execute(LispObject arg)
327        {
328            if (arg instanceof SingleFloat) {
329                double d = Math.tanh(((SingleFloat)arg).value);
330                return new SingleFloat((float)d);
331            } else if (arg instanceof DoubleFloat) {
332                double d = Math.tanh(((DoubleFloat)arg).value);
333                return new DoubleFloat(d);
334            }
335            return sinh(arg).divideBy(cosh(arg));
336        }
337    };
338
339    // ### asinh
340    private static final Primitive ASINH = new Primitive("asinh", "number")
341    {
342        @Override
343        public LispObject execute(LispObject arg)
344        {
345            return asinh(arg);
346        }
347    };
348
349    static LispObject asinh(LispObject arg)
350    {
351        if (arg instanceof Complex) {
352            LispObject im = ((Complex)arg).getImaginaryPart();
353            if (im.zerop())
354                return Complex.getInstance(asinh(((Complex)arg).getRealPart()),
355                                           im);
356        }
357        LispObject result = arg.multiplyBy(arg);
358        result = Fixnum.ONE.add(result);
359        result = sqrt(result);
360        result = result.add(arg);
361        result = log(result);
362
363        return complexToRealFixup(result, arg);
364    }
365
366    // ### acosh
367    private static final Primitive ACOSH = new Primitive("acosh", "number")
368    {
369        @Override
370        public LispObject execute(LispObject arg)
371        {
372            return acosh(arg);
373        }
374    };
375
376    static LispObject acosh(LispObject arg)
377    {
378        if (arg instanceof Complex) {
379            LispObject im = ((Complex)arg).getImaginaryPart();
380            if (im.zerop())
381                return Complex.getInstance(acosh(((Complex)arg).getRealPart()),
382                                           im);
383        }
384        LispObject n1 = arg.add(Fixnum.ONE);
385        n1 = n1.divideBy(Fixnum.TWO);
386        n1 = sqrt(n1);
387        LispObject n2 = arg.subtract(Fixnum.ONE);
388        n2 = n2.divideBy(Fixnum.TWO);
389        n2 = sqrt(n2);
390        LispObject result = n1.add(n2);
391        result = log(result);
392        result = result.multiplyBy(Fixnum.TWO);
393
394        return complexToRealFixup(result, arg);
395    }
396
397    // ### atanh
398    private static final Primitive ATANH = new Primitive("atanh", "number")
399    {
400        @Override
401        public LispObject execute(LispObject arg)
402        {
403            return atanh(arg);
404        }
405    };
406
407    static LispObject atanh(LispObject arg)
408    {
409        if (arg instanceof Complex) {
410            LispObject im = ((Complex)arg).getImaginaryPart();
411            if (im.zerop())
412                return Complex.getInstance(atanh(((Complex)arg).getRealPart()),
413                                           im);
414        }
415        LispObject n1 = log(Fixnum.ONE.add(arg));
416        LispObject n2 = log(Fixnum.ONE.subtract(arg));
417        LispObject result = n1.subtract(n2);
418        result = result.divideBy(Fixnum.TWO);
419
420        return complexToRealFixup(result, arg);
421    }
422
423    // ### cis
424    private static final Primitive CIS = new Primitive("cis", "radians")
425    {
426        @Override
427        public LispObject execute(LispObject arg)
428        {
429            return cis(arg);
430        }
431    };
432
433    static LispObject cis(LispObject arg)
434    {
435        if (arg.realp())
436            return Complex.getInstance(cos(arg), sin(arg));
437        return type_error(arg, Symbol.REAL);
438    }
439
440    // ### exp
441    private static final Primitive EXP = new Primitive("exp", "number")
442    {
443        @Override
444        public LispObject execute(LispObject arg)
445        {
446            return exp(arg);
447        }
448    };
449
450    static LispObject exp(LispObject arg)
451    {
452        if (arg.realp()) {
453            if (arg instanceof DoubleFloat) {
454                double d = Math.pow(Math.E, ((DoubleFloat)arg).value);
455                return OverUnderFlowCheck(new DoubleFloat(d));
456            } else {
457                float f = (float) Math.pow(Math.E, SingleFloat.coerceToFloat(arg).value);
458                return OverUnderFlowCheck(new SingleFloat(f));
459            }
460        }
461        if (arg instanceof Complex) {
462            Complex c = (Complex) arg;
463            return exp(c.getRealPart()).multiplyBy(cis(c.getImaginaryPart()));
464        }
465        return type_error(arg, Symbol.NUMBER);
466    }
467
468    // ### sqrt
469    private static final Primitive SQRT = new Primitive("sqrt", "number")
470    {
471        @Override
472        public LispObject execute(LispObject arg)
473        {
474            return sqrt(arg);
475        }
476    };
477
478    static final LispObject sqrt(LispObject obj)
479    {
480        if (obj instanceof DoubleFloat) {
481            if (obj.minusp())
482                return Complex.getInstance(new DoubleFloat(0), sqrt(obj.negate()));
483            return new DoubleFloat(Math.sqrt(DoubleFloat.coerceToFloat(obj).value));
484        }
485        if (obj.realp()) {
486            if (obj.minusp())
487                return Complex.getInstance(new SingleFloat(0), sqrt(obj.negate()));
488            return new SingleFloat((float)Math.sqrt(SingleFloat.coerceToFloat(obj).value));
489        }
490        if (obj instanceof Complex) {
491            LispObject imagpart = ((Complex)obj).imagpart;
492            if (imagpart.zerop()) {
493                LispObject realpart = ((Complex)obj).realpart;
494                if (realpart.minusp())
495                    return Complex.getInstance(imagpart, sqrt(realpart.negate()));
496                else
497                    return Complex.getInstance(sqrt(realpart), imagpart);
498            }
499            return exp(log(obj).divideBy(Fixnum.TWO));
500        }
501        return type_error(obj, Symbol.NUMBER);
502    }
503
504    // ### log
505    private static final Primitive LOG =
506        new Primitive("log", "number &optional base")
507    {
508        @Override
509        public LispObject execute(LispObject arg)
510        {
511            return log(arg);
512        }
513        @Override
514        public LispObject execute(LispObject number, LispObject base)
515
516        {
517            if (number.realp() && !number.minusp()
518                && base.isEqualTo(Fixnum.getInstance(10))) {
519                double d =
520                    Math.log10(DoubleFloat.coerceToFloat(number).value);
521                if (number instanceof DoubleFloat
522                    || base instanceof DoubleFloat)
523                    return new DoubleFloat(d);
524                else
525                    return new SingleFloat((float)d);
526            }
527            return log(number).divideBy(log(base));
528        }
529    };
530
531    static final LispObject log(LispObject obj)
532    {
533        if (obj.realp() && !obj.minusp()) {
534            // Result is real.
535            if (obj instanceof Fixnum)
536                return new SingleFloat((float)Math.log(((Fixnum)obj).value));
537            if (obj instanceof Bignum)
538                return new SingleFloat((float)Math.log(((Bignum)obj).doubleValue()));
539            if (obj instanceof Ratio)
540                return new SingleFloat((float)Math.log(((Ratio)obj).doubleValue()));
541            if (obj instanceof SingleFloat)
542                return new SingleFloat((float)Math.log(((SingleFloat)obj).value));
543            if (obj instanceof DoubleFloat)
544                return new DoubleFloat(Math.log(((DoubleFloat)obj).value));
545        } else {
546            // Result is complex.
547            if (obj.realp() && obj.minusp()) {
548                if (obj instanceof DoubleFloat) {
549                    DoubleFloat re = DoubleFloat.coerceToFloat(obj);
550                    DoubleFloat abs = new DoubleFloat(Math.abs(re.value));
551                    DoubleFloat phase = new DoubleFloat(Math.PI);
552                    return Complex.getInstance(new DoubleFloat(Math.log(abs.getValue())), phase);
553                } else {
554                    SingleFloat re = SingleFloat.coerceToFloat(obj);
555                    SingleFloat abs = new SingleFloat(Math.abs(re.value));
556                    SingleFloat phase = new SingleFloat((float)Math.PI);
557                    return Complex.getInstance(new SingleFloat((float)Math.log(abs.value)), phase);
558                }
559            } else if (obj instanceof Complex) {
560                if (((Complex)obj).getRealPart() instanceof DoubleFloat) {
561                    DoubleFloat re = DoubleFloat.coerceToFloat(((Complex)obj).getRealPart());
562                    DoubleFloat im = DoubleFloat.coerceToFloat(((Complex)obj).getImaginaryPart());
563                    DoubleFloat phase =
564                        new DoubleFloat(Math.atan2(im.getValue(), re.getValue()));  // atan(y/x)
565                    DoubleFloat abs = DoubleFloat.coerceToFloat(obj.ABS());
566                    return Complex.getInstance(new DoubleFloat(Math.log(abs.getValue())), phase);
567                } else {
568                    SingleFloat re = SingleFloat.coerceToFloat(((Complex)obj).getRealPart());
569                    SingleFloat im = SingleFloat.coerceToFloat(((Complex)obj).getImaginaryPart());
570                    SingleFloat phase =
571                        new SingleFloat((float)Math.atan2(im.value, re.value));  // atan(y/x)
572                    SingleFloat abs = SingleFloat.coerceToFloat(obj.ABS());
573                    return Complex.getInstance(new SingleFloat((float)Math.log(abs.value)), phase);
574                }
575            }
576        }
577        type_error(obj, Symbol.NUMBER);
578        return NIL;
579    }
580
581    // ### expt base-number power-number => result
582    public static final Primitive EXPT =
583        new Primitive("expt", "base-number power-number")
584    {
585        @Override
586        public LispObject execute(LispObject base, LispObject power)
587
588        {
589            if (power.zerop()) {
590                if (power instanceof Fixnum) {
591                    if (base instanceof SingleFloat)
592                        return SingleFloat.ONE;
593                    if (base instanceof DoubleFloat)
594                        return DoubleFloat.ONE;
595                    if (base instanceof Complex) {
596                        if (((Complex)base).realpart instanceof SingleFloat)
597                            return Complex.getInstance(SingleFloat.ONE,
598                                                       SingleFloat.ZERO);
599                        if (((Complex)base).realpart instanceof DoubleFloat)
600                            return Complex.getInstance(DoubleFloat.ONE,
601                                                       DoubleFloat.ZERO);
602                    }
603                    return Fixnum.ONE;
604                }
605                if (power instanceof DoubleFloat)
606                    return DoubleFloat.ONE;
607                if (base instanceof DoubleFloat)
608                    return DoubleFloat.ONE;
609                return SingleFloat.ONE;
610            }
611            if (base.zerop())
612                return base;
613            if (base.isEqualTo(1))
614                return base;
615           
616            if ((power instanceof Fixnum
617                 || power instanceof Bignum)
618                 && (base.rationalp()
619                     || (base instanceof Complex
620                         && ((Complex)base).realpart.rationalp()))) {
621                // exact math version
622                return intexp(base, power);
623            }
624            // for anything not a rational or complex rational, use
625            // float approximation.
626            boolean wantDoubleFloat = false;
627            if (base instanceof DoubleFloat)
628                wantDoubleFloat = true;
629            else if (power instanceof DoubleFloat)
630                wantDoubleFloat = true;
631            else if (base instanceof Complex
632                     && (((Complex)base).getRealPart() instanceof DoubleFloat
633                         || ((Complex)base).getImaginaryPart() instanceof DoubleFloat))
634                wantDoubleFloat = true;
635            else if (power instanceof Complex
636                    && (((Complex)power).getRealPart() instanceof DoubleFloat
637                         || ((Complex)power).getImaginaryPart() instanceof DoubleFloat))
638                wantDoubleFloat = true;
639
640            if (wantDoubleFloat) {
641                if (power instanceof Complex)
642                    power = ((Complex)power).coerceToDoubleFloat();
643                else
644                    power = DoubleFloat.coerceToFloat(power);
645
646                if (base instanceof Complex)
647                    base = ((Complex)base).coerceToDoubleFloat();
648                else
649                    base = DoubleFloat.coerceToFloat(base);
650            }
651
652
653
654            if (base instanceof Complex || power instanceof Complex)
655                return exp(power.multiplyBy(log(base)));
656            final double x; // base
657            final double y; // power
658            if (base instanceof Fixnum)
659                x = ((Fixnum)base).value;
660            else if (base instanceof Bignum)
661                x = ((Bignum)base).doubleValue();
662            else if (base instanceof Ratio)
663                x = ((Ratio)base).doubleValue();
664            else if (base instanceof SingleFloat)
665                x = ((SingleFloat)base).value;
666            else if (base instanceof DoubleFloat)
667                x = ((DoubleFloat)base).value;
668            else
669                return error(new LispError("EXPT: unsupported case: base is of type " +
670                                            base.typeOf().writeToString()));
671
672            if (power instanceof Fixnum)
673                y = ((Fixnum)power).value;
674            else if (power instanceof Bignum)
675                y = ((Bignum)power).doubleValue();
676            else if (power instanceof Ratio)
677                y = ((Ratio)power).doubleValue();
678            else if (power instanceof SingleFloat)
679                y = ((SingleFloat)power).value;
680            else if (power instanceof DoubleFloat)
681                y = ((DoubleFloat)power).value;
682            else
683                return error(new LispError("EXPT: unsupported case: power is of type " +
684                                            power.typeOf().writeToString()));
685            double r = Math.pow(x, y);
686            if (Double.isNaN(r)) {
687                if (x < 0) {
688                    r = Math.pow(-x, y);
689                    double realPart = r * Math.cos(y * Math.PI);
690                    double imagPart = r * Math.sin(y * Math.PI);
691                    if (base instanceof DoubleFloat || power instanceof DoubleFloat)
692                        return Complex
693                            .getInstance(OverUnderFlowCheck(new DoubleFloat(realPart)),
694                                         OverUnderFlowCheck(new DoubleFloat(imagPart)));
695                    else
696                        return Complex
697                            .getInstance(OverUnderFlowCheck(new SingleFloat((float)realPart)),
698                                         OverUnderFlowCheck(new SingleFloat((float)imagPart)));
699                }
700            }
701            if (base instanceof DoubleFloat || power instanceof DoubleFloat)
702                return OverUnderFlowCheck(new DoubleFloat(r));
703            else
704                return OverUnderFlowCheck(new SingleFloat((float)r));
705        }
706    };
707
708    /** Checks number for over- or underflow values.
709     *
710     * @param number
711     * @return number or signals an appropriate error
712     */
713    final static LispObject OverUnderFlowCheck(LispObject number)
714
715    {
716        if (number instanceof Complex) {
717            OverUnderFlowCheck(((Complex)number).realpart);
718            OverUnderFlowCheck(((Complex)number).imagpart);
719            return number;
720        }
721
722        if (TRAP_OVERFLOW) {
723            if (number instanceof SingleFloat)
724                if (Float.isInfinite(((SingleFloat)number).value))
725                    return error(new FloatingPointOverflow(NIL));
726            if (number instanceof DoubleFloat)
727                if (Double.isInfinite(((DoubleFloat)number).value))
728                    return error(new FloatingPointOverflow(NIL));
729        }
730        if (TRAP_UNDERFLOW) {
731            if (number.zerop())
732                return error(new FloatingPointUnderflow(NIL));
733        }
734        return number;
735    }
736
737    /** Checks number for over- or underflow values.
738     *
739     * @param number
740     * @return number or signals an appropriate error
741     */
742    final static float OverUnderFlowCheck(float number)
743
744    {
745        if (TRAP_OVERFLOW) {
746            if (Float.isInfinite(number))
747                error(new FloatingPointOverflow(NIL));
748        }
749        if (TRAP_UNDERFLOW) {
750            if (number == 0)
751                error(new FloatingPointUnderflow(NIL));
752        }
753        return number;
754    }
755
756    /** Checks number for over- or underflow values.
757     *
758     * @param number
759     * @return number or signals an appropriate error
760     */
761    public final static double OverUnderFlowCheck(double number)
762
763    {
764        if (TRAP_OVERFLOW) {
765            if (Double.isInfinite(number))
766                error(new FloatingPointOverflow(NIL));
767        }
768        if (TRAP_UNDERFLOW) {
769            if (number == 0)
770                error(new FloatingPointUnderflow(NIL));
771        }
772        return number;
773    }
774    // Adapted from SBCL.
775    /** Return the exponent of base taken to the integer exponent power
776     *
777     * @param base A value of any type
778     * @param power An integer (fixnum or bignum) value
779     */
780    static final LispObject intexp(LispObject base, LispObject power)
781
782    {
783        if (power.isEqualTo(0))
784            return Fixnum.ONE;
785        if (base.isEqualTo(1))
786            return base;
787        if (base.isEqualTo(0))
788            return base;
789
790        if (power.minusp()) {
791            power = Fixnum.ZERO.subtract(power);
792            return Fixnum.ONE.divideBy(intexp(base, power));
793        }
794        if (base.eql(Fixnum.TWO))
795            return Fixnum.ONE.ash(power);
796
797        LispObject nextn = power.ash(Fixnum.MINUS_ONE);
798        LispObject total;
799        if (power.oddp())
800            total = base;
801        else
802            total = Fixnum.ONE;
803        while (true) {
804            if (nextn.zerop())
805                return total;
806            base = base.multiplyBy(base);
807
808            if (nextn.oddp())
809                total = base.multiplyBy(total);
810            nextn = nextn.ash(Fixnum.MINUS_ONE);
811        }
812    }
813}
Note: See TracBrowser for help on using the repository browser.