ref: d7e016b552696b8f5581f5093f82a598bc0033d6
dir: /mathlib.c/
#include <u.h> #include <libc.h> #include "dat.h" #include "quakedef.h" #include "fns.h" vec3_t vec3_origin = {0,0,0}; int nanmask = 255<<23; /*-----------------------------------------------------------------*/ #define DEG2RAD( a ) ( a * M_PI ) / 180.0F void ProjectPointOnPlane( vec3_t dst, const vec3_t p, const vec3_t normal ) { float d; vec3_t n; float inv_denom; inv_denom = 1.0F / DotProduct( normal, normal ); d = DotProduct( normal, p ) * inv_denom; n[0] = normal[0] * inv_denom; n[1] = normal[1] * inv_denom; n[2] = normal[2] * inv_denom; dst[0] = p[0] - d * n[0]; dst[1] = p[1] - d * n[1]; dst[2] = p[2] - d * n[2]; } /* ** assumes "src" is normalized */ void PerpendicularVector( vec3_t dst, const vec3_t src ) { int pos; int i; float minelem = 1.0F; vec3_t tempvec; /* ** find the smallest magnitude axially aligned vector */ for ( pos = 0, i = 0; i < 3; i++ ) { if ( fabs( src[i] ) < minelem ) { pos = i; minelem = fabs( src[i] ); } } tempvec[0] = tempvec[1] = tempvec[2] = 0.0F; tempvec[pos] = 1.0F; /* ** project the point onto the plane defined by src */ ProjectPointOnPlane( dst, tempvec, src ); /* ** normalize the result */ VectorNormalize( dst ); } void RotatePointAroundVector( vec3_t dst, const vec3_t dir, const vec3_t point, float degrees ) { float m[3][3]; float im[3][3]; float zrot[3][3]; float tmpmat[3][3]; float rot[3][3]; int i; vec3_t vr, vup, vf; vf[0] = dir[0]; vf[1] = dir[1]; vf[2] = dir[2]; PerpendicularVector( vr, dir ); CrossProduct( vr, vf, vup ); m[0][0] = vr[0]; m[1][0] = vr[1]; m[2][0] = vr[2]; m[0][1] = vup[0]; m[1][1] = vup[1]; m[2][1] = vup[2]; m[0][2] = vf[0]; m[1][2] = vf[1]; m[2][2] = vf[2]; memcpy(im, m, sizeof im); im[0][1] = m[1][0]; im[0][2] = m[2][0]; im[1][0] = m[0][1]; im[1][2] = m[2][1]; im[2][0] = m[0][2]; im[2][1] = m[1][2]; memset(zrot, 0, sizeof zrot); zrot[0][0] = zrot[1][1] = zrot[2][2] = 1.0F; zrot[0][0] = cos( DEG2RAD( degrees ) ); zrot[0][1] = sin( DEG2RAD( degrees ) ); zrot[1][0] = -sin( DEG2RAD( degrees ) ); zrot[1][1] = cos( DEG2RAD( degrees ) ); R_ConcatRotations( m, zrot, tmpmat ); R_ConcatRotations( tmpmat, im, rot ); for ( i = 0; i < 3; i++ ) { dst[i] = rot[i][0] * point[0] + rot[i][1] * point[1] + rot[i][2] * point[2]; } } /*-----------------------------------------------------------------*/ float anglemod(float a) { /* if (a >= 0) a -= 360*(int)(a/360); else a += 360*( 1 + (int)(-a/360) ); */ a = (360.0/65536) * ((int)(a*(65536/360.0)) & 65535); return a; } /* ================== BOPS_Error Split out like this for ASM to call. ================== */ void BOPS_Error (void) { fatal ("BoxOnPlaneSide: Bad signbits"); } /* ================== BoxOnPlaneSide Returns 1, 2, or 1 + 2 ================== */ int BoxOnPlaneSide (vec3_t emins, vec3_t emaxs, mplane_t *p) { float dist1, dist2; int sides; /* // this is done by the BOX_ON_PLANE_SIDE macro before calling this function // fast axial cases if (p->type < 3) { if (p->dist <= emins[p->type]) return 1; if (p->dist >= emaxs[p->type]) return 2; return 3; } */ // general case switch (p->signbits) { case 0: dist1 = p->normal[0]*emaxs[0] + p->normal[1]*emaxs[1] + p->normal[2]*emaxs[2]; dist2 = p->normal[0]*emins[0] + p->normal[1]*emins[1] + p->normal[2]*emins[2]; break; case 1: dist1 = p->normal[0]*emins[0] + p->normal[1]*emaxs[1] + p->normal[2]*emaxs[2]; dist2 = p->normal[0]*emaxs[0] + p->normal[1]*emins[1] + p->normal[2]*emins[2]; break; case 2: dist1 = p->normal[0]*emaxs[0] + p->normal[1]*emins[1] + p->normal[2]*emaxs[2]; dist2 = p->normal[0]*emins[0] + p->normal[1]*emaxs[1] + p->normal[2]*emins[2]; break; case 3: dist1 = p->normal[0]*emins[0] + p->normal[1]*emins[1] + p->normal[2]*emaxs[2]; dist2 = p->normal[0]*emaxs[0] + p->normal[1]*emaxs[1] + p->normal[2]*emins[2]; break; case 4: dist1 = p->normal[0]*emaxs[0] + p->normal[1]*emaxs[1] + p->normal[2]*emins[2]; dist2 = p->normal[0]*emins[0] + p->normal[1]*emins[1] + p->normal[2]*emaxs[2]; break; case 5: dist1 = p->normal[0]*emins[0] + p->normal[1]*emaxs[1] + p->normal[2]*emins[2]; dist2 = p->normal[0]*emaxs[0] + p->normal[1]*emins[1] + p->normal[2]*emaxs[2]; break; case 6: dist1 = p->normal[0]*emaxs[0] + p->normal[1]*emins[1] + p->normal[2]*emins[2]; dist2 = p->normal[0]*emins[0] + p->normal[1]*emaxs[1] + p->normal[2]*emaxs[2]; break; case 7: dist1 = p->normal[0]*emins[0] + p->normal[1]*emins[1] + p->normal[2]*emins[2]; dist2 = p->normal[0]*emaxs[0] + p->normal[1]*emaxs[1] + p->normal[2]*emaxs[2]; break; default: dist1 = dist2 = 0; // shut up compiler BOPS_Error (); break; } /* int i; vec3_t corners[2]; for (i=0 ; i<3 ; i++) { if (plane->normal[i] < 0) { corners[0][i] = emins[i]; corners[1][i] = emaxs[i]; } else { corners[1][i] = emins[i]; corners[0][i] = emaxs[i]; } } dist = DotProduct (plane->normal, corners[0]) - plane->dist; dist2 = DotProduct (plane->normal, corners[1]) - plane->dist; sides = 0; if (dist1 >= 0) sides = 1; if (dist2 < 0) sides |= 2; */ sides = 0; if (dist1 >= p->dist) sides = 1; if (dist2 < p->dist) sides |= 2; #ifdef PARANOID if (sides == 0) fatal ("BoxOnPlaneSide: sides==0"); #endif return sides; } void AngleVectors (vec3_t angles, vec3_t forward, vec3_t right, vec3_t up) { float angle; float sr, sp, sy, cr, cp, cy; angle = angles[YAW] * (M_PI*2 / 360); sy = sin(angle); cy = cos(angle); angle = angles[PITCH] * (M_PI*2 / 360); sp = sin(angle); cp = cos(angle); angle = angles[ROLL] * (M_PI*2 / 360); sr = sin(angle); cr = cos(angle); forward[0] = cp*cy; forward[1] = cp*sy; forward[2] = -sp; right[0] = (-1*sr*sp*cy+-1*cr*-sy); right[1] = (-1*sr*sp*sy+-1*cr*cy); right[2] = -1*sr*cp; up[0] = (cr*sp*cy+-sr*-sy); up[1] = (cr*sp*sy+-sr*cy); up[2] = cr*cp; } int VectorCompare (vec3_t v1, vec3_t v2) { int i; for (i=0 ; i<3 ; i++) if (v1[i] != v2[i]) return 0; return 1; } void VectorMA (vec3_t veca, float scale, vec3_t vecb, vec3_t vecc) { vecc[0] = veca[0] + scale*vecb[0]; vecc[1] = veca[1] + scale*vecb[1]; vecc[2] = veca[2] + scale*vecb[2]; } vec_t _DotProduct (vec3_t v1, vec3_t v2) { return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]; } void _VectorSubtract (vec3_t veca, vec3_t vecb, vec3_t out) { out[0] = veca[0]-vecb[0]; out[1] = veca[1]-vecb[1]; out[2] = veca[2]-vecb[2]; } void _VectorAdd (vec3_t veca, vec3_t vecb, vec3_t out) { out[0] = veca[0]+vecb[0]; out[1] = veca[1]+vecb[1]; out[2] = veca[2]+vecb[2]; } void _VectorCopy (vec3_t in, vec3_t out) { out[0] = in[0]; out[1] = in[1]; out[2] = in[2]; } void CrossProduct (vec3_t v1, vec3_t v2, vec3_t cross) { cross[0] = v1[1]*v2[2] - v1[2]*v2[1]; cross[1] = v1[2]*v2[0] - v1[0]*v2[2]; cross[2] = v1[0]*v2[1] - v1[1]*v2[0]; } double sqrt(double x); vec_t Length(vec3_t v) { int i; float length; length = 0; for (i=0 ; i< 3 ; i++) length += v[i]*v[i]; length = sqrt (length); // FIXME return length; } float VectorNormalize (vec3_t v) { float length, ilength; length = v[0]*v[0] + v[1]*v[1] + v[2]*v[2]; length = sqrt (length); // FIXME if (length) { ilength = 1/length; v[0] *= ilength; v[1] *= ilength; v[2] *= ilength; } return length; } void VectorInverse (vec3_t v) { v[0] = -v[0]; v[1] = -v[1]; v[2] = -v[2]; } void VectorScale (vec3_t in, vec_t scale, vec3_t out) { out[0] = in[0]*scale; out[1] = in[1]*scale; out[2] = in[2]*scale; } /* ================ R_ConcatRotations ================ */ void R_ConcatRotations (float in1[3][3], float in2[3][3], float out[3][3]) { out[0][0] = in1[0][0] * in2[0][0] + in1[0][1] * in2[1][0] + in1[0][2] * in2[2][0]; out[0][1] = in1[0][0] * in2[0][1] + in1[0][1] * in2[1][1] + in1[0][2] * in2[2][1]; out[0][2] = in1[0][0] * in2[0][2] + in1[0][1] * in2[1][2] + in1[0][2] * in2[2][2]; out[1][0] = in1[1][0] * in2[0][0] + in1[1][1] * in2[1][0] + in1[1][2] * in2[2][0]; out[1][1] = in1[1][0] * in2[0][1] + in1[1][1] * in2[1][1] + in1[1][2] * in2[2][1]; out[1][2] = in1[1][0] * in2[0][2] + in1[1][1] * in2[1][2] + in1[1][2] * in2[2][2]; out[2][0] = in1[2][0] * in2[0][0] + in1[2][1] * in2[1][0] + in1[2][2] * in2[2][0]; out[2][1] = in1[2][0] * in2[0][1] + in1[2][1] * in2[1][1] + in1[2][2] * in2[2][1]; out[2][2] = in1[2][0] * in2[0][2] + in1[2][1] * in2[1][2] + in1[2][2] * in2[2][2]; } /* ================ R_ConcatTransforms ================ */ void R_ConcatTransforms (float in1[3][4], float in2[3][4], float out[3][4]) { out[0][0] = in1[0][0] * in2[0][0] + in1[0][1] * in2[1][0] + in1[0][2] * in2[2][0]; out[0][1] = in1[0][0] * in2[0][1] + in1[0][1] * in2[1][1] + in1[0][2] * in2[2][1]; out[0][2] = in1[0][0] * in2[0][2] + in1[0][1] * in2[1][2] + in1[0][2] * in2[2][2]; out[0][3] = in1[0][0] * in2[0][3] + in1[0][1] * in2[1][3] + in1[0][2] * in2[2][3] + in1[0][3]; out[1][0] = in1[1][0] * in2[0][0] + in1[1][1] * in2[1][0] + in1[1][2] * in2[2][0]; out[1][1] = in1[1][0] * in2[0][1] + in1[1][1] * in2[1][1] + in1[1][2] * in2[2][1]; out[1][2] = in1[1][0] * in2[0][2] + in1[1][1] * in2[1][2] + in1[1][2] * in2[2][2]; out[1][3] = in1[1][0] * in2[0][3] + in1[1][1] * in2[1][3] + in1[1][2] * in2[2][3] + in1[1][3]; out[2][0] = in1[2][0] * in2[0][0] + in1[2][1] * in2[1][0] + in1[2][2] * in2[2][0]; out[2][1] = in1[2][0] * in2[0][1] + in1[2][1] * in2[1][1] + in1[2][2] * in2[2][1]; out[2][2] = in1[2][0] * in2[0][2] + in1[2][1] * in2[1][2] + in1[2][2] * in2[2][2]; out[2][3] = in1[2][0] * in2[0][3] + in1[2][1] * in2[1][3] + in1[2][2] * in2[2][3] + in1[2][3]; } /* =================== FloorDivMod Returns mathematically correct (floor-based) quotient and remainder for numer and denom, both of which should contain no fractional part. The quotient must fit in 32 bits. ==================== */ void FloorDivMod (double numer, double denom, int *quotient, int *rem) { int q, r; double x; #ifdef PARANOID if (denom <= 0.0) fatal ("FloorDivMod: bad denominator %d\n", denom); if (floor(numer) != numer || floor(denom) != denom) fatal ("FloorDivMod: non-integer numer or denom %f %f\n", numer, denom); #endif if (numer >= 0.0) { x = floor(numer / denom); q = (int)x; r = (int)floor(numer - (x * denom)); } else { // // perform operations with positive values, and fix mod to make floor-based // x = floor(-numer / denom); q = -(int)x; r = (int)floor(-numer - (x * denom)); if (r != 0) { q--; r = (int)denom - r; } } *quotient = q; *rem = r; } /* =================== GreatestCommonDivisor ==================== */ int GreatestCommonDivisor (int i1, int i2) { if (i1 > i2) { if (i2 == 0) return (i1); return GreatestCommonDivisor (i2, i1 % i2); } else { if (i1 == 0) return (i2); return GreatestCommonDivisor (i1, i2 % i1); } } // TODO: move to nonintel.c /* =================== Invert24To16 Inverts an 8.24 value to a 16.16 value ==================== */ fixed16_t Invert24To16(fixed16_t val) { if (val < 256) return (0xFFFFFFFF); return (fixed16_t) (((double)0x10000 * (double)0x1000000 / (double)val) + 0.5); }