jruby/core/src/main/java/org/jruby/RubyMath.java at master · jruby/jruby

739 lines (634 loc) · 29 KB
/***** BEGIN LICENSE BLOCK *****
 * Version: EPL 2.0/GPL 2.0/LGPL 2.1
 * The contents of this file are subject to the Eclipse Public
 * License Version 2.0 (the "License"); you may not use this file
 * except in compliance with the License. You may obtain a copy of
 * the License at http://www.eclipse.org/legal/epl-v20.html
 * Software distributed under the License is distributed on an "AS
 * IS" basis, WITHOUT WARRANTY OF ANY KIND, either express or
 * implied. See the License for the specific language governing
 * rights and limitations under the License.
 * Copyright (C) 2001 Chad Fowler <chadfowler@chadfowler.com>
 * Copyright (C) 2001-2002 Jan Arne Petersen <jpetersen@uni-bonn.de>
 * Copyright (C) 2002 Benoit Cerrina <b.cerrina@wanadoo.fr>
 * Copyright (C) 2002-2004 Thomas E Enebo <enebo@acm.org>
 * Copyright (C) 2004 Stefan Matthias Aust <sma@3plus4.de>
 * Alternatively, the contents of this file may be used under the terms of
 * either of the GNU General Public License Version 2 or later (the "GPL"),
 * or the GNU Lesser General Public License Version 2.1 or later (the "LGPL"),
 * in which case the provisions of the GPL or the LGPL are applicable instead
 * of those above. If you wish to allow use of your version of this file only
 * under the terms of either the GPL or the LGPL, and not to allow others to
 * use your version of this file under the terms of the EPL, indicate your
 * decision by deleting the provisions above and replace them with the notice
 * and other provisions required by the GPL or the LGPL. If you do not delete
 * the provisions above, a recipient may use your version of this file under
 * the terms of any one of the EPL, the GPL or the LGPL.
 ***** END LICENSE BLOCK *****/
package org.jruby;
import org.jruby.anno.JRubyMethod;
import org.jruby.anno.JRubyModule;
import org.jruby.runtime.ThreadContext;
import org.jruby.runtime.Visibility;
import org.jruby.runtime.builtin.IRubyObject;
import static org.jruby.api.Convert.asFixnum;
import static org.jruby.api.Convert.asFloat;
import static org.jruby.api.Convert.numToDouble;
import static org.jruby.api.Convert.toInt;
import static org.jruby.api.Create.newArray;
import static org.jruby.api.Define.defineModule;
@JRubyModule(name="Math")
public class RubyMath {
    /** Create the Math module and add it to the Ruby runtime.
    public static RubyModule createMathModule(ThreadContext context) {
        return defineModule(context, "Math").
                defineConstant(context, "E", asFloat(context, Math.E)).
                defineConstant(context, "PI", asFloat(context, Math.PI)).
                defineMethods(context, RubyMath.class);
    private static void domainCheck(ThreadContext context, double value, String msg) {
        if (Double.isNaN(value)) throw context.runtime.newMathDomainError(msg);
    public static double chebylevSerie(double x, double coef[]) {
        double  b0, b1, b2, twox;
        int i;
        b1 = 0.0;
        b0 = 0.0;
        b2 = 0.0;
        twox = 2.0 * x;
        for (i = coef.length-1; i >= 0; i--) {
            b0 = twox * b1 - b2 + coef[i];
        return 0.5*(b0 - b2);
    public static double sign(double x, double y) {
        double abs = ((x < 0) ? -x : x);
        return (y < 0.0) ? -abs : abs;
    @JRubyMethod(name = "atan2", module = true, visibility = Visibility.PRIVATE)
    public static RubyFloat atan2(ThreadContext context, IRubyObject recv, IRubyObject x, IRubyObject y) {
        return asFloat(context, Math.atan2(numToDouble(context, x), numToDouble(context, y)));
    @JRubyMethod(name = "cos", module = true, visibility = Visibility.PRIVATE)
    public static RubyFloat cos(ThreadContext context, IRubyObject recv, IRubyObject x) {
        return asFloat(context, Math.cos(numToDouble(context, x)));
    @JRubyMethod(name = "sin", module = true, visibility = Visibility.PRIVATE)
    public static RubyFloat sin(ThreadContext context, IRubyObject recv, IRubyObject x) {
        return asFloat(context, Math.sin(numToDouble(context, x)));
    @JRubyMethod(name = "tan", module = true, visibility = Visibility.PRIVATE)
    public static RubyFloat tan(ThreadContext context, IRubyObject recv,  IRubyObject x) {
        return asFloat(context, Math.tan(numToDouble(context, x)));
    @JRubyMethod(name = "asin", module = true, visibility = Visibility.PRIVATE)
    public static RubyFloat asin(ThreadContext context, IRubyObject recv, IRubyObject x) {
        double value = numToDouble(context, x);
        if (value < -1.0 || value > 1.0) throw context.runtime.newMathDomainError("asin");
        return asFloat(context, Math.asin(value));
    @JRubyMethod(name = "acos", module = true, visibility = Visibility.PRIVATE)
    public static RubyFloat acos(ThreadContext context, IRubyObject recv, IRubyObject x) {
        double value = numToDouble(context, x);
        if (value < -1.0 || value > 1.0) throw context.runtime.newMathDomainError("acos");
        return asFloat(context, Math.acos(value));
    @JRubyMethod(name = "atan", module = true, visibility = Visibility.PRIVATE)
    public static RubyFloat atan(ThreadContext context, IRubyObject recv, IRubyObject x) {
        return asFloat(context, Math.atan(numToDouble(context, x)));
    @JRubyMethod(name = "cosh", module = true, visibility = Visibility.PRIVATE)
    public static RubyFloat cosh(ThreadContext context, IRubyObject recv, IRubyObject x) {
        double value = numToDouble(context, x);
        return asFloat(context, (Math.exp(value) + Math.exp(-value)) / 2.0);
    @JRubyMethod(name = "sinh", module = true, visibility = Visibility.PRIVATE)
    public static RubyFloat sinh(ThreadContext context, IRubyObject recv, IRubyObject x) {
        double value = numToDouble(context, x);
        return asFloat(context, (Math.exp(value) - Math.exp(-value)) / 2.0);
    @JRubyMethod(name = "tanh", module = true, visibility = Visibility.PRIVATE)
    public static RubyFloat tanh(ThreadContext context, IRubyObject recv, IRubyObject x) {
        return asFloat(context, Math.tanh(numToDouble(context, x)));
    @JRubyMethod(name = "acosh", module = true, visibility = Visibility.PRIVATE)
    public static RubyFloat acosh(ThreadContext context, IRubyObject recv, IRubyObject x) {
        double value = numToDouble(context, x);
        double result;
        if (Double.isNaN(value)) {
            result = Double.NaN;
        } else if (value < 1) {
            throw context.runtime.newMathDomainError("acosh");
        } else if (value < 94906265.62) {
            result = Math.log(value + Math.sqrt(value * value - 1.0));
        } else{
            result = 0.69314718055994530941723212145818 + Math.log(value);
        return asFloat(context, result);
    public static final double[] ASINH_COEF = {
    @JRubyMethod(name = "asinh", module = true, visibility = Visibility.PRIVATE)
    public static RubyFloat asinh(ThreadContext context, IRubyObject recv, IRubyObject x) {
        double value = numToDouble(context, x);
        double  y = Math.abs(value);
        double result;
        if (Double.isNaN(value)) {
            result = Double.NaN;
        } else if (y <= 1.05367e-08) {
            result = value;
        } else if (y <= 1.0) {          
            result = value * (1.0 + chebylevSerie(2.0 * value * value - 1.0, ASINH_COEF));
        } else if (y < 94906265.62) {
            result = Math.log(value + Math.sqrt(value * value + 1.0));
        } else {    
            result = 0.69314718055994530941723212145818 + Math.log(y);
            if (value < 0) result *= -1;
        return asFloat(context, result);
    public static final double[] ATANH_COEF = {
    @JRubyMethod(name = "atanh", module = true, visibility = Visibility.PRIVATE)
    public static RubyFloat atanh(ThreadContext context, IRubyObject recv, IRubyObject x) {
        double value = numToDouble(context, x);
        if (value < -1.0 || value > 1.0) throw context.runtime.newMathDomainError("atanh");
        final double y = Math.abs(value);
        final double result;
        if (Double.isNaN(value)) {
            result = Double.NaN;
        } else if (y < 1.82501e-08) {
            result = value;
        } else if (y <= 0.5) {
            result = value * (1.0 + chebylevSerie(8.0 * value * value - 1.0, ATANH_COEF));
        } else if (y < 1.0) {
            result = 0.5 * Math.log((1.0 + value) / (1.0 - value));
        } else if (y == 1.0) {
            result = value * Double.POSITIVE_INFINITY;
        } else {
            result = Double.NaN;
        return asFloat(context, result);
    @JRubyMethod(name = "exp", module = true, visibility = Visibility.PRIVATE)
    public static RubyFloat exp(ThreadContext context, IRubyObject recv, IRubyObject exponent) {
        return exp(context, exponent);
    public static RubyFloat exp(ThreadContext context, IRubyObject exponent) {
        return asFloat(context, Math.exp(numToDouble(context, exponent)));
    // MRI : get_double_rshift
    private static double[] get_double_rshift(ThreadContext context, IRubyObject x) {
        int numbits = 0;
        if (x instanceof RubyBignum) {
            RubyBignum bignum = (RubyBignum)x;
            numbits = bignum.getValue().abs().bitLength();
            if (bignum.getValue().signum() > 0 && numbits >= DBL_MAX_EXP) {
                numbits -= DBL_MANT_DIG;
                x = bignum.op_rshift(context, numbits);
            } else {
                numbits = 0;
        return new double[]{numToDouble(context, x), numbits};
    private static int DBL_MANT_DIG = 53;
    private static int DBL_MAX_EXP = 1024;
    private static double LOG_E_2 = Math.log(2);
    private static double LOG_10_2 = Math.log10(2);
    /** Returns the natural logarithm of x.
    @JRubyMethod(name = "log", module = true, visibility = Visibility.PRIVATE)
    public static RubyFloat log(ThreadContext context, IRubyObject recv, IRubyObject val) {
        return log(context, val);
    public static RubyFloat log(ThreadContext context, IRubyObject val) {
        double[] ret = get_double_rshift(context, val);
        if (ret[0] < 0) throw context.runtime.newMathDomainError("log");
        /* log(d * 2 ** numbits) */
        return asFloat(context, Math.log(ret[0]) + ret[1] * LOG_E_2);
    @JRubyMethod(name = "log", module = true, visibility = Visibility.PRIVATE)
    public static RubyFloat log(ThreadContext context, IRubyObject recv, IRubyObject val, IRubyObject base) {
        double [] ret = get_double_rshift(context, val);
        double _base = numToDouble(context, base);
        if (ret[0] < 0 || _base < 0) throw context.runtime.newMathDomainError("log");
        /* log(d * 2 ** numbits) / log(base) */
        return asFloat(context, Math.log(ret[0]) / Math.log(_base) + ret[1]);
    public static RubyFloat log(ThreadContext context, IRubyObject recv, IRubyObject... args) {
        return args.length == 2 ? log(context, recv, args[0], args[1]) : log(context, recv, args[0]);
    /** Returns the base 10 logarithm of x.
    @JRubyMethod(name = "log10", module = true, visibility = Visibility.PRIVATE)
    public static RubyFloat log10(ThreadContext context, IRubyObject recv, IRubyObject x) {
        double [] ret = get_double_rshift(context, x);
        if (ret[0] < 0) throw context.runtime.newMathDomainError("log10");
        /* log10(d * 2 ** numbits) */
        return asFloat(context, Math.log10(ret[0]) + ret[1] * LOG_10_2);
    /** Returns the base 2 logarithm of x.
    @JRubyMethod(name = "log2", module = true, visibility = Visibility.PRIVATE)
    public static RubyFloat log2(ThreadContext context, IRubyObject recv, IRubyObject x) {
        double [] ret = get_double_rshift(context, x);
        if (ret[0] < 0) throw context.runtime.newMathDomainError("log2");
        /* log2(d * 2 ** numbits) */
        return asFloat(context, Math.log(ret[0]) / LOG_E_2 + ret[1]);
    @JRubyMethod(name = "sqrt", module = true, visibility = Visibility.PRIVATE)
    public static RubyFloat sqrt(ThreadContext context, IRubyObject recv, IRubyObject x) {
        double value = numToDouble(context, x);
        if (value < 0) throw context.runtime.newMathDomainError("sqrt");
        return asFloat(context, value == 0.0 ? 0.0 : Math.sqrt(value));
    @JRubyMethod(name = "cbrt", module = true, visibility = Visibility.PRIVATE)
    public static RubyFloat cbrt(ThreadContext context, IRubyObject recv, IRubyObject x) {
        return asFloat(context, Math.cbrt(numToDouble(context, x)));
    @Deprecated(since = "10.0.0.0")
    public static RubyFloat hypot19(ThreadContext context, IRubyObject recv, IRubyObject x, IRubyObject y) {
        return hypot(context, recv, x, y);
    @JRubyMethod(name = "hypot", module = true, visibility = Visibility.PRIVATE)
    public static RubyFloat hypot(ThreadContext context, IRubyObject recv, IRubyObject x, IRubyObject y) {
        double valuea = numToDouble(context, x);
        double valueb = numToDouble(context, y);
        double result;
        if (Math.abs(valuea) > Math.abs(valueb)) {
            result = valueb / valuea;
            result = Math.abs(valuea) * Math.sqrt(1 + result * result);
        } else if (valueb != 0) {
            result = valuea / valueb;
            result = Math.abs(valueb) * Math.sqrt(1 + result * result);
        } else if (Double.isNaN(valuea) || Double.isNaN(valueb)) {
            result = Double.NaN;
        } else {
            result = 0;
        return asFloat(context, result);
     * x = mantissa * 2 ** exponent
     * Where mantissa is in the range of [.5, 1)
    @JRubyMethod(name = "frexp", module = true, visibility = Visibility.PRIVATE)
    public static RubyArray frexp(ThreadContext context, IRubyObject recv, IRubyObject other) {
        double mantissa = numToDouble(context, other);
        short sign = 1;
        long exponent = 0;
        if (!Double.isInfinite(mantissa) && mantissa != 0.0) {
            // Make mantissa same sign so we only have one code path.
            if (mantissa < 0) {
                mantissa = -mantissa;
                sign = -1;
            // Increase value to hit lower range.
            for (; mantissa < 0.5; mantissa *= 2.0, exponent -=1) { }
            // Decrease value to hit upper range.  
            for (; mantissa >= 1.0; mantissa *= 0.5, exponent +=1) { }
        return newArray(context, asFloat(context, sign * mantissa), asFixnum(context, exponent));
     * r = x * 2 ** y
    @JRubyMethod(name = "ldexp", module = true, visibility = Visibility.PRIVATE)
    public static RubyFloat ldexp(ThreadContext context, IRubyObject recv, IRubyObject mantissa, IRubyObject exponent) {
        double m = numToDouble(context, mantissa);
        int e = toInt(context, exponent);
        return e > 1023 ? // avoid overflow. Math.power(2.0, 1024) is greater than Math.MAX_VALUE.
                asFloat(context, m * Math.pow(2.0, e - 1023) * Math.pow(2.0, 1023)) :
                asFloat(context, m * Math.pow(2.0, e));
    public static final double[] ERFC_COEF = {
    @JRubyMethod(name = "erf", module = true, visibility = Visibility.PRIVATE)
    public static RubyFloat erf(ThreadContext context, IRubyObject recv, IRubyObject x) {
        double value = numToDouble(context, x);
        double  y = Math.abs(value);
        double  result;
        if (y <= 1.49012e-08) {
            result = 2 * value / 1.77245385090551602729816748334;
        } else if (y <= 1) {
            result = value * (1 + chebylevSerie(2 * value * value - 1, ERFC_COEF));
        } else if (y < 6.013687357) {
            result = sign(1 - erfc(context, recv, asFloat(context, y)).value, value);
        } else if (Double.isNaN(y)) {
            result = Double.NaN;
        } else {
            result = sign(1, value);
        return asFloat(context,result);
    public static final double[] ERFC2_COEF = {
    public static final double[] ERFCC_COEF = {
    @JRubyMethod(name = "erfc", module = true, visibility = Visibility.PRIVATE)
    public static RubyFloat erfc(ThreadContext context, IRubyObject recv, IRubyObject x) {
        double value = numToDouble(context, x);
        double  result;
        double  y = Math.abs(value);
        if (value <= -6.013687357) {
            result = 2;
        } else if (y < 1.49012e-08) {
            result = 1 - 2 * value / 1.77245385090551602729816748334;
        } else {
            double ysq = y*y;
            if (y < 1) {
                result = 1 - value * (1 + chebylevSerie(2 * ysq - 1, ERFC_COEF));
            } else if (y <= 4.0) {
                result = Math.exp(-ysq)/y*(0.5+chebylevSerie((8.0 / ysq - 5.0) / 3.0, ERFC2_COEF));
                if (value < 0) result = 2.0 - result;
                if (value < 0) result = 2.0 - result;
                if (value < 0) result = 2.0 - result;
            } else {
                result = Math.exp(-ysq) / y * (0.5 + chebylevSerie(8.0 / ysq - 1, ERFCC_COEF));
                if (value < 0) result = 2.0 - result;
        return asFloat(context, result);
    private static final double FACTORIAL[] = {
        /*  0! */ 1.0,
        /*  1! */ 1.0,
        /*  2! */ 2.0,
        /*  3! */ 6.0,
        /*  4! */ 24.0,
        /*  5! */ 120.0,
        /*  6! */ 720.0,
        /*  7! */ 5040.0,
        /*  8! */ 40320.0,
        /*  9! */ 362880.0,
        /* 10! */ 3628800.0,
        /* 11! */ 39916800.0,
        /* 12! */ 479001600.0,
        /* 13! */ 6227020800.0,
        /* 14! */ 87178291200.0,
        /* 15! */ 1307674368000.0,
    private static final double NEMES_GAMMA_COEFF[] = {
     * Based on Gerg&#337; Nemes's Gamma Function approximation formula, we compute
     * approximate value of Gamma function of x.
     * @param recv Math module
     * @param x a real number
     * @return &Gamma;(x) for real number x
     * @see <a href="http://www.ebyte.it/library/downloads/2008_MTH_Nemes_GammaApproximationUpdate.pdf">
     * New asymptotic expansion for the &Gamma;(x) function</a>
    @JRubyMethod(name = "gamma", module = true, visibility = Visibility.PRIVATE)
    public static RubyFloat gamma(ThreadContext context, IRubyObject recv, IRubyObject x) {
        double value = RubyKernel.new_float(context, x).value;
        double result = nemes_gamma(value);
        /* note nemes_gamma can return Double.POSITIVE_INFINITY or Double.NEGATIVE_INFINITY
         * when value is an integer less than 1.
         * We treat 0 as a special case to avoid Domain error.
        if (Double.isInfinite(result)) {
            if (value < 0) {
                result = Double.NaN;
            } else {
                if (value == 0 && 1 / value < 0) {
                    result = Double.NEGATIVE_INFINITY;
                } else {
                    result = Double.POSITIVE_INFINITY;
        if (Double.isNaN(value)) return asFloat(context, Double.NaN);
        domainCheck(context, result, "gamma");
        return asFloat(context, result);
     * Based on Gerg&#337; Nemes's Gamma Function approximation formula, we compute
     * Log Gamma function for real number x.
     * @param recv Math module
     * @param x a real number
     * @return 2-element array [ln(&Gamma;(x)), sgn] for real number x,
     *  where sgn is the sign of &Gamma;(x) when exponentiated
     * @see #gamma(ThreadContext, org.jruby.runtime.builtin.IRubyObject, org.jruby.runtime.builtin.IRubyObject)
    @JRubyMethod(name = "lgamma", module = true, visibility = Visibility.PRIVATE)
    public static RubyArray lgamma(ThreadContext context, IRubyObject recv, IRubyObject x) {
        double value = RubyKernel.new_float(context, x).value;
        // JRUBY-4653: Could this error checking done more elegantly?
        if (value < 0 && Double.isInfinite(value)) throw context.runtime.newMathDomainError("lgamma");
        NemesLogGamma l = new NemesLogGamma(value);
        return newArray(context, asFloat(context, l.value), asFixnum(context, (int) l.sign));
    public static double nemes_gamma(double x) {
        double int_part = (int) x;
        if ((x - int_part) == 0.0 && 0 < int_part && int_part <= FACTORIAL.length) {
            return FACTORIAL[(int) int_part - 1];
        NemesLogGamma l = new NemesLogGamma(x);
        return l.sign * Math.exp(l.value);
     * Inner class to help with &Gamma; functions
    public static class NemesLogGamma {
        public final double value;
        public final double sign;
        public NemesLogGamma(double x) {
            if (Double.isInfinite(x)) {
                value = Double.POSITIVE_INFINITY; sign = 1;
                return;
            if (Double.isNaN(x)) {
                value = Double.NaN; sign = 1;
                return;
            double int_part = (int) x;
            sign = signum(x, int_part);
            if ((x - int_part) == 0.0 && 0 < int_part && int_part <= FACTORIAL.length) {
                value = Math.log(FACTORIAL[(int) int_part - 1]);
            else if (x < 10) {
                double rising_factorial = 1;
                for (int i = 0; i < (int) Math.abs(x) - int_part + 10; i++) {
                    rising_factorial *= (x + i);
                NemesLogGamma l = new NemesLogGamma(x + (int) Math.abs(x) - int_part + 10);
                value = l.value - Math.log(Math.abs(rising_factorial));
            } else {
                double temp = 0.0;
                for (int i = 0; i < NEMES_GAMMA_COEFF.length; i++) {
                    temp += NEMES_GAMMA_COEFF[i] * 1.0 / Math.pow(x, i);
                value = x * (Math.log(x) - 1 + Math.log(temp)) +
                        (Math.log(2) + Math.log(Math.PI) - Math.log(x)) / 2.0;
        private static int signum(final double x, final double int_part) {
            return ( (int_part % 2 == 0 && (x - int_part) != 0.0 && (x < 0)) || negZero(x) ) ? -1 : 1;
        private static boolean negZero(final double x) {
            return x == 0.0 && Double.doubleToRawLongBits(x) != 0; // detect -0.0 (since in Java: `0.0 == -0.0`)
    @JRubyMethod(module = true, visibility = Visibility.PRIVATE)
    public static IRubyObject expm1(ThreadContext context, IRubyObject unused, IRubyObject x) {
        return asFloat(context, Math.expm1(numToDouble(context, x)));
    @JRubyMethod(module = true, visibility = Visibility.PRIVATE)
    public static IRubyObject log1p(ThreadContext context, IRubyObject unused, IRubyObject x) {
        double xDouble = numToDouble(context, x);
        if (xDouble < -1.0) {
            throw context.runtime.newMathDomainError("log1p");
        return asFloat(context, Math.log1p(xDouble));
pFad - Phone/Frame/Anonymizer/Declutterfier! Saves Data!

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Pfad - The Proxy pFad © 2024 Your Company Name. All rights reserved.

pFad - Phone/Frame/Anonymizer/Declutterfier! Saves Data!

Uh oh!

FilesExpand file tree

RubyMath.java

Latest commit

History

RubyMath.java

File metadata and controls

Pfad - The Proxy pFad © 2024 Your Company Name. All rights reserved.
