1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.apache.commons.math4.legacy.analysis.integration.gauss;
18
19 import java.math.BigDecimal;
20 import java.math.MathContext;
21
22 import org.apache.commons.math4.legacy.core.Pair;
23
24
25
26
27
28
29
30
31
32
33
34 public class LegendreHighPrecisionRuleFactory extends BaseRuleFactory<BigDecimal> {
35
36 private final MathContext mContext;
37
38 private final BigDecimal two;
39
40 private final BigDecimal minusOne;
41
42 private final BigDecimal oneHalf;
43
44
45
46
47 public LegendreHighPrecisionRuleFactory() {
48 this(MathContext.DECIMAL128);
49 }
50
51
52
53
54 public LegendreHighPrecisionRuleFactory(MathContext mContext) {
55 this.mContext = mContext;
56 two = new BigDecimal("2", mContext);
57 minusOne = new BigDecimal("-1", mContext);
58 oneHalf = new BigDecimal("0.5", mContext);
59 }
60
61
62 @Override
63 protected Pair<BigDecimal[], BigDecimal[]> computeRule(int numberOfPoints) {
64 if (numberOfPoints == 1) {
65
66 return new Pair<>(new BigDecimal[] { BigDecimal.ZERO },
67 new BigDecimal[] { two });
68 }
69
70
71
72
73 final BigDecimal[] previousPoints = getRuleInternal(numberOfPoints - 1).getFirst();
74
75
76 final BigDecimal[] points = new BigDecimal[numberOfPoints];
77 final BigDecimal[] weights = new BigDecimal[numberOfPoints];
78
79
80 final int iMax = numberOfPoints / 2;
81 for (int i = 0; i < iMax; i++) {
82
83 BigDecimal a = (i == 0) ? minusOne : previousPoints[i - 1];
84
85 BigDecimal b = (iMax == 1) ? BigDecimal.ONE : previousPoints[i];
86
87 BigDecimal pma = BigDecimal.ONE;
88
89 BigDecimal pa = a;
90
91 BigDecimal pmb = BigDecimal.ONE;
92
93 BigDecimal pb = b;
94 for (int j = 1; j < numberOfPoints; j++) {
95 final BigDecimal b_two_j_p_1 = new BigDecimal(2 * j + 1, mContext);
96 final BigDecimal b_j = new BigDecimal(j, mContext);
97 final BigDecimal b_j_p_1 = new BigDecimal(j + 1, mContext);
98
99
100
101
102 BigDecimal tmp1 = a.multiply(b_two_j_p_1, mContext);
103 tmp1 = pa.multiply(tmp1, mContext);
104 BigDecimal tmp2 = pma.multiply(b_j, mContext);
105
106 BigDecimal ppa = tmp1.subtract(tmp2, mContext);
107 ppa = ppa.divide(b_j_p_1, mContext);
108
109
110
111
112 tmp1 = b.multiply(b_two_j_p_1, mContext);
113 tmp1 = pb.multiply(tmp1, mContext);
114 tmp2 = pmb.multiply(b_j, mContext);
115
116 BigDecimal ppb = tmp1.subtract(tmp2, mContext);
117 ppb = ppb.divide(b_j_p_1, mContext);
118
119 pma = pa;
120 pa = ppa;
121 pmb = pb;
122 pb = ppb;
123 }
124
125
126 BigDecimal c = a.add(b, mContext).multiply(oneHalf, mContext);
127
128 BigDecimal pmc = BigDecimal.ONE;
129
130 BigDecimal pc = c;
131 boolean done = false;
132 while (!done) {
133 BigDecimal tmp1 = b.subtract(a, mContext);
134 BigDecimal tmp2 = c.ulp().multiply(BigDecimal.TEN, mContext);
135 done = tmp1.compareTo(tmp2) <= 0;
136 pmc = BigDecimal.ONE;
137 pc = c;
138 for (int j = 1; j < numberOfPoints; j++) {
139 final BigDecimal b_two_j_p_1 = new BigDecimal(2 * j + 1, mContext);
140 final BigDecimal b_j = new BigDecimal(j, mContext);
141 final BigDecimal b_j_p_1 = new BigDecimal(j + 1, mContext);
142
143
144 tmp1 = c.multiply(b_two_j_p_1, mContext);
145 tmp1 = pc.multiply(tmp1, mContext);
146 tmp2 = pmc.multiply(b_j, mContext);
147
148 BigDecimal ppc = tmp1.subtract(tmp2, mContext);
149 ppc = ppc.divide(b_j_p_1, mContext);
150
151 pmc = pc;
152 pc = ppc;
153 }
154
155 if (!done) {
156 if (pa.signum() * pc.signum() <= 0) {
157 b = c;
158 pmb = pmc;
159 pb = pc;
160 } else {
161 a = c;
162 pma = pmc;
163 pa = pc;
164 }
165 c = a.add(b, mContext).multiply(oneHalf, mContext);
166 }
167 }
168 final BigDecimal nP = new BigDecimal(numberOfPoints, mContext);
169 BigDecimal tmp1 = pmc.subtract(c.multiply(pc, mContext), mContext);
170 tmp1 = tmp1.multiply(nP);
171 tmp1 = tmp1.pow(2, mContext);
172 BigDecimal tmp2 = c.pow(2, mContext);
173 tmp2 = BigDecimal.ONE.subtract(tmp2, mContext);
174 tmp2 = tmp2.multiply(two, mContext);
175 tmp2 = tmp2.divide(tmp1, mContext);
176
177 points[i] = c;
178 weights[i] = tmp2;
179
180 final int idx = numberOfPoints - i - 1;
181 points[idx] = c.negate(mContext);
182 weights[idx] = tmp2;
183 }
184
185
186
187
188 if ((numberOfPoints & 1) != 0) {
189 BigDecimal pmc = BigDecimal.ONE;
190 for (int j = 1; j < numberOfPoints; j += 2) {
191 final BigDecimal b_j = new BigDecimal(j, mContext);
192 final BigDecimal b_j_p_1 = new BigDecimal(j + 1, mContext);
193
194
195 pmc = pmc.multiply(b_j, mContext);
196 pmc = pmc.divide(b_j_p_1, mContext);
197 pmc = pmc.negate(mContext);
198 }
199
200
201 final BigDecimal nP = new BigDecimal(numberOfPoints, mContext);
202 BigDecimal tmp1 = pmc.multiply(nP, mContext);
203 tmp1 = tmp1.pow(2, mContext);
204 BigDecimal tmp2 = two.divide(tmp1, mContext);
205
206 points[iMax] = BigDecimal.ZERO;
207 weights[iMax] = tmp2;
208 }
209
210 return new Pair<>(points, weights);
211 }
212 }