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}