001/*
002 * $Id: Roots.java 4865 2014-07-05 19:27:53Z 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        if ( A.signum() < 0 ) {
160            throw new ArithmeticException("root of negative not defined");
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"); //e-13"); // TODO
173        int p = Math.max(mc.getPrecision(),java.math.MathContext.DECIMAL64.getPrecision());
174        //java.math.MathContext.UNLIMITED.getPrecision() == 0
175        eps = Power.positivePower(eps,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");
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"); //e-10"); // TODO
232        int p = Math.max(mc.getPrecision(),java.math.MathContext.DECIMAL64.getPrecision());
233        //java.math.MathContext.UNLIMITED.getPrecision() == 0
234        eps = Power.positivePower(eps,(p*2)/3);
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 = 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        BigDecimal r = a.re.abs().sum(a.abs().re);
273        BigDecimal t = new BigDecimal(2);
274        BigDecimal ti = new BigDecimal("0.5");
275        //BigDecimal u = r.divide(t);
276        BigDecimal u = r.multiply(ti);
277        BigDecimal v = Roots.sqrt(u);
278        //System.out.println("r = " + r + ", a = " + a);
279        //System.out.println("v = " + v + ", u = " + u);
280        if (a.re.signum() >= 0) {
281            return new BigDecimalComplex(v, a.im.divide(v.multiply(t)));
282        }
283        u = v;
284        if (a.im.signum() < 0) {
285            u = u.negate();
286        }
287        return new BigDecimalComplex(a.im.abs().divide(v.multiply(t)), u);
288    }
289
290}