source: CGBLisp/src/RCS/modular.lisp,v@ 1

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
@