/************************************************************************
 *                                                                      *
 * A replacement log() routine for PML library.  It really computes     *
 * base 2 logarithm which can be multiplied by a proper constant        *
 * to provide final answer.  A rational approximation of an original    *
 * is replaced by a polynomial one.  In practice, with a software       *
 * floating point, this gives a better precision.                       *
 *                                                                      *
 * Michal Jaegermann, December 1990                                     *
 *                                                                      *
 * If GCC_HACK is defined then we are folding log and log10 routines    *
 * by making in assembler an extra entry point.  Do not define that     *
 * for portable routine!!                                               *
 *                                                                      *
 ************************************************************************
 */
/* #define GCC_HACK */
/************************************************************************
 *									*
 *				N O T I C E				*
 *									*
 *			Copyright Abandoned, 1987, Fred Fish		*
 *									*
 *	This previously copyrighted work has been placed into the	*
 *	public domain by the author (Fred Fish) and may be freely used	*
 *	for any purpose, private or commercial.  I would appreciate	*
 *	it, as a courtesy, if this notice is left in all copies and	*
 *	derivative works.  Thank you, and enjoy...			*
 *									*
 *	The author makes no warranty of any kind with respect to this	*
 *	product and explicitly disclaims any implied warranties of	*
 *	merchantability or fitness for any particular purpose.		*
 *									*
 ************************************************************************
 */


/*
 *  FUNCTION
 *
 *	log   double precision natural log
 *
 *  KEY WORDS
 *
 *	log
 *	machine independent routines
 *	math libraries
 *
 *  DESCRIPTION
 *
 *	Returns double precision natural log of double precision
 *	floating point argument.
 *
 *  USAGE
 *
 *	double log (x)
 *	double x;
 *
 *  REFERENCES
 *
 *	Computer Approximations, J.F. Hart et al, John Wiley & Sons,
 *	1968, pp. 105-111.
 *
 *  RESTRICTIONS
 *
 *	The absolute error in the approximating rational function is
 *	10**(-19.38).  Note that contrary to DEC's assertion
 *	in the F4P user's guide, the error is absolute and not
 *	relative.
 *      ** modified ** theoretical precision is only 10**(-16.5);
 *      it works better in practice. 
 *
 *	This error bound assumes exact arithmetic
 *	in the polynomial evaluation.  Additional rounding and
 *	truncation errors may occur as the argument is reduced
 *	to the range over which the polynomial approximation
 *	is valid, and as the polynomial is evaluated using
 *	finite-precision arithmetic.
 *	
 *  PROGRAMMER
 *
 *	Fred Fish
 *      Modifications - Michal Jaegermann
 *
 *  INTERNALS
 *
 *	Computes log(X) from:
 *
 *	  (1)	If argument is zero then flag an error
 *		and return minus infinity (or rather its
 *		machine representation).
 *
 *	  (2)	If argument is negative then flag an
 *		error and return minus infinity.
 *
 *	  (3)	Given that x = m * 2**k then extract
 *		mantissa m and exponent k.
 *        (3a)  If m = 0.5 then we know what is the final
 *              result and calculations can be shortened. 
 *
 *	  (4)	Transform m which is in range [0.5,1.0]
 *		to range [1/sqrt(2), sqrt(2)] by:
 *
 *			s = m * sqrt(2)
 *
 *	  (5)	Compute z = (s - 1) / (s + 1)
 *        (5a)  Modified - combine steps (4) and (5) computing
 *              z = (m - 1/sqrt(2)) / (m + 1/sqrt(2))
 *
 *	  (6)	Now use the approximation from HART
 *		page 111 to find log(s):
 *
 *		log(s) = z * ( P(z**2) / Q(z**2) )
 *
 *		Where:
 *
 *		P(z**2) =  SUM [ Pj * (z**(2*j)) ]
 *		over j = {0,1,2,3}
 *
 *		Q(z**2) =  SUM [ Qj * (z**(2*j)) ]
 *		over j = {0,1,2,3}
 *
 *		P0 =  -0.240139179559210509868484e2
 *		P1 =  0.30957292821537650062264e2
 *		P2 =  -0.96376909336868659324e1
 *		P3 =  0.4210873712179797145
 *		Q0 =  -0.120069589779605254717525e2
 *		Q1 =  0.19480966070088973051623e2
 *		Q2 =  -0.89111090279378312337e1
 *		Q3 =  1.0000
 *
 *		(coefficients from HART table #2705 pg 229)
  *	  (6a)	Modified - compute, as a polynomial of z, an
 *              approximation of log2(m * sqrt(2)). Coefficients
 *              for this polynomial are not given explicitely in HART
 *              but can be computed from table #2665, for example.
 *
 *	  (7)	Finally, compute log(x) from:
 *
 *		log(x) = (k * log(2)) - log(sqrt(2)) + log(s)
 *
 *        (7a)  log2(x) = k - 1/2 + log2(m * sqrt(2)).  Multiply
 *              by a constant to adjust a logarithm base.
 *
 */

#include <stdio.h>
#include <pmluser.h>
#include "pml.h"

#define HALFSQRT2		0.70710678118654752440

static double log_coeffs[] = {
	2.88539008177793058026,
	0.9617966939212138784,
	0.577078018095541030,
	0.4121983030496934,
	0.32062199711811,
	0.2612917312170,
	0.24451102108
};

#ifdef GCC_HACK
static double log_of2[] = {
	LN2, 0.30102999566398119521 /* log10(2) */
};

static char funcname[] = "log";
#else
static char funcname[] = "log2";
#endif



#ifdef GCC_HACK
/*
 * This fake function header may change even from version to a version
 * of gcc compiler.  Ensure that first four assembler instructions in
 * log and log10 are the same!
 */
__asm__("
	.text
	.even
	.globl _log10
_log10:");
#ifdef __MSHORT__
__asm__("
	link a6,#-24
	moveml #0x3e30,sp@-");
#else
__asm__("
	link a6,#-28
	moveml #0x3f30,sp@-");
#endif  /* __MSHORT__ */
__asm__("
	movel a6@(8),d2
	movel a6@(12),d3
	moveq #1,d6
	jra   lgentry");

double log (x)
#else
double log2 (x)
#endif  /* GCC_HACK */

double x;
{
    auto int k;
    auto double s;
    auto struct exception xcpt;
    auto char * xcptstr;
#ifdef GCC_HACK
    auto int index;
#endif
    extern double frexp ();
    extern double poly ();


    DBUG_ENTER (funcname);
    DBUG_3 ("login", "arg %le", x);
#ifdef GCC_HACK
    index = 0;

__asm__("
lgentry:");
#endif

    if (x <= 0.0) {
	xcpt.retval = -MAXDOUBLE;
	xcpt.name = funcname;
	xcpt.arg1 = x;
    	if (x == 0.0) {
	    xcpt.type = SING;
	    xcptstr = "SINGULARITY";
	}
	else {
	    xcpt.type = DOMAIN;
	    xcptstr = "DOMAIN";
	}
	if (!matherr (&xcpt)) {
	    fprintf (stderr, "%s: %s error\n", xcptstr, funcname);
	    errno = EDOM;
	}
    }
    else {
	s = frexp (x, &k);
	if (0.5 == s ) {
		s = -1.0;
	}
	else {
	    /* range reduction with s multiplied by SQRT2 */
	    s = (s - HALFSQRT2) / (s + HALFSQRT2);
	    /* polynomial approximation */
	    s *= poly ((sizeof(log_coeffs)/sizeof(double)) - 1, 
	               log_coeffs, s * s);
	    /* and subtract log2 of SQRT2 */
	    s -= 0.5;
	}
	/* add the original binary exponent */
	s += k;
	/* multiply to get a requested logarithm */
#ifdef GCC_HACK
	xcpt.retval = s * log_of2[index];
#else
	xcpt.retval = s;
#endif
    }
    DBUG_3 ("logout", "result %le", xcpt.retval);
    DBUG_RETURN (xcpt.retval);
}

#ifndef GCC_HACK
double log (x)
double x;
{
    return (LN2 * log2(x));
}

double log10 (x)
double x;
{
    return (0.30102999566398119521 * log2(x));
}

#endif /* GCC_HACK */