The Tor Browser / annotate

media/libopus/celt/mathops.c@ac0c01689b40 (annotated)

media/libopus/celt/mathops.c

Thu, 15 Jan 2015 15:59:08 +0100

author
Michael Schloh von Bennewitz <michael@schloh.com>
date
Thu, 15 Jan 2015 15:59:08 +0100
branch
TOR_BUG_9701
changeset 10
ac0c01689b40
permissions
-rw-r--r--

Implement a real Private Browsing Mode condition by changing the API/ABI;
This solves Tor bug #9701, complying with disk avoidance documented in
https://www.torproject.org/projects/torbrowser/design/#disk-avoidance.

michael@0 1 /* Copyright (c) 2002-2008 Jean-Marc Valin
michael@0 2 Copyright (c) 2007-2008 CSIRO
michael@0 3 Copyright (c) 2007-2009 Xiph.Org Foundation
michael@0 4 Written by Jean-Marc Valin */
michael@0 5 /**
michael@0 6 @file mathops.h
michael@0 7 @brief Various math functions
michael@0 8 */
michael@0 9 /*
michael@0 10 Redistribution and use in source and binary forms, with or without
michael@0 11 modification, are permitted provided that the following conditions
michael@0 12 are met:
michael@0 13
michael@0 14 - Redistributions of source code must retain the above copyright
michael@0 15 notice, this list of conditions and the following disclaimer.
michael@0 16
michael@0 17 - Redistributions in binary form must reproduce the above copyright
michael@0 18 notice, this list of conditions and the following disclaimer in the
michael@0 19 documentation and/or other materials provided with the distribution.
michael@0 20
michael@0 21 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
michael@0 22 ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
michael@0 23 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
michael@0 24 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
michael@0 25 OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
michael@0 26 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
michael@0 27 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
michael@0 28 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
michael@0 29 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
michael@0 30 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
michael@0 31 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
michael@0 32 */
michael@0 33
michael@0 34 #ifdef HAVE_CONFIG_H
michael@0 35 #include "config.h"
michael@0 36 #endif
michael@0 37
michael@0 38 #include "mathops.h"
michael@0 39
michael@0 40 /*Compute floor(sqrt(_val)) with exact arithmetic.
michael@0 41 This has been tested on all possible 32-bit inputs.*/
michael@0 42 unsigned isqrt32(opus_uint32 _val){
michael@0 43 unsigned b;
michael@0 44 unsigned g;
michael@0 45 int bshift;
michael@0 46 /*Uses the second method from
michael@0 47 http://www.azillionmonkeys.com/qed/sqroot.html
michael@0 48 The main idea is to search for the largest binary digit b such that
michael@0 49 (g+b)*(g+b) <= _val, and add it to the solution g.*/
michael@0 50 g=0;
michael@0 51 bshift=(EC_ILOG(_val)-1)>>1;
michael@0 52 b=1U<<bshift;
michael@0 53 do{
michael@0 54 opus_uint32 t;
michael@0 55 t=(((opus_uint32)g<<1)+b)<<bshift;
michael@0 56 if(t<=_val){
michael@0 57 g+=b;
michael@0 58 _val-=t;
michael@0 59 }
michael@0 60 b>>=1;
michael@0 61 bshift--;
michael@0 62 }
michael@0 63 while(bshift>=0);
michael@0 64 return g;
michael@0 65 }
michael@0 66
michael@0 67 #ifdef FIXED_POINT
michael@0 68
michael@0 69 opus_val32 frac_div32(opus_val32 a, opus_val32 b)
michael@0 70 {
michael@0 71 opus_val16 rcp;
michael@0 72 opus_val32 result, rem;
michael@0 73 int shift = celt_ilog2(b)-29;
michael@0 74 a = VSHR32(a,shift);
michael@0 75 b = VSHR32(b,shift);
michael@0 76 /* 16-bit reciprocal */
michael@0 77 rcp = ROUND16(celt_rcp(ROUND16(b,16)),3);
michael@0 78 result = MULT16_32_Q15(rcp, a);
michael@0 79 rem = PSHR32(a,2)-MULT32_32_Q31(result, b);
michael@0 80 result = ADD32(result, SHL32(MULT16_32_Q15(rcp, rem),2));
michael@0 81 if (result >= 536870912) /* 2^29 */
michael@0 82 return 2147483647; /* 2^31 - 1 */
michael@0 83 else if (result <= -536870912) /* -2^29 */
michael@0 84 return -2147483647; /* -2^31 */
michael@0 85 else
michael@0 86 return SHL32(result, 2);
michael@0 87 }
michael@0 88
michael@0 89 /** Reciprocal sqrt approximation in the range [0.25,1) (Q16 in, Q14 out) */
michael@0 90 opus_val16 celt_rsqrt_norm(opus_val32 x)
michael@0 91 {
michael@0 92 opus_val16 n;
michael@0 93 opus_val16 r;
michael@0 94 opus_val16 r2;
michael@0 95 opus_val16 y;
michael@0 96 /* Range of n is [-16384,32767] ([-0.5,1) in Q15). */
michael@0 97 n = x-32768;
michael@0 98 /* Get a rough initial guess for the root.
michael@0 99 The optimal minimax quadratic approximation (using relative error) is
michael@0 100 r = 1.437799046117536+n*(-0.823394375837328+n*0.4096419668459485).
michael@0 101 Coefficients here, and the final result r, are Q14.*/
michael@0 102 r = ADD16(23557, MULT16_16_Q15(n, ADD16(-13490, MULT16_16_Q15(n, 6713))));
michael@0 103 /* We want y = x*r*r-1 in Q15, but x is 32-bit Q16 and r is Q14.
michael@0 104 We can compute the result from n and r using Q15 multiplies with some
michael@0 105 adjustment, carefully done to avoid overflow.
michael@0 106 Range of y is [-1564,1594]. */
michael@0 107 r2 = MULT16_16_Q15(r, r);
michael@0 108 y = SHL16(SUB16(ADD16(MULT16_16_Q15(r2, n), r2), 16384), 1);
michael@0 109 /* Apply a 2nd-order Householder iteration: r += r*y*(y*0.375-0.5).
michael@0 110 This yields the Q14 reciprocal square root of the Q16 x, with a maximum
michael@0 111 relative error of 1.04956E-4, a (relative) RMSE of 2.80979E-5, and a
michael@0 112 peak absolute error of 2.26591/16384. */
michael@0 113 return ADD16(r, MULT16_16_Q15(r, MULT16_16_Q15(y,
michael@0 114 SUB16(MULT16_16_Q15(y, 12288), 16384))));
michael@0 115 }
michael@0 116
michael@0 117 /** Sqrt approximation (QX input, QX/2 output) */
michael@0 118 opus_val32 celt_sqrt(opus_val32 x)
michael@0 119 {
michael@0 120 int k;
michael@0 121 opus_val16 n;
michael@0 122 opus_val32 rt;
michael@0 123 static const opus_val16 C[5] = {23175, 11561, -3011, 1699, -664};
michael@0 124 if (x==0)
michael@0 125 return 0;
michael@0 126 else if (x>=1073741824)
michael@0 127 return 32767;
michael@0 128 k = (celt_ilog2(x)>>1)-7;
michael@0 129 x = VSHR32(x, 2*k);
michael@0 130 n = x-32768;
michael@0 131 rt = ADD16(C[0], MULT16_16_Q15(n, ADD16(C[1], MULT16_16_Q15(n, ADD16(C[2],
michael@0 132 MULT16_16_Q15(n, ADD16(C[3], MULT16_16_Q15(n, (C[4])))))))));
michael@0 133 rt = VSHR32(rt,7-k);
michael@0 134 return rt;
michael@0 135 }
michael@0 136
michael@0 137 #define L1 32767
michael@0 138 #define L2 -7651
michael@0 139 #define L3 8277
michael@0 140 #define L4 -626
michael@0 141
michael@0 142 static OPUS_INLINE opus_val16 _celt_cos_pi_2(opus_val16 x)
michael@0 143 {
michael@0 144 opus_val16 x2;
michael@0 145
michael@0 146 x2 = MULT16_16_P15(x,x);
michael@0 147 return ADD16(1,MIN16(32766,ADD32(SUB16(L1,x2), MULT16_16_P15(x2, ADD32(L2, MULT16_16_P15(x2, ADD32(L3, MULT16_16_P15(L4, x2
michael@0 148 ))))))));
michael@0 149 }
michael@0 150
michael@0 151 #undef L1
michael@0 152 #undef L2
michael@0 153 #undef L3
michael@0 154 #undef L4
michael@0 155
michael@0 156 opus_val16 celt_cos_norm(opus_val32 x)
michael@0 157 {
michael@0 158 x = x&0x0001ffff;
michael@0 159 if (x>SHL32(EXTEND32(1), 16))
michael@0 160 x = SUB32(SHL32(EXTEND32(1), 17),x);
michael@0 161 if (x&0x00007fff)
michael@0 162 {
michael@0 163 if (x<SHL32(EXTEND32(1), 15))
michael@0 164 {
michael@0 165 return _celt_cos_pi_2(EXTRACT16(x));
michael@0 166 } else {
michael@0 167 return NEG32(_celt_cos_pi_2(EXTRACT16(65536-x)));
michael@0 168 }
michael@0 169 } else {
michael@0 170 if (x&0x0000ffff)
michael@0 171 return 0;
michael@0 172 else if (x&0x0001ffff)
michael@0 173 return -32767;
michael@0 174 else
michael@0 175 return 32767;
michael@0 176 }
michael@0 177 }
michael@0 178
michael@0 179 /** Reciprocal approximation (Q15 input, Q16 output) */
michael@0 180 opus_val32 celt_rcp(opus_val32 x)
michael@0 181 {
michael@0 182 int i;
michael@0 183 opus_val16 n;
michael@0 184 opus_val16 r;
michael@0 185 celt_assert2(x>0, "celt_rcp() only defined for positive values");
michael@0 186 i = celt_ilog2(x);
michael@0 187 /* n is Q15 with range [0,1). */
michael@0 188 n = VSHR32(x,i-15)-32768;
michael@0 189 /* Start with a linear approximation:
michael@0 190 r = 1.8823529411764706-0.9411764705882353*n.
michael@0 191 The coefficients and the result are Q14 in the range [15420,30840].*/
michael@0 192 r = ADD16(30840, MULT16_16_Q15(-15420, n));
michael@0 193 /* Perform two Newton iterations:
michael@0 194 r -= r*((r*n)-1.Q15)
michael@0 195 = r*((r*n)+(r-1.Q15)). */
michael@0 196 r = SUB16(r, MULT16_16_Q15(r,
michael@0 197 ADD16(MULT16_16_Q15(r, n), ADD16(r, -32768))));
michael@0 198 /* We subtract an extra 1 in the second iteration to avoid overflow; it also
michael@0 199 neatly compensates for truncation error in the rest of the process. */
michael@0 200 r = SUB16(r, ADD16(1, MULT16_16_Q15(r,
michael@0 201 ADD16(MULT16_16_Q15(r, n), ADD16(r, -32768)))));
michael@0 202 /* r is now the Q15 solution to 2/(n+1), with a maximum relative error
michael@0 203 of 7.05346E-5, a (relative) RMSE of 2.14418E-5, and a peak absolute
michael@0 204 error of 1.24665/32768. */
michael@0 205 return VSHR32(EXTEND32(r),i-16);
michael@0 206 }
michael@0 207
michael@0 208 #endif

Mercurial Repository: The Tor Browser

Europalab SCM Repository Server

mercurial