// Shave and a Haircut // (c) 2019 Epic Games // US Patent 6720962 typedef float Vector[3]; typedef float Vector4[4]; /*N*/ static float vcl( Vector v ); /*N*/ static void vnml( Vector out, Vector in ); /*N*/ static void vcross( Vector out, Vector a, Vector c ); /*N*/ static void mcp( Matrix matOut, /* matrix (destination) */ Matrix matIn ); /*N*/ /*N*/ /*N*/ static void matinit( Matrix mat /* matrix to initialize (in-out) */ ); /*N*/ static void matrot( Matrix mat, /* resulting rotation matrix (out) */ VERT vrot /* rotation vector (in) */ ); static void rotx( Matrix mat, /* resulting rotation matrix (out) */ float vrot /* rotation vector (in) */ ); static void roty( Matrix mat, /* resulting rotation matrix (out) */ float vrot /* rotation vector (in) */ ); static void rotz( Matrix mat, /* resulting rotation matrix (out) */ float vrot /* rotation vector (in) */ ); /* vector_stuff.c a collection of vector manipulation routines (some adapted from Gems II) 1993 joe alter */ static VERT Vcross( VERT a, VERT b ) { VERT c; c.x = ( a.y * b.z ) - ( a.z * b.y ); c.y = ( a.z * b.x ) - ( a.x * b.z ); c.z = ( a.x * b.y ) - ( a.y * b.x ); return ( c ); } static void MatMul( Matrix a, Matrix b, Matrix c ) { int i, j, k; for( i = 0; i < 4; i++ ) { for( j = 0; j < 4; j++ ) { c[i][j] = 0; for( k = 0; k < 4; k++ ) c[i][j] += a[i][k] * b[k][j]; } } // return(c); } static float vsql( Vector v ) { return ( ( v[0] * v[0] ) + ( v[1] * v[1] ) + ( v[2] * v[2] ) ); } static float vcl( Vector v ) { return ( sqrt( vsql( v ) ) ); } static void vcp( Vector outVector, Vector inVector ) { outVector[0] = inVector[0]; outVector[1] = inVector[1]; outVector[2] = inVector[2]; } /* vcross product c = a cross b */ static void vcross( Vector out, Vector a, Vector b ) { out[0] = ( a[1] * b[2] ) - ( a[2] * b[1] ); out[1] = ( a[2] * b[0] ) - ( a[0] * b[2] ); out[2] = ( a[0] * b[1] ) - ( a[1] * b[0] ); } static void vnml( Vector out, Vector in ) { float len; vcp( out, in ); len = vcl( out ); if( len != 0.0 ) { out[0] /= len; out[1] /= len; out[2] /= len; } } static void mcp( Matrix matOut, /* matrix (destination) */ Matrix matIn ) { int i, j; for( i = 0; i < 4; i++ ) { for( j = 0; j < 4; j++ ) { matOut[i][j] = matIn[i][j]; } } } /* vxm multiply the vector 'v times the matrix m, result in vector r */ static void mkmatrix( Matrix m, VERT xx, VERT yy, VERT zz, VERT b ) { // xx.x+=.00000001; // yy.y+=.00000001; // zz.z+=.00000001; //#ifdef crap if( ( xx.x == 0 ) && ( xx.y == 0 ) && ( xx.z == 0 ) ) { xx.x = 1.0f; xx.y = 0.0f; xx.z = 0.0f; } if( ( yy.x == 0 ) && ( yy.y == 0 ) && ( yy.z == 0 ) ) { yy.x = 0.0f; yy.y = 1.0f; yy.z = 0.0f; } if( ( zz.x == 0 ) && ( zz.y == 0 ) && ( zz.z == 0 ) ) { zz.x = 0.0f; zz.y = 0.0f; zz.z = 1.0f; } xx = Vnorm( xx ); yy = Vnorm( yy ); zz = Vnorm( zz ); //#endif //xx=Vcross(zz,yy); //yy=Vcross(zz,xx); m[0][0] = xx.x; m[0][1] = xx.y; m[0][2] = xx.z; m[1][0] = yy.x; m[1][1] = yy.y; m[1][2] = yy.z; // m[2][0] = zz.x; m[2][1] = zz.y; m[2][2] = zz.z; m[0][3] = 1.0f; m[1][3] = 1.0f; m[2][3] = 1.0f; m[3][3] = 1.0f; m[3][0] = 0.0f; m[3][1] = 0.0f; m[3][2] = 0.0f; } static void mk_polymat( Matrix m, int n ) { VERT xv, yv, zv, junk; junk.x = 0; junk.y = 0; junk.z = 0; xv.x = v[facelist[face_start[n] + 1]].x - v[facelist[face_start[n]]].x; xv.y = v[facelist[face_start[n] + 1]].y - v[facelist[face_start[n]]].y; xv.z = v[facelist[face_start[n] + 1]].z - v[facelist[face_start[n]]].z; zv.x = v[facelist[face_end[n] - 1]].x - v[facelist[face_start[n]]].x; zv.y = v[facelist[face_end[n] - 1]].y - v[facelist[face_start[n]]].y; zv.z = v[facelist[face_end[n] - 1]].z - v[facelist[face_start[n]]].z; xv = Vnorm( xv ); zv = Vnorm( zv ); yv = Vcross( xv, zv ); zv = Vcross (yv, xv); mkmatrix( m, xv, yv, zv, junk ); } static void mk_polymat2( Matrix m, int n,VERT handle ) { VERT xv, yv, zv, junk; junk.x = 0; junk.y = 0; junk.z = 0; xv.x = v[facelist[face_start[n] + 1]].x - v[facelist[face_start[n]]].x; xv.y = v[facelist[face_start[n] + 1]].y - v[facelist[face_start[n]]].y; xv.z = v[facelist[face_start[n] + 1]].z - v[facelist[face_start[n]]].z; zv.x = v[facelist[face_end[n] - 1]].x - v[facelist[face_start[n]]].x; zv.y = v[facelist[face_end[n] - 1]].y - v[facelist[face_start[n]]].y; zv.z = v[facelist[face_end[n] - 1]].z - v[facelist[face_start[n]]].z; xv=handle; xv = Vnorm( xv ); zv = Vnorm( zv ); yv = Vcross( xv, zv ); zv = Vcross (yv, xv); mkmatrix( m, xv, yv, zv, junk ); } static void Smk_polymat( Matrix m, int n ) { VERT xv, yv, zv, junk; junk.x = 0; junk.y = 0; junk.z = 0; xv.x = Shair[n + 1].hv[0].x - Shair[n].hv[0].x; xv.y = Shair[n + 1].hv[0].y - Shair[n].hv[0].y; xv.z = Shair[n + 1].hv[0].z - Shair[n].hv[0].z; yv.x = Shair[n].hv[1].x - Shair[n].hv[0].x; yv.y = Shair[n].hv[1].y - Shair[n].hv[0].y; yv.z = Shair[n].hv[1].z - Shair[n].hv[0].z; xv = Vnorm( xv ); yv = Vnorm( yv ); zv = Vcross( xv, yv ); zv = Vnorm( zv ); mkmatrix( m, xv, yv, zv, junk ); } static void mkmatrix2( Matrix m, VERT xx, VERT yy, VERT zz, VERT b ) { int x,y; // xx.x+=.00000001; // yy.y+=.00000001; // zz.z+=.00000001; //xx=Vnorm(xx); //yy=Vnorm(yy); //zz=Vnorm(zz); //xx=Vcross(zz,yy); //yy=Vcross(zz,xx); for (x=0;x<4;x++) for (y=0;y<4;y++) m[x][y]=0.0f; m[0][0] = xx.x; m[0][1] = xx.y; m[0][2] = xx.z; m[1][0] = yy.x; m[1][1] = yy.y; m[1][2] = yy.z; // m[2][0] = zz.x; m[2][1] = zz.y; m[2][2] = zz.z; m[0][3] = 1.0f; m[1][3] = 1.0f; m[2][3] = 1.0f; m[3][3] = 1.0f; m[3][0] = 0.0f; m[3][1] = 0.0f; m[3][2] = 0.0f; // m[0][3]=0; // m[1][3]=0; // m[2][3]=0; // m[3][3]=1.0f; // m[3][0]=b.x; // m[3][1]=b.y; // m[3][2]=b.z; } static VERT vxm( Matrix m, /* input matrix (in) */ VERT v /* input vector (in) */ ) { VERT r; r.x = v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0] + m[3][0]; r.y = v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1] + m[3][1]; r.z = v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2] + m[3][2]; return ( r ); } static VERT vxm3x3( Matrix m, /* input matrix (in) */ VERT v /* input vector (in) */ ) { VERT r; r.x = v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0] ; r.y = v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1] ; r.z = v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2] ; return ( r ); } static void printmatrix( Matrix m ) { printf( "%f %f %f %f\n", m[0][0], m[1][0], m[2][0], m[3][0] ); printf( "%f %f %f %f\n", m[0][1], m[1][1], m[2][1], m[3][1] ); printf( "%f %f %f %f\n", m[0][2], m[1][2], m[2][2], m[3][2] ); printf( "%f %f %f %f\n", m[0][3], m[1][3], m[2][3], m[3][3] ); } static void rotz( Matrix m, /* resulting rotation matrix (out) */ float th /* rotation vector (in) */ ) { m[0][0] = cos( th ); m[1][0] = -sin( th ); m[2][0] = 0.0f; m[3][0] = 0.0f; m[0][1] = sin( th ); m[1][1] = cos( th ); m[2][1] = 0.0f; m[3][1] = 0.0f; m[0][2] = 0.0f; m[1][2] = 0.0f; m[2][2] = 1.0f; m[3][2] = 0.0f; m[0][3] = 0.0f; m[1][3] = 0.0f; m[2][3] = 0.0f; m[3][3] = 1.0f; } static void roty( Matrix m, /* resulting rotation matrix (out) */ float th /* rotation vector (in) */ ) { m[0][0] = cos( th ); m[1][0] = 0.0f; m[2][0] = sin( th ); m[3][0] = 0.0f; m[0][1] = 0.0f; m[1][1] = 1.0f; m[2][1] = 0.0f; m[3][1] = 0.0f; m[0][2] = -sin( th ); m[1][2] = 0; m[2][2] = cos( th ); m[3][2] = 0.0f; m[0][3] = 0.0f; m[1][3] = 0.0f; m[2][3] = 0.0f; m[3][3] = 1.0f; } static void rotx( Matrix m, /* resulting rotation matrix (out) */ float th /* rotation vector (in) */ ) { m[0][0] = 1.0f; m[1][0] = 0.0f; m[2][0] = 0.0f; m[3][0] = 0.0f; m[0][1] = 0.0f; m[1][1] = cos( th ); m[2][1] = -sin( th ); m[3][1] = 0.0f; m[0][2] = 0.0f; m[1][2] = sin( th ); m[2][2] = cos( th ); m[3][2] = 0.0f; m[0][3] = 0.0f; m[1][3] = 0.0f; m[2][3] = 0.0f; m[3][3] = 1.0f; } /* matrot builds the roation matrix mat given the SI vrot rotation vector */ static void matrot( Matrix mat, /* resulting rotation matrix (out) */ VERT vrot /* rotation vector (in) */ ) { float sx, cx, sy, cy, sz, cz, tmp; VERT tv; sx = sin( vrot.x ); cx = cos( vrot.x ); sy = sin( vrot.y ); cy = cos( vrot.y ); sz = sin( vrot.z ); cz = cos( vrot.z ); mat[0][0] = cy * cz; mat[0][1] = cy * sz; mat[0][2] = -sy; tmp = sx * sy; mat[1][0] = tmp * cz - cx * sz; mat[1][1] = tmp * sz + cx * cz; mat[1][2] = sx * cy; tmp = cx * sy; mat[2][0] = tmp * cz + sx * sz; mat[2][1] = tmp * sz - sx * cz; mat[2][2] = cx * cy; mat[0][3] = mat[1][3] = mat[2][3] = 0.0f; mat[3][0] = mat[3][1] = mat[3][2] = 0.0f; mat[3][3] = 1.0f; #ifdef SLOW tv.x = mat[0][0]; tv.y = mat[0][1]; tv.z = mat[0][2]; tv = Vnorm( tv ); mat[0][0] = tv.x; mat[0][1] = tv.y; mat[0][2] = tv.z; tv.x = mat[1][0]; tv.y = mat[1][1]; tv.z = mat[1][2]; tv = Vnorm( tv ); mat[1][0] = tv.x; mat[1][1] = tv.y; mat[1][2] = tv.z; tv.x = mat[2][0]; tv.y = mat[2][1]; tv.z = mat[2][2]; tv = Vnorm( tv ); mat[2][0] = tv.x; mat[2][1] = tv.y; mat[2][2] = tv.z; #endif } /* matinit initialize the matrix 'mat' to identity */ static void matinit( Matrix mat /* matrix to initialize (in-out) */ ) { register int i, j; for( i = 0; i < 4; i++ ) for( j = 0; j < 4; j++ ) mat[i][j] = ( i == j ? 1.0 : 0.0 ); } static float VDot( VERT a, VERT b ) { return ( ( float ) ( a.x * b.x ) + ( a.y * b.y ) + ( a.z * b.z ) ); } /* return the cross product c = a cross b */ /* return vector sum c = a+b */ static VERT Vadd( VERT a, VERT b ) { VERT c; c.x = a.x + b.x; c.y = a.y + b.y; c.z = a.z + b.z; return ( c ); } /* return vector difference c = a-b */ static VERT Vsub( VERT a, VERT b ) { VERT c; c.x = a.x - b.x; c.y = a.y - b.y; c.z = a.z - b.z; return ( c ); } /* multiply two vectors together component-wise and return the result */ static VERT VMult( VERT a, VERT b ) { VERT result; result.x = a.x * b.x; result.y = a.y * b.y; result.z = a.z * b.z; return ( result ); } /* return the dot product of vectors a and b */ static float VSquaredLength( VERT a ) { return ( ( float ) ( ( a.x * a.x ) + ( a.y * a.y ) + ( a.z * a.z ) ) ); } /* returns length of input vector */ static float Vlength( VERT a ) { return ( ( float ) sqrt( ( double ) VSquaredLength( a ) ) ); } /* scales the input vector to the new length and returns it */ static VERT Vscale( VERT v, float newlen ) { VERT ret; float len = Vlength( v ); ret.x = 0; ret.y = 0; ret.z = 0; if( len != 0.0 ) { ret.x *= newlen / len; ret.y *= newlen / len; ret.z *= newlen / len; } return ( ret ); } /* normalizes the input vector and returns it */ // qqqqq /* Matrix Inversion by Richard Carling from "Graphics Gems", Academic Press, 1990 */ #define SMALL_NUMBER 1.e-8 /* * inverse( original_matrix, inverse_matrix ) * * calculate the inverse of a 4x4 matrix * * -1 * A = ___1__ adjoint A * det A */ //#include "joemath.h" #include static double det4x4( Matrix m ); static double det2x2( double a, double b, double c, double d ); static double det3x3( double, double, double, double, double, double, double, double, double ); static void adjoint( Matrix in, Matrix out ); static void inverse( Matrix in, Matrix out ) { int i, j; double det; /* calculate the adjoint matrix */ adjoint( in, out ); /* calculate the 4x4 determinant * if the determinant is zero, * then the inverse matrix is not unique. */ det = det4x4( in ); if( fabs( det ) < SMALL_NUMBER ) { // printf("Non-singular matrix, no inverse!\n"); /*exit(1); */ matinit( out ); } /* scale the adjoint matrix to get the inverse */ for( i = 0; i < 4; i++ ) for( j = 0; j < 4; j++ ) out[i][j] = out[i][j] / det; } /* * adjoint( original_matrix, inverse_matrix ) * * calculate the adjoint of a 4x4 matrix * * Let a denote the minor determinant of matrix A obtained by * ij * * deleting the ith row and jth column from A. * * i+j * Let b = (-1) a * ij ji * * The matrix B = (b ) is the adjoint of A * ij */ void adjoint( Matrix in, Matrix out ) //Matrix in; Matrix out; { double a1, a2, a3, a4, b1, b2, b3, b4; double c1, c2, c3, c4, d1, d2, d3, d4; /* assign to individual variable names to aid */ /* selecting correct values */ a1 = in[0][0]; b1 = in[0][1]; c1 = in[0][2]; d1 = in[0][3]; a2 = in[1][0]; b2 = in[1][1]; c2 = in[1][2]; d2 = in[1][3]; a3 = in[2][0]; b3 = in[2][1]; c3 = in[2][2]; d3 = in[2][3]; a4 = in[3][0]; b4 = in[3][1]; c4 = in[3][2]; d4 = in[3][3]; /* row column labeling reversed since we transpose rows & columns */ out[0][0] = det3x3( b2, b3, b4, c2, c3, c4, d2, d3, d4 ); out[1][0] = -det3x3( a2, a3, a4, c2, c3, c4, d2, d3, d4 ); out[2][0] = det3x3( a2, a3, a4, b2, b3, b4, d2, d3, d4 ); out[3][0] = -det3x3( a2, a3, a4, b2, b3, b4, c2, c3, c4 ); out[0][1] = -det3x3( b1, b3, b4, c1, c3, c4, d1, d3, d4 ); out[1][1] = det3x3( a1, a3, a4, c1, c3, c4, d1, d3, d4 ); out[2][1] = -det3x3( a1, a3, a4, b1, b3, b4, d1, d3, d4 ); out[3][1] = det3x3( a1, a3, a4, b1, b3, b4, c1, c3, c4 ); out[0][2] = det3x3( b1, b2, b4, c1, c2, c4, d1, d2, d4 ); out[1][2] = -det3x3( a1, a2, a4, c1, c2, c4, d1, d2, d4 ); out[2][2] = det3x3( a1, a2, a4, b1, b2, b4, d1, d2, d4 ); out[3][2] = -det3x3( a1, a2, a4, b1, b2, b4, c1, c2, c4 ); out[0][3] = -det3x3( b1, b2, b3, c1, c2, c3, d1, d2, d3 ); out[1][3] = det3x3( a1, a2, a3, c1, c2, c3, d1, d2, d3 ); out[2][3] = -det3x3( a1, a2, a3, b1, b2, b3, d1, d2, d3 ); out[3][3] = det3x3( a1, a2, a3, b1, b2, b3, c1, c2, c3 ); } /* * double = det4x4( matrix ) * * calculate the determinant of a 4x4 matrix. */ double det4x4( Matrix m ) { double ans; double a1, a2, a3, a4, b1, b2, b3, b4, c1, c2, c3, c4, d1, d2, d3, d4; /* assign to individual variable names to aid selecting */ /* correct elements */ a1 = m[0][0]; b1 = m[0][1]; c1 = m[0][2]; d1 = m[0][3]; a2 = m[1][0]; b2 = m[1][1]; c2 = m[1][2]; d2 = m[1][3]; a3 = m[2][0]; b3 = m[2][1]; c3 = m[2][2]; d3 = m[2][3]; a4 = m[3][0]; b4 = m[3][1]; c4 = m[3][2]; d4 = m[3][3]; ans = a1 * det3x3( b2, b3, b4, c2, c3, c4, d2, d3, d4 ) - b1 * det3x3( a2, a3, a4, c2, c3, c4, d2, d3, d4 ) + c1 * det3x3( a2, a3, a4, b2, b3, b4, d2, d3, d4 ) - d1 * det3x3( a2, a3, a4, b2, b3, b4, c2, c3, c4 ); return ans; } /* * double = det3x3( a1, a2, a3, b1, b2, b3, c1, c2, c3 ) * * calculate the determinant of a 3x3 matrix * in the form * * | a1, b1, c1 | * | a2, b2, c2 | * | a3, b3, c3 | */ double det3x3( double a1, double a2, double a3, double b1, double b2, double b3, double c1, double c2, double c3 ) { double ans; // double det2x2(); ans = a1 * det2x2( b2, b3, c2, c3 ) - b1 * det2x2( a2, a3, c2, c3 ) + c1 * det2x2( a2, a3, b2, b3 ); return ans; } /* * double = det2x2( double a, double b, double c, double d ) * * calculate the determinant of a 2x2 matrix. */ double det2x2( double a, double b, double c, double d ) { double ans; ans = a * d - b * c; return ans; } static int Fround( float f ) { int ret; if( f > 0 ) ret = ( int ) ( f + 0.5 ); // Posotive number else ret = ( int ) ( f - 0.5 ); // Negative number return ( ret ); } static float getnoise( float u, float v, float z ) //float u,v,z; { int ui, vi, zi, ui2, vi2, zi2; float b1, a, b, c, a1, b2, a2, ur, vr, zr, iur, ivr, izr; //#ifdef crap //if (GmultON==0) { //u+=(float)Gcurpass*21.5; //v+=(float)Gcurpass*21.1; //z+=(float)Gcurpass*22.323; } //#endif if( u < 0 ) u = -u; if( v < 0 ) v = -v; if( z < 0 ) z = -z; // while (u>18.0f) u-=18.0f; // while (v>18.0f) v-=18.0f; // while (z>18.0f) z-=18.0f; ui = ( int ) floor( u ); vi = ( int ) floor( v ); zi = ( int ) floor( z ); zr = z - zi; ur = u - ui; vr = v - vi; iur = 1.0f-ur; ivr = 1.0f-vr; izr = 1.0f-zr; ui = ui % 17; vi = vi % 17; zi = zi % 17; ui2 = ui + 1; vi2 = vi + 1; zi2 = zi + 1; ui2 = ui2 % 17; vi2 = vi2 % 17; zi2 = zi2 % 17; a = noise[ui][vi][zi]; b = noise[ui2][vi][zi]; a1 = a * iur + b * ur; a = noise[ui][vi2][zi]; b = noise[ui2][vi2][zi]; a2 = a * iur + b * ur; b1 = a1 * ivr + a2 * vr; a = noise[ui][vi][zi2]; b = noise[ui2][vi][zi2]; a1 = a * iur + b * ur; a = noise[ui][vi2][zi2]; b = noise[ui2][vi2][zi2]; a2 = a * iur + b * ur; b2 = a1 * ivr + a2 * vr; c = b1 * izr + b2 * zr; // if( Gfast_eval == 0 ) // c = smoothstep( c ); // c=(float)rand()/(float)RAND_MAX; if( c > 1 ) c = 1; if( c < 0 ) c = 0; return ( c ); } static VERT coord_cross( VERT a, VERT b ) { VERT out; out.x = ( a.y * b.z ) - ( a.z * b.y ); out.y = ( a.z * b.x ) - ( a.x * b.z ); out.z = ( a.x * b.y ) - ( a.y * b.x ); return ( out ); } static VERT M_normalize( VERT m ) //VERT m; { float d; d = fsqrt( m.x * m.x + m.y * m.y + m.z * m.z ); if( d > 0.0f ) { m.x /= d; m.y /= d; m.z /= d; } return ( m ); } static VERT Vnorm( VERT v ) //VERT v; { VERT ret; float len; len = ( float ) fsqrt( ( double ) ( v.x * v.x + v.y * v.y + v.z * v.z ) ); ret = v; if( len != 0.0f ) { ret.x /= len; ret.y /= len; ret.z /= len; } if( len == 0.0 ) { ret.x = .0000001f; ret.y = 0; ret.x = 0; } return ( ret ); } /* fsqrt.c * * A fast square root program adapted from the code of * Paul Lalonde and Robert Dawson in Graphics Gems I. * The format of IEEE double precision floating point numbers is: * * SEEEEEEEEEEEMMMM MMMMMMMMMMMMMMMM MMMMMMMMMMMMMMMM MMMMMMMMMMMMMMMM * * S = Sign bit for whole number * E = Exponent bit (exponent in excess 1023 form) * M = Mantissa bit */ /* MOST_SIG_OFFSET gives the (int *) offset from the address of the double * to the part of the number containing the sign and exponent. * You will need to find the relevant offset for your architecture. */ #define MOST_SIG_OFFSET 1 /* SQRT_TAB_SIZE - the size of the lookup table - must be a power of four. */ #define SQRT_TAB_SIZE 16384 /* MANT_SHIFTS is the number of shifts to move mantissa into position. * If you quadruple the table size subtract two from this constant, * if you quarter the table size then add two. * Valid values are: (16384, 7) (4096, 9) (1024, 11) (256, 13) */ #define MANT_SHIFTS 7 #define EXP_BIAS 1023 /* Exponents are always positive */ #define EXP_SHIFTS 20 /* Shifs exponent to least sig. bits */ #define EXP_LSB 0x00100000 /* 1 << EXP_SHIFTS */ #define MANT_MASK 0x000FFFFF /* Mask to extract mantissa */ static int sqrt_tab[SQRT_TAB_SIZE]; static void init_sqrt_tab( void ) { int i; double f; unsigned int *fi = ( unsigned int * ) &f + MOST_SIG_OFFSET; for( i = 0; i < SQRT_TAB_SIZE / 2; i++ ) { f = 0; /* Clears least sig part */ *fi = ( i << MANT_SHIFTS ) | ( EXP_BIAS << EXP_SHIFTS ); f = sqrt( f ); sqrt_tab[i] = *fi & MANT_MASK; f = 0; /* Clears least sig part */ *fi = ( i << MANT_SHIFTS ) | ( ( EXP_BIAS + 1 ) << EXP_SHIFTS ); f = sqrt( f ); sqrt_tab[i + SQRT_TAB_SIZE / 2] = *fi & MANT_MASK; } } #ifndef IRIX static float fsqrt( float ff ) { ff = sqrt( ff ); return ( ff ); } #endif static float fsqrt2( float ff ) { double f; unsigned int e; unsigned int *fi = ( unsigned int * ) &f + MOST_SIG_OFFSET; f = ( double ) ff; if( f == 0.0f ) return ( 0.0f ); e = ( *fi >> EXP_SHIFTS ) - EXP_BIAS; *fi &= MANT_MASK; if( e & 1 ) *fi |= EXP_LSB; e >>= 1; *fi = ( sqrt_tab[*fi >> MANT_SHIFTS] ) | ( ( e + EXP_BIAS ) << EXP_SHIFTS ); ff = ( float ) f; return ( ff ); } static void dump_sqrt_tab( void ) { int i, nl = 0; printf( "unsigned int sqrt_tab[] = {\n" ); for( i = 0; i < SQRT_TAB_SIZE - 1; i++ ) { printf( "0x%x,", sqrt_tab[i] ); nl++; if( nl > 8 ) { nl = 0; putchar( '\n' ); } } printf( "0x%x\n", sqrt_tab[SQRT_TAB_SIZE - 1] ); printf( "};\n" ); }