util.c 3.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
  1. /*Copyright (C) 2008-2009 Timothy B. Terriberry (tterribe@xiph.org)
  2. You can redistribute this library and/or modify it under the terms of the
  3. GNU Lesser General Public License as published by the Free Software
  4. Foundation; either version 2.1 of the License, or (at your option) any later
  5. version.*/
  6. //#include <stdlib.h>
  7. #include "util.h"
  8. /*Computes floor(sqrt(_val)) exactly.*/
  9. unsigned qr_isqrt(unsigned _val){
  10. unsigned g;
  11. unsigned b;
  12. int bshift;
  13. /*Uses the second method from
  14. http://www.azillionmonkeys.com/qed/sqroot.html
  15. The main idea is to search for the largest binary digit b such that
  16. (g+b)*(g+b) <= _val, and add it to the solution g.*/
  17. g=0;
  18. b=0x8000;
  19. for(bshift=16;bshift-->0;){
  20. unsigned t;
  21. t=(g<<1)+b<<bshift;
  22. if(t<=_val){
  23. g+=b;
  24. _val-=t;
  25. }
  26. b>>=1;
  27. }
  28. return g;
  29. }
  30. /*Computes sqrt(_x*_x+_y*_y) using CORDIC.
  31. This implementation is valid for all 32-bit inputs and returns a result
  32. accurate to about 27 bits of precision.
  33. It has been tested for all postiive 16-bit inputs, where it returns correctly
  34. rounded results in 99.998% of cases and the maximum error is
  35. 0.500137134862598032 (for _x=48140, _y=63018).
  36. Very nearly all results less than (1<<16) are correctly rounded.
  37. All Pythagorean triples with a hypotenuse of less than ((1<<27)-1) evaluate
  38. correctly, and the total bias over all Pythagorean triples is -0.04579, with
  39. a relative RMS error of 7.2864E-10 and a relative peak error of 7.4387E-9.*/
  40. unsigned qr_ihypot(int _x,int _y){
  41. unsigned x;
  42. unsigned y;
  43. int mask;
  44. int shift;
  45. int u;
  46. int v;
  47. int i;
  48. x=_x=abs(_x);
  49. y=_y=abs(_y);
  50. mask=-(x>y)&(_x^_y);
  51. x^=mask;
  52. y^=mask;
  53. _y^=mask;
  54. shift=31-qr_ilog(y);
  55. shift=QR_MAXI(shift,0);
  56. x=(unsigned)((x<<shift)*0x9B74EDAAULL>>32);
  57. _y=(int)((_y<<shift)*0x9B74EDA9LL>>32);
  58. u=x;
  59. mask=-(_y<0);
  60. x+=_y+mask^mask;
  61. _y-=u+mask^mask;
  62. u=x+1>>1;
  63. v=_y+1>>1;
  64. mask=-(_y<0);
  65. x+=v+mask^mask;
  66. _y-=u+mask^mask;
  67. for(i=1;i<16;i++){
  68. int r;
  69. u=x+1>>2;
  70. r=(1<<2*i)>>1;
  71. v=_y+r>>2*i;
  72. mask=-(_y<0);
  73. x+=v+mask^mask;
  74. _y=_y-(u+mask^mask)<<1;
  75. }
  76. return x+((1U<<shift)>>1)>>shift;
  77. }
  78. #if defined(__GNUC__) && defined(HAVE_FEATURES_H)
  79. # include <features.h>
  80. # if __GNUC_PREREQ(3,4)
  81. # include <limits.h>
  82. # if INT_MAX>=2147483647
  83. # define QR_CLZ0 sizeof(unsigned)*CHAR_BIT
  84. # define QR_CLZ(_x) (__builtin_clz(_x))
  85. # elif LONG_MAX>=2147483647L
  86. # define QR_CLZ0 sizeof(unsigned long)*CHAR_BIT
  87. # define QR_CLZ(_x) (__builtin_clzl(_x))
  88. # endif
  89. # endif
  90. #endif
  91. int qr_ilog(unsigned _v){
  92. #if defined(QR_CLZ)
  93. /*Note that __builtin_clz is not defined when _x==0, according to the gcc
  94. documentation (and that of the x86 BSR instruction that implements it), so
  95. we have to special-case it.*/
  96. return QR_CLZ0-QR_CLZ(_v)&-!!_v;
  97. #else
  98. int ret;
  99. int m;
  100. m=!!(_v&0xFFFF0000)<<4;
  101. _v>>=m;
  102. ret=m;
  103. m=!!(_v&0xFF00)<<3;
  104. _v>>=m;
  105. ret|=m;
  106. m=!!(_v&0xF0)<<2;
  107. _v>>=m;
  108. ret|=m;
  109. m=!!(_v&0xC)<<1;
  110. _v>>=m;
  111. ret|=m;
  112. ret|=!!(_v&0x2);
  113. return ret + !!_v;
  114. #endif
  115. }
  116. #if defined(QR_TEST_SQRT)
  117. #include <math.h>
  118. #include <stdio.h>
  119. /*Exhaustively test the integer square root function.*/
  120. int main(void){
  121. unsigned u;
  122. u=0;
  123. do{
  124. unsigned r;
  125. unsigned s;
  126. r=qr_isqrt(u);
  127. s=(int)sqrt(u);
  128. if(r!=s)printf("%u: %u!=%u\n",u,r,s);
  129. u++;
  130. }
  131. while(u);
  132. return 0;
  133. }
  134. #endif