package edu.northwestern.at.utils.math;

/*  Please see the license information in the header below. */

/** Basic numeric operations not included in standard Java run-time library.
 *
 *  <p>
 *  The binomial methods are modified from those in the Colt library.
 *  </p>
 *
 *  <p>
 *  The following methods for trigonometric functions come from the
 *  Sfun class written by Visual Numerics.
 *  </p>
 *
 *  <ul>
 *  <li>acosh</li>
 *  <li>asinh</li>
 *  <li>atanh</li>
 *  <li>cot</li>
 *  <li>cosh</li>
 *  <li>sinh</li>
 *  <li>tanh</li>
 *  </ul>
 *
 *  <p>
 *  These methods are covered by the following license.
 *  </p>
 *
 * -------------------------------------------------------------------------
 *  $Id: Sfun.java,v 1.1.1.1 1999/03/05 21:43:39 brophy Exp $
 * -------------------------------------------------------------------------
 * Copyright (c) 1997 - 1998 by Visual Numerics, Inc. All rights reserved.
 *
 * Permission to use, copy, modify, and distribute this software is freely
 * granted by Visual Numerics, Inc., provided that the copyright notice
 * above and the following warranty disclaimer are preserved in human
 * readable form.
 *
 * Because this software is licensed free of charge, it is provided
 * "AS IS", with NO WARRANTY.  TO THE EXTENT PERMITTED BY LAW, VNI
 * DISCLAIMS ALL WARRANTIES, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED
 * TO ITS PERFORMANCE, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 * VNI WILL NOT BE LIABLE FOR ANY DAMAGES WHATSOEVER ARISING OUT OF THE USE
 * OF OR INABILITY TO USE THIS SOFTWARE, INCLUDING BUT NOT LIMITED TO DIRECT,
 * INDIRECT, SPECIAL, CONSEQUENTIAL, PUNITIVE, AND EXEMPLARY DAMAGES, EVEN
 * IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGES.
 *
 * -------------------------------------------------------------------------
 */

public class ArithUtils
{
    /** The smallest relative spacing for doubles.*/

    public final static double EPSILON_SMALL = Constants.MACHEPS / 2.0D;

    /** The largest relative spacing for doubles. */

    public final static double EPSILON_LARGE = Constants.MACHEPS;

    // Series on [0,0.0625]

    private static final double COT_COEF[] =
    {
        .240259160982956302509553617744970e+0,
        -.165330316015002278454746025255758e-1,
        -.429983919317240189356476228239895e-4,
        -.159283223327541046023490851122445e-6,
        -.619109313512934872588620579343187e-9,
        -.243019741507264604331702590579575e-11,
        -.956093675880008098427062083100000e-14,
        -.376353798194580580416291539706666e-16,
        -.148166574646746578852176794666666e-18
    };

    // Series on the interval [0,1]

    private static final double SINH_COEF[] =
    {
        0.1730421940471796,
        0.08759422192276048,
        0.00107947777456713,
        0.00000637484926075,
        0.00000002202366404,
        0.00000000004987940,
        0.00000000000007973,
        0.00000000000000009
    };

    // Series on [0,1]

    private static final double TANH_COEF[] =
    {
        -.25828756643634710,
        -.11836106330053497,
        .009869442648006398,
        -.000835798662344582,
        .000070904321198943,
        -.000006016424318120,
        .000000510524190800,
        -.000000043320729077,
        .000000003675999055,
        -.000000000311928496,
        .000000000026468828,
        -.000000000002246023,
        .000000000000190587,
        -.000000000000016172,
        .000000000000001372,
        -.000000000000000116,
        .000000000000000009
    };

    // Series on the interval [0,1]

    private static final double ASINH_COEF[] =
    {
         -.12820039911738186343372127359268e+0,
        -.58811761189951767565211757138362e-1,
        .47274654322124815640725249756029e-2,
        -.49383631626536172101360174790273e-3,
        .58506207058557412287494835259321e-4,
        -.74669983289313681354755069217188e-5,
        .10011693583558199265966192015812e-5,
        -.13903543858708333608616472258886e-6,
        .19823169483172793547317360237148e-7,
        -.28847468417848843612747272800317e-8,
        .42672965467159937953457514995907e-9,
        -.63976084654366357868752632309681e-10,
        .96991686089064704147878293131179e-11,
        -.14844276972043770830246658365696e-11,
        .22903737939027447988040184378983e-12,
        -.35588395132732645159978942651310e-13,
        .55639694080056789953374539088554e-14,
        -.87462509599624678045666593520162e-15,
        .13815248844526692155868802298129e-15,
        -.21916688282900363984955142264149e-16,
        .34904658524827565638313923706880e-17
    };

    // Series on the interval [0,0.25]

    private static final double ATANH_COEF[] =
    {
        .9439510239319549230842892218633e-1,
        .4919843705578615947200034576668e-1,
        .2102593522455432763479327331752e-2,
        .1073554449776116584640731045276e-3,
        .5978267249293031478642787517872e-5,
        .3505062030889134845966834886200e-6,
        .2126374343765340350896219314431e-7,
        .1321694535715527192129801723055e-8,
        .8365875501178070364623604052959e-10,
        .5370503749311002163881434587772e-11,
        .3486659470157107922971245784290e-12,
        .2284549509603433015524024119722e-13,
        .1508407105944793044874229067558e-14,
        .1002418816804109126136995722837e-15,
        .6698674738165069539715526882986e-17,
        .4497954546494931083083327624533e-18
    };

    /** Return the inverse (arc) hyperbolic cosine of a double.
     *
     *  @param  x   Double whose inverse hyperbolic cosine is desired.
     *
     *  @return     The inverse hyperbolic cosine of x.
     *
     *  <p>
     *  If x is NaN or less than one, the result is NaN.
     *  </p>
     *
     *  <p>
     *  This method is a modified version of the one in the
     *  Visual Numerics Sfun class.
     *  </p>
     */

    public static double acosh( double x )
    {
        double ans;

        if ( Double.isNaN( x ) || ( x < 1 ) )
        {
            ans = Double.NaN;
        }
                                // 94906265.62 = 1.0/Math.sqrt(EPSILON_SMALL)

        else if ( x < 94906265.62 )
        {
            ans = safeLog( x + Math.sqrt( x * x - 1.0D ) );
        }
        else
        {
            ans = 0.69314718055994530941723212145818D + safeLog( x );
        }

        return ans;
    }

    /** Check if two doubles are equal to machine precision.
     *
     *  @param  a   First double.
     *  @param  b   Second double.
     *
     *  @return     True if a and b are equal to machine precision.
     */

    public static boolean areEqual( double a , double b )
    {
        return areEqual( a , b , Constants.MACHEPS );
    }

    /** Check if two doubles are equal to specified tolerance.
     *
     *  @param  a           First double.
     *  @param  b           Second double.
     *  @param  tolerance   Tolerance.
     *
     *  @return             True if a and b are equal to specified tolerance.
     */

    public static boolean areEqual( double a , double b , double tolerance )
    {
        boolean result  = ( a == b );

        if ( !result )
        {
            if ( a == 0.0D )
            {
                result  = ( Math.abs( b ) <= tolerance );
            }
            else if ( b == 0.0D )
            {
                result  = ( Math.abs( a ) <= tolerance );
            }
            else
            {
                result  = ( fuzzyCompare( a , b , tolerance ) == 0 );
            }
        }

        return result;
    }

    /** Return the inverse (arc) hyperbolic sine of a double.
     *
     *  @param  x   The value whose inverse hyperbolic sine is desired.
     *
     *  @return     The inverse hyperbolic sine of x.
     *
     *  <p>
     *  If x is NaN, the result is NaN.
     *  </p>
     *
     *  <p>
     *  This method is a modified version of the one in the
     *  Visual Numerics Sfun class.
     *  </p>
     */

    public static double asinh( double x )
    {
        double  ans;

        double  y   = Math.abs( x );

        if ( Double.isNaN( x ) )
        {
            ans = Double.NaN;
        }
                                // 1.05367e-08 = Math.sqrt(EPSILON_SMALL)

        else if ( y <= 1.05367e-08 )
        {
            ans = x;
        }
        else if ( y <= 1.0D )
        {
            ans =
                x *
                ( 1.0D +
                    Polynomial.evaluateChebyschev(
                        ASINH_COEF , 2.0D * x * x - 1.0D ) );
        }
                                // 94906265.62 = 1/Math.sqrt(EPSILON_SMALL)

        else if ( y < 94906265.62D )
        {
            ans = safeLog( y + Math.sqrt( y * y + 1.0D ) );
        }
        else
        {
            ans = 0.69314718055994530941723212145818D + safeLog( y );
        }

        if ( x < 0.0D ) ans = -ans;

        return ans;
    }

    /** Returns the inverse (arc) hyperbolic tangent of a double.
     *
     *  @param  x   The value whose inverse hyperbolic tangent is desired.
     *
     *  @return      The arc hyperbolic tangent of x.
     *
     *  <p>
     *  If x is NaN or |x|>1, the result is NaN.
     *  </p>
     *
     *  <p>
     *  This method is a modified version of the one in the
     *  Visual Numerics Sfun class.
     *  </p>
     */

    public static double atanh( double x )
    {
        double  ans;

        double  y   = Math.abs( x );

        if ( Double.isNaN( x ) )
        {
            ans = Double.NaN;
        }
                                // 1.82501e-08 = Math.sqrt(3.0*EPSILON_SMALL)

        else if ( y < 1.82501e-08 )
        {
            ans = x;
        }
        else if ( y <= 0.5D )
        {
            ans =
                x *
                ( 1.0D +
                    Polynomial.evaluateChebyschev(
                        ATANH_COEF , 8.0D * x * x - 1.0D ) );
        }
        else if ( y < 1.0D )
        {
            ans = 0.5D * safeLog( ( 1.0D + x ) / ( 1.0D - x ) );
        }
        else if ( y == 1.0D )
        {
            ans = x * Double.POSITIVE_INFINITY;
        }
        else
        {
            ans = Double.NaN;
        }

        return ans;
    }

    /** Efficiently returns the binomial coefficient, often also referred to as "n over k" or "n choose k".
     *
     *  <p>
     *  The binomial coefficient is defined as <tt>(n * n-1 * ... * n-k+1 ) / ( 1 * 2 * ... * k )</tt>.
     *  </p>
     *
     *  <ul>
     *  <li>k<0<tt>: <tt>0</tt>.</li>
     *  <li>k==0<tt>: <tt>1</tt>.</li>
     *  <li>k==1<tt>: <tt>n</tt>.</li>
     *  <li>else: <tt>(n * n-1 * ... * n-k+1 ) / ( 1 * 2 * ... * k )</tt>.</li>
     *  </ul>
     *
     *  @return     The binomial coefficient.
     */

    public static double binomial( double n , long k )
    {
        if ( k < 0 ) return 0.0D;
        if ( k == 0 ) return 1.0D;
        if ( k == 1 ) return (double)n;

                                // Compute binomial(n,k) =
                                // (n * n-1 * ... * n-k+1 ) / ( 1 * 2 * ... * k )

        double a        = n - k + 1;
        double b        = 1.0D;
        double binomial = 1.0D;

        for ( long i = k; i-- > 0; )
        {
            binomial *= (a++) / (b++);
        }

        return binomial;
    }

    /** Efficiently returns the binomial coefficient, often also referred to as "n over k" or "n choose k".
     *
     *  <p>
     *  The binomial coefficient is defined as
     *  </p>
     *  <ul>
     *  <li>k<0<tt>: <tt>0</tt>.</li>
     *  <li>k==0 || k==n<tt>: <tt>1</tt>.</li>
     *  <li>k==1 || k==n-1<tt>: <tt>n</tt>.</li>
     *  <li>else: <tt>(n * n-1 * ... * n-k+1 ) / ( 1 * 2 * ... * k )</tt>.</li>
     *  </ul>
     *
     *  @return         The binomial coefficient.
     */

    public static double binomial( long n , long k )
    {
        if ( k < 0 ) return 0.0D;
        if ( ( k == 0 ) || ( k == n ) ) return 1.0D;
        if ( ( k == 1 ) || ( k == ( n - 1 ) ) ) return n;

        // try quick version and see whether we get numeric overflows.
        // factorial(..) is O(1); requires no loop; only a table lookup.
        if (n > k)
        {
            int max = Factorial.longFactorials.length +
                Factorial.doubleFactorials.length;

            if (n < max)
            { // if (n! < inf && k! < inf)
                double n_fac = Factorial.factorial((int) n);
                double k_fac = Factorial.factorial((int) k);
                double n_minus_k_fac = Factorial.factorial((int) (n - k));
                double nk = n_minus_k_fac * k_fac;

                if (nk != Double.POSITIVE_INFINITY)
                { // no numeric overflow?
                    // now this is completely safe and accurate
                    return n_fac / nk;
                }
            }
            if ( k > ( n / 2 ) ) k = n - k; // quicker
        }

        // binomial(n,k) = (n * n-1 * ... * n-k+1 ) / ( 1 * 2 * ... * k )

        long a          = n - k + 1;
        long b          = 1;
        double binomial = 1.0D;

        for ( long i = k ; i-- > 0; )
        {
            binomial    *= ( (double)( a++)) / (b++);
        }

        return binomial;
    }

    /** Return the hyperbolic cosine of a double.
     *
     *  @param  x   The value whose hyperbolic cosine is desired.
     *
     *  @return     The hyperbolic cosine of x.
     *
     *  <p>
     *  If x is NaN, the result is NaN.
     *  </p>
     *
     *  <p>
     *  This method is a modified version of the one in the
     *  Visual Numerics Sfun class.
     *  </p>
     */

    static public double cosh( double x )
    {
        double  ans;

        double  y   = Math.exp( Math.abs( x ) );

        if ( Double.isNaN( x ) )
        {
            ans = Double.NaN;
        }
        else if ( Double.isInfinite( x ) )
        {
            ans = x;
        }
                                // 94906265.62 = 1.0/Math.sqrt(EPSILON_SMALL)

        else if ( y < 94906265.62D )
        {
            ans = 0.5D * ( y + 1.0D / y );
        }
        else
        {
            ans = 0.5D * y;
        }

        return ans;
    }

    /** Return the contangent of a double.
     *
     *  @param  x   The number whose cotangent is desired.
     *
     *  @return     The cotangent.
     *
     *  <p>
     *  This method is a modified version of the one in the
     *  Visual Numerics Sfun class.
     *  </p>
     */

    public static double cot( double x )
    {
        double ans, ainty, ainty2, prodbg, y, yrem;
        double pi2rec = 0.011619772367581343075535053490057; //  2/PI - 0.625

        y   = Math.abs( x );

                                // 4.5036e+15 = 1.0/EPSILON_LARGE

        if ( y > 4.5036e+15 )
        {
            return Double.NaN;
        }

        // Carefully compute
        // Y * (2/PI) = (AINT(Y) + REM(Y)) * (.625 + PI2REC)
        //      = AINT(.625*Y) + REM(.625*Y) + Y*PI2REC  =  AINT(.625*Y) + Z
        //      = AINT(.625*Y) + AINT(Z) + REM(Z)

        ainty  = (int)y;
        yrem   = y - ainty;
        prodbg = 0.625D * ainty;
        ainty  = (int)prodbg;
        y      = ( prodbg - ainty ) + 0.625D * yrem + y * pi2rec;
        ainty2 = (int)y;
        ainty  = ainty + ainty2;
        y      = y - ainty2;

        int ifn = (int)( ainty % 2.0 );

        if ( ifn == 1 ) y   = 1.0D - y;

        if ( y == 0.0D )
        {
            ans = Double.POSITIVE_INFINITY;
        }
                                // 1.82501e-08 = Math.sqrt(3.0*EPSILON_SMALL)

        else if ( y <= 1.82501e-08 )
        {
            ans = 1.0D / y;
        }

        else if ( y <= 0.25D )
        {
            ans =
                ( 0.5D +
                    Polynomial.evaluateChebyschev(
                        COT_COEF , 32.0D * y * y - 1.0D ) ) / y;
        }

        else if ( y <= 0.5D )
        {
            ans =
                (   0.5D +
                    Polynomial.evaluateChebyschev(
                        COT_COEF , 8.0D * y * y - 1.0D ) ) / ( 0.5D * y );

            ans = ( ans * ans - 1.0D ) * 0.5D / ans;
        }

        else
        {
            ans =
                (   0.5D +
                    Polynomial.evaluateChebyschev(
                        COT_COEF , 2.0D * y * y - 1.0D ) ) / ( 0.25D * y );

//$$$PIB$$$ Is one of the following two lines bogus?

            ans = ( ans * ans - 1.0D ) * 0.5D / ans;
            ans = ( ans * ans - 1.0D ) * 0.5D / ans;
        }

        if ( x != 0.0D ) ans = sign( ans , x );

        if ( ifn == 1 ) ans = -ans;

        return ans;
    }

    /** Get exp( x ) - 1.
     *
     *  @param  x   The number for which to find exp( x ) - 1.
     *
     *  @return     exp( x ) - 1.
     *
     *  <p>
     *  Example:    expm1( 9.995003330835334E-4 ) is 0.001 .
     *  </p>
     *
     *  <p>
     *  Implements a method suggested by William Kahan.
     *  </p>
     */

    public static double expm1( double x )
    {
        double result;

        double u    = Math.exp( x );

        if ( u == 1.0D )
        {
            result  = x;
        }

        else if ( ( u - 1.0D )  == -1.0D )
        {
            result  = -1.0D;
        }

        else
        {
            result  = ( u - 1.0D ) * ( x / Math.log( u ) );
        }

        return result;
    }

    /** Perform fuzzy comparison of two doubles with specified tolerance.
     *
     *  @param  a           First double.
     *  @param  b           Second double.
     *  @param  tolerance   Tolerance value.
     *
     *  @return             1 if a > b.
     *                      0 if a ~= b.
     *                      -1 if a < b.
     *
     *  <p>
     *  This is an implementation of an algorithm suggested by
     *  Donald E. Knuth in Section 4.2.2 of
     *  <em>Seminumerical Algorithms (3rd edition)</em>.
     *  </p>
     */

    public static int fuzzyCompare( double a , double b , double tolerance )
    {
                                //  Check for exacty equality first
                                //  to handle NaNs.

        if ( a == b ) return 0;

                                //  Compute difference of a and b.

        double difference   = a - b;

                                //  Find exponent for whichever of a or b
                                //  has largest absolute value.

        double maxAbs       = Math.max( Math.abs( a ) , Math.abs( b ) );

                                //  Get exponent.  Round up to next
                                //  power of two.

        int exponent        = new SplitDouble( maxAbs ).exponent;

                                //  Form neighborhood of size  2 * delta.

        double delta        = tolerance * Math.pow( 2 , exponent );

                                //  Assume a and b at least approximately
                                //  equal.

        int result          = 0;

                                //  Return 1 if a > b and difference
                                //  is outside delta neighborhood.

        if ( difference > delta )
        {
            result  = 1;
        }
                                //  Return -1 if a < b and difference is
                                //  outside delta neighborhood.

        else if ( difference < -delta )
        {
            result  = -1;
        }
                                //  If difference lies between
                                //  -delta and delta, a is exactly
                                //  or approximately equal to b .
                                //  Return 0 in this case.
        return result;
    }

    /** Safely calculate hypotenuse value.
     *
     *  @param  a   One leg of triangle.
     *  @param  b   Second leg of triangle.
     *
     *  @return     Hypotenuse value.
     *
     *  <p>
     *  The hypotenuse value is given mathematically as
     *  the sqrt( a^2 + b^2 ).  The method implemented
     *  here reduces the chances of cancellation and
     *  roundoff error.  If the |a| > |b|, we compute
     *  the hypotenuse as:
     *  </p>
     *
     *  <p>
     *  hypotenuse  = |a| * sqrt( 1 + (b/a) * (b/a) )
     *  </p>
     *
     *  <p>
     *  Otherwise b != 0 compute the hypotenuse as:
     *  </p>
     *
     *  <p>
     *  hypotenuse  = |b| * sqrt( 1 + (a/b) * (a/b) )
     *  </p>
     *
     *  <p>
     *  If b is zero, the hypotenuse is zero.
     *  </p>
     */

    public static double hypot( double a , double b )
    {
        double r;

        if ( Math.abs( a ) > Math.abs( b ) )
        {
            r   = b / a;
            r   = Math.abs( a ) * Math.sqrt( 1.0D + r * r );
        }
        else if ( b != 0 )
        {
            r   = a / b;
            r   = Math.abs( b ) * Math.sqrt( 1.0D + r * r );
        }
        else
        {
            r   = 0.0D;
        }

        return r;
    }

    /** Check if number is negative zero.
     *
     *  @param  x   The number to check.
     *
     *  @return     true if x is negative zero, false otherwise.
     */

    public static boolean isNegativeZero( double x )
    {
        return ( ( x == 0.0D ) && ( ( 1.0D / x ) < 0.0D ) );
    }

    /** Get log base 2 of a double.
     *
     *  @param  x   The number whose log base 2 value is desired.
     *
     *  @return     The log base 2 of x.
     *              If x is <= 0, 0 is returned.
     *
     *  <p>
     *  Example:    log2( 32 ) is 5.0D .
     *  </p>
     */

    public static double log2( double x )
    {
        double result   = 0.0D;

        if ( x > 0.0D )
        {
            result  = safeLog( x ) * Constants.LN2INV;
        }

        return result;
    }

    /** Get log base 10 of a double.
     *
     *  @param  x   The number whose log base 10 value is desired.
     *
     *  @return     The log base 10 of x.
     *              If x is <= 0, 0 is returned.
     *
     *  <p>
     *  Example:    log10( 100.0D ) is 2.0D .
     *  </p>
     */

    public static double log10( double x )
    {
        double result   = 0.0D;

        if ( x > 0.0D )
        {
            result  = safeLog( x ) * Constants.LN10INV;
        }

        return result;
    }

    /** Get log of a double + 1.
     *
     *  @param  x   The number for which we want log( x + 1 ).
     *
     *  @return     The log of ( x + 1 ).
     *
     *  <p>
     *  Example:    log1p( 0.001D ) is 9.995003330835334E-4
     *  </p>
     *
     *  <p>
     *  Implements a method suggested by William Kahan.
     *  </p>
     */

    public static double log1p( double x )
    {
        double result;
        double u;
                                //  Use log(), corrected to first order
                                //  for truncation loss.
        u   = 1.0D + x;

        if ( u == 1.0D )
        {
            result  = x;
        }
        else
        {
            result  = ( Math.log( u ) * ( x / ( u - 1.0D ) ) );
        }

        return result;
    }

    /** Get the log( exp( logX ) + exp( logY ) ).
     *
     *  @param  logX    Log( x )
     *  @param  logY    Log( y )
     *
     *  @return         log( x + y ) .
     */

    public static double logSumLogs( double logX , double logY )
    {
        return Math.max( logX , logY ) +
            log1p( Math.exp( -Math.abs( logX - logY ) ) );
    }

    /** Get the log( exp( logX ) - exp( logY ) ).
     *
     *  @param  logX    Log( x )
     *  @param  logY    Log( y )
     *
     *  @return         log( x - y ) .
     */

    public static double logDiffLogs( double logX , double logY )
    {
        return logX + log1p( -Math.exp( logY - logX ) );
    }

    /** Round double to specified number of decimal places.
     *
     *  @param  x   The double to round.
     *  @param  n   The number of decimal places to round to.
     *
     *  @return     x rounded to n decimal places.
     */

    public static double round( double x , int n )
    {
        double  pow10   = Math.pow( 10 , n );

        return Math.round( x * pow10 ) / pow10;
    }

    /** Return natural log of a double.
     *
     *  @param  x   The number whose natural log is desired.
     *
     *  @return     The natural log of x.  If x is zero,
     *              returns 0.
     */

    public static double safeLog( double x )
    {
        if ( x == 0.0D )
        {
            return 0.0D;
        }
        else
        {
            return Math.log( x );
        }
    }

    /** Return sign of an integer.
     *
     *  @param  n   Number whose sign is desired.
     *
     *  @return     -1 if n < 0, 0 if n isn 0, 1 if n > 0.
     */

    public static int sign( int n )
    {
        if ( n > 0 )
        {
            return 1;
        }
        else if ( n < 0 )
        {
            return -1;
        }

        return 0;
    }

    /** Return sign of a double.
     *
     *  @param  d   double whose sign is desired.
     *
     *  @return     -1 if d < 0, 0 if d is 0, 1 if d > 0.
     */

    public static int sign( double d )
    {
        if ( d > 0.0D )
        {
            return 1;
        }
        else if ( d < 0.0D )
        {
            return -1;
        }

        return 0;
    }

    /** Return the value of a double with the sign of another double.
     *
     *  @param  x   First double.
     *  @param  y   Second double.
     *
     *  @return     x with the sign of y.
     */

    public static double sign( double x , double y )
    {
        double abs_x    = ( ( x < 0 ) ? -x : x );

        return ( y < 0.0 ) ? -abs_x : abs_x;
    }

    /** Compute hyperbolic sine of a double.
     *
     *  @param  x   The number whose hyperbolic sine is desired.
     *
     *  @return     The hyperbolic sine of x.
     *
     *  <p>
     *  This method is a modified version of the one in the
     *  Visual Numerics Sfun class.
     *  </p>
     */

    public static double sinh( double x )
    {
        double  ans;

        double  y   = Math.abs( x );

        if ( Double.isNaN( x ) )
        {
            ans = Double.NaN;
        }

        else if ( Double.isInfinite( y ) )
        {
            return x;
        }
                                // 2.58096e-08 = Math.sqrt( 6.0 * EPSILON_SMAL L)

        else if ( y < 2.58096e-08 )
        {
            ans = x;
        }

        else if ( y <= 1.0D )
        {
            ans =
                x *
                ( 1.0D +
                    Polynomial.evaluateChebyschev(
                        SINH_COEF , 2.0D * x * x - 1.0D ) );
        }

        else
        {
            y   = Math.exp( y );

                                // 94906265.62 = 1.0/Math.sqrt(EPSILON_SMALL)

            if ( y >= 94906265.62D )
            {
                ans = sign( 0.5D * y , x );
            }
            else
            {
                ans = sign( 0.5D * ( y - 1.0D / y ) , x );
            }
        }

        return ans;
    }

    /** Return the integer portion of a double as a double.
     *
     *  @param  x   The double whose integer portion is to be found.
     *
     *  @return     The integer portion of x.
     *
     *  <p>
     *  Example:    trunc( 30.12345D ) is 30.0D .
     *  </p>
     */

    public static double trunc( double x )
    {
        long lx = (long)x;
        return (double)lx;
    }


    /** Return the hyperbolic tangent of a double.
     *
     *  @param  x   The value whose hyperbolic tangent is desired.
     *
     *  @return     The hyperbolic tangent of x.
     *
     *  <p>
     *  This method is a modified version of the one in the
     *  Visual Numerics Sfun class.
     *  </p>
     */

    public static double tanh( double x )
    {
        double  ans, y;

        y   = Math.abs( x );

        if ( Double.isNaN( x ) )
        {
            ans = Double.NaN;
        }
                                // 1.82501e-08 = Math.sqrt(3.0*EPSILON_SMALL)

        else if ( y < 1.82501e-08 )
        {
            ans = x;
        }
        else if ( y <= 1.0D )
        {
            ans =
                x *
                ( 1.0D +
                    Polynomial.evaluateChebyschev(
                        TANH_COEF , 2.0D * x * x - 1.0D ) );
        }
                                // 7.977294885 = -0.5*Math.log(EPSILON_SMALL)

        else if ( y < 7.977294885 )
        {
            y   = Math.exp( y );
            ans = sign( ( y - 1.0D / y ) / ( y + 1.0D / y ) , x );
        }
        else
        {
            ans = sign( 1.0D , x );
        }

        return ans;
    }

    /* Don't allow instantiation but allow overrides. */

    protected ArithUtils()
    {
    }
}

