How to build the sqrt(2) for many digits? Use the Newton approximation:
Newton's method is an extremely powerful technique—in general the convergence is quadratic: as the method converges on the root, the difference between the root and the approximation is squared (the number of accurate digits roughly doubles) at each step.
The idea of the method is as follows:
one starts with an initial guess which is reasonably close to the true root,
then the function is approximated by its tangent line,
and one computes the x-intercept of this tangent line.
This x-intercept will typically be a better approximation to the function's root than the original guess,
and the method can be iterated.
x = sqrt(2)
x²-2 = 0
f(x) = x²-2 =0
f'(x) = 2x
according to Newton:
x' = x - f(x)/f'(x)
= x - (x² -2) / 2x
= (2x² - x² +2) / 2x
= (x +2/x ) /2
Another approach:
Or let's think we start with an rectangle with area size A = 2.
+---------------------+
| |
a A |
| |
+---------x-----------+
There are 2 sides called a and x. One called x is initialized with 1.
So the other side a is defined as A / x.
Next we try to get a and x closer to the same size, so the rectangle will migrate to a square.
So an improved value for x can be setup as (a + x) / 2.
Finally the result will be x==a==sqrt(A)
x = 1 // init value
a = A / x // the other side
x'= (x+a) /2 // improved value for x, just as arithmetic average
x'= (x + A/x) /2
Code listing
/*
calc sqrt(2) with a lot of digits
trick:
start calcuation with low scale and improve scale with each loop
BigDecimal is based on BigIntger.
BigInteger use a BIT-Array of size 2^31 -1
So max value stored is 2^2.147.483.646
idx binary decimal
---------------------------------------------------
1 max bin=2^1 max decimal=10^0
2 max bin=2^3 max decimal=10^0
3 max bin=2^7 max decimal=10^2
4 max bin=2^15 max decimal=10^4
5 max bin=2^31 max decimal=10^9
6 max bin=2^63 max decimal=10^18
7 max bin=2^127 max decimal=10^38
8 max bin=2^255 max decimal=10^76
9 max bin=2^511 max decimal=10^153
10 max bin=2^1023 max decimal=10^307
11 max bin=2^2047 max decimal=10^616
12 max bin=2^4095 max decimal=10^1232
13 max bin=2^8191 max decimal=10^2465
14 max bin=2^16383 max decimal=10^4931
15 max bin=2^32767 max decimal=10^9863
16 max bin=2^65535 max decimal=10^19728
17 max bin=2^131071 max decimal=10^39456
18 max bin=2^262143 max decimal=10^78912
19 max bin=2^524287 max decimal=10^157826
20 max bin=2^1048575 max decimal=10^315652
21 max bin=2^2097151 max decimal=10^631305
22 max bin=2^4194303 max decimal=10^1262611
23 max bin=2^8388607 max decimal=10^2525222
24 max bin=2^16777215 max decimal=10^5050444
25 max bin=2^33554431 max decimal=10^10100890
26 max bin=2^67108863 max decimal=10^20201780
27 max bin=2^134217727 max decimal=10^40403561
28 max bin=2^268435455 max decimal=10^80807123
29 max bin=2^536870911 max decimal=10^161614248
30 max bin=2^1073741823 max decimal=10^323228496
31 max bin=2^2147483646 max decimal=10^646456992
wich means log(10) / log(2) * (2 ^ 31) -1 digits in decimal notation
mda = 646.456.992 (max digit amount)
if we want x*x=A, digit amount for x = sqrt(max digit amount)
= 323.228.496
*/
package sqrt;
import java.io.BufferedReader;
import java.io.BufferedWriter;
import java.io.FileInputStream;
import java.io.FileNotFoundException;
import java.io.FileOutputStream;
import java.io.IOException;
import java.io.InputStreamReader;
import java.io.OutputStreamWriter;
import java.math.BigDecimal;
import java.math.BigInteger;
import java.util.Map;
import java.util.Set;
/**
*
* @author LangR
*/
public class Sqrt {
static boolean verbose = false; // more output
static boolean onScreen = false;
static boolean chk = false; // compare solution with reference value
static boolean niceOutput = true; // use some nice spaces and linefeeds
static int maxDigit = 1000;
//--------------------------------------------------------------------------
private static String addXML(String tag, String value, String comment) {
if (comment.length() > 0) {
return "<" + tag + ">"
+ "<comment>" + comment + "</comment>"
+ value
+ "</" + tag + ">\n";
} else {
return "<" + tag + ">"
+ value
+ "</" + tag + ">\n";
}
}
//--------------------------------------------------------------------------
public static String check() throws FileNotFoundException, IOException {
FileInputStream fis = new FileInputStream("src/sqrt/reference");
BufferedReader buf = new BufferedReader(new InputStreamReader(fis));
String line = buf.readLine();
StringBuilder sb = new StringBuilder();
while (line != null) {
sb.append(line);
line = buf.readLine();
}
fis.close();
return sb.toString();
}
//--------------------------------------------------------------------------
public static void out(String fn, String x) throws FileNotFoundException, IOException {
StringBuffer sb = new StringBuffer();
int dotPos = x.indexOf(".") + 1;
for (int i = 0; i < x.length(); i++) {
if (niceOutput) {
if (i % 100 == dotPos && i > dotPos) {
sb.append("\n");
}
if (i % 1000 == dotPos && i > dotPos) {
sb.append("\n");
}
if (i % 10 == dotPos) {
sb.append(" ");
}
if (i == dotPos) {
sb.append("\n ");
}
}
sb.append(x.charAt(i));
}
if (fn.length() > 0) {
FileOutputStream fos = new FileOutputStream(fn);
BufferedWriter buf = new BufferedWriter(new OutputStreamWriter(fos));
buf.write(sb.toString());
buf.flush();
fos.close();
} else {
System.out.println(sb);
}
// output
if (chk) {
System.out.println("start to compare with reference string");
System.out.println(x.replace(".", ""));
if (check().contains(x.replace(".", ""))) {
System.out.println("ok");
} else {
System.out.println("Failed");
System.out.println(check());
}
System.exit(1);
}
System.out.println();
}
//--------------------------------------------------------------------------
public static void chkOverflow(int max) {
BigDecimal one = BigDecimal.ONE;
BigDecimal two = one.add(one);
System.out.println("Try to find a/b can handle " + max + " digits");
do {
System.out.println("try " + max + " next");
try {
one.add(one).divide(two, max, BigDecimal.ROUND_HALF_UP);
} catch (Exception ex) {
System.out.println("Error:\t" + max + "\t" + ex.getMessage());
chkOverflow(max - 1);
}
} while (max < maxDigit);
System.out.println("Try to find a/b can handle " + maxDigit + " digits");
System.out.println("ok");
System.exit(1);
}
//--------------------------------------------------------------------------
public static BigDecimal mySqrt(BigDecimal A, BigDecimal a) {
/*
*-----------------*
| A a
*---------------b-*
A = a * b
a'= (a + A/a) /2
*/
BigDecimal TWO = new BigDecimal("2.0");
BigDecimal aOld;
int iter = 0;
int scale = 1;
do {
if (verbose) {
System.out.println("\titeration \t" + (++iter)+"\tscale = "+scale);
}
aOld = a;
a = (a.add(
A.divide(a, scale, BigDecimal.ROUND_HALF_UP)))
.divide(TWO, scale, BigDecimal.ROUND_HALF_UP);
scale += scale; // inc precision
if (scale > maxDigit) {
scale = maxDigit;
}
} while (!aOld.equals(a)); // stop on no change
return a;
}
//--------------------------------------------------------------------------
private static void usage() {
System.out.println("optional arguments\n"
+ "-chk to compare sqrt(2) against reference 1E06\n"
+ "-v to verbose\n"
+ "-len to define the length of digits used\n"
+ "-screen to report on screen\n"
+ "-test to chk overflow condition\n"
+ "and a number n e.g. 3 for scale digits, where n as exponent of 10^n is used.");
System.exit(1);
}
//--------------------------------------------------------------------------
private static String Z(long n) {
// return a 2 digit number with leading 0
if (n < 10) {
return "0" + n;
} else {
return n + "";
}
}
//--------------------------------------------------------------------------
private static String outXML(String[] args, long timeMs) throws IOException {
StringBuffer sb = new StringBuffer();
sb.append("<root>\n");
sb.append(addXML("author", "R.LANG", ""));
StackTraceElement[] stack = Thread.currentThread().getStackTrace();
StackTraceElement main = stack[stack.length - 1];
String mainClass = main.getClassName();
sb.append(addXML("program", mainClass, ""));
for (int i = 0; i < args.length; i++) {
sb.append(addXML("argument_" + i, args[i], ""));
}
long secondsInMilli = 1000;
long minutesInMilli = secondsInMilli * 60;
long hoursInMilli = minutesInMilli * 60;
long daysInMilli = hoursInMilli * 24;
long elapsedDays = timeMs / daysInMilli;
timeMs = timeMs % daysInMilli;
long elapsedHours = timeMs / hoursInMilli;
timeMs = timeMs % hoursInMilli;
long elapsedMinutes = timeMs / minutesInMilli;
timeMs = timeMs % minutesInMilli;
long elapsedSeconds = timeMs / secondsInMilli;
timeMs = timeMs % secondsInMilli;
String s;
if (elapsedDays > 0) {
s = elapsedDays + " days " + Z(elapsedHours) + ":" + Z(elapsedMinutes) + ":" + Z(elapsedSeconds);
} else if (elapsedHours > 0) {
s = Z(elapsedHours) + ":" + Z(elapsedMinutes) + ":" + Z(elapsedSeconds);
} else if (elapsedMinutes > 0) {
s = Z(elapsedMinutes) + ":" + Z(elapsedSeconds);
} else if (elapsedSeconds > 0) {
s = Z(elapsedSeconds) + "s " + timeMs + "ms";;
} else {
s = timeMs + "ms";
}
sb.append(addXML("execution_time", s, ""));
String version = System.getProperty("java.version");
sb.append(addXML("jre_version", version, "JRE version used"));
sb.append("<hw>\n");
Runtime runtime = Runtime.getRuntime();
sb.append(addXML("available_processors", runtime.availableProcessors() + "", ""));
sb.append(addXML("max_memory", runtime.maxMemory() / 1024 / 1024 + " MByte", "Returns the maximum amount of memory that the Java virtual machine will attempt to use"));
sb.append(addXML("free_memory", runtime.freeMemory() / 1024 / 1024 + " MByte", "Returns the amount of free memory in the Java Virtual Machine"));
sb.append(addXML("total_memory_in_use", runtime.totalMemory() / 1024 / 1024 + " MByte", "Returns the total amount of memory in the Java virtual machine"));
sb.append(addXML("operating_system", System.getProperty("os.name").toLowerCase(), ""));
Map<String, String> map = System.getenv();
Set<String> keys = map.keySet();
for (String key : keys) {
if (key.startsWith("PROC")) {
sb.append(addXML(key.toLowerCase(), map.get(key), ""));
}
}
sb.append("</hw>\n");
sb.append("</root>\n");
return sb.toString();
}
//--------------------------------------------------------------------------
public static void main(String[] args) throws FileNotFoundException, IOException {
double f = Math.log(10.0) / Math.log(2.0);
/*
for (int i = 1; i < 32; i++) {
int exp= (int) Math.pow(2, i)-1;
BigInteger b = BigInteger.valueOf(2).pow(exp);
//BigDecimal d = BigDecimal.valueOf(10).pow((int) (exp/f)) ;
System.out.println(i + "\tmax bin=2^" + (long) exp +"\tmax decimal=10^"+(long) (exp/f) );
}
*/
long start = System.currentTimeMillis();
String mantisseArg = "2";
for (int i = 0; i < args.length; i++) {
if (args[i].contains("-v")) {
verbose = true;
} else if (args[i].contains("-root")) {
i++;
mantisseArg = args[i];
} else if (args[i].contains("-screen")) {
onScreen = true;
} else if (args[i].contains("-chk")) {
chk = true;
} else if (args[i].contains("-len")) {
i++;
maxDigit = new Integer(args[i]);
} else if (args[i].contains("-test")) {
chkOverflow(536870919);
} else {
usage();
}
}
int exp = (int) (maxDigit * f);
try {
BigDecimal b = BigDecimal.valueOf(2).pow(exp);
} catch (Exception ex) {
System.out.println(ex.getMessage());
System.exit(1);
}
BigDecimal A = new BigDecimal(mantisseArg);
BigDecimal x = mySqrt(A, BigDecimal.ONE);
long finish = System.currentTimeMillis();
String fn = "sqrt(" + mantisseArg + ")_" + maxDigit;
String sqrt = x.toString(); // this conversion needs time!!
if (onScreen) {
out("", sqrt); // on screen
} else {
out(fn + ".txt", sqrt);
niceOutput = false;
out(fn + ".log", sqrt);
String s = outXML(args, finish - start);
out(fn + ".xml", s);
}
System.out.println("Digit length = " + sqrt.length() + "\tcalc time = " + (finish - start) + "ms");
}
}