Delete files that were added to the repository by mistake.
This commit is contained in:
Gilles Sadowski
parent
a5e0b2e7f6
commit
44d1f03669
File diff suppressed because it is too large
Load Diff
|
|
@ -1,658 +0,0 @@
|
|||
/*
|
||||
* 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.legacy.util;
|
||||
|
||||
import java.io.PrintStream;
|
||||
|
||||
import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
|
||||
|
||||
/** Class used to compute the classical functions tables.
|
||||
* @since 3.0
|
||||
*/
|
||||
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 (d != 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
|
||||
* @exception DimensionMismatchException if the two lengths are not equal
|
||||
*/
|
||||
private static void checkLen(int expectedLen, int actual)
|
||||
throws DimensionMismatchException {
|
||||
if (expectedLen != actual) {
|
||||
throw new DimensionMismatchException(actual, expectedLen);
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
|
|
@ -1,262 +0,0 @@
|
|||
/*
|
||||
* 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.legacy.util;
|
||||
|
||||
import java.lang.reflect.InvocationTargetException;
|
||||
import java.lang.reflect.Method;
|
||||
import java.lang.reflect.Modifier;
|
||||
import java.lang.reflect.Type;
|
||||
import java.util.ArrayList;
|
||||
import java.util.Arrays;
|
||||
import java.util.List;
|
||||
|
||||
import org.apache.commons.math4.legacy.exception.MathArithmeticException;
|
||||
import org.apache.commons.numbers.core.Precision;
|
||||
import org.junit.Assert;
|
||||
import org.junit.Test;
|
||||
import org.junit.runner.RunWith;
|
||||
import org.junit.runners.Parameterized;
|
||||
import org.junit.runners.Parameterized.Parameters;
|
||||
|
||||
/**
|
||||
* Test to compare AccurateMath results against StrictMath results for boundary values.
|
||||
* <p>
|
||||
* Running all tests independently: <br>
|
||||
* {@code mvn test -Dtest=AccurateMathStrictComparisonTest}<br>
|
||||
* or just run tests against a single method (e.g. scalb):<br>
|
||||
* {@code mvn test -Dtest=AccurateMathStrictComparisonTest -DargLine="-DtestMethod=scalb"}
|
||||
*/
|
||||
@SuppressWarnings("boxing")
|
||||
@RunWith(Parameterized.class)
|
||||
public class AccurateMathStrictComparisonTest {
|
||||
|
||||
// Values which often need special handling
|
||||
private static final Double[] DOUBLE_SPECIAL_VALUES = {
|
||||
-0.0, +0.0, // 1,2
|
||||
Double.NaN, // 3
|
||||
Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY, // 4,5
|
||||
-Double.MAX_VALUE, Double.MAX_VALUE, // 6,7
|
||||
// decreasing order of absolute value to help catch first failure
|
||||
-Precision.EPSILON, Precision.EPSILON, // 8,9
|
||||
-Precision.SAFE_MIN, Precision.SAFE_MIN, // 10,11
|
||||
-Double.MIN_VALUE, Double.MIN_VALUE, // 12,13
|
||||
};
|
||||
|
||||
private static final Float [] FLOAT_SPECIAL_VALUES = {
|
||||
-0.0f, +0.0f, // 1,2
|
||||
Float.NaN, // 3
|
||||
Float.NEGATIVE_INFINITY, Float.POSITIVE_INFINITY, // 4,5
|
||||
Float.MIN_VALUE, Float.MAX_VALUE, // 6,7
|
||||
-Float.MIN_VALUE, -Float.MAX_VALUE, // 8,9
|
||||
};
|
||||
|
||||
private static final Object [] LONG_SPECIAL_VALUES = {
|
||||
-1,0,1, // 1,2,3
|
||||
Long.MIN_VALUE, Long.MAX_VALUE, // 4,5
|
||||
};
|
||||
|
||||
private static final Object[] INT_SPECIAL_VALUES = {
|
||||
-1,0,1, // 1,2,3
|
||||
Integer.MIN_VALUE, Integer.MAX_VALUE, // 4,5
|
||||
};
|
||||
|
||||
private final Method mathMethod;
|
||||
private final Method fastMethod;
|
||||
private final Type[] types;
|
||||
private final Object[][] valueArrays;
|
||||
|
||||
public AccurateMathStrictComparisonTest(Method m, Method f, Type[] types, Object[][] data) throws Exception{
|
||||
this.mathMethod=m;
|
||||
this.fastMethod=f;
|
||||
this.types=types;
|
||||
this.valueArrays=data;
|
||||
}
|
||||
|
||||
@Test
|
||||
public void test1() throws Exception{
|
||||
setupMethodCall(mathMethod, fastMethod, types, valueArrays);
|
||||
}
|
||||
private static boolean isNumber(Double d) {
|
||||
return !(d.isInfinite() || d.isNaN());
|
||||
}
|
||||
|
||||
private static boolean isNumber(Float f) {
|
||||
return !(f.isInfinite() || f.isNaN());
|
||||
}
|
||||
|
||||
private static void reportFailedResults(Method mathMethod, Object[] params, Object expected, Object actual, int[] entries){
|
||||
final String methodName = mathMethod.getName();
|
||||
String format = null;
|
||||
long actL=0;
|
||||
long expL=0;
|
||||
if (expected instanceof Double) {
|
||||
Double exp = (Double) expected;
|
||||
Double act = (Double) actual;
|
||||
if (isNumber(exp) && isNumber(act) && exp != 0) { // show difference as hex
|
||||
actL = Double.doubleToLongBits(act);
|
||||
expL = Double.doubleToLongBits(exp);
|
||||
if (Math.abs(actL-expL)==1) {
|
||||
// Not 100% sure off-by-one errors are allowed everywhere, so only allow for these methods
|
||||
if (methodName.equals("toRadians") || methodName.equals("atan2")) {
|
||||
return;
|
||||
}
|
||||
}
|
||||
format = "%016x";
|
||||
}
|
||||
} else if (expected instanceof Float ){
|
||||
Float exp = (Float) expected;
|
||||
Float act = (Float) actual;
|
||||
if (isNumber(exp) && isNumber(act) && exp != 0) { // show difference as hex
|
||||
actL = Float.floatToIntBits(act);
|
||||
expL = Float.floatToIntBits(exp);
|
||||
format = "%08x";
|
||||
}
|
||||
}
|
||||
StringBuilder sb = new StringBuilder();
|
||||
sb.append(mathMethod.getReturnType().getSimpleName());
|
||||
sb.append(" ");
|
||||
sb.append(methodName);
|
||||
sb.append("(");
|
||||
String sep = "";
|
||||
for(Object o : params){
|
||||
sb.append(sep);
|
||||
sb.append(o);
|
||||
sep=", ";
|
||||
}
|
||||
sb.append(") expected ");
|
||||
if (format != null){
|
||||
sb.append(String.format(format, expL));
|
||||
} else {
|
||||
sb.append(expected);
|
||||
}
|
||||
sb.append(" actual ");
|
||||
if (format != null){
|
||||
sb.append(String.format(format, actL));
|
||||
} else {
|
||||
sb.append(actual);
|
||||
}
|
||||
sb.append(" entries ");
|
||||
sb.append(Arrays.toString(entries));
|
||||
String message = sb.toString();
|
||||
final boolean fatal = true;
|
||||
if (fatal) {
|
||||
Assert.fail(message);
|
||||
} else {
|
||||
System.out.println(message);
|
||||
}
|
||||
}
|
||||
|
||||
private static void callMethods(Method mathMethod, Method fastMethod,
|
||||
Object[] params, int[] entries) throws IllegalAccessException {
|
||||
try {
|
||||
Object expected;
|
||||
try {
|
||||
expected = mathMethod.invoke(mathMethod, params);
|
||||
} catch (InvocationTargetException ite) {
|
||||
expected = ite.getCause();
|
||||
}
|
||||
Object actual;
|
||||
try {
|
||||
actual = fastMethod.invoke(mathMethod, params);
|
||||
} catch (InvocationTargetException ite) {
|
||||
actual = ite.getCause();
|
||||
}
|
||||
if (expected instanceof ArithmeticException) {
|
||||
Assert.assertEquals(MathArithmeticException.class, actual.getClass());
|
||||
} else if (!expected.equals(actual)) {
|
||||
reportFailedResults(mathMethod, params, expected, actual, entries);
|
||||
}
|
||||
} catch (IllegalArgumentException e) {
|
||||
Assert.fail(mathMethod+" "+e);
|
||||
}
|
||||
}
|
||||
|
||||
private static void setupMethodCall(Method mathMethod, Method fastMethod,
|
||||
Type[] types, Object[][] valueArrays) throws Exception {
|
||||
Object[] params = new Object[types.length];
|
||||
int entry1 = 0;
|
||||
int[] entries = new int[types.length];
|
||||
for(Object d : valueArrays[0]) {
|
||||
entry1++;
|
||||
params[0] = d;
|
||||
entries[0] = entry1;
|
||||
if (params.length > 1){
|
||||
int entry2 = 0;
|
||||
for(Object d1 : valueArrays[1]) {
|
||||
entry2++;
|
||||
params[1] = d1;
|
||||
entries[1] = entry2;
|
||||
callMethods(mathMethod, fastMethod, params, entries);
|
||||
}
|
||||
} else {
|
||||
callMethods(mathMethod, fastMethod, params, entries);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@Parameters
|
||||
public static List<Object[]> data() throws Exception {
|
||||
String singleMethod = System.getProperty("testMethod");
|
||||
List<Object[]> list = new ArrayList<>();
|
||||
for(Method mathMethod : StrictMath.class.getDeclaredMethods()) {
|
||||
method:
|
||||
if (Modifier.isPublic(mathMethod.getModifiers())){// Only test public methods
|
||||
Type []types = mathMethod.getGenericParameterTypes();
|
||||
if (types.length >=1) { // Only check methods with at least one parameter
|
||||
try {
|
||||
// Get the corresponding AccurateMath method
|
||||
Method fastMethod = AccurateMath.class.getDeclaredMethod(mathMethod.getName(), (Class[]) types);
|
||||
if (Modifier.isPublic(fastMethod.getModifiers())) { // It must be public too
|
||||
if (singleMethod != null && !fastMethod.getName().equals(singleMethod)) {
|
||||
break method;
|
||||
}
|
||||
Object [][] values = new Object[types.length][];
|
||||
int index = 0;
|
||||
for(Type t : types) {
|
||||
if (t.equals(double.class)){
|
||||
values[index]=DOUBLE_SPECIAL_VALUES;
|
||||
} else if (t.equals(float.class)) {
|
||||
values[index]=FLOAT_SPECIAL_VALUES;
|
||||
} else if (t.equals(long.class)) {
|
||||
values[index]=LONG_SPECIAL_VALUES;
|
||||
} else if (t.equals(int.class)) {
|
||||
values[index]=INT_SPECIAL_VALUES;
|
||||
} else {
|
||||
System.out.println("Cannot handle class "+t+" for "+mathMethod);
|
||||
break method;
|
||||
}
|
||||
index++;
|
||||
}
|
||||
// System.out.println(fastMethod);
|
||||
/*
|
||||
* The current implementation runs each method as a separate test.
|
||||
* Could be amended to run each value as a separate test
|
||||
*/
|
||||
list.add(new Object[]{mathMethod, fastMethod, types, values});
|
||||
// setupMethodCall(mathMethod, fastMethod, params, data);
|
||||
} else {
|
||||
System.out.println("Cannot find public AccurateMath method corresponding to: "+mathMethod);
|
||||
}
|
||||
} catch (NoSuchMethodException e) {
|
||||
System.out.println("Cannot find AccurateMath method corresponding to: "+mathMethod);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
return list;
|
||||
}
|
||||
}
|
||||
File diff suppressed because it is too large
Load Diff
Loading…
Reference in New Issue