Fixed Point Library


goobers

Member
Joined
Jan 22, 2007
Messages
344
Location
The UK
Website
Visit site
Hey, I'm working on a Mode7 style engine for the GP2X at the moment, and I think I'm gonna have to convert all of the calculations to fixed point to get a resonable speed out of it. I've got working Mode7 code that I wrote in a different language (with floats) that I'm rewriting in C++

What is the best way to do this?

I've had a look at math-sll.h, but I'm not entirely sure how I'd use the functions to replace what i've already got...

for instance, I've got a line as follows:

float flrPointDist = (strtFlrPointDist / cos(convRad((countX * 0.1875) -30)));

(convRad is a function I wrote)

Not been coding in C++ very long, hoping one of you guys knows what I need to do :p
 
Any reason not to use opengl and the gpu940 library?

Check out how cool egoboo is.
 
As much as I hate to come up with answers that do not quite match to what was asked for, I think the best is not to use a library or anything but to do it all by hand using nothing but good old integers. You can get much better results, considering the following assumptions.

Let's say that all you really need for the coordinates of your mode 7 engine is a 8.8 format, meaning your coodinates go from -128 to +127.996, with a step of 0.003906 between each possible value. Well the reason why it's much faster than using a 16.16 format, is that with a 16.16 format, a multiplication will look like this :

CODE
inline int QUNSIGN(fixed *a)
{
if (*a >= 0)
return 0;
*a = -*a;
return 1;
}

fixed QMUL(fixed x, fixed y)
{
int sign = QUNSIGN(&x) ^ QUNSIGN(&y);
fixed xf = QFRAC(x), yf = QFRAC(y), res;
Uint32 r = xf*yf;
r >>= 16;

x = QINT(x);
y = QINT(y);
r += x*yf + y*xf;
r += ((x*y)<<16);

res = (r & 0x7FFFFFFF);
return (sign ? -res : res);
}


as with a 8.8 format, a multiplication will look like this :

CODE
r = x * y >> 16;


I let you count how many cycles each would take. My point is, while a fixed point library is nice and all, the best thing, performance wise and in order to suit your needs, is to define precisely what you're going to need in terms of range or precision, and to not use more than needed. That's how you gotta do if you want to achieve the best performance.

Not too complicated, that's what I did for the physics of Sponge Blob Tennis, and it worked out great! I used a 10.6 format there, and did something special for that format, and now my game runs at 32 MHz at 80 FPS, and yet its main bottleneck has to do with copying stuff!

But yeah, the downside of this is that you have to think. Not too complicated, but yeah, you have to think.
 
This is the most recent I have. I haven't used it on the GP2X yet, but the ARM ASM parts should work anyway. If yo know how to use functions this is pretty easy to use...

CODE

/*
* $Id: math-sll.c,v 1.15 2002/08/20 18:01:54 andrewm Exp $
*
* Changed by CHUI
*
* Purpose
* A fixed point (31.32 bit) math library.
*
* Description
* Floating point packs the most accuracy in the available bits, but it
* often provides more accuracy than is required. It is time consuming to
* carry the extra precision around, particularly on platforms that don't
* have a dedicated floating point processor.
*
* This library is a compromise. All math is done using the 64 bit signed
* "long long" format (sll), and is not intended to be portable, just as
* fast as possible. Since "long long" is a elementary type, it can be
* passed around without resorting to the use of pointers. Since the
* format used is fixed point, there is never a need to do time consuming
* checks and adjustments to maintain normalized numbers, as is the case
* in floating point.
*
* Simply put, this library is limited to handling numbers with a whole
* part of up to 2^31 - 1 = 2.147483647e9 in magnitude, and fractional
* parts down to 2^-32 = 2.3283064365e-10 in magnitude. This yields a
* decent range and accuracy for many applications.
*
* IMPORTANT
* No checking for arguments out of range (error).
* No checking for divide by zero (error).
* No checking for overflow (error).
* No checking for underflow (warning).
* Chops, doesn't round.
*
* Functions
* sll dbl2sll(double x) double -> sll
* double slldbl(sll x) sll -> double
*
* sll slladd(sll x, sll y) x + y
* sll sllsub(sll x, sll y) x - y
* sll sllmul(sll x, sll y) x * y
* sll slldiv(sll x, sll y) x / y
*
* sll sllinv(sll v) 1 / x
* sll sllmul2(sll x) x * 2
* sll sllmul4(sll x) x * 4
* sll sllmul2n(sll x, int n) x * 2^n, 0 <= n <= 31
* sll slldiv2(sll x) x / 2
* sll slldiv4(sll x) x / 4
* sll slldiv2n(sll x, int n) x / 2^n, 0 <= n <= 31
*
* sll sllcos(sll x) cos x
* sll sllsin(sll x) sin x
* sll slltan(sll x) tan x
* sll sllatan(sll x) atan x
*
* sll sllexp(sll x) e^x
* sll slllog(sll x) ln x
*
* sll sllpow(sll x, sll y) x^y
* sll sllsqrt(sll x) x^(1 / 2)
*
*/

#ifndef MATHSLL_H
#define MATHSLL_H


#ifdef __cplusplus
extern "C" {
#endif

#ifndef USE_FIXED_POINT

#include <math.h>

/*
typedef float sll;
typedef float ull;
*/

typedef double sll;
typedef double ull;

#define int2sll(X) ((sll) (X))
#define sllvalue(X) (X)
#define sll2int(X) ((int) (X))
#define sll_abs(X) ((sll)fabs((double)X))
#define sllabs(X) ((sll)fabs((double)X))
#define sllint(X) ((sll)floor((double)X))
#define sllfrac(X) ((X) - ((sll)sll2int(X)))
#define sllneg(X) (-(X))
#define slladd(X,Y) ((X) + (Y))
#define sllsub(X,Y) ((X) - (Y))
#define sllmul(X,Y) ((X) * (Y))
#define slldiv(X,Y) ((X) / (Y))
#define sll2dbl(X) ((double)(X))
#define dbl2sll(X) ((sll)(X))
#define sllinv(X) (((sll)1.0) / (X))
#define sllmul2(X) (((sll)2.0) * (X))
#define sllmul4(X) (((sll)4.0) * (X))
#define sllmul2n(X,N) (((sll)(1<<N)) * (X))
#define slldiv2(X) ((X) / ((sll)2.0))
#define slldiv4(X) ((X) / ((sll)4.0))
#define slldiv2n(X,N) ((X) / ((sll)(1<<N)))
#define sllcos(X) ((sll)cos((double)(X)))
#define sllsin(X) ((sll)sin((double)(X)))
#define slltan(X) ((sll)tan((double)(X)))
#define sllatan(X) ((sll)atan((double)(X)))
#define sllexp(X) ((sll)exp((double)(X)))
#define slllog(X) ((sll)log((double)(X)))
#define sllpow(X) ((sll)pow((double)(X)))
#define sllsqrt(X) ((sll)sqrt((double)(X)))



#define sllrotr(X,N) ((X) / ((sll)(1<<N)))
#define sllrotl(X,N) ((X) * ((sll)(1<<N)))

#define SLL_CONST_0 ((sll)0.0)
#define SLL_CONST_1 ((sll)1.0)
#define SLL_CONST_2 ((sll)2.0)
#define SLL_CONST_3 ((sll)3.0)
#define SLL_CONST_4 ((sll)4.0)
#define SLL_CONST_10 ((sll)10.0)
#define SLL_CONST_15 ((sll)15.0)
#define SLL_CONST_1_2 ((sll)0.5)
#define SLL_CONST_1_3 ((sll)0.33333333333333333333333333333333)
#define SLL_CONST_1_4 ((sll)0.25)
#define SLL_CONST_1_5 ((sll)0.2)
#define SLL_CONST_1_6 ((sll)0.16666666666666666666666666666667)
#define SLL_CONST_1_7 ((sll)0.14285714285714285714285714285714)
#define SLL_CONST_1_8 ((sll)0.125)
#define SLL_CONST_1_9 ((sll)0.11111111111111111111111111111111)
#define SLL_CONST_1_10 ((sll)0.10)
#define SLL_CONST_1_11 ((sll)0.090909090909090909090909090909091)
#define SLL_CONST_1_12 ((sll)0.083333333333333333333333333333333)
#define SLL_CONST_1_20 ((sll)0.05)
#define SLL_CONST_1_30 ((sll)0.033333333333333333333333333333333)
#define SLL_CONST_1_42 ((sll)0.023809523809523809523809523809524)
#define SLL_CONST_1_56 ((sll)0.017857142857142857142857142857143)
#define SLL_CONST_1_72 ((sll)0.013888888888888888888888888888889)
#define SLL_CONST_1_90 ((sll)0.011111111111111111111111111111111)
#define SLL_CONST_1_110 ((sll)0.0090909090909090909090909090909091)
#define SLL_CONST_1_132 ((sll)0.0075757575757575757575757575757576)
#define SLL_CONST_1_156 ((sll)0.0064102564102564102564102564102564)
#define SLL_CONST_E ((sll)2.7182818284590452354)
#define SLL_CONST_1_E ((sll)0.36787944117144232159014926384893)
#define SLL_CONST_1_SQRTE ((sll)0.13533528323661269189004515421424)
#define SLL_CONST_LOG10_E ((sll)0.43429448190325182765747371064254)
#define SLL_CONST_PI ((sll)3.1415926535897932384626433832795)
#define SLL_CONST_2PI ((sll)6.283185307179586476925286766559)
#define SLL_CONST_4PI ((sll)12.566370614359172953850573533118)
#define SLL_CONST_PI_2 ((sll)1.5707963267948966192313216916398)
#define SLL_CONST_PI_4 ((sll)0.78539816339744830961566084581988)
#define SLL_CONST_PI_8 ((sll)0.3926990816987242)
#define SLL_CONST_PI_12 ((sll)0.26179938779914943653855361527329)
#define SLL_CONST_1_PI ((sll)0.31830988618379067153776752674503)
#define SLL_CONST_2_PI ((sll)3.1415926535897932384626433832795)
#define SLL_CONST_ATAN_1 ((sll)0.785398)


#else


/* Data types */
typedef signed long long sll;
typedef unsigned long long ull;

/* Macros */
#define int2sll(X) (((sll) (X)) << 32)
// #define sllvalue(X) (sll2dbl(X))
#define sllvalue(X) (X)
#define sll2int(X) ((int) ((X) >> 32))
#define sll_abs(X) ((X) & 0xefffffffffffffffLL)
#define sllabs(A) (((A)<SLL_CONST_0)? -(A):(A))
#define sllint(X) ((X) & 0xffffffff00000000LL)
#define sllfrac(X) ((X) & 0x00000000ffffffffLL)
#define sllneg(X) (-(X))
#define _slladd(X,Y) ((X) + (Y))
#define _sllsub(X,Y) ((X) - (Y))

#define sllrotl(X,N) ((X)<<N)
#define sllrotr(X,N) ((X)>>N)

/* Constants (converted from double) */
#define SLL_CONST_0 0x0000000000000000LL
#define SLL_CONST_1 0x0000000100000000LL
#define SLL_CONST_2 0x0000000200000000LL
#define SLL_CONST_3 0x0000000300000000LL
#define SLL_CONST_4 0x0000000400000000LL
#define SLL_CONST_10 0x0000000a00000000LL
#define SLL_CONST_15 0x0000000f00000000LL
#define SLL_CONST_1_2 0x0000000080000000LL
#define SLL_CONST_1_3 0x0000000055555555LL
#define SLL_CONST_1_4 0x0000000040000000LL
#define SLL_CONST_1_5 0x0000000033333333LL
#define SLL_CONST_1_6 0x000000002aaaaaaaLL
#define SLL_CONST_1_7 0x0000000024924924LL
#define SLL_CONST_1_8 0x0000000020000000LL
#define SLL_CONST_1_9 0x000000001c71c71cLL
#define SLL_CONST_1_10 0x0000000019999999LL
#define SLL_CONST_1_11 0x000000001745d174LL
#define SLL_CONST_1_12 0x0000000015555555LL
#define SLL_CONST_1_20 0x000000000cccccccLL
#define SLL_CONST_1_30 0x0000000008888888LL
#define SLL_CONST_1_42 0x0000000006186186LL
#define SLL_CONST_1_56 0x0000000004924924LL
#define SLL_CONST_1_72 0x00000000038e38e3LL
#define SLL_CONST_1_90 0x0000000002d82d82LL
#define SLL_CONST_1_110 0x000000000253c825LL
#define SLL_CONST_1_132 0x0000000001f07c1fLL
#define SLL_CONST_1_156 0x0000000001a41a41LL
#define SLL_CONST_E 0x00000002b7e15162LL
#define SLL_CONST_1_E 0x000000005e2d58d8LL
#define SLL_CONST_SQRTE 0x00000001a61298e1LL
#define SLL_CONST_1_SQRTE 0x000000009b4597e3LL
#define SLL_CONST_LOG2_E 0x0000000171547652LL
#define SLL_CONST_LOG10_E 0x000000006f2dec54LL
#define SLL_CONST_LN2 0x00000000b17217f7LL
#define SLL_CONST_LN10 0x000000024d763776LL
#define SLL_CONST_PI 0x00000003243f6a88LL
#define SLL_CONST_2PI 0x00000006487ED510LL
#define SLL_CONST_4PI 0x0000000C90FDAA22LL
#define SLL_CONST_PI_2 0x00000001921fb544LL
#define SLL_CONST_PI_4 0x00000000c90fdaa2LL
#define SLL_CONST_PI_8 0x000000006487ED51LL
#define SLL_CONST_PI_12 0x00000000430548E0LL
#define SLL_CONST_1_PI 0x00000000517cc1b7LL
#define SLL_CONST_2_PI 0x00000000a2f9836eLL
#define SLL_CONST_2_SQRTPI 0x0000000120dd7504LL
#define SLL_CONST_SQRT2 0x000000016a09e667LL
#define SLL_CONST_1_SQRT2 0x00000000b504f333LL
#define SLL_CONST_ATAN_1 0x00000000C90FD7E4LL

static __inline__ double sll2dbl(sll s)
{

union {
double d;
unsigned u[2];
ull _ull;
sll _sll;
} in, retval;
register unsigned exp;
register unsigned flag;

if (s == 0)
return 0.0;

/* Move into memory as args might be passed in regs */
in._sll = s;

/* Handle the negative flag */
if (in._sll < 1) {
flag = 0x80000000;
in._ull = sllneg(in._sll);
} else
flag = 0x00000000;

/* Normalize */
for (exp = 1053; in._ull && (in.u[1] & 0x80000000) == 0; exp--) {
in._ull <<= 1;
}
in._ull <<= 1;
exp++;
in._ull >>= 12;
retval._ull = in._ull;
retval.u[1] |= flag | (exp << 20);

#if defined(__arm__)

/* ARM architecture has a big-endian double */
exp = retval.u[0];
retval.u[0] = retval.u[1];
retval.u[1] = exp;

#endif /* defined(__arm__) */

return retval.d;
}

static __inline__ sll slladd(sll x, sll y)
{
return (x + y);
}

static __inline__ sll sllsub(sll x, sll y)
{
return (x - y);
}

/*
* Let a = A * 2^32 + a_hi * 2^0 + a_lo * 2^(-32)
* Let b = B * 2^32 + b_hi * 2^0 + b_lo * 2^(-32)
*
* Where:
* *_hi is the integer part
* *_lo the fractional part
* A and B are the sign (0 for positive, -1 for negative).
*
* a * b = (A * 2^32 + a_hi * 2^0 + a_lo * 2^-32)
* * (B * 2^32 + b_hi * 2^0 + b_lo * 2^-32)
*
* Expanding the terms, we get:
*
* = A * B * 2^64 + A * b_h * 2^32 + A * b_l * 2^0
* + a_h * B * 2^32 + a_h * b_h * 2^0 + a_h * b_l * 2^-32
* + a_l * B * 2^0 + a_l * b_h * 2^-32 + a_l * b_l * 2^-64
*
* Grouping by powers of 2, we get:
*
* = A * B * 2^64
* Meaningless overflow from sign extension - ignore
*
* + (A * b_h + a_h * B) * 2^32
* Overflow which we can't handle - ignore
*
* + (A * b_l + a_h * b_h + a_l * B) * 2^0
* We only need the low 32 bits of this term, as the rest is overflow
*
* + (a_h * b_l + a_l * b_h) * 2^-32
* We need all 64 bits of this term
*
* + a_l * b_l * 2^-64
* We only need the high 32 bits of this term, as the rest is underflow
*
* Note that:
* a > 0 && b > 0: A = 0, B = 0 and the third term is a_h * b_h
* a < 0 && b > 0: A = -1, B = 0 and the third term is a_h * b_h - b_l
* a > 0 && b < 0: A = 0, B = -1 and the third term is a_h * b_h - a_l
* a < 0 && b < 0: A = -1, B = -1 and the third term is a_h * b_h - a_l - b_l
*/
#if defined(__arm__)
static __inline__ sll sllmul(sll left, sll right)
{
/*
* From gcc/config/arm/arm.h:
* In a pair of registers containing a DI or DF value the 'Q'
* operand returns the register number of the register containing
* the least significant part of the value. The 'R' operand returns
* the register number of the register containing the most
* significant part of the value.
*/
sll retval;

__asm__ (
"@ sllmul\n\t"
"umull %R0, %Q0, %Q1, %Q2\n\t"
"mul %R0, %R1, %R2\n\t"
"umlal %Q0, %R0, %Q1, %R2\n\t"
"umlal %Q0, %R0, %Q2, %R1\n\t"
"tst %R1, #0x80000000\n\t"
"subne %R0, %R0, %Q2\n\t"
"tst %R2, #0x80000000\n\t"
"subne %R0, %R0, %Q1\n\t"
: "=&r" (retval)
: "%r" (left), "r" (right)
: "cc"
);

return retval;
}
#elif defined(__i386__)
static __inline__ sll sllmul(sll left, sll right)
{
register sll retval;
__asm__(
"# sllmul\n\t"
" movl %1, %%eax\n\t"
" mull %3\n\t"
" movl %%edx, %%ebx\n\t"
"\n\t"
" movl %2, %%eax\n\t"
" mull %4\n\t"
" movl %%eax, %%ecx\n\t"
"\n\t"
" movl %1, %%eax\n\t"
" mull %4\n\t"
" addl %%eax, %%ebx\n\t"
" adcl %%edx, %%ecx\n\t"
"\n\t"
" movl %2, %%eax\n\t"
" mull %3\n\t"
" addl %%ebx, %%eax\n\t"
" adcl %%ecx, %%edx\n\t"
"\n\t"
" btl $31, %2\n\t"
" jnc 1f\n\t"
" subl %3, %%edx\n\t"
"1: btl $31, %4\n\t"
" jnc 1f\n\t"
" subl %1, %%edx\n\t"
"1:\n\t"
: "=&A" (retval)
: "m" (left), "m" (((unsigned *) &left)[1]),
"m" (right), "m" (((unsigned *) &right)[1])
: "ebx", "ecx", "cc"
);
return retval;
}
#else
/* Plain C version: not optimal but portable. */
#warning Fixed Point no optimal
static __inline__ sll sllmul(sll a, sll b)
{
unsigned int a_lo, b_lo;
signed int a_hi, b_hi;
sll x;

a_lo = a;
a_hi = (ull) a >> 32;
b_lo = b;
b_hi = (ull) b >> 32;

x = ((ull) (a_hi * b_hi) << 32)
+ (((ull) a_lo * b_lo) >> 32)
+ (sll) a_lo * b_hi
+ (sll) b_lo * a_hi;

return x;
}
#endif

static __inline__ sll sllinv(sll v)
{
int sgn = 0;
sll u;
ull s = -1; //0xFFFFFFFFFFFFFFFF; //-1;

/* Use positive numbers, or the approximation won't work */
if (v < SLL_CONST_0) {
v = sllneg(v);
sgn = 1;
}

/* An approximation - must be larger than the actual value */
for (u = v; u; ((ull)u) >>= 1)
s >>= 1;

/* Newton's Method */
u = sllmul(s, _sllsub(SLL_CONST_2, sllmul(v, s)));
u = sllmul(u, _sllsub(SLL_CONST_2, sllmul(v, u)));
u = sllmul(u, _sllsub(SLL_CONST_2, sllmul(v, u)));
u = sllmul(u, _sllsub(SLL_CONST_2, sllmul(v, u)));
u = sllmul(u, _sllsub(SLL_CONST_2, sllmul(v, u)));
u = sllmul(u, _sllsub(SLL_CONST_2, sllmul(v, u)));

return ((sgn) ? sllneg(u): u);
}

static __inline__ sll slldiv(sll left, sll right)
{
return sllmul(left, sllinv(right));
}

static __inline__ sll sllmul2(sll x)
{
return x << 1;
}

static __inline__ sll sllmul4(sll x)
{
return x << 2;
}

static __inline__ sll sllmul2n(sll x, int n)
{
sll y;

#if defined(__arm__)
/*
* On ARM we need to do explicit assembly since the compiler
* doesn't know the range of n is limited and decides to call
* a library function instead.
*/
__asm__ (
"@ sllmul2n\n\t"
"mov %R0, %R1, lsl %2\n\t"
"orr %R0, %R0, %Q1, lsr %3\n\t"
"mov %Q0, %Q1, lsl %2\n\t"
: "=r" (y)
: "r" (x), "rM" (n), "rM" (32 - n)
);
#else
y = x << n;
#endif

return y;
}

static __inline__ sll slldiv2(sll x)
{
return x >> 1;
}

static __inline__ sll slldiv4(sll x)
{
return x >> 2;
}

static __inline__ sll slldiv2n(sll x, int n)
{
sll y;

#if defined(__arm__)
/*
* On ARM we need to do explicit assembly since the compiler
* doesn't know the range of n is limited and decides to call
* a library function instead.
*/
__asm__ (
"@ slldiv2n\n\t"
"mov %Q0, %Q1, lsr %2\n\t"
"orr %Q0, %Q0, %R1, lsl %3\n\t"
"mov %R0, %R1, asr %2\n\t"
: "=r" (y)
: "r" (x), "rM" (n), "rM" (32 - n)
);
#else
y = x >> n;
#endif

return y;
}

/*
* Unpack the IEEE floating point double format and put it in fixed point
* sll format.
*/
static __inline__ sll dbl2sll(double dbl)
{
union {
double d;
unsigned u[2];
ull _ull;
sll _sll;
} in, retval;
register unsigned exp;

/* Move into memory as args might be passed in regs */
in.d = dbl;

#if defined(__arm__)

/* ARM architecture has a big-endian double */
exp = in.u[0];
in.u[0] = in.u[1];
in.u[1] = exp;

#endif /* defined(__arm__) */

/* Leading 1 is assumed by IEEE */
retval.u[1] = 0x40000000;

/* Unpack the mantissa into the unsigned long */
retval.u[1] |= (in.u[1] << 10) & 0x3ffffc00;
retval.u[1] |= (in.u[0] >> 22) & 0x000003ff;
retval.u[0] = in.u[0] << 10;

/* Extract the exponent and align the decimals */
exp = (in.u[1] >> 20) & 0x7ff;
if (exp)
retval._ull >>= 1053 - exp;
else
return 0L;

/* Negate if negative flag set */
if (in.u[1] & 0x80000000)
retval._sll = -retval._sll;

return retval._sll;
}

static __inline__ sll float2sll(float f)
{
return dbl2sll((double)f);
}

static __inline__ float sll2float(sll s)
{
return ((float)sll2dbl(s));
}

/*
* Calculate cos x where -pi/4 <= x <= pi/4
*
* Description:
* cos x = 1 - x^2 / 2! + x^4 / 4! - ... + x^(2N) / (2N)!
* Note that (pi/4)^12 / 12! < 2^-32 which is the smallest possible number.
*/
static __inline__ sll _sllcos(sll x)
{
sll retval, x2;
x2 = sllmul(x, x);
/*
* cos x = t0 + t1 + t2 + t3 + t4 + t5 + t6
*
* f0 = 0! = 1
* f1 = 2! = 2 * 1 * f0 = 2 * f0
* f2 = 4! = 4 * 3 * f1 = 12 x f1
* f3 = 6! = 6 * 5 * f2 = 30 * f2
* f4 = 8! = 8 * 7 * f3 = 56 * f3
* f5 = 10! = 10 * 9 * f4 = 90 * f4
* f6 = 12! = 12 * 11 * f5 = 132 * f5
*
* t0 = 1
* t1 = -t0 * x2 / 2 = -t0 * x2 * SLL_CONST_1_2
* t2 = -t1 * x2 / 12 = -t1 * x2 * SLL_CONST_1_12
* t3 = -t2 * x2 / 30 = -t2 * x2 * SLL_CONST_1_30
* t4 = -t3 * x2 / 56 = -t3 * x2 * SLL_CONST_1_56
* t5 = -t4 * x2 / 90 = -t4 * x2 * SLL_CONST_1_90
* t6 = -t5 * x2 / 132 = -t5 * x2 * SLL_CONST_1_132
*/
retval = _sllsub(SLL_CONST_1, sllmul(x2, SLL_CONST_1_132));
retval = _sllsub(SLL_CONST_1, sllmul(sllmul(x2, retval), SLL_CONST_1_90));
retval = _sllsub(SLL_CONST_1, sllmul(sllmul(x2, retval), SLL_CONST_1_56));
retval = _sllsub(SLL_CONST_1, sllmul(sllmul(x2, retval), SLL_CONST_1_30));
retval = _sllsub(SLL_CONST_1, sllmul(sllmul(x2, retval), SLL_CONST_1_12));
retval = _sllsub(SLL_CONST_1, slldiv2(sllmul(x2, retval)));
return retval;
}

/*
* Calculate sin x where -pi/4 <= x <= pi/4
*
* Description:
* sin x = x - x^3 / 3! + x^5 / 5! - ... + x^(2N+1) / (2N+1)!
* Note that (pi/4)^13 / 13! < 2^-32 which is the smallest possible number.
*/
static __inline__ sll _sllsin(sll x)
{
sll retval, x2;
x2 = sllmul(x, x);
/*
* sin x = t0 + t1 + t2 + t3 + t4 + t5 + t6
*
* f0 = 0! = 1
* f1 = 3! = 3 * 2 * f0 = 6 * f0
* f2 = 5! = 5 * 4 * f1 = 20 x f1
* f3 = 7! = 7 * 6 * f2 = 42 * f2
* f4 = 9! = 9 * 8 * f3 = 72 * f3
* f5 = 11! = 11 * 10 * f4 = 110 * f4
* f6 = 13! = 13 * 12 * f5 = 156 * f5
*
* t0 = 1
* t1 = -t0 * x2 / 6 = -t0 * x2 * SLL_CONST_1_6
* t2 = -t1 * x2 / 20 = -t1 * x2 * SLL_CONST_1_20
* t3 = -t2 * x2 / 42 = -t2 * x2 * SLL_CONST_1_42
* t4 = -t3 * x2 / 72 = -t3 * x2 * SLL_CONST_1_72
* t5 = -t4 * x2 / 110 = -t4 * x2 * SLL_CONST_1_110
* t6 = -t5 * x2 / 156 = -t5 * x2 * SLL_CONST_1_156
*/
retval = _sllsub(x, sllmul(x2, SLL_CONST_1_156));
retval = _sllsub(x, sllmul(sllmul(x2, retval), SLL_CONST_1_110));
retval = _sllsub(x, sllmul(sllmul(x2, retval), SLL_CONST_1_72));
retval = _sllsub(x, sllmul(sllmul(x2, retval), SLL_CONST_1_42));
retval = _sllsub(x, sllmul(sllmul(x2, retval), SLL_CONST_1_20));
retval = _sllsub(x, sllmul(sllmul(x2, retval), SLL_CONST_1_6));
return retval;
}

static __inline__ sll sllcos(sll x)
{
int i;
sll retval;

/* Calculate cos (x - i * pi/2), where -pi/4 <= x - i * pi/2 <= pi/4 */
i = sll2int(_slladd(sllmul(x, SLL_CONST_2_PI), SLL_CONST_1_2));
x = _sllsub(x, sllmul(int2sll(i), SLL_CONST_PI_2));

switch (i & 3) {
default:
case 0:
retval = _sllcos(x);
break;
case 1:
retval = sllneg(_sllsin(x));
break;
case 2:
retval = sllneg(_sllcos(x));
break;
case 3:
retval = _sllsin(x);
break;
}
return retval;
}

static __inline__ sll sllsin(sll x)
{
int i;
sll retval;

/* Calculate sin (x - n * pi/2), where -pi/4 <= x - i * pi/2 <= pi/4 */
i = sll2int(_slladd(sllmul(x, SLL_CONST_2_PI), SLL_CONST_1_2));
x = _sllsub(x, sllmul(int2sll(i), SLL_CONST_PI_2));

switch (i & 3) {
default:
case 0:
retval = _sllsin(x);
break;
case 1:
retval = _sllcos(x);
break;
case 2:
retval = sllneg(_sllsin(x));
break;
case 3:
retval = sllneg(_sllcos(x));
break;
}
return retval;
}

static __inline__ sll slltan(sll x)
{
int i;
sll retval;

i = sll2int(_slladd(sllmul(x, SLL_CONST_2_PI), SLL_CONST_1_2));
x = _sllsub(x, sllmul(int2sll(i), SLL_CONST_PI_2));
switch (i & 3) {
default:
case 0:
case 2:
retval = slldiv(_sllsin(x), _sllcos(x));
break;
case 1:
case 3:
retval = sllneg(slldiv(_sllcos(x), _sllsin(x)));
break;
}
return retval;
}

/*
* atan x = SUM[n=0,) (-1)^n * x^(2n + 1)/(2n + 1), |x| < 1
*
* Two term approximation
* a = x - x^3/3
* Gives us
* atan x = a + ??
* Let ?? = arctan ?
* atan x = a + arctan ?
* Rearrange
* atan x - a = arctan ?
* Apply tan to both sides
* tan (atan x - a) = tan arctan ?
* tan (atan x - a) = ?
* Applying the standard formula
* tan (u - v) = (tan u - tan v) / (1 + tan u * tan v)
* Gives us
* tan (atan x - a) = (tan atan x - tan a) / (1 + tan arctan x * tan a)
* Let t = tan a
* tan (atan x - a) = (x - t) / (1 + x * t)
* So finally
* arctan x = a + arctan ((tan x - t) / (1 + x * t))
* And the typical worst case is x = 1.0 which converges in 3 iterations.
*/
static __inline__ sll _sllatan(sll x)
{
sll a, t, retval;

/* First iteration */
a = sllmul(x, _sllsub(SLL_CONST_1, sllmul(x, sllmul(x, SLL_CONST_1_3))));
retval = a;

/* Second iteration */
t = slldiv(_sllsin(a), _sllcos(a));
x = slldiv(_sllsub(x, t), _slladd(SLL_CONST_1, sllmul(t, x)));
a = sllmul(x, _sllsub(SLL_CONST_1, sllmul(x, sllmul(x, SLL_CONST_1_3))));
retval = _slladd(retval, a);

/* Third iteration */
t = slldiv(_sllsin(a), _sllcos(a));
x = slldiv(_sllsub(x, t), _slladd(SLL_CONST_1, sllmul(t, x)));
a = sllmul(x, _sllsub(SLL_CONST_1, sllmul(x, sllmul(x, SLL_CONST_1_3))));
return _slladd(retval, a);
}

static __inline__ sll sllatan(sll x)
{
sll retval;

if (x < -sllneg(SLL_CONST_1))
retval = sllneg(SLL_CONST_PI_2);
else if (x > SLL_CONST_1)
retval = SLL_CONST_PI_2;
else
return _sllatan(x);
return _sllsub(retval, _sllatan(sllinv(x)));
}

/*
* Calculate e^x where -0.5 <= x <= 0.5
*
* Description:
* e^x = x^0 / 0! + x^1 / 1! + ... + x^N / N!
* Note that 0.5^11 / 11! < 2^-32 which is the smallest possible number.
*/
static __inline__ sll _sllexp(sll x)
{
sll retval;
retval = _slladd(SLL_CONST_1, sllmul(0, sllmul(x, SLL_CONST_1_11)));
retval = _slladd(SLL_CONST_1, sllmul(retval, sllmul(x, SLL_CONST_1_11)));
retval = _slladd(SLL_CONST_1, sllmul(retval, sllmul(x, SLL_CONST_1_10)));
retval = _slladd(SLL_CONST_1, sllmul(retval, sllmul(x, SLL_CONST_1_9)));
retval = _slladd(SLL_CONST_1, sllmul(retval, slldiv2n(x, 3)));
retval = _slladd(SLL_CONST_1, sllmul(retval, sllmul(x, SLL_CONST_1_7)));
retval = _slladd(SLL_CONST_1, sllmul(retval, sllmul(x, SLL_CONST_1_6)));
retval = _slladd(SLL_CONST_1, sllmul(retval, sllmul(x, SLL_CONST_1_5)));
retval = _slladd(SLL_CONST_1, sllmul(retval, slldiv4(x)));
retval = _slladd(SLL_CONST_1, sllmul(retval, sllmul(x, SLL_CONST_1_3)));
retval = _slladd(SLL_CONST_1, sllmul(retval, slldiv2(x)));
return retval;
}

/*
* Calculate e^x where x is arbitrary
*/
static __inline__ sll sllexp(sll x)
{
int i;
sll e, retval;

e = SLL_CONST_E;

/* -0.5 <= x <= 0.5 */
i = sll2int(_slladd(x, SLL_CONST_1_2));
retval = _sllexp(_sllsub(x, int2sll(i)));

/* i >= 0 */
if (i < 0) {
i = -i;
e = SLL_CONST_1_E;
}

/* Scale the result */
for (;i; i >>= 1) {
if (i & 1)
retval = sllmul(retval, e);
e = sllmul(e, e);
}
return retval;
}

/*
* Calculate natural logarithm using Netwton-Raphson method
*/
static __inline__ sll slllog(sll x)
{
sll x1, ln = 0;

/* Scale: e^(-1/2) <= x <= e^(1/2) */
while (x < SLL_CONST_1_SQRTE) {
ln = _sllsub(ln, SLL_CONST_1);
x = sllmul(x, SLL_CONST_E);
}
while (x > SLL_CONST_SQRTE) {
ln = _slladd(ln, SLL_CONST_1);
x = sllmul(x, SLL_CONST_1_E);
}

/* First iteration */
x1 = sllmul(_sllsub(x, SLL_CONST_1), slldiv2(_sllsub(x, SLL_CONST_3)));
ln = _sllsub(ln, x1);
x = sllmul(x, _sllexp(x1));

/* Second iteration */
x1 = sllmul(_sllsub(x, SLL_CONST_1), slldiv2(_sllsub(x, SLL_CONST_3)));
ln = _sllsub(ln, x1);
x = sllmul(x, _sllexp(x1));

/* Third iteration */
x1 = sllmul(_sllsub(x, SLL_CONST_1), slldiv2(_sllsub(x, SLL_CONST_3)));
ln = _sllsub(ln, x1);

return ln;
}

/*
* ln x^y = y * log x
* e^(ln x^y) = e^(y * log x)
* x^y = e^(y * ln x)
*/
static __inline__ sll sllpow(sll x, sll y)
{
if (y == SLL_CONST_0)
return SLL_CONST_1;
return sllexp(sllmul(y, slllog(x)));
}

/*
* Consider a parabola centered on the y-axis
* y = a * x^2 + b
* Has zeros (y = 0) at
* a * x^2 + b = 0
* a * x^2 = -b
* x^2 = -b / a
* x = +- (-b / a)^(1 / 2)
* Letting a = 1 and b = -X
* y = x^2 - X
* x = +- X^(1 / 2)
* Which is convenient since we want to find the square root of X, and we can
* use Newton's Method to find the zeros of any f(x)
* xn = x - f(x) / f'(x)
* Applied Newton's Method to our parabola
* f(x) = x^2 - X
* xn = x - (x^2 - X) / (2 * x)
* xn = x - (x - X / x) / 2
* To make this converge quickly, we scale X so that
* X = 4^N * z
* Taking the roots of both sides
* X^(1 / 2) = (4^n * z)^(1 / 2)
* X^(1 / 2) = 2^n * z^(1 / 2)
* Let N = 2^n
* x^(1 / 2) = N * z^(1 / 2)
* We want this to converge to the positive root, so we must start at a point
* 0 < start <= x^(1 / 2)
* or
* x^(1/2) <= start <= infinity
* since
* (1/2)^(1/2) = 0.707
* 2^(1/2) = 1.414
* A good choice is 1 which lies in the middle, and takes 4 iterations to
* converge from either extreme.
*/
static __inline__ sll sllsqrt(sll x)
{
sll n, xn;

/* Start with a scaling factor of 1 */
n = SLL_CONST_1;

/* Quick solutions for the simple cases */
if (x <= SLL_CONST_0 || x == SLL_CONST_1)
return x;

/* Scale x so that 0.5 <= x < 2 */
while (x >= SLL_CONST_2) {
x = slldiv4(x);
n = sllmul2(n);
}
while (x < SLL_CONST_1_2) {
x = sllmul4(x);
n = slldiv2(n);
}

/* Simple solution if x = 4^n */
if (x == SLL_CONST_1)
return n;

/* The starting point */
xn = SLL_CONST_1;

/* Four iterations will be enough */
xn = _sllsub(xn, slldiv2(_sllsub(xn, slldiv(x, xn))));
xn = _sllsub(xn, slldiv2(_sllsub(xn, slldiv(x, xn))));
xn = _sllsub(xn, slldiv2(_sllsub(xn, slldiv(x, xn))));
xn = _sllsub(xn, slldiv2(_sllsub(xn, slldiv(x, xn))));

/* Scale the result */
return sllmul(n, xn);
}


#endif

static inline sll sllatan2(sll x, sll y)
{
register sll ret;
register sll d=slldiv(y,x);

if (d <= SLL_CONST_1 && d >= sllneg(SLL_CONST_1))
ret=sllmul(SLL_CONST_ATAN_1, d);
else
ret=sllatan(d);

if (x > SLL_CONST_0)
ret= slladd(ret, SLL_CONST_PI);

return ret;
}


#ifdef __cplusplus
}
#endif

#endif /* MATHSLL_H */
 
For not being portable, you sure made an effort to make it portable.
31.32 is probably far over the top for most games.
 
Thanks for the advice- I already had a quick look at math-sll but I'll try and use that in my code... I'll post up the results when I'm done. :)
 
Back
Top