;;;;;;;;;
;;;
;;; Exercise 1.
;;;
;;; An O(log n)-time (counting multiplication as a single step) and
;;; O(1)-space algorithm for computing Fibonacci numbers
;;;
;;;;;;;;;

(load-option 'format)

;;; O(log n)-time and O(1)-space implementation of exponentiation
(define (expt x n)
  (define (aux x n res)
;   (format #t "expt: aux: n = ~A\n" n)
    (cond ((= n 0) res)
          ((= n 1) (* x res))
          (else (let ((qr (integer-divide n 2)))
                    (if (eq? (cdr qr) 0)
                        (aux (* x x) (car qr) res)
                        (aux (* x x) (car qr) (* x res)))))))
  (aux x n 1))

;;; O(log n)-time and O(1)-space implementation of Fibnacci function
(define (fastfib-aux x11 x12 x21 x22 n r11 r12 r21 r22)
;  (format #t "fastfib-aux: n = ~A\n" n)
   (cond ((= n 0) r21)
         ((= n 1) (+ (* x21 r11) (* x22 r21)))
         (else (let ((qr (integer-divide n 2)))
                 (if (eq? (cdr qr) 0)
                     (fastfib-aux (+ (* x11 x11) (* x12 x21))
                                (+ (* x11 x21) (* x12 x22))
                                (+ (* x21 x11) (* x22 x21))
                                (+ (* x21 x21) (* x22 x22))
                                (car qr) r11 r12 r21 r22)
                     (fastfib-aux (+ (* x11 x11) (* x12 x21))
                                (+ (* x11 x21) (* x12 x22))
                                (+ (* x21 x11) (* x22 x21))
                                (+ (* x21 x21) (* x22 x22))
                                (car qr)
                                (+ (* x11 r11) (* x12 r21))
                                (+ (* x11 r12) (* x12 r22))
                                (+ (* x21 r11) (* x22 r21))
                                (+ (* r21 r12) (* x22 r22))))))))

(define (fastfib n) (fastfib-aux 0 1 1 1 n 1 0 0 1))

(define base (expt 10 8))
(define (+mod x y) (remainder (+ x y) base))
(define (*mod x y) (remainder (* x y) base))

;;; a modified version of 'fastfib-aux' that uses modulo arithmetic
(define (fastmodfib-aux n x11 x12 x21 x22 r11 r12 r21 r22)
;  (format #t "fastfibmod-aux: n = ~A\n" n)
   (cond ((= n 0) r21)
         ((= n 1) (+mod (*mod x21 r11) (*mod x22 r21)))
         (else (let ((qr (integer-divide n 2)))
                 (if (eq? (cdr qr) 0)
                     (fastmodfib-aux (car qr)
                                (+mod (*mod x11 x11) (*mod x12 x21))
                                (+mod (*mod x11 x21) (*mod x12 x22))
                                (+mod (*mod x21 x11) (*mod x22 x21))
                                (+mod (*mod x21 x21) (*mod x22 x22))
                                r11 r12 r21 r22)
                     (fastmodfib-aux (car qr)
                                (+mod (*mod x11 x11) (*mod x12 x21))
                                (+mod (*mod x11 x21) (*mod x12 x22))
                                (+mod (*mod x21 x11) (*mod x22 x21))
                                (+mod (*mod x21 x21) (*mod x22 x22))
                                (+mod (*mod x11 r11) (*mod x12 r21))
                                (+mod (*mod x11 r12) (*mod x12 r22))
                                (+mod (*mod x21 r11) (*mod x22 r21))
                                (+mod (*mod r21 r12) (*mod x22 r22))))))))

(define (fastmodfib n) (fastmodfib-aux n 0 1 1 1 1 0 0 1))


;;;;;;;;;
;;;
;;; Exercise 2: smoothing a function
;;;
;;;;;;;;;

(define (smooth f dx)
  (lambda (x) (/ (+ (f (+ x dx)) (f x) (f (- x dx))) 3)))

;;;;;;;;;
;;;
;;; Exercise 3: repeated application
;;;
;;;;;;;;;

(define (repapp f n)
  (define (aux f n x)
    (if (= n 0) x (aux f (- n 1) (f x))))
  (lambda (x) (aux f n x)))

;;; Here is another implementation

(define (compose f g) (lambda (x) (g (f x))))

(define (repapp f n)
  (if (= n 0)
      (lambda (x) x)
      (compose f (repapp f (- n 1)))))

;;;;;;;;;
;;;
;;; Scheme code for finding congruent numbers
;;;
;;;;;;;;;

(load-option 'format)

(define unit '())
(define true #t)
(define false #f)
(define (int-rat-div x y) (/ x y))
(define (square x) (* x x))

;;; the following function checks whether the number 'area'
;;; can be a congruent number with a witnessing triangle in
;;; which one side is px/qx.
(define (check px qx area)
  (if (= 1 (gcd px qx))
      (let* ((X (int-rat-div px qx))
             (Y (int-rat-div (* 2 (* qx area)) px))
             (py (numerator Y))
             (qy (denominator Y))
             (pqx (+ px qx))
             (pqy (+ py qy)))
        (if (or (< pqy pqx)
                (and (<= pqy pqx) (< py px)))
            unit
            (let* ((pxsq (* px px))
                   (qxsq (* qx qx))
                   (pzsq (+ (square pxsq)
                            (* 4 (* (square qxsq) (square area)))))
                   (pz (floor (sqrt pzsq))))
              (if (= pzsq (* pz pz))
                  (begin (format true
                                 "area = ~S\nX = ~S\nY = ~S\nZ = ~S\n\n"
                                 area X Y (int-rat-div pz (* px qx)))
                         true)
                  false))))
      false))

;;; the function 'find-congruent-number' finds congruent numbers between
;;; 'lower' and 'upper'
(define (find-congruent-numbers lower upper)
  ;;; (loop px sum congs):
  ;;; px: is the numerator
  ;;; px + qx = sum
  ;;; congs stores the congruent numbers that have been found
  (define (loop px sum congs)
    ;;; 'loop-area' tries all numbers between 'lower' and 'upper' that
    ;;; have not been shown to be congruent
    (define (loop-area area congs)
      (if (> area upper)
          (loop (1+ px) sum congs)
          (if (member area congs)
              (loop-area (1+ area) congs)
              (if (check px (- sum px) area)
                  (loop-area (1+ area) (cons area congs))
                  (loop-area (1+ area) congs)))))
    (if (= px sum) (loop 1 (1+ sum) congs) (loop-area lower congs)))
  (loop 1 2 '()))

; (find-congruent-numbers 10 50) quickly finds the the following
; 14 congruent numbers

;area = 20
;X = 3
;Y = 40/3
;Z = 41/3

;area = 15
;X = 4
;Y = 15/2
;Z = 17/2

;area = 30
;X = 5
;Y = 12
;Z = 13

;area = 39
;X = 5/2
;Y = 156/5
;Z = 313/10

;area = 24
;X = 6
;Y = 8
;Z = 10

;area = 21
;X = 7/2
;Y = 12
;Z = 25/2

;area = 14
;X = 8/3
;Y = 21/2
;Z = 65/6

;area = 45
;X = 9/2
;Y = 20
;Z = 41/2

;area = 34
;X = 15/2
;Y = 136/15
;Z = 353/30

;area = 28
;X = 35/6
;Y = 48/5
;Z = 337/30

;area = 41
;X = 40/3
;Y = 123/20
;Z = 881/60

;area = 22
;X = 33/35
;Y = 140/3
;Z = 4901/105

;area = 46
;X = 168/11
;Y = 253/42
;Z = 7585/462

;area = 13
;X = 323/30
;Y = 780/323
;Z = 106921/9690