Source code

001/*
002 * $Id$
003 */
004
005package edu.jas.arith;
006
007
008import java.math.MathContext;
009
010
011/**
012 * Root computation algorithms. Roots for BigInteger and BigDecimals.
013 * @author Heinz Kredel
014 */
015public class Roots {
016
017
018    /**
019     * Integer n-th root. Uses BigDecimal and newton iteration. R is the n-th
020     * root of A.
021     * @param A big integer.
022     * @param n long.
023     * @return the n-th root of A.
024     */
025    public static BigInteger root(BigInteger A, int n) {
026        if (n == 1) {
027            return A;
028        }
029        if (n == 2) {
030            return sqrt(A);
031        }
032        if (n < 1) {
033            throw new IllegalArgumentException("negative root not defined");
034        }
035        if (A == null || A.isZERO() || A.isONE()) {
036            return A;
037        }
038        if (A.signum() < 0) {
039            throw new ArithmeticException("root of negative not defined: " + 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 = R.power(n); //Power.positivePower(R, n);
053            //System.out.println("P = " + P);
054            if (A.compareTo(P) >= 0) {
055                break;
056            }
057            R = R.subtract(BigInteger.ONE);
058        }
059        return R;
060    }
061
062
063    /**
064     * Integer square root. Uses BigDecimal and newton iteration. R is the
065     * square root of A.
066     * @param A big integer.
067     * @return the square root of A.
068     */
069    public static BigInteger sqrt(BigInteger A) {
070        if (A == null || A.isZERO() || A.isONE()) {
071            return A;
072        }
073        if (A.signum() < 0) {
074            throw new ArithmeticException("root of negative not defined: " + A);
075        }
076        // ensure enough precision
077        int s = A.val.bitLength() + 2;
078        MathContext mc = new MathContext(s);
079        //System.out.println("mc = " + mc);
080        // newton iteration
081        BigDecimal Ap = new BigDecimal(A.val, mc);
082        //System.out.println("Ap = " + Ap);
083        BigDecimal Ar = sqrt(Ap);
084        //System.out.println("Ar = " + Ar);
085        java.math.BigInteger RP = Ar.val.toBigInteger();
086        BigInteger R = new BigInteger(RP);
087        while (true) {
088            BigInteger P = R.multiply(R);
089            //System.out.println("P = " + P);
090            if (A.compareTo(P) >= 0) {
091                break;
092            }
093            R = R.subtract(BigInteger.ONE);
094        }
095        return R;
096    }
097
098
099    /**
100     * Integer square root. Uses BigInteger only. R is the square root of A.
101     * @param A big integer.
102     * @return the square root of A.
103     */
104    public static BigInteger sqrtInt(BigInteger A) {
105        if (A == null || A.isZERO() || A.isONE()) {
106            return A;
107        }
108        int s = A.signum();
109        if (s < 0) {
110            throw new ArithmeticException("root of negative not defined: " + A);
111        }
112        BigInteger R, R1, d;
113        int log2 = A.val.bitLength();
114        //System.out.println("A = " + A + ", log2 = " + log2);
115        int rootlog2 = log2 - log2 / 2;
116        R = new BigInteger(A.val.shiftRight(rootlog2));
117        //System.out.println("R = " + R + ", rootlog2 = " + rootlog2);
118        d = R;
119        while (!d.isZERO()) {
120            d = new BigInteger(d.val.shiftRight(1)); // div 2
121            R1 = R.sum(d);
122            s = A.compareTo(R1.multiply(R1));
123            if (s == 0) {
124                return R1;
125            }
126            if (s > 0) {
127                R = R1;
128            }
129            //System.out.println("R1 = " + R1);
130            //System.out.println("d  = " + d);
131        }
132        while (true) {
133            R1 = R.sum(BigInteger.ONE);
134            //System.out.println("R1 = " + R1);
135            s = A.compareTo(R1.multiply(R1));
136            if (s == 0) {
137                return R1;
138            }
139            if (s > 0) {
140                R = R1;
141            }
142            if (s < 0) {
143                return R;
144            }
145        }
146        //return R;
147    }
148
149
150    /**
151     * Square root. R is the square root of A.
152     * @param A big decimal.
153     * @return the square root of A.
154     */
155    public static BigDecimal sqrt(BigDecimal A) {
156        if (A == null || A.isZERO() || A.isONE()) {
157            return A;
158        }
159        if (A.signum() < 0) {
160            throw new ArithmeticException("root of negative not defined: " + A);
161        }
162        // for small A use root of inverse
163        if (A.abs().val.compareTo(BigDecimal.ONE.val) < 0) {
164            BigDecimal Ap = A.inverse();
165            Ap = sqrt(Ap);
166            Ap = Ap.inverse();
167            //System.out.println("sqrt(A).inverse() = " + Ap);
168            return Ap;
169        }
170        // ensure enough precision
171        MathContext mc = A.context;
172        BigDecimal eps = new BigDecimal("0.1");
173        int p = Math.max(mc.getPrecision(), java.math.MathContext.DECIMAL64.getPrecision());
174        //java.math.MathContext.UNLIMITED.getPrecision() == 0
175        eps = eps.power(p / 2);
176        // newton iteration
177        BigDecimal Ap = new BigDecimal(A.val, mc);
178        BigDecimal ninv = new BigDecimal(0.5, mc);
179        BigDecimal R1, R = Ap.multiply(ninv); // initial guess
180        BigDecimal d;
181        int i = 0;
182        while (true) {
183            R1 = R.sum(Ap.divide(R));
184            R1 = R1.multiply(ninv); // div n
185            d = R.subtract(R1).abs();
186            R = R1;
187            if (d.val.compareTo(eps.val) <= 0) {
188                //System.out.println("d  = " + d + ", R = " + R);
189                break;
190            }
191            if (i++ % 11 == 0) {
192                eps = eps.sum(eps);
193            }
194            //System.out.println("eps  = " + eps + ", d = " + d);
195        }
196        return R;
197    }
198
199
200    /**
201     * N-th root. R is the n-th root of A.
202     * @param A big decimal.
203     * @param n long.
204     * @return the n-th root of A.
205     */
206    public static BigDecimal root(BigDecimal A, int n) {
207        if (n == 1) {
208            return A;
209        }
210        if (n == 2) {
211            return sqrt(A);
212        }
213        if (n < 1) {
214            throw new IllegalArgumentException("negative root not defined");
215        }
216        if (A == null || A.isZERO() || A.isONE()) {
217            return A;
218        }
219        if (A.signum() < 0) {
220            throw new ArithmeticException("root of negative not defined: " + A);
221        }
222        // for small A use root of inverse
223        if (A.abs().val.compareTo(BigDecimal.ONE.val) < 0) {
224            BigDecimal Ap = A.inverse();
225            //System.out.println("A.inverse() = " + Ap);
226            Ap = root(Ap, n);
227            return Ap.inverse();
228        }
229        // ensure enough precision
230        MathContext mc = A.context;
231        BigDecimal eps = new BigDecimal("0.1");
232        int p = Math.max(mc.getPrecision(), java.math.MathContext.DECIMAL64.getPrecision());
233        //java.math.MathContext.UNLIMITED.getPrecision() == 0
234        eps = eps.power((p / 3) * 2);
235        // newton iteration
236        BigDecimal Ap = A;
237        BigDecimal N = new BigDecimal(n, mc);
238        BigDecimal ninv = new BigDecimal(1.0 / n, mc);
239        BigDecimal nsub = new BigDecimal(1.0, mc); // because of precision
240        nsub = nsub.subtract(ninv);
241        BigDecimal P, R1, R = Ap.multiply(ninv); // initial guess
242        BigDecimal d;
243        int i = 0;
244        while (true) {
245            P = R.power(n - 1); //Power.positivePower(R, n - 1);
246            R1 = Ap.divide(P.multiply(N));
247            R1 = R.multiply(nsub).sum(R1);
248            d = R.subtract(R1).abs();
249            R = R1;
250            if (d.val.compareTo(eps.val) <= 0) {
251                //System.out.println("d.val  = " + d.val);
252                break;
253            }
254            if (i++ % 11 == 0) {
255                eps = eps.sum(eps);
256            }
257        }
258        // System.out.println("eps  = " + eps + ", d = " + d);
259        return R;
260    }
261
262
263    /**
264     * Complex decimal number square root.
265     * @param a big decimal complex.
266     * @return sqrt(a).
267     */
268    public static BigDecimalComplex sqrt(BigDecimalComplex a) {
269        if (a.isZERO() || a.isONE()) {
270            return a;
271        }
272        if (a.im.isZERO() && a.re.signum() > 0) {
273            BigDecimal v = Roots.sqrt(a.re);
274            return new BigDecimalComplex(v);
275        }
276        BigDecimal r = a.re.abs().sum(a.abs().re);
277        BigDecimal t = new BigDecimal(2);
278        BigDecimal ti = new BigDecimal("0.5");
279        //BigDecimal u = r.divide(t);
280        BigDecimal u = r.multiply(ti);
281        BigDecimal v = Roots.sqrt(u);
282        //System.out.println("r = " + r + ", a = " + a);
283        //System.out.println("v = " + v + ", u = " + u);
284        if (a.re.signum() >= 0) {
285            return new BigDecimalComplex(v, a.im.divide(v.multiply(t)));
286        }
287        u = v;
288        if (a.im.signum() < 0) {
289            u = u.negate();
290        }
291        return new BigDecimalComplex(a.im.abs().divide(v.multiply(t)), u);
292    }
293
294
295    /**
296     * Square root. R is the square root approximation of A.
297     * Convert to BigDecimal and compute square root.
298     * @param A big rational.
299     * @return the square root approximation of A.
300     */
301    public static BigRational sqrt(BigRational A) {
302        if (A == null || A.isZERO() || A.isONE()) {
303            return A;
304        }
305        if (A.signum() < 0) {
306            throw new ArithmeticException("root of negative not defined: " + A);
307        }
308        BigDecimal ad = new BigDecimal(A);
309        BigDecimal s = sqrt(ad);
310        //System.out.println("s = " + s);
311        BigRational S = new BigRational(s.toString());
312        return S;
313    }
314
315
316    /**
317     * Square root. R is the square root approximation of A.
318     * Convert to BigDecimalComplex and compute square root.
319     * @param A big complex rational.
320     * @return the square root approximation of A.
321     */
322    public static BigComplex sqrt(BigComplex A) {
323        if (A == null || A.isZERO() || A.isONE()) {
324            return A;
325        }
326        BigDecimalComplex ad = new BigDecimalComplex(A);
327        BigDecimalComplex s = sqrt(ad);
328        //System.out.println("s = " + s);
329        BigComplex S = new BigComplex(s.toString());
330        return S;
331    }
332
333}