[VA] SRC #009 - Pi Day 2021 Special

[robve]

(03-18-2021 02:44 PM)Albert Chan Wrote:  Thanks for spending the time to add code challenges.
Both codes are calculating pi by 2/(sin(x)/x) at x=pi/2, in 2 different ways.

You're very welcome!

Because this is going to be a long post, I will post the two parts a and b separately.

a.

10 P=SQR(2)
20 Q=P/2
30 DISP 2/Q
40 P=SQR(2+P)
50 Q=Q*P/2
60 IF P<2 GOTO 30

As Albert replied correctly, the code computes Viète's formula, published in 1593:

$$ \frac{2}{\pi} = \frac{\sqrt{2}}{2} \cdot \frac{\sqrt{2+\sqrt{2}}}{2} \cdot \frac{\sqrt{2+\sqrt{2+\sqrt{2}}}}{2} \cdots $$

In recurrence form:

$$ \lim_{n\rightarrow\infty} \prod_{i=1}^n \frac{a_i}{2} = \frac{2}{\pi};\qquad a_1 = \sqrt{2},\quad a_n = \sqrt{2+a_{n-1}} $$

In functional form:

$$ \pi \approx \mathrm{viete}(n) = 2/\prod_{i=1}^n v(n)/2;\qquad v(n) = \begin{cases} \sqrt{2} & n=1 \\ \sqrt{2+v(n-1)} & n>1 \end{cases} $$

Viète obtained his formula by comparing the areas of regular polygons with \( 2^n \) and \( 2^{n+1} \) sides inscribed in a circle. The first term in the product, \( \frac{\sqrt{2}}{2} \), is the ratio of areas of a square and an octagon, the second term is the ratio of areas of an octagon and a hexadecagon, etc. - Wikipedia

See also "Mathematics: From the Birth of Numbers" by Jan Gullberg.

This directly relates to Archimedes's famous work (ca.225BC) on approximating the area of a circle by polygons inside and outside the circle squeezing the circle: "At each stage, he needed to approximate sophisticated square roots, yet he never faltered. When he reached the 96-gon, his estimate was \( \frac{6336}{2017\frac{1}{4}} > 3\frac{10}{71} \) ." - "Journey Through Genius" by William Dunham.

Note that Euclid's Elements does not prove the ratio of the radius of the circle to its circumference \( 2\pi \). Sometimes confused with Euclid VI.33 that proves an important property of angles and arcs but does not compare them to circles with different radii "in equal circles [emphasis mine], angles, whether at the center or circumference, are in the same ratio to one another as the arcs on which they stand." - Oliver Byrne's Euclid 1847 VI.33

   
[image credit: Oliver Byrne's Euclid 1847]

Viète's formula can also be written as

$$ \frac{\pi}{2} = \frac{1}{\sqrt{\frac{1}{2}}} \cdot \frac{1}{\sqrt{\frac{1}{2}+\frac{1}{2}\sqrt{\frac{1}{2}}}} \cdot \frac{1}{\sqrt{\frac{1}{2}+\frac{1}{2}\sqrt{\frac{1}{2}+\frac{1}{2}\sqrt{\frac{1​}{2}}}}} \cdots $$

In recurrence form:

$$ \lim_{n\rightarrow\infty} \prod_{i=1}^n \frac{1}{a_i} = \frac{\pi}{2};\qquad a_1 = \sqrt{\frac{1}{2}},\quad a_n = \sqrt{\frac{1+a_{n-1}}{2}} $$

In functional form:

$$ \pi \approx \mathrm{viete}(n) = 2/\prod_{i=1}^n v(n);\qquad v(n) = \begin{cases} \sqrt{1/2} & n=1 \\ \sqrt{(1+v(n-1))/2} & n>1 \end{cases} $$

While not based on Viète's formula, with our calculators we can quickly integrate the unit quarter circle defined by \( 1=x^2+y^2 \), i.e. integrate \( y = f(x) = \sqrt{1-x^2} \) to get \( \pi \), using a square root:

$$ \pi = 4\int_0^1 \sqrt{1-x^2} \,dx $$

Proof (I've simplified this somewhat):

$$ 4 \int_0^1 \sqrt{1-x^2} \,dx = 4 \int_0^{\frac{\pi}{2}} \sqrt{1-\sin^2 \theta} \cos \theta \,d\theta = 4 \int_0^{\frac{\pi}{2}} \sqrt{\cos^2 \theta} \cos \theta \,d\theta = 4 \int_0^{\frac{\pi}{2}} \cos^2 \theta \,d\theta = \pi $$

Programs I wrote today to illustrate Viete's formula

My own functional programming language Husky:

v(1) := (r,r) where r := sqrt(2)/2;
v(n) := (t,q) where t := p*q
              where q := sqrt(2+2*r)/2
              where (p,r) := v(n-1).
viete(n) := 2/p where (p,r) := v(n).
> viete(20).

and a list-based version:

prod := foldr(*, 1).
v(1) := [sqrt(2)/2];
v(n) := [sqrt(2+2*x)/2, x | xs] where x.xs := v(n-1).
viete(n) := 2/prod(v(n)).
> viete(20).

Haskell:

v 1 = (r,r) where r = (sqrt 2)/2
v n = (t,q) where t = p*q
                  q = (sqrt (2+2*r))/2
                  (p,r) = v (n-1)
viete n = 2/r where (r,_) = v n
main = putStrLn (show (viete 20))

and a list-based version:

prod = foldr (*) 1
viete n = 2/(prod (v n))
v 1 = [(sqrt 2)/2]
v n = (sqrt (2+2*x))/2 : x : xs where x:xs = v (n-1)
main = putStrLn (show (viete 20))

Prolog:

viete(N, P) :- v(N, R, _), P is 2/R.
v(1, R, R) :- R is sqrt(2)/2, !.
v(N, T, Q) :- M is N-1, v(M, P, R), Q is sqrt(2+2*R)/2, T is P*Q.
?- viete(20, P).

C (version with convergence check)

#include <stdio.h>
#include <math.h>
int main() {
  double p = sqrt(2), q = p/2;
  while (p < 2)
    q *= (p = sqrt(2+p))/2;
  printf("%.17g\n", 2/q);
}

My own MiniC C-like language:

int main() {
  float p, q;
  p = sqrt(2.0);
  q = p/2;
  while (p < 2.0)
    q *= (p = sqrt(2+p))/2;
  print 2/q;
}
$ ./minic viete.c
$ java viete

Java (version with convergence check):

import java.lang.*;
public class Viete
{
  public static void main(String[] arg)
  {
    double p = Math.sqrt(2), q = p/2;
    while (p < 2)
      q *= (p = Math.sqrt(2+p))/2;
    System.out.println(2/q);
  }
}

Python (version with convergence check):

from math import sqrt
def viete():
   p = sqrt(2)
   q = p/2
   while p < 2:
      p = sqrt(2+p)
      q *= p/2
   print(2/q)
if __name__ == "__main__":
   viete()

HPPL (version with convergence check):

EXPORT viete()
BEGIN
  LOCAL p = √2, q = p/2;
  REPEAT
    p := √(2+p);
    q := q*p/2;
  UNTIL p >= 2;
  RETURN 2/q;
END;

HP-71B FORTH with MultiMod (forthcoming: where can I find the pdf manual?)

Edit: Thanks to rprosperi's help to locate the missing FORTH manual, here is my HP-71B FORTH program:

FVARIABLE P
FVARIABLE Q
: VIETE
  2. SQRT P STO
  2. F/ Q STO
  BEGIN
    2. P RCL F+ SQRT P STO
    2. F/ Q RCL F* Q STO
    2. P RCL X>=Y?
  UNTIL
  2. Q RCL F/ F.
;
VIETE

HP-71B FORTH significantly extends ANS FORTH. Couldn't be happier with MultiMod installed on my HP-71B!

- Rob

Minor edit: fixed typo and added list-based functional versions of Husky and Haskell programs and HP-71B FORTH.