;;; -*- Mode: Lisp; Syntax: Common-Lisp; Package: Grobner; Base: 10 -*-
#|
	$Id: modular.lisp,v 1.4 2009/01/22 04:03:50 marek Exp $
  *--------------------------------------------------------------------------*
  |  Copyright (C) 1994, Marek Rychlik (e-mail: rychlik@math.arizona.edu)    |
  |    Department of Mathematics, University of Arizona, Tucson, AZ 85721    |
  |                                                                          |
  | Everyone is permitted to copy, distribute and modify the code in this    |
  | directory, as long as this copyright note is preserved verbatim.         |
  *--------------------------------------------------------------------------*
|#

(defpackage "MODULAR"
  (:export modular-division make-modular-division)
  (:use "XGCD" "COMMON-LISP"))
    
(in-package "MODULAR")

#+debug(proclaim '(optimize (speed 0) (debug 3)))
#-debug(proclaim '(optimize (speed 3) (debug 0)))

(defun modular-inverse (x p)
  "Find the inverse of X modulo prime P, using Euclid algorithm."
  (multiple-value-bind (gcd u v)
      (xgcd x p)
    (declare (ignore gcd v))
    (mod u p)))

(defun modular-division (x y p)
  "Divide X by Y modulo prime P."
  (mod (* x (modular-inverse y p)) p))

(defvar *inverse-by-lookup-limit* 100000
  "If prime modulus is < this number then the division algorithm
will use a lookup table of inverses created at the time when field-modulo-prime is called.")
  
(defun make-inverse-table (modulus &aux (table (list 0)))
  "Make a vector of length MODULUS containing all inverses modulo MODULUS,
which should be a prime number. The inverse of 0 is 0."
  (do ((x 1 (1+ x))) ((>= x modulus) (apply #'vector (nreverse table)))
    (push (modular-inverse x modulus) table)))

(defun make-modular-division (modulus)
  "Return a function of two arguments which will perform division
modulo MODULUS. Currently, if MODULUS is < *INVERSE-BY-LOOKUP-LIMIT*
then the returned function does table lookup, otherwise it uses
the Euclid algorithm to find the inverse."
  (cond
   ((>= modulus *inverse-by-lookup-limit*)
    #'(lambda (x y) (modular-division x y modulus)))
   (t
    (let ((table (make-inverse-table modulus)))
      #'(lambda (x y)
	  (mod (* x (svref table y)) modulus))))))