(module ma158 mzscheme (require (lib "etc.ss") (only (lib "1.ss" "srfi") any every append-map fold find filter take) (only (lib "43.ss" "srfi") vector-map vector-fold)) (define currently-deriving (make-parameter #t)) (define (derive . args) (when (currently-deriving) (apply printf args))) (define pairwise (opt-lambda (pred? ls [error-msg #f]) (or (null? ls) (and (every (lambda (x) (or (pred? (car ls) x) (and error-msg (error (format "~a and ~a are not ~a" (car ls) x error-msg))))) (cdr ls)) (pairwise pred? (cdr ls) error-msg))))) ;;====================================================================== ;; generally useful mod arithmetic stuff (define div quotient) (define mod remainder) (define ^ expt) (define square-and-multiply (opt-lambda (base power modulus [acc 1]) (cond [(zero? power) acc] [(odd? power) (square-and-multiply base (sub1 power) modulus ((* acc base) . mod . modulus))] [else (square-and-multiply ((* base base) . mod . modulus) (/ power 2) modulus acc)]))) (define naive-^-mod (opt-lambda (base power modulus [acc 1]) (if (zero? power) acc (naive-^-mod base (sub1 power) modulus ((* acc base) . mod . modulus))))) (define ^-mod square-and-multiply) (define (normalized-mod n p) (let ([n (n . mod . p)]) (if (negative? n) (+ p n) n))) ; generates a random number between 0 and limit; the maximum limit is ; 2147483647 (define (gen-random limit) (random (min 2147483647 limit))) ;;====================================================================== ;; stuff for primitive roots (define (powers-in-group n m) (let loop ([ls (list n)]) (let ([next (mod (* (car ls) n) m)]) (if (memq next ls) ls (loop (cons next ls)))))) (define (primitive-root? n m) (let ([all-elements (build-list (sub1 m) add1)] [powers (powers-in-group n m)]) (and (every (lambda (n) (memq n powers)) all-elements) #t))) ; #t so we get a truth value, rather than the last memq result ;;====================================================================== ;; an exponential-time brute force discrete-log solver (define (naive-dlp base value modulus) (let loop ([exp 0]) (cond [(= exp modulus) (error "never found an exponent...huh")] [(= (mod (expt base exp) modulus) value) exp] [else (loop (add1 exp))]))) ;;====================================================================== ;; a . divides? . d <=> a | b (define (divides? a b) (zero? (mod b a))) (define (relatively-prime? a b) (= 1 (gcd a b))) ;;====================================================================== ;; ord_p(n) operator (define (ord p n) (if (divides? p n) (add1 (ord p (/ n p))) 0)) ;;====================================================================== ;; Extended GCD (solves ua + vb = 1); used for inverse-mod (define (extended-gcd a b) (if (zero? (mod a b)) (values 0 1) (let-values ([(x y) (extended-gcd b (mod a b))]) (values y (- x (* y (div a b))))))) (define (inverse-mod n p) (if (= n 1) 1 (let-values ([(inv invertible?) (extended-gcd n p)]) (if (or (= inv 0) (= invertible? 0)) (error 'inverse-mod "~a is not invertible mod ~a" n p) (inv . normalized-mod . p))))) ;;====================================================================== ;; Affine encryption algorithm (define (affine-e p k1 k2 m) (mod (+ (* k1 m) k2) p)) (define (affine-d p k1 k2 c) (mod (* (- c k2) (k1 . inverse-mod . p)) p)) ;;====================================================================== ;; ElGamal encryption algorithm (define (elgamal-pubkey privkey base p) (^-mod base privkey p)) (define (elgamal-e pubkey k m base p) (let ([c1 (^-mod base k p)] [c2 ((* m (^ pubkey k)) . mod . p)]) (values c1 c2))) (define (elgamal-d c1 c2 privkey p) (let ([c1-inv ((^-mod c1 privkey p) . inverse-mod . p)]) ((* c1-inv c2) . mod . p))) ;;====================================================================== ;; a congruence solver, via Chinese remainder theorem ;; ;; USE: ;; ;; (x-congruent-to (a mod b) ;; (c mod d) ;; ..... ;; (y mod z)) ;; ;; will return the lowest such x (mod bdf...z) if one exists, or will signal ;; an error if on of b, d, f, ..., z aren't pairwise relatively prime. (define (solve-2-congruence a amod b bmod) (unless (relatively-prime? amod bmod) (error 'solve-2-congruence "~a and ~a aren't relatively prime; gcd is ~a" amod bmod (gcd amod bmod))) ; x = a (mod amod) ; x = b (mod bmod) ; ; so x = y*amod + a = b (mod bmod) ; y*amod = b - a (mod bmod) ; y = (b - a)*amod^(-1) ; ; then we solve for x and return the smallest possible such value (mod amod*bmod) (let* ([y (* (- b a) (amod . inverse-mod . bmod))] [x (+ (* y amod) a)]) (x . normalized-mod . (* amod bmod)))) (define (solve-congruences a amod congs) (if (null? congs) ; single congruence case a (let* ([b-cong (car congs)] [b (car b-cong)] [bmod (cadr b-cong)]) (solve-congruences (solve-2-congruence a amod b bmod) (* amod bmod) (cdr congs))))) (define-syntax x-congruent-to (syntax-rules (mod) [(_ (a mod amod) (a* mod amod*) ...) (begin (pairwise relatively-prime? (list amod amod* ...) "relatively prime") (solve-congruences a amod (list (list a* amod*) ...)))])) ;;====================================================================== ;; a square root finder (define (naive-sqrt a p) (let ([a (a . normalized-mod . p)]) (let loop ([root 1]) (cond [(= root p) (error 'naive-sqrt "no root found for ~a in ~a" a p)] [(= ((* root root) . mod . p) a) root] [else (loop (add1 root))])))) (define (sqrt-mod-composite a . factors) (when (null? factors) (error 'sqrt-mod-composite "no factors given")) (display (cons 'x-congruent-to (map (lambda (factor) (list (a . mod . factor) 'mod factor)) factors))) (newline) (let* ([sqrts (map (lambda (p) (if ((p . mod . 4) . = . 3) (^-mod (a . mod . p) ((add1 p) . / . 4) p) (naive-sqrt (a . mod . p) p))) factors)] [a (car sqrts)] [amod (car factors)] [congs (map list (cdr sqrts) (cdr factors))]) (display (cons 'sqrt-x-congruent-to (map (lambda (sqrt factor) (list sqrt 'and (- factor sqrt) 'mod factor)) sqrts factors))) (newline) (solve-congruences a amod congs))) ;;====================================================================== ;; Pohlig-Hellman DLP solver (define (ordered-dlp h b p order) (display (format "We solve the DLP $~a^x \\equiv ~a (\\mod ~a)$ -- the order here is ~a.~n" h b p order)) (let loop ([exp 0]) (cond [(= exp order) (error "never found an exponent...huh")] [(= (^-mod h exp p) b) exp] [else (loop (add1 exp))]))) (define (next-x-for-factored-dlp h b p q e x i) (let* ([q^e-1 (^ q (sub1 e))] [q^e-i (^ q (- e 1 i))] [h^-x ((^ h x) . inverse-mod . p)] [bh^-x ((* b h^-x) . mod . p)] [b (^-mod bh^-x q^e-i p)]) (let ([xi (* (ordered-dlp (^-mod h q^e-1 p) b p q) (^ q i))]) (display (format "Given x = ~a, we compute the $x_{~a}$ to be the solution to $(~a^{~a^{~a}})^{x_~a} = (~a \\cdot {~a^{~a}}^{-1})^{~a^{~a}}$; it is ~a.~n" x i h q e i b h x q (- e 1 i) xi)) (+ x xi)))) (define (factored-dlp h b p q e) (let loop ([i 0] [x 0]) (if (= i e) x (loop (add1 i) (next-x-for-factored-dlp h b p q e x i))))) (define (pohlig-hellman-dlp base value . factors) (let* ([p-1 (fold (lambda (factor acc) (* (^ (car factor) (cadr factor)) acc)) 1 factors)] [p (add1 p-1)] [mod-logs (map (lambda (factor) (let* ([fac-base (car factor)] [fac-expt (cadr factor)] [factor (^ fac-base fac-expt)] [exponent (/ p-1 factor)] [h (^-mod base exponent p)] [b (^-mod value exponent p)]) (display (format "~nThen we must solve for the root $~a^{~a}$.~n" fac-base fac-expt)) (let ([log (factored-dlp h b p fac-base fac-expt)]) (display (format "Our final x for the root $~a^{~a}$ is thus ~a.~n" fac-base fac-expt log)) log))) factors)] [factors (map (lambda (factor) (^ (car factor) (cadr factor))) factors)] [a (car mod-logs)] [amod (car factors)] [congs (map list (cdr mod-logs) (cdr factors))]) (display (format "~nAnd now we can simply solve the congruence~n\\begin{eqnarray*}~nx & \\equiv & ~a (\\mod ~a) ~a~n\\end{eqnarray*}~nwith the result " a amod (apply string-append (map (lambda (equiv) (let ([val (car equiv)] [modulus (cadr equiv)]) (format "\\\\~n& \\equiv & ~a (\\mod ~a) " val modulus))) congs)))) (solve-congruences a amod congs))) ;;====================================================================== ;; Miller-Rabin test ; returns values 2^k q, such that 2^k*q = n and q is odd (define split-odd-factor (opt-lambda (n [k 0]) (if (even? n) (split-odd-factor (/ n 2) (add1 k)) (values k n)))) (define (miller-rabin-witness? n a) (let ([a-orig a]) (cond [(and (even? n) (not (= n 2))) (derive "$~a$ needs no witness -- it's even!~n" n) #t] ; even numbers are obviously composite, excepting 2 ;) [(< 1 (gcd n a) n) (derive "~a and ~a have gcd $1 < ~a < ~a$, so ~a is a witness.~n" n a (gcd n a) n a) #t] ; if there's a common divisor, we're also composite [else ; the interesting part: first, get k and q s.t. 2^k*q = n-1, q is odd (let-values ([(k q) (split-odd-factor (sub1 n))]) (derive "$~a = 2^{~a} \\cdot ~a$. " n k q) (let ([a^q (^-mod a q n)]) (derive "First, we take $a = ~a^{~a} = ~a$. " a q a^q) (cond [(= 1 a^q) (derive "Now $a \\equiv -1 (\\mod ~a)$, so $~a$ may be prime; in any case, $~a$ isn't a witness. " n n a-orig) #f] [else (let loop ([a a^q] [i 0]) (cond [(= (sub1 n) a) (derive "Now $a \\equiv -1 (\\mod ~a)$, so $~a$ may be prime; in any case, $~a$ isn't a witness. " n n a-orig) #f] [(= i k) (derive "We've gone through ~a iterations, but ~a still didn't invert -- that is, $~a^{2^{~a} \\cdot ~a} \\equiv (\\mod ~a)$, so ~a is a witness that ~a is composite." k a-orig a-orig k q n a-orig n) #t] [else (let ([next-a (^-mod a 2 n)]) (derive "$~a^{2^{~a} \\cdot ~a} \\not\\equiv 1~~\\text{or}~~-1 (\\mod ~a)$, so we continue to the next iteration, and our new $a = ~a^2 = ~a$. " a-orig k q n a next-a) (loop next-a (add1 i)))]))])))]))) ;;====================================================================== ;; Probable prime generation ; we start with a prime base and a relative prime to that base -- this'll speed up our search (define 2*3*5*7*11 (* 2 3 5 7 11)) (define rp-2*3*5*7*11 1139) (define (random-witness limit) (+ 2 (gen-random limit))) (define (passes-n-miller-rabin-tests? p n) (let ([witnesses (build-list n (lambda (i) (random-witness (sub1 p))))]) (cond [(find (lambda (a) (parameterize ((currently-deriving #f)) (miller-rabin-witness? p a))) witnesses) => (lambda (a) (derive "~a is a witness for ~a.~n" a p) #f)] [else (derive "None of ~a are witnesses for ~a.~n" witnesses p) #t]))) (define (miller-rabin-probable-prime scale prime-base rel-prime) (let loop ([j 0]) (let ([p (+ (* (+ scale j) prime-base) rel-prime)]) (if (passes-n-miller-rabin-tests? p 10) (if (passes-n-miller-rabin-tests? p 100) ; can't hurt to be sure ;) p (loop (add1 j))) (loop (add1 j)))))) ;;====================================================================== ;; Difference of squares (define (perfect-square? n) (exact? (sqrt n))) (define (perfect-squares-up-to N lim) (map (lambda (n) (list n (perfect-square? n))) (build-list 30 (lambda (x) (+ N (^ (add1 x) 2)))))) ;;====================================================================== ;; Quadratic reciprocity (define (jacobi q p) (unless (> p 1) (error 'jacobi "expected p > 1, got (jacobi ~a ~a)" q p)) (cond [(= 2 q) (case (p . mod . 8) [(1 7) 1] [(3 5) -1] [else (error 'jacobi "expected p odd prime, got ~a" p)])] [(= 1 q) 1] ; pretty trivial... [(= -1 q) (case (p . mod . 4) [(1) 1] [(3) -1] [else (error 'jacobi "expected p odd prime, got ~a" p)])] [(< q 0) (* (jacobi -1 p) (jacobi (* -1 q) p))] [(even? q) (derive "$~a = 2 \\cdot ~a$, so we compute $\\jacobi{2}{~a}\\cdot\\jacobi{~a}{~a}$. " q (/ q 2) p (/ q 2) p) (let* ([j2 (jacobi 2 p)] [jq/2 (jacobi (/ q 2) p)] [result (* j2 jq/2)]) (derive "We thus find that $\\jacobi{~a}{~a} = \\jacobi{2}{~a}\\cdot\\jacobi{~a}{~a} = ~a \\cdot ~a = ~a.$ " q p p (/ q 2) p j2 jq/2 result) result)] [(even? p) (derive "$~a = 2 \\cdot ~a$, so we compute $\\jacobi{~a}{2}\\cdot\\jacobi{~a}{~a}$. " p (/ p 2) q q (/ p 2)) (let* ([j2 (jacobi q 2)] [jp/2 (jacobi q (/ p 2))] [result (* j2 jp/2)]) (derive "We thus find that $\\jacobi{~a}{~a} = \\jacobi{~a}{2}\\cdot\\jacobi{~a}{~a} = ~a \\cdot ~a = ~a.$ " q p q q (/ p 2) j2 jp/2 result) result)] [(or (= (p . mod . 4) 1) (= (q . mod . 4) 1)) (derive "$~a \\mod 4 \\equiv 1$, so $\\jacobi{~a}{~a} = \\jacobi{~a}{~a}$. " (if (= (p . mod . 4) 1) p q) q p (p . mod . q) q) (jacobi (p . mod . q) q)] [(and (= (p . mod . 4) (q . mod . 4) 3)) (derive "$~a \\mod 4 \\equiv ~a \\mod 4 \\equiv 1$, so $\\jacobi{~a}{~a} = -\\jacobi{~a}{~a}$. " q p q p (p . mod . q) q) (let* ([pmq (jacobi (p . mod . q) q)] [result (* -1 pmq)]) (derive "We find that $\\jacobi{~a}{~a} = ~a$; so we find that $\\jacobi{~a}{~a} = ~a.$ " (p . mod . q) q pmq q p result) result)] [else (error 'jacobi "can't further reduce (jacobi ~a ~a) without factoring ~a" q p q)])) ;;====================================================================== ;; Goldwasser-Micali encryption (define goldwasser-micali-e (opt-lambda (N a m [r #f]) (let* ([r (or r (add1 (gen-random (sub1 N))))] [r^2 (^-mod r 2 N)]) ((* (if m a 1) r^2) . mod . N)))) (define (goldwasser-micali-d p c) (= (jacobi c p) -1)) ;;====================================================================== ;; Continued fractions (define (continued-fraction->latex ls) (cond [(null? ls) ""] [(null? (cdr ls)) (format "~a" (car ls))] [else (format "~a + \\frac{1}{~a}" (car ls) (continued-fraction->latex (cdr ls)))])) (define (evaluate-continued-fraction . ls) (cond [(null? ls) 0] [(null? (cdr ls)) (car ls)] [else (+ (car ls) (/ 1 (apply evaluate-continued-fraction (cdr ls))))])) (define (sublists ls) (map (lambda (i) (take ls i)) (build-list (length ls) add1))) (define (evaluate-each-continued-fraction ls) (map (lambda (ls) (apply evaluate-continued-fraction ls)) (sublists ls))) (define (expand-continued-fraction n) (cond [(zero? n) (values 0 0)] [(integer? n) (values n 0)] [(negative? n) (raise-type-error 'expand-continued-fraction "expected n>=0, got ~a" n)] [else (let* ([fl (floor n)] [fr (- n fl)] [recip (/ 1 fr)]) (derive "$~a = ~a + ~a = ~a + \\frac{1}{~a}$. " n fl fr fl recip) (values fl recip))])) (define (continued-fraction n) (unless (rational? n) (raise-type-error 'continued-fraction "expected rational n, got ~a" n)) (let ([n-pos (max n (- n))]) (let-values ([(val frac) (expand-continued-fraction n-pos)]) (let ([val (* val (/ n n-pos))]) (cond [(zero? frac) `(,val)] [else (cons val (continued-fraction frac))]))))) ;;====================================================================== ;; Vectors (define (vec-length v) (sqrt (vector-fold (lambda (i sum n) (+ sum (* n n))) 0 v))) (define (dot-product v1 v2) (vector-fold (lambda (i sum n) (+ sum n)) 0 (vector-map (lambda (i x1 x2) (* x1 x2)) v1 v2))) )