Mathematica has a function called PowerMod that's incredibly useful in
just about all of the currently interesting cryptographic protocols:
The modular power function PowerMod[a, b, n] gives exactly the same
results as Mod[a^b, n] for \!\(TraditionalForm\`b > 0\). PowerMod is
much more efficient, however, because it avoids generating the full
form of a^b.
You can use PowerMod not only to find positive modular powers, but
also to find modular inverses. For negative b, PowerMod[a, b, n]
gives, if possible, an integer \!\(TraditionalForm\`k\) such that
\!\(TraditionalForm\`k a\^\(-b\) \[Congruent] 1 mod n\). (Whenever
such an integer exists, it is guaranteed to be unique modulo
\!\(TraditionalForm\`n\).) If no such integer
\!\(TraditionalForm\`k\) exists, Mathematica leaves PowerMod
unevaluated.
Has anyone (who's able and willing to share) implemented this in
Common Lisp?
-- Kaelin
Kaelin Colclasure wrote:
> Mathematica has a function called PowerMod that's incredibly useful in
> just about all of the currently interesting cryptographic protocols:
>
> The modular power function PowerMod[a, b, n] gives exactly the same
> results as Mod[a^b, n] for \!\(TraditionalForm\`b > 0\). PowerMod is
> much more efficient, however, because it avoids generating the full
> form of a^b.
>
> You can use PowerMod not only to find positive modular powers, but
> also to find modular inverses. For negative b, PowerMod[a, b, n]
> gives, if possible, an integer \!\(TraditionalForm\`k\) such that
> \!\(TraditionalForm\`k a\^\(-b\) \[Congruent] 1 mod n\). (Whenever
> such an integer exists, it is guaranteed to be unique modulo
> \!\(TraditionalForm\`n\).) If no such integer
> \!\(TraditionalForm\`k\) exists, Mathematica leaves PowerMod
> unevaluated.
>
> Has anyone (who's able and willing to share) implemented this in
> Common Lisp?
Yes sure,
;;; Fast (mod (expt number exponent) modulus)
(defun expt-mod (number exponent modulus)
(loop with result = 1
for exp = exponent then (floor exp 2)
for sqr = number then (mod (* sqr sqr) modulus)
until (zerop exp)
when (oddp exp) do
(setf result (mod (* result sqr) modulus))
finally (return result)))
Nothing special so far...
Here's some more...
;;; if d = gcd(a,b) then there is a x and b y so that d = ax + by
(defun extended-gcd (a b)
(loop for x = 1 then x2
and y = 0 then y2
and d = a then d2
and x2 = 0 then (- x (* k x2))
and y2 = 1 then (- y (* k y2))
and d2 = b then (- d (* k d2))
until (zerop d2)
for k = (floor d d2)
finally (return (values d x y))))
;;; a^-1 mod n , provided that it exists
(defun multiplicative-inverse-mod (a n)
(multiple-value-bind (d x y)
(extended-gcd a n)
(if (> d 1)
(error "INVERSE-MOD: no inverse~
=> modulus is not prime.")
(loop while (minusp x) do
(incf x n)
finally (return x)))))
I've used it for implementing a miller-rabin pseudo-prime test and an
implementation of the RSA public key cipher.
If I find the time I will place this code (and the rest like miller-rabin
and RSA) on my Webpage http://www.dataheaven.de.
Meanwhile permission is hereby granted to use the above code if you include
the following copyright statement and disclaimer. (I know it's not very
much code but it's part of some more code that I worked some time to get it
working as it does)
;;;; ;
;;;; (c) 2001 by Jochen Schmidt.
;;;;
;;;; 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.
;;;;
;;;; THIS SOFTWARE IS PROVIDED "AS IS" AND THERE ARE NEITHER
;;;; EXPRESSED NOR IMPLIED WARRANTIES - THIS INCLUDES, BUT
;;;; IS NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY
;;;; AND FITNESS FOR A PARTICULAR PURPOSE.IN NO WAY ARE THE
;;;; AUTHORS 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)
;;;;
;;;; For further details contact the authors of this software.
;;;;
;;;; Jochen Schmidt
;;;; Zuckmantelstr. 11
;;;; 91616 Neusitz
;;;; GERMANY
;;;;
;;;; N�rnberg, 16.Apr.2001 Jochen Schmidt
;;;;
Regards,
Jochen Schmidt
Jochen Schmidt <···@dataheaven.de> writes:
> Yes sure,
>
> ;;; Fast (mod (expt number exponent) modulus)
> (defun expt-mod (number exponent modulus)
> (loop with result = 1
> for exp = exponent then (floor exp 2)
> for sqr = number then (mod (* sqr sqr) modulus)
> until (zerop exp)
> when (oddp exp) do
> (setf result (mod (* result sqr) modulus))
> finally (return result)))
>
> Nothing special so far...
>
> Here's some more...
>
> ;;; if d = gcd(a,b) then there is a x and b y so that d = ax + by
> (defun extended-gcd (a b)
> (loop for x = 1 then x2
> and y = 0 then y2
> and d = a then d2
> and x2 = 0 then (- x (* k x2))
> and y2 = 1 then (- y (* k y2))
> and d2 = b then (- d (* k d2))
> until (zerop d2)
> for k = (floor d d2)
> finally (return (values d x y))))
>
> ;;; a^-1 mod n , provided that it exists
> (defun multiplicative-inverse-mod (a n)
> (multiple-value-bind (d x y)
> (extended-gcd a n)
> (if (> d 1)
> (error "INVERSE-MOD: no inverse~
> => modulus is not prime.")
> (loop while (minusp x) do
> (incf x n)
> finally (return x)))))
Perfect! Well, almost... ;-) I might suggest that the error message
above is slightly misleading to the mathematically unsophisticated --
a group for whom I speak with some confidence. Rather than ``modulus is
not prime'' I would suggest ``<a> is not primitive modulo <n>''.
Thank you! And I do hope you find the time to post the rest of this code.
I've added your site to my bookmarks and will keep an eye out for it. :-)
> I've used it for implementing a miller-rabin pseudo-prime test and an
> implementation of the RSA public key cipher.
> If I find the time I will place this code (and the rest like miller-rabin
> and RSA) on my Webpage http://www.dataheaven.de.
> Meanwhile permission is hereby granted to use the above code if you include
> the following copyright statement and disclaimer. (I know it's not very
> much code but it's part of some more code that I worked some time to get it
> working as it does)
>
> ;;;; ;
> ;;;; (c) 2001 by Jochen Schmidt.
> ;;;;
> ;;;; 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.
> ;;;;
> ;;;; THIS SOFTWARE IS PROVIDED "AS IS" AND THERE ARE NEITHER
> ;;;; EXPRESSED NOR IMPLIED WARRANTIES - THIS INCLUDES, BUT
> ;;;; IS NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY
> ;;;; AND FITNESS FOR A PARTICULAR PURPOSE.IN NO WAY ARE THE
> ;;;; AUTHORS 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)
> ;;;;
> ;;;; For further details contact the authors of this software.
> ;;;;
> ;;;; Jochen Schmidt
> ;;;; Zuckmantelstr. 11
> ;;;; 91616 Neusitz
> ;;;; GERMANY
> ;;;;
> ;;;; N�rnberg, 16.Apr.2001 Jochen Schmidt
> ;;;;
>
> Regards,
> Jochen Schmidt
From: F. Xavier Noria
Subject: Re: CL implementation of Mathematica's PowerMod?
Date:
Message-ID: <3adac518.678328@news.iddeo.es>
On Mon, 16 Apr 2001 01:27:25 +0200, Jochen Schmidt <···@dataheaven.de> wrote:
Just some comments (I am not a Lisp programmer, so comments
on my comments or code will be very welcomed!)
: ;;; Fast (mod (expt number exponent) modulus)
: (defun expt-mod (number exponent modulus)
: (loop with result = 1
: for exp = exponent then (floor exp 2)
: for sqr = number then (mod (* sqr sqr) modulus)
: until (zerop exp)
: when (oddp exp) do
: (setf result (mod (* result sqr) modulus))
: finally (return result)))
There is a way to write this using bit operators, if I can
write it I'll send it. Also, for almost the same price we can
allow negative exponents.
: ;;; if d = gcd(a,b) then there is a x and b y so that d = ax + by
: (defun extended-gcd (a b)
: (loop for x = 1 then x2
: and y = 0 then y2
: and d = a then d2
: and x2 = 0 then (- x (* k x2))
: and y2 = 1 then (- y (* k y2))
: and d2 = b then (- d (* k d2))
: until (zerop d2)
: for k = (floor d d2)
: finally (return (values d x y))))
There is no need to compute both x2 and y2 along the loops,
provided x2 and the gcd you can obtain y2.
Also, since the Euclidean Algorithm is naturally expressed
using recursion I implement a function that uses an inner
tail-recursive routine:
(defun euclidean-div (a b)
"Given integers a, b, being b nonzero, (euclidean-div a b) returns
the unique integers q, r, such that a = bq + r and 0 <= r < |b|."
;; a/b = q + r/b
(if (< 0 b)
(floor a b)
(ceiling a b)))
;;; Extended Euclidean Algorithm.
(defun bezout (a b)
"Given nonzero integers a, b, (bezout a b) returns integers x, y, d,
such that d = gcd(a, b) and ax + by = d."
;; x[-2] = 1, x[-1] = 0, x[i] = x[i-2] - q[i]x[i-1]
;; r[-2] = a, r[-1] = b, r[i-2] = r[i-1]q[i] + r[i]
(labels ((rec (x[i-2] x[i-1] r[i-2] r[i-1])
(multiple-value-bind (q[i] r[i]) (euclidean-div r[i-2] r[i-1])
(if (zerop r[i])
(values x[i-1] (/ (- r[i-1] (* a x[i-1])) b) r[i-1])
(rec x[i-1] (- x[i-2] (* q[i] x[i-1])) r[i-1] r[i])))))
(rec 1 0 a b)))
: ;;; a^-1 mod n , provided that it exists
: (defun multiplicative-inverse-mod (a n)
: (multiple-value-bind (d x y)
: (extended-gcd a n)
: (if (> d 1)
: (error "INVERSE-MOD: no inverse~
: => modulus is not prime.")
From d > 1 doesn't follow that the modulus is not prime. For
instance 2 has no inverse mod 6 (and 6 is not prime). An integer
n is invertible mod m iff gcd(n, m) is 1, so the error message
could say they are not coprime.
: (loop while (minusp x) do
: (incf x n)
: finally (return x)))))
Now you have no need to compute y2 at all, which is done inside
EXTENDED-GCD. An alternative:
;;; Extended Euclidean Algorithm.
(defun inv-mod (g m)
"Given integers g, m, being m > 0, if [g] is invertible in Z/Zm
(inv-mod g m) returns a representant of its inverse class, otherwise
returns nil."
;; x[-2] = 1, x[-1] = 0, x[i] = x[i-2] - q[i]x[i-1]
;; r[-2] = a, r[-1] = b, r[i-2] = r[i-1]q[i] + r[i]
(labels ((rec (x[i-2] x[i-1] r[i-2] r[i-1])
(multiple-value-bind (q[i] r[i]) (euclidean-div r[i-2] r[i-1])
(if (zerop r[i])
(if (= r[i-1] 1)
(mod x[i-1] m))
(rec x[i-1] (- x[i-2] (* q[i] x[i-1])) r[i-1] r[i])))))
(rec 1 0 g m)))
-- fxn
F. Xavier Noria wrote:
> On Mon, 16 Apr 2001 01:27:25 +0200, Jochen Schmidt <···@dataheaven.de>
> wrote:
>
> Just some comments (I am not a Lisp programmer, so comments
> on my comments or code will be very welcomed!)
>
> : ;;; Fast (mod (expt number exponent) modulus)
> : (defun expt-mod (number exponent modulus)
> : (loop with result = 1
> : for exp = exponent then (floor exp 2)
> : for sqr = number then (mod (* sqr sqr) modulus)
> : until (zerop exp)
> : when (oddp exp) do
> : (setf result (mod (* result sqr) modulus))
> : finally (return result)))
>
> There is a way to write this using bit operators, if I can
> write it I'll send it. Also, for almost the same price we can
> allow negative exponents.
I've already done that (bit-operator-stuff) and it doesn't get you much.
I've done some tests in the past and it fully depends on the Lispsystems
implementation of FLOOR or LOGBITP/INTEGER-LENGTH (for the bit-operator
variant).
The most costly part of the whole thing are the bignum multiplications.
It is maybe possible to implement an own bignum-multiplication but I like
readable and understandabe code until the performance is really a problem.
Here the same function using bit-operators:
(defun expt-mod (n exponent modulus)
"Same as (mod (expt n exponent) modulus) just more efficient."
(declare (optimize (speed 3) (safety 0)))
(loop with result = 1
for i of-type fixnum from 0 below (integer-length exponent)
for sqr = n then (mod (* sqr sqr) modulus)
when (logbitp i exponent) do
(setf result (mod (* result sqr) modulus))
finally (return result)))
There's no measurable performance difference between this and the former
variant (I tested in the past)
the "handbook of applied cryptography" describes some more specialized
variants of modular exponentiation.
> : ;;; if d = gcd(a,b) then there is a x and b y so that d = ax + by
> : (defun extended-gcd (a b)
> : (loop for x = 1 then x2
> : and y = 0 then y2
> : and d = a then d2
> : and x2 = 0 then (- x (* k x2))
> : and y2 = 1 then (- y (* k y2))
> : and d2 = b then (- d (* k d2))
> : until (zerop d2)
> : for k = (floor d d2)
> : finally (return (values d x y))))
>
> There is no need to compute both x2 and y2 along the loops,
> provided x2 and the gcd you can obtain y2.
I have to admit that I have not looked deeply into this function it is more
of a 1:1 conversion from the pseudo-code of the "Handbook of Applied
Cryptography"
> Also, since the Euclidean Algorithm is naturally expressed
> using recursion I implement a function that uses an inner
> tail-recursive routine:
<snipped some code>
> : ;;; a^-1 mod n , provided that it exists
> : (defun multiplicative-inverse-mod (a n)
> : (multiple-value-bind (d x y)
> : (extended-gcd a n)
> : (if (> d 1)
> : (error "INVERSE-MOD: no inverse~
> : => modulus is not prime.")
>
> From d > 1 doesn't follow that the modulus is not prime. For
> instance 2 has no inverse mod 6 (and 6 is not prime). An integer
> n is invertible mod m iff gcd(n, m) is 1, so the error message
> could say they are not coprime.
Yes your're right. The message is nonsense. If I remember right it was used
in some cryptopackage sources I looked into.
>
> : (loop while (minusp x) do
> : (incf x n)
> : finally (return x)))))
>
> Now you have no need to compute y2 at all, which is done inside
> EXTENDED-GCD. An alternative:
Yes that's true too - but I thought the performance is good enough to allow
the generalization aspect of using the more general extended-gcd.
> ;;; Extended Euclidean Algorithm.
> (defun inv-mod (g m)
> "Given integers g, m, being m > 0, if [g] is invertible in Z/Zm
> (inv-mod g m) returns a representant of its inverse class, otherwise
> returns nil."
> ;; x[-2] = 1, x[-1] = 0, x[i] = x[i-2] - q[i]x[i-1]
> ;; r[-2] = a, r[-1] = b, r[i-2] = r[i-1]q[i] + r[i]
> (labels ((rec (x[i-2] x[i-1] r[i-2] r[i-1])
> (multiple-value-bind (q[i] r[i]) (euclidean-div r[i-2]
> r[i-1])
> (if (zerop r[i])
> (if (= r[i-1] 1)
> (mod x[i-1] m))
> (rec x[i-1] (- x[i-2] (* q[i] x[i-1])) r[i-1] r[i])))))
> (rec 1 0 g m)))
Looks nice I will test it if I find time.
You may take a look into the complete package available on
http://www.dataheaven.de
there's an implementation of the Miller-Rabin Pseudoprimetest and
the RSA algorithm.
Regards,
Jochen