(import "java.math.BigInteger")
(import "java.net.URL")
(import "java.io.BufferedReader")
(import "java.io.InputStreamReader")
#{
This is about solving the google puzzle, posted on a billboard
http://www.google.com/googleblog/2004/07/warning-we-brake-for-number-theory.html

The billboard says:

 {first 10-digit prime found in consecutive digits of e}.com

When i first saw the billboard from the url above i could not read the
"of e" so i didn't under stand the puzzle until the other day.  By
now, you can solve this problem, and the next one, easily by googling.

However, i was interested in learning something from the problem.

By googling, you can find 5 million digits of e.  Here's a procedure that returns a generator that produces one digit of e each time it is called.

}#
(define  (echar)
  ;; Returns a thunk that returns a digit in the first million digits
  ;; of e each time it is called.d
  (let ((r (BufferedReader.
            (InputStreamReader.
             (.getContent
              (.openConnection
               (URL.
                "http://antwrp.gsfc.nasa.gov/htmltest/gifcity/e.5mil")))))))
    (let loop0 ((line (.readLine r)))
      (if (not (equal? line "e = ")) (loop0 (.readLine r))
          (begin (.readLine r)
                 (lambda ()
                   (let loop1 ((it (.read r)))
                     (if (= it -1) #f
                         (let ((it (integer->char it)))
                           (if (char-numeric? it)
                               (- (char->integer it) (char->integer #\0))
                               (loop1 (.read r))))))))))))
{

It will be good for testing to be able to compare this generator
against another.

}
(define (etest n generator)
  ;; Test that generator produces digits of e.
  (let ((g1 (echar))
        (g2 generator))
    (let loop ((i 0))
      (if (< i n)
          (let ((c1 (g1))
                (c2 (g2)))
            (if (= c1 c2)
                (begin
                  (print `(,i ,c1))
                  (loop (+ i 1)))
                (print `(,i ,c1 ,c2))))))))

{
I used http://www.hulver.com/scoop/story/2004/7/22/153549/352 who references
http://web.comlab.ox.ac.uk/oucl/work/jeremy.gibbons/publications/spigot.pdf
which is about "spigot algorithms" that let you generate one digit of a number
at a time.

Here we use a computation of e-1:
e-1 = 1 + 1/2*(1 + 1/3* (1 + 1/4* ...))

Notice that this can be written as a composition of function producer (gen):
}
(define (compose f g) (lambda (x) (f (g x))))

(define (gen n) (lambda (x) (+ 1.0 (/ (+ x 0.0) n))))

{
So we can write an expression for e as:
   (+ 1 (compose (compose (compose (gen 2) (gen 3))
                          (gen 4))
                 (gen 5)) x)
x is the remainder of the computation.

What's interesting is that x is between 1 and 2 so you can do branch
and bound to determine an accurate value of e.
}

(define (e)
  ;; Compute e to double precision.
  ;; This gives one digit less than (exp 1). 
  (let loop ((i 2)
             (g (gen 2)))
    (let ((low (g 1.0))
          (high (g 2.0)))
      (if (= high low) (list (+ low 1.0) i)
          (loop (+ i 1) (compose g (gen (+ i 1))))))))

{
A "spigot algorithm" generates one digit at a time.

Say g is the current composition of functions, as above.  When low and
high have the same integer value, say d, we can output d 
make a new g by (compose (lambda (x) (* (- x d) 10)) g).

}
(define (egen)
  ;; Generate digits of e (good for 16 digits).
  (let ((firstTime? #t)
        (g (gen 2))
        (n 2))
    (define (grow)
      (let ((low (g 1.0))
            (hi  (g 2.0)))
        (if (= (inexact->exact low) (inexact->exact hi))
            (let* ((d (inexact->exact low)))
              (set! g (compose (lambda (x) (* (- x d) 10.0)) g))
              (if firstTime? (begin (set! firstTime? #f)
                                    (+ d 1))
                  d))
            (begin 
              (set! n (+ n 1))
              (set! g (compose g (gen n)))
              (grow)))))))

{

The interesting thing is we can replace function composition by matrix
multiplication, using a "linear fractional transform", lft.  An lft is
a function of the form:

f(x) =(a*x + b)/(c*x + d)

If we represent this function as the matrix:

a b
c d

function composition become matrix multiplication.

We need matrices of the form

1 n  for the function (lambda (x) (+ 1 (/ x n)))
0 n

and 

10 -10*d for the function (lambda (x) (* (- x d) 10))
 0     1

Using long integers with this approach you can generate 25 digits of e
correctly.  To do more, use BigIntegers.

We'll represent the matrix as (vector a b c d) and the rational number
as (vector numerator denomintor), both with BigIntegers.
}

(define zero (BigInteger. "0"))
(define one (BigInteger. "1"))

(define-method (big (n BigInteger)) n)
(define-method (big (n String)) (BigInteger. n))
(define-method (big n) (big (.toString n)))

(define (b+ a b) (.add (big a) (big b)))
(define (b* a b) (.multiply (big a) (big b)))
(define (b/ a b) (.divide (big a) (big b)))
(define (bgcd a b) (.gcd (big a) (big b)))
(define (b- a . b) (if (pair? b) (.subtract (big a) (big b))
                       (.subtract zero (big a))))

(define (bvector . xs) (apply vector (map big xs)))

(define vr vector-ref)

(define (dot a b c d) (b+ (b* a b) (b* c d)))

(define (mmul a b)
  ;; matrix - matrix multiplication.
  ;; 12 operations.
  (msimplify (bvector (dot (vr a 0) (vr b 0) (vr a 1) (vr b 2))
                     (dot (vr a 0) (vr b 1) (vr a 1) (vr b 3))
                     (dot (vr a 2) (vr b 0) (vr a 3) (vr b 2))
                     (dot (vr a 2) (vr b 1) (vr a 3) (vr b 3)))))
  
(define (mmul a b)
  ;; matrix - matrix multiplication.
  ;; This version uses the fact that a[2] and b[2] are 0.
  ;; 5 operations
  (msimplify (bvector (b* (vr a 0) (vr b 0))
                     (dot (vr a 0) (vr b 1) (vr a 1) (vr b 3))
                     zero
                     (b* (vr a 3) (vr b 3)))))

(define (msimplify m)
  ;; Remove any common factor.
  (let ((it (bgcd (bgcd (vr m 0) (vr m 1))
                  (bgcd (vr m 2) (vr m 3)))))
    (if (.equals it one) m
        (bvector (b/ (vr m 0) it)
                 (b/ (vr m 1) it)
                 (b/ (vr m 2) it)
                 (b/ (vr m 3) it)))))

(define (meval m x)
  ;; Returns a rational number as #(numerator denominator).
  (vsimplify
             (vector (b+ (b* (vr m 0) x) (vr m 1))
                     (b+ (b* (vr m 2) x) (vr m 3)))))

(define (intpart x) (b/ (vr x 0) (vr x 1)))

(define (vsimplify v)
  (let ((it (bgcd (vr v 0) (vr v 1))))
    (if (= it one) v
        (vector (b/ (vr v 0) it) (b/ (vr v 1) it)))))

(define (genm n)(bvector 1 n 0 n))
(define (dig d) (bvector 10 (b* -10 d) 0 1))

(define (egen)
  (let ((firstTime? #t)
        (m (genm 2))
        (n 2))
    (define (grow)
      (let ((low (meval m 1))
            (hi  (meval m 2)))
        (print (list n m))
        (if (= (intpart low) (intpart hi))
            (let* ((d (intpart low)))
              (set! m (mmul (dig d) m))
              (if firstTime? (begin (set! firstTime? #f)
                                    (b+ d 1))
                  d))
            (begin 
              (set! n (b+ n 1))
              (set! m (mmul m (genm n)))
              (grow)))))))


#{
main(){int N=9009,n=N,a[9009],x;while(--n)a[n]=1+1/n;
for(;N>9;printf("%d",x))
for(n=N--;--n;a[n]=x%n,x=10*a[n-1]+x/n);}
}#

;;; http://www.7427466391.com/
;;;http://www.google.com/labjobs/index.html
;;; http://antwrp.gsfc.nasa.gov/htmltest/gifcity/e.5mil
;;; http://numbers.computation.free.fr/Constants/TinyPrograms/tinycodes.html