fp.java

国外人写的操作感挺好的一款滑雪游戏.结构也比较清楚

// Copyright (c) 2004 Venan Entertainment, Inc. All rights reserved.
//
// Venan Entertainment, Inc.,  Middletown, Connecticut 06457
// http://www.venan.com

// Original Source:
// Copyright (c) 2001 Beartronics 
// Author: Henry Minsky ([e-mail])
// Licensed under terms "Artistic License" http://www.opensource.org/licenses/artistic-license.html
//
// Numerical algorithms based on
//	http://www.cs.clemson.edu/html_docs/SUNWspro/common-tools/numerical_comp_guide/ncg_examples.doc.html
// Trig routines based on numerical algorithms described in 
//	http://www.magic-software.com/MgcNumerics.html
//	http://www.dattalo.com/technical/theory/logs.html

public class FP
{
	public final static String toString( int x )
	{
		String s = new String( );
		if ( x<0 )
		{
			x = -x;
			s = "-";
		}
		int whole = x >> 16;
		long frac = (((long)(x & 0x0000ffff)) * 100000) >> 16;

		Integer b = new Integer( whole );
		s += b.toString( );
		s += ".";
		String f = new String( );
		f = Long.toString( frac );
		for ( int i=0; i<5-f.length( ); i++ )
			s += "0";
		s += f;

		return s;
	}

	/** Convert a 16:16 fixed-point to an int
	*/
	public final static int toInt( int x )
	{
		return x>>16;
	}

	/** Convert an int to a 16:16 fixed-point 
	*/

	public final static int intToFP( int x )
	{
		return x<<16;
	}

	public final static int fracPart( int x )
	{
		return x&0xffff;
	}

	/** Multiply two fixed-point numbers */
	public final static int mul( int x, int y )
	{
		long z = (long) x * (long) y;
		return ( (int) (z >> 16) );
	}

	/** Divides two fixed-point numbers */
	public final static int div( int x, int y )
	{
		long z = (((long) x) << 32);
		return (int) ((z / y) >> 16);
	}


	/** Compute square-root of a 16:16 fixed point number */
	public final static int sqrt( long n )
	{
		long s = (n + 0x10000) >> 1;
		for ( int i = 0; i < 8; i++ )
		{
			//converge six times
			s = (s + (n << 16)/s) >> 1;
		}
		return (int)s;
	}

	/** Vector length */
	public final static int getMag( int x, int y )
	{
		long n = ((long) x * (long) x + (long) y * (long) y) >> 16;
		return sqrt(n);
	}

	/** Round to nearest fixed point integer */
	public final static int round( int n )
	{
		if ( n > 0 )
		{
			if ( (n & 0x8000) != 0 )
			{
				return ( ((n+0x10000)>>16)<<16 );
			}
			else
			{
				return ( ((n)>>16)<<16 );
			}
		}
		else
		{
			int k;
			n = -n;
			if ( (n & 0x8000) != 0 )
			{
				k = (((n+0x10000)>>16)<<16);
			}
			else
			{
				k = (((n)>>16)<<16);
			}
			return -k;
		}
	}

	public final static int TWO_PI = 411775;
	public final static int PI = 205887;
	public final static int PI_OVER_2 = 102944;
	public final static int E = 178145;
	public final static int HALF = (1<<15);
	public final static int ONE = (1<<16);
	public final static int TWO = (2<<16);
	public final static int MILLI = 66;			// 1/1000
	public final static int ONE_OVER_SQRT_2 = 46341;
	public final static int SIXTH = 10923;

	private final static int INVFACT7 = 13;		// 1/7!
	private final static int INVFACT6 = 91;		// 1/6!
	private final static int INVFACT5 = 546;	// 1/5!
	private final static int INVFACT4 = 2731;	// 1/4!
	private final static int INVFACT3 = 10923;	// 1/3!
	private final static int INVFACT2 = 32768;	// 1/2!

	private final static int TK9 = 1433;		// 62/2835
	private final static int TK7 = 3537;		// 17/315
	private final static int TK5 = 8738;		// 2/15
	private final static int TK3 = 21845;		// 1/3

	private final static int AT9 = 7282;		// 1/9
	private final static int AT7 = 9362;		// 1/7
	private final static int AT5 = 13107;		// 1/5
	private final static int AT3 = 21845;		// 1/3

	private final static int AS9 = 1991;		// 1/2 * 3/4 * 5/6 * 7/8 * 1/9
	private final static int AS7 = 2926;		// 1/2 * 3/4 * 5/6 * 1/7
	private final static int AS5 = 4915;		// 1/2 * 3/4 * 1/5
	private final static int AS3 = 10923;		// 1/2 * 1/3

	private final static int log2arr[] =
	{
		26573,	// ln(3/2)
		14624,	// ln(5/4)
		7719,	// ln(9/8)
		3973,	// ln(17/16)
		2017,	// etc.
		1016,
		510,
		256,
		128,
		64,
		32,
		16,
		8,
		4,
		2,
		1,
		0,
		0
    };

	private final static int lnscale[] =
	{
		0,			// ln(2) * 0
		45426,		// ln(2) * 1
		90852,		// ln(2) * 2
		136278,		// ln(2) * 3
		181704,		// etc.
		227130,
		272557,
		317983,
		363409,
		408835,
		454261,
		499687,
		545113,
		590539,
		635965,
		681391,
		626817,
		772244
    };

	/**
	* For the inverse tangent calls, all approximations are valid for |t| <= 1.
	* To compute ATAN(t) for t > 1, use ATAN(t) = PI/2 - ATAN(1/t).  For t < -1,
	* use ATAN(t) = -PI/2 - ATAN(1/t).
	*/

	/** Computes SIN(f), f is a fixed point number in radians.
	* 0 <= f < 2PI
	*/
	public final static int sin( int f )
	{
		// Taylor series expansion
		// If in range [0,pi/2): nothing needs to be done.
		// otherwise, we need to get f into that range and account for
		// sign change.
		if ( f < 0 )
			f += FP.TWO_PI;
		else if ( f >= FP.TWO_PI )
			f -= FP.TWO_PI;

		int sign = 1;
		if ( (f >= PI_OVER_2) && (f < PI) )
		{
			f = PI - f;
		}

		if ( (f >= PI) && (f < (PI + PI_OVER_2)) )
		{
			f = f - PI;
			sign = -1;
		}

		if ( f >= (PI + PI_OVER_2) )
		{
			f = TWO_PI - f;
			sign = -1;
		}

		int sqr = mul( f,f );
		int result = -INVFACT7;
		result = mul( result, sqr );
		result += INVFACT5;
		result = mul( result, sqr );
		result -= INVFACT3;
		result = mul( result, sqr );
		result += ONE;
		result = mul( result, f );
		return sign * result;
	}

	/** Computes COS(f), f is a fixed point number in radians.
	* 0 <= f < 2PI
	*/
	public final static int cos( int f )
	{
		// Use an alternate (smaller) version of cos
		// sin^2 + cos^2 = 1
		int iSin = sin(f);
		return FP.sqrt( ONE - FP.mul(iSin,iSin));
/*
		// If in range [0,pi/2): nothing needs to be done.
		// otherwise, we need to get f into that range and account for
		// sign change.
		if ( f < 0 )
			f += FP.TWO_PI;
		else if ( f >= FP.TWO_PI )
			f -= FP.TWO_PI;

		int sign = 1;
		if ( (f >= PI_OVER_2) && (f < PI) )
		{
			f = PI - f;
			sign = -1;
		}
		else if ( (f >= PI_OVER_2) && (f < (PI + PI_OVER_2)) )
		{
			f = f - PI;
			sign = -1;
		}
		else if ( f >= (PI + PI_OVER_2) )
		{
			f = TWO_PI - f;
		}

		int sqr = mul( f,f );
		int result = -INVFACT6;
		result = mul( result, sqr );
		result += INVFACT4;
		result = mul( result, sqr );
		result -= INVFACT2;
		result = mul( result, sqr );
		result += ONE;
		return result * sign;
*/
	}


	/** Computes Tan(f), f is a fixed point number in radians.
	* -PI/2 <= f < PI/2
	*/

	public final static int tan( int f )
	{
		if ( f < -FP.PI_OVER_2 )
			f += FP.PI;
		else if ( f >= FP.PI_OVER_2 )
			f -= FP.PI;
		
		int sqr = mul( f,f );
		int result = TK9;
		result = mul( result, sqr );
		result += TK7;
		result = mul( result, sqr );
		result += TK5;
		result = mul( result, sqr );
		result += TK3;
		result = mul( result, sqr );
		result += ONE;
		result = mul( result, f );
		return result;
	}


	public final static int atan2( int y, int x )
	{
		if (x != 0)
			return FP.atan( FP.div( y, x ));
		else if (y != 0)
			return PI_OVER_2;
		else
			return 0;
	}

	/** Computes ArcTan(f), f is a fixed point number
	* |f| <= 1
	* <p>
	* For the inverse tangent calls, all approximations are valid for |t| <= 1.
	* To compute ATAN(t) for t > 1, use ATAN(t) = PI/2 - ATAN(1/t).  For t < -1,
	* use ATAN(t) = -PI/2 - ATAN(1/t).
	*/

	public final static int atan( int f )
	{
		int sqr = mul( f,f );

		int result = AT9;
		if (f >= FP.ONE)
		{
			result = div( result, sqr );
			result -= AT7;
			result = div( result, sqr );
			result += AT5;
			result = div( result, sqr );
			result -= AT3;
			result = div( result, sqr );
			result += ONE;
			result = div( result, -f );
			result += PI_OVER_2;
		}
		else if (f <= -FP.ONE)
		{
			result = div( result, sqr );
			result -= AT7;
			result = div( result, sqr );
			result += AT5;
			result = div( result, sqr );
			result -= AT3;
			result = div( result, sqr );
			result += ONE;
			result = div( result, -f );
			result -= PI_OVER_2;
		}
		else
		{
			result = mul( result, sqr );
			result -= AT7;
			result = mul( result, sqr );
			result += AT5;
			result = mul( result, sqr );
			result -= AT3;
			result = mul( result, sqr );
			result += ONE;
			result = mul( result, f );
		}
		return result;
	}

	/** Compute ArcSin(f), 0 <= f <= 1
	*/

	public final static int asin( int f )
	{
		int sqr = mul( f,f );
		int result = AS9;
		result = mul( result, sqr );
		result += AS7;
		result = mul( result, sqr );
		result += AS5;
		result = mul( result, sqr );
		result += AS3;
		result = mul( result, sqr );
		result += ONE;
		result = mul( result, f );
		return result;
	}


	/** Compute ArcCos(f), 0 <= f <= 1
	*/

	public final static int acos( int f )
	{
		return PI_OVER_2 - asin( f );
	}


	/** Power
		x^y
	*/
	public final static int pow( int x, int y )
	{
		return exp(mul(y,ln( x )));
	}	

	/** Exponential
	*/
	public final static int exp( int x )
	{
		// Taylor series approximation
		int s = ONE;
		int g = s;
		for ( int i=1; i<12; i++ )
		{
			int f = FP.div( x,FP.intToFP( i ) );
			s = FP.mul( s,f );
			g += s;
		}

		return g;
	}


	/** Logarithms: 
	* 
	* (2) Knuth, Donald E., "The Art of Computer Programming Vol 1",
	* Addison-Wesley Publishing Company, ISBN 0-201-03822-6 ( this
	* comes from Knuth (2), section 1.2.3, exercise 25).
	*
	* http://www.dattalo.com/technical/theory/logs.html
	*

	*/

	/** This table is created using base of e. 

	(defun fixedpoint (z) 
	  (round (* z (lsh 1 16))))

	(loop for k from 0 to 16 do
	      (setq z (log (+ 1 (expt 2.0 (- (+ k 1)))))) 
	      (insert (format "%d\n"  (fixedpoint z))))


	*/

	/*
	Binary Logarithm:

	case is very similar to the previous one. The only difference is
	in how the input is factored. Like before we are given:

	Input: 16 bit unsigned integer x; 0 < x < 65536

	(or 8 bit unsigned integer...)

	Output: g, lg(x) the logarithm of x with respect to base 2.

	3(b).i) Create a table of logarithms of the following constants:
	log2arr[i] = lg(1 + 2^(-(i+1))) 

	i = 0..M, M == desired size of the table.

	The first few values of the array are

	lg(3/2), lg(5/4), lg(9/8), lg(17/16),...

	Recall that in the previous case the factors were

	lg(2/1), lg(4/3), lg(8/7), lg(16/15),...

	Again, if you wish to compute logarithms to a different base,
	then substitute the lg() function with the appropriate based
	logarithm function.

	3(b).ii)Scale y to a value between 1 and 2.
	This is identical to 3(a).ii.

	3(b).iii) Changing Perspective
	Again, this is identical to 3(a).iii.

	3(b).iv) Factor y.

	This is very similar to step (iv) above. However, we now have
	different factors. Using the same example, x = 1.9, we can find
	the factors for this case.

	a) 1.9 > 1.5,
	x = x/1.5 ==> 1.266666
	b) 1.26 >  1.25
	x ==> 1.0133333
	c) 1.0133 < 1.125 so don't divide
	d) 1.0133 < 1.0625  "     "
	e) etc.


	So, x ~= 1.5 * 1.25 * etc.


	Like the previouse case, these factors are not perfect. Also,
	they're somewhat redundant in the sense that
	1.5*1.25*1.125*1.0625*... spans a range that is larger than 2 (
	~2.38423 for i<=22). So unlike the previous factoring method,
	this one will not have repeated factors. Here's some psuedo
	code:

	for(i=1,d=0.5; i<M; i++, d/=2)
	if( x > 1+d)
	{
	x /= (1+d);
	g += log2arr[i-1];   // log2arr[i-1] = log2(1+d);
	}


	Here, d takes on the values of 0.5, 0.25, 0.125, ... , 2^(-i). Then
	1+d is the trial factor at each step. If x is greater than this trial
	factor, then we divide the trial factor out and add to g (ultimately
	the logarithm of x) the partial logarithm of the factor.

	*/

	/*
	(loop for k from 0 to 16 do
	(setq z (log (expt 2 k)))
	(insert (format "%d,\n" (fixedpoint z))))
	*/



	public final static int ln( int x )
	{
		// prescale so x is between 1 and 2
		int shift = 0;

		while ( x > 1<<17 )
		{
			shift++;
			x >>= 1;
		}

		int g = 0;
		int d = HALF;
		for ( int i = 1; i < 16; i++ )
		{
			if ( x > (ONE + d) )
			{
				x = div( x, ( ONE + d) );
				g += log2arr[i-1];   // log2arr[i-1] = ln(1+d);
			}
			d >>= 1;
		}
		return g + lnscale[shift];
	}

	public final static int splint( final int xa[], final int ya[], final int y2a[], int n, int x )
	{
		int klo,khi,k;
		int h,b,a;
		
		klo = 0;
		khi = n-1;
		
		while (khi-klo > 1)
		{
			k = (khi+klo) >> 1;
			if (xa[k] > x)
				khi=k;
			else
				klo=k;
		}
	
		h = xa[khi]-xa[klo];
		a = FP.div( xa[khi]-x, h );
		b = FP.div( x-xa[klo], h );
	
		// evaluate the cubic spline polynomial
		return FP.mul( a, ya[klo] ) +
			FP.mul( b, ya[khi] ) +
			FP.mul(
	        	FP.mul( FP.mul( FP.mul( a, a ) - FP.ONE, xa[khi]-x), y2a[klo]) +
				FP.mul( FP.mul( FP.mul( b, b ) - FP.ONE, x-xa[klo]), y2a[khi]),
				FP.mul( h, FP.SIXTH )
			);
	}
	
}