| /* |
| * Licensed to the Apache Software Foundation (ASF) under one or more |
| * contributor license agreements. See the NOTICE file distributed with |
| * this work for additional information regarding copyright ownership. |
| * The ASF licenses this file to You under the Apache License, Version 2.0 |
| * (the "License"); you may not use this file except in compliance with |
| * the License. You may obtain a copy of the License at |
| * |
| * http://www.apache.org/licenses/LICENSE-2.0 |
| * |
| * Unless required by applicable law or agreed to in writing, software |
| * distributed under the License is distributed on an "AS IS" BASIS, |
| * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| * See the License for the specific language governing permissions and |
| * limitations under the License. |
| */ |
| package org.apache.commons.math4.core.jdkmath; |
| |
| import java.io.PrintStream; |
| |
| |
| /** Class used to compute the classical functions tables. |
| * @since 3.0 |
| */ |
| final class AccurateMathCalc { |
| |
| /** |
| * 0x40000000 - used to split a double into two parts, both with the low order bits cleared. |
| * Equivalent to 2^30. |
| */ |
| private static final long HEX_40000000 = 0x40000000L; // 1073741824L |
| |
| /** Factorial table, for Taylor series expansions. 0!, 1!, 2!, ... 19! */ |
| private static final double[] FACT = new double[] |
| { |
| +1.0d, // 0 |
| +1.0d, // 1 |
| +2.0d, // 2 |
| +6.0d, // 3 |
| +24.0d, // 4 |
| +120.0d, // 5 |
| +720.0d, // 6 |
| +5040.0d, // 7 |
| +40320.0d, // 8 |
| +362880.0d, // 9 |
| +3628800.0d, // 10 |
| +39916800.0d, // 11 |
| +479001600.0d, // 12 |
| +6227020800.0d, // 13 |
| +87178291200.0d, // 14 |
| +1307674368000.0d, // 15 |
| +20922789888000.0d, // 16 |
| +355687428096000.0d, // 17 |
| +6402373705728000.0d, // 18 |
| +121645100408832000.0d, // 19 |
| }; |
| |
| /** Coefficients for slowLog. */ |
| private static final double[][] LN_SPLIT_COEF = { |
| {2.0, 0.0}, |
| {0.6666666269302368, 3.9736429850260626E-8}, |
| {0.3999999761581421, 2.3841857910019882E-8}, |
| {0.2857142686843872, 1.7029898543501842E-8}, |
| {0.2222222089767456, 1.3245471311735498E-8}, |
| {0.1818181574344635, 2.4384203044354907E-8}, |
| {0.1538461446762085, 9.140260083262505E-9}, |
| {0.13333332538604736, 9.220590270857665E-9}, |
| {0.11764700710773468, 1.2393345855018391E-8}, |
| {0.10526403784751892, 8.251545029714408E-9}, |
| {0.0952233225107193, 1.2675934823758863E-8}, |
| {0.08713622391223907, 1.1430250008909141E-8}, |
| {0.07842259109020233, 2.404307984052299E-9}, |
| {0.08371849358081818, 1.176342548272881E-8}, |
| {0.030589580535888672, 1.2958646899018938E-9}, |
| {0.14982303977012634, 1.225743062930824E-8}, |
| }; |
| |
| /** Table start declaration. */ |
| private static final String TABLE_START_DECL = " {"; |
| |
| /** Table end declaration. */ |
| private static final String TABLE_END_DECL = " };"; |
| |
| /** |
| * Private Constructor. |
| */ |
| private AccurateMathCalc() { |
| } |
| |
| /** Build the sine and cosine tables. |
| * @param SINE_TABLE_A table of the most significant part of the sines |
| * @param SINE_TABLE_B table of the least significant part of the sines |
| * @param COSINE_TABLE_A table of the most significant part of the cosines |
| * @param COSINE_TABLE_B table of the most significant part of the cosines |
| * @param SINE_TABLE_LEN length of the tables |
| * @param TANGENT_TABLE_A table of the most significant part of the tangents |
| * @param TANGENT_TABLE_B table of the most significant part of the tangents |
| */ |
| @SuppressWarnings("unused") |
| private static void buildSinCosTables(double[] SINE_TABLE_A, double[] SINE_TABLE_B, |
| double[] COSINE_TABLE_A, double[] COSINE_TABLE_B, |
| int SINE_TABLE_LEN, double[] TANGENT_TABLE_A, double[] TANGENT_TABLE_B) { |
| final double[] result = new double[2]; |
| |
| /* Use taylor series for 0 <= x <= 6/8 */ |
| for (int i = 0; i < 7; i++) { |
| double x = i / 8.0; |
| |
| slowSin(x, result); |
| SINE_TABLE_A[i] = result[0]; |
| SINE_TABLE_B[i] = result[1]; |
| |
| slowCos(x, result); |
| COSINE_TABLE_A[i] = result[0]; |
| COSINE_TABLE_B[i] = result[1]; |
| } |
| |
| /* Use angle addition formula to complete table to 13/8, just beyond pi/2 */ |
| for (int i = 7; i < SINE_TABLE_LEN; i++) { |
| double[] xs = new double[2]; |
| double[] ys = new double[2]; |
| double[] as = new double[2]; |
| double[] bs = new double[2]; |
| double[] temps = new double[2]; |
| |
| if ((i & 1) == 0) { |
| // Even, use double angle |
| xs[0] = SINE_TABLE_A[i / 2]; |
| xs[1] = SINE_TABLE_B[i / 2]; |
| ys[0] = COSINE_TABLE_A[i / 2]; |
| ys[1] = COSINE_TABLE_B[i / 2]; |
| |
| /* compute sine */ |
| splitMult(xs, ys, result); |
| SINE_TABLE_A[i] = result[0] * 2.0; |
| SINE_TABLE_B[i] = result[1] * 2.0; |
| |
| /* Compute cosine */ |
| splitMult(ys, ys, as); |
| splitMult(xs, xs, temps); |
| temps[0] = -temps[0]; |
| temps[1] = -temps[1]; |
| splitAdd(as, temps, result); |
| COSINE_TABLE_A[i] = result[0]; |
| COSINE_TABLE_B[i] = result[1]; |
| } else { |
| xs[0] = SINE_TABLE_A[i / 2]; |
| xs[1] = SINE_TABLE_B[i / 2]; |
| ys[0] = COSINE_TABLE_A[i / 2]; |
| ys[1] = COSINE_TABLE_B[i / 2]; |
| as[0] = SINE_TABLE_A[i / 2 + 1]; |
| as[1] = SINE_TABLE_B[i / 2 + 1]; |
| bs[0] = COSINE_TABLE_A[i / 2 + 1]; |
| bs[1] = COSINE_TABLE_B[i / 2 + 1]; |
| |
| /* compute sine */ |
| splitMult(xs, bs, temps); |
| splitMult(ys, as, result); |
| splitAdd(result, temps, result); |
| SINE_TABLE_A[i] = result[0]; |
| SINE_TABLE_B[i] = result[1]; |
| |
| /* Compute cosine */ |
| splitMult(ys, bs, result); |
| splitMult(xs, as, temps); |
| temps[0] = -temps[0]; |
| temps[1] = -temps[1]; |
| splitAdd(result, temps, result); |
| COSINE_TABLE_A[i] = result[0]; |
| COSINE_TABLE_B[i] = result[1]; |
| } |
| } |
| |
| /* Compute tangent = sine/cosine */ |
| for (int i = 0; i < SINE_TABLE_LEN; i++) { |
| double[] xs = new double[2]; |
| double[] ys = new double[2]; |
| double[] as = new double[2]; |
| |
| as[0] = COSINE_TABLE_A[i]; |
| as[1] = COSINE_TABLE_B[i]; |
| |
| splitReciprocal(as, ys); |
| |
| xs[0] = SINE_TABLE_A[i]; |
| xs[1] = SINE_TABLE_B[i]; |
| |
| splitMult(xs, ys, as); |
| |
| TANGENT_TABLE_A[i] = as[0]; |
| TANGENT_TABLE_B[i] = as[1]; |
| } |
| } |
| |
| /** |
| * For x between 0 and pi/4 compute cosine using Talor series |
| * cos(x) = 1 - x^2/2! + x^4/4! ... |
| * @param x number from which cosine is requested |
| * @param result placeholder where to put the result in extended precision |
| * (may be null) |
| * @return cos(x) |
| */ |
| static double slowCos(final double x, final double[] result) { |
| |
| final double[] xs = new double[2]; |
| final double[] ys = new double[2]; |
| final double[] facts = new double[2]; |
| final double[] as = new double[2]; |
| split(x, xs); |
| ys[0] = ys[1] = 0.0; |
| |
| for (int i = FACT.length - 1; i >= 0; i--) { |
| splitMult(xs, ys, as); |
| ys[0] = as[0]; |
| ys[1] = as[1]; |
| |
| if ((i & 1) != 0) { // skip odd entries |
| continue; |
| } |
| |
| split(FACT[i], as); |
| splitReciprocal(as, facts); |
| |
| if ((i & 2) != 0) { // alternate terms are negative |
| facts[0] = -facts[0]; |
| facts[1] = -facts[1]; |
| } |
| |
| splitAdd(ys, facts, as); |
| ys[0] = as[0]; ys[1] = as[1]; |
| } |
| |
| if (result != null) { |
| result[0] = ys[0]; |
| result[1] = ys[1]; |
| } |
| |
| return ys[0] + ys[1]; |
| } |
| |
| /** |
| * For x between 0 and pi/4 compute sine using Taylor expansion: |
| * sin(x) = x - x^3/3! + x^5/5! - x^7/7! ... |
| * @param x number from which sine is requested |
| * @param result placeholder where to put the result in extended precision |
| * (may be null) |
| * @return sin(x) |
| */ |
| static double slowSin(final double x, final double[] result) { |
| final double[] xs = new double[2]; |
| final double[] ys = new double[2]; |
| final double[] facts = new double[2]; |
| final double[] as = new double[2]; |
| split(x, xs); |
| ys[0] = ys[1] = 0.0; |
| |
| for (int i = FACT.length - 1; i >= 0; i--) { |
| splitMult(xs, ys, as); |
| ys[0] = as[0]; |
| ys[1] = as[1]; |
| |
| if ((i & 1) == 0) { // Ignore even numbers |
| continue; |
| } |
| |
| split(FACT[i], as); |
| splitReciprocal(as, facts); |
| |
| if ((i & 2) != 0) { // alternate terms are negative |
| facts[0] = -facts[0]; |
| facts[1] = -facts[1]; |
| } |
| |
| splitAdd(ys, facts, as); |
| ys[0] = as[0]; ys[1] = as[1]; |
| } |
| |
| if (result != null) { |
| result[0] = ys[0]; |
| result[1] = ys[1]; |
| } |
| |
| return ys[0] + ys[1]; |
| } |
| |
| |
| /** |
| * For x between 0 and 1, returns exp(x), uses extended precision. |
| * @param x argument of exponential |
| * @param result placeholder where to place exp(x) split in two terms |
| * for extra precision (i.e. exp(x) = result[0] + result[1] |
| * @return exp(x) |
| */ |
| static double slowexp(final double x, final double[] result) { |
| final double[] xs = new double[2]; |
| final double[] ys = new double[2]; |
| final double[] facts = new double[2]; |
| final double[] as = new double[2]; |
| split(x, xs); |
| ys[0] = ys[1] = 0.0; |
| |
| for (int i = FACT.length - 1; i >= 0; i--) { |
| splitMult(xs, ys, as); |
| ys[0] = as[0]; |
| ys[1] = as[1]; |
| |
| split(FACT[i], as); |
| splitReciprocal(as, facts); |
| |
| splitAdd(ys, facts, as); |
| ys[0] = as[0]; |
| ys[1] = as[1]; |
| } |
| |
| if (result != null) { |
| result[0] = ys[0]; |
| result[1] = ys[1]; |
| } |
| |
| return ys[0] + ys[1]; |
| } |
| |
| /** Compute split[0], split[1] such that their sum is equal to d, |
| * and split[0] has its 30 least significant bits as zero. |
| * @param d number to split |
| * @param split placeholder where to place the result |
| */ |
| private static void split(final double d, final double[] split) { |
| if (d < 8e298 && d > -8e298) { |
| final double a = d * HEX_40000000; |
| split[0] = (d + a) - a; |
| split[1] = d - split[0]; |
| } else { |
| final double a = d * 9.31322574615478515625E-10; |
| split[0] = (d + a - d) * HEX_40000000; |
| split[1] = d - split[0]; |
| } |
| } |
| |
| /** Recompute a split. |
| * @param a input/out array containing the split, changed |
| * on output |
| */ |
| private static void resplit(final double[] a) { |
| final double c = a[0] + a[1]; |
| final double d = -(c - a[0] - a[1]); |
| |
| if (c < 8e298 && c > -8e298) { // MAGIC NUMBER |
| double z = c * HEX_40000000; |
| a[0] = (c + z) - z; |
| a[1] = c - a[0] + d; |
| } else { |
| double z = c * 9.31322574615478515625E-10; |
| a[0] = (c + z - c) * HEX_40000000; |
| a[1] = c - a[0] + d; |
| } |
| } |
| |
| /** Multiply two numbers in split form. |
| * @param a first term of multiplication |
| * @param b second term of multiplication |
| * @param ans placeholder where to put the result |
| */ |
| private static void splitMult(double[] a, double[] b, double[] ans) { |
| ans[0] = a[0] * b[0]; |
| ans[1] = a[0] * b[1] + a[1] * b[0] + a[1] * b[1]; |
| |
| /* Resplit */ |
| resplit(ans); |
| } |
| |
| /** Add two numbers in split form. |
| * @param a first term of addition |
| * @param b second term of addition |
| * @param ans placeholder where to put the result |
| */ |
| private static void splitAdd(final double[] a, final double[] b, final double[] ans) { |
| ans[0] = a[0] + b[0]; |
| ans[1] = a[1] + b[1]; |
| |
| resplit(ans); |
| } |
| |
| /** Compute the reciprocal of in. Use the following algorithm. |
| * in = c + d. |
| * want to find x + y such that x+y = 1/(c+d) and x is much |
| * larger than y and x has several zero bits on the right. |
| * |
| * Set b = 1/(2^22), a = 1 - b. Thus (a+b) = 1. |
| * Use following identity to compute (a+b)/(c+d) |
| * |
| * (a+b)/(c+d) = a/c + (bc - ad) / (c^2 + cd) |
| * set x = a/c and y = (bc - ad) / (c^2 + cd) |
| * This will be close to the right answer, but there will be |
| * some rounding in the calculation of X. So by carefully |
| * computing 1 - (c+d)(x+y) we can compute an error and |
| * add that back in. This is done carefully so that terms |
| * of similar size are subtracted first. |
| * @param in initial number, in split form |
| * @param result placeholder where to put the result |
| */ |
| static void splitReciprocal(final double[] in, final double[] result) { |
| final double b = 1.0 / 4194304.0; |
| final double a = 1.0 - b; |
| |
| if (in[0] == 0.0) { |
| in[0] = in[1]; |
| in[1] = 0.0; |
| } |
| |
| result[0] = a / in[0]; |
| result[1] = (b * in[0] - a * in[1]) / (in[0] * in[0] + in[0] * in[1]); |
| |
| if (result[1] != result[1]) { // can happen if result[1] is NAN |
| result[1] = 0.0; |
| } |
| |
| /* Resplit */ |
| resplit(result); |
| |
| for (int i = 0; i < 2; i++) { |
| /* this may be overkill, probably once is enough */ |
| double err = 1.0 - result[0] * in[0] - result[0] * in[1] - |
| result[1] * in[0] - result[1] * in[1]; |
| /*err = 1.0 - err; */ |
| err *= result[0] + result[1]; |
| /*printf("err = %16e\n", err); */ |
| result[1] += err; |
| } |
| } |
| |
| /** Compute (a[0] + a[1]) * (b[0] + b[1]) in extended precision. |
| * @param a first term of the multiplication |
| * @param b second term of the multiplication |
| * @param result placeholder where to put the result |
| */ |
| private static void quadMult(final double[] a, final double[] b, final double[] result) { |
| final double[] xs = new double[2]; |
| final double[] ys = new double[2]; |
| final double[] zs = new double[2]; |
| |
| /* a[0] * b[0] */ |
| split(a[0], xs); |
| split(b[0], ys); |
| splitMult(xs, ys, zs); |
| |
| result[0] = zs[0]; |
| result[1] = zs[1]; |
| |
| /* a[0] * b[1] */ |
| split(b[1], ys); |
| splitMult(xs, ys, zs); |
| |
| double tmp = result[0] + zs[0]; |
| result[1] -= tmp - result[0] - zs[0]; |
| result[0] = tmp; |
| tmp = result[0] + zs[1]; |
| result[1] -= tmp - result[0] - zs[1]; |
| result[0] = tmp; |
| |
| /* a[1] * b[0] */ |
| split(a[1], xs); |
| split(b[0], ys); |
| splitMult(xs, ys, zs); |
| |
| tmp = result[0] + zs[0]; |
| result[1] -= tmp - result[0] - zs[0]; |
| result[0] = tmp; |
| tmp = result[0] + zs[1]; |
| result[1] -= tmp - result[0] - zs[1]; |
| result[0] = tmp; |
| |
| /* a[1] * b[0] */ |
| split(a[1], xs); |
| split(b[1], ys); |
| splitMult(xs, ys, zs); |
| |
| tmp = result[0] + zs[0]; |
| result[1] -= tmp - result[0] - zs[0]; |
| result[0] = tmp; |
| tmp = result[0] + zs[1]; |
| result[1] -= tmp - result[0] - zs[1]; |
| result[0] = tmp; |
| } |
| |
| /** Compute exp(p) for a integer p in extended precision. |
| * @param p integer whose exponential is requested |
| * @param result placeholder where to put the result in extended precision |
| * @return exp(p) in standard precision (equal to result[0] + result[1]) |
| */ |
| static double expint(int p, final double[] result) { |
| //double x = M_E; |
| final double[] xs = new double[2]; |
| final double[] as = new double[2]; |
| final double[] ys = new double[2]; |
| //split(x, xs); |
| //xs[1] = (double)(2.7182818284590452353602874713526625L - xs[0]); |
| //xs[0] = 2.71827697753906250000; |
| //xs[1] = 4.85091998273542816811e-06; |
| //xs[0] = Double.longBitsToDouble(0x4005bf0800000000L); |
| //xs[1] = Double.longBitsToDouble(0x3ed458a2bb4a9b00L); |
| |
| /* E */ |
| xs[0] = 2.718281828459045; |
| xs[1] = 1.4456468917292502E-16; |
| |
| split(1.0, ys); |
| |
| while (p > 0) { |
| if ((p & 1) != 0) { |
| quadMult(ys, xs, as); |
| ys[0] = as[0]; ys[1] = as[1]; |
| } |
| |
| quadMult(xs, xs, as); |
| xs[0] = as[0]; xs[1] = as[1]; |
| |
| p >>= 1; |
| } |
| |
| if (result != null) { |
| result[0] = ys[0]; |
| result[1] = ys[1]; |
| |
| resplit(result); |
| } |
| |
| return ys[0] + ys[1]; |
| } |
| /** xi in the range of [1, 2]. |
| * 3 5 7 |
| * x+1 / x x x \ |
| * ln ----- = 2 * | x + ---- + ---- + ---- + ... | |
| * 1-x \ 3 5 7 / |
| * |
| * So, compute a Remez approximation of the following function |
| * |
| * ln ((sqrt(x)+1)/(1-sqrt(x))) / x |
| * |
| * This will be an even function with only positive coefficents. |
| * x is in the range [0 - 1/3]. |
| * |
| * Transform xi for input to the above function by setting |
| * x = (xi-1)/(xi+1). Input to the polynomial is x^2, then |
| * the result is multiplied by x. |
| * @param xi number from which log is requested |
| * @return log(xi) |
| */ |
| static double[] slowLog(double xi) { |
| double[] x = new double[2]; |
| double[] x2 = new double[2]; |
| double[] y = new double[2]; |
| double[] a = new double[2]; |
| |
| split(xi, x); |
| |
| /* Set X = (x-1)/(x+1) */ |
| x[0] += 1.0; |
| resplit(x); |
| splitReciprocal(x, a); |
| x[0] -= 2.0; |
| resplit(x); |
| splitMult(x, a, y); |
| x[0] = y[0]; |
| x[1] = y[1]; |
| |
| /* Square X -> X2*/ |
| splitMult(x, x, x2); |
| |
| |
| //x[0] -= 1.0; |
| //resplit(x); |
| |
| y[0] = LN_SPLIT_COEF[LN_SPLIT_COEF.length - 1][0]; |
| y[1] = LN_SPLIT_COEF[LN_SPLIT_COEF.length - 1][1]; |
| |
| for (int i = LN_SPLIT_COEF.length - 2; i >= 0; i--) { |
| splitMult(y, x2, a); |
| y[0] = a[0]; |
| y[1] = a[1]; |
| splitAdd(y, LN_SPLIT_COEF[i], a); |
| y[0] = a[0]; |
| y[1] = a[1]; |
| } |
| |
| splitMult(y, x, a); |
| y[0] = a[0]; |
| y[1] = a[1]; |
| |
| return y; |
| } |
| |
| |
| /** |
| * Print an array. |
| * @param out text output stream where output should be printed |
| * @param name array name |
| * @param expectedLen expected length of the array |
| * @param array2d array data |
| */ |
| static void printarray(PrintStream out, String name, int expectedLen, double[][] array2d) { |
| out.println(name); |
| checkLen(expectedLen, array2d.length); |
| out.println(TABLE_START_DECL + " "); |
| int i = 0; |
| for (double[] array : array2d) { // "double array[]" causes PMD parsing error |
| out.print(" {"); |
| for (double d : array) { // assume inner array has very few entries |
| out.printf("%-25.25s", format(d)); // multiple entries per line |
| } |
| out.println("}, // " + i++); |
| } |
| out.println(TABLE_END_DECL); |
| } |
| |
| /** |
| * Print an array. |
| * @param out text output stream where output should be printed |
| * @param name array name |
| * @param expectedLen expected length of the array |
| * @param array array data |
| */ |
| static void printarray(PrintStream out, String name, int expectedLen, double[] array) { |
| out.println(name + "="); |
| checkLen(expectedLen, array.length); |
| out.println(TABLE_START_DECL); |
| for (double d : array) { |
| out.printf(" %s%n", format(d)); // one entry per line |
| } |
| out.println(TABLE_END_DECL); |
| } |
| |
| /** Format a double. |
| * @param d double number to format |
| * @return formatted number |
| */ |
| static String format(double d) { |
| if (Double.isNaN(d)) { |
| return "Double.NaN,"; |
| } else { |
| return ((d >= 0) ? "+" : "") + Double.toString(d) + "d,"; |
| } |
| } |
| |
| /** |
| * Check two lengths are equal. |
| * @param expectedLen expected length |
| * @param actual actual length |
| * @throws IllegalStateException if the two lengths are not equal |
| */ |
| private static void checkLen(int expectedLen, int actual) { |
| if (expectedLen != actual) { |
| throw new IllegalStateException(actual + " != " + expectedLen); |
| } |
| } |
| } |
