;; --------------------------------------------
;; This file contains a bunch of lisp (scheme) code that was found by
;; studying the first chapter of the book "Structure and
;; Interpretation of Computer Programs".

;; Date: 2020-09-17

;; --------------------------------------------

(define (even? n) (= (remainder n 2) 0))
(define (divides? a b) (= (remainder b a) 0))

(define (average x y)
  (/ (+ x y) 2))

(define (average-damp f)
  (lambda (x) (average x (f x))))

(define (abs x)
  (if (< x 0) (- x)
      x))

(define (square n) (* n n))
(define (cube x) (* x x x))

;; computes n!
(define (factorial n)
  (define (iter product counter)
    (if (> counter n)
	product
	(iter (* counter product)
	      (+ counter 1))))
  (iter 1 1))

;; computes the n-th fibonacci number F_n
(define (fib n)
  (define (fib-iter a b count)
    (cond ((= count n) a)
          (else (fib-iter b (+ a b) (+ count 1)))))
  (fib-iter 0 1 0))

;; logarithmic fibonacci
(define (log-fib n)
  (define (fib-iter a b p q count )
    (cond ((= count 0) b)
	  ((even? count)
	   (fib-iter a
		     b
		     (+ (* q q) (* p p))
		     (+ (* 2 p q) (* q q))
		     (/ count 2)))
	  (else (fib-iter (+ (* b q) (* a q) (* a p))
			  (+ (* b p) (* a q))
			  p 
			  q
			  (- count 1)))))  
  (fib-iter 1 0 0 1 n))

;; computes b^n 
(define (expt b n)
  (define (expt-iter a b n)
    (cond ((= n 0) a)
	  ((even? n) (expt-iter a (square b) (/ n 2)))
	  (else (expt-iter (* a b) b (- n 1)))))
  (expt-iter 1 b n))

;; computes base^exp mod m
(define (expmod base exp m)
  (cond ((= exp 0) 1)
	((= exp 1) base)
	((even? exp)
	 (remainder
	  (square (expmod base (/ exp 2) m))
	  m))
	(else
	 (remainder
	  (* base (expmod base (- exp 1) m))
	  m))))

;; computes GCD(a, b)
(define (gcd a b)
  (if (= b 0) a
      (gcd b (remainder a b))))

;; -------

;; basic primality testing
(define (find-divisor n test-divisor)
  (cond ((> (square test-divisor) n) n)
	((divides? test-divisor n) test-divisor)
	(else (find-divisor n (+ test-divisor 1)))))

(define (smallest-divisor n) (find-divisor n 2))

(define (prime? n)
  (= n (smallest-divisor n)))

;; implements Fermat's test for primality testing
(define (fermat-test n)
  (define (try-it a)
    (= (expmod a n n) a))
  (try-it (+ 1 (random (- n 1)))))

(define (fast-prime? n times)
  (cond ((= times 0) #t)
	((fermat-test n) (fast-prime? n (- times 1)))
	(else #f)))

;; check if n is a carmicahel number
(define (carmichael? n)
  (define (carmichael-iter a n)
    (cond ((= a n) #t)
	  ((not (= (expmod a n n) a)) #f)
	  (else (carmichael-iter (+ a 1) n))))
  (carmichael-iter 2 n))

;; Implements Miller-Rabin primality test
(define (decompose n)
  (define (decompose-iter n count)
    (if (and (not (= n 0))
	     (= (remainder n 2) 0))
	(decompose-iter (/ n 2) (+ count 1))
	(list n count)))
  (decompose-iter n 0))

(define (witness a n)
  (define (witness-iter count base)
    (cond ((and (= count 0)
		(not (= base 1))) #t)
	  ((and (= count 0)
		(= base 1)) #f)
	  ((and (not (= base 1))
		(not (or (= base -1) (= base (- n 1))))
		(= (expmod base 2 n) 1)) #t)
	  (else (witness-iter (- count 1) (expmod base 2 n)))))
  (witness-iter (cadr (decompose (- n 1)))
		(expmod a (car (decompose (- n 1))) n)))

(define (miller-rabin n)
  (witness (+ 1 (random (- n 1))) n))

;; -------

;; Approximates sin(x) using the following facts:
;; 
;; 1. For small values of x we have sin(x) \approx x
;; 2. To reduce the size of x we can use the following identity
;;
;;    \sin{x} = 3 \cdot \sin{\frac{x}{3}} - 4 \cdot \sin^3{\frac{x}{3}}
(define (sine angle)
  (define (p x) 
    (- (* 3 x) (* 4 (cube x))))  
  (if (not (> (abs angle) 0.1))
      angle
      (p (sine (/ angle 3.0)))))


;; approximate integrals using simpson-rule
(define (simpson-rule f a b n)
  (define (simpson-term k h) (f (+ a (* k h))))
  (define (simpson-rule-sum total count h)
    (cond ((= count 0) (simpson-rule-sum (f a) (+ count 1) h))
	  ((= count n) (* (/ h 3)
			  (+ total (simpson-term n h))))

	  ((even? count) (simpson-rule-sum (+ total (* 2 (simpson-term count h)))
					   (+ count 1)
					   h))
	  (else (simpson-rule-sum (+ total (* 4 (simpson-term count h)))
				  (+ count 1)
				  h))))
  (simpson-rule-sum 0 0 (/ (- b a) n)))

;; computes the sequence of numbers that converges to pi/8
(define (pi-sum a b)
  (if (> a b)
      0
      (+ (/ 1.0 (* a (+ a 2)))
	 (pi-sum (+ a 4) b))))

;; approximates e - 2 with a continued fraction expansion
(define (cont-frac n d k)
  (define (cont-frac-iter result i)
    (if (= i 0) result
	(cont-frac-iter (/ (n i) (+ (d i) result)) (- i 1))))
  (cont-frac-iter (/ (n k) (d k)) (- k 1)))

(define (d i)
  (cond ((or (= i 1)
	     (= (remainder (- i 1) 3) 2)
	     (= (remainder (- i 1) 3) 0)) 1)
	((= i 2) 2)
	(else
	 (+ 2 (* 2 (floor (/ (- i 1) 3)))))))

;; approximate tan(x) using Lambert's formula
(define (tan-cf x k)
  (define (d i) (- (* 2 i) 1))
  (define (n i)
    (if (= i 1) x
	(- (* x x))))
  (cont-frac n d k))

;; --------

;; computes square root using netwon's method
(define (sqrt x)
  (define (good-enough? guess)
    (< (abs (- (* guess guess) x)) 0.001))
  (define (improve guess) (average guess (/ x guess)))
  (define (sqrt-iter guess)
    (if (good-enough? guess)
	guess
	(sqrt-iter (improve guess))))
  (sqrt-iter 1.0))

;; half-interval method
(define (search f neg-point pos-point)
  (define (close-enough? x y) (< (abs (- x y) 0.001)))
  (let ((midpoint (average neg-point pos-point)))
    (if (close-enough? neg-point pos-point)
	midpoint
	(let ((test-value (f midpoint)))
	  (cond ((positive? test-value)
		 (search f neg-point midpoint)
		 ((negative? test-value)
		  (search f midpoint pos-point)
		  (else midpoint))))))))

(define (half-interval-method f a b)
  (let ((a-value (f a))
	(b-value (f b)))
    (cond ((and (negative? a-value) (positive? b-value)) 
	   (search f a b))
	  ((and (negative? b-value) (positive? a-value))
	   (search f b a))
	  (else
	   (error "Values are not of opposite sign" a b)))))

;; fixed-point method
(define (fixed-point f first-guess)
  (define tolerance 0.00001)
  (define (close-enough? v1 v2) (< (abs (- v1 v2)) tolerance))
  (define (try guess)
    (let ((next (f guess)))
      (if (close-enough? guess next)
	  next
	  (try next))))
  (try first-guess))

;; Newton's method
(define dx 0.00001)

(define (deriv g)
  (lambda (x) (/ (- (g (+ x dx)) (g x)) dx)))

(define (netwon-transform g)
  (lambda (x) (- x (/ (g x) ((deriv g) x)))))

(define (newtons-method g guess)
  (fixed-point (newton-transform g) guess))

;; approximates the golden ratio (1.6180...)
;; (fixed-point (lambda (x) (+ 1 (/ 1 x))) 1.0)

;; computes the n-th root
(define (compose f g)
  (lambda (x) (f (g x))))

(define (repeated f n)
  (define (repeated-iter comp count)
    (if (= count n)
	comp
	(repeated-iter (compose f comp) (+ count 1))))
  (repeated-iter f 1))

(define (nth-root n x)
  (fixed-point
   ((repeated average-damp (floor (base-log 2 n))) (lambda (y) (/ x (expt y (- n 1)))))
   1.0))

;; -------------------
;; Misc
;; -------------------

;; Compute number of ways you can change a given amount of money
(define (count-change amount) (cc amount 5))

(define (cc amount kinds-of-coins)
  (cond ((= amount 0) 1)
	((or (< amount 0) (= kinds-of-coins 0)) 0)
	(else (+ (cc amount
		     (- kinds-of-coins 1))
		 (cc (- amount (first-denomination kinds-of-coins))
		     kinds-of-coins)))))

(define (first-denomination kinds-of-coins)
  (cond ((= kinds-of-coins 1) 1)
	((= kinds-of-coins 2) 5)
	((= kinds-of-coins 3) 10)
	((= kinds-of-coins 4) 25)
	((= kinds-of-coins 5) 50)))

;; Ackermann's function
(define (A x y)
  (cond ((= y 0) 0)
	((= x 0) (* 2 y))
	((= y 1) 2)
	(else (A (- x 1) (A x (- y 1))))))

;; Pascal triangle
(define (pascal-triangle row col)
  (define (pascal-triangle-rec row col)
    (cond  ((= row 1) 1)
	   ((= col 1) 1)
	   ((= col row) 1)
	   (else 
	    (+ 
	     (pascal-triangle-rec (- row 1) (- col 1))
	     (pascal-triangle-rec (- row 1) col)))))
  (if (> col row) -1
      (pascal-triangle-rec row col)))