001/*
002 * $Id: Roots.java 4054 2012-07-26 17:34:57Z kredel $
003 */
004
005package edu.jas.arith;
006
007
008// import java.util.Random;
009import java.math.MathContext;
010
011import edu.jas.structure.Power;
012
013
014/**
015 * Root computation algorithms. Roots for BigInteger and BigDecimals.
016 * @author Heinz Kredel
017 */
018public 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            }
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        // ensure enough precision
074        int s = A.val.bitLength() + 2;
075        MathContext mc = new MathContext(s);
076        //System.out.println("mc = " + mc);
077        // newton iteration
078        BigDecimal Ap = new BigDecimal(A.val, mc);
079        //System.out.println("Ap = " + Ap);
080        BigDecimal Ar = sqrt(Ap);
081        //System.out.println("Ar = " + Ar);
082        java.math.BigInteger RP = Ar.val.toBigInteger();
083        BigInteger R = new BigInteger(RP);
084        while (true) {
085            BigInteger P = R.multiply(R);
086            //System.out.println("P = " + P);
087            if (A.compareTo(P) >= 0) {
088                break;
089            }
090            R = R.subtract(BigInteger.ONE);
091        }
092        return R;
093    }
094
095
096    /**
097     * Integer square root. Uses BigInteger only. R is the square root of A.
098     * @param A big integer.
099     * @return the square root of A.
100     */
101    public static BigInteger sqrtInt(BigInteger A) {
102        if (A == null || A.isZERO() || A.isONE()) {
103            return A;
104        }
105        int s = A.signum();
106        if (s < 0) {
107            throw new ArithmeticException("root of negative not defined");
108        }
109        if (s == 0) {
110            return 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        // for small A use root of inverse
160        if (A.val.compareTo(BigDecimal.ONE.val) < 0) {
161            BigDecimal Ap = A.inverse();
162            //System.out.println("A.inverse() = " + Ap);
163            Ap = sqrt(Ap);
164            return Ap.inverse();
165        }
166        MathContext mc = A.context;
167        // newton iteration
168        BigDecimal Ap = new BigDecimal(A.val, mc);
169        BigDecimal ninv = new BigDecimal(0.5, mc);
170        BigDecimal R1, R = Ap.multiply(ninv); // initial guess
171        BigDecimal d;
172        while (true) {
173            R1 = R.sum(Ap.divide(R));
174            R1 = R1.multiply(ninv); // div n
175            d = R.subtract(R1).abs();
176            R = R1;
177            if (d.val.compareTo(BigDecimal.ONE.val) <= 0) {
178                //System.out.println("d  = " + d);
179                break;
180            }
181        }
182        return R;
183    }
184
185
186    /**
187     * N-th root. R is the n-th root of A.
188     * @param A big decimal.
189     * @param n long.
190     * @return the n-th root of A.
191     */
192    public static BigDecimal root(BigDecimal A, int n) {
193        if (n == 1) {
194            return A;
195        }
196        if (n == 2) {
197            return sqrt(A);
198        }
199        if (n < 1) {
200            throw new IllegalArgumentException("negative root not defined");
201        }
202        if (A == null || A.isZERO() || A.isONE()) {
203            return A;
204        }
205        // for small A use root of inverse
206        if (A.val.compareTo(BigDecimal.ONE.val) < 0) {
207            BigDecimal Ap = A.inverse();
208            //System.out.println("A.inverse() = " + Ap);
209            Ap = root(Ap, n);
210            return Ap.inverse();
211        }
212        // ensure enough precision
213        MathContext mc = A.context;
214        // newton iteration
215        BigDecimal Ap = A;
216        BigDecimal N = new BigDecimal(n, mc);
217        BigDecimal ninv = new BigDecimal(1.0 / n, mc);
218        BigDecimal nsub = new BigDecimal(1.0, mc); // because of precision
219        nsub = nsub.subtract(ninv);
220        //BigDecimal half = BigDecimal.ONE.sum(BigDecimal.ONE).inverse();
221        BigDecimal half = new BigDecimal(BigDecimal.ONE.val.divide(java.math.BigDecimal.TEN));
222        BigDecimal P, R1, R = Ap.multiply(ninv); // initial guess
223        BigDecimal d;
224        while (true) {
225            P = Power.positivePower(R, n - 1);
226            R1 = Ap.divide(P.multiply(N));
227            R1 = R.multiply(nsub).sum(R1);
228            d = R.subtract(R1).abs();
229            R = R1;
230            //if ( d.compareTo( BigDecimal.ONE ) <= 0 ) {
231            //    System.out.println("d  = " + d);
232            if (d.val.compareTo(half.val) <= 0) {
233                //System.out.println("d.val  = " + d.val);
234                break;
235            }
236            //}
237        }
238        return R;
239    }
240
241}