;;; ursa.lisp.txt — Zoë Trout's favorite programs.
;;;
;;; Name kept from Zoë's canonical file at https://wedgewack.org/ursa.lisp.txt —
;;; same programs, ported to idiomatic Scheme for lumbda. Every iteration
;;; here is a named-let tail call, so TCO keeps the stack flat and portal
;;; determinism stays intact. See ticket 0004 and whitepaper §9.2 for the
;;; CL compat path (loop macro + setf/defun shim) that runs Zoë's original
;;; CL source unchanged.
;;;
;;; Load: (load "examples/ursa.lisp.txt")

;;; ── modular exponentiation ───────────────────────────────────────
;;; (expt-mod b e m) = b^e mod m, in O(log e) multiplies.

(define (expt-mod base exponent modulus)
  (cond ((= modulus 1) 0)
        (else
         (let loop ((b (modulo base modulus))
                    (e exponent)
                    (acc 1))
           (cond ((zero? e) acc)
                 ((odd? e)
                  (loop (modulo (* b b) modulus)
                        (quotient e 2)
                        (modulo (* acc b) modulus)))
                 (else
                  (loop (modulo (* b b) modulus)
                        (quotient e 2)
                        acc)))))))

;;; ── Miller-Rabin primality ───────────────────────────────────────

(define (miller-rabin-base n a d s)
  (let ((x (expt-mod a d n)))
    (if (or (= x 1) (= x (- n 1)))
        #t
        (let loop ((i 0) (x x))
          (cond ((>= i s) #f)
                (else
                 (let ((x2 (expt-mod x 2 n)))
                   (if (= x2 (- n 1))
                       #t
                       (loop (+ i 1) x2)))))))))

;;; n-1 = 2^s * d with d odd, then k witness trials.
(define (miller-rabin-k n k)
  (cond ((<= n 1) #f)
        ((<= n 3) #t)
        ((even? n) #f)
        (else
         (let decompose ((d (- n 1)) (s 0))
           (if (even? d)
               (decompose (quotient d 2) (+ s 1))
               (let witness ((i 0))
                 (cond ((>= i k) #t)
                       (else
                        (let ((a (+ 2 (random-int (- n 2)))))
                          (if (miller-rabin-base n a d s)
                              (witness (+ i 1))
                              #f))))))))))

;;; Variadic wrapper — k defaults to 10.
(define (miller-rabin n . opt)
  (miller-rabin-k n (if (null? opt) 10 (car opt))))

;;; (primep n) returns n if prime, #f otherwise — matches Zoë's CL shape.
(define (primep n . opt)
  (let ((k (if (null? opt) 10 (car opt))))
    (if (miller-rabin-k n k) n #f)))

;;; ── Pollard rho factorization ────────────────────────────────────

;;; rhoff: one factor-finding attempt on n. Returns a non-trivial factor
;;; or #f if the tortoise/hare cycled before separating.
(define (rhoff n)
  (cond ((even? n) 2)
        ((primep n) n)
        (else
         (let ((c (+ 1 (random-int (- n 1)))))
           (define (f z) (modulo (+ (* z z) c) n))
           ;; Birthday-bound iteration cap. Pollard rho expects ≈√n
           ;; iterations before a collision. Capping at 4·√n + 32
           ;; lets honest runs finish while rejecting pathological c
           ;; values quickly; rho's outer retry draws a new c and
           ;; keeps the total work bounded. Without this, a bad c
           ;; can loop for ~n iters allocating let* bindings each
           ;; time, which on the non-GC asm tier builds up heap.
           (let ((cap (+ 32 (* 4 (isqrt n)))))
             (let loop ((x 2) (y 2) (d 1) (k 0))
               (cond ((> k cap) #f)
                     ((not (= d 1))
                      (if (and (< 1 d) (< d n)) d #f))
                     (else
                      (let* ((x2 (f x))
                             (y2 (f (f y)))
                             (d2 (gcd (abs (- x2 y2)) n)))
                        (loop x2 y2 d2 (+ k 1)))))))))))

;;; rho: full factorization via a work-list. Zoë's CL uses an adjustable
;;; vector with vector-push-extend; the Scheme version uses a list as a
;;; stack (cons = push, car+cdr = pop). Same algorithm, no new mutable
;;; type needed — so lumbda's portal continues to serialize cleanly.
(define (rho n)
  (and (> n 1)
       (let loop ((pending (list n)) (results '()))
         (cond ((null? pending) (sort results))
               (else
                (let ((m (car pending)) (rest (cdr pending)))
                  (cond ((= m 1)
                         (loop rest results))
                        ((primep m)
                         (loop rest (cons m results)))
                        ((even? m)
                         (loop (cons (quotient m 2) rest)
                               (cons 2 results)))
                        (else
                         (let retry ((d (rhoff m)))
                           (if d
                               (loop (cons d (cons (quotient m d) rest))
                                     results)
                               (retry (rhoff m))))))))))))

(define (factor n) (rho n))

;;; ── repunit-value ────────────────────────────────────────────────
;;; digit * sum_{i=0..n-1} base^i — a repunit in the given base.

(define (repunit-value n . opt)
  (let* ((digit (if (null? opt) 1 (car opt)))
         (base  (if (or (null? opt) (null? (cdr opt))) 2 (cadr opt))))
    (and (> base digit)
         (let loop ((i 1) (j digit))
           (if (> i (- n 1))
               j
               (loop (+ i 1) (+ j (* digit (expt base i)))))))))

;;; ── Mersenne + Lucas-Lehmer ──────────────────────────────────────

(define (mersenne-number p) (- (expt 2 p) 1))

;;; Lucas-Lehmer recurrence: s_0=4, s_{i+1} = s_i^2 - 2 (mod M_p).
;;; Returns s_{p-2}. M_p is prime iff the residue is zero.
(define (lucas-lehmer-residue p)
  (let ((m (mersenne-number p)))
    (let loop ((i 0) (s 4))
      (if (>= i (- p 2))
          s
          (loop (+ i 1) (modulo (- (* s s) 2) m))))))

(define (lucas-lehmer-primep p)
  (and (primep p)
       (zero? (lucas-lehmer-residue p))))

;;; ── digit conversion ─────────────────────────────────────────────
;;; Zoë's CL version uses defgeneric/defmethod to dispatch on type
;;; (integer → vector, vector → integer, sequence → integer via coerce).
;;; In lumbda we dispatch by predicate — same behavior, no CLOS needed,
;;; no new namespace to serialize. Portal stays clean.

(define (digits->integer seq base)
  (let* ((vec (if (vector? seq) seq (list->vector seq)))
         (len (vector-length vec)))
    (let loop ((i 0) (result 0))
      (if (>= i len)
          result
          (loop (+ i 1) (+ (* result base) (vector-ref vec i)))))))

(define (integer->digits number base)
  (cond ((< number base) (vector number))
        (else
         (let loop ((n number) (acc '()))
           (if (zero? n)
               (list->vector acc)
               (loop (quotient n base)
                     (cons (modulo n base) acc)))))))

(define (digits number base)
  (cond ((integer? number) (integer->digits number base))
        ((vector? number)  (digits->integer number base))
        ((list? number)    (digits->integer number base))
        (else (error "digits: expected integer, vector, or list"))))

;;; ── random-n-bit integer primitives ──────────────────────────────

;;; Uniformly random integer with the top bit set — always in
;;; [2^{n-1}, 2^n). Matches Zoë's CL `of-n-bits`.
(define (of-n-bits n)
  (and (>= n 2)
       (let loop ((i (- n 2)) (sum (expt 2 (- n 1))))
         (if (< i 0)
             sum
             (loop (- i 1)
                   (+ sum (* (random-int 2) (expt 2 i))))))))

;;; Rejection sample: draw random n-bit integers until one is prime.
(define (prime-of-n-bits n)
  (and (>= n 2)
       (let try ()
         (let ((candidate (of-n-bits n)))
           (if (primep candidate) candidate (try))))))