2 THE COMPUTER CODE CONTAINED HEREIN IS THE SOLE PROPERTY OF PARALLAX
3 SOFTWARE CORPORATION ("PARALLAX"). PARALLAX, IN DISTRIBUTING THE CODE TO
4 END-USERS, AND SUBJECT TO ALL OF THE TERMS AND CONDITIONS HEREIN, GRANTS A
5 ROYALTY-FREE, PERPETUAL LICENSE TO SUCH END-USERS FOR USE BY SUCH END-USERS
6 IN USING, DISPLAYING, AND CREATING DERIVATIVE WORKS THEREOF, SO LONG AS
7 SUCH USE, DISPLAY OR CREATION IS FOR NON-COMMERCIAL, ROYALTY OR REVENUE
8 FREE PURPOSES. IN NO EVENT SHALL THE END-USER USE THE COMPUTER CODE
9 CONTAINED HEREIN FOR REVENUE-BEARING PURPOSES. THE END-USER UNDERSTANDS
10 AND AGREES TO THE TERMS HEREIN AND ACCEPTS THE SAME BY USE OF THIS FILE.
11 COPYRIGHT 1993-1998 PARALLAX SOFTWARE CORPORATION. ALL RIGHTS RESERVED.
14 * $Source: /cvs/cvsroot/d2x/maths/fixc.c,v $
17 * $Date: 2001-10-31 07:41:54 $
19 * C version of fixed point library
21 * $Log: not supported by cvs2svn $
22 * Revision 1.2 2001/01/31 15:18:04 bradleyb
23 * Makefile and conf.h fixes
25 * Revision 1.1.1.1 2001/01/19 03:29:58 bradleyb
28 * Revision 1.3 1999/10/18 00:31:01 donut
29 * more alpha fixes from Falk Hueffner
31 * Revision 1.2 1999/08/05 22:53:41 sekmu
33 * D3D patch(es) from ADB
35 * Revision 1.1.1.1 1999/06/14 22:13:35 donut
36 * Import of d1x 1.37 source.
38 * Revision 1.7 1995/09/22 14:08:16 allender
39 * fixed fix_atan2 to work correctly with doubles
41 * Revision 1.6 1995/08/31 15:43:49 allender
42 * *** empty log message ***
44 * Revision 1.5 1995/07/05 16:15:15 allender
45 * make fixmuldiv use doubles for PPC implementation
47 * Revision 1.4 1995/05/15 13:57:36 allender
48 * make fixmuldiv compile when compiling under 68k
50 * Revision 1.3 1995/05/11 13:02:59 allender
51 * some routines are now in assembly
53 * Revision 1.2 1995/05/04 20:04:45 allender
54 * use MPW fixdiv if compiling with MPW (why did I do this?)
56 * Revision 1.1 1995/04/17 11:37:54 allender
61 * Revision 1.1 1995/03/08 18:55:09 matt
72 static char rcsid[] = "$Id: fixc.c,v 1.3 2001-10-31 07:41:54 bradleyb Exp $";
83 #pragma message ("warning: FIX NOT INLINED")
85 #warning "FIX NOT INLINED"
89 extern ubyte guess_table[];
90 extern short sincos_table[];
91 extern ushort asin_table[];
92 extern ushort acos_table[];
93 extern fix isqrt_guess_table[];
96 void fixquadnegate(quad *q)
99 q->high = 0 - q->high - (q->low != 0);
102 //multiply two ints & add 64-bit result to 64-bit sum
103 void fixmulaccum(quad *q,fix a,fix b)
106 u_int32_t ah,al,bh,bl;
112 aa = labs(a); bb = labs(b);
114 ah = aa>>16; al = aa&0xffff;
115 bh = bb>>16; bl = bb&0xffff;
124 if (q->low < old) q->high++;
128 if (q->low < old) q->high++;
130 q->high += ah*bh + (t>>16) + c;
137 //extract a fix from a quad product
138 fix fixquadadjust(quad *q)
140 return (q->high<<16) + (q->low>>16);
144 #define EPSILON (F1_0/100)
147 #define QLONG __int64
149 #define QLONG long long
153 fix fixmul(fix a, fix b) {
154 /* return (fix)(((double)a*(double)b)/65536.0);*/
156 asm("imul %%edx; shrd $16,%%edx,%%eax" : "=a" (ret) : "a" (a), "d" (b) : "%edx");
158 return (fix)((((QLONG)a)*b) >> 16);
161 fix fixdiv(fix a, fix b)
163 /* return (fix)(((double)a * 65536.0) / (double)b);*/
164 return (fix)((((QLONG)a) << 16)/b);
166 asm("mov %%eax,%%edx; sar $16,%%edx; shl $16,%%eax; idiv %%ebx" : "=a" (ret) : "a" (a), "b" (b) : "%edx");
170 fix fixmuldiv(fix a, fix b, fix c)
173 asm("imul %%edx; idiv %%ebx" : "=a" (ret) : "a" (a), "d" (b), "b" (c) : "%edx");
178 d = (double)a * (double) b;
179 return (fix)(d / (double) c);
181 return (fix)((((QLONG)a)*b)/c);
185 //given cos & sin of an angle, return that angle.
186 //parms need not be normalized, that is, the ratio of the parms cos/sin must
187 //equal the ratio of the actual cos & sin for the result angle, but the parms
188 //need not be the actual cos & sin.
189 //NOTE: this is different from the standard C atan2, since it is left-handed.
191 fixang fix_atan2(fix cos,fix sin)
193 double d, dsin, dcos;
196 //Assert(!(cos==0 && sin==0));
198 //find smaller of two
202 d = sqrt((dsin * dsin) + (dcos * dcos));
207 if (labs(sin) < labs(cos)) { //sin is smaller, use arcsin
208 t = fix_asin((fix)((dsin / d) * 65536.0));
214 t = fix_acos((fix)((dcos / d) * 65536.0));
222 //divide a quad by a fix, returning a fix
223 int32_t fixdivquadlong(u_int32_t nl,u_int32_t nh,u_int32_t d)
233 Q = ((nh&0x80000000)!=0);
234 M = ((d&0x80000000)!=0);
239 for (i=0; i<32; i++ ) {
243 T = ((nl&0x80000000L)!=0);
249 Q = (unsigned char)((0x80000000L & nh) != 0 );
250 nh = (nh << 1) | (u_int32_t)T;
258 Q = (unsigned char)(tmp1 == 0);
261 Q = (unsigned char)((0x80000000L & nh) != 0 );
262 nh = (nh << 1) | (u_int32_t)T;
270 Q = (unsigned char)(tmp1 == 0);
278 for (i=0; i<32; i++ ) {
282 T = ((nl&0x80000000L)!=0);
288 Q = (unsigned char)((0x80000000L & nh) != 0 );
289 nh = (nh << 1) | (u_int32_t)T;
297 Q = (unsigned char)(tmp1 == 0);
300 Q = (unsigned char)((0x80000000L & nh) != 0 );
301 nh = (nh << 1) | (u_int32_t)T;
309 Q = (unsigned char)(tmp1 == 0);
321 unsigned int fixdivquadlongu(uint nl, uint nh, uint d)
323 return fixdivquadlong((u_int32_t) nl,(u_int32_t) nh,(u_int32_t) d);
326 int32_t fixdivquadlong(u_int32_t nl,u_int32_t nh,u_int32_t d) {
330 :"a" (nl), "d" (nh), "r" (d)
335 static inline u_int32_t fixdivquadlongu(u_int32_t nl,u_int32_t nh,u_int32_t d) {
339 :"a" (nl), "d" (nh), "r" (d)
346 u_int32_t quad_sqrt(u_int32_t low,int32_t high)
355 if (high==0 && (int32_t)low>=0)
356 return long_sqrt((int32_t)low);
358 if (high & 0xff000000) {
359 cnt=12+16; i = high >> 24;
360 } else if (high & 0xff0000) {
361 cnt=8+16; i = high >> 16;
362 } else if (high & 0xff00) {
363 cnt=4+16; i = high >> 8;
368 r = guess_table[i]<<cnt;
370 //quad loop usually executed 4 times
372 r = (fixdivquadlongu(low,high,r)+r)/2;
373 r = (fixdivquadlongu(low,high,r)+r)/2;
374 r = (fixdivquadlongu(low,high,r)+r)/2;
379 t = fixdivquadlongu(low,high,r);
386 } while (!(r==t || r==old_r));
388 t = fixdivquadlongu(low,high,r);
389 //edited 05/17/99 Matt Mueller - tq.high is undefined here.. so set them to = 0
392 fixmulaccum(&tq,r,t);
393 if (tq.low!=low || tq.high!=high)
399 //computes the square root of a long, returning a short
400 ushort long_sqrt(int32_t a)
409 else if (a & 0xff0000)
416 r = guess_table[(a>>cnt)&0xff]<<cnt;
418 //the loop nearly always executes 3 times, so we'll unroll it 2 times and
419 //not do any checking until after the third time. By my calcutations, the
420 //loop is executed 2 times in 99.97% of cases, 3 times in 93.65% of cases,
421 //four times in 16.18% of cases, and five times in 0.44% of cases. It never
422 //executes more than five times. By timing, I determined that is is faster
423 //to always execute three times and not check for termination the first two
424 //times through. This means that in 93.65% of cases, we save 6 cmp/jcc pairs,
425 //and in 6.35% of cases we do an extra divide. In real life, these numbers
426 //might not be the same.
441 } while (!(r==t || r==old_r));
449 //computes the square root of a fix, returning a fix
452 return ((fix) long_sqrt(a)) << 8;
456 //compute sine and cosine of an angle, filling in the variables
457 //either of the pointers can be NULL
459 void fix_sincos(fix a,fix *s,fix *c)
467 ss = sincos_table[i];
468 if (s) *s = (ss + (((sincos_table[i+1] - ss) * f)>>8))<<2;
470 cc = sincos_table[i+64];
471 if (c) *c = (cc + (((sincos_table[i+64+1] - cc) * f)>>8))<<2;
474 //compute sine and cosine of an angle, filling in the variables
475 //either of the pointers can be NULL
477 void fix_fastsincos(fix a,fix *s,fix *c)
483 if (s) *s = sincos_table[i] << 2;
484 if (c) *c = sincos_table[i+64] << 2;
487 //compute inverse sine
488 fixang fix_asin(fix v)
495 if (vv >= f1_0) //check for out of range
502 aa = aa + (((asin_table[i+1] - aa) * f)>>8);
510 //compute inverse cosine
511 fixang fix_acos(fix v)
518 if (vv >= f1_0) //check for out of range
525 aa = aa + (((acos_table[i+1] - aa) * f)>>8);
533 #define TABLE_SIZE 1024
535 //for passed value a, returns 1/sqrt(a)
536 fix fix_isqrt( fix a )
542 if ( a == 0 ) return 0;
544 while( b >= TABLE_SIZE ) {
549 //printf( "Count = %d (%d>>%d)\n", cnt, b, (cnt+1)/2 );
550 r = isqrt_guess_table[b] >> ((cnt+1)/2);
552 //printf( "Initial r = %d\n", r );
554 for (i=0; i<3; i++ ) {
556 r = fixmul( ( (3*65536) - fixmul(fixmul(r,r),a) ), r) / 2;
557 //printf( "r %d = %d\n", i, r );
558 if ( old_r >= r ) return (r+old_r)/2;