home *** CD-ROM | disk | FTP | other *** search
/ SuperHack / SuperHack CD.bin / CODING / CPP / FASTINT.ZIP / FASTINT.CPP
Encoding:
C/C++ Source or Header  |  1997-01-09  |  7.1 KB  |  330 lines
  1. // integer trig functions by Karl Lager
  2. // 76752,2243@compuserve.com
  3.  
  4. // The following code is freeware.  You may prune and modify it for
  5. // your own use, but if you share the code please include the origional
  6. // fastint.cpp file.
  7.  
  8. // 9/12/96
  9. // a slight error in the fastasin and fastacos has been corrected.
  10. // When averaging two numbers, they should shift the sum right, not
  11. // left.  They should work now.
  12. // The high precision trig functions are present, but have not all
  13. // been thouroughly tested.  Caveat hacker.
  14.  
  15. #include <math.h>
  16.  
  17. /**************************************************************/
  18. //                 high precision trig
  19.  
  20. long hsinof[1025]; // sin of 0-90 deg
  21. int htanhashtable[1025];
  22. //
  23.  
  24.  
  25. int isinof[360],icosof[360]; //sin&cos <<14
  26. int tanhashtable[129];
  27.  
  28.  
  29. void init_trig()
  30. { double t;int i;
  31.   for (i = 0;i<=128;i++)
  32.    { t= i;
  33.     tanhashtable[i] = atan(t/128)*180/M_PI;
  34.    }
  35.   for (i = 0;i<=1025;i++)
  36.    { t= i;
  37.     htanhashtable[i] = atan(t/1024)*32768/M_PI;
  38.    }
  39.   for (i=0;i<360;i++)
  40.    { t=M_PI/180*i;
  41.      isinof[i] = sin(t)*(1<<14);
  42.      icosof[i] = cos(t)*(1<<14);
  43.    }
  44.   for (i=0;i<1025;i++)
  45.    { t =M_PI/2048L*i;
  46.      hsinof[i] = sin(t)*65536L;
  47.    }
  48. }
  49.  
  50.  
  51. int fastatan2(long y, long x)
  52. {int t,d,a;
  53.  long tan;
  54.  if (x == 0 && y == 0) return 0;//bulletproof
  55.  if (x >= 0)
  56.  { if (y >= 0)
  57.    { if (x>y) { tan = (y<<7)/x; d = 1; a = 0;}
  58.      else { tan = (x<<7)/y; d = 0; a = 90;}
  59.    }
  60.    else
  61.    { if (x>-y) { tan = -(y<<7)/x; d = 0; a = 360;}
  62.      else { tan = (x<<7)/-y; d = 1; a = 270;}
  63.    }
  64.  }
  65.  else
  66.  { if (y >= 0)
  67.    { if (-x>y) { tan = (y<<7)/-x; d = 0; a = 180; }
  68.      else { tan = (-x<<7)/y; d = 1; a = 90;}
  69.    }
  70.    else
  71.    { if (-x>-y) { tan = (y<<7)/x; d = 1; a = 180; }
  72.      else { tan = (x<<7)/y; d = 0; a = 270;}
  73.    }
  74.  }
  75.  t = tanhashtable[(int)tan];
  76.  if(d) t = a+t;
  77.  else t = a-t;
  78.  if (t==360) t=0;
  79.  return t;
  80.  }
  81.  
  82.  
  83. int fastasin(int x)
  84. { int a,b,i;
  85.  if (x<0)
  86.   {a=270;b=359;}
  87.  else
  88.   {a=0;b=90;}
  89.  while(a<b-1)      //binary search
  90.  { i= (a+b)>>1;
  91.    if (isinof[i]<x)
  92.      a=i;
  93.    else
  94.      b=i;
  95.  }
  96.  if (x-isinof[a] > isinof[b]-x) return b;
  97.  else return a;
  98. }
  99.  
  100. int fastacos(int x) // 90-fastasin(x)
  101. { int a,b,i;
  102.   a=0;b=180;
  103.  while(a<b-1)      //binary search
  104.  { i= (a+b)>>1;
  105.    if (icosof[i]>x)
  106.      a=i;
  107.    else
  108.      b=i;
  109.  }
  110.  if (icosof[a]-x > x-icosof[b]) return b; //interpolation would go here for
  111.  else return a;                       //higher precision angles
  112. }
  113.  
  114.  
  115. unsigned sqrti(unsigned long x)//{float f=x;unsigned i=sqrt(f);return i;}
  116. { unsigned long rt;
  117.   long i;
  118.   if (x==0) return 0;
  119.   if (x<4) return 1;
  120. //  for(i=2;(x>>i)>0;i+=2);  // find top bit
  121.   asm { mov cx,0
  122.         mov eax, x
  123.       }
  124. sqrtlooptop:
  125.   asm { inc cx
  126.         shr eax,2
  127.         jnz sqrtlooptop
  128.       }
  129.  
  130.   //rt = x>>(i>>1);
  131.   asm { mov eax,x;
  132.         shr eax,cl
  133.         mov rt, eax
  134.       }
  135.   if (rt == 0xFFFF) return (unsigned)rt;
  136.   while (x>=rt*rt+(rt<<1)+1)
  137.  
  138.    { i = (x/rt-rt)>>1;          //(a²+b²)=a²+2ab+b²; c²-a²=2ab+b²
  139.      rt += i;                                //   (b ~= c²-a²)/2a
  140.      while(x<=rt*rt-(rt<<1)+1)
  141.       { i = (rt-x/rt)>>1;
  142.         rt -= i;
  143.       }
  144.    }
  145.   if (rt*rt<x) {if (x-rt*rt>rt*rt+(rt<<1)+1-x) rt++;}
  146.   else { if (rt*rt-x>x-(rt*rt-(rt<<1)+1)) rt--;}
  147.   return (unsigned)rt;
  148. }
  149.  
  150. // high precision trig functions
  151. // 64k = 1 rev  (65536 points around a circle)
  152.  
  153. long hsin(unsigned x) // 16.16 fixed point
  154. { int s;
  155.   int index;
  156.   long sin1;
  157.  
  158.   if (x > 32768)     // if x > 180 sign is negative
  159.     { s = 1;
  160.       x -=32768;
  161.     }
  162.     else s = 0;
  163.  
  164.   if (x > 16384)
  165.     x = 32768 - x;   // x is now between 0 and 90
  166.   // return hsinof[x>>4] for no interpolation
  167.   index = x>>4;
  168.   sin1 = hsinof[index];
  169.   if ((x&15)>0)        // that's x%16
  170.   sin1 += (hsinof[index+1]-sin1)*(x&15) >> 4;
  171.   if (s == 0) return sin1;
  172.   else return - sin1;
  173.  
  174. }
  175.  
  176. long hcos(unsigned x) // 16.16 fixed point
  177. { int s;
  178.   int index;
  179.   long sin1;
  180.  
  181.   x = 16384 - x;     // cos(x) = sin(90-x)
  182.  
  183.   if (x > 32768)
  184.     { s = 1;
  185.       x -=32768;
  186.     }
  187.     else s = 0;
  188.  
  189.   if (x > 16384)
  190.     x = 32768 - x;   // x is now between 0 and 90° (0-16k)
  191.   index = x>>4;
  192.   sin1 = hsinof[index];
  193.   if ((x&15)>0)
  194.   sin1 += (hsinof[index+1]-sin1)*(x&15) >> 4;
  195.   if (s == 0) return sin1;
  196.   else return -sin1;
  197.  
  198. }
  199.  
  200. unsigned hfastatan2(long y, long x)
  201. {int d;
  202.  unsigned t,a;
  203.  unsigned f;
  204.  long tan;
  205.  if (x == 0 && y == 0) return 0;
  206.  if (x >= 0)
  207.  { if (y >= 0)
  208.    { if (x>y)
  209.      { tan = (y<<10)/x;
  210.        f = (((y<<10)-tan*x)<<4)/x;
  211.        d = 1; a = 0;
  212.      }
  213.      else
  214.      { tan = (x<<10)/y;
  215.        f = (((x<<10)-tan*y)<<4)/y;
  216.        d = 0; a = 16384;}
  217.    }
  218.    else
  219.    { if (x>-y)
  220.      { tan = -(y<<10)/x;
  221.        f = (((-y<<10)-tan*x)<<4)/x;
  222.        d = 0; a = 0; //== 65536;
  223.      }
  224.      else
  225.      { tan = (x<<10)/-y;
  226.        f = (((x<<10)+tan*y)<<4)/-y;
  227.        d = 1; a = 49152;
  228.      }
  229.    }
  230.  }
  231.  else
  232.  { if (y >= 0)
  233.    { if (-x>y)
  234.      { tan = (y<<10)/-x;
  235.        f = (((y<<10)+tan*x)<<4)/-x;
  236.        d = 0; a = 32768;
  237.      }
  238.      else
  239.      { tan = (-x<<10)/y;
  240.        f = (((-x<<10)-tan*y)<<4)/y;
  241.        d = 1; a = 16384;
  242.      }
  243.    }
  244.    else
  245.    { if (-x>-y)
  246.      { tan = (y<<10)/x;
  247.        f = (((y<<10)-tan*x)<<4)/x;
  248.        d = 1; a = 32768;
  249.      }
  250.      else
  251.      { tan = (x<<10)/y;
  252.        f = (((x<<10)-tan*y)<<4)/y;
  253.        d = 0; a = 49152;}
  254.    }
  255.  }
  256.  t = tanhashtable[(int)tan];
  257.                                // add fraction for interpolation
  258.  if (f>0) t += (tanhashtable[(int)tan+1]-t)*f >> 4;
  259.  if(d) t = a+t;
  260.  else t = a-t;
  261. // if (t==65536) t=0; // handled by wraparound
  262.  return t;
  263.  }
  264.  
  265.  
  266. unsigned hfastasin(long x)   // x is in the range -64k to 64k
  267. { int a,b,i,s,result;
  268.  if (x<0)
  269.   {s = 1;x = -x;} //{a=270;b=359;}
  270.  else
  271.   s = 0;
  272.  
  273.   a=0;b=1024;
  274.  while(a<b-1)      //binary search
  275.  { i= (a+b)>>1;
  276.    if (hsinof[i]<x)
  277.      a=i;
  278.    else
  279.      b=i;
  280.  }
  281.  // if (x-isinof[a] > isinof[b]-x) return b;
  282.                    // interpolate
  283.  result = (a<<4) + ((x-hsinof[a])<<4) / (hsinof[b]-hsinof[a]);
  284.  
  285.  if (s==0) return result;
  286.  else return -result;
  287. }
  288.  
  289.  
  290. unsigned hfastacos(long x)   // x is in the range -64k to 64k
  291. { int a,b,i,s,result;
  292.  if (x<0)
  293.   {s = 1;x = -x;} //{a=270;b=359;}
  294.  else
  295.   s = 0;
  296.  
  297.   a=0;b=1024;
  298.  while(a<b-1)      //binary search
  299.  { i= (a+b)>>1;
  300.    if (hsinof[i]<x)
  301.      a=i;
  302.    else
  303.      b=i;
  304.  }
  305.  // if (x-isinof[a] > isinof[b]-x) return b;
  306.                    // interpolate
  307.  result = (a<<4) + ((x-hsinof[a])<<4) / (hsinof[b]-hsinof[a]);
  308.  
  309.  if (s==0) result = 16384 - result;
  310.  else result = 16384 + result;
  311.  return result;
  312. }
  313.  
  314.  
  315. // And finally, a quick way to get integers out of the coprocessor..
  316.  
  317. const double d2i = 6755399441055744.0;
  318. const float f2i = 12582912.0;
  319. double tempdub;
  320. float tempfl;
  321. #define fl2int(x,i)asm{fld dword ptr x; fadd f2i; fstp tempfl; fwait; mov eax,dword ptr tempfl;sub eax,f2i;mov i,eax}
  322. #define dub2int(x,i) asm{fld qword ptr x; fadd d2i; fstp tempdub; fwait; mov eax,dword ptr tempdub;mov i,eax}
  323.  
  324. long float2int(double x)
  325.  { long i;
  326.    dub2int(x,i);
  327.    return i;
  328.  }
  329.  
  330.