Artifact 7633662194bd2c73fd71e1942a547eea49d4ae74df02b378f86fe20faf6e60f2:

• Executable file r37/lisp/csl/jlisp/MyMath.java — part of check-in [f2fda60abd] at 2011-09-02 18:13:33 on branch master — Some historical releases purely for archival purposes

  git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 3832) [annotate] [blame] [check-ins using] [more...]

//
// This file is part of the Jlisp implementation of Standard Lisp
// Copyright \u00a9 (C) Codemist Ltd, 1998-2000.
//

import java.math.*;

// The implementations here are intended to avoid total stupidity about
// overflow and accuracy, but are not striving for the very best 
// least-significant-bit results.

class MyMath
{

static double acosh(double a)
{
    double a1 = a - 1.0;
    if (a1 == 0.0) return 1.0;
    else if (a1 < 0.0) return (0.0/0.0); // a NaN
// The series shown here was developed as a Chebychev approximation
// to the function acosh(x)/sqrt(x-1) near x=1. The cut-off at x-1=0.1
// is somewhat arbitrary, but at that point the simple formula used for
// large x only (?) loses a few bits when computing a*a-1 and the argument
// for the logarithm is 1.55 (far enough from 1.0 that I think that is
// not a danger).
    else if (a1 < 0.1)
    {   double w = 
           (((((((-0.00013247403176210292*a1+
           0.00038027261161096525)*a1-
           0.00098848023694688735)*a1+
           0.0026854001707801832)*a1-
           0.0078918165513171978)*a1+
           0.026516504292612476)*a1-
           0.1178511301977528)*a1+
           1.4142135623730950);
        return Math.sqrt(a1)*w;
    }
    else return Math.log(a + Math.sqrt(a*a - 1.0));
}

static double acoth(double a)
{
    return atanh(1.0/a);
}

static double acsch(double a)
{
    return asinh(1.0/a);
}

static double asech(double a)
{
    return acosh(1.0/a);
}

static double asinh(double a)
{
    if (a >= 0.0) return Math.log(a + Math.sqrt(1.0 + a*a));
    else return -Math.log(-a + Math.sqrt(1.0 + a*a));
}

static double atanh(double a)
{
    return 0.5*Math.log((1.0+a)/(1.0-a)); 
}

static double cosh(double a)
{
    if (a < 0.0) a = -a;
// no cancellation worries for small argument
    if (a < 20.0)
    {   double ea = Math.exp(a);
        return (ea + 1.0/ea)/2.0;
    }
// if exp(-a) is tiny compared with exp(a) I can ignore it
    else if (a < 700) return Math.exp(a)/2.0;
// avoid premature overflow in extreme cases
    else return Math.exp(a - Math.log(2.0));
}

static double coth(double a)
{
    return 1.0/tanh(a);
}

static double csch(double a)
{
    return 1.0/sinh(a);
}

static double sech(double a)
{
    return 1.0/cosh(a);
}

static double sinh(double a)
{
// for small arguments I use the series expansion to avoid cancellation
    double aa = Math.abs(a);
    if (aa < 0.35)
    {   double a2 = a*a;
        double r = (((((a2/110.0 + 1.0)*a2/72.0 + 1.0)*a2/42.0 +
	           1.0)*a2/20.0 + 1.0)*a2/6.0 + 1.0)*a;
	return r;  
    }
// for medium arguments the full formula can be used
    else if (aa < 20.0) return (Math.exp(a) - Math.exp(-a))/2.0;
// for |a| > 20 I can use a simplified version
    else if (aa < 700.0) aa = Math.exp(aa)/2.0;
// fially for very large args I must avoid premature overflow
    else aa = Math.exp(aa - Math.log(2.0));
    if (a < 0.0) return -aa;
    else return aa;
}

static double tanh(double a)
{
    double aa = Math.abs(a);
// for small argument I will first range reduce and then use a simple
// powert series.
    if (aa < 0.40)
    {   int n = 0;
        while (aa >= 0.05)
	{   aa = aa/2.0;
	    n++;
	}
	double a2 = aa*aa;
	double r = ((((62.0/2835.0*a2 - 17.0/315)*a2 + 2.0/15.0)*a2 -
	               1.0/3.0)*a2 + 1.0)*aa;
	while (n != 0)
	{   r = 2.0*r/(1 + r*r);
	    n--;
	}
	if (a < 0) return -r;
        else return r;
    }
// for large enough argument the value will be +1 or -1 to within accuracy
// limits
    else if (aa >= 20.0)
    {   if (a < 0.0) return -1.0;
        else return 1.0;
    }
// for intermediate ranges I can use the normal formula, secure that there
// will be no overflow, underflow or serious cancellation
    double ea = Math.exp(a);
    double ema = 1.0/ea;
    return (ea-ema)/(ea+ema);
}

}

// end of MyMath.java