source: trunk/abcl/src/org/armedbear/lisp/MathFunctions.java @ 11729

Last change on this file since 11729 was 11729, checked in by ehuelsmann, 16 years ago

Fix EXPT with a Bignum exponent.
At the same time, simplify EXPT and add some documentation.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 32.8 KB
Line 
1/*
2 * MathFunctions.java
3 *
4 * Copyright (C) 2004-2006 Peter Graves
5 * $Id: MathFunctions.java 11729 2009-04-04 21:57:58Z 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 java.lang.reflect.Method;
37
38public final class MathFunctions extends Lisp
39{
40    // ### sin
41    private static final Primitive SIN = new Primitive("sin", "radians")
42    {
43        @Override
44        public LispObject execute(LispObject arg) throws ConditionThrowable
45        {
46            return sin(arg);
47        }
48    };
49
50    private static LispObject sin(LispObject arg) throws ConditionThrowable
51    {
52        if (arg instanceof DoubleFloat)
53            return new DoubleFloat(Math.sin(((DoubleFloat)arg).value));
54        if (arg.realp())
55            return new SingleFloat((float)Math.sin(SingleFloat.coerceToFloat(arg).value));
56        if (arg instanceof Complex) {
57            LispObject n = arg.multiplyBy(Complex.getInstance(Fixnum.ZERO,
58                                                              Fixnum.ONE));
59            LispObject result = exp(n);
60            result = result.subtract(exp(n.multiplyBy(Fixnum.MINUS_ONE)));
61            return result.divideBy(Fixnum.TWO.multiplyBy(Complex.getInstance(Fixnum.ZERO,
62                                                                             Fixnum.ONE)));
63        }
64        return type_error(arg, Symbol.NUMBER);
65    }
66
67    // ### cos
68    private static final Primitive COS = new Primitive("cos", "radians")
69    {
70        @Override
71        public LispObject execute(LispObject arg) throws ConditionThrowable
72        {
73            return cos(arg);
74        }
75    };
76
77    private static LispObject cos(LispObject arg) throws ConditionThrowable
78    {
79        if (arg instanceof DoubleFloat)
80            return new DoubleFloat(Math.cos(((DoubleFloat)arg).value));
81        if (arg.realp())
82            return new SingleFloat((float)Math.cos(SingleFloat.coerceToFloat(arg).value));
83        if (arg instanceof Complex) {
84            LispObject n = arg.multiplyBy(Complex.getInstance(Fixnum.ZERO,
85                                                              Fixnum.ONE));
86            LispObject result = exp(n);
87            result = result.add(exp(n.multiplyBy(Fixnum.MINUS_ONE)));
88            return result.divideBy(Fixnum.TWO);
89        }
90        return type_error(arg, Symbol.NUMBER);
91    }
92
93    // ### tan
94    private static final Primitive TAN = new Primitive("tan", "radians")
95    {
96        @Override
97        public LispObject execute(LispObject arg) throws ConditionThrowable
98        {
99            if (arg instanceof DoubleFloat)
100                return new DoubleFloat(Math.tan(((DoubleFloat)arg).value));
101            if (arg.realp())
102                return new SingleFloat((float)Math.tan(SingleFloat.coerceToFloat(arg).value));
103            return sin(arg).divideBy(cos(arg));
104        }
105    };
106
107    // ### asin
108    private static final Primitive ASIN = new Primitive("asin", "number")
109    {
110        @Override
111        public LispObject execute(LispObject arg) throws ConditionThrowable
112        {
113            return asin(arg);
114        }
115    };
116
117    private static LispObject asin(LispObject arg) throws ConditionThrowable
118    {
119        if (arg instanceof SingleFloat) {
120            float f = ((SingleFloat)arg).value;
121            if (Math.abs(f) <= 1)
122                return new SingleFloat((float)Math.asin(f));
123        }
124        if (arg instanceof DoubleFloat) {
125            double d = ((DoubleFloat)arg).value;
126            if (Math.abs(d) <= 1)
127                return new DoubleFloat(Math.asin(d));
128        }
129        LispObject result = arg.multiplyBy(arg);
130        result = Fixnum.ONE.subtract(result);
131        result = sqrt(result);
132        LispObject n = Complex.getInstance(Fixnum.ZERO, Fixnum.ONE);
133        n = n.multiplyBy(arg);
134        result = n.add(result);
135        result = log(result);
136        result = result.multiplyBy(Complex.getInstance(Fixnum.ZERO,
137                                                       Fixnum.MINUS_ONE));
138        if (result instanceof Complex) {
139            if (arg instanceof Complex)
140                return result;
141            LispObject im = ((Complex)result).getImaginaryPart();
142            if (im.zerop())
143                return ((Complex)result).getRealPart();
144        }
145        return result;
146    }
147
148    // ### acos
149    private static final Primitive ACOS = new Primitive("acos", "number")
150    {
151        @Override
152        public LispObject execute(LispObject arg) throws ConditionThrowable
153        {
154            return acos(arg);
155        }
156    };
157
158    private static LispObject acos(LispObject arg) throws ConditionThrowable
159    {
160        if (arg instanceof DoubleFloat) {
161            double d = ((DoubleFloat)arg).value;
162            if (Math.abs(d) <= 1)
163                return new DoubleFloat(Math.acos(d));
164        }
165        if (arg instanceof SingleFloat) {
166            float f = ((SingleFloat)arg).value;
167            if (Math.abs(f) <= 1)
168                return new SingleFloat((float)Math.acos(f));
169        }
170        LispObject result = new DoubleFloat(Math.PI/2);
171        if (!(arg instanceof DoubleFloat)) {
172            if (arg instanceof Complex &&
173                    ((Complex)arg).getRealPart() instanceof DoubleFloat) {
174                    // do nothing; we want to keep the double float value
175            }
176            else
177                result = new SingleFloat((float)((DoubleFloat)result).value);
178        }
179        result = result.subtract(asin(arg));
180        if (result instanceof Complex) {
181            if (arg instanceof Complex)
182                return result;
183            LispObject im = ((Complex)result).getImaginaryPart();
184            if (im.zerop())
185                return ((Complex)result).getRealPart();
186        }
187        return result;
188    }
189
190    // ### atan
191    private static final Primitive ATAN =
192        new Primitive("atan", "number1 &optional number2")
193    {
194        @Override
195        public LispObject execute(LispObject arg) throws ConditionThrowable
196        {
197            if (arg.numberp())
198                return atan(arg);
199            return type_error(arg, Symbol.NUMBER);
200        }
201
202        // "If both number1 and number2 are supplied for atan, the result is
203        // the arc tangent of number1/number2."
204
205        // y = +0     x = +0       +0
206        // y = -0     x = +0       -0
207        // y = +0     x = -0       +<PI>
208        // y = -0     x = -0       -<PI>
209        @Override
210        public LispObject execute(LispObject y, LispObject x)
211            throws ConditionThrowable
212        {
213            if (!y.realp())
214                return type_error(y, Symbol.REAL);
215            if (!x.realp())
216                return type_error(x, Symbol.REAL);
217            double d1, d2;
218            d1 = DoubleFloat.coerceToFloat(y).value;
219            d2 = DoubleFloat.coerceToFloat(x).value;
220            double result = Math.atan2(d1, d2);
221            if (y instanceof DoubleFloat || x instanceof DoubleFloat)
222                return new DoubleFloat(result);
223            else
224                return new SingleFloat((float)result);
225        }
226    };
227
228    private static LispObject atan(LispObject arg) throws ConditionThrowable
229    {
230        if (arg instanceof Complex) {
231            LispObject im = ((Complex)arg).imagpart;
232            if (im.zerop())
233                return Complex.getInstance(atan(((Complex)arg).realpart),
234                                           im);
235            LispObject result = arg.multiplyBy(arg);
236            result = result.add(Fixnum.ONE);
237            result = Fixnum.ONE.divideBy(result);
238            result = sqrt(result);
239            LispObject n = Complex.getInstance(Fixnum.ZERO, Fixnum.ONE);
240            n = n.multiplyBy(arg);
241            n = n.add(Fixnum.ONE);
242            result = n.multiplyBy(result);
243            result = log(result);
244            result = result.multiplyBy(Complex.getInstance(Fixnum.ZERO, Fixnum.MINUS_ONE));
245            return result;
246        }
247        if (arg instanceof DoubleFloat)
248            return new DoubleFloat(Math.atan(((DoubleFloat)arg).value));
249        return new SingleFloat((float)Math.atan(SingleFloat.coerceToFloat(arg).value));
250    }
251
252    // ### sinh
253    private static final Primitive SINH = new Primitive("sinh", "number")
254    {
255        @Override
256        public LispObject execute(LispObject arg) throws ConditionThrowable
257        {
258            return sinh(arg);
259        }
260    };
261
262    private static Method sinhMethod = null;
263    static {
264        try {
265            sinhMethod = Class.forName("java.lang.Math")
266                    .getMethod("sinh", new Class[] { Double.TYPE });
267        }
268        catch (Throwable t) {
269            Debug.trace(t);
270        }
271    }
272
273    private static LispObject sinh(LispObject arg) throws ConditionThrowable
274    {
275        if (arg instanceof Complex) {
276            LispObject im = ((Complex)arg).getImaginaryPart();
277            if (im.zerop())
278                return Complex.getInstance(sinh(((Complex)arg).getRealPart()),
279                                           im);
280        }
281        if (arg instanceof SingleFloat) {
282            try {
283                if (sinhMethod != null) {
284                    Object[] args;
285                    args = new Object[1];
286                    args[0] = new Double(((SingleFloat)arg).value);
287                    Double d = (Double) sinhMethod.invoke(null, args);
288                    return new SingleFloat((float)d.doubleValue());
289                }
290            }
291            catch (Throwable t) {
292                Debug.trace(t);
293                // Fall through...
294            }
295        } else if (arg instanceof DoubleFloat) {
296            try {
297                if (sinhMethod != null) {
298                    Object[] args;
299                    args = new Object[1];
300                    args[0] = new Double(((DoubleFloat)arg).value);
301                    Double d = (Double) sinhMethod.invoke(null, args);
302                    return new DoubleFloat(d.doubleValue());
303                }
304            }
305            catch (Throwable t) {
306                Debug.trace(t);
307                // Fall through...
308            }
309        }
310        LispObject result = exp(arg);
311        result = result.subtract(exp(arg.multiplyBy(Fixnum.MINUS_ONE)));
312        result = result.divideBy(Fixnum.TWO);
313        if (result instanceof Complex) {
314            if (arg instanceof Complex)
315                return result;
316            LispObject im = ((Complex)result).getImaginaryPart();
317            if (im.zerop())
318                return ((Complex)result).getRealPart();
319        }
320        return result;
321    }
322
323    // ### cosh
324    private static final Primitive COSH = new Primitive("cosh", "number")
325    {
326        @Override
327        public LispObject execute(LispObject arg) throws ConditionThrowable
328        {
329            return cosh(arg);
330        }
331    };
332
333    private static Method coshMethod = null;
334    static {
335        try {
336            coshMethod = Class.forName("java.lang.Math")
337                    .getMethod("cosh", new Class[] { Double.TYPE });
338        }
339        catch (Throwable t) {
340            Debug.trace(t);
341        }
342    }
343
344    private static LispObject cosh(LispObject arg) throws ConditionThrowable
345    {
346        if (arg instanceof Complex) {
347            LispObject im = ((Complex)arg).getImaginaryPart();
348            if (im.zerop())
349                return Complex.getInstance(cosh(((Complex)arg).getRealPart()),
350                                           im);
351        }
352        if (arg instanceof SingleFloat) {
353            try {
354                if (coshMethod != null) {
355                    Object[] args;
356                    args = new Object[1];
357                    args[0] = new Double(((SingleFloat)arg).value);
358                    Double d = (Double) coshMethod.invoke(null, args);
359                    return new SingleFloat((float)d.doubleValue());
360                }
361            }
362            catch (Throwable t) {
363                Debug.trace(t);
364                // Fall through...
365            }
366        } else if (arg instanceof DoubleFloat) {
367            try {
368                if (coshMethod != null) {
369                    Object[] args;
370                    args = new Object[1];
371                    args[0] = new Double(((DoubleFloat)arg).value);
372                    Double d = (Double) coshMethod.invoke(null, args);
373                    return new DoubleFloat(d.doubleValue());
374                }
375            }
376            catch (Throwable t) {
377                Debug.trace(t);
378                // Fall through...
379            }
380        }
381        LispObject result = exp(arg);
382        result = result.add(exp(arg.multiplyBy(Fixnum.MINUS_ONE)));
383        result = result.divideBy(Fixnum.TWO);
384        if (result instanceof Complex) {
385            if (arg instanceof Complex)
386                return result;
387            LispObject im = ((Complex)result).getImaginaryPart();
388            if (im.zerop())
389                return ((Complex)result).getRealPart();
390        }
391        return result;
392    }
393
394    private static Method tanhMethod = null;
395    static {
396        try {
397            tanhMethod = Class.forName("java.lang.Math")
398                    .getMethod("tanh", new Class[] { Double.TYPE });
399        }
400        catch (Throwable t) {
401            Debug.trace(t);
402        }
403    }
404
405    // ### tanh
406    private static final Primitive TANH = new Primitive("tanh", "number")
407    {
408        @Override
409        public LispObject execute(LispObject arg) throws ConditionThrowable
410        {
411            if (arg instanceof SingleFloat) {
412                try {
413                    if (tanhMethod != null) {
414                        Object[] args;
415                        args = new Object[1];
416                        args[0] = new Double(((SingleFloat)arg).value);
417                        Double d = (Double) tanhMethod.invoke(null, args);
418                        return new SingleFloat((float)d.doubleValue());
419                    }
420                }
421                catch (Throwable t) {
422                    Debug.trace(t);
423                    // Fall through...
424                }
425            } else if (arg instanceof DoubleFloat) {
426                try {
427                    if (tanhMethod != null) {
428                        Object[] args;
429                        args = new Object[1];
430                        args[0] = new Double(((DoubleFloat)arg).value);
431                        Double d = (Double) tanhMethod.invoke(null, args);
432                        return new DoubleFloat(d.doubleValue());
433                    }
434                }
435                catch (Throwable t) {
436                    Debug.trace(t);
437                    // Fall through...
438                }
439            }
440            return sinh(arg).divideBy(cosh(arg));
441        }
442    };
443
444    // ### asinh
445    private static final Primitive ASINH = new Primitive("asinh", "number")
446    {
447        @Override
448        public LispObject execute(LispObject arg) throws ConditionThrowable
449        {
450            return asinh(arg);
451        }
452    };
453
454    private static LispObject asinh(LispObject arg) throws ConditionThrowable
455    {
456        if (arg instanceof Complex) {
457            LispObject im = ((Complex)arg).getImaginaryPart();
458            if (im.zerop())
459                return Complex.getInstance(asinh(((Complex)arg).getRealPart()),
460                                           im);
461        }
462        LispObject result = arg.multiplyBy(arg);
463        result = Fixnum.ONE.add(result);
464        result = sqrt(result);
465        result = result.add(arg);
466        result = log(result);
467        if (result instanceof Complex) {
468            if (arg instanceof Complex)
469                return result;
470            LispObject im = ((Complex)result).getImaginaryPart();
471            if (im.zerop())
472                return ((Complex)result).getRealPart();
473        }
474        return result;
475    }
476
477    // ### acosh
478    private static final Primitive ACOSH = new Primitive("acosh", "number")
479    {
480        @Override
481        public LispObject execute(LispObject arg) throws ConditionThrowable
482        {
483            return acosh(arg);
484        }
485    };
486
487    private static LispObject acosh(LispObject arg) throws ConditionThrowable
488    {
489        if (arg instanceof Complex) {
490            LispObject im = ((Complex)arg).getImaginaryPart();
491            if (im.zerop())
492                return Complex.getInstance(acosh(((Complex)arg).getRealPart()),
493                                           im);
494        }
495        LispObject n1 = arg.add(Fixnum.ONE);
496        n1 = n1.divideBy(Fixnum.TWO);
497        n1 = sqrt(n1);
498        LispObject n2 = arg.subtract(Fixnum.ONE);
499        n2 = n2.divideBy(Fixnum.TWO);
500        n2 = sqrt(n2);
501        LispObject result = n1.add(n2);
502        result = log(result);
503        result = result.multiplyBy(Fixnum.TWO);
504        if (result instanceof Complex) {
505            if (arg instanceof Complex)
506                return result;
507            LispObject im = ((Complex)result).getImaginaryPart();
508            if (im.zerop())
509                return ((Complex)result).getRealPart();
510        }
511        return result;
512    }
513
514    // ### atanh
515    private static final Primitive ATANH = new Primitive("atanh", "number")
516    {
517        @Override
518        public LispObject execute(LispObject arg) throws ConditionThrowable
519        {
520            return atanh(arg);
521        }
522    };
523
524    private static LispObject atanh(LispObject arg) throws ConditionThrowable
525    {
526        if (arg instanceof Complex) {
527            LispObject im = ((Complex)arg).getImaginaryPart();
528            if (im.zerop())
529                return Complex.getInstance(atanh(((Complex)arg).getRealPart()),
530                                           im);
531        }
532        LispObject n1 = log(Fixnum.ONE.add(arg));
533        LispObject n2 = log(Fixnum.ONE.subtract(arg));
534        LispObject result = n1.subtract(n2);
535        result = result.divideBy(Fixnum.TWO);
536        if (result instanceof Complex) {
537            if (arg instanceof Complex)
538                return result;
539            LispObject im = ((Complex)result).getImaginaryPart();
540            if (im.zerop())
541                return ((Complex)result).getRealPart();
542        }
543        return result;
544    }
545
546    // ### cis
547    private static final Primitive CIS = new Primitive("cis", "radians")
548    {
549        @Override
550        public LispObject execute(LispObject arg) throws ConditionThrowable
551        {
552            return cis(arg);
553        }
554    };
555
556    private static LispObject cis(LispObject arg) throws ConditionThrowable
557    {
558        if (arg.realp())
559            return Complex.getInstance(cos(arg), sin(arg));
560        return type_error(arg, Symbol.REAL);
561    }
562
563    // ### exp
564    private static final Primitive EXP = new Primitive("exp", "number")
565    {
566        @Override
567        public LispObject execute(LispObject arg) throws ConditionThrowable
568        {
569            return exp(arg);
570        }
571    };
572
573    private static LispObject exp(LispObject arg) throws ConditionThrowable
574    {
575        if (arg.realp()) {
576            if (arg instanceof DoubleFloat) {
577                double d = Math.pow(Math.E, ((DoubleFloat)arg).value);
578                if (TRAP_OVERFLOW && Double.isInfinite(d))
579                    return error(new FloatingPointOverflow(NIL));
580                if (d == 0 && TRAP_UNDERFLOW)
581                    return error(new FloatingPointUnderflow(NIL));
582                return new DoubleFloat(d);
583            } else {
584                float f = (float) Math.pow(Math.E, SingleFloat.coerceToFloat(arg).value);
585                if (TRAP_OVERFLOW && Float.isInfinite(f))
586                    return error(new FloatingPointOverflow(NIL));
587                if (f == 0 && TRAP_UNDERFLOW)
588                    return error(new FloatingPointUnderflow(NIL));
589                return new SingleFloat(f);
590            }
591        }
592        if (arg instanceof Complex) {
593            Complex c = (Complex) arg;
594            return exp(c.getRealPart()).multiplyBy(cis(c.getImaginaryPart()));
595        }
596        return type_error(arg, Symbol.NUMBER);
597    }
598
599    // ### sqrt
600    private static final Primitive SQRT = new Primitive("sqrt", "number")
601    {
602        @Override
603        public LispObject execute(LispObject arg) throws ConditionThrowable
604        {
605            return sqrt(arg);
606        }
607    };
608
609    private static final LispObject sqrt(LispObject obj) throws ConditionThrowable
610    {
611        if (obj instanceof DoubleFloat) {
612            if (obj.minusp())
613                return Complex.getInstance(new DoubleFloat(0), sqrt(obj.negate()));
614            return new DoubleFloat(Math.sqrt(DoubleFloat.coerceToFloat(obj).value));
615        }
616        if (obj.realp()) {
617            if (obj.minusp())
618                return Complex.getInstance(new SingleFloat(0), sqrt(obj.negate()));
619            return new SingleFloat((float)Math.sqrt(SingleFloat.coerceToFloat(obj).value));
620        }
621        if (obj instanceof Complex) {
622            LispObject imagpart = ((Complex)obj).imagpart;
623            if (imagpart.zerop()) {
624                LispObject realpart = ((Complex)obj).realpart;
625                if (realpart.minusp())
626                    return Complex.getInstance(imagpart, sqrt(realpart.negate()));
627                else
628                    return Complex.getInstance(sqrt(realpart), imagpart);
629            }
630            return exp(log(obj).divideBy(Fixnum.TWO));
631        }
632        return type_error(obj, Symbol.NUMBER);
633    }
634
635    // ### log
636    private static final Primitive LOG =
637        new Primitive("log", "number &optional base")
638    {
639        @Override
640        public LispObject execute(LispObject arg) throws ConditionThrowable
641        {
642            return log(arg);
643        }
644        @Override
645        public LispObject execute(LispObject number, LispObject base)
646            throws ConditionThrowable
647        {
648            if (number.realp() && !number.minusp()
649                && base.isEqualTo(Fixnum.getInstance(10))) {
650                try {
651                    double d =
652                        Math.log10(DoubleFloat.coerceToFloat(number).value);
653                    if (number instanceof DoubleFloat
654                        || base instanceof DoubleFloat)
655                        return new DoubleFloat(d);
656                    else
657                        return new SingleFloat((float)d);
658                }
659                catch (Throwable t) {
660                    Debug.trace(t);
661                    // Fall through...
662                }
663            }
664            return log(number).divideBy(log(base));
665        }
666    };
667
668    private static final LispObject log(LispObject obj) throws ConditionThrowable
669    {
670        if (obj.realp() && !obj.minusp()) {
671            // Result is real.
672            if (obj instanceof Fixnum)
673                return new SingleFloat((float)Math.log(((Fixnum)obj).value));
674            if (obj instanceof Bignum)
675                return new SingleFloat((float)Math.log(((Bignum)obj).doubleValue()));
676            if (obj instanceof Ratio)
677                return new SingleFloat((float)Math.log(((Ratio)obj).doubleValue()));
678            if (obj instanceof SingleFloat)
679                return new SingleFloat((float)Math.log(((SingleFloat)obj).value));
680            if (obj instanceof DoubleFloat)
681                return new DoubleFloat(Math.log(((DoubleFloat)obj).value));
682        } else {
683            // Result is complex.
684            if (obj.realp() && obj.minusp()) {
685                if (obj instanceof DoubleFloat) {
686                    DoubleFloat re = DoubleFloat.coerceToFloat(obj);
687                    DoubleFloat abs = new DoubleFloat(Math.abs(re.value));
688                    DoubleFloat phase = new DoubleFloat(Math.PI);
689                    return Complex.getInstance(new DoubleFloat(Math.log(abs.getValue())), phase);
690                } else {
691                    SingleFloat re = SingleFloat.coerceToFloat(obj);
692                    SingleFloat abs = new SingleFloat(Math.abs(re.value));
693                    SingleFloat phase = new SingleFloat((float)Math.PI);
694                    return Complex.getInstance(new SingleFloat((float)Math.log(abs.value)), phase);
695                }
696            } else if (obj instanceof Complex) {
697                if (((Complex)obj).getRealPart() instanceof DoubleFloat) {
698                    DoubleFloat re = DoubleFloat.coerceToFloat(((Complex)obj).getRealPart());
699                    DoubleFloat im = DoubleFloat.coerceToFloat(((Complex)obj).getImaginaryPart());
700                    DoubleFloat phase =
701                        new DoubleFloat(Math.atan2(im.getValue(), re.getValue()));  // atan(y/x)
702                    DoubleFloat abs = DoubleFloat.coerceToFloat(obj.ABS());
703                    return Complex.getInstance(new DoubleFloat(Math.log(abs.getValue())), phase);
704                } else {
705                    SingleFloat re = SingleFloat.coerceToFloat(((Complex)obj).getRealPart());
706                    SingleFloat im = SingleFloat.coerceToFloat(((Complex)obj).getImaginaryPart());
707                    SingleFloat phase =
708                        new SingleFloat((float)Math.atan2(im.value, re.value));  // atan(y/x)
709                    SingleFloat abs = SingleFloat.coerceToFloat(obj.ABS());
710                    return Complex.getInstance(new SingleFloat((float)Math.log(abs.value)), phase);
711                }
712            }
713        }
714        type_error(obj, Symbol.NUMBER);
715        return NIL;
716    }
717
718    // ### expt base-number power-number => result
719    public static final Primitive EXPT =
720        new Primitive("expt", "base-number power-number")
721    {
722        @Override
723        public LispObject execute(LispObject base, LispObject power)
724            throws ConditionThrowable
725        {
726            if (power.zerop()) {
727                if (power instanceof Fixnum) {
728                    if (base instanceof SingleFloat)
729                        return SingleFloat.ONE;
730                    if (base instanceof DoubleFloat)
731                        return DoubleFloat.ONE;
732                    if (base instanceof Complex) {
733                        if (((Complex)base).realpart instanceof SingleFloat)
734                            return Complex.getInstance(SingleFloat.ONE,
735                                                       SingleFloat.ZERO);
736                        if (((Complex)base).realpart instanceof DoubleFloat)
737                            return Complex.getInstance(DoubleFloat.ONE,
738                                                       DoubleFloat.ZERO);
739                    }
740                    return Fixnum.ONE;
741                }
742                if (power instanceof DoubleFloat)
743                    return DoubleFloat.ONE;
744                if (base instanceof DoubleFloat)
745                    return DoubleFloat.ONE;
746                return SingleFloat.ONE;
747            }
748            if (base.zerop())
749                return base;
750            if (base.isEqualTo(1))
751                return base;
752           
753            if ((power instanceof Fixnum
754                 || power instanceof Bignum)
755                 && (base.rationalp()
756                     || (base instanceof Complex
757                         && ((Complex)base).realpart.rationalp()))) {
758                // exact math version
759                return intexp(base, power);
760            }
761            // for anything not a rational or complex rational, use
762            // float approximation.
763            if (base instanceof Complex || power instanceof Complex)
764                return exp(power.multiplyBy(log(base)));
765            final double x; // base
766            final double y; // power
767            if (base instanceof Fixnum)
768                x = ((Fixnum)base).value;
769            else if (base instanceof Bignum)
770                x = ((Bignum)base).doubleValue();
771            else if (base instanceof Ratio)
772                x = ((Ratio)base).doubleValue();
773            else if (base instanceof SingleFloat)
774                x = ((SingleFloat)base).value;
775            else if (base instanceof DoubleFloat)
776                x = ((DoubleFloat)base).value;
777            else
778                return error(new LispError("EXPT: unsupported case: base is of type " +
779                                            base.typeOf().writeToString()));
780
781            if (power instanceof Fixnum)
782                y = ((Fixnum)power).value;
783            else if (power instanceof Bignum)
784                y = ((Bignum)power).doubleValue();
785            else if (power instanceof Ratio)
786                y = ((Ratio)power).doubleValue();
787            else if (power instanceof SingleFloat)
788                y = ((SingleFloat)power).value;
789            else if (power instanceof DoubleFloat)
790                y = ((DoubleFloat)power).value;
791            else
792                return error(new LispError("EXPT: unsupported case: power is of type " +
793                                            power.typeOf().writeToString()));
794            double r = Math.pow(x, y);
795            if (Double.isNaN(r)) {
796                if (x < 0) {
797                    r = Math.pow(-x, y);
798                    double realPart = r * Math.cos(y * Math.PI);
799                    double imagPart = r * Math.sin(y * Math.PI);
800                    if (base instanceof DoubleFloat || power instanceof DoubleFloat)
801                        return Complex
802                            .getInstance(OverUnderFlowCheck(new DoubleFloat(realPart)),
803                                         OverUnderFlowCheck(new DoubleFloat(imagPart)));
804                    else
805                        return Complex
806                            .getInstance(OverUnderFlowCheck(new SingleFloat((float)realPart)),
807                                         OverUnderFlowCheck(new SingleFloat((float)imagPart)));
808                }
809            }
810            if (base instanceof DoubleFloat || power instanceof DoubleFloat)
811                return OverUnderFlowCheck(new DoubleFloat(r));
812            else
813                return OverUnderFlowCheck(new SingleFloat((float)r));
814        }
815    };
816
817    /** Checks number for over- or underflow values.
818     *
819     * @param number
820     * @return number or signals an appropriate error
821     * @throws org.armedbear.lisp.ConditionThrowable
822     */
823    private final static LispObject OverUnderFlowCheck(LispObject number)
824            throws ConditionThrowable
825    {
826        if (number instanceof Complex) {
827            OverUnderFlowCheck(((Complex)number).realpart);
828            OverUnderFlowCheck(((Complex)number).imagpart);
829            return number;
830        }
831
832        if (TRAP_OVERFLOW) {
833            if (number instanceof SingleFloat)
834                if (Float.isInfinite(((SingleFloat)number).value))
835                    return error(new FloatingPointOverflow(NIL));
836            if (number instanceof DoubleFloat)
837                if (Double.isInfinite(((DoubleFloat)number).value))
838                    return error(new FloatingPointOverflow(NIL));
839        }
840        if (TRAP_UNDERFLOW) {
841            if (number.zerop())
842                return error(new FloatingPointUnderflow(NIL));
843        }
844        return number;
845    }
846
847    // Adapted from SBCL.
848    /** Return the exponent of base taken to the integer exponent power
849     *
850     * @param base A value of any type
851     * @param power An integer (fixnum or bignum) value
852     * @throws org.armedbear.lisp.ConditionThrowable
853     */
854    private static final LispObject intexp(LispObject base, LispObject power)
855        throws ConditionThrowable
856    {
857        if (power.isEqualTo(0))
858            return Fixnum.ONE;
859        if (base.isEqualTo(1))
860            return base;
861        if (base.isEqualTo(0))
862            return base;
863
864        if (power.minusp()) {
865            power = Fixnum.ZERO.subtract(power);
866            return Fixnum.ONE.divideBy(intexp(base, power));
867        }
868        if (base.eql(Fixnum.TWO))
869            return Fixnum.ONE.ash(power);
870
871        LispObject nextn = power.ash(Fixnum.MINUS_ONE);
872        LispObject total;
873        if (power.oddp())
874            total = base;
875        else
876            total = Fixnum.ONE;
877        while (true) {
878            if (nextn.zerop())
879                return total;
880            base = base.multiplyBy(base);
881
882            if (nextn.oddp())
883                total = base.multiplyBy(total);
884            nextn = nextn.ash(Fixnum.MINUS_ONE);
885        }
886    }
887}
Note: See TracBrowser for help on using the repository browser.