\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 2014: Homework 7}
\date{\vspace{-2ex} due 10am Tuesday May 20th}
\maketitle

{\em
Last one! Shorter one to let you focus on projects. The last question has a little more
coding, but I give you lots of clues, and encourage collaboration.}

\ben

%\item Code fast exp via strassen integers eg $2^{1639}$ first?

\item BBP algorithm for arbitrary binary digits of $\ln 2$.
\ben
\item Prove that $\ln 2 = \sum_{k=1}^\infty (1/k2^k)$.
[Hint: take $\frac{1}{1-x^{-1}} - 1 - \frac{1}{x}$, expand, and
integrate over $(a,+\infty)$, for some $a$, in two different ways.]

\item 
Write a Matlab (using standard doubles to represent integers),
or python (not arbitrary precision) function {\tt r = expmodk(b,n,k)}
that implements fast binary exponentiation (mod $k$) to compute
$r = b^n (\mbox{mod} k)$.
Include a driver that tests it on small numbers for which the answer
can be computed directly, including annoying cases like $n=0$, $k=1$, etc.

\item
Use this to code up the formula for $\{ 2^d \ln 2\}$ using the two sums.
You'll want Matlab's {\tt mod(x,1)} or python's {\tt x \% 1}
applied to {\em each} term in the 1st sum.
Include 50 or so terms in the exponentially convergent 2nd sum.
Output the fractional part in binary, e.g.\ using {\tt dec2bin} in Matlab.
[Hint: test your code by changing $d$ by e.g.\ 20 and checking
that the overlapping digits are the same.]

\item
Compute the 50 binary digits starting at the $10^7$th, for $\ln 2$,
and state how long your code took to run. [Hint: warm up on
smaller $d$ values, and make sure to check the last digits
using the above overlapping method.]
\een


\item
In the file {\tt threekeys.txt} are three
1024-bit RSA public keys (numbers of size around $2^{1024}$ of the form $N=qp$
with $p$, $q$ prime), in decimal. Your goal is to crack (factor) them. You'll want to work
in sage/python.
\ben
\item It turns out that two of them have a factor in common. Find it, and hence crack two
of the three, i.e. give the factors. (This is a toy version of the following:
by mining thousands of RSA keys for common factors
it was discovered in February 2012 that, due to faulty random number generator,
such factors are way more common than expected!\footnote{This required
Dan Bernstein's idea of {\em batch} searching.})
\item
Factor the remaining key by writing a short code for Fermat's method. How many steps
did you need? (This shows the danger of having $q-p$ not much bigger than $\sqrt{p}$.)
\een



\item
Build a python function for Kraitchik's algorithm from class; this is a baby
``Quadratic Sieve'' (without an actual sieve). This will be good practise handling python
{\em lists}, and satisfying when done.
Your interface should be {\tt kraitchik(N,y,r)} where
$N$ is the integer to factor, $y$ sets the maximum size of prime to include in
your factor base, and $r$ is the number of successive $x$ values to try.
Use your code to factor the numbers:
%(which are both way more expensive via trial division or Fermat!):
\ben
\item
1180591624032052314157 (factors are too far apart for Fermat, way too large
for trial division)
\item
$2^{2^6}+1 = 18446744073709551617$, the 6th Fermat number (although more specialized
methods exist for Fermat numbers\ldots)
\een
Tuning the parameters: increase $r$ until you get around as many $y$-smooth
candidates $x^2-N$ as the size of your factor base. Increase $y$ a bunch then repeat.
Eventually you'll start finding useful vectors in the kernel.
(a) and (b) shouldn't exceed 1 minute runtime (seek help if they do).

Hints:

\bi
\item I will let you use \verb?fb = prime_range(y)? to set the factor base
\item I will let you use the built-in {\tt f = factor(...)} to extract the small
prime factors of the ``small'' (compared to $N$) numbers $x^2-N$.
This is wasteful but saves you coding a sieve (or batch trial division).

\item To build a matrix over $\{0,1\}$ from a list of lists, eg, {\tt a = [[1,2],[3,4],[5,6]]},
use {\tt A = matrix(GF(2),3,2,a)}

\item To find the set of null-space row vectors use \verb?A.left_kernel().basis()?.
Each vector is a {\em tuple}.

\item Python tricks: elementwise eg multiplication of two lists {\tt x} and {\tt y}
done by: {\tt [a*b for a,b in zip(x,y)]}.
You will only need this if you decide to construct $v$ from its prime factors\footnote{Actually this is more efficient for huge cases---see Crandall--Pomerance book p.268 (6)---but you'll find it easier to get $v$ from $\sqrt{(x_1^2-N)(x_2^2-N)\cdots}$}.
Product of a list is done by {\tt prod}.
Add {\tt y} to a list {\tt x} via {\tt x.append(y)}. Ask if stuck.

\item
Here's a routine that does the messy job of converting the output of {\tt f=factor(...)}
into a list of exponents (0,1,2, etc) given {\tt fb} the factor base list:
\begin{verbatim}
def factorbaseexponents(fb,f):
    """convert factor tuple list f into a list of factor base exponents
    for the factor base list fb
    """
    fbex = [0]*len(fb)       # list of zeros
    for t in list(f):        # loop over tuples in factors
        i = fb.index(t[0])   # 1st el of tuple has to be in fb list since smooth
        fbex[i] = t[1]
    return fbex
\end{verbatim}

\item Debug your code on the worksheet example, printing out everything. Then debug
on\\
{\tt kraitchik(1098413,25,200)},
which should find the six $x$ values 1051, 1063, 1077, 1119, 1142, 1237 for which
$x^2-N$ is $25$-smooth, and
a kernel vector (0, 1, 1, 0, 0, 1) which leads to the factorization $563 \times 1951$.
The latter example comes from Brent's 2010 slides on the course website.
\ei


\een
\end{document}

%log 2 digits via BBP ?  % few sec for 1e6; 77 sec for 1e7.
