Lab 4

Get the code

Find the Lab 4 code in the class repository.

$ ls
README.md  classes/   labs/      scripts/
$ ls labs/lab4
README    cpp/      julia/    python/   rust/
$ ls labs/lab4/cpp
collect.bash*   pi.cpp          series.cpp      view-pi.gp
makefile        ratpow.cpp      view-conv.gp    view-series.gp

Be sure that your copy of the bash script is flagged as executable. (Meta data in the form of file flags is not always preserved when copying.) If necessary, you can use chmod +x to make the change:

$ ls -l collect.bash
-rw-------@ 1 kbeach  staff  174 Aug 27  2020 collect.bash
$ chmod +x collect.bash
$ ls -l collect.bash
-rwx--x--x@ 1 kbeach  staff  174 Aug 27  2020 collect.bash*

The instructions here focus on C++, but you may carry out the comparable steps in Julia, Python, or Rust.

Raising to rational powers

Earlier in these notes, we mentioned Heron’s iterative method for computing the square root of a number. (One of the most ancient numerical algorithms—it’s been known for nearly two thousands years.) In Exercise 2, you derived a generalization of this method that computes \(\mathsf{a^{p/q}}\) with \(\mathsf{p}\) and \(\mathsf{q}\) integer:

\[\begin{aligned} \mathsf{x_0} & := \mathsf{1}\\ \mathsf{x_{n+1}} & := \mathsf{\Bigl(1-\frac{1}{q}\Bigr)x_n} + \mathsf{\frac{a^p}{q}x_n^{1-q}} \end{aligned}\]

If the sequence \(\mathsf{(x_n)}\) converges, it must converge to \(\mathsf{a^{p/q}}\).

Fill in the body of the program ratpow.cpp so that the recurrence relation shown above is computed starting from \(\mathsf{x = 1}\). Arrange for values of \(\mathsf{a}\), \(\mathsf{p}\), and \(\mathsf{q}\) to be read in from the command line. The formally correct stopping condition is that \(\mathsf{(x_n)^q}\) becomes sufficiently close to \(\mathsf{a^p}\). Instead, let’s just have the loop terminate when the next value in the sequence is sufficiently close to the last one. That comparison should be carried out not with the == operator but with a function nearly_equal that checks whether two numbers \(\mathsf{x}\) and \(\mathsf{y}\) satisfy \(\mathsf{|x-y|} < \mathsf{\epsilon |x|}\) and \(\mathsf{|x-y|} < \mathsf{\epsilon |y|}\).

$ make ratpow
g++ -o ratpow ratpow.cpp -O2 -ansi -pedantic -Wall -lm
$ ./ratpow
Returns the rational power a^(p/q) and compares to the
value computed with std::pow in the cmath library
Usage: ratpow a p q
        a=nonnegative fp number
      p,q=nonzero integers
$ ./ratpow 1.0 1 1
Newton's method value: 1
 C Math library value: 1
$ ./ratpow 2.0 1 2
Newton's method value: 1.414213562373095
 C Math library value: 1.414213562373095
$ ./ratpow 3.5 4 5
Newton's method value: 2.724296895429098
 C Math library value: 2.724296895429098

Series expansion for atanh

Version 1

The file series.cpp defines a function atan_series that computes the first N terms of the series expansion

\[\begin{aligned} \mathsf{arctan}\,\mathsf{x} &= \mathsf{x} - \mathsf{\frac{1}{3}x^3} + \mathsf{\frac{1}{5}x^5} - \mathsf{\frac{1}{7}x^7} + \cdots\\ &= \mathsf{x}\biggl[ \biggl(\biggl(\biggl(\biggl(\cdots + \mathsf{\frac{1}{9}\biggr)x^2} - \mathsf{\frac{1}{7}\biggr)x^2} + \mathsf{\frac{1}{5}\biggr)x^2} - \mathsf{\frac{1}{3}\biggr)x^2} + \mathsf{1}\biggr]. \end{aligned}\]

The function makes use of Horner’s scheme (we encountered this in Lab 2) to minimize the number of arithmetic operations required. Change the code in main so that the output is in scientific notation and uses 10 decimal digits of precision. You should be able to reproduce exactly the session shown below.

$ make series
g++ -o series series.cpp -O2 -ansi -pedantic -Wall -lm
$ ./series
Calculates the sequence of partial series expansions

            x
            x - x^3/3
            x - x^3/3 + x^5/5
            x - x^3/3 + x^5/5 - x^7/7
            etc.

approximating atan(x) = x - x^3/3 + x^5/5 - x^7/7 + ...

Usage: series x-value
$ ./series 0.5
For x = 0.5, the successive partial expansions have the following values:
     1    4.5833333333e-01
     2    4.6458333333e-01
     3    4.6346726190e-01
     4    4.6368427579e-01
     5    4.6363988659e-01
     .
     .
     .
    28    4.6364760900e-01
    29    4.6364760900e-01
    30    4.6364760900e-01
--------------------------
   inf    4.6364760900e-01
$ gnuplot -persist view-series.gp

The view-series.gp script plots the results for x = 0.5, 0.75, 1.0, 1.05 as a function of the number of terms included. The radius of convergence for this expansion is 1.

Alternate series for \(\textsf{arctan}\,\mathsf{x}\) with better convergence properties can be constructed by expanding with respect to new variable \(\mathsf{y = x^2/(1+x^2)}\) that maps \(\mathsf{x} \in \mathbb{R}\) onto \(\mathsf{y(x)} \in \mathsf{[0,1)}\). Keep in mind that \(\mathsf{y}\) itself has an expansion \(\mathsf{y} = \mathsf{x^2} - \mathsf{x^4} + \mathsf{x^6} - \mathsf{x^8} + \cdots\), so a series truncated at finite order in \(\mathsf{y}\) actually includes terms at all orders in \(\mathsf{x}\).

Version 2

Since the inverse of \(\mathsf{y(x)}\) is

\[\mathsf{x(y)} = \pm \mathsf{\frac{\sqrt{(1-y)y}}{1-y}},\]

we can expand \(\mathsf{(1/x)\mathsf{arctan}(x)}\) as

\[\mathsf{\frac{1-y}{\sqrt{(1-y)y}}} \textsf{arctan}(\mathsf{x(y)}) = \mathsf{1} - \mathsf{\frac{1}{3}y - \frac{2}{3\cdot 5}y^2} - \mathsf{\frac{2\cdot 4}{3\cdot 5\cdot 7}y^3} - \mathsf{\frac{2\cdot 4\cdot 6} {3\cdot 5\cdot 7\cdot 9}y^4} - \cdots,\]

which leads to

\[\begin{aligned} \textsf{arctan}(\mathsf{x}) &= \mathsf{x}\biggl(\cdots + \mathsf{\frac{6\cdot 4\cdot 2\cdot(-1)} {9\cdot 7\cdot 5\cdot 3}y^4} + \mathsf{\frac{4\cdot 2\cdot(-1)}{7\cdot 5\cdot 3}y^3} + \mathsf{\frac{2\cdot(-1)}{5\cdot 3}y^2} + \mathsf{\frac{(-1)}{3}y} + \mathsf{1} \biggr)\\ &= \mathsf{x}\biggl[ \biggl(\biggl(\biggl(\biggl( \cdots + \mathsf{1\biggr)\frac{6}{9}y} + \mathsf{1\biggr)\frac{4}{7}y} + \mathsf{1\biggr)\frac{2}{5}y} + \mathsf{1\biggr)\frac{(-1)}{3}y} + \mathsf{1} \biggr]. \end{aligned}\]

Version 3

Similarly, \((1/x+x)\arctan(x)\) has an expansion

\[\mathsf{\frac{x(y)}{y}}\textsf{arctan}(\mathsf{x(y)}) = \mathsf{1} + \mathsf{\frac{2}{3}y} + \mathsf{\frac{2\cdot 4}{3 \cdot 5}y^2} + \mathsf{\frac{2\cdot 4\cdot 6}{3\cdot 5\cdot 7}y^3} + \mathsf{\frac{2\cdot 4\cdot 6\cdot 8}{3\cdot 5\cdot 7\cdot 9}y^4} + \cdots\]

Rearranging the terms gives

\[\begin{aligned} \textsf{arctan}(\mathsf{x}) &= \mathsf{\frac{y}{x}}\biggl( \cdots + \mathsf{\frac{2\cdot 4\cdot 6\cdot 8}{3\cdot 5\cdot 7\cdot 9}y^4} + \mathsf{\frac{2\cdot 4\cdot 6}{3\cdot 5\cdot 7}y^3} + \mathsf{\frac{2\cdot 4}{3 \cdot 5}y^2} + \mathsf{\frac{2}{3}y} + \mathsf{1} \biggr)\\ &= \mathsf{\frac{y}{x}}\biggl[ \biggl(\biggl(\biggl(\biggl( \cdots + \mathsf{1\biggr)\frac{8}{9}y} + \mathsf{1\biggr)\frac{6}{7}y} + \mathsf{1\biggr)\frac{4}{5}y} + \mathsf{1\biggr)\frac{2}{3}y} + \mathsf{1} \biggr]. \end{aligned} \]

With these results in mind, make the following changes to series.cpp. Implement versions 2 and 3 of the expansion in the functions named atan_series2 and atan_series3. Have the program output each version in its own column.

$ ./series 0.99
For x = 0.99, the successive partial expansions have the following values:
     1    6.6656700000e-01    9.9000000000e-01    4.9997474875e-01
     2    8.5676500998e-01    8.2665824958e-01    6.6495808335e-01
     3    7.2361281742e-01    7.9431821634e-01    7.3028818435e-01
     4    8.2511473381e-01    7.8517106891e-01    7.5800541305e-01
     5    7.4372034706e-01    7.8215266270e-01    7.7020037778e-01
     .
     .
     .
    28    7.8518321535e-01    7.8037308008e-01    7.8037307961e-01
    29    7.7581569879e-01    7.8037308007e-01    7.8037307985e-01
    30    7.8469578200e-01    7.8037308007e-01    7.8037307996e-01
--------------------------
   inf    7.8037308007e-01
$ gnuplot -persist view-conv.gp

The view-conv.gp script shows the rate of convergence graphically. Convince yourself that version 2 of the expansion converges the fastest.

Elliptic integral

The complete elliptic integral of the first kind,

\[\mathsf{I(a,b)} = \mathsf{\frac{1}{2}} \int_{-\infty}^{\infty} \mathsf{\frac{dt}{\sqrt{(a^2+t^2)(b^2+t^2)}}},\]

has the special property that

\[\mathsf{I(a,b)} = \mathsf{I\biggl( \frac{a+b}{2},\sqrt{ab}\biggr)}.\]

Amazingly, this simple result has rather deep implications (related to the so-called arithmetic-geometric mean) and can be used to derive [Borwein and Borwein, SIAM Review 26, 351 (1984)] an iterative scheme for computing \(\pi\) that converges exponentially: viz., one starts with the initial conditions \(\alpha_\mathsf{0} = \sqrt{\mathsf{2}}, \ \beta_0 = \mathsf{0}, \ \pi_\mathsf{0} = \mathsf{2} + \sqrt{\mathsf{2}}\) and generates successive approximations recursively via

\[\begin{aligned} \alpha_{\mathsf{n+1}} &= \mathsf{\frac{1}{2}}\bigl(\mathsf{\alpha_n^{1/2} + \alpha_n^{-1/2}}\bigr),\\ \mathsf{\beta_{n+1}} &= \mathsf{\alpha_n^{1/2}}\biggl(\mathsf{\frac{1+\beta_n}{\alpha_n+\beta_n}}\biggr),\\ \pi_{\mathsf{n+1}} &= \mathsf{\pi_n\beta_{n+1}}\biggl(\mathsf{\frac{1+\alpha_{n+1}}{1+\beta_{n+1}}}\biggr). \end{aligned}\]

There are also several exact identities relating linear combinations of \(\textsf{arctan}\) to \(\pi\).

\[\begin{aligned} \pi &= \mathsf{4}\,\textsf{arctan}(\mathsf{1/2}) + \mathsf{4}\,\textsf{arctan}(\mathsf{1/5}) + \mathsf{4}\textsf{arctan}(\mathsf{1/8})\\ &= \mathsf{16}\,\textsf{arctan}(\mathsf{1/5}) - \mathsf{4}\,\textsf{arctan}(\mathsf{1/239})\\ &= \mathsf{24}\,\textsf{arctan}(\mathsf{1/8}) + \mathsf{8}\,\textsf{arctan}(\mathsf{1/57}) + \mathsf{4}\,\textsf{arctan}(\mathsf{1/239}) \end{aligned}\]

Modify the file pi.cpp so that the program outputs in the first column an integer count, in the second column the AGM expansion result, and in the three other columns the numerical estimates of \(\pi\) based on the arctan formulas shown above. The first column refers either to the number of iterations in the AGM scheme or to the number of terms included in the series expansion of arctan. Be sure to use the function atan_series2 that you wrote previously. The final output should look like this:

$ make pi
g++ -o pi pi.cpp -O2 -ansi -pedantic -Wall -lm
$ ./pi
Usage: pi #iterations
$ ./pi 7
  1     3.1426067539416                 3.3     3.1832635983264     3.1570872788666
  2      3.141592660966     3.1538461538462     3.1422380549654     3.1416881708556
  3     3.1415926535898     3.1430059171598      3.141606891258      3.141593494527
  4     3.1415926535898      3.141783262891     3.1415930195281     3.1415926622229
  5     3.1415926535898     3.1416206328643     3.1415926638427     3.1415926536865
  6     3.1415926535898     3.1415969882255     3.1415926538935     3.1415926535909
  7     3.1415926535898        3.1415933509     3.1415926535991     3.1415926535898
-----------------------------------------------------------------------------------
inf     3.1415926535898
$ gnuplot -persist view-pi.gp

If you’ve done everything correctly, the script view-pi.gp should give you a plot of the discrepancy between the approximation and exact value in each of the four cases. Be sure you understand what is implied by the shapes of the various curves.