\documentclass[10pt,letterpaper]{article}
\usepackage{graphicx,amssymb}
\usepackage[colorlinks=true]{hyperref}

\topmargin -0.5in
\textheight 9in
\oddsidemargin=-0.05in
\evensidemargin=-0.05in
\textwidth 6.5in
\pagestyle{empty}

\newcommand{\bi}{\begin{itemize}}
\newcommand{\ei}{\end{itemize}}
\newcommand\ben{\begin{enumerate}}
\newcommand\een{\end{enumerate}}
\newcommand\vg{\vspace{2ex}}
\newcommand{\bc}{\begin{center}}
\newcommand{\ec}{\end{center}}
\newcommand{\ie}{{\it i.e.\ }}
\newcommand{\eg}{{\it e.g.\ }}
\newcommand{\eps}{\varepsilon}
\newcommand{\mbf}[1]{{\mathbf #1}}
\newcommand{\tbox}[1]{{\mbox{\tiny #1}}}


\begin{document}
\title{Math 56 Compu \& Expt Math, Spring 2013: Homework 6}
\date{\vspace{-2ex} due 10am Thursday May 9th}
\maketitle

{\em
Options for python/SAGE setup:
1)
install python and the {\tt mpmath} and {\tt gmpy} packages, or
2) install SAGE, or
3) work ``in the cloud''
by creating a notebook account at {\tt sage.dartmouth.edu}

This lets you do arbitrary precision; for now you use only the four basic
arithmetic operations
(although of course they have many more available).

For {\tt mpmath} usage see: \verb?http://mpmath.googlecode.com/svn/trunk/doc/build/basics.html?
}

\ben

\item Cementing the concept of super-algebraic convergence.
\ben
\item Write down a function $f(n)$ exhibiting super-algebraic convergence
to zero as $n\to\infty$,
but such that $e^{-c\sqrt{n}} = o(f(n))$ for each $c>0$.

\item Fix $c>0$. Let $k$ be an algebraic convergence order.
Find the smallest $k$-dependent ``constant'' $C_k$ you can, such that,
for each order $k=1,2,\ldots$, we have
$e^{-c\sqrt{n}} \le C_k n^{-k}$ for all $n=1,2,\ldots$.
What is $\lim_{k\to\infty}C_k$ ?
This shows that the implied big-O constant\footnote{or this could
instead be $N_0$ as in $n\ge N_0$, although here $N_0=1$ hence the constant
has to take the hit.}
in the definition
of super-algebraic need not be uniformly bounded with respect to order.
[Hint: Fix $k$, put all upstairs in exp then maximize vs $n$.]
%Skip any proof that the maximum you find is the only one.]
\een

\item Fast multiplication, and I mean fast.
\ben
\item
Write a Matlab function for Strassen's fast multiplication of
arbitrary-precision integers in base 10,
with an accompanying driver/test code.
Test that the output is a vector of integers in the range 0 to 9.
You may use the built-in FFT.
[Hint: this is pretty quick if you adapt my code for addition of integers.
Make sure to do something to the output of the convolution so that the carrying
works reliably!]
\item
Use your code to compute the ridiculously big number
$2^{2^{22}}$ (note the top number is the decimal number
22). Taking 1 and doubling it $2^{22}$ times is impractical---instead devise a
much faster scheme.
How many decimal digits does this number have?!
Print the last 10 digits. Dial them into a phone.\footnote{only joking.}
How long does the calculation take? Estimate (roughly, ie
within a factor of two) how long it would take using
naive arbitrary-precision multiplication at $10^9$ flop/sec.
[Hint: if {\tt x} is a base-10
vector, \verb?fprintf('%d',x)? is a useful way to display]
\item [BONUS] Make conjectures about the last couple of decimal digits
of $2^{2^j}$ for $j=2,3,\ldots$ and prove them.
\een


\item Fast division.
\ben
\item
In python (or SAGE) make a function {\tt x = fastrecip(z,tol,xguess)}
which implements the Newton method for the reciprocal of {\tt z},
with initial guess {\tt xguess},
stopping at relative tolerance {\tt tol}.
This being python, you can have the test driver and the function in
the same file.
Note {\tt xguess} needs to be reasonably close to the answer (how near?)

\item
Import {\tt mpmath} and set working precision to 1000 decimal digits.
Compute $1/97$ using your above routine (don't use any division!).
What is the repetition period of its decimal expansion?
Hints:
\begin{verbatim}
from mpmath import *
z = mpf('97')                # the arb-precision representation of 97
\end{verbatim}

% https://en.wikipedia.org/wiki/Repeating_decimal


\item Assuming constant effort per flop%
\footnote{This is not strictly true, since as $N$ grows, more floating-point accuracy is needed in the FFT so that it rounds to the correct integers; see Crandall--Pomerance book Sec 9.5.} in the FFT,
deduce the complexity in big-O notation for
your fast reciprocal of an $N$-digit number.

\item
The wikipedia page
\verb?http://en.wikipedia.org/wiki/Division_algorithm?
in the section {\bf Newton--Raphson division}, in the paragraph starting
``From a computation point of view\ldots''
claims
that $x_{n+1} = x_n(2 - z x_n)$ is (much) less accurate
than the mathematically equivalent $x_{n+1} = x_n + x_n(1-z x_n)$.
Evaluate this claim by using your skills in forward error analysis
(rules of floating point). Ignore initial rounding, and you may assume
$x_n$ is near its final value.
\een


\item
Algorithms for digits of $\pi$. Please work in python/SAGE/mpmath/gmpy
and make use of all basic arbitrary-precision arithmetic and sqrt.
\ben
\item
Make a function {\tt y = atantaylor(x,n)} which sums the {\tt n}-term
Taylor series about the origin (i.e. Maclaurin series)
to approximate $y = \tan^{-1} x$,
at whatever the current working precision is.
Test against {\tt mpmath}'s {\tt atan} to check.
[Hint: make sure all constants are multi-precision, e.g. {\tt mpf('0')}]
\item
Use the above to evaluate any of the Machin--Euler type formulae for
$\pi$ to 10000 digits. Give the last 10 of those
digits, state your convergence
rate, how many terms $n$ you used, and how many seconds it took.

Python tips: {\tt print str(x)[-10:]} for printing, and the following
for timing some chunk of code,

{\tt import time

t = time.time()

... (do something) ...

print time.time()-t 'secs'}

\item Consider Gauss' arithmetic-geometric mean iteration
that lies at the heart of the Brent--Salamin algorithm for $\pi$:
\begin{eqnarray}
x_{x+1} & = & (x_n+y_n)/2 \nonumber \\
y_{n+1}& = & \sqrt{x_n y_n} \nonumber
\end{eqnarray}
Prove that it is quadratically convergent to something.
[Hints: bound $|x_{n+1} - y_{n+1}|$ in terms of its previous value.
Play with the difference of two squares $x_{n+1}^2 - y_{n+1}^2$.]

\item Code up Brent--Salamin using {\tt mpmath} for all arithmetic
including the square-root (to make your life easier).
[Hint: debug at a ``low'' precision of only 100 digits,
while printing $2x_n^2/\alpha_n$ each iteration.]

How many times faster is it for $N=10^4$ digits than in (b)?
Quote the 10 digits of $\pi$ starting at the millionth
(since the last few digits in your answer will be wrong, you'll need
to go a few digits beyond what you need).

\een


\een
\end{document}

%log 2 digits via BBP ?
