1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.apache.commons.math3.analysis.differentiation;
18
19 import java.util.ArrayList;
20 import java.util.Arrays;
21 import java.util.List;
22 import java.util.concurrent.atomic.AtomicReference;
23
24 import org.apache.commons.math3.exception.DimensionMismatchException;
25 import org.apache.commons.math3.exception.MathArithmeticException;
26 import org.apache.commons.math3.exception.MathInternalError;
27 import org.apache.commons.math3.exception.NotPositiveException;
28 import org.apache.commons.math3.exception.NumberIsTooLargeException;
29 import org.apache.commons.math3.util.CombinatoricsUtils;
30 import org.apache.commons.math3.util.FastMath;
31 import org.apache.commons.math3.util.MathArrays;
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127 public class DSCompiler {
128
129
130 private static AtomicReference<DSCompiler[][]> compilers =
131 new AtomicReference<DSCompiler[][]>(null);
132
133
134 private final int parameters;
135
136
137 private final int order;
138
139
140 private final int[][] sizes;
141
142
143 private final int[][] derivativesIndirection;
144
145
146 private final int[] lowerIndirection;
147
148
149 private final int[][][] multIndirection;
150
151
152 private final int[][][] compIndirection;
153
154
155
156
157
158
159
160
161 private DSCompiler(final int parameters, final int order,
162 final DSCompiler valueCompiler, final DSCompiler derivativeCompiler)
163 throws NumberIsTooLargeException {
164
165 this.parameters = parameters;
166 this.order = order;
167 this.sizes = compileSizes(parameters, order, valueCompiler);
168 this.derivativesIndirection =
169 compileDerivativesIndirection(parameters, order,
170 valueCompiler, derivativeCompiler);
171 this.lowerIndirection =
172 compileLowerIndirection(parameters, order,
173 valueCompiler, derivativeCompiler);
174 this.multIndirection =
175 compileMultiplicationIndirection(parameters, order,
176 valueCompiler, derivativeCompiler, lowerIndirection);
177 this.compIndirection =
178 compileCompositionIndirection(parameters, order,
179 valueCompiler, derivativeCompiler,
180 sizes, derivativesIndirection);
181
182 }
183
184
185
186
187
188
189
190 public static DSCompiler getCompiler(int parameters, int order)
191 throws NumberIsTooLargeException {
192
193
194 final DSCompiler[][] cache = compilers.get();
195 if (cache != null && cache.length > parameters &&
196 cache[parameters].length > order && cache[parameters][order] != null) {
197
198 return cache[parameters][order];
199 }
200
201
202 final int maxParameters = FastMath.max(parameters, cache == null ? 0 : cache.length);
203 final int maxOrder = FastMath.max(order, cache == null ? 0 : cache[0].length);
204 final DSCompiler[][] newCache = new DSCompiler[maxParameters + 1][maxOrder + 1];
205
206 if (cache != null) {
207
208 for (int i = 0; i < cache.length; ++i) {
209 System.arraycopy(cache[i], 0, newCache[i], 0, cache[i].length);
210 }
211 }
212
213
214 for (int diag = 0; diag <= parameters + order; ++diag) {
215 for (int o = FastMath.max(0, diag - parameters); o <= FastMath.min(order, diag); ++o) {
216 final int p = diag - o;
217 if (newCache[p][o] == null) {
218 final DSCompiler valueCompiler = (p == 0) ? null : newCache[p - 1][o];
219 final DSCompiler derivativeCompiler = (o == 0) ? null : newCache[p][o - 1];
220 newCache[p][o] = new DSCompiler(p, o, valueCompiler, derivativeCompiler);
221 }
222 }
223 }
224
225
226 compilers.compareAndSet(cache, newCache);
227
228 return newCache[parameters][order];
229
230 }
231
232
233
234
235
236
237
238 private static int[][] compileSizes(final int parameters, final int order,
239 final DSCompiler valueCompiler) {
240
241 final int[][] sizes = new int[parameters + 1][order + 1];
242 if (parameters == 0) {
243 Arrays.fill(sizes[0], 1);
244 } else {
245 System.arraycopy(valueCompiler.sizes, 0, sizes, 0, parameters);
246 sizes[parameters][0] = 1;
247 for (int i = 0; i < order; ++i) {
248 sizes[parameters][i + 1] = sizes[parameters][i] + sizes[parameters - 1][i + 1];
249 }
250 }
251
252 return sizes;
253
254 }
255
256
257
258
259
260
261
262
263 private static int[][] compileDerivativesIndirection(final int parameters, final int order,
264 final DSCompiler valueCompiler,
265 final DSCompiler derivativeCompiler) {
266
267 if (parameters == 0 || order == 0) {
268 return new int[1][parameters];
269 }
270
271 final int vSize = valueCompiler.derivativesIndirection.length;
272 final int dSize = derivativeCompiler.derivativesIndirection.length;
273 final int[][] derivativesIndirection = new int[vSize + dSize][parameters];
274
275
276 for (int i = 0; i < vSize; ++i) {
277
278 System.arraycopy(valueCompiler.derivativesIndirection[i], 0,
279 derivativesIndirection[i], 0,
280 parameters - 1);
281 }
282
283
284 for (int i = 0; i < dSize; ++i) {
285
286
287 System.arraycopy(derivativeCompiler.derivativesIndirection[i], 0,
288 derivativesIndirection[vSize + i], 0,
289 parameters);
290
291
292 derivativesIndirection[vSize + i][parameters - 1]++;
293
294 }
295
296 return derivativesIndirection;
297
298 }
299
300
301
302
303
304
305
306
307
308
309
310
311 private static int[] compileLowerIndirection(final int parameters, final int order,
312 final DSCompiler valueCompiler,
313 final DSCompiler derivativeCompiler) {
314
315 if (parameters == 0 || order <= 1) {
316 return new int[] { 0 };
317 }
318
319
320 final int vSize = valueCompiler.lowerIndirection.length;
321 final int dSize = derivativeCompiler.lowerIndirection.length;
322 final int[] lowerIndirection = new int[vSize + dSize];
323 System.arraycopy(valueCompiler.lowerIndirection, 0, lowerIndirection, 0, vSize);
324 for (int i = 0; i < dSize; ++i) {
325 lowerIndirection[vSize + i] = valueCompiler.getSize() + derivativeCompiler.lowerIndirection[i];
326 }
327
328 return lowerIndirection;
329
330 }
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345 private static int[][][] compileMultiplicationIndirection(final int parameters, final int order,
346 final DSCompiler valueCompiler,
347 final DSCompiler derivativeCompiler,
348 final int[] lowerIndirection) {
349
350 if ((parameters == 0) || (order == 0)) {
351 return new int[][][] { { { 1, 0, 0 } } };
352 }
353
354
355 final int vSize = valueCompiler.multIndirection.length;
356 final int dSize = derivativeCompiler.multIndirection.length;
357 final int[][][] multIndirection = new int[vSize + dSize][][];
358
359 System.arraycopy(valueCompiler.multIndirection, 0, multIndirection, 0, vSize);
360
361 for (int i = 0; i < dSize; ++i) {
362 final int[][] dRow = derivativeCompiler.multIndirection[i];
363 List<int[]> row = new ArrayList<int[]>(dRow.length * 2);
364 for (int j = 0; j < dRow.length; ++j) {
365 row.add(new int[] { dRow[j][0], lowerIndirection[dRow[j][1]], vSize + dRow[j][2] });
366 row.add(new int[] { dRow[j][0], vSize + dRow[j][1], lowerIndirection[dRow[j][2]] });
367 }
368
369
370 final List<int[]> combined = new ArrayList<int[]>(row.size());
371 for (int j = 0; j < row.size(); ++j) {
372 final int[] termJ = row.get(j);
373 if (termJ[0] > 0) {
374 for (int k = j + 1; k < row.size(); ++k) {
375 final int[] termK = row.get(k);
376 if (termJ[1] == termK[1] && termJ[2] == termK[2]) {
377
378 termJ[0] += termK[0];
379
380 termK[0] = 0;
381 }
382 }
383 combined.add(termJ);
384 }
385 }
386
387 multIndirection[vSize + i] = combined.toArray(new int[combined.size()][]);
388
389 }
390
391 return multIndirection;
392
393 }
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410 private static int[][][] compileCompositionIndirection(final int parameters, final int order,
411 final DSCompiler valueCompiler,
412 final DSCompiler derivativeCompiler,
413 final int[][] sizes,
414 final int[][] derivativesIndirection)
415 throws NumberIsTooLargeException {
416
417 if ((parameters == 0) || (order == 0)) {
418 return new int[][][] { { { 1, 0 } } };
419 }
420
421 final int vSize = valueCompiler.compIndirection.length;
422 final int dSize = derivativeCompiler.compIndirection.length;
423 final int[][][] compIndirection = new int[vSize + dSize][][];
424
425
426 System.arraycopy(valueCompiler.compIndirection, 0, compIndirection, 0, vSize);
427
428
429
430
431
432 for (int i = 0; i < dSize; ++i) {
433 List<int[]> row = new ArrayList<int[]>();
434 for (int[] term : derivativeCompiler.compIndirection[i]) {
435
436
437
438
439 int[] derivedTermF = new int[term.length + 1];
440 derivedTermF[0] = term[0];
441 derivedTermF[1] = term[1] + 1;
442 int[] orders = new int[parameters];
443 orders[parameters - 1] = 1;
444 derivedTermF[term.length] = getPartialDerivativeIndex(parameters, order, sizes, orders);
445 for (int j = 2; j < term.length; ++j) {
446
447
448 derivedTermF[j] = convertIndex(term[j], parameters,
449 derivativeCompiler.derivativesIndirection,
450 parameters, order, sizes);
451 }
452 Arrays.sort(derivedTermF, 2, derivedTermF.length);
453 row.add(derivedTermF);
454
455
456 for (int l = 2; l < term.length; ++l) {
457 int[] derivedTermG = new int[term.length];
458 derivedTermG[0] = term[0];
459 derivedTermG[1] = term[1];
460 for (int j = 2; j < term.length; ++j) {
461
462
463 derivedTermG[j] = convertIndex(term[j], parameters,
464 derivativeCompiler.derivativesIndirection,
465 parameters, order, sizes);
466 if (j == l) {
467
468 System.arraycopy(derivativesIndirection[derivedTermG[j]], 0, orders, 0, parameters);
469 orders[parameters - 1]++;
470 derivedTermG[j] = getPartialDerivativeIndex(parameters, order, sizes, orders);
471 }
472 }
473 Arrays.sort(derivedTermG, 2, derivedTermG.length);
474 row.add(derivedTermG);
475 }
476
477 }
478
479
480 final List<int[]> combined = new ArrayList<int[]>(row.size());
481 for (int j = 0; j < row.size(); ++j) {
482 final int[] termJ = row.get(j);
483 if (termJ[0] > 0) {
484 for (int k = j + 1; k < row.size(); ++k) {
485 final int[] termK = row.get(k);
486 boolean equals = termJ.length == termK.length;
487 for (int l = 1; equals && l < termJ.length; ++l) {
488 equals &= termJ[l] == termK[l];
489 }
490 if (equals) {
491
492 termJ[0] += termK[0];
493
494 termK[0] = 0;
495 }
496 }
497 combined.add(termJ);
498 }
499 }
500
501 compIndirection[vSize + i] = combined.toArray(new int[combined.size()][]);
502
503 }
504
505 return compIndirection;
506
507 }
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541 public int getPartialDerivativeIndex(final int ... orders)
542 throws DimensionMismatchException, NumberIsTooLargeException {
543
544
545 if (orders.length != getFreeParameters()) {
546 throw new DimensionMismatchException(orders.length, getFreeParameters());
547 }
548
549 return getPartialDerivativeIndex(parameters, order, sizes, orders);
550
551 }
552
553
554
555
556
557
558
559
560
561
562
563 private static int getPartialDerivativeIndex(final int parameters, final int order,
564 final int[][] sizes, final int ... orders)
565 throws NumberIsTooLargeException {
566
567
568
569 int index = 0;
570 int m = order;
571 int ordersSum = 0;
572 for (int i = parameters - 1; i >= 0; --i) {
573
574
575 int derivativeOrder = orders[i];
576
577
578 ordersSum += derivativeOrder;
579 if (ordersSum > order) {
580 throw new NumberIsTooLargeException(ordersSum, order, true);
581 }
582
583 while (derivativeOrder-- > 0) {
584
585
586
587 index += sizes[i][m--];
588 }
589
590 }
591
592 return index;
593
594 }
595
596
597
598
599
600
601
602
603
604
605
606
607
608 private static int convertIndex(final int index,
609 final int srcP, final int[][] srcDerivativesIndirection,
610 final int destP, final int destO, final int[][] destSizes)
611 throws NumberIsTooLargeException {
612 int[] orders = new int[destP];
613 System.arraycopy(srcDerivativesIndirection[index], 0, orders, 0, FastMath.min(srcP, destP));
614 return getPartialDerivativeIndex(destP, destO, destSizes, orders);
615 }
616
617
618
619
620
621
622
623
624
625 public int[] getPartialDerivativeOrders(final int index) {
626 return derivativesIndirection[index];
627 }
628
629
630
631
632 public int getFreeParameters() {
633 return parameters;
634 }
635
636
637
638
639 public int getOrder() {
640 return order;
641 }
642
643
644
645
646
647
648
649
650 public int getSize() {
651 return sizes[parameters][order];
652 }
653
654
655
656
657
658
659
660
661
662
663
664
665
666 public void linearCombination(final double a1, final double[] c1, final int offset1,
667 final double a2, final double[] c2, final int offset2,
668 final double[] result, final int resultOffset) {
669 for (int i = 0; i < getSize(); ++i) {
670 result[resultOffset + i] =
671 MathArrays.linearCombination(a1, c1[offset1 + i], a2, c2[offset2 + i]);
672 }
673 }
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690 public void linearCombination(final double a1, final double[] c1, final int offset1,
691 final double a2, final double[] c2, final int offset2,
692 final double a3, final double[] c3, final int offset3,
693 final double[] result, final int resultOffset) {
694 for (int i = 0; i < getSize(); ++i) {
695 result[resultOffset + i] =
696 MathArrays.linearCombination(a1, c1[offset1 + i],
697 a2, c2[offset2 + i],
698 a3, c3[offset3 + i]);
699 }
700 }
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720 public void linearCombination(final double a1, final double[] c1, final int offset1,
721 final double a2, final double[] c2, final int offset2,
722 final double a3, final double[] c3, final int offset3,
723 final double a4, final double[] c4, final int offset4,
724 final double[] result, final int resultOffset) {
725 for (int i = 0; i < getSize(); ++i) {
726 result[resultOffset + i] =
727 MathArrays.linearCombination(a1, c1[offset1 + i],
728 a2, c2[offset2 + i],
729 a3, c3[offset3 + i],
730 a4, c4[offset4 + i]);
731 }
732 }
733
734
735
736
737
738
739
740
741
742
743 public void add(final double[] lhs, final int lhsOffset,
744 final double[] rhs, final int rhsOffset,
745 final double[] result, final int resultOffset) {
746 for (int i = 0; i < getSize(); ++i) {
747 result[resultOffset + i] = lhs[lhsOffset + i] + rhs[rhsOffset + i];
748 }
749 }
750
751
752
753
754
755
756
757
758
759 public void subtract(final double[] lhs, final int lhsOffset,
760 final double[] rhs, final int rhsOffset,
761 final double[] result, final int resultOffset) {
762 for (int i = 0; i < getSize(); ++i) {
763 result[resultOffset + i] = lhs[lhsOffset + i] - rhs[rhsOffset + i];
764 }
765 }
766
767
768
769
770
771
772
773
774
775
776
777 public void multiply(final double[] lhs, final int lhsOffset,
778 final double[] rhs, final int rhsOffset,
779 final double[] result, final int resultOffset) {
780 for (int i = 0; i < multIndirection.length; ++i) {
781 final int[][] mappingI = multIndirection[i];
782 double r = 0;
783 for (int j = 0; j < mappingI.length; ++j) {
784 r += mappingI[j][0] *
785 lhs[lhsOffset + mappingI[j][1]] *
786 rhs[rhsOffset + mappingI[j][2]];
787 }
788 result[resultOffset + i] = r;
789 }
790 }
791
792
793
794
795
796
797
798
799
800
801
802 public void divide(final double[] lhs, final int lhsOffset,
803 final double[] rhs, final int rhsOffset,
804 final double[] result, final int resultOffset) {
805 final double[] reciprocal = new double[getSize()];
806 pow(rhs, lhsOffset, -1, reciprocal, 0);
807 multiply(lhs, lhsOffset, reciprocal, 0, result, resultOffset);
808 }
809
810
811
812
813
814
815
816
817
818
819 public void remainder(final double[] lhs, final int lhsOffset,
820 final double[] rhs, final int rhsOffset,
821 final double[] result, final int resultOffset) {
822
823
824 final double rem = FastMath.IEEEremainder(lhs[lhsOffset], rhs[rhsOffset]);
825 final double k = FastMath.rint((lhs[lhsOffset] - rem) / rhs[rhsOffset]);
826
827
828 result[resultOffset] = rem;
829
830
831 for (int i = 1; i < getSize(); ++i) {
832 result[resultOffset + i] = lhs[lhsOffset + i] - k * rhs[rhsOffset + i];
833 }
834
835 }
836
837
838
839
840
841
842
843
844
845
846
847 public void pow(final double a,
848 final double[] operand, final int operandOffset,
849 final double[] result, final int resultOffset) {
850
851
852
853 final double[] function = new double[1 + order];
854 if (a == 0) {
855 if (operand[operandOffset] == 0) {
856 function[0] = 1;
857 double infinity = Double.POSITIVE_INFINITY;
858 for (int i = 1; i < function.length; ++i) {
859 infinity = -infinity;
860 function[i] = infinity;
861 }
862 } else if (operand[operandOffset] < 0) {
863 Arrays.fill(function, Double.NaN);
864 }
865 } else {
866 function[0] = FastMath.pow(a, operand[operandOffset]);
867 final double lnA = FastMath.log(a);
868 for (int i = 1; i < function.length; ++i) {
869 function[i] = lnA * function[i - 1];
870 }
871 }
872
873
874
875 compose(operand, operandOffset, function, result, resultOffset);
876
877 }
878
879
880
881
882
883
884
885
886
887
888 public void pow(final double[] operand, final int operandOffset, final double p,
889 final double[] result, final int resultOffset) {
890
891
892
893 double[] function = new double[1 + order];
894 double xk = FastMath.pow(operand[operandOffset], p - order);
895 for (int i = order; i > 0; --i) {
896 function[i] = xk;
897 xk *= operand[operandOffset];
898 }
899 function[0] = xk;
900 double coefficient = p;
901 for (int i = 1; i <= order; ++i) {
902 function[i] *= coefficient;
903 coefficient *= p - i;
904 }
905
906
907 compose(operand, operandOffset, function, result, resultOffset);
908
909 }
910
911
912
913
914
915
916
917
918
919
920 public void pow(final double[] operand, final int operandOffset, final int n,
921 final double[] result, final int resultOffset) {
922
923 if (n == 0) {
924
925 result[resultOffset] = 1.0;
926 Arrays.fill(result, resultOffset + 1, resultOffset + getSize(), 0);
927 return;
928 }
929
930
931
932 double[] function = new double[1 + order];
933
934 if (n > 0) {
935
936 final int maxOrder = FastMath.min(order, n);
937 double xk = FastMath.pow(operand[operandOffset], n - maxOrder);
938 for (int i = maxOrder; i > 0; --i) {
939 function[i] = xk;
940 xk *= operand[operandOffset];
941 }
942 function[0] = xk;
943 } else {
944
945 final double inv = 1.0 / operand[operandOffset];
946 double xk = FastMath.pow(inv, -n);
947 for (int i = 0; i <= order; ++i) {
948 function[i] = xk;
949 xk *= inv;
950 }
951 }
952
953 double coefficient = n;
954 for (int i = 1; i <= order; ++i) {
955 function[i] *= coefficient;
956 coefficient *= n - i;
957 }
958
959
960 compose(operand, operandOffset, function, result, resultOffset);
961
962 }
963
964
965
966
967
968
969
970
971
972
973
974 public void pow(final double[] x, final int xOffset,
975 final double[] y, final int yOffset,
976 final double[] result, final int resultOffset) {
977 final double[] logX = new double[getSize()];
978 log(x, xOffset, logX, 0);
979 final double[] yLogX = new double[getSize()];
980 multiply(logX, 0, y, yOffset, yLogX, 0);
981 exp(yLogX, 0, result, resultOffset);
982 }
983
984
985
986
987
988
989
990
991
992
993 public void rootN(final double[] operand, final int operandOffset, final int n,
994 final double[] result, final int resultOffset) {
995
996
997
998 double[] function = new double[1 + order];
999 double xk;
1000 if (n == 2) {
1001 function[0] = FastMath.sqrt(operand[operandOffset]);
1002 xk = 0.5 / function[0];
1003 } else if (n == 3) {
1004 function[0] = FastMath.cbrt(operand[operandOffset]);
1005 xk = 1.0 / (3.0 * function[0] * function[0]);
1006 } else {
1007 function[0] = FastMath.pow(operand[operandOffset], 1.0 / n);
1008 xk = 1.0 / (n * FastMath.pow(function[0], n - 1));
1009 }
1010 final double nReciprocal = 1.0 / n;
1011 final double xReciprocal = 1.0 / operand[operandOffset];
1012 for (int i = 1; i <= order; ++i) {
1013 function[i] = xk;
1014 xk *= xReciprocal * (nReciprocal - i);
1015 }
1016
1017
1018 compose(operand, operandOffset, function, result, resultOffset);
1019
1020 }
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030 public void exp(final double[] operand, final int operandOffset,
1031 final double[] result, final int resultOffset) {
1032
1033
1034 double[] function = new double[1 + order];
1035 Arrays.fill(function, FastMath.exp(operand[operandOffset]));
1036
1037
1038 compose(operand, operandOffset, function, result, resultOffset);
1039
1040 }
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050 public void expm1(final double[] operand, final int operandOffset,
1051 final double[] result, final int resultOffset) {
1052
1053
1054 double[] function = new double[1 + order];
1055 function[0] = FastMath.expm1(operand[operandOffset]);
1056 Arrays.fill(function, 1, 1 + order, FastMath.exp(operand[operandOffset]));
1057
1058
1059 compose(operand, operandOffset, function, result, resultOffset);
1060
1061 }
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071 public void log(final double[] operand, final int operandOffset,
1072 final double[] result, final int resultOffset) {
1073
1074
1075 double[] function = new double[1 + order];
1076 function[0] = FastMath.log(operand[operandOffset]);
1077 if (order > 0) {
1078 double inv = 1.0 / operand[operandOffset];
1079 double xk = inv;
1080 for (int i = 1; i <= order; ++i) {
1081 function[i] = xk;
1082 xk *= -i * inv;
1083 }
1084 }
1085
1086
1087 compose(operand, operandOffset, function, result, resultOffset);
1088
1089 }
1090
1091
1092
1093
1094
1095
1096
1097
1098 public void log1p(final double[] operand, final int operandOffset,
1099 final double[] result, final int resultOffset) {
1100
1101
1102 double[] function = new double[1 + order];
1103 function[0] = FastMath.log1p(operand[operandOffset]);
1104 if (order > 0) {
1105 double inv = 1.0 / (1.0 + operand[operandOffset]);
1106 double xk = inv;
1107 for (int i = 1; i <= order; ++i) {
1108 function[i] = xk;
1109 xk *= -i * inv;
1110 }
1111 }
1112
1113
1114 compose(operand, operandOffset, function, result, resultOffset);
1115
1116 }
1117
1118
1119
1120
1121
1122
1123
1124
1125 public void log10(final double[] operand, final int operandOffset,
1126 final double[] result, final int resultOffset) {
1127
1128
1129 double[] function = new double[1 + order];
1130 function[0] = FastMath.log10(operand[operandOffset]);
1131 if (order > 0) {
1132 double inv = 1.0 / operand[operandOffset];
1133 double xk = inv / FastMath.log(10.0);
1134 for (int i = 1; i <= order; ++i) {
1135 function[i] = xk;
1136 xk *= -i * inv;
1137 }
1138 }
1139
1140
1141 compose(operand, operandOffset, function, result, resultOffset);
1142
1143 }
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153 public void cos(final double[] operand, final int operandOffset,
1154 final double[] result, final int resultOffset) {
1155
1156
1157 double[] function = new double[1 + order];
1158 function[0] = FastMath.cos(operand[operandOffset]);
1159 if (order > 0) {
1160 function[1] = -FastMath.sin(operand[operandOffset]);
1161 for (int i = 2; i <= order; ++i) {
1162 function[i] = -function[i - 2];
1163 }
1164 }
1165
1166
1167 compose(operand, operandOffset, function, result, resultOffset);
1168
1169 }
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179 public void sin(final double[] operand, final int operandOffset,
1180 final double[] result, final int resultOffset) {
1181
1182
1183 double[] function = new double[1 + order];
1184 function[0] = FastMath.sin(operand[operandOffset]);
1185 if (order > 0) {
1186 function[1] = FastMath.cos(operand[operandOffset]);
1187 for (int i = 2; i <= order; ++i) {
1188 function[i] = -function[i - 2];
1189 }
1190 }
1191
1192
1193 compose(operand, operandOffset, function, result, resultOffset);
1194
1195 }
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205 public void tan(final double[] operand, final int operandOffset,
1206 final double[] result, final int resultOffset) {
1207
1208
1209 final double[] function = new double[1 + order];
1210 final double t = FastMath.tan(operand[operandOffset]);
1211 function[0] = t;
1212
1213 if (order > 0) {
1214
1215
1216
1217
1218
1219
1220
1221
1222 final double[] p = new double[order + 2];
1223 p[1] = 1;
1224 final double t2 = t * t;
1225 for (int n = 1; n <= order; ++n) {
1226
1227
1228 double v = 0;
1229 p[n + 1] = n * p[n];
1230 for (int k = n + 1; k >= 0; k -= 2) {
1231 v = v * t2 + p[k];
1232 if (k > 2) {
1233 p[k - 2] = (k - 1) * p[k - 1] + (k - 3) * p[k - 3];
1234 } else if (k == 2) {
1235 p[0] = p[1];
1236 }
1237 }
1238 if ((n & 0x1) == 0) {
1239 v *= t;
1240 }
1241
1242 function[n] = v;
1243
1244 }
1245 }
1246
1247
1248 compose(operand, operandOffset, function, result, resultOffset);
1249
1250 }
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260 public void acos(final double[] operand, final int operandOffset,
1261 final double[] result, final int resultOffset) {
1262
1263
1264 double[] function = new double[1 + order];
1265 final double x = operand[operandOffset];
1266 function[0] = FastMath.acos(x);
1267 if (order > 0) {
1268
1269
1270
1271
1272
1273
1274
1275 final double[] p = new double[order];
1276 p[0] = -1;
1277 final double x2 = x * x;
1278 final double f = 1.0 / (1 - x2);
1279 double coeff = FastMath.sqrt(f);
1280 function[1] = coeff * p[0];
1281 for (int n = 2; n <= order; ++n) {
1282
1283
1284 double v = 0;
1285 p[n - 1] = (n - 1) * p[n - 2];
1286 for (int k = n - 1; k >= 0; k -= 2) {
1287 v = v * x2 + p[k];
1288 if (k > 2) {
1289 p[k - 2] = (k - 1) * p[k - 1] + (2 * n - k) * p[k - 3];
1290 } else if (k == 2) {
1291 p[0] = p[1];
1292 }
1293 }
1294 if ((n & 0x1) == 0) {
1295 v *= x;
1296 }
1297
1298 coeff *= f;
1299 function[n] = coeff * v;
1300
1301 }
1302 }
1303
1304
1305 compose(operand, operandOffset, function, result, resultOffset);
1306
1307 }
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317 public void asin(final double[] operand, final int operandOffset,
1318 final double[] result, final int resultOffset) {
1319
1320
1321 double[] function = new double[1 + order];
1322 final double x = operand[operandOffset];
1323 function[0] = FastMath.asin(x);
1324 if (order > 0) {
1325
1326
1327
1328
1329
1330
1331
1332 final double[] p = new double[order];
1333 p[0] = 1;
1334 final double x2 = x * x;
1335 final double f = 1.0 / (1 - x2);
1336 double coeff = FastMath.sqrt(f);
1337 function[1] = coeff * p[0];
1338 for (int n = 2; n <= order; ++n) {
1339
1340
1341 double v = 0;
1342 p[n - 1] = (n - 1) * p[n - 2];
1343 for (int k = n - 1; k >= 0; k -= 2) {
1344 v = v * x2 + p[k];
1345 if (k > 2) {
1346 p[k - 2] = (k - 1) * p[k - 1] + (2 * n - k) * p[k - 3];
1347 } else if (k == 2) {
1348 p[0] = p[1];
1349 }
1350 }
1351 if ((n & 0x1) == 0) {
1352 v *= x;
1353 }
1354
1355 coeff *= f;
1356 function[n] = coeff * v;
1357
1358 }
1359 }
1360
1361
1362 compose(operand, operandOffset, function, result, resultOffset);
1363
1364 }
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374 public void atan(final double[] operand, final int operandOffset,
1375 final double[] result, final int resultOffset) {
1376
1377
1378 double[] function = new double[1 + order];
1379 final double x = operand[operandOffset];
1380 function[0] = FastMath.atan(x);
1381 if (order > 0) {
1382
1383
1384
1385
1386
1387
1388
1389 final double[] q = new double[order];
1390 q[0] = 1;
1391 final double x2 = x * x;
1392 final double f = 1.0 / (1 + x2);
1393 double coeff = f;
1394 function[1] = coeff * q[0];
1395 for (int n = 2; n <= order; ++n) {
1396
1397
1398 double v = 0;
1399 q[n - 1] = -n * q[n - 2];
1400 for (int k = n - 1; k >= 0; k -= 2) {
1401 v = v * x2 + q[k];
1402 if (k > 2) {
1403 q[k - 2] = (k - 1) * q[k - 1] + (k - 1 - 2 * n) * q[k - 3];
1404 } else if (k == 2) {
1405 q[0] = q[1];
1406 }
1407 }
1408 if ((n & 0x1) == 0) {
1409 v *= x;
1410 }
1411
1412 coeff *= f;
1413 function[n] = coeff * v;
1414
1415 }
1416 }
1417
1418
1419 compose(operand, operandOffset, function, result, resultOffset);
1420
1421 }
1422
1423
1424
1425
1426
1427
1428
1429
1430
1431
1432
1433 public void atan2(final double[] y, final int yOffset,
1434 final double[] x, final int xOffset,
1435 final double[] result, final int resultOffset) {
1436
1437
1438 double[] tmp1 = new double[getSize()];
1439 multiply(x, xOffset, x, xOffset, tmp1, 0);
1440 double[] tmp2 = new double[getSize()];
1441 multiply(y, yOffset, y, yOffset, tmp2, 0);
1442 add(tmp1, 0, tmp2, 0, tmp2, 0);
1443 rootN(tmp2, 0, 2, tmp1, 0);
1444
1445 if (x[xOffset] >= 0) {
1446
1447
1448 add(tmp1, 0, x, xOffset, tmp2, 0);
1449 divide(y, yOffset, tmp2, 0, tmp1, 0);
1450 atan(tmp1, 0, tmp2, 0);
1451 for (int i = 0; i < tmp2.length; ++i) {
1452 result[resultOffset + i] = 2 * tmp2[i];
1453 }
1454
1455 } else {
1456
1457
1458 subtract(tmp1, 0, x, xOffset, tmp2, 0);
1459 divide(y, yOffset, tmp2, 0, tmp1, 0);
1460 atan(tmp1, 0, tmp2, 0);
1461 result[resultOffset] =
1462 ((tmp2[0] <= 0) ? -FastMath.PI : FastMath.PI) - 2 * tmp2[0];
1463 for (int i = 1; i < tmp2.length; ++i) {
1464 result[resultOffset + i] = -2 * tmp2[i];
1465 }
1466
1467 }
1468
1469
1470 result[resultOffset] = FastMath.atan2(y[yOffset], x[xOffset]);
1471
1472 }
1473
1474
1475
1476
1477
1478
1479
1480
1481
1482 public void cosh(final double[] operand, final int operandOffset,
1483 final double[] result, final int resultOffset) {
1484
1485
1486 double[] function = new double[1 + order];
1487 function[0] = FastMath.cosh(operand[operandOffset]);
1488 if (order > 0) {
1489 function[1] = FastMath.sinh(operand[operandOffset]);
1490 for (int i = 2; i <= order; ++i) {
1491 function[i] = function[i - 2];
1492 }
1493 }
1494
1495
1496 compose(operand, operandOffset, function, result, resultOffset);
1497
1498 }
1499
1500
1501
1502
1503
1504
1505
1506
1507
1508 public void sinh(final double[] operand, final int operandOffset,
1509 final double[] result, final int resultOffset) {
1510
1511
1512 double[] function = new double[1 + order];
1513 function[0] = FastMath.sinh(operand[operandOffset]);
1514 if (order > 0) {
1515 function[1] = FastMath.cosh(operand[operandOffset]);
1516 for (int i = 2; i <= order; ++i) {
1517 function[i] = function[i - 2];
1518 }
1519 }
1520
1521
1522 compose(operand, operandOffset, function, result, resultOffset);
1523
1524 }
1525
1526
1527
1528
1529
1530
1531
1532
1533
1534 public void tanh(final double[] operand, final int operandOffset,
1535 final double[] result, final int resultOffset) {
1536
1537
1538 final double[] function = new double[1 + order];
1539 final double t = FastMath.tanh(operand[operandOffset]);
1540 function[0] = t;
1541
1542 if (order > 0) {
1543
1544
1545
1546
1547
1548
1549
1550
1551 final double[] p = new double[order + 2];
1552 p[1] = 1;
1553 final double t2 = t * t;
1554 for (int n = 1; n <= order; ++n) {
1555
1556
1557 double v = 0;
1558 p[n + 1] = -n * p[n];
1559 for (int k = n + 1; k >= 0; k -= 2) {
1560 v = v * t2 + p[k];
1561 if (k > 2) {
1562 p[k - 2] = (k - 1) * p[k - 1] - (k - 3) * p[k - 3];
1563 } else if (k == 2) {
1564 p[0] = p[1];
1565 }
1566 }
1567 if ((n & 0x1) == 0) {
1568 v *= t;
1569 }
1570
1571 function[n] = v;
1572
1573 }
1574 }
1575
1576
1577 compose(operand, operandOffset, function, result, resultOffset);
1578
1579 }
1580
1581
1582
1583
1584
1585
1586
1587
1588
1589 public void acosh(final double[] operand, final int operandOffset,
1590 final double[] result, final int resultOffset) {
1591
1592
1593 double[] function = new double[1 + order];
1594 final double x = operand[operandOffset];
1595 function[0] = FastMath.acosh(x);
1596 if (order > 0) {
1597
1598
1599
1600
1601
1602
1603
1604 final double[] p = new double[order];
1605 p[0] = 1;
1606 final double x2 = x * x;
1607 final double f = 1.0 / (x2 - 1);
1608 double coeff = FastMath.sqrt(f);
1609 function[1] = coeff * p[0];
1610 for (int n = 2; n <= order; ++n) {
1611
1612
1613 double v = 0;
1614 p[n - 1] = (1 - n) * p[n - 2];
1615 for (int k = n - 1; k >= 0; k -= 2) {
1616 v = v * x2 + p[k];
1617 if (k > 2) {
1618 p[k - 2] = (1 - k) * p[k - 1] + (k - 2 * n) * p[k - 3];
1619 } else if (k == 2) {
1620 p[0] = -p[1];
1621 }
1622 }
1623 if ((n & 0x1) == 0) {
1624 v *= x;
1625 }
1626
1627 coeff *= f;
1628 function[n] = coeff * v;
1629
1630 }
1631 }
1632
1633
1634 compose(operand, operandOffset, function, result, resultOffset);
1635
1636 }
1637
1638
1639
1640
1641
1642
1643
1644
1645
1646 public void asinh(final double[] operand, final int operandOffset,
1647 final double[] result, final int resultOffset) {
1648
1649
1650 double[] function = new double[1 + order];
1651 final double x = operand[operandOffset];
1652 function[0] = FastMath.asinh(x);
1653 if (order > 0) {
1654
1655
1656
1657
1658
1659
1660
1661 final double[] p = new double[order];
1662 p[0] = 1;
1663 final double x2 = x * x;
1664 final double f = 1.0 / (1 + x2);
1665 double coeff = FastMath.sqrt(f);
1666 function[1] = coeff * p[0];
1667 for (int n = 2; n <= order; ++n) {
1668
1669
1670 double v = 0;
1671 p[n - 1] = (1 - n) * p[n - 2];
1672 for (int k = n - 1; k >= 0; k -= 2) {
1673 v = v * x2 + p[k];
1674 if (k > 2) {
1675 p[k - 2] = (k - 1) * p[k - 1] + (k - 2 * n) * p[k - 3];
1676 } else if (k == 2) {
1677 p[0] = p[1];
1678 }
1679 }
1680 if ((n & 0x1) == 0) {
1681 v *= x;
1682 }
1683
1684 coeff *= f;
1685 function[n] = coeff * v;
1686
1687 }
1688 }
1689
1690
1691 compose(operand, operandOffset, function, result, resultOffset);
1692
1693 }
1694
1695
1696
1697
1698
1699
1700
1701
1702
1703 public void atanh(final double[] operand, final int operandOffset,
1704 final double[] result, final int resultOffset) {
1705
1706
1707 double[] function = new double[1 + order];
1708 final double x = operand[operandOffset];
1709 function[0] = FastMath.atanh(x);
1710 if (order > 0) {
1711
1712
1713
1714
1715
1716
1717
1718 final double[] q = new double[order];
1719 q[0] = 1;
1720 final double x2 = x * x;
1721 final double f = 1.0 / (1 - x2);
1722 double coeff = f;
1723 function[1] = coeff * q[0];
1724 for (int n = 2; n <= order; ++n) {
1725
1726
1727 double v = 0;
1728 q[n - 1] = n * q[n - 2];
1729 for (int k = n - 1; k >= 0; k -= 2) {
1730 v = v * x2 + q[k];
1731 if (k > 2) {
1732 q[k - 2] = (k - 1) * q[k - 1] + (2 * n - k + 1) * q[k - 3];
1733 } else if (k == 2) {
1734 q[0] = q[1];
1735 }
1736 }
1737 if ((n & 0x1) == 0) {
1738 v *= x;
1739 }
1740
1741 coeff *= f;
1742 function[n] = coeff * v;
1743
1744 }
1745 }
1746
1747
1748 compose(operand, operandOffset, function, result, resultOffset);
1749
1750 }
1751
1752
1753
1754
1755
1756
1757
1758
1759
1760
1761
1762 public void compose(final double[] operand, final int operandOffset, final double[] f,
1763 final double[] result, final int resultOffset) {
1764 for (int i = 0; i < compIndirection.length; ++i) {
1765 final int[][] mappingI = compIndirection[i];
1766 double r = 0;
1767 for (int j = 0; j < mappingI.length; ++j) {
1768 final int[] mappingIJ = mappingI[j];
1769 double product = mappingIJ[0] * f[mappingIJ[1]];
1770 for (int k = 2; k < mappingIJ.length; ++k) {
1771 product *= operand[operandOffset + mappingIJ[k]];
1772 }
1773 r += product;
1774 }
1775 result[resultOffset + i] = r;
1776 }
1777 }
1778
1779
1780
1781
1782
1783
1784
1785
1786 public double taylor(final double[] ds, final int dsOffset, final double ... delta)
1787 throws MathArithmeticException {
1788 double value = 0;
1789 for (int i = getSize() - 1; i >= 0; --i) {
1790 final int[] orders = getPartialDerivativeOrders(i);
1791 double term = ds[dsOffset + i];
1792 for (int k = 0; k < orders.length; ++k) {
1793 if (orders[k] > 0) {
1794 try {
1795 term *= FastMath.pow(delta[k], orders[k]) /
1796 CombinatoricsUtils.factorial(orders[k]);
1797 } catch (NotPositiveException e) {
1798
1799 throw new MathInternalError(e);
1800 }
1801 }
1802 }
1803 value += term;
1804 }
1805 return value;
1806 }
1807
1808
1809
1810
1811
1812 public void checkCompatibility(final DSCompiler compiler)
1813 throws DimensionMismatchException {
1814 if (parameters != compiler.parameters) {
1815 throw new DimensionMismatchException(parameters, compiler.parameters);
1816 }
1817 if (order != compiler.order) {
1818 throw new DimensionMismatchException(order, compiler.order);
1819 }
1820 }
1821
1822 }