CS3343/3341 Analysis of Algorithms Weird Topic |
Newton's Method for
Arbitrary Square Roots |
The program below incorporates normalization and calculates the square root of an arbitrary positive double.
For convenience the input double is given as a command line argument. In the several example outputs, there is one very small positive number and one very large one. The output shown is from the Java program. The output from the C program differs only in trivial ways. In the data below, black digits are correct, while red digits are incorrect.| Newton's Method, Square Root of the double r | |
|---|---|
| Root.java | root.c |
// Root.java: use Newton's method for sqrt(r)
// includes normalization.
public class Root {
private double mul; // denormalizing factor
// norm: normalize to 0.25 <= r <= 1.0
public double norm(double r) {
// divide d by even power of 2
mul = 1.0; // power of 2
double r0 = r;
while (r0 > 1) {
// divide by 4, will
// later multiply by 2
r0 = r0/4.0;
mul = mul*2.0;
}
while (r0 < 0.25) {
r0 = r0*4.0;
mul = mul/2.0;
}
// r0 now normalized
return r0;
}
// denorm: reverse the normalization
public double denorm(double sqrtR0) {
return sqrtR0*mul;
}
// calculate sqrt(r0), 0.25 <= d0 <= 1
public double sqrt0(double r0) {
// initial approx: x0 = (2/3)r0 + 0.3506
double x0 = 0.66666*r0 + 0.3506;
double x1;
System.out.println(
"1st approx to sqrt(r0): " + x0);
for ( ; ; ) {
x1 = x0/2.0 + r0/(2.0*x0);
System.out.println(x1);
if (Math.abs(x1 - x0) < 1.0e-7) break;
x0 = x1;
}
return x1;
}
public static void main(String[] args) {
// calculate sqrt(r) for any r > 0
Root root = new Root();
double r = Double.parseDouble(args[0]);
System.out.println("input is r: " + r +
", want sqrt(r)");
// normalize to 0.25 <= r<= 1
double r0 = root.norm(r);
// r0 now normalized
System.out.println("normalized is r0: " +
r0 + ", want sqrt(r0)");
double sqrtR0 = root.sqrt0(r0);
// denoramlize: multiply sqrtR0 by mul
double sqrtR = root.denorm(sqrtR0);
System.out.println("sqrt(" + r + ") = " +
sqrtR);
System.out.println("check on ans: " +
Math.sqrt(r));
}
}
| // root.c: use Newton's method for sqrt(r)
// includes normalization.
#include <stdio.h>
#include <math.h>
#include<stdlib.h>
double mul; // denormalizing factor
// norm: normalize to 0.25 <= r <= 1.0
double norm(double r) {
// divide d by even power of 2
mul = 1.0; // power of 2
double r0 = r;
while (r0 > 1) {
// divide by 4, will
// later multiply by 2
r0 = r0/4.0;
mul = mul*2.0;
}
while (r0 < 0.25) {
r0 = r0*4.0;
mul = mul/2.0;
}
// r0 now normalized
return r0;
}
// denorm: reverse the normalization
double denorm(double sqrtr0) {
return sqrtr0*mul;
}
// calculate sqrt(r0), 0.25 <= d0 <= 1
double sqrt0(double r0) {
// initial approx: x0 = (2/3)r0 + 0.3506
double x0 = 0.66666*r0 + 0.3506;
double x1;
printf("1st approx to sqrt(r0): %18.16f\n",
x0);
for ( ; ; ) {
x1 = x0/2.0 + r0/(2.0*x0);
printf("%18.16f\n", x1);
if (fabs(x1 - x0) < 1.0e-7) break;
x0 = x1;
}
return x1;
}
void main(int argc, char *argv[]) {
// calculate sqrt(r) for any r > 0
double r, r0, sqrtr0, sqrtr;
if (argc != 2) {
printf("Usage: \"./root 2\"\n");
exit(1);
}
sscanf(argv[1], "%lf", &r);
printf("input is r: %18.16f, want sqrt(r)\n",
r);
// normalize to 0.25 <= r<= 1
r0 = norm(r);
// r0 now normalized
printf("normalized is r0: %18.16f", r0);
printf(", want sqrt(r0)\n");
sqrtr0 = sqrt0(r0);
// denoramlize: multipy sqrtR0 by mul
sqrtr = denorm(sqrtr0);
printf("sqrt(%18.16f) = %18.16f\n",
r, sqrtr);
printf("check on ans: %18.16f\n", sqrt(r));
}
|
% javac Root.java % java Root 2.0 input is r: 2.0, want sqrt(r) normalized is r0: 0.5, want sqrt(r0) 1st approx to sqrt(r0): 0.68393 0.7074994845232699 0.7071068901731362 0.7071067811865559 0.7071067811865475 sqrt(2.0) = 1.414213562373095 check on ans: 1.4142135623730951 | % cc -o root -lm root.c % ./root 10 input is r: 2.012E-7, want sqrt(r) normalized is r0: 0.8438939648, want sqrt(r0) 1st approx to sqrt(r0): 0.913190350573568 0.918653257848627 0.9186370148561589 0.9186370147125578 sqrt(2.012E-7) = 4.485532298401161E-4 check on ans: 4.485532298401161E-4 |