From 2036c57b8ab86067c598cee6c9d4726fb383a0f5 Mon Sep 17 00:00:00 2001 From: "Timothy B. Terriberry" Date: Sun, 14 Dec 2008 14:40:34 -0500 Subject: [PATCH] Ensure that log2_frac() is _really_ an upper bound. This version has actually been tested for all 32-bit inputs. --- libcelt/cwrs.c | 52 +++++++++++++++++++++++++++------------------------- 1 file changed, 27 insertions(+), 25 deletions(-) diff --git a/libcelt/cwrs.c b/libcelt/cwrs.c index a39c1a1d..fa5af87d 100644 --- 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> 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); + else val<<=16-l; + l=l-1<>16); + l+=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<