head	1.4;
access;
symbols;
locks; strict;
comment	@;;; @;


1.4
date	2009.01.22.04.03.50;	author marek;	state Exp;
branches;
next	1.3;

1.3
date	2009.01.19.09.26.32;	author marek;	state Exp;
branches;
next	1.2;

1.2
date	2009.01.19.07.39.06;	author marek;	state Exp;
branches;
next	1.1;

1.1
date	2009.01.19.06.44.13;	author marek;	state Exp;
branches;
next	;


desc
@@


1.4
log
@*** empty log message ***
@
text
@;;; -*- Mode: Lisp; Syntax: Common-Lisp; Package: Grobner; Base: 10 -*-
#|
	$Id$
  *--------------------------------------------------------------------------*
  |  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))))))@


1.3
log
@*** empty log message ***
@
text
@d19 2
a20 2
;;(proclaim '(optimize (speed 0) (debug 3)))
(proclaim '(optimize (speed 3) (debug 0)))
@


1.2
log
@*** empty log message ***
@
text
@d19 2
a20 1
(proclaim '(optimize (speed 0) (debug 3)))
@


1.1
log
@Initial revision
@
text
@d3 1
a3 1
	$Id: modular.lisp,v 1.6 1997/12/13 15:55:32 marek Exp $
d19 2
@