Eric <ericbelcastro@gmail.com> wrote:
+---------------
| Say I want to create a function that fills a list of length n with a
| sequence that counts from 1 to n, at a given interval r and returns
| the list. And this process being generalized for any given length or
| interval. Of course the length and interval have to be relatively prime.
| For example: If the interval is 3 and the length is 17, the list
| returned would be...
| (1 7 13 2 8 14 3 9 15 4 10 16 5 11 17 6 12)
| You can see that, since the interval is 3, if you read
| every 3 in the final list you will see 1 2 3 4...etc.
+---------------
Spending a little time on theory up front often makes the resulting
program simpler in the end. First, let's do a small coordinate
transformation on your requirement into a domain that's a little
easier to deal with, specifially, seeking the desired permutation
of the numbers 0..N-1 rather than 1..N. [We can do a simple "fixup"
at the end.] Counting from zero is usually easier to deal with in such
matters, see <http://www.cs.utexas.edu/users/EWD/ewd08xx/EWD831.PDF>.
I use my version of IOTA a lot [doesn't everyone have one in
their toolbox]:
> (defun iota (count &optional (start 0) (step 1))
(loop repeat count for i upfrom start by step collect i))
IOTA
> (iota 17)
(0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16)
> (iota 17 1)
(1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17)
> (iota 17 1 3)
(1 4 7 10 13 16 19 22 25 28 31 34 37 40 43 46 49)
>
Let's rename & extend this into an IOTA with a STRIDE [what you called
INTERVAL] modulo LENGTH using Ken Tilton's temporary array method
[shown in a parallel response]:
> (defun iota/mod/stride (length &key (start 0) (step 1) (stride 1))
(loop with temp-array = (make-array length)
repeat length
for val = start then (+ val step)
and index = 0 then (mod (+ index stride) length)
do (setf (svref temp-array index) val)
finally (return (coerce temp-array 'list))))
IOTA/MOD/STRIDE
> (iota/mod/stride 17 :stride 3 :start 1)
(1 7 13 2 8 14 3 9 15 4 10 16 5 11 17 6 12)
> (iota/mod/stride 27 :stride 5 :start 1)
(1 12 23 7 18 2 13 24 8 19 3 14 25 9 20 4 15 26 10 21 5 16 27 11 22 6 17)
>
That's what we want, yes? Now let's see how to get rid of that ugly
temporary array... First, let's go back to counting from zero:
(iota/mod/stride 17 :stride 3)
(0 6 12 1 7 13 2 8 14 3 9 15 4 10 16 5 11)
>
Observe that now the array INDEX is always (mod (* VAL STRIDE) LENGTH),
and *that* means that at the instant we store into array index 1, VAL
is the multiplicative inverse of STRIDE (mod LENGTH), by the definition
of multiplicative inverses:
X * X^(-1) = 1 (mod N)
So to get rid of the array, all we have to do is to find the
multiplicative inverse of STRIDE modulo LENGTH, and then we
simply generate the desired elements in order, this way:
> (defun iota/mod/stride (length &key (start 0) (stride 1))
(loop with stride^-1 = (reciprocal/mod stride length)
repeat length
for x = 0 then (mod (+ x stride^-1) length)
collect (+ x start)))
IOTA/MOD/STRIDE
> (iota/mod/stride 17 :stride 3 :start 1)
(1 7 13 2 8 14 3 9 15 4 10 16 5 11 17 6 12)
> (iota/mod/stride 27 :stride 5 :start 1)
(1 12 23 7 18 2 13 24 8 19 3 14 25 9 20 4 15 26 10 21 5 16 27 11 22 6 17)
>
Whoops! But where did that RECIPROCAL/MOD function come from?
That is, how do we find the multiplicative inverse of an integer
modulo another integer? After some helpful sidebar discussions
with Kenny, and rummaging around in Knuth's & Schneier's books,
it hit me. (D'Oh!) "Use the Web, Luke!":
http://en.wikipedia.org/wiki/Modular_multiplicative_inverse
This gives two methods, the extended Euclidean algorithm and
direct modular exponentiation using Euler's totient function.
Since the latter is a bit messy, let's look at the former:
http://en.wikipedia.org/wiki/Extended_Euclidean_algorithm
The simplest[1] version of this to code in Lisp is the recursive
method, which we can directly transliterate from the pseudocode
in the article like this:
> (defun extended-gcd (a b)
(if (zerop (mod a b))
(values 0 1)
(multiple-value-bind (x y)
(extended-gcd b (mod a b))
(values y (- x (* y (floor a b)))))))
EXTENDED-GCD
> (extended-gcd 3 17)
6
-1
> (extended-gcd 5 27)
11
-2
>
The primary value in each case is the desired reciprocol, that is,
3^-1 (mod 17) = 6, and 5^-1 (mod 27) = 11. Now all we need is a
simple wrapper that verifies that the LENGTH and STRIDE are co-prime:
> (defun reciprocal/mod (r n)
(multiple-value-bind (x y)
(extended-gcd r n)
(let ((gcd (+ (* r x) (* n y))))
(assert (= 1 gcd) () "~s and ~s are not co-prime!" r n)
(if (minusp x)
(+ x n) ; not necessary, but esthetically pleasing
x))))
RECIPROCAL/MOD
> (reciprocal/mod 3 17)
6
> (reciprocal/mod 5 27)
11
> (reciprocal/mod 9 27)
Error in function LISP::ASSERT-ERROR: 9 and 27 are not co-prime!
[Condition of type SIMPLE-ERROR]
...
With those two helpful functions, IOTA/MOD/STRIDE will work as shown.
QED.
[Note that in CMUCL, when all the above functions are compiled,
no bytes are consed for R, N less than ~10000, and only 48 bytes
[due to the "(+ (* r x) (* n y)" term in RECIPROCAL/MOD] for
R, N that are large fixnums. So the above will definitely cons
less than the temporary array version.]
-Rob
[1] I said, "The simplest method in Lisp is the recursive method",
but the iterative method might be faster. Exercise for the reader:
Recode EXTENDED-GCD using the iterative method and compare with
the recursive method for execution time & bytes consed.
-----
Rob Warnock <rpw3@rpw3.org>
627 26th Avenue <URL:http://rpw3.org/>
San Mateo, CA 94403 (650)572-2607