/*
** Command & Conquer Renegade(tm)
** Copyright 2025 Electronic Arts Inc.
**
** 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, either version 3 of the License, or
** (at your option) any later version.
**
** This program 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 General Public License for more details.
**
** You should have received a copy of the GNU General Public License
** along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
/*
This file contains some math functions that Jan worked on while visiting westwood
between Meltdown and Siggraph 2001.
*/
#include <iostream>
#include <iomanip>
#include <cmath>
#include "p_timer.h"
#include "wwmath.h"
using namespace std;
#undef for
#define for if(false); else for
/* (1-u)(1-u)(1-u)
--------------- 1/6 0/6
6
3*u*u*u - 6u*u +4
--------------- 4/6 1/6
6
-3*u*u*u + 3*u*u + 3*u + 1
--------------- 1/6 4/6
6
u*u*u
--------------- 0/6 1/6
6
*/
const double pi = 3.1415926535897932384626433832795;
const int SINTABLESIZE = 16;
float sinTable[SINTABLESIZE+1][3];
/*
float sinTable[SINTABLESIZE+1][3] =
{
0.f, 0.0327249329938f, 0.0654498139838f,
0.0980171403296f, 0.130584524375f, 0.16299423692f,
0.195090322016f, 0.22718650856f, 0.258968953094f,
0.290284677254f, 0.321600570763f, 0.352449632798f,
0.382683432365f, 0.412917427848f, 0.442536062284f,
0.471396736826f, 0.5002577f, 0.528360611467f,
0.55557023302f, 0.582780177536f, 0.609096734804f,
0.634393284164f, 0.659690191419f, 0.683966978964f,
0.707106781187f, 0.730247032377f, 0.752250161646f,
0.773010453363f, 0.79377112016f, 0.813288873188f,
0.831469612303f, 0.849650815039f, 0.866495135964f,
0.881921264348f, 0.897347904711f, 0.911356511164f,
0.923879532511f, 0.936403079986f, 0.947441092993f,
0.956940335732f, 0.966440181792f, 0.974401217255f,
0.980785280403f, 0.987169903631f, 0.991977366409f,
0.995184726672f, 0.998392578191f, 1.0000003005f,
1.f, 1.0000003005f, 0.998392578191f,
};
*/
void init_bez_sin()
{
for(int i=0; i<SINTABLESIZE+1; i++)
{
double x = sin(pi/2*(i*3)/(SINTABLESIZE*3));
double y = sin(pi/2*(i*3+1)/(SINTABLESIZE*3)); // * 1.00053f;
double z = sin(pi/2*(i*3+2)/(SINTABLESIZE*3)); // * 1.00053f;
double w = sin(pi/2*(i*3+3)/(SINTABLESIZE*3));
double mse = 4.0;
double my = 1.0;
double mz = 1.0;
// * nu2 * nu;
// nu2;
// u * nu;
// * u2 * u;
double wy = 4.0;
for(double dy=1.000178; dy<1.000537; dy+=0.0000001)
{
bool b = false;
double wz = 4.0;
for(double dz=1.000534; dz<1.000715; dz+=0.0000001)
{
double by = x * 8.0 / 27.0 +
y * dy * 12.0 / 27.0 +
z * dz * 6.0 / 27.0 +
w / 27.0;
double bz = x / 27.0 +
y * dy * 6.0 / 27.0 +
z * dz * 12.0 / 27.0 +
w * 8.0 / 27.0;
if(y*dy < 0.0 || y*dy > 1.0)
continue;
if(z*dz < 0.0 || z*dz > 1.0)
continue;
double yerr, zerr;
// if(by > y)
// yerr = by/y;
// else
yerr = y/by;
// if(bz > z)
// zerr = bz/z;
// else
zerr = z/bz;
yerr -= 1.0;
zerr -= 1.0;
double se = yerr*yerr + zerr*zerr;
if(mse > se)
{
wy = se;
wz = se;
mse = se;
my = dy;
mz = dz;
b = false;
}
else
{
// if(se > wy)
// b = true;
// if(se > wz)
// break;
}
}
// if(b)
// break;
}
sinTable[i][0] = float(x);
sinTable[i][1] = float(y * my);
sinTable[i][2] = float(z * mz);
cout << setprecision(12);
cout << setw(16) << x << "f, " << setw(16) << y*my << "f, " << setw(16) << z*mz << "f," << endl;
}
}
const int ACOSTABLESIZE = 32;
//float acosTable[ACOSTABLESIZE+1][3];
float acosTable[ACOSTABLESIZE+1][3] =
{
0, 0, 0,
0, 0, 0,
0, 0, 0,
0, 0, 0,
0, 0, 0,
0, 0, 0,
0, 0, 0,
0, 0, 0,
2.09439510239f, 2.07033956061f, 2.04678164613f,
2.02361292154f, 2.00044498842f, 1.97766674853f,
1.95519310129f, 1.93271994669f, 1.91055213001f,
1.88862003072f, 1.86668845355f, 1.84499333386f,
1.82347658194f, 1.80196011523f, 1.7806223216f,
1.75941271297f, 1.73820326172f, 1.71712229397f,
1.69612415796f, 1.67512622765f, 1.65421121018f,
1.63333708859f, 1.61246303408f, 1.59162963139f,
1.57079632679f, 1.54996304584f, 1.52912960772f,
1.508255565f, 1.48738138962f, 1.46646655946f,
1.44546849563f, 1.42447034954f, 1.40338938542f,
1.38217994062f, 1.36097042076f, 1.33963248005f,
1.31811607165f, 1.29659929615f, 1.27490425438f,
1.25297262287f, 1.23104063033f, 1.20887262684f,
1.1863995523f, 1.16392587108f, 1.14114765337f,
1.11797973205f, 1.09481095108f, 1.07125317716f,
1.0471975512f, 1.0231405f, 0.998586616212f,
0, 0, 0,
0, 0, 0,
0, 0, 0,
0, 0, 0,
0, 0, 0,
0, 0, 0,
0, 0, 0,
0, 0, 0
};
double calc_cos_plus_45(double cx)
{
double sx = sqrt(1-cx*cx);
cx = cx*cos(pi/4) - sx*sin(pi/4);
return cx;
}
void init_bez_acos()
{
for(int i=0; i<ACOSTABLESIZE; i++)
{
// if(i < ACOSTABLESIZE/4 || i > ACOSTABLESIZE/4*3)
// continue;
double x = acos(float(i*3)/(ACOSTABLESIZE*3)*2-1);
double y = acos(float(i*3+1)/(ACOSTABLESIZE*3)*2-1); // * 1.00053f;
double z = acos(float(i*3+2)/(ACOSTABLESIZE*3)*2-1); // * 1.00053f;
double w;
if(i != ACOSTABLESIZE-1)
w = acos(float(i*3+3)/(ACOSTABLESIZE*3)*2-1);
else
w = 0.0;//-acos(float(i*3+1)/(ACOSTABLESIZE*3)*2-1);
double mse = 4.0;
double my = 1.0;
double mz = 1.0;
double wy = 4.0;
for(double dy=0.9996944; dy<=1.0006926; dy+=16.0/8388608)
{
bool b = false;
double wz = 4.0;
for(double dz=0.9997872; dz<=1.0011092; dz+=16.0/8388608)
{
double by = x * 8.0 / 27.0 +
y * dy * 12.0 / 27.0 +
z * dz * 6.0 / 27.0 +
w / 27.0;
double bz = x / 27.0 +
y * dy * 6.0 / 27.0 +
z * dz * 12.0 / 27.0 +
w * 8.0 / 27.0;
if(y*dy < 0.0 || y*dy > pi)
continue;
if(z*dz < 0.0 || z*dz > pi)
continue;
double yerr, zerr;
// if(by > y)
// yerr = by/y;
// else
yerr = y/by;
// if(bz > z)
// zerr = bz/z;
// else
zerr = z/bz;
yerr -= 1.0;
zerr -= 1.0;
double se = yerr*yerr + zerr*zerr;
if(mse > se)
{
wy = se;
wz = se;
mse = se;
my = dy;
mz = dz;
b = false;
}
else
{
// if(se > wy)
// b = true;
// if(se > wz)
// break;
}
}
// if(b)
// break;
}
acosTable[i][0] = float(x);
acosTable[i][1] = float(y * my);
acosTable[i][2] = float(z * mz);
cout << setprecision(12);
cout << setw(16) << x << "f, " << setw(16) << y*my << "f, " << setw(16) << z*mz << "f," << endl;
// cout << my << " " << mz << endl;
}
acosTable[ACOSTABLESIZE][0] = 0.f;
}
float bez_sin(float x)
{
x *= (SINTABLESIZE*4)/(pi*2);
int ix0 = int(floor(x));
float u = x - ix0;
float nu = 1.f-u;
float sign = 1.f;
ix0 &= SINTABLESIZE*4-1;
if(ix0 >= SINTABLESIZE*2)
{
sign = -1.f;
ix0 -= SINTABLESIZE*2;
}
if(ix0 >= SINTABLESIZE)
{
ix0 = (SINTABLESIZE-1)-(ix0-SINTABLESIZE);
float t = u;
u = nu;
nu = t;
}
int ix1 = ix0+1;
// ix0 &= SINTABLESIZE-1;
// ix1 &= SINTABLESIZE-1;
float u2 = u*u;
float nu2 = nu*nu;
float r0 = sinTable[ix0][0] * nu2 * nu;
float r1 = sinTable[ix0][1] * 3 * u * nu2;
float r2 = sinTable[ix0][2] * 3 * u2 * nu;
float r3 = sinTable[ix1][0] * u2 * u;
r0 += r1;
r2 += r3;
return (r0+r2) * sign;
}
float bez_acos(float x)
{
// double sx = sqrt(1-x*x);
// x = x*cos(pi/4) - sx*sin(pi/4);
x = calc_cos_plus_45(x);
x += 1.f;
x *= ACOSTABLESIZE/2;
int ix0 = int(floor(x));
float u = x - ix0;
float nu = 1.f-u;
if(ix0 < 0)
return pi-pi/4; // -infinite...
if(ix0 >= ACOSTABLESIZE)
return 0.f;
int ix1 = ix0+1;
float u2 = u*u;
float nu2 = nu*nu;
float r0 = acosTable[ix0][0] * nu2 * nu;
float r1 = acosTable[ix0][1] * 3 * u * nu2;
float r2 = acosTable[ix0][2] * 3 * u2 * nu;
float r3 = acosTable[ix1][0] * u2 * u;
r0 += r1;
r2 += r3;
return (r0+r2)-pi/4;
}
float mcos(float c, float m)
{
float ac = acos(c);
float c1 = cos(ac*m);
float c2 = cos(ac*(1-m));
float c3 = c * cos(m);
return c1;
}
float cheb_asin(float x)
{
float table[] =
{
1.101460576f,
+.5764093225f,
-.09735417608f,
+.003083774398f,
-.000844493744f,
+.005722363752f,
-.005962135984f
};
float r = table[0] + table[1]*x;
r -= sqrt(1-x*x);
x = 2.f*x-1.f;
float x2 = x;
x *= x;
for(int i=2; i<sizeof(table)/sizeof(table[0]); i++)
{
r += x * table[i];
x *= x2;
}
return r;
}
float taylor_acos(float x)
{
/*
(1) /(2*3),
(3) /(2*4*5),
(3*5) /(2*4*6*7),
(3*5*7) /(2*4*6*8*9),
(3*5*7*9) /(2*4*6*8*10*11),
(3*5*7*9*11) /(2*4*6*8*10*12*13),
(3*5*7*9*11*13) /(2*4*6*8*10*12*14*15),
(3*5*7*9*11*13*15) /(2*4*6*8*10*12*14*16*17),
(3*5*7*9*11*13*15*17) /(2*4*6*8*10*12*14*16*18*19),
*/
float table[] =
{
1.f,
float(1) /(2*3),
float(3) /(2*4*5),
float(5) /(2*4*2*7),
float(5*7) /(2*4*2*8*9),
float(7*9) /(2*4*2*8*2*11),
float(7*9*11) /(2*4*2*8*2*12*13),
float(9*11*13) /(2*4*2*8*2*12*2*15),
float(9*11*13*15) /(2*4*2*8*2*12*2*16*17),
float(11*13*15*17) /(2*4*2*8*2*12*2*16*2*19),
};
// x = sqrt((1.f+x)/2.f);
bool flip = false;
if(x > 0.707106781186547524400844362104849f)
{
x = 0.707106781186547524400844362104849f - (x-0.707106781186547524400844362104849f);
flip = true;
}
float r = 1.57079632679489661923132169163975f;
float x2 = x*x;
for(int i=0; i<sizeof(table)/sizeof(table[0]); i++)
{
r -= x * table[i];
x *= x2;
}
if(flip)
{
r = 1.57079632679489661923132169163975f-r;
}
return r;
}
/*
taylor acos
3 5 7 9 11
%PI X 3 X 5 X 35 X 63 X
--- - X - -- - ---- - ---- - ----- - ------
2 6 40 112 1152 2816
taylor(sin(x),x,0,9);
3 5 7 9
X X X X
X - -- + --- - ---- + ------ + . . .
3! 5! 7! 9!
taylor cos
2 4 6 8
X X X X
1 - -- + -- - --- + ----- + . . .
2! 4! 6! 8!
*/
int intChop(const float& f)
{
int a = *reinterpret_cast<const int*>(&f); // take bit pattern of float into a register
int sign = (a>>31); // sign = 0xFFFFFFFF if original value is negative, 0 if positive
int mantissa = (a&((1<<23)-1))|(1<<23); // extract mantissa and add the hidden bit
int exponent = ((a&0x7fffffff)>>23)-127; // extract the exponent
int r = (unsigned int(mantissa)<<8)>>(31-exponent); // ((1<<exponent)*mantissa)>>24 -- (we know that mantissa > (1<<24))
return ((r ^ (sign)) - sign ) &~ (exponent>>31); // add original sign. If exponent was negative, make return value 0.
}
int intFloor (const float& f)
{
int a = *reinterpret_cast<const int*>(&f); // take bit pattern of float into a register
int sign = (a>>31); // sign = 0xFFFFFFFF if original value is negative, 0 if positive
a&=0x7fffffff; // we don't need the sign any more
int exponent = (a>>23)-127; // extract the exponent
int expsign = ~(exponent>>31); // 0xFFFFFFFF if exponent is positive, 0 otherwise
int imask = ( (1<<(31-(exponent))))-1; // mask for true integer values
int mantissa = (a&((1<<23)-1)); // extract mantissa (without the hidden bit)
int r = (unsigned int(mantissa|(1<<23))<<8)>>(31-exponent); // ((1<<exponent)*(mantissa|hidden bit))>>24 -- (we know that mantissa > (1<<24))
r = ((r & expsign) ^ (sign)) + ((!((mantissa<<8)&imask)&(expsign^((a-1)>>31)))&sign); // if (fabs(value)<1.0) value = 0; copy sign; if (value < 0 && value==(int)(value)) value++;
return r;
}
const int CACHE_TRASH_BUFFER_SIZE = 4024000;
int * CacheTrashBuffer = NULL;
volatile int ReadVar;
void trash_the_cache(void)
{
// read a million random bytes and overwrite them with a new value
for (int j=0; j<1024000; j++) {
int address = rand() % CACHE_TRASH_BUFFER_SIZE;
ReadVar = CacheTrashBuffer[address];
CacheTrashBuffer[address] = ReadVar+1;
//cout<<ReadVar;
}
}
int jan_main(int argc, char* argv[])
{
// init_bez_sin();
init_bez_acos();
CacheTrashBuffer = new int[CACHE_TRASH_BUFFER_SIZE];
for (int i=0; i<CACHE_TRASH_BUFFER_SIZE; i++) {
CacheTrashBuffer[i] = rand();
}
for(int i=0; i<10; i++)
{
// float foo = i/float(47) * pi * 2;
float foo = float(i)/64.f;
float r0 = 0.0f;
float r1 = 0.0f;
float r2 = 0.0f;
unsigned long acos_cycles = 0;
unsigned long bez_cycles = 0;
unsigned long table_cycles = 0;
unsigned long acos_sum = 0;
unsigned long bez_sum = 0;
unsigned long table_sum = 0;
const int SAMPLE_COUNT = 20;
{
for (int j=0; j<SAMPLE_COUNT; j++) {
foo = WWMath::Random_Float();
acos_cycles = Get_CPU_Clock();
r0 = acos(foo);
acos_sum += Get_CPU_Clock() - acos_cycles;
trash_the_cache();
}
}
{
for (int j=0; j<SAMPLE_COUNT; j++) {
foo = WWMath::Random_Float();
bez_cycles = Get_CPU_Clock();
r1 = bez_acos(foo);
bez_sum += Get_CPU_Clock() - bez_cycles;
trash_the_cache();
}
}
{
for (int j=0; j<SAMPLE_COUNT; j++) {
foo = WWMath::Random_Float();
table_cycles = Get_CPU_Clock();
r2 = WWMath::Fast_Acos(foo);
table_sum += Get_CPU_Clock() - table_cycles;
trash_the_cache();
}
}
// cout << "x: " << setw(8) << foo << " sin(x): " << setw(8) << r0 << " bez_sin(x): " << setw(8) << r1 << " ratio: " << setw(8) << r0-r1 << endl;
// cout << "x: " << setw(8) << foo << " acos(x): " << setw(8) << r0 << " bez_acos(x): " << setw(8) << r1 << " ratio: " << setw(8) << r0/r1 << endl;
cout << "acos clocks: "<<acos_sum<<" bez clocks: "<<bez_sum<<" table clocks: "<<table_sum << endl;
}
delete CacheTrashBuffer;
return 0;
}