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}