Xref: utzoo comp.theory.dynamic-sys:186 sci.math:15784 sci.physics:17509
Path: utzoo!news-server.csri.toronto.edu!cs.utexas.edu!swrinde!elroy.jpl.nasa.gov!decwrl!pa.dec.com!shlump.nac.dec.com!ryn.mro4.dec.com!allvax.enet.dec.com!jroth
From: jroth@allvax.enet.dec.com (Jim Roth)
Newsgroups: comp.theory.dynamic-sys,sci.math,sci.physics
Subject: Re: square roots by hand or computer (or by spigot)
Message-ID: <4141@ryn.mro4.dec.com>
Date: 15 Mar 91 08:09:17 GMT
Sender: guest@ryn.mro4.dec.com
Followup-To: comp.theory.dynamic-sys
Organization: Digital Equipment Corporation
Lines: 85


In article <1991Mar14.113902.97@bsu-ucs.uucp>, 00lhramer@bsu-ucs.uucp (Leslie Ramer) writes...
>In article <1991Mar1.001103.6341@cec1.wustl.edu>, amc4919@cec_ss1.wustl.edu (Adam M. Costello) writes:
>> ... pencil and paper method
> ... newtons method

If you just want to see a decimal expansion to many significant figures
a simple way to do it is with a spigot algorithm.

Approximate the root of n with a convenient fraction p/q such that
n = p^2/q^2*(1+s/p^2), where s/p^2 is small.

For example sqrt(5) is approximately 9/4 and 5 = 81/16*(1-1/81).

Then the root can be approximated by a binomial expansion

    sqrt(n) = p/q (1 + s/p^2)^(1/2)

            = 1/q*(p + s/p^2*(p - s/(2*p^2)*(p - 3*s/(4*p^2)*(p - ...

            = 1/q*(a[0] + b[1]/c[1]*(a[1] - b[2]/c[2]*(a[2] - b[3]/c[3]*(a[3]...

where the a[k]'s are initially set to p.

Now, multiply all a[k]'s by your radix (a power of 10) and reduce them
modulo the denominators, c[k], passing the quotients times the numerators
b[k] to the previous a[k]'s.

At the end of the series a new digit will pop out.

Here is a simple program for sqrt(5); for roots of other numbers
initialize p, q, and s accordingly.

Practical?  Of course not!  But it does generate lots of digits easily :-)

- Jim

/*
 *  Spigot algorithm for sqrt(n)
 */

#include <stdio.h>

main()
{
    long
	*a,
	nterms, ndigits,
	p, q, r, s, t, d, n,
	i, j, k, tprev;

    nterms = 210;	/* use this many terms in the binomial expansion */
    ndigits = 100;	/* generate 100*4 digits */
    r = 10000;		/* 4 digits at a time */

    a = (long *)malloc(sizeof(long)*(nterms+1));

    n = 5;
    p = 9;
    q = 4;
    s = -1;

    for (i = 0; i <= nterms; i++)
	a[i] = p;

    for (i = 0; i < ndigits; i++) {
	for (d = 0, k = nterms; k; k--) {
	    a[k] = a[k]*r - d;
	    t = a[k]/(2*k*p*p);
	    a[k] -= t*(2*k*p*p);
	    d = t*(2*k-3)*s;
	}
        a[0] = a[0]*r - d;
        t = a[0]/q;
        a[0] = a[0] - t*q;
	if (t < 0) { 
	    t += r;
	    tprev--;
	}
        if (i) printf("%8d", tprev);
	tprev = t;
        if (!(i % 8)) printf("\n");
    }
    printf("\n");
}