; primitive roots
(define (factors n) ; 2,3,5-wheel
(let ((wheel (vector 1 2 2 4 2 4 2 4 6 2 6)))
(let loop ((n n) (f 2) (w 0) (fs (list)))
(if (< n (* f f)) (reverse (cons n fs))
(if (zero? (modulo n f))
(loop (/ n f) f w (cons f fs))
(loop n (+ f (vector-ref wheel w))
(if (= w 10) 3 (+ w 1)) fs))))))
(define (unique eql? xs)
(cond ((null? xs) '())
((null? (cdr xs)) xs)
((eql? (car xs) (cadr xs))
(unique eql? (cdr xs)))
(else (cons (car xs) (unique eql? (cdr xs))))))
(define (expm b e m)
(define (m* x y) (modulo (* x y) m))
(cond ((zero? e) 1)
((even? e) (expm (m* b b) (/ e 2) m))
(else (m* b (expm (m* b b) (/ (- e 1) 2) m)))))
(define (prim-roots p)
(let ((fs (unique = (factors (- p 1)))))
(let loop ((a 1) (ps fs))
(if (pair? ps)
(if (= (expm a (/ (- p 1) (car ps)) p) 1)
(loop (+ a 1) fs) ; try next a
(loop a (cdr ps)))
(let loop ((m 1) (rs (list)))
(if (= m p) (sort rs <)
(if (= (gcd m (- p 1)) 1)
(loop (+ m 1) (cons (expm a m p) rs))
(loop (+ m 1) rs))))))))
(display (prim-roots 47)) (newline)