; modular multiplication without overflow
(define (mod-mul a b m) ; assumes a,b < 2^64 and m < 2^63
(define (shift x k)
(let loop ((k k) (x x))
(if (zero? k) x
(let* ((x (* x 2))
(x (if (< x m) x (- x m))))
(loop (- k 1) x)))))
(let* ((two32 (expt 2 32))
(a1 (remainder a two32))
(a2 (quotient a two32))
(b1 (remainder b two32))
(b2 (quotient b two32))
(p11 (modulo (* a1 b1) m))
(p12 (shift (modulo (* a1 b2) m) 32))
(p21 (shift (modulo (* a2 b1) m) 32))
(p22 (shift (modulo (* a2 b2) m) 64))
(s (+ p11 p12))
(s (if (< s m) s (- s m)))
(s (+ s p21))
(s (if (< s m) s (- s m)))
(s (+ s p22))
(s (if (< s m) s (- s m))))
s))
(define a 17259738289493410580109721600)
(define b 10327523882682224844906430464)
(define m 8631375519702822467321987072)
(display (mod-mul a b m))