ref: 2036c57b8ab86067c598cee6c9d4726fb383a0f5
parent: c7ace558fe3faa30701011fce37aba522963ffa9
author: Timothy B. Terriberry <[email protected]>
date: Sun Dec 14 09:40:34 EST 2008
Ensure that log2_frac() is _really_ an upper bound. This version has actually been tested for all 32-bit inputs.
--- a/libcelt/cwrs.c
+++ b/libcelt/cwrs.c
@@ -49,33 +49,35 @@
#include "mathops.h"
#include "arch.h"
+/*Guaranteed to return a conservatively large estimate of the binary logarithm
+ with frac bits of fractional precision.
+ Tested for all possible 32-bit inputs with frac=4, where the maximum
+ overestimation is 0.06254243 bits.*/
int log2_frac(ec_uint32 val, int frac)
{- int i;
- /* EC_ILOG() actually returns log2()+1, go figure */
- int L = EC_ILOG(val)-1;
- int pow2;
- pow2=!(val&val-1);
- /*printf ("in: %d %d ", val, L);*/- if (L>14)
- val >>= L-14;
- else if (L<14)
- val <<= 14-L;
- L <<= frac;
- /*printf ("%d\n", val);*/- for (i=0;i<frac;i++)
-{- val = (val*val) >> 15;
- /*printf ("%d\n", val);*/- if (val > 16384)
- L |= (1<<(frac-i-1));
- else
- val <<= 1;
-}
- /*The previous loop returns a conservatively low estimate.
- If val wasn't a power of two, we should round up.*/
- L += !pow2;
- return L;
+ int l;
+ l=EC_ILOG(val);
+ if(val&val-1){+ /*This is (val>>l-16), but guaranteed to round up, even if adding a bias
+ before the shift would cause overflow (e.g., for 0xFFFFxxxx).*/
+ if(l>16)val=(val>>l-16)+((val&(1<<l-16)-1)+(1<<l-16)-1>>l-16);
+ else val<<=16-l;
+ l=l-1<<frac;
+ /*Note that we always need one iteration, since the rounding up above means
+ that we might need to adjust the integer part of the logarithm.*/
+ do{+ int b;
+ b=(int)(val>>16);
+ l+=b<<frac;
+ val=val+b>>b;
+ val=val*val+0x7FFF>>15;
+ }
+ while(frac-->0);
+ /*If val is not exactly 0x8000, then we have to round up the remainder.*/
+ return l+(val>0x8000);
+ }
+ /*Exact powers of two require no rounding.*/
+ else return l-1<<frac;
}
int fits_in32(int _n, int _m)