1*a58d3d2aSXin Li /* Copyright (c) 2008-2011 Xiph.Org Foundation, Mozilla Corporation,
2*a58d3d2aSXin Li Gregory Maxwell
3*a58d3d2aSXin Li Written by Jean-Marc Valin, Gregory Maxwell, and Timothy B. Terriberry */
4*a58d3d2aSXin Li /*
5*a58d3d2aSXin Li Redistribution and use in source and binary forms, with or without
6*a58d3d2aSXin Li modification, are permitted provided that the following conditions
7*a58d3d2aSXin Li are met:
8*a58d3d2aSXin Li
9*a58d3d2aSXin Li - Redistributions of source code must retain the above copyright
10*a58d3d2aSXin Li notice, this list of conditions and the following disclaimer.
11*a58d3d2aSXin Li
12*a58d3d2aSXin Li - Redistributions in binary form must reproduce the above copyright
13*a58d3d2aSXin Li notice, this list of conditions and the following disclaimer in the
14*a58d3d2aSXin Li documentation and/or other materials provided with the distribution.
15*a58d3d2aSXin Li
16*a58d3d2aSXin Li THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
17*a58d3d2aSXin Li ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
18*a58d3d2aSXin Li LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
19*a58d3d2aSXin Li A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
20*a58d3d2aSXin Li OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
21*a58d3d2aSXin Li EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
22*a58d3d2aSXin Li PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
23*a58d3d2aSXin Li PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
24*a58d3d2aSXin Li LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
25*a58d3d2aSXin Li NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
26*a58d3d2aSXin Li SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27*a58d3d2aSXin Li */
28*a58d3d2aSXin Li
29*a58d3d2aSXin Li #ifdef HAVE_CONFIG_H
30*a58d3d2aSXin Li #include "config.h"
31*a58d3d2aSXin Li #endif
32*a58d3d2aSXin Li
33*a58d3d2aSXin Li #ifndef CUSTOM_MODES
34*a58d3d2aSXin Li #define CUSTOM_MODES
35*a58d3d2aSXin Li #endif
36*a58d3d2aSXin Li
37*a58d3d2aSXin Li #include <stdio.h>
38*a58d3d2aSXin Li #include <math.h>
39*a58d3d2aSXin Li #include "mathops.h"
40*a58d3d2aSXin Li #include "bands.h"
41*a58d3d2aSXin Li
42*a58d3d2aSXin Li #ifdef FIXED_POINT
43*a58d3d2aSXin Li #define WORD "%d"
44*a58d3d2aSXin Li #else
45*a58d3d2aSXin Li #define WORD "%f"
46*a58d3d2aSXin Li #endif
47*a58d3d2aSXin Li
48*a58d3d2aSXin Li int ret = 0;
49*a58d3d2aSXin Li
testdiv(void)50*a58d3d2aSXin Li void testdiv(void)
51*a58d3d2aSXin Li {
52*a58d3d2aSXin Li opus_int32 i;
53*a58d3d2aSXin Li for (i=1;i<=327670;i++)
54*a58d3d2aSXin Li {
55*a58d3d2aSXin Li double prod;
56*a58d3d2aSXin Li opus_val32 val;
57*a58d3d2aSXin Li val = celt_rcp(i);
58*a58d3d2aSXin Li #ifdef FIXED_POINT
59*a58d3d2aSXin Li prod = (1./32768./65526.)*val*i;
60*a58d3d2aSXin Li #else
61*a58d3d2aSXin Li prod = val*i;
62*a58d3d2aSXin Li #endif
63*a58d3d2aSXin Li if (fabs(prod-1) > .00025)
64*a58d3d2aSXin Li {
65*a58d3d2aSXin Li fprintf (stderr, "div failed: 1/%d="WORD" (product = %f)\n", i, val, prod);
66*a58d3d2aSXin Li ret = 1;
67*a58d3d2aSXin Li }
68*a58d3d2aSXin Li }
69*a58d3d2aSXin Li }
70*a58d3d2aSXin Li
testsqrt(void)71*a58d3d2aSXin Li void testsqrt(void)
72*a58d3d2aSXin Li {
73*a58d3d2aSXin Li opus_int32 i;
74*a58d3d2aSXin Li for (i=1;i<=1000000000;i++)
75*a58d3d2aSXin Li {
76*a58d3d2aSXin Li double ratio;
77*a58d3d2aSXin Li opus_val16 val;
78*a58d3d2aSXin Li val = celt_sqrt(i);
79*a58d3d2aSXin Li ratio = val/sqrt(i);
80*a58d3d2aSXin Li if (fabs(ratio - 1) > .0005 && fabs(val-sqrt(i)) > 2)
81*a58d3d2aSXin Li {
82*a58d3d2aSXin Li fprintf (stderr, "sqrt failed: sqrt(%d)="WORD" (ratio = %f)\n", i, val, ratio);
83*a58d3d2aSXin Li ret = 1;
84*a58d3d2aSXin Li }
85*a58d3d2aSXin Li i+= i>>10;
86*a58d3d2aSXin Li }
87*a58d3d2aSXin Li }
88*a58d3d2aSXin Li
testbitexactcos(void)89*a58d3d2aSXin Li void testbitexactcos(void)
90*a58d3d2aSXin Li {
91*a58d3d2aSXin Li int i;
92*a58d3d2aSXin Li opus_int32 min_d,max_d,last,chk;
93*a58d3d2aSXin Li chk=max_d=0;
94*a58d3d2aSXin Li last=min_d=32767;
95*a58d3d2aSXin Li for(i=64;i<=16320;i++)
96*a58d3d2aSXin Li {
97*a58d3d2aSXin Li opus_int32 d;
98*a58d3d2aSXin Li opus_int32 q=bitexact_cos(i);
99*a58d3d2aSXin Li chk ^= q*i;
100*a58d3d2aSXin Li d = last - q;
101*a58d3d2aSXin Li if (d>max_d)max_d=d;
102*a58d3d2aSXin Li if (d<min_d)min_d=d;
103*a58d3d2aSXin Li last = q;
104*a58d3d2aSXin Li }
105*a58d3d2aSXin Li if ((chk!=89408644)||(max_d!=5)||(min_d!=0)||(bitexact_cos(64)!=32767)||
106*a58d3d2aSXin Li (bitexact_cos(16320)!=200)||(bitexact_cos(8192)!=23171))
107*a58d3d2aSXin Li {
108*a58d3d2aSXin Li fprintf (stderr, "bitexact_cos failed\n");
109*a58d3d2aSXin Li ret = 1;
110*a58d3d2aSXin Li }
111*a58d3d2aSXin Li }
112*a58d3d2aSXin Li
testbitexactlog2tan(void)113*a58d3d2aSXin Li void testbitexactlog2tan(void)
114*a58d3d2aSXin Li {
115*a58d3d2aSXin Li int i,fail;
116*a58d3d2aSXin Li opus_int32 min_d,max_d,last,chk;
117*a58d3d2aSXin Li fail=chk=max_d=0;
118*a58d3d2aSXin Li last=min_d=15059;
119*a58d3d2aSXin Li for(i=64;i<8193;i++)
120*a58d3d2aSXin Li {
121*a58d3d2aSXin Li opus_int32 d;
122*a58d3d2aSXin Li opus_int32 mid=bitexact_cos(i);
123*a58d3d2aSXin Li opus_int32 side=bitexact_cos(16384-i);
124*a58d3d2aSXin Li opus_int32 q=bitexact_log2tan(mid,side);
125*a58d3d2aSXin Li chk ^= q*i;
126*a58d3d2aSXin Li d = last - q;
127*a58d3d2aSXin Li if (q!=-1*bitexact_log2tan(side,mid))
128*a58d3d2aSXin Li fail = 1;
129*a58d3d2aSXin Li if (d>max_d)max_d=d;
130*a58d3d2aSXin Li if (d<min_d)min_d=d;
131*a58d3d2aSXin Li last = q;
132*a58d3d2aSXin Li }
133*a58d3d2aSXin Li if ((chk!=15821257)||(max_d!=61)||(min_d!=-2)||fail||
134*a58d3d2aSXin Li (bitexact_log2tan(32767,200)!=15059)||(bitexact_log2tan(30274,12540)!=2611)||
135*a58d3d2aSXin Li (bitexact_log2tan(23171,23171)!=0))
136*a58d3d2aSXin Li {
137*a58d3d2aSXin Li fprintf (stderr, "bitexact_log2tan failed\n");
138*a58d3d2aSXin Li ret = 1;
139*a58d3d2aSXin Li }
140*a58d3d2aSXin Li }
141*a58d3d2aSXin Li
142*a58d3d2aSXin Li #ifndef FIXED_POINT
testlog2(void)143*a58d3d2aSXin Li void testlog2(void)
144*a58d3d2aSXin Li {
145*a58d3d2aSXin Li float x;
146*a58d3d2aSXin Li for (x=0.001f;x<1677700.0;x+=(x/8.0))
147*a58d3d2aSXin Li {
148*a58d3d2aSXin Li float error = fabs((1.442695040888963387*log(x))-celt_log2(x));
149*a58d3d2aSXin Li if (error>0.0009)
150*a58d3d2aSXin Li {
151*a58d3d2aSXin Li fprintf (stderr, "celt_log2 failed: fabs((1.442695040888963387*log(x))-celt_log2(x))>0.001 (x = %f, error = %f)\n", x,error);
152*a58d3d2aSXin Li ret = 1;
153*a58d3d2aSXin Li }
154*a58d3d2aSXin Li }
155*a58d3d2aSXin Li }
156*a58d3d2aSXin Li
testexp2(void)157*a58d3d2aSXin Li void testexp2(void)
158*a58d3d2aSXin Li {
159*a58d3d2aSXin Li float x;
160*a58d3d2aSXin Li for (x=-11.0;x<24.0;x+=0.0007f)
161*a58d3d2aSXin Li {
162*a58d3d2aSXin Li float error = fabs(x-(1.442695040888963387*log(celt_exp2(x))));
163*a58d3d2aSXin Li if (error>0.0002)
164*a58d3d2aSXin Li {
165*a58d3d2aSXin Li fprintf (stderr, "celt_exp2 failed: fabs(x-(1.442695040888963387*log(celt_exp2(x))))>0.0005 (x = %f, error = %f)\n", x,error);
166*a58d3d2aSXin Li ret = 1;
167*a58d3d2aSXin Li }
168*a58d3d2aSXin Li }
169*a58d3d2aSXin Li }
170*a58d3d2aSXin Li
testexp2log2(void)171*a58d3d2aSXin Li void testexp2log2(void)
172*a58d3d2aSXin Li {
173*a58d3d2aSXin Li float x;
174*a58d3d2aSXin Li for (x=-11.0;x<24.0;x+=0.0007f)
175*a58d3d2aSXin Li {
176*a58d3d2aSXin Li float error = fabs(x-(celt_log2(celt_exp2(x))));
177*a58d3d2aSXin Li if (error>0.001)
178*a58d3d2aSXin Li {
179*a58d3d2aSXin Li fprintf (stderr, "celt_log2/celt_exp2 failed: fabs(x-(celt_log2(celt_exp2(x))))>0.001 (x = %f, error = %f)\n", x,error);
180*a58d3d2aSXin Li ret = 1;
181*a58d3d2aSXin Li }
182*a58d3d2aSXin Li }
183*a58d3d2aSXin Li }
184*a58d3d2aSXin Li #else
testlog2(void)185*a58d3d2aSXin Li void testlog2(void)
186*a58d3d2aSXin Li {
187*a58d3d2aSXin Li opus_val32 x;
188*a58d3d2aSXin Li for (x=8;x<1073741824;x+=(x>>3))
189*a58d3d2aSXin Li {
190*a58d3d2aSXin Li float error = fabs((1.442695040888963387*log(x/16384.0))-celt_log2(x)/1024.0);
191*a58d3d2aSXin Li if (error>0.003)
192*a58d3d2aSXin Li {
193*a58d3d2aSXin Li fprintf (stderr, "celt_log2 failed: x = %ld, error = %f\n", (long)x,error);
194*a58d3d2aSXin Li ret = 1;
195*a58d3d2aSXin Li }
196*a58d3d2aSXin Li }
197*a58d3d2aSXin Li }
198*a58d3d2aSXin Li
testexp2(void)199*a58d3d2aSXin Li void testexp2(void)
200*a58d3d2aSXin Li {
201*a58d3d2aSXin Li opus_val16 x;
202*a58d3d2aSXin Li for (x=-32768;x<15360;x++)
203*a58d3d2aSXin Li {
204*a58d3d2aSXin Li float error1 = fabs(x/1024.0-(1.442695040888963387*log(celt_exp2(x)/65536.0)));
205*a58d3d2aSXin Li float error2 = fabs(exp(0.6931471805599453094*x/1024.0)-celt_exp2(x)/65536.0);
206*a58d3d2aSXin Li if (error1>0.0002&&error2>0.00004)
207*a58d3d2aSXin Li {
208*a58d3d2aSXin Li fprintf (stderr, "celt_exp2 failed: x = "WORD", error1 = %f, error2 = %f\n", x,error1,error2);
209*a58d3d2aSXin Li ret = 1;
210*a58d3d2aSXin Li }
211*a58d3d2aSXin Li }
212*a58d3d2aSXin Li }
213*a58d3d2aSXin Li
testexp2log2(void)214*a58d3d2aSXin Li void testexp2log2(void)
215*a58d3d2aSXin Li {
216*a58d3d2aSXin Li opus_val32 x;
217*a58d3d2aSXin Li for (x=8;x<65536;x+=(x>>3))
218*a58d3d2aSXin Li {
219*a58d3d2aSXin Li float error = fabs(x-0.25*celt_exp2(celt_log2(x)))/16384;
220*a58d3d2aSXin Li if (error>0.004)
221*a58d3d2aSXin Li {
222*a58d3d2aSXin Li fprintf (stderr, "celt_log2/celt_exp2 failed: fabs(x-(celt_exp2(celt_log2(x))))>0.001 (x = %ld, error = %f)\n", (long)x,error);
223*a58d3d2aSXin Li ret = 1;
224*a58d3d2aSXin Li }
225*a58d3d2aSXin Li }
226*a58d3d2aSXin Li }
227*a58d3d2aSXin Li
testilog2(void)228*a58d3d2aSXin Li void testilog2(void)
229*a58d3d2aSXin Li {
230*a58d3d2aSXin Li opus_val32 x;
231*a58d3d2aSXin Li for (x=1;x<=268435455;x+=127)
232*a58d3d2aSXin Li {
233*a58d3d2aSXin Li opus_val32 lg;
234*a58d3d2aSXin Li opus_val32 y;
235*a58d3d2aSXin Li
236*a58d3d2aSXin Li lg = celt_ilog2(x);
237*a58d3d2aSXin Li if (lg<0 || lg>=31)
238*a58d3d2aSXin Li {
239*a58d3d2aSXin Li printf("celt_ilog2 failed: 0<=celt_ilog2(x)<31 (x = %d, celt_ilog2(x) = %d)\n",x,lg);
240*a58d3d2aSXin Li ret = 1;
241*a58d3d2aSXin Li }
242*a58d3d2aSXin Li y = 1<<lg;
243*a58d3d2aSXin Li
244*a58d3d2aSXin Li if (x<y || (x>>1)>=y)
245*a58d3d2aSXin Li {
246*a58d3d2aSXin Li printf("celt_ilog2 failed: 2**celt_ilog2(x)<=x<2**(celt_ilog2(x)+1) (x = %d, 2**celt_ilog2(x) = %d)\n",x,y);
247*a58d3d2aSXin Li ret = 1;
248*a58d3d2aSXin Li }
249*a58d3d2aSXin Li }
250*a58d3d2aSXin Li }
251*a58d3d2aSXin Li #endif
252*a58d3d2aSXin Li
main(void)253*a58d3d2aSXin Li int main(void)
254*a58d3d2aSXin Li {
255*a58d3d2aSXin Li testbitexactcos();
256*a58d3d2aSXin Li testbitexactlog2tan();
257*a58d3d2aSXin Li testdiv();
258*a58d3d2aSXin Li testsqrt();
259*a58d3d2aSXin Li testlog2();
260*a58d3d2aSXin Li testexp2();
261*a58d3d2aSXin Li testexp2log2();
262*a58d3d2aSXin Li #ifdef FIXED_POINT
263*a58d3d2aSXin Li testilog2();
264*a58d3d2aSXin Li #endif
265*a58d3d2aSXin Li return ret;
266*a58d3d2aSXin Li }
267