Re: Pi with bad calculator

In sci.physics, Jim Black
<tramspap@xxxxxxxxx>
wrote
on 18 Mar 2007 21:58:45 -0700
<1174280325.623898.235740@xxxxxxxxxxxxxxxxxxxxxxxxxxx>:
On Mar 18, 7:45 pm, BioFreak <BioFr...@xxxxxxxxxxxxxxx> wrote:
While vacationing in the Ozarks you're stuck with a
little piece of grocery store calculator whose
extravagance doesn't go beyond a mere square root
button. But you're doing a back of the envelope
calculation of something that begs for a value for Pi
way better than 3.1415 which you can remember for sure.
Other than cursing yourself for not remembering the
next few digits what else can you do to manage?

Well, there are nice algorithms for doing this, but the only one I
remember off the top of my head is

pi = 4/1 - 4/3 + 4/5 - 4/7 + ...

and that converges way to slow. And if I could look up an algorithm,
I could look up pi.

So it's back to first principles. We can use the old trick of
doubling the number of sides in a regular polygon repeatedly to
approximate a circle. So in my head I draw a little diagram
containing half of a side of an n-gon, and a side of a 2n-gon, and
lines connecting everything to the center. I set the radius of the
circle to one, and I'll call the length of the half-side of the n-gon
y. By two applications of the Pythagorean theorem, the length of the
2n-gon side is:

sqrt( y^2 + (1 - sqrt(1-y^2))^2 )

Simplifying, we get:

sqrt( 2 - 2 sqrt(1-y^2) )

So what we need to iterate is

s' = sqrt( 2 - 2 sqrt(1 - s^2/4) ),

starting from s=1 in the case of the hexagon. But notice that we
would be taking a square root and then immediately squaring it, which
is clearly wasteful. So let's simplify the operations we have to do
between two of the square roots that can't be eliminated:

1 - (sqrt(2-2x))^2 / 4

= (x + 1)/2

That isn't so painful, after all. First, we do the calculation up to
the square root from s=1, getting

sqrt(3/4).

Then we iterate

x' = sqrt((x + 1)/2)

N times, where N is some large number. Then we finish off the
calculation of the side length with

s = sqrt(2 - 2x).

Last, we multiply by 3 * 2^(N+1) to get pi.

I tried this for N = 20 on my PC's simulacrum of a stupid calculator.

I got: 3.1415926535897605995224090804006
pi = 3.1415926535897932384626433832795

On a real stupid calculator, though, I'd be a bit worried about loss
of precision. Notice that at the end of the calculation,

sqrt(2 - 2x) ~ pi / 3 / 2^(N+1)

so

x ~ 1 - pi^2 / 9 / 2^(2N+3).

That's going to limit how far you can take the calculation out, unless
you want to spend even more time thinking and come up with a clever
way to get around the problem.

--
Jim E. Black


This isn't too bad. I might approach this from a slightly
different perspective, by noting that the perimeter is
N times the side, for some value of N, then taking that
chord as an abstract value, and computing the perimeter directly.

So let's assume s is the side of a regular polygon inscribed in the
circle, and p = N * s is our starting estimate for 2*pi. Therefore,
s/2 is the side of a right triangle whose hypotenuse is 1 (the radius
of a unit circle); the other side is therefore sqrt(1 - s^2/4).
Another right triangle has 1 - sqrt(1 - s^2/4) as one side (this is
the radius between chord and the unit circle), s/2 as the other side,
and the new hypotenuse s'; this new hypotenuse is of course

s' = sqrt((1 - sqrt(1 - s^2/4))^2 + s^2/4)
= sqrt((1 + 1 - s^2/4 - 2*sqrt(1 - s^2/4)) + s^2/4)
= sqrt(2 - 2*sqrt(1 - s^2/4))

Since p = N * s and p' = 2 * N * s', we can multiply both sides
by 2 * N and get:

p' = 2 * N * sqrt(2 - 2*sqrt(1 - s^2/4))
= 2 * sqrt(2*N^2 - 2*N^2*sqrt(1 - s^2/4))
= 2 * sqrt(2*N^2 - N*sqrt(4 * N^2 - N^2*s^2))
= 2 * sqrt(2*N^2 - N*sqrt(4 * N^2 - p^2))

Not the best of expressions -- I was hoping to simplify
N out entirely -- but one can start with N = 6, p = 6
(a regular hexagon) and hope for the best. I get:

N p pp(N,p)
6 6.000 6.211657082460498296373572103
12 6.211657082460498296373572103 6.265257226562476394323498940
24 6.265257226562476394323498940 6.278700406093734414270293643
48 6.278700406093734414270293643 6.282063901781019276222705853
96 6.282063901781019276222705853 6.282904944570924150901218618
192 6.282904944570924150901218618 6.283115215823715291032926690
384 6.283115215823715291032926690 6.283167784296636817337939207
....
6.283185307179586476925286767

I can't say I'm thrilled with this but at least it's
converging about as fast, or faster than, 1/N. The main
problem with my variant is the outer subtraction, which
will lose precision for large N.

If one takes q = (4*N^2 - p^2), q' = (16*N^2 - p'^2), then

q' = 16*N^2 - (4 * (2*N^2 - N*sqrt(q)))
= 8*N^2 + 4*N*sqrt(q)

This expression will not lose precision because of subtraction, except
at the final step where we calculate p = sqrt(4*N^2 - q). However,
q starts out rather large -- for N = 6, q = 108 -- and gets larger.
A non-exponential calculator might have problems. :-)

--
#191, ewill3@xxxxxxxxxxxxx
Useless C++ Programming Idea #10239993:
char * f(char *p) {char *q = malloc(strlen(p)); strcpy(q,p); return q; }

--
Posted via a free Usenet account from http://www.teranews.com

.

Quantcast