;;; Arithmetic modulo 2^31 - 1.  Assumes a Scheme with >48-bit fixnums,
;;; i.e. MIT Scheme on amd64.
;;;
;;; Public domain, or if you can't handle that then cc0
;;; <https://creativecommons.org/publicdomain/zero/1.0/>.

(declare (usual-integrations))

(define-integrable (reduce311 x)
  (fix:+ (fix:lsh x -31) (fix:and x #x7fffffff)))

(define-integrable (add311 a b)
  (reduce311 (fix:+ a b)))

(define-integrable (neg311 x)
  (fix:- #x7fffffff x))

(define-integrable (mul311 a b)
  #;
  (let ((p (* a b)))
    (add311 (shift-right p 31) (bitwise-and p #x7fffffff)))
  (let ((ah (fix:lsh a -16))
        (al (fix:and a #xffff))
        (bh (fix:lsh b -16))
        (bl (fix:and b #xffff)))
    (declare (integrate ah al bh bl))
    (let ((p0 (fix:lsh (fix:* ah bh) 1))        ;2^32 \equiv 2^1 (mod 2^31 - 1)
          (p1 (fix:lsh (fix:* ah bl) 16))
          (p2 (fix:lsh (fix:* al bh) 16))
          (p3 (fix:* al bl)))
      (declare (integrate p0 p1 p2 p3))
      (add311 (add311 p0 p1) (add311 p2 p3)))))

(define (order311 x)
  (let loop ((n 1) (xn x))
    (if (fix:= xn 1)
        n
        (loop (fix:+ n 1) (mul311 xn x)))))