Reorders some operations in anti-collapse to reuse values
[opus.git] / libcelt / mathops.c
index 2a9ab5e..494649e 100644 (file)
 
 #include "mathops.h"
 
+/*Compute floor(sqrt(_val)) with exact arithmetic.
+  This has been tested on all possible 32-bit inputs.*/
+unsigned isqrt32(celt_uint32 _val){
+  unsigned b;
+  unsigned g;
+  int      bshift;
+  /*Uses the second method from
+     http://www.azillionmonkeys.com/qed/sqroot.html
+    The main idea is to search for the largest binary digit b such that
+     (g+b)*(g+b) <= _val, and add it to the solution g.*/
+  g=0;
+  bshift=EC_ILOG(_val)-1>>1;
+  b=1U<<bshift;
+  do{
+    celt_uint32 t;
+    t=((celt_uint32)g<<1)+b<<bshift;
+    if(t<=_val){
+      g+=b;
+      _val-=t;
+    }
+    b>>=1;
+    bshift--;
+  }
+  while(bshift>=0);
+  return g;
+}
+
 #ifdef FIXED_POINT
 
 celt_word32 frac_div32(celt_word32 a, celt_word32 b)
 {
    celt_word16 rcp;
    celt_word32 result, rem;
-   int shift = 30-celt_ilog2(b);
-   a = SHL32(a,shift);
-   b = SHL32(b,shift);
-
+   int shift = celt_ilog2(b)-29;
+   a = VSHR32(a,shift);
+   b = VSHR32(b,shift);
    /* 16-bit reciprocal */
-   rcp = ROUND16(celt_rcp(ROUND16(b,16)),2);
-   result = SHL32(MULT16_32_Q15(rcp, a),1);
+   rcp = ROUND16(celt_rcp(ROUND16(b,16)),3);
+   result = SHL32(MULT16_32_Q15(rcp, a),2);
    rem = a-MULT32_32_Q31(result, b);
-   result += SHL32(MULT16_32_Q15(rcp, rem),1);
+   result += SHL32(MULT16_32_Q15(rcp, rem),2);
    return result;
 }