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 RfFieldDuplication<T extends CalculusFieldElement<T>> extends FieldDuplication<T> {
29
30
31
32
33
34
35 RfFieldDuplication(final T x, final T y, final T z) {
36 super(x, y, z);
37 }
38
39
40 @Override
41 protected void initialMeanPoint(final T[] va) {
42 va[3] = va[0].add(va[1]).add(va[2]).divide(3.0);
43 }
44
45
46 @Override
47 protected T convergenceCriterion(final T r, final T max) {
48 return max.divide(FastMath.sqrt(FastMath.sqrt(FastMath.sqrt(r.multiply(3.0)))));
49 }
50
51
52 @Override
53 protected void update(final int m, final T[] vaM, final T[] sqrtM, final double fourM) {
54
55
56 final T lambdaA = sqrtM[0].multiply(sqrtM[1]);
57 final T lambdaB = sqrtM[0].multiply(sqrtM[2]);
58 final T lambdaC = sqrtM[1].multiply(sqrtM[2]);
59
60
61 vaM[0] = vaM[0].linearCombination(0.25, vaM[0], 0.25, lambdaA, 0.25, lambdaB, 0.25, lambdaC);
62 vaM[1] = vaM[1].linearCombination(0.25, vaM[1], 0.25, lambdaA, 0.25, lambdaB, 0.25, lambdaC);
63 vaM[2] = vaM[2].linearCombination(0.25, vaM[2], 0.25, lambdaA, 0.25, lambdaB, 0.25, lambdaC);
64 vaM[3] = vaM[3].linearCombination(0.25, vaM[3], 0.25, lambdaA, 0.25, lambdaB, 0.25, lambdaC);
65
66 }
67
68
69 @Override
70 protected T evaluate(final T[] va0, final T aM, final double fourM) {
71
72
73 final T inv = aM.multiply(fourM).reciprocal();
74 final T bigX = va0[3].subtract(va0[0]).multiply(inv);
75 final T bigY = va0[3].subtract(va0[1]).multiply(inv);
76 final T bigZ = bigX.add(bigY).negate();
77
78
79 final T e2 = bigX.multiply(bigY).subtract(bigZ.multiply(bigZ));
80 final T e3 = bigX.multiply(bigY).multiply(bigZ);
81
82 final T e2e2 = e2.multiply(e2);
83 final T e2e3 = e2.multiply(e3);
84 final T e3e3 = e3.multiply(e3);
85 final T e2e2e2 = e2e2.multiply(e2);
86
87
88
89 final T poly = e2e2e2.multiply(RfRealDuplication.E2_E2_E2).
90 add(e3e3.multiply(RfRealDuplication.E3_E3)).
91 add(e2e3.multiply(RfRealDuplication.E2_E3)).
92 add(e2e2.multiply(RfRealDuplication.E2_E2)).
93 add(e3.multiply(RfRealDuplication.E3)).
94 add(e2.multiply(RfRealDuplication.E2)).
95 add(RfRealDuplication.CONSTANT).
96 divide(RfRealDuplication.DENOMINATOR);
97 return poly.divide(FastMath.sqrt(aM));
98
99 }
100
101
102 @Override
103 public T integral() {
104 final T x = getVi(0);
105 final T y = getVi(1);
106 final T z = getVi(2);
107 if (x.isZero()) {
108 return completeIntegral(y, z);
109 } else if (y.isZero()) {
110 return completeIntegral(x, z);
111 } else if (z.isZero()) {
112 return completeIntegral(x, y);
113 } else {
114 return super.integral();
115 }
116 }
117
118
119
120
121
122
123 private T completeIntegral(final T x, final T y) {
124
125 T xM = x.sqrt();
126 T yM = y.sqrt();
127
128
129 for (int i = 1; i < RfRealDuplication.AGM_MAX; ++i) {
130
131 final T xM1 = xM;
132 final T yM1 = yM;
133
134
135 xM = xM1.add(yM1).multiply(0.5);
136
137
138 yM = xM1.multiply(yM1).sqrt();
139
140
141 if (xM.subtract(yM).norm() <= 4 * FastMath.ulp(xM).getReal()) {
142
143 break;
144 }
145
146 }
147
148 return xM.add(yM).reciprocal().multiply(xM.getPi());
149
150 }
151
152 }