source: src/main/java/agents/anac/y2019/harddealer/math3/util/FastMathCalc.java

Last change on this file was 204, checked in by Katsuhide Fujita, 5 years ago

Fixed errors of ANAC2019 agents
  • Property svn:executable set to *
File size: 21.2 KB
Line 
1/*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17package agents.anac.y2019.harddealer.math3.util;
18
19import java.io.PrintStream;
20
21import agents.anac.y2019.harddealer.math3.exception.DimensionMismatchException;
22
23/** Class used to compute the classical functions tables.
24 * @since 3.0
25 */
26class FastMathCalc {
27
28 /**
29 * 0x40000000 - used to split a double into two parts, both with the low order bits cleared.
30 * Equivalent to 2^30.
31 */
32 private static final long HEX_40000000 = 0x40000000L; // 1073741824L
33
34 /** Factorial table, for Taylor series expansions. 0!, 1!, 2!, ... 19! */
35 private static final double FACT[] = new double[]
36 {
37 +1.0d, // 0
38 +1.0d, // 1
39 +2.0d, // 2
40 +6.0d, // 3
41 +24.0d, // 4
42 +120.0d, // 5
43 +720.0d, // 6
44 +5040.0d, // 7
45 +40320.0d, // 8
46 +362880.0d, // 9
47 +3628800.0d, // 10
48 +39916800.0d, // 11
49 +479001600.0d, // 12
50 +6227020800.0d, // 13
51 +87178291200.0d, // 14
52 +1307674368000.0d, // 15
53 +20922789888000.0d, // 16
54 +355687428096000.0d, // 17
55 +6402373705728000.0d, // 18
56 +121645100408832000.0d, // 19
57 };
58
59 /** Coefficients for slowLog. */
60 private static final double LN_SPLIT_COEF[][] = {
61 {2.0, 0.0},
62 {0.6666666269302368, 3.9736429850260626E-8},
63 {0.3999999761581421, 2.3841857910019882E-8},
64 {0.2857142686843872, 1.7029898543501842E-8},
65 {0.2222222089767456, 1.3245471311735498E-8},
66 {0.1818181574344635, 2.4384203044354907E-8},
67 {0.1538461446762085, 9.140260083262505E-9},
68 {0.13333332538604736, 9.220590270857665E-9},
69 {0.11764700710773468, 1.2393345855018391E-8},
70 {0.10526403784751892, 8.251545029714408E-9},
71 {0.0952233225107193, 1.2675934823758863E-8},
72 {0.08713622391223907, 1.1430250008909141E-8},
73 {0.07842259109020233, 2.404307984052299E-9},
74 {0.08371849358081818, 1.176342548272881E-8},
75 {0.030589580535888672, 1.2958646899018938E-9},
76 {0.14982303977012634, 1.225743062930824E-8},
77 };
78
79 /** Table start declaration. */
80 private static final String TABLE_START_DECL = " {";
81
82 /** Table end declaration. */
83 private static final String TABLE_END_DECL = " };";
84
85 /**
86 * Private Constructor.
87 */
88 private FastMathCalc() {
89 }
90
91 /** Build the sine and cosine tables.
92 * @param SINE_TABLE_A table of the most significant part of the sines
93 * @param SINE_TABLE_B table of the least significant part of the sines
94 * @param COSINE_TABLE_A table of the most significant part of the cosines
95 * @param COSINE_TABLE_B table of the most significant part of the cosines
96 * @param SINE_TABLE_LEN length of the tables
97 * @param TANGENT_TABLE_A table of the most significant part of the tangents
98 * @param TANGENT_TABLE_B table of the most significant part of the tangents
99 */
100 @SuppressWarnings("unused")
101 private static void buildSinCosTables(double[] SINE_TABLE_A, double[] SINE_TABLE_B,
102 double[] COSINE_TABLE_A, double[] COSINE_TABLE_B,
103 int SINE_TABLE_LEN, double[] TANGENT_TABLE_A, double[] TANGENT_TABLE_B) {
104 final double result[] = new double[2];
105
106 /* Use taylor series for 0 <= x <= 6/8 */
107 for (int i = 0; i < 7; i++) {
108 double x = i / 8.0;
109
110 slowSin(x, result);
111 SINE_TABLE_A[i] = result[0];
112 SINE_TABLE_B[i] = result[1];
113
114 slowCos(x, result);
115 COSINE_TABLE_A[i] = result[0];
116 COSINE_TABLE_B[i] = result[1];
117 }
118
119 /* Use angle addition formula to complete table to 13/8, just beyond pi/2 */
120 for (int i = 7; i < SINE_TABLE_LEN; i++) {
121 double xs[] = new double[2];
122 double ys[] = new double[2];
123 double as[] = new double[2];
124 double bs[] = new double[2];
125 double temps[] = new double[2];
126
127 if ( (i & 1) == 0) {
128 // Even, use double angle
129 xs[0] = SINE_TABLE_A[i/2];
130 xs[1] = SINE_TABLE_B[i/2];
131 ys[0] = COSINE_TABLE_A[i/2];
132 ys[1] = COSINE_TABLE_B[i/2];
133
134 /* compute sine */
135 splitMult(xs, ys, result);
136 SINE_TABLE_A[i] = result[0] * 2.0;
137 SINE_TABLE_B[i] = result[1] * 2.0;
138
139 /* Compute cosine */
140 splitMult(ys, ys, as);
141 splitMult(xs, xs, temps);
142 temps[0] = -temps[0];
143 temps[1] = -temps[1];
144 splitAdd(as, temps, result);
145 COSINE_TABLE_A[i] = result[0];
146 COSINE_TABLE_B[i] = result[1];
147 } else {
148 xs[0] = SINE_TABLE_A[i/2];
149 xs[1] = SINE_TABLE_B[i/2];
150 ys[0] = COSINE_TABLE_A[i/2];
151 ys[1] = COSINE_TABLE_B[i/2];
152 as[0] = SINE_TABLE_A[i/2+1];
153 as[1] = SINE_TABLE_B[i/2+1];
154 bs[0] = COSINE_TABLE_A[i/2+1];
155 bs[1] = COSINE_TABLE_B[i/2+1];
156
157 /* compute sine */
158 splitMult(xs, bs, temps);
159 splitMult(ys, as, result);
160 splitAdd(result, temps, result);
161 SINE_TABLE_A[i] = result[0];
162 SINE_TABLE_B[i] = result[1];
163
164 /* Compute cosine */
165 splitMult(ys, bs, result);
166 splitMult(xs, as, temps);
167 temps[0] = -temps[0];
168 temps[1] = -temps[1];
169 splitAdd(result, temps, result);
170 COSINE_TABLE_A[i] = result[0];
171 COSINE_TABLE_B[i] = result[1];
172 }
173 }
174
175 /* Compute tangent = sine/cosine */
176 for (int i = 0; i < SINE_TABLE_LEN; i++) {
177 double xs[] = new double[2];
178 double ys[] = new double[2];
179 double as[] = new double[2];
180
181 as[0] = COSINE_TABLE_A[i];
182 as[1] = COSINE_TABLE_B[i];
183
184 splitReciprocal(as, ys);
185
186 xs[0] = SINE_TABLE_A[i];
187 xs[1] = SINE_TABLE_B[i];
188
189 splitMult(xs, ys, as);
190
191 TANGENT_TABLE_A[i] = as[0];
192 TANGENT_TABLE_B[i] = as[1];
193 }
194
195 }
196
197 /**
198 * For x between 0 and pi/4 compute cosine using Talor series
199 * cos(x) = 1 - x^2/2! + x^4/4! ...
200 * @param x number from which cosine is requested
201 * @param result placeholder where to put the result in extended precision
202 * (may be null)
203 * @return cos(x)
204 */
205 static double slowCos(final double x, final double result[]) {
206
207 final double xs[] = new double[2];
208 final double ys[] = new double[2];
209 final double facts[] = new double[2];
210 final double as[] = new double[2];
211 split(x, xs);
212 ys[0] = ys[1] = 0.0;
213
214 for (int i = FACT.length-1; i >= 0; i--) {
215 splitMult(xs, ys, as);
216 ys[0] = as[0]; ys[1] = as[1];
217
218 if ( (i & 1) != 0) { // skip odd entries
219 continue;
220 }
221
222 split(FACT[i], as);
223 splitReciprocal(as, facts);
224
225 if ( (i & 2) != 0 ) { // alternate terms are negative
226 facts[0] = -facts[0];
227 facts[1] = -facts[1];
228 }
229
230 splitAdd(ys, facts, as);
231 ys[0] = as[0]; ys[1] = as[1];
232 }
233
234 if (result != null) {
235 result[0] = ys[0];
236 result[1] = ys[1];
237 }
238
239 return ys[0] + ys[1];
240 }
241
242 /**
243 * For x between 0 and pi/4 compute sine using Taylor expansion:
244 * sin(x) = x - x^3/3! + x^5/5! - x^7/7! ...
245 * @param x number from which sine is requested
246 * @param result placeholder where to put the result in extended precision
247 * (may be null)
248 * @return sin(x)
249 */
250 static double slowSin(final double x, final double result[]) {
251 final double xs[] = new double[2];
252 final double ys[] = new double[2];
253 final double facts[] = new double[2];
254 final double as[] = new double[2];
255 split(x, xs);
256 ys[0] = ys[1] = 0.0;
257
258 for (int i = FACT.length-1; i >= 0; i--) {
259 splitMult(xs, ys, as);
260 ys[0] = as[0]; ys[1] = as[1];
261
262 if ( (i & 1) == 0) { // Ignore even numbers
263 continue;
264 }
265
266 split(FACT[i], as);
267 splitReciprocal(as, facts);
268
269 if ( (i & 2) != 0 ) { // alternate terms are negative
270 facts[0] = -facts[0];
271 facts[1] = -facts[1];
272 }
273
274 splitAdd(ys, facts, as);
275 ys[0] = as[0]; ys[1] = as[1];
276 }
277
278 if (result != null) {
279 result[0] = ys[0];
280 result[1] = ys[1];
281 }
282
283 return ys[0] + ys[1];
284 }
285
286
287 /**
288 * For x between 0 and 1, returns exp(x), uses extended precision
289 * @param x argument of exponential
290 * @param result placeholder where to place exp(x) split in two terms
291 * for extra precision (i.e. exp(x) = result[0] + result[1]
292 * @return exp(x)
293 */
294 static double slowexp(final double x, final double result[]) {
295 final double xs[] = new double[2];
296 final double ys[] = new double[2];
297 final double facts[] = new double[2];
298 final double as[] = new double[2];
299 split(x, xs);
300 ys[0] = ys[1] = 0.0;
301
302 for (int i = FACT.length-1; i >= 0; i--) {
303 splitMult(xs, ys, as);
304 ys[0] = as[0];
305 ys[1] = as[1];
306
307 split(FACT[i], as);
308 splitReciprocal(as, facts);
309
310 splitAdd(ys, facts, as);
311 ys[0] = as[0];
312 ys[1] = as[1];
313 }
314
315 if (result != null) {
316 result[0] = ys[0];
317 result[1] = ys[1];
318 }
319
320 return ys[0] + ys[1];
321 }
322
323 /** Compute split[0], split[1] such that their sum is equal to d,
324 * and split[0] has its 30 least significant bits as zero.
325 * @param d number to split
326 * @param split placeholder where to place the result
327 */
328 private static void split(final double d, final double split[]) {
329 if (d < 8e298 && d > -8e298) {
330 final double a = d * HEX_40000000;
331 split[0] = (d + a) - a;
332 split[1] = d - split[0];
333 } else {
334 final double a = d * 9.31322574615478515625E-10;
335 split[0] = (d + a - d) * HEX_40000000;
336 split[1] = d - split[0];
337 }
338 }
339
340 /** Recompute a split.
341 * @param a input/out array containing the split, changed
342 * on output
343 */
344 private static void resplit(final double a[]) {
345 final double c = a[0] + a[1];
346 final double d = -(c - a[0] - a[1]);
347
348 if (c < 8e298 && c > -8e298) { // MAGIC NUMBER
349 double z = c * HEX_40000000;
350 a[0] = (c + z) - z;
351 a[1] = c - a[0] + d;
352 } else {
353 double z = c * 9.31322574615478515625E-10;
354 a[0] = (c + z - c) * HEX_40000000;
355 a[1] = c - a[0] + d;
356 }
357 }
358
359 /** Multiply two numbers in split form.
360 * @param a first term of multiplication
361 * @param b second term of multiplication
362 * @param ans placeholder where to put the result
363 */
364 private static void splitMult(double a[], double b[], double ans[]) {
365 ans[0] = a[0] * b[0];
366 ans[1] = a[0] * b[1] + a[1] * b[0] + a[1] * b[1];
367
368 /* Resplit */
369 resplit(ans);
370 }
371
372 /** Add two numbers in split form.
373 * @param a first term of addition
374 * @param b second term of addition
375 * @param ans placeholder where to put the result
376 */
377 private static void splitAdd(final double a[], final double b[], final double ans[]) {
378 ans[0] = a[0] + b[0];
379 ans[1] = a[1] + b[1];
380
381 resplit(ans);
382 }
383
384 /** Compute the reciprocal of in. Use the following algorithm.
385 * in = c + d.
386 * want to find x + y such that x+y = 1/(c+d) and x is much
387 * larger than y and x has several zero bits on the right.
388 *
389 * Set b = 1/(2^22), a = 1 - b. Thus (a+b) = 1.
390 * Use following identity to compute (a+b)/(c+d)
391 *
392 * (a+b)/(c+d) = a/c + (bc - ad) / (c^2 + cd)
393 * set x = a/c and y = (bc - ad) / (c^2 + cd)
394 * This will be close to the right answer, but there will be
395 * some rounding in the calculation of X. So by carefully
396 * computing 1 - (c+d)(x+y) we can compute an error and
397 * add that back in. This is done carefully so that terms
398 * of similar size are subtracted first.
399 * @param in initial number, in split form
400 * @param result placeholder where to put the result
401 */
402 static void splitReciprocal(final double in[], final double result[]) {
403 final double b = 1.0/4194304.0;
404 final double a = 1.0 - b;
405
406 if (in[0] == 0.0) {
407 in[0] = in[1];
408 in[1] = 0.0;
409 }
410
411 result[0] = a / in[0];
412 result[1] = (b*in[0]-a*in[1]) / (in[0]*in[0] + in[0]*in[1]);
413
414 if (result[1] != result[1]) { // can happen if result[1] is NAN
415 result[1] = 0.0;
416 }
417
418 /* Resplit */
419 resplit(result);
420
421 for (int i = 0; i < 2; i++) {
422 /* this may be overkill, probably once is enough */
423 double err = 1.0 - result[0] * in[0] - result[0] * in[1] -
424 result[1] * in[0] - result[1] * in[1];
425 /*err = 1.0 - err; */
426 err *= result[0] + result[1];
427 /*printf("err = %16e\n", err); */
428 result[1] += err;
429 }
430 }
431
432 /** Compute (a[0] + a[1]) * (b[0] + b[1]) in extended precision.
433 * @param a first term of the multiplication
434 * @param b second term of the multiplication
435 * @param result placeholder where to put the result
436 */
437 private static void quadMult(final double a[], final double b[], final double result[]) {
438 final double xs[] = new double[2];
439 final double ys[] = new double[2];
440 final double zs[] = new double[2];
441
442 /* a[0] * b[0] */
443 split(a[0], xs);
444 split(b[0], ys);
445 splitMult(xs, ys, zs);
446
447 result[0] = zs[0];
448 result[1] = zs[1];
449
450 /* a[0] * b[1] */
451 split(b[1], ys);
452 splitMult(xs, ys, zs);
453
454 double tmp = result[0] + zs[0];
455 result[1] -= tmp - result[0] - zs[0];
456 result[0] = tmp;
457 tmp = result[0] + zs[1];
458 result[1] -= tmp - result[0] - zs[1];
459 result[0] = tmp;
460
461 /* a[1] * b[0] */
462 split(a[1], xs);
463 split(b[0], ys);
464 splitMult(xs, ys, zs);
465
466 tmp = result[0] + zs[0];
467 result[1] -= tmp - result[0] - zs[0];
468 result[0] = tmp;
469 tmp = result[0] + zs[1];
470 result[1] -= tmp - result[0] - zs[1];
471 result[0] = tmp;
472
473 /* a[1] * b[0] */
474 split(a[1], xs);
475 split(b[1], ys);
476 splitMult(xs, ys, zs);
477
478 tmp = result[0] + zs[0];
479 result[1] -= tmp - result[0] - zs[0];
480 result[0] = tmp;
481 tmp = result[0] + zs[1];
482 result[1] -= tmp - result[0] - zs[1];
483 result[0] = tmp;
484 }
485
486 /** Compute exp(p) for a integer p in extended precision.
487 * @param p integer whose exponential is requested
488 * @param result placeholder where to put the result in extended precision
489 * @return exp(p) in standard precision (equal to result[0] + result[1])
490 */
491 static double expint(int p, final double result[]) {
492 //double x = M_E;
493 final double xs[] = new double[2];
494 final double as[] = new double[2];
495 final double ys[] = new double[2];
496 //split(x, xs);
497 //xs[1] = (double)(2.7182818284590452353602874713526625L - xs[0]);
498 //xs[0] = 2.71827697753906250000;
499 //xs[1] = 4.85091998273542816811e-06;
500 //xs[0] = Double.longBitsToDouble(0x4005bf0800000000L);
501 //xs[1] = Double.longBitsToDouble(0x3ed458a2bb4a9b00L);
502
503 /* E */
504 xs[0] = 2.718281828459045;
505 xs[1] = 1.4456468917292502E-16;
506
507 split(1.0, ys);
508
509 while (p > 0) {
510 if ((p & 1) != 0) {
511 quadMult(ys, xs, as);
512 ys[0] = as[0]; ys[1] = as[1];
513 }
514
515 quadMult(xs, xs, as);
516 xs[0] = as[0]; xs[1] = as[1];
517
518 p >>= 1;
519 }
520
521 if (result != null) {
522 result[0] = ys[0];
523 result[1] = ys[1];
524
525 resplit(result);
526 }
527
528 return ys[0] + ys[1];
529 }
530 /** xi in the range of [1, 2].
531 * 3 5 7
532 * x+1 / x x x \
533 * ln ----- = 2 * | x + ---- + ---- + ---- + ... |
534 * 1-x \ 3 5 7 /
535 *
536 * So, compute a Remez approximation of the following function
537 *
538 * ln ((sqrt(x)+1)/(1-sqrt(x))) / x
539 *
540 * This will be an even function with only positive coefficents.
541 * x is in the range [0 - 1/3].
542 *
543 * Transform xi for input to the above function by setting
544 * x = (xi-1)/(xi+1). Input to the polynomial is x^2, then
545 * the result is multiplied by x.
546 * @param xi number from which log is requested
547 * @return log(xi)
548 */
549 static double[] slowLog(double xi) {
550 double x[] = new double[2];
551 double x2[] = new double[2];
552 double y[] = new double[2];
553 double a[] = new double[2];
554
555 split(xi, x);
556
557 /* Set X = (x-1)/(x+1) */
558 x[0] += 1.0;
559 resplit(x);
560 splitReciprocal(x, a);
561 x[0] -= 2.0;
562 resplit(x);
563 splitMult(x, a, y);
564 x[0] = y[0];
565 x[1] = y[1];
566
567 /* Square X -> X2*/
568 splitMult(x, x, x2);
569
570
571 //x[0] -= 1.0;
572 //resplit(x);
573
574 y[0] = LN_SPLIT_COEF[LN_SPLIT_COEF.length-1][0];
575 y[1] = LN_SPLIT_COEF[LN_SPLIT_COEF.length-1][1];
576
577 for (int i = LN_SPLIT_COEF.length-2; i >= 0; i--) {
578 splitMult(y, x2, a);
579 y[0] = a[0];
580 y[1] = a[1];
581 splitAdd(y, LN_SPLIT_COEF[i], a);
582 y[0] = a[0];
583 y[1] = a[1];
584 }
585
586 splitMult(y, x, a);
587 y[0] = a[0];
588 y[1] = a[1];
589
590 return y;
591 }
592
593
594 /**
595 * Print an array.
596 * @param out text output stream where output should be printed
597 * @param name array name
598 * @param expectedLen expected length of the array
599 * @param array2d array data
600 */
601 static void printarray(PrintStream out, String name, int expectedLen, double[][] array2d) {
602 out.println(name);
603 checkLen(expectedLen, array2d.length);
604 out.println(TABLE_START_DECL + " ");
605 int i = 0;
606 for(double[] array : array2d) { // "double array[]" causes PMD parsing error
607 out.print(" {");
608 for(double d : array) { // assume inner array has very few entries
609 out.printf("%-25.25s", format(d)); // multiple entries per line
610 }
611 out.println("}, // " + i++);
612 }
613 out.println(TABLE_END_DECL);
614 }
615
616 /**
617 * Print an array.
618 * @param out text output stream where output should be printed
619 * @param name array name
620 * @param expectedLen expected length of the array
621 * @param array array data
622 */
623 static void printarray(PrintStream out, String name, int expectedLen, double[] array) {
624 out.println(name + "=");
625 checkLen(expectedLen, array.length);
626 out.println(TABLE_START_DECL);
627 for(double d : array){
628 out.printf(" %s%n", format(d)); // one entry per line
629 }
630 out.println(TABLE_END_DECL);
631 }
632
633 /** Format a double.
634 * @param d double number to format
635 * @return formatted number
636 */
637 static String format(double d) {
638 if (d != d) {
639 return "Double.NaN,";
640 } else {
641 return ((d >= 0) ? "+" : "") + Double.toString(d) + "d,";
642 }
643 }
644
645 /**
646 * Check two lengths are equal.
647 * @param expectedLen expected length
648 * @param actual actual length
649 * @exception DimensionMismatchException if the two lengths are not equal
650 */
651 private static void checkLen(int expectedLen, int actual)
652 throws DimensionMismatchException {
653 if (expectedLen != actual) {
654 throw new DimensionMismatchException(actual, expectedLen);
655 }
656 }
657
658}
Note: See TracBrowser for help on using the repository browser.