-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathlib-gfp.lisp
178 lines (153 loc) · 10.1 KB
/
lib-gfp.lisp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
;; -*- Mode:Lisp; Syntax:ANSI-Common-LISP; Coding:us-ascii-unix; fill-column:158 -*-
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;; @file lib-gfp.lisp
;; @author Mitch Richling <https://www.mitchr.me>
;; @brief Interactive GF(p) library -- modular arithmatic.@EOL
;; @std Common Lisp
;; @see tst-gfp.lisp
;; @copyright
;; @parblock
;; Copyright (c) 1997,1998,2004,2010,2013,2015, Mitchell Jay Richling <https://www.mitchr.me> All rights reserved.
;;
;; Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
;;
;; 1. Redistributions of source code must retain the above copyright notice, this list of conditions, and the following disclaimer.
;;
;; 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions, and the following disclaimer in the documentation
;; and/or other materials provided with the distribution.
;;
;; 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software
;; without specific prior written permission.
;;
;; THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
;; IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
;; LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
;; OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
;; LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
;; DAMAGE.
;; @endparblock
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defpackage :MJR_GFP
(:USE :COMMON-LISP
:MJR_INTU)
(:DOCUMENTATION "Brief: Interactive GF(p) library -- modular arithmatic.;")
(:EXPORT #:mjr_gfp_help
#:mjr_gfp_simplify
#:mjr_gfp_+ #:mjr_gfp_- #:mjr_gfp_* #:mjr_gfp_/
#:mjr_gfp_iexpt #:mjr_gfp_imul
#:mjr_gfp_< #:mjr_gfp_> #:mjr_gfp_=
#:mjr_gfp_zerop #:mjr_gfp_onep #:mjr_gfp_gfpp
#:mjr_gfp_divides?
))
(in-package :MJR_GFP)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_gfp_help ()
"Basic arithmetic over prime order finite fields (i.e. known as prime order modular arithmetic in computer science circles).
I should have said 'in computer science FIELDS' -- HA. No. Really. That was funny! Well, OK... Fine...
While this library may be used interactively, it proimarly provides a support roll for other libraries like MJR_POLYGFP and MJR_PRIME."
(documentation 'mjr_gfp_help 'function))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defconstant +mjr_gfp_zero+ 0
"Additive unit in the field GF(p)")
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defconstant +mjr_gfp_one+ 1
"Multiplicative unit in the field GF(p)")
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_gfp_gfpp (p a)
"non-NIL if object is, or can be safely coerced into an element of GF(p).
In particular, the return will be one of the following:
* 1 ...... An integer with no need of simplification
* 2 ...... An integer in need of simplification
* 3 ...... A non-complex number (it will be truncated by mjr_gfp_simplify)
* nil .... All other cases."
(if (integerp a)
(if (and (<= 0 a) (< a p))
1
2)
(if (and (numberp a) (not (complexp a)))
3)))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_gfp_simplify (p a)
"Truncate lisp numbers to integers if required and then simplify the result mod P"
(if (integerp a)
(mod a p)
(mjr_gfp_simplify p (truncate a))))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_gfp_= (p a b)
"non-NIL if A=B mod P"
(= (mjr_gfp_simplify p a) (mjr_gfp_simplify p b)))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_gfp_onep (p a)
"non-NIL if A=1 mod P"
(mjr_gfp_= p a +mjr_gfp_one+))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_gfp_zerop (p a)
"non-NIL if A=0 mod P"
(mjr_gfp_= p a +mjr_gfp_zero+))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_gfp_divides? (p a b)
"non-NIL if a divides b --- $a|b$ if $ak=b$ for some integer $k$. $\\mathrm{GF}(p)$ is a FIELD, so $a|b$ unless $a=0$ and $b\\neq0$."
(if (mjr_gfp_zerop p a)
(mjr_gfp_zerop p b)
't))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_gfp_< (p a b)
"non-NIL if A<B mod P"
(< (mjr_gfp_simplify p a) (mjr_gfp_simplify p b)))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_gfp_> (p a b)
"non-NIL if A>B mod P"
(> (mjr_gfp_simplify p a) (mjr_gfp_simplify p b)))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_gfp_+ (p &rest rest)
"Compute (X+Y mod P). Return additive identity when given no arguments."
(reduce (lambda (x y) (mjr_gfp_simplify p (+ x (mjr_gfp_simplify p y)))) rest :initial-value +mjr_gfp_zero+))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_gfp_- (p &rest rest)
"Compute (X-Y mod P). Return additive inverse when given a single argument."
(if rest
(reduce (lambda (x y) (mjr_gfp_simplify p (- x (mjr_gfp_simplify p y)))) rest :initial-value +mjr_gfp_zero+)
(error "mjr_gfp_-: Missing argument")))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_gfp_* (p &rest rest)
"Compute (X*Y mod P). Return multiplicative identity when given no arguments."
(reduce (lambda (x y) (mjr_gfp_simplify p (* x (mjr_gfp_simplify p y)))) rest :initial-value +mjr_gfp_one+))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_gfp_/ (p &rest rest)
"Compute (X/Y mod P). Returns multiplicative inverse when given a single argument."
(if (cdr rest)
(reduce (lambda (x y) (mjr_gfp_simplify p (* x (mjr_gfp_/ p y)))) (cdr rest) :initial-value (mjr_gfp_simplify p (car rest)))
(if rest
(let ((tmp (mjr_gfp_simplify p (car rest))))
(if (= tmp +mjr_gfp_zero+)
(error "mjr_gfp_/: DIVISION-BY-ZERO!")
(multiple-value-bind (s1 s2 gcd) (mjr_intu_extended-gcd tmp p)
(declare (ignore s2))
(if (= 1 gcd)
(mjr_gfp_simplify p s1)))))
(error "mjr_gfp_/: Missing argument"))))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_gfp_iexpt (p x n)
"Compute (X^N mod P)
Uses the right-to-left binary algorithm (also known as 'exponentiation by squaring' and 'binary exponentiation'). Also uses a trick to shrink exponents
mod (p-1) -- making it faster exponents much larger than p, but a tiny bit slower for smaller exponents.
References:
Bruce Schneier (1996); Applied Cryptography: Protocols, Algorithms, and Source Code in C 2nd; ISBN: 978-0471117094; pp224
David Bressoud (1989); Factorization and primality testing; ISBN: 0-387-97040-1; pp33-34"
(cond ((not (integerp n)) (error "mjr_gfp_imul: n must be an integer!")))
(let ((result +mjr_gfp_one+)
(x (mjr_gfp_simplify p x))
(n (mjr_gfp_simplify (1- p) n))) ;; This will always work because p is prime...
(loop while (not (zerop n))
do (if (logbitp 0 n)
(setq result (mjr_gfp_simplify p (* result x))))
do (setq n (ash n -1)
x (mjr_gfp_simplify p (* x x))))
result))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_gfp_imul (p x n)
"Compute (X*N mod P) -- for GF(p) this is easy, for other rings it must be implemented as repeated addition.."
(cond ((not (integerp n)) (error "mjr_gfp_imul: n must be an integer!")))
(mjr_gfp_* p x n))