Lab 4
=====

Get the code
------------

Find the Lab 4 code in the class repository.

.. code-block:: console

   $ 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:

.. code-block:: console

   $ 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
:math:`\mathsf{a^{p/q}}` with :math:`\mathsf{p}` and :math:`\mathsf{q}`
integer:

.. math::

   \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 :math:`\mathsf{(x_n)}` converges, it must converge to :math:`\mathsf{a^{p/q}}`.

Fill in the body of the program ``ratpow.cpp`` so that the recurrence
relation shown above is computed starting from :math:`\mathsf{x = 1}`.
Arrange for values of :math:`\mathsf{a}`, :math:`\mathsf{p}`, and
:math:`\mathsf{q}` to be read in from the command line. The formally
correct stopping condition is that :math:`\mathsf{(x_n)^q}` becomes
sufficiently close to :math:`\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 :math:`\mathsf{x}` and :math:`\mathsf{y}` satisfy
:math:`\mathsf{|x-y|} < \mathsf{\epsilon |x|}` and :math:`\mathsf{|x-y|}
< \mathsf{\epsilon |y|}`.

.. code-block:: console

   $ 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 

.. math::

   \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.

.. code-block:: console

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

Version 2
~~~~~~~~~

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

.. math::

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

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

.. math::

   \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

.. math::

   \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, :math:`(1/x+x)\arctan(x)` has an expansion

.. math::

   \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

.. math::
   \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.

.. code-block:: console

   $ ./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*,

.. math::

   \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

.. math::

   \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 :math:`\pi` that converges exponentially: viz., one
starts with the initial conditions :math:`\alpha_\mathsf{0}
= \sqrt{\mathsf{2}}, \ \beta_0 = \mathsf{0}, \ \pi_\mathsf{0} = \mathsf{2}
+ \sqrt{\mathsf{2}}` and generates successive approximations recursively
via

.. math::

   \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 :math:`\textsf{arctan}` to :math:`\pi`.

.. math::

   \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 :math:`\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:

.. code-block:: console

   $ 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.