#lang racket

;;;;;;;;;;
;; 1.29 ;;
;;;;;;;;;;

(define (sum term a next b)
  (if (> a b)
    0
    (+ (term a)
       (sum term (next a) next b))))

(define (integral f a b dx)
  (define (add-dx x) (+ x dx))
  (* (sum f (+ a (/ dx 2.0)) add-dx b)
     dx))

(define (cube x) (* x x x))

(integral cube 0 1 0.01)
;; 0.24998750000000042
(integral cube 0 1 0.001)
;; 0.249999875000001

;; Simpson's Rule:  Estimate the definite integral of f(x) from x = a to x = b as
;; [(b - a) / (* 3 n)] times
;; [(sum all terms) + (sum all inner terms) + 2 * (sum all odd inner terms)]
;; where n is an even number

(define (simpson-integral f a b even-n)
  (define h (/ (- b a) even-n))
  (define (add-dx x) (+ x h))
  (define (add-two-dx x) (+ x h h))
  (* (/ h 3.0)
     (+ (sum f a add-dx b)
        (sum f (+ a h) add-dx (- b h))
        (* 2 (sum f (+ a h) add-two-dx (- b h))))))

(simpson-integral cube 0 1 100)
;; 0.25
(simpson-integral cube 0 1 1000)
;; 0.25

;;;;;;;;;;
;; 1.30 ;;
;;;;;;;;;;

;; iterative-sum is tail-recursive hence iterative

(define (iterative-sum term a next b)
  (define (iter a result)
    (if (> a b)
      result
      (iter (next a) (+ result (term a)))))
  (iter a 0))

(define (identity x) x)

(iterative-sum identity 1 add1 100) ; sum numbers from 1 to 100
;; 5050

;;;;;;;;;;
;; 1.31 ;;
;;;;;;;;;;

(define (product term a next b)
  (if (> a b)
    1
    (* (term a) (product term (next a) next b))))

(define (factorial n)
  (product identity 1 add1 n))

(factorial 3)
;; 6
(factorial 4)
;; 24
(factorial 5)
;; 120

(define (iterative-product term a next b)
  (define (iter a result)
    (if (> a b)
      result
      (iter (next a) (* result (term a)))))
  (iter a 1))

(iterative-product identity 1 add1 5) ; should be factorial(5)
;; 120

;; Wallis' formula:
;; pi / 4 = (2 / 3) * (4 / 3) * (4 / 5) * (6 / 5) * (6 / 7) * (8 / 7) ...

(define (wallis-term i)
  (define j (* 2 (add1 (quotient i 2))))
  (if (odd? i)
    (/ j (add1 j))
    (/ j (sub1 j))))


(for ([i (in-range 1 7)]) (displayln (wallis-term i)))
;; 2/3
;; 4/3
;; 4/5
;; 6/5
;; 6/7
;; 8/7

(* 4 (product wallis-term 1.0 add1 10))
;; 3.275101041334807
(* 4 (product wallis-term 1.0 add1 100))
;; 3.1570301764551645
(* 4 (product wallis-term 1.0 add1 1000))
;; 3.143160705532257
(* 4 (product wallis-term 1.0 add1 10000))
;; 3.1417497057379635
(* 4 (product wallis-term 1.0 add1 100000))
;; 3.1416083612780903

;;;;;;;;;;
;; 1.32 ;;
;;;;;;;;;;

(define (accumulate combiner null-value term a next b)
  (if (> a b)
    null-value
    (combiner (term a)
              (accumulate combiner null-value term (next a) next b))))

(define (new-sum term a next b)
  (accumulate + 0 term a next b))

(new-sum identity 1 add1 100) ; add numbers from 1 to 100
;; 5050

(define (new-product term a next b)
  (accumulate * 1 term a next b))

(new-product identity 1 add1 5) ; should be factorial(5)
;; 120

(define (iterative-accumulate combiner null-value term a next b)
  (define (iter a result)
    (if (> a b)
      result
      (iter (next a) (combiner result (term a)))))
  (iter a null-value))

(iterative-accumulate + 0 identity 0 add1 100) ; add numbers from 1 to 100
;; 5050

(iterative-accumulate * 1 identity 1 add1 5) ; should be factorial(5)
;; 120

;;;;;;;;;;
;; 1.33 ;;
;;;;;;;;;;

(define (filtered-accumulate filter? combiner null-value term a next b)
  (define (iter a result)
    (cond [(> a b) result]
          [(filter? a)
           (iter (next a) (combiner result (term a)))]
          [else (iter (next a) result)]))
  (iter a null-value))

(define (prime? n [times 10]) ; miller-rabin test from 1.26
  (define (expmod base expo m)
    (cond [(zero? expo) 1]
          [(even? expo)
           (define u (expmod base (/ expo 2) m))
           (define u-squared (remainder (sqr u) m))
           (if (and (= u-squared 1)
                    (> u 1)
                    (< u (sub1 m)))
             0
             u-squared)]
          [else (remainder (* base (expmod base (sub1 expo) m)) m)]))
  (define (try-it a)
    (= (expmod a (sub1 n) n) 1))
  (if (< n 10)
    (member n '(2 3 5 7))
    (for/and ([i (in-range times)])
             ; do not test 0, 1, or n - 1
             (try-it (+ 2 (random (min 4294967087 (- n 3))))))))

(define (sum-of-squares-of-primes a b)
  (filtered-accumulate prime? + 0 sqr a add1 b))

(sum-of-squares-of-primes 1 4) ; 2 ^ 2 + 3 ^ 2 = 13
;; 13
(sum-of-squares-of-primes 1 10) ; 4 + 9 + 25 + 49 = 87
;; 87

(define (product-of-relative-primes n)
  (define (relatively-prime? i)
    (= (gcd i n) 1))
  (filtered-accumulate relatively-prime? * 1 identity 1 add1 n))

(product-of-relative-primes 6) ; 1 * 5 = 5
;; 5
(product-of-relative-primes 7) ; 1 * 2 * 3 * 4 * 5 * 6 = 720 
;; 720
(product-of-relative-primes 8) ; 1 * 3 * 5 * 7 = 105
;; 105
(product-of-relative-primes 20)
;; 8729721