1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.hipparchus.special.elliptic.carlson;
18
19 import org.hipparchus.CalculusFieldElement;
20 import org.hipparchus.complex.Complex;
21 import org.hipparchus.complex.FieldComplex;
22 import org.hipparchus.util.FastMath;
23
24
25
26
27
28 class RjFieldDuplication<T extends CalculusFieldElement<T>> extends FieldDuplication<T> {
29
30
31 private T delta;
32
33
34 private T sM;
35
36
37
38
39
40
41
42
43 RjFieldDuplication(final T x, final T y, final T z, final T p, final T delta) {
44 super(x, y, z, p);
45 this.delta = delta;
46 }
47
48
49 @Override
50 protected void initialMeanPoint(final T[] va) {
51 va[4] = va[0].add(va[1]).add(va[2]).add(va[3].multiply(2)).divide(5.0);
52 }
53
54
55 @Override
56 protected T convergenceCriterion(final T r, final T max) {
57 return max.divide(FastMath.sqrt(FastMath.sqrt(FastMath.sqrt(r.multiply(0.25)))));
58 }
59
60
61 @Override
62 protected void update(final int m, final T[] vaM, final T[] sqrtM, final double fourM) {
63 final T dM = sqrtM[3].add(sqrtM[0]).
64 multiply(sqrtM[3].add(sqrtM[1])).
65 multiply(sqrtM[3].add(sqrtM[2]));
66 if (m == 0) {
67 sM = dM.multiply(0.5);
68 } else {
69
70 final T rM = sM.multiply(delta.divide(sM.multiply(sM).multiply(fourM)).add(1.0).sqrt().add(1.0));
71 sM = dM.multiply(rM).subtract(delta.divide(fourM * fourM)).
72 divide(dM.add(rM.divide(fourM)).multiply(2));
73 }
74
75
76 final T lambdaA = sqrtM[0].multiply(sqrtM[1]);
77 final T lambdaB = sqrtM[0].multiply(sqrtM[2]);
78 final T lambdaC = sqrtM[1].multiply(sqrtM[2]);
79
80
81 vaM[0] = vaM[0].linearCombination(0.25, vaM[0], 0.25, lambdaA, 0.25, lambdaB, 0.25, lambdaC);
82 vaM[1] = vaM[1].linearCombination(0.25, vaM[1], 0.25, lambdaA, 0.25, lambdaB, 0.25, lambdaC);
83 vaM[2] = vaM[2].linearCombination(0.25, vaM[2], 0.25, lambdaA, 0.25, lambdaB, 0.25, lambdaC);
84 vaM[3] = vaM[3].linearCombination(0.25, vaM[3], 0.25, lambdaA, 0.25, lambdaB, 0.25, lambdaC);
85 vaM[4] = vaM[4].linearCombination(0.25, vaM[4], 0.25, lambdaA, 0.25, lambdaB, 0.25, lambdaC);
86
87 }
88
89
90 @Override
91 protected T evaluate(final T[] va0, final T aM, final double fourM) {
92
93
94 final T inv = aM.multiply(fourM).reciprocal();
95 final T bigX = va0[4].subtract(va0[0]).multiply(inv);
96 final T bigY = va0[4].subtract(va0[1]).multiply(inv);
97 final T bigZ = va0[4].subtract(va0[2]).multiply(inv);
98 final T bigP = bigX.add(bigY).add(bigZ).multiply(-0.5);
99 final T bigP2 = bigP.multiply(bigP);
100
101
102 final T xyz = bigX.multiply(bigY).multiply(bigZ);
103 final T e2 = bigX.multiply(bigY.add(bigZ)).add(bigY.multiply(bigZ)).
104 subtract(bigP.multiply(bigP).multiply(3));
105 final T e3 = xyz.add(bigP.multiply(2).multiply(e2.add(bigP2.multiply(2))));
106 final T e4 = xyz.multiply(2).add(bigP.multiply(e2.add(bigP2.multiply(3)))).multiply(bigP);
107 final T e5 = xyz.multiply(bigP2);
108
109 final T e2e2 = e2.multiply(e2);
110 final T e2e3 = e2.multiply(e3);
111 final T e2e4 = e2.multiply(e4);
112 final T e2e5 = e2.multiply(e5);
113 final T e3e3 = e3.multiply(e3);
114 final T e3e4 = e3.multiply(e4);
115 final T e2e2e2 = e2e2.multiply(e2);
116 final T e2e2e3 = e2e2.multiply(e3);
117
118
119
120 final T poly = e3e4.add(e2e5).multiply(RdRealDuplication.E3_E4_P_E2_E5).
121 add(e2e2e3.multiply(RdRealDuplication.E2_E2_E3)).
122 add(e2e4.multiply(RdRealDuplication.E2_E4)).
123 add(e3e3.multiply(RdRealDuplication.E3_E3)).
124 add(e2e2e2.multiply(RdRealDuplication.E2_E2_E2)).
125 add(e5.multiply(RdRealDuplication.E5)).
126 add(e2e3.multiply(RdRealDuplication.E2_E3)).
127 add(e4.multiply(RdRealDuplication.E4)).
128 add(e2e2.multiply(RdRealDuplication.E2_E2)).
129 add(e3.multiply(RdRealDuplication.E3)).
130 add(e2.multiply(RdRealDuplication.E2)).
131 add(RdRealDuplication.CONSTANT).
132 divide(RdRealDuplication.DENOMINATOR);
133 final T polyTerm = poly.divide(aM.multiply(FastMath.sqrt(aM)).multiply(fourM));
134
135
136 final T rcTerm = new RcFieldDuplication<>(poly.getField().getOne(), delta.divide(sM.multiply(sM).multiply(fourM)).add(1)).
137 integral().
138 multiply(3).divide(sM);
139
140 return polyTerm.add(rcTerm);
141
142 }
143
144 }