/*************************************************************************
 *
 *  $RCSfile: solmath.hxx,v $
 *
 *  $Revision: 1.7.2.2 $
 *
 *  last change: $Author: armin $ $Date: 2002/10/29 16:59:55 $
 *
 *  The Contents of this file are made available subject to the terms of
 *  either of the following licenses
 *
 *         - GNU Lesser General Public License Version 2.1
 *         - Sun Industry Standards Source License Version 1.1
 *
 *  Sun Microsystems Inc., October, 2000
 *
 *  GNU Lesser General Public License Version 2.1
 *  =============================================
 *  Copyright 2000 by Sun Microsystems, Inc.
 *  901 San Antonio Road, Palo Alto, CA 94303, USA
 *
 *  This library is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU Lesser General Public
 *  License version 2.1, as published by the Free Software Foundation.
 *
 *  This library is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *  Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public
 *  License along with this library; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston,
 *  MA  02111-1307  USA
 *
 *
 *  Sun Industry Standards Source License Version 1.1
 *  =================================================
 *  The contents of this file are subject to the Sun Industry Standards
 *  Source License Version 1.1 (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.openoffice.org/license.html.
 *
 *  Software provided under this License is provided on an "AS IS" basis,
 *  WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING,
 *  WITHOUT LIMITATION, WARRANTIES THAT THE SOFTWARE IS FREE OF DEFECTS,
 *  MERCHANTABLE, FIT FOR A PARTICULAR PURPOSE, OR NON-INFRINGING.
 *  See the License for the specific provisions governing your rights and
 *  obligations concerning the Software.
 *
 *  The Initial Developer of the Original Code is: Sun Microsystems, Inc.
 *
 *  Copyright: 2000 by Sun Microsystems, Inc.
 *
 *  All Rights Reserved.
 *
 *  Contributor(s): _______________________________________
 *
 *
 ************************************************************************/

#ifndef _TOOLS_SOLMATH_HXX
#define _TOOLS_SOLMATH_HXX

#ifdef MAC
// MAC braucht natuerlich 'ne Extrawurst
#include <mac_start.h>
#include <fp.h>
#include <mac_end.h>
#else
#include <float.h>
#endif

#ifdef SOLARIS
#include <ieeefp.h>
#endif

#include <math.h>

#ifndef _STRING_HXX
#include <tools/string.hxx>
#endif
#ifndef _SOLAR_H
#include <tools/solar.h>
#endif

/// 2**-48 or 3.5527136788005e-015 for ApproxEqual scaling
#define SOMA_MIN_PRECISION	(1.0 / (16777216.0 * 16777216.0))

/// Rounding modes for SolarMath::Round
enum SolarMathRoundingMode
{
	SolarMathRoundCorrected = 0,	/// like HalfUp, but corrects roundoff errors, preferred
	SolarMathRoundDown,		   		/// floor of absolute value, signed return (commercial)
	SolarMathRoundUp,				/// ceil of absolute value, signed return (commercial)
	SolarMathRoundFloor,			/// floor of signed value
	SolarMathRoundCeiling,			/// ceil of signed value
	SolarMathRoundHalfDown,			/// frac <= 0.5 ? floor of abs : ceil of abs, signed return
	SolarMathRoundHalfUp,  			/// frac < 0.5 ? floor of abs : ceil of abs, signed return (mathematical)
	SolarMathRoundHalfEven			/// IEEE rounding mode (statistical)
};

/**
    Methods to work with IEEE 754 double values.
    Special Approx...() methods to iron out roundoff errors.
    Example: 0.1+0.2-0.3 results in 5.55112e-017 instead of 0.0, calling ceil()
    with that value results in 1.0 instead of expected 0.0.
    ApproxSub( 0.1+0.2, 0.3 ) returns 0.0, ApproxCeil( 0.1+0.2-0.3 ) returns 0.0.
 */
class SolarMath
{
public:

    /** Test equalness of two values with an accuracy of the magnitude of the
        given values scaled by 2^-48 (4 bits roundoff stripped).
        @ATTENTION ApproxEqual( value!=0.0, 0.0 ) _never_ yields TRUE */
	static	inline	BOOL	ApproxEqual( double a, double b );

    /** Add two values. If signs differ and the absolute values are equal
        according to ApproxEqual() the method returns 0.0 instead of
        calculating the sum.
        If you wanted to sum up multiple values it would be convenient not to
        call ApproxAdd() for each value but instead remember the first value
        not equal to 0.0, add all other values using normal + operator, and
        with the result and the remembered value call ApproxAdd(). */
	static	inline	double	ApproxAdd( double a, double b );

    /** Substract two values (a-b). If signs are identical and the values are
        equal according to ApproxEqual() the method returns 0.0 instead of
        calculating the substraction. */
	static	inline	double	ApproxSub( double a, double b );

    /** floor() method taking ApproxEqual() into account.
        Use for expected integer values being calculated by double functions.
        @ATTENTION The threshhold value is 3.55271e-015 */
	static	inline	double	ApproxFloor( double a );

    /** ceil() method taking ApproxEqual() into account.
        Use for expected integer values being calculated by double functions.
        @ATTENTION The threshhold value is 3.55271e-015 */
	static	inline	double	ApproxCeil( double a );

    /** For compiler libs not having a finite() function, i.e. OS2BLCI, WINMSCI.
        May be used by the define SOMA_FINITE(double) which should be used
        instead of calling finite() directly and normally maps to finite().
        SOMA_FINITE(d)==0 => value is +INF, -INF or NAN */
	static	inline	int		Finite( double d );

    /** If a value represents +INF or -INF. The sign bit may be queried with
        IsSignBitSet().
        If SOMA_FINITE(d)==0 and IsINF(d)==FALSE then NAN. */
	static	inline	BOOL	IsINF( double d );

	/** Test on any QNAN or SNAN */
	static	inline	BOOL	IsNAN( double d );

    /** Test on a special value 0xFFF8000000000000 used by Microsoft to
        indicate 0/0. This is also a QNAN.
        DO NOT USE IN PORTABLE CODE!
        Could only be interesting for document importing code. */
    static  inline  BOOL    IsIND( double d );

    /** If the sign bit is set */
	static	inline	BOOL	IsSignBitSet( double d );

    /** Set to +INF if bNegative==FALSE or -INF if bNegative==TRUE. */
	static	inline	void	SetINF( double& rd, BOOL bNegative );

	/** Set a QNAN */
    static  inline  void    SetNAN( double& rd, BOOL bNegative = FALSE );

    /** Conversions analogous to sprintf() using internal rounding.
        +/-HUGE_VAL are converted to "1.#INF" and "-1.#INF",
        NAN values are converted to "1.#NAN" and "-1.#NAN",
        of course using cDec instead of '.'

        @param cType
            'F': like sprintf() %f
            'E': like sprintf() %E
            'G': like sprintf() %G, 'F' or 'E' format is used depending on
                 which one is more compact. If nDec==INT_MAX the default value
                 (currently 6) is used.
            'A': Automatic, 'F' or 'E' format is used depending on the double
                 value. If nDec==INT_MAX the highest number of significant
                 decimals possible is generated (bEraseTrailingDecZero=TRUE
                 may be passed where applicable). If nDec!=INT_MAX the
                 specified number of decimal places is generated.

        @param bEraseTrailingDecZero
            Trailing zeros in decimal places are erased.
     */
	static	void			DoubleToString(
								XubString& sStr,
								double fVal,
								xub_Unicode cType,			/// 'F'  'E'  'G'  'A'
                                int nDec,                   /// decimal places
                                xub_Unicode cDec = '.',     /// decimal separator
								BOOL bEraseTrailingDecZero = FALSE
								);

    /** Conversion analogous to strtod(), convert a string representing a
        decimal number into a double value.
        Leading blanks are eaten. If ppEnd!=NULL *ppEnd is set behind the last
        character parsed away. Thus if pStr just contains the numerical string
        to be parsed, **ppEnd=='\0' and *ppEnd-pStr==strlen(pStr) on success.
        Overflow returns +/-HUGE_VAL, underflow 0. In both cases nErrno is set
        to ERANGE, otherwise 0.
        "+/-1.#INF" is recognized as +/-HUGE_VAL, nErrno is set to ERANGE.
        "+/-1.#NAN" is recognized and the value is set to +/-NAN, nErrno is not
        set. */
	static	double			StringToDouble(
								const xub_Unicode* pStr,
								const xub_Unicode cGrpSep,	///	group (AKA thousands) separator
								const xub_Unicode cDecSep,	/// decimal separator
								int& nErrno,
								const xub_Unicode** ppEnd = NULL
								);

	/** Rounds a double value.

		@param fVal
			specifies the value to be rounded.

		@param nDec
			specifies the decimal place where rounding occurs in the range
			-20 to 20. <0 if position is before the decimal point.

		@param eMode
			specifies the rounding mode

	 */
	static	double			Round(
								double fVal,
								short nDec = 0,
								SolarMathRoundingMode eMode = SolarMathRoundCorrected
								);

    /** Scales fVal to a power of 10 without calling pow() or div() for nExp
        values between -16 and +16, providing a faster method.

        @return fVal * pow( 10.0, nExp )
     */
    static  double          Pow10Exp( double fVal, int nExp );

    /** If a value is a valid argument for sin(), cos(), tan().
        IEEE 754 specifies that absolute values up to 2^64 (=1.844e19) for the
        radian must be supported by trigonometric functions.
        Unfortunately, at least on x86 architectures, the FPU doesn't generate
        an error pattern for values >2^64 but produces erroneous results
        instead and sets only the "invalid operation" (IM) flag in the status
        word :-(  Thus the application has to handle it itself. */
    static  inline  BOOL    IsValidArcArg( double d );

    /** Safe sin(), returns NAN if not valid. */
    static  inline  double  Sin( double d );

    /** Safe cos(), returns NAN if not valid. */
    static  inline  double  Cos( double d );

    /** Safe tan(), returns NAN if not valid. */
    static  inline  double  Tan( double d );
};

/** @descr
    Handle FP errors using exceptions and setjmp/longjmp? DON'T in C++!
    Switch off with SOMA_FPCONTROL() and just calculate, test results with
    SOMA_FINITE(d), 0==error (INFinitiy, NaN), else ok.
 */
#if 0
#define SOMA_FPSIGNAL_JUMP 1
#else
#define SOMA_FPSIGNAL_JUMP 0
#endif

// signal() ID for FP exceptions
#ifdef RS6000
#define SOMA_SIGFPE SIGTRAP
#else
#define SOMA_SIGFPE SIGFPE
#endif

// SOMA_FPCONTROL(): switch FP exceptions on/off, depending on SOMA_FPSIGNAL_JUMP
// SOMA_FPRESET(): reinitialize math package
#if defined(WNT) || defined(WIN)
#define SOMA_FPEXCEPTIONS_ON()	_control87( _MCW_EM, 0 )
#define SOMA_FPEXCEPTIONS_OFF()	_control87( _MCW_EM, _MCW_EM )
#define SOMA_FPRESET()			_fpreset()

#elif defined(OS2)
#define SOMA_FPEXCEPTIONS_ON()	_control87( MCW_EM, 0 )
#define SOMA_FPEXCEPTIONS_OFF()	_control87( MCW_EM, MCW_EM )
#define SOMA_FPRESET()			_fpreset()

#elif defined(RS6000)
#define SOMA_FPEXCEPTIONS_ON()	fp_enable_all()
#define SOMA_FPEXCEPTIONS_OFF()	fp_disable_all()
#define SOMA_FPRESET()

#elif defined(LINUX)
#include <fpu_control.h>

#if defined(X86_64)
#include <fenv.h>
#endif

#if defined(POWERPC)

// we need to work around a compiler issue for ppc linux in gcc 3.1

#define _MYFPU_GETCW(cw) ( { \
  union { double d; fpu_control_t cw[2]; } tmp __attribute__ ((__aligned__(8))); \
  __asm__ ("mffs 0; fmr %0, 0" : "=f" (tmp.d) : : "fr0"); \
  (cw)=tmp.cw[1]; \
  tmp.cw[1]; } )

#define _MYFPU_SETCW(cw) { \
  union { double d; fpu_control_t cw[2]; } tmp __attribute__ ((__aligned__(8))); \
  tmp.cw[0] = 0xFFF80000; /* More-or-less arbitrary; this is a QNaN. */ \
  tmp.cw[1] = cw; \
  __asm__ ("lfd%U0 0,%0; mtfsf 255,0" : : "m" (tmp.d) : "fr0"); \
}

#endif

extern inline void glibc_setfpucw( fpu_control_t set )
{
	fpu_control_t cw;

	/* Fetch the current control word.  */

#if defined(POWERPC)
	_MYFPU_GETCW(cw);
#else
	_FPU_GETCW(cw);
#endif

	/* Preserve the reserved bits, and set the rest as the user
	   specified (or the default, if the user gave zero).  */
	cw &= _FPU_RESERVED;
	cw |= set & ~_FPU_RESERVED;

#if defined(POWERPC)
	_MYFPU_SETCW(cw);
#else
	_FPU_SETCW(cw);
#endif

}

#define __setfpucw( control_word ) glibc_setfpucw( control_word )

#if defined(POWERPC)
/* set bit to 1 to enable that exception */
/* _FPU_MASK_ZM | _FPU_MASK_OM | _FPU_MASK_UM | _FPU_MASK_IM */
#define SOMA_FPEXCEPTIONS_ON()	__setfpucw( _FPU_DEFAULT | 0x000000F0 )
#define SOMA_FPEXCEPTIONS_OFF()	__setfpucw( _FPU_DEFAULT )
#elif defined(X86_64)
#define SOMA_FPEXCEPTIONS_ON()	feenableexcept ( FE_ALLEXCEPT )
#define SOMA_FPEXCEPTIONS_OFF()	fedisableexcept( FE_ALLEXCEPT )
#else
#define SOMA_FPEXCEPTIONS_ON()	__setfpucw( _FPU_DEFAULT & ~0x001F )
#define SOMA_FPEXCEPTIONS_OFF()	__setfpucw( _FPU_IEEE )
#endif
#define SOMA_FPRESET()

#else
#define SOMA_FPEXCEPTIONS_ON()
#define SOMA_FPEXCEPTIONS_OFF()
#define SOMA_FPRESET()
#endif

#if SOMA_FPSIGNAL_JUMP
#define SOMA_FPCONTROL() SOMA_FPEXCEPTIONS_ON()
#else
#define SOMA_FPCONTROL() SOMA_FPEXCEPTIONS_OFF()
#endif

/// SOMA_FINITE(d) : test double d on INFINITY, NaN et al.
#if defined(WNT)
#define SOMA_FINITE(d) _finite(d)
#elif defined(MAC)
#define SOMA_FINITE(d) isfinite(d)
#elif defined(LINUX)
#define SOMA_FINITE(d) finite(d)
#elif defined(UNX)
#define SOMA_FINITE(d) finite(d)
#else
// for compiler libs not having finite(), i.e. OS2BLCI, WINMSCI
#define SOMA_FINITE(d) SolarMath::Finite(d)
#endif

#if defined(__IEEEDOUBLE) && (__SIZEOFDOUBLE == 8)

#if defined(__BIGENDIAN)

/// IEEE 754 double structures for BigEndian
typedef union
{
	struct
	{
		UINT32 sign			: 1;
		UINT32 exponent		:11;
		UINT32 fraction_hi	:20;
		UINT32 fraction_lo	:32;
	} inf_parts;
	struct
	{
		UINT32 sign			: 1;
		UINT32 exponent		:11;
		UINT32 qnan_bit		: 1;
		UINT32 bits			:19;
		UINT32 fraction_lo	:32;
	} nan_parts;
	struct
	{
		UINT32 msw			:32;
		UINT32 lsw			:32;
	} w32_parts;
	double value;
} SolarMathDouble;

#elif defined(__LITTLEENDIAN)

/// IEEE 754 double structures for LittleEndian
typedef union
{
	struct {
		UINT32 fraction_lo	:32;
		UINT32 fraction_hi	:20;
		UINT32 exponent		:11;
		UINT32 sign			: 1;
	} inf_parts;
	struct {
		UINT32 fraction_lo	:32;
		UINT32 bits			:19;
		UINT32 qnan_bit		: 1;
		UINT32 exponent		:11;
		UINT32 sign			: 1;
	} nan_parts;
	struct
	{
		UINT32 lsw			:32;
		UINT32 msw			:32;
	} w32_parts;
	double value;
} SolarMathDouble;

#else
#error __FILE__ ": unsupported ENDIAN for SolarMath IEEE 754"
#endif

#else	// !#if defined(__IEEEDOUBLE) && (__SIZEOFDOUBLE == 8)
#error __FILE__ ": don't know how to handle SolarMath IEEE 754"
#endif

inline BOOL SolarMath::ApproxEqual( double a, double b )
{
	if ( a == b )
		return TRUE;
	double x = a - b;
	return (x < 0.0 ? -x : x)  <  ((a < 0.0 ? -a : a) * SOMA_MIN_PRECISION);
}

inline double SolarMath::ApproxAdd( double a, double b )
{
	if ( ((a < 0.0 && b > 0.0) || (b < 0.0 && a > 0.0)) && ApproxEqual( a, -b ) )
		return 0.0;
	return a + b;
}

inline double SolarMath::ApproxSub( double a, double b )
{
	if ( ((a < 0.0 && b < 0.0) || (a > 0.0 && b > 0.0)) && ApproxEqual( a, b ) )
		return 0.0;
	return a - b;
}

inline double SolarMath::ApproxFloor( double a )
{
	double b = floor( a );
	if ( ApproxEqual( a - 1.0, b ) )
		return b + 1.0;
	return b;
}

inline double SolarMath::ApproxCeil( double a )
{
	double b = ceil( a );
	if ( ApproxEqual( a + 1.0, b ) )
		return b - 1.0;
	return b;
}

inline int SolarMath::Finite( double d )
{   // test double on INFINITY, NaN et al.
    // Use with define SOMA_FINITE(d)
    // bit mask in exponent is 0x7FF == 2047
	// sign, exponent, qnan_bit:
	// 0x7FF0 ==  INF
	// 0x7FF1 ==  SNAN  bis 0x7FF7
	// 0x7FF8 ==  QNAN  bis 0x7FFF
	// 0xFFF0 == -INF
	// 0xFFF1 == -SNAN  bis 0xFFF7
    // 0xFFF8 == -IND   (0/0 M$ extension, also a QNAN)
    // 0xFFF9 == -QNAN  bis 0xFFFF
    return ((SolarMathDouble*)&d)->inf_parts.exponent != 0x7ff;
}

inline BOOL SolarMath::IsINF( double d )
{	// exponent==0x7ff fraction==0
	return (SOMA_FINITE(d) == 0) &&
        (((SolarMathDouble*)&d)->inf_parts.fraction_hi == 0)
        && (((SolarMathDouble*)&d)->inf_parts.fraction_lo == 0);
}

inline BOOL SolarMath::IsNAN( double d )
{	// exponent==0x7ff fraction!=0
	return (SOMA_FINITE(d) == 0) && (
        (((SolarMathDouble*)&d)->inf_parts.fraction_hi != 0)
        || (((SolarMathDouble*)&d)->inf_parts.fraction_lo != 0) );
}

// The special M$ extension
inline BOOL SolarMath::IsIND( double d )
{   // sign bit set  exponent==0x7ff  QNAN bit set  remaining fraction==0
    return (((SolarMathDouble*)&d)->w32_parts.msw == 0xfff80000)
        && (((SolarMathDouble*)&d)->w32_parts.lsw == 0);
}

inline BOOL SolarMath::IsSignBitSet( double d )
{
    return ((SolarMathDouble*)&d)->inf_parts.sign != 0;
}

inline void SolarMath::SetINF( double& rd, BOOL bNegative )
{
	((SolarMathDouble&)rd).w32_parts.msw = ( bNegative ? 0xfff00000 : 0x7ff00000 );
	((SolarMathDouble&)rd).w32_parts.lsw = 0;
}

inline void SolarMath::SetNAN( double& rd, BOOL bNegative )
{
	((SolarMathDouble&)rd).w32_parts.msw = ( bNegative ? 0xffffffff : 0x7fffffff );
	((SolarMathDouble&)rd).w32_parts.lsw = 0xffffffff ;
}

inline BOOL SolarMath::IsValidArcArg( double d )
{
    const double fMax = (double)(unsigned long)0x80000000 * (double)(unsigned long)0x80000000 * 2;
    return fabs(d) <= fMax;
}

// We could just call sin() and test the status word for IM flag instead of
// calling IsValidArcArg(), would be less overhead but also platform dependent.
// Similar for cos(), tan()
inline double SolarMath::Sin( double d )
{
    if ( IsValidArcArg( d ) )
        return sin( d );
    SetNAN( d );
    return d;
}

inline double SolarMath::Cos( double d )
{
    if ( IsValidArcArg( d ) )
        return cos( d );
    SetNAN( d );
    return d;
}

inline double SolarMath::Tan( double d )
{
    if ( IsValidArcArg( d ) )
        return tan( d );
    SetNAN( d );
    return d;
}

#endif	// _TOOLS_SOLMATH_HXX

