343 lines
14 KiB
C++
343 lines
14 KiB
C++
/***************************************************************************
|
|
* Math.C *
|
|
* *
|
|
* Some basic math functions. *
|
|
* *
|
|
* HISTORY *
|
|
* Name Date Description *
|
|
* *
|
|
* arvo 06/21/93 Initial coding. *
|
|
* *
|
|
*--------------------------------------------------------------------------*
|
|
* Copyright (C) 1999, James Arvo *
|
|
* *
|
|
* This program is free software; you can redistribute it and/or modify it *
|
|
* under the terms of the GNU General Public License as published by the *
|
|
* Free Software Foundation. See http://www.fsf.org/copyleft/gpl.html *
|
|
* *
|
|
* This program is distributed in the hope that it will be useful, but *
|
|
* WITHOUT EXPRESS OR IMPLIED WARRANTY of merchantability or fitness for *
|
|
* any particular purpose. See the GNU General Public License for more *
|
|
* details. *
|
|
* *
|
|
***************************************************************************/
|
|
#include <math.h>
|
|
#include <stdlib.h>
|
|
#include <iostream>
|
|
#include <assert.h>
|
|
#include "ArvoMath.h"
|
|
#include "form.h"
|
|
|
|
namespace ArvoMath {
|
|
static const float Epsilon = 1.0E-5;
|
|
static const double LogTwo = log( 2.0 );
|
|
|
|
#define BinCoeffMax 500
|
|
|
|
double RelErr( double x, double y )
|
|
{
|
|
double z = x - y;
|
|
if( x < 0.0 ) x = -x;
|
|
if( y < 0.0 ) y = -y;
|
|
return z / ( x > y ? x : y );
|
|
}
|
|
|
|
/***************************************************************************
|
|
* A R C Q U A D *
|
|
* *
|
|
* Returns the theta / ( 2*PI ) where the input variables x and y are *
|
|
* such that x == COS( theta ) and y == SIN( theta ). *
|
|
* *
|
|
***************************************************************************/
|
|
float ArcQuad( float x, float y )
|
|
{
|
|
if( Abs( x ) > Epsilon )
|
|
{
|
|
float temp = OverTwoPi * atan( Abs( y ) / Abs( x ) );
|
|
if( x < 0.0 ) temp = 0.5 - temp;
|
|
if( y < 0.0 ) temp = 1.0 - temp;
|
|
return( temp );
|
|
}
|
|
else if( y > Epsilon ) return( 0.25 );
|
|
else if( y < -Epsilon ) return( 0.75 );
|
|
else return( 0.0 );
|
|
}
|
|
|
|
/***************************************************************************
|
|
* A R C T A N *
|
|
* *
|
|
* Returns the angle theta such that x = COS( theta ) & y = SIN( theta ). *
|
|
* *
|
|
***************************************************************************/
|
|
float ArcTan( float x, float y )
|
|
{
|
|
if( Abs( x ) > Epsilon )
|
|
{
|
|
float temp = atan( Abs( y ) / Abs( x ) );
|
|
if( x < 0.0 ) temp = Pi - temp;
|
|
if( y < 0.0 ) temp = TwoPi - temp;
|
|
return( temp );
|
|
}
|
|
else if( y > Epsilon ) return( PiOverTwo );
|
|
else if( y < -Epsilon ) return( 3 * PiOverTwo );
|
|
else return( 0.0 );
|
|
}
|
|
|
|
/***************************************************************************
|
|
* M A C H I N E E P S I L O N *
|
|
* *
|
|
* Returns the machine epsilon. *
|
|
* *
|
|
***************************************************************************/
|
|
float MachineEpsilon()
|
|
{
|
|
float x = 1.0;
|
|
float y;
|
|
float z = 1.0 + x;
|
|
while( z > 1.0 )
|
|
{
|
|
y = x;
|
|
x /= 2.0;
|
|
z = (float)( 1.0 + (float)x ); // Avoid double precision!
|
|
}
|
|
return (float)y;
|
|
}
|
|
|
|
/***************************************************************************
|
|
* L O G G A M M A *
|
|
* *
|
|
* Computes the natural log of the gamma function using the Lanczos *
|
|
* approximation formula. Gamma is defined by *
|
|
* *
|
|
* ( z - 1 ) -t *
|
|
* gamma( z ) = Integral[ t e dt ] *
|
|
* *
|
|
* *
|
|
* where the integral ranges from 0 to infinity. The gamma function *
|
|
* satisfies *
|
|
* gamma( n + 1 ) = n! *
|
|
* *
|
|
* This algorithm has been adapted from "Numerical Recipes", p. 157. *
|
|
* *
|
|
***************************************************************************/
|
|
double LogGamma( double x )
|
|
{
|
|
static const double
|
|
coeff0 = 7.61800917300E+1,
|
|
coeff1 = -8.65053203300E+1,
|
|
coeff2 = 2.40140982200E+1,
|
|
coeff3 = -1.23173951600E+0,
|
|
coeff4 = 1.20858003000E-3,
|
|
coeff5 = -5.36382000000E-6,
|
|
stp = 2.50662827465E+0,
|
|
half = 5.00000000000E-1,
|
|
fourpf = 4.50000000000E+0,
|
|
one = 1.00000000000E+0,
|
|
two = 2.00000000000E+0,
|
|
three = 3.00000000000E+0,
|
|
four = 4.00000000000E+0,
|
|
five = 5.00000000000E+0;
|
|
double r = coeff0 / ( x ) + coeff1 / ( x + one ) +
|
|
coeff2 / ( x + two ) + coeff3 / ( x + three ) +
|
|
coeff4 / ( x + four ) + coeff5 / ( x + five ) ;
|
|
double s = x + fourpf;
|
|
double t = ( x - half ) * log( s ) - s;
|
|
return t + log( stp * ( r + one ) );
|
|
}
|
|
|
|
/***************************************************************************
|
|
* L O G F A C T *
|
|
* *
|
|
* Returns the natural logarithm of n factorial. For efficiency, some *
|
|
* of the values are cached, so they need be computed only once. *
|
|
* *
|
|
***************************************************************************/
|
|
double LogFact( int n )
|
|
{
|
|
static const int Cache_Size = 100;
|
|
static double c[ Cache_Size ] = { 0.0 }; // Cache some of the values.
|
|
if( n <= 1 ) return 0.0;
|
|
if( n < Cache_Size )
|
|
{
|
|
if( c[n] == 0.0 ) c[n] = LogGamma((double)(n+1));
|
|
return c[n];
|
|
}
|
|
return LogGamma((double)(n+1)); // gamma(n+1) == n!
|
|
}
|
|
|
|
/***************************************************************************
|
|
* M U L T I N O M I A L C O E F F *
|
|
* *
|
|
* Returns the multinomial coefficient ( n; X1 X2 ... Xk ) which is *
|
|
* defined to be n! / ( X1! X2! ... Xk! ). This is done by computing *
|
|
* exp( log(n!) - log(X1!) - log(X2!) - ... - log(Xk!) ). The value of *
|
|
* n is obtained by summing the Xi's. *
|
|
* *
|
|
***************************************************************************/
|
|
double MultinomialCoeff( int k, int X[] )
|
|
{
|
|
int i;
|
|
// Find n by summing the coefficients.
|
|
|
|
int n = X[0];
|
|
for( i = 1; i < k; i++ ) n += X[i];
|
|
|
|
// Compute log(n!) then subtract log(X!) for each X.
|
|
|
|
double LogCoeff = LogFact( n );
|
|
for( i = 0; i < k; i++ ) LogCoeff -= LogFact( X[i] );
|
|
|
|
// Round the exponential of the result to the nearest integer.
|
|
|
|
return floor( exp( LogCoeff ) + 0.5 );
|
|
}
|
|
|
|
|
|
double MultinomialCoeff( int i, int j, int k )
|
|
{
|
|
int n = i + j + k;
|
|
double x = LogFact( n ) - LogFact( i ) - LogFact( j ) - LogFact( k );
|
|
return floor( exp( x ) + 0.5 );
|
|
}
|
|
|
|
/***************************************************************************
|
|
* B I N O M I A L C O E F F S *
|
|
* *
|
|
* Generate all n+1 binomial coefficents for a given n. This is done by *
|
|
* computing the n'th row of Pascal's triange, starting from the top. *
|
|
* No additional storage is required. *
|
|
* *
|
|
***************************************************************************/
|
|
void BinomialCoeffs( int n, long *coeff )
|
|
{
|
|
coeff[0] = 1;
|
|
for( int i = 1; i <= n; i++ )
|
|
{
|
|
long a = coeff[0];
|
|
long b = coeff[1];
|
|
for( int j = 1; j < i; j++ ) // Make next row of Pascal's triangle.
|
|
{
|
|
coeff[j] = a + b; // Overwrite the old row.
|
|
a = b;
|
|
b = coeff[j+1];
|
|
}
|
|
coeff[i] = 1; // The last entry in any row is always 1.
|
|
}
|
|
}
|
|
|
|
void BinomialCoeffs( int n, double *coeff )
|
|
{
|
|
coeff[0] = 1.0;
|
|
for( int i = 1; i <= n; i++ )
|
|
{
|
|
double a = coeff[0];
|
|
double b = coeff[1];
|
|
for( int j = 1; j < i; j++ ) // Make next row of Pascal's triangle.
|
|
{
|
|
coeff[j] = a + b; // Overwrite the old row.
|
|
a = b;
|
|
b = coeff[j+1];
|
|
}
|
|
coeff[i] = 1.0; // The last entry in any row is always 1.
|
|
}
|
|
}
|
|
|
|
const double *BinomialCoeffs( int n )
|
|
{
|
|
static double *coeff[ BinCoeffMax + 1 ] = { 0 };
|
|
if( n > BinCoeffMax || n < 0 )
|
|
{
|
|
std::cerr << form( "%d is outside of (0,%d) in BinomialCoeffs", n, BinCoeffMax );
|
|
return NULL;
|
|
}
|
|
if( coeff[n] == NULL ) // Fill in this entry.
|
|
{
|
|
double *c = new double[ n + 1 ];
|
|
if( c == NULL )
|
|
{
|
|
std::cerr << form( "Could not allocate for BinomialCoeffs(%d)", n );
|
|
return NULL;
|
|
}
|
|
BinomialCoeffs( n, c );
|
|
coeff[n] = c;
|
|
}
|
|
return coeff[n];
|
|
}
|
|
|
|
/***************************************************************************
|
|
* B I N O M I A L C O E F F *
|
|
* *
|
|
* Compute a given binomial coefficient. Several rows of Pascal's *
|
|
* triangle are stored for efficiently computing the small coefficients. *
|
|
* Higher-order terms are computed using LogFact. *
|
|
* *
|
|
***************************************************************************/
|
|
double BinomialCoeff( int n, int k )
|
|
{
|
|
double b;
|
|
int p = n - k;
|
|
if( k <= 1 || p <= 1 ) // Check for errors and special cases.
|
|
{
|
|
if( k == 0 || p == 0 ) return 1;
|
|
if( k == 1 || p == 1 ) return n;
|
|
std::cerr << form( "BinomialCoeff(%d,%d) is undefined", n, k );
|
|
return 0;
|
|
}
|
|
static const int // Store part of Pascal's triange for small coeffs.
|
|
n0[] = { 1 },
|
|
n1[] = { 1, 1 },
|
|
n2[] = { 1, 2, 1 },
|
|
n3[] = { 1, 3, 3, 1 },
|
|
n4[] = { 1, 4, 6, 4, 1 },
|
|
n5[] = { 1, 5, 10, 10, 5, 1 },
|
|
n6[] = { 1, 6, 15, 20, 15, 6, 1 },
|
|
n7[] = { 1, 7, 21, 35, 35, 21, 7, 1 },
|
|
n8[] = { 1, 8, 28, 56, 70, 56, 28, 8, 1 },
|
|
n9[] = { 1, 9, 36, 84, 126, 126, 84, 36, 9, 1 };
|
|
switch( n )
|
|
{
|
|
case 0 : b = n0[k]; break;
|
|
case 1 : b = n1[k]; break;
|
|
case 2 : b = n2[k]; break;
|
|
case 3 : b = n3[k]; break;
|
|
case 4 : b = n4[k]; break;
|
|
case 5 : b = n5[k]; break;
|
|
case 6 : b = n6[k]; break;
|
|
case 7 : b = n7[k]; break;
|
|
case 8 : b = n8[k]; break;
|
|
case 9 : b = n9[k]; break;
|
|
default:
|
|
{
|
|
double x = LogFact( n ) - LogFact( p ) - LogFact( k );
|
|
b = floor( exp( x ) + 0.5 );
|
|
}
|
|
}
|
|
return b;
|
|
}
|
|
|
|
|
|
/***************************************************************************
|
|
* L O G D O U B L E F A C T (Log of double factorial) *
|
|
* *
|
|
* Return log( n!! ) where the double factorial is defined by *
|
|
* *
|
|
* (2 n + 1)!! = 1 * 3 * 5 * ... * (2n + 1) (Odd integers) *
|
|
* *
|
|
* (2 n)!! = 2 * 4 * 6 * ... * 2n (Even integers) *
|
|
* *
|
|
* and is related to the single factorial via *
|
|
* *
|
|
* (2 n + 1)!! = (2 n + 1)! / ( 2^n n! ) (Odd integers) *
|
|
* *
|
|
* (2 n)!! = 2^n n! (Even integers) *
|
|
* *
|
|
***************************************************************************/
|
|
double LogDoubleFact( int n ) // log( n!! )
|
|
{
|
|
int k = n / 2;
|
|
double f = LogFact( k ) + k * LogTwo;
|
|
if( Odd(n) ) f = LogFact( n ) - f;
|
|
return f;
|
|
}
|
|
};
|