1/*
2 * Copyright (C) 2015 The Android Open Source Project
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 *   http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 */
16
17// THIS TYPICALLY TAKES > 4 MINUTES TO RUN!
18// It should generate no output during that time.
19
20package com.hp.creals;
21
22import junit.framework.AssertionFailedError;
23import junit.framework.TestCase;
24
25import java.math.BigInteger;
26import java.util.Random;
27
28public class SlowCRTest extends TestCase {
29    private static void check(boolean x, String s) {
30        if (!x) throw new AssertionFailedError(s);
31    }
32    final static int TEST_PREC = -200; // 200 bits to the right of
33                                       // binary point.
34    final static int NRANDOM = 100;    // Number of random values to
35                                       // test.  Bigger ==> slower
36    private static void checkEq(CR x, CR y, String s) {
37        check(x.compareTo(y, TEST_PREC) == 0, s);
38    }
39    private static void checkApprEq(double x, double y, String s) {
40        check(Math.abs(x - y) <= 0.000001, s);
41    }
42    final static BigInteger MASK =
43            BigInteger.ONE.shiftLeft(-TEST_PREC).subtract(BigInteger.ONE);
44    private static boolean isApprInt(CR x) {
45        BigInteger appr = x.get_appr(TEST_PREC);
46        return appr.and(MASK).signum() == 0;
47    }
48
49    final static CR ZERO = CR.valueOf(0);
50    final static CR ONE = CR.valueOf(1);
51    final static CR TWO = CR.valueOf(2);
52    final static CR BIG = CR.valueOf(200).exp();
53    final static CR SMALL = BIG.inverse();
54    final static CR HALF_PI = CR.PI.divide(CR.valueOf(2));
55
56    final static UnaryCRFunction ATAN = UnaryCRFunction.atanFunction;
57    final static UnaryCRFunction TAN = UnaryCRFunction.tanFunction;
58    final static UnaryCRFunction COSINE = UnaryCRFunction.sinFunction
59                                             .monotoneDerivative(ZERO, CR.PI);
60    final static UnaryCRFunction ARCSINE =
61                 UnaryCRFunction.sinFunction.inverseMonotone(HALF_PI.negate(),
62                                                             HALF_PI);
63
64    // Perform some consistency checks on trig functions at x
65    // We assume that x is within floating point range.
66    private static void checkTrig(CR x) {
67        double xAsDouble = x.doubleValue();
68        if (Math.abs(xAsDouble) < 1000000.0) {
69            checkApprEq(x.sin().doubleValue(), Math.sin(xAsDouble),
70                          "sin float compare:" + xAsDouble);
71            checkApprEq(x.cos().doubleValue(), Math.cos(xAsDouble),
72                          "cos float compare:" + xAsDouble);
73            checkApprEq(TAN.execute(x).doubleValue(), Math.tan(xAsDouble),
74                          "tan float compare:" + xAsDouble);
75            checkApprEq(ATAN.execute(x).doubleValue(), Math.atan(xAsDouble),
76                          "atan float compare:" + xAsDouble);
77        }
78        if (Math.abs(xAsDouble) < 1.0) {
79            checkApprEq(x.asin().doubleValue(), Math.asin(xAsDouble),
80                          "asin float compare:" + xAsDouble);
81            checkApprEq(x.acos().doubleValue(), Math.acos(xAsDouble),
82                          "acos float compare:" + xAsDouble);
83            checkEq(ARCSINE.execute(x), x.asin(),
84                          "inverse(sin) compare:" + xAsDouble);
85        }
86        if (xAsDouble < 3.1415926535 && xAsDouble > 0.0) {
87            checkApprEq(COSINE.execute(x).doubleValue(), Math.cos(xAsDouble),
88                          "deriv(sin) float compare:"  + xAsDouble);
89            checkEq(COSINE.execute(x), x.cos(),
90                          "deriv(sin) float compare:"  + xAsDouble);
91        }
92        // Check that sin(x+v) = sin(x)cos(v) + cos(x)sin(v)
93        // for a couple of different values of v.
94        for (int i = 1; i <= 5; ++i) {
95            CR v = CR.valueOf(i);
96            checkEq(
97                x.add(v).sin(),
98                x.sin().multiply(v.cos()).add(x.cos().multiply(v.sin())),
99                "Angle sum formula failed for " + xAsDouble + " + " + i);
100        }
101        checkEq(x.cos().multiply(x.cos()).add(x.sin().multiply(x.sin())),
102                    CR.valueOf(1),
103                    "sin(x)^2 + cos(x)^2 != 1:" + xAsDouble);
104        // Check that inverses are consistent
105        checkEq(x, TAN.execute(ATAN.execute(x)),
106                      "tan(atan(" + xAsDouble + ")");
107        CR xcos = x.cos();
108        CR tmp = xcos.acos();
109        // Result or its inverse should differ from x by an
110        // exact multiple of pi.
111        check(isApprInt(tmp.subtract(x).divide(CR.PI))
112              || isApprInt(tmp.add(x).divide(CR.PI)),
113              "acos(cos):" + xAsDouble);
114        CR xsin = x.sin();
115        tmp = ARCSINE.execute(xsin);
116        CR tmp2 = xsin.asin();
117        checkEq(tmp, tmp2, "Asin(sin) computations differ:" + xAsDouble);
118        // Result or its inverse should differ from x by an
119        // exact multiple of pi.
120        check(isApprInt(tmp.subtract(x).divide(CR.PI))
121              || isApprInt(tmp.add(x).divide(CR.PI)),
122              "acos(cos):" + xAsDouble);
123    }
124
125    private static void checkExpLn(CR x) {
126        double xAsDouble = x.doubleValue();
127        if (Math.abs(xAsDouble) < 10.0) {
128            checkApprEq(x.exp().doubleValue(), Math.exp(xAsDouble),
129                          "exp float compare:" + xAsDouble);
130        }
131        if (Math.abs(xAsDouble) <= 1000.0) {
132            checkEq(x, x.exp().ln(), "ln(exp) failed:" + xAsDouble);
133            checkEq(x.multiply(CR.valueOf(2)).exp(),
134
135                        x.exp().multiply(x.exp()),
136                        "exp^2 failed:" + xAsDouble);
137        }
138        if (xAsDouble > 0.000000001) {
139            checkApprEq(x.ln().doubleValue(), Math.log(xAsDouble),
140                          "exp float compare:" + xAsDouble);
141            checkEq(x, x.ln().exp(), "exp(ln) failed:" + xAsDouble);
142            checkEq(x.ln().divide(CR.valueOf(2)), x.sqrt().ln(),
143                        "ln(sqrt) failed:" + xAsDouble);
144            // Check that ln(xv) = ln(x) + ln(v) for various v
145            for (int i = 1; i <= 5; ++i) {
146                CR v = CR.valueOf(i);
147                checkEq(
148                    x.ln().add(v.ln()),
149                    x.multiply(v).ln(),
150                    "ln(product) formula failed for:" + xAsDouble + "," + i);
151            }
152        }
153    }
154
155    private static void checkBasic(CR x) {
156        checkEq(x.abs().sqrt().multiply(x.abs().sqrt()), x.abs(),
157                    "sqrt*sqrt:" + x.doubleValue());
158        if (!x.get_appr(TEST_PREC).equals(BigInteger.ZERO)) {
159            checkEq(x.inverse().inverse(), x,
160                        "inverse(inverse):" + x.doubleValue());
161        }
162    }
163
164    public void testSlowTrig() {
165        checkEq(ZERO.acos(), CR.PI.divide(TWO), "acos(0)");
166        checkEq(ONE.acos(), ZERO, "acos(1)");
167        checkEq(ONE.negate().acos(), CR.PI, "acos(-1)");
168        checkEq(ZERO.asin(), ZERO, "asin(0)");
169        checkEq(ONE.asin(), CR.PI.divide(TWO), "asin(1)");
170        checkEq(ONE.negate().asin(), CR.PI.divide(TWO).negate(), "asin(-1)");
171        checkTrig(ZERO);
172        CR BIG = CR.valueOf(200).exp();
173        checkTrig(BIG);
174        checkTrig(BIG.negate());
175        checkTrig(SMALL);
176        checkTrig(SMALL.negate());
177        checkTrig(CR.PI);
178        checkTrig(CR.PI.subtract(SMALL));
179        checkTrig(CR.PI.add(SMALL));
180        checkTrig(CR.PI.negate());
181        checkTrig(CR.PI.negate().subtract(SMALL));
182        checkTrig(CR.PI.negate().add(SMALL));
183        Random r = new Random();  // Random seed!
184        for (int i = 0; i < NRANDOM; ++i) {
185            double d = Math.exp(2.0 * r.nextDouble() - 1.0);
186            if (r.nextBoolean()) d = -d;
187            final CR x = CR.valueOf(d);
188            checkTrig(x);
189        }
190        // And a few big ones
191        for (int i = 0; i < 10; ++i) {
192            double d = Math.exp(200.0 * r.nextDouble());
193            if (r.nextBoolean()) d = -d;
194            final CR x = CR.valueOf(d);
195            checkTrig(x);
196        }
197    }
198
199    public void testSlowExpLn() {
200        checkEq(CR.valueOf(1).ln(), CR.valueOf(0), "ln(1) != 0");
201        checkExpLn(CR.valueOf(0));
202        CR BIG = CR.valueOf(200).exp();
203        checkExpLn(BIG);
204        checkExpLn(BIG.negate());
205        checkExpLn(SMALL);
206        checkExpLn(SMALL.negate());
207        checkExpLn(CR.PI);
208        checkExpLn(ONE);
209        checkExpLn(ONE.subtract(SMALL));
210        checkExpLn(ONE.negate().subtract(SMALL));
211        Random r = new Random();  // Random seed!
212        for (int i = 0; i < NRANDOM; ++i) {
213            double d = Math.exp(10.0 * r.nextDouble() - 1.0);
214            if (r.nextBoolean()) d = -d;
215            final CR x = CR.valueOf(d);
216            checkExpLn(x);
217        }
218        // And a few big ones
219        for (int i = 0; i < 10; ++i) {
220            double d = Math.exp(200.0 * r.nextDouble());
221            if (r.nextBoolean()) d = -d;
222            final CR x = CR.valueOf(d);
223            checkExpLn(x);
224        }
225    }
226
227    public void testSlowBasic() {
228        checkEq(ZERO.sqrt(), ZERO, "sqrt(0)");
229        checkEq(ZERO.abs(), ZERO, "abs(0)");
230        Random r = new Random();  // Random seed!
231        for (int i = 0; i < NRANDOM; ++i) {
232            double d = Math.exp(10.0 * r.nextDouble() - 1.0);
233            if (r.nextBoolean()) d = -d;
234            final CR x = CR.valueOf(d);
235            checkBasic(x);
236        }
237        // And a few very big ones, but within IEEE double range
238        for (int i = 0; i < 10; ++i) {
239            double d = Math.exp(600.0 * r.nextDouble());
240            if (r.nextBoolean()) d = -d;
241            final CR x = CR.valueOf(d);
242            checkBasic(x);
243        }
244    }
245}
246