001 /*
002 * $Id: Roots.java 3296 2010-08-26 17:30:55Z kredel $
003 */
004
005 package edu.jas.arith;
006
007
008 // import java.util.Random;
009 import java.math.MathContext;
010
011 import edu.jas.structure.Power;
012
013
014 /**
015 * Root computation algorithms. Roots for BigInteger and BigDecimals.
016 * @author Heinz Kredel
017 */
018 public class Roots {
019
020
021 /**
022 * Integer n-th root. Uses BigDecimal and newton iteration. R is the n-th
023 * root of A.
024 * @param A big integer.
025 * @param n long.
026 * @return the n-th root of A.
027 */
028 public static BigInteger root(BigInteger A, int n) {
029 if (n == 1) {
030 return A;
031 }
032 if (n == 2) {
033 return sqrt(A);
034 }
035 if (n < 1) {
036 throw new IllegalArgumentException("negative root not defined");
037 }
038 if (A == null || A.isZERO() || A.isONE()) {
039 return A;
040 }
041 // ensure enough precision
042 int s = A.val.bitLength() + 2;
043 MathContext mc = new MathContext(s);
044 //System.out.println("mc = " + mc);
045 BigDecimal Ap = new BigDecimal(A.val, mc);
046 //System.out.println("Ap = " + Ap);
047 BigDecimal Ar = root(Ap, n);
048 //System.out.println("Ar = " + Ar);
049 java.math.BigInteger RP = Ar.val.toBigInteger();
050 BigInteger R = new BigInteger(RP);
051 while (true) {
052 BigInteger P = Power.positivePower(R, n);
053 //System.out.println("P = " + P);
054 if (A.compareTo(P) >= 0) {
055 break;
056 } else {
057 R = R.subtract(BigInteger.ONE);
058 }
059 }
060 return R;
061 }
062
063
064 /**
065 * Integer square root. Uses BigDecimal and newton iteration. R is the
066 * square root of A.
067 * @param A big integer.
068 * @return the square root of A.
069 */
070 public static BigInteger sqrt(BigInteger A) {
071 if (A == null || A.isZERO() || A.isONE()) {
072 return A;
073 }
074 // ensure enough precision
075 int s = A.val.bitLength() + 2;
076 MathContext mc = new MathContext(s);
077 //System.out.println("mc = " + mc);
078 // newton iteration
079 BigDecimal Ap = new BigDecimal(A.val, mc);
080 //System.out.println("Ap = " + Ap);
081 BigDecimal Ar = sqrt(Ap);
082 //System.out.println("Ar = " + Ar);
083 java.math.BigInteger RP = Ar.val.toBigInteger();
084 BigInteger R = new BigInteger(RP);
085 while (true) {
086 BigInteger P = R.multiply(R);
087 //System.out.println("P = " + P);
088 if (A.compareTo(P) >= 0) {
089 break;
090 } else {
091 R = R.subtract(BigInteger.ONE);
092 }
093 }
094 return R;
095 }
096
097
098 /**
099 * Integer square root. Uses BigInteger only. R is the square root of A.
100 * @param A big integer.
101 * @return the square root of A.
102 */
103 public static BigInteger sqrtInt(BigInteger A) {
104 if (A == null || A.isZERO() || A.isONE()) {
105 return A;
106 }
107 int s = A.signum();
108 if (s < 0) {
109 throw new ArithmeticException("root of negative not defined");
110 }
111 if (s == 0) {
112 return A;
113 }
114 BigInteger R, R1, d;
115 int log2 = A.val.bitLength();
116 //System.out.println("A = " + A + ", log2 = " + log2);
117 int rootlog2 = log2 - log2 / 2;
118 R = new BigInteger(A.val.shiftRight(rootlog2));
119 //System.out.println("R = " + R + ", rootlog2 = " + rootlog2);
120 d = R;
121 while (!d.isZERO()) {
122 d = new BigInteger(d.val.shiftRight(1)); // div 2
123 R1 = R.sum(d);
124 s = A.compareTo(R1.multiply(R1));
125 if (s == 0) {
126 return R1;
127 }
128 if (s > 0) {
129 R = R1;
130 }
131 //System.out.println("R1 = " + R1);
132 //System.out.println("d = " + d);
133 }
134 while (true) {
135 R1 = R.sum(BigInteger.ONE);
136 //System.out.println("R1 = " + R1);
137 s = A.compareTo(R1.multiply(R1));
138 if (s == 0) {
139 return R1;
140 }
141 if (s > 0) {
142 R = R1;
143 }
144 if (s < 0) {
145 return R;
146 }
147 }
148 //return R;
149 }
150
151
152 /**
153 * Square root. R is the square root of A.
154 * @param A big decimal.
155 * @return the square root of A.
156 */
157 public static BigDecimal sqrt(BigDecimal A) {
158 if (A == null || A.isZERO() || A.isONE()) {
159 return A;
160 }
161 // for small A use root of inverse
162 if (A.val.compareTo(BigDecimal.ONE.val) < 0) {
163 BigDecimal Ap = A.inverse();
164 //System.out.println("A.inverse() = " + Ap);
165 Ap = sqrt(Ap);
166 return Ap.inverse();
167 }
168 MathContext mc = A.context;
169 // newton iteration
170 BigDecimal Ap = new BigDecimal(A.val, mc);
171 BigDecimal ninv = new BigDecimal(0.5, mc);
172 BigDecimal R1, R = Ap.multiply(ninv); // initial guess
173 BigDecimal d;
174 while (true) {
175 R1 = R.sum(Ap.divide(R));
176 R1 = R1.multiply(ninv); // div n
177 d = R.subtract(R1).abs();
178 R = R1;
179 if (d.val.compareTo(BigDecimal.ONE.val) <= 0) {
180 //System.out.println("d = " + d);
181 break;
182 }
183 }
184 return R;
185 }
186
187
188 /**
189 * N-th root. R is the n-th root of A.
190 * @param A big decimal.
191 * @param n long.
192 * @return the n-th root of A.
193 */
194 public static BigDecimal root(BigDecimal A, int n) {
195 if (n == 1) {
196 return A;
197 }
198 if (n == 2) {
199 return sqrt(A);
200 }
201 if (n < 1) {
202 throw new IllegalArgumentException("negative root not defined");
203 }
204 if (A == null || A.isZERO() || A.isONE()) {
205 return A;
206 }
207 // for small A use root of inverse
208 if (A.val.compareTo(BigDecimal.ONE.val) < 0) {
209 BigDecimal Ap = A.inverse();
210 //System.out.println("A.inverse() = " + Ap);
211 Ap = root(Ap, n);
212 return Ap.inverse();
213 }
214 // ensure enough precision
215 MathContext mc = A.context;
216 // newton iteration
217 BigDecimal Ap = A;
218 BigDecimal N = new BigDecimal(n, mc);
219 BigDecimal ninv = new BigDecimal(1.0 / n, mc);
220 BigDecimal nsub = new BigDecimal(1.0, mc); // because of precision
221 nsub = nsub.subtract(ninv);
222 //BigDecimal half = BigDecimal.ONE.sum(BigDecimal.ONE).inverse();
223 BigDecimal half = new BigDecimal(BigDecimal.ONE.val.divide(java.math.BigDecimal.TEN));
224 BigDecimal P, R1, R = Ap.multiply(ninv); // initial guess
225 BigDecimal d;
226 while (true) {
227 P = Power.positivePower(R, n - 1);
228 R1 = Ap.divide(P.multiply(N));
229 R1 = R.multiply(nsub).sum(R1);
230 d = R.subtract(R1).abs();
231 R = R1;
232 //if ( d.compareTo( BigDecimal.ONE ) <= 0 ) {
233 // System.out.println("d = " + d);
234 if (d.val.compareTo(half.val) <= 0) {
235 //System.out.println("d.val = " + d.val);
236 break;
237 }
238 //}
239 }
240 return R;
241 }
242
243 }