Need help with Lucas Sequences...

[WraithX]

Hello everyone,

I'm trying to write some code that will calculate the value of the k-th step of the Lucas V sequence modulo some number. (ie, I'm trying to find out V_k(P,Q) modulo n ) I've searched online for some pseudo-code and haven't found anything I can really use.

I can see how to calculate it in O(k) time just by using the recurrence relation:
V₀ = 2
V₁ = P
V_k(P,Q) = P * V_k-1(P,Q) - Q * V_k-2(P,Q), k > 1

However, for some problems I'll be working on, this will be way too slow. I've found mention of a way to calculate it in O(lg(k)) time by using Lucas Chains. I've even found some code that was close to this, but it isn't completely generic. It assumes that Q=1, and I need it for any Q. I've tried to modify it, but I can't seem to get it right. Can someone help me figure this out?

Here is the code I found online that is close:
(this calculate V_e(p,1) modulo n )

Integer Lucas(const Integer &e, const Integer &p, const Integer &n)
{
        unsigned i = e.BitCount();

        if (i==0)
                return 2;

        Integer v=p, v1=(p*p-2)%n;
        i--;
        while (i--)
        {
                if (e[i])       // if i-th bit of e is 1
                {
                        v = (v*v1 - p) % n;
                        v1 = (v1*v1 - 2) % n;
                }
                else
                {
                        v1 = (v*v1 - p) % n;
                        v = (v*v - 2) % n;
                }
        }
        return v;
}

I've tried to modify it to include Q like so:

Lucas (k, p, q, n)
  i = bitcount(k)
  if i == 0 return 2
  v = p
  qn = q
  v1 = (p*p - 2*qn) % n
  qn = (qn*qn) % n
  i--
  while (i--) {
    if k[i] {
      v  = (v*v1 - p*qn) % n
      v1 = (v1*v1 - 2*qn) % n
      qn = (qn*qn*q) % n
    }
    else {
      v1 = (v*v1 - p*qn) % n
      v  = (v*v - 2*qn) % n
      qn = (qn*qn) % n
    }
  }
  return v

However, I think this is giving me wrong answers. Can anyone see what I might have done wrong, or does anyone just happen to have the correct pseudo-code lying around? ;)

There are a couple of recurrence relations that allow the O(lg(k)) speedup to be realized, and they are:
V_2k(P,Q) = V_k(P,Q)² - 2*Q^k
V_2k+1(P,Q) = V_k+1(P,Q) * V_k(P,Q) - P*Q^k

Also, to make sure I'm calculating some of this correctly I've tried to work an example with k=5, p=1, q=2, n=349. From the recurrence relation, I believe the answer should be 7 (ie, V₅(1,2) % 349 = 7 ) If anyone can confirm or deny this, I'd appreciate it.

[wblipp]

Originally Posted by WraithX

V₀ = 2
V₁ = P
V_k(P,Q) = P * V_k-1(P,Q) - Q * V_k-2(P,Q), k > 1

Also, to make sure I'm calculating some of this correctly I've tried to work an example with k=5, p=1, q=2, n=349. From the recurrence relation, I believe the answer should be 7 (ie, V₅(1,2) % 349 = 7 ) If anyone can confirm or deny this, I'd appreciate it.

assuming p = P and q = Q, I get

V₀(1,2) = 2
V₁(1,2) = 1
V₂(1,2) = 1 * V₁(1,2) - 2 * V₀(1,2) = 1*1-2*2 = -3
V₃(1,2) = 1 * V₂(1,2) - 2 * V₁(1,2) = 1*(-3)-2*1 = -5
V₄(1,2) = 1 * V₃(1,2) - 2 * V₂(1,2) = 1*(-5)-2*(-3) = 1
V₅(1,2) = 1 * V₄(1,2) - 2 * V₃(1,2) = 1*(1)-2*(-5) = 11

[science_man_88]

http://en.wikipedia.org/wiki/Lucas_sequence

might help in one sense if i knew anything useful I could try to make Pari Code and get someone like CRG to translate into C.

[WraithX]

Originally Posted by wblipp

assuming p = P and q = Q, I get

V₀(1,2) = 2
V₁(1,2) = 1
V₂(1,2) = 1 * V₁(1,2) - 2 * V₀(1,2) = 1*1-2*2 = -3
V₃(1,2) = 1 * V₂(1,2) - 2 * V₁(1,2) = 1*(-3)-2*1 = -5
V₄(1,2) = 1 * V₃(1,2) - 2 * V₂(1,2) = 1*(-5)-2*(-3) = 1
V₅(1,2) = 1 * V₄(1,2) - 2 * V₃(1,2) = 1*(1)-2*(-5) = 11

Thank you for this, somehow I (mistakenly) got -3 in that last place to come up with 7.

Originally Posted by science_man_88

http://en.wikipedia.org/wiki/Lucas_sequence

might help in one sense if i knew anything useful I could try to make Pari Code and get someone like CRG to translate into C.

Yeah, I've looked there but they don't provide enough info for me to write a Lucas Chain version of a Lucas Sequence solver.

If you stumble upon a web site that does give pseudo code (or pari or C) I'd love to hear about it. I can take an actual algorithm and code it up how I need it. Thank you for looking, hopefully we can find an implementation on the internet or find out how to fix my implementation up above.

[CRGreathouse]

Originally Posted by WraithX

Thank you for looking, hopefully we can find an implementation on the internet or find out how to fix my implementation up above.

There's basic code here:
http://pari.math.u-bordeaux.fr/cgi-b...2541&root=pari

It has reasonable performance: the 208,988-digit F_1,000,000 is calculated in 78 milliseconds on a slow PC. If you need a speed demon, though, this isn't it.

[science_man_88]

Originally Posted by http://www.newsletter-library.com/Science/Math/Number_Theory/Prime_Numbers/Primality_Tests

Lucas Sequences in Cryptography
Focus on their use in primality testing, with pseudo-code. Includes an explanation of the strong Lucas PRP test.
site exerpt
Lucas Sequences in Cryptography This is a short note on the practical usefulness of Lucas sequences in applied cryptography. A Lucas sequence is a sequence of integers characterized by two parameters, P and Q. In practice Q is always...
http://www.eskimo.com/~weidai/lucas.html

does this fit the bill ?

[science_man_88]

Originally Posted by http://www.numericana.com/answer/factoring.htm

We can't resist making a few incidental remarks (inspired by D.E. Knuth):
A fairly useless expression for the index n thus computed is:
n = l max ( 1, é m/ l ù )

Once n is known, we observe that:

m is the smallest i such that Un+i = Ui
l is the smallest i such that Un+i = Un

We need not check both of these conditions beyond the point where the first one is met (which gives directly the preperiod m): If the second one has already been met, we're done with the period (l) as well. Otherwise, we have m < l and must have l = n (since l divides n < m+l < 2l). All told, we may obtain m and l by calling the function f only 2m additional times, using the following pseudocode (the above has "initialized" n and y = Un ).

i := 0 ; x := U0 ; z := y ; l := 0
while x ¹ z
i := i+1 ; x := f (x) ; z := f (z)
if z = y and l = 0 then l := i
wend
m := i
if l = 0 then l := n

I think this looks closer lol

[WraithX]

Originally Posted by CRGreathouse

There's basic code here:
http://pari.math.u-bordeaux.fr/cgi-b...2541&root=pari

It has reasonable performance: the 208,988-digit F_1,000,000 is calculated in 78 milliseconds on a slow PC. If you need a speed demon, though, this isn't it.

Hmmm, looking at that page, the only reference I see to Lucas is the following code:

/*******************************************************************/
/**                                                               **/
/**                      LUCAS & FIBONACCI                        **/
/**                                                               **/
/*******************************************************************/
static void
lucas(ulong n, GEN *a, GEN *b)
{
  GEN z, t, zt;
  if (!n) { *a = gen_2; *b = gen_1; return; }
  lucas(n >> 1, &z, &t); zt = mulii(z, t);
  switch(n & 3) {
    case  0: *a = addsi(-2,sqri(z)); *b = addsi(-1,zt); break;
    case  1: *a = addsi(-1,zt);      *b = addsi(2,sqri(t)); break;
    case  2: *a = addsi(2,sqri(z));  *b = addsi(1,zt); break;
    case  3: *a = addsi(1,zt);       *b = addsi(-2,sqri(t));
  }
}

GEN
fibo(long n)
{
  pari_sp av = avma;
  GEN a, b;
  if (!n) return gen_0;
  lucas((ulong)(labs(n)-1), &a, &b);
  a = diviuexact(addii(shifti(a,1),b), 5);
  if (n < 0 && !odd(n)) setsigne(a, -1);
  return gerepileuptoint(av, a);
}

This only looks like it calculates Fibonacci and Lucas numbers, which are special cases of the generic Lucas Sequence. With F(n) = U_n(1,-1) and L(n) = V_n(1,-1).

I'm going to keep searching and/or trying to correct my original code from the first post, but if anyone knows (where I can find) (pseudo-) code to generate the k-th value of a generic Lucas V Sequence, I'd love to hear about it.

[WraithX]

Ok, you'll all be happy to know I now have a working implementation of a generic Lucas V Sequence. I found some reference code at:
http://home.netsurf.de/wolfgang.ehrhardt/mp_intro.html
He seems to have a ton of math functions programmed in Pascal in his own math library. I was able to adapt his:

mp_lucasvmod Calculate v[k] mod n of Lucas V sequence for p,q

into a C function. I'm including it below in case anyone may need this in the future.

int lucasvmod(int p, int q, int n, int k)
{
  /* adapted from the math package at http://home.netsurf.de/wolfgang.ehrhardt/mp_intro.html */
  /* calculate v[k] mod n of Lucas V sequence for p,q. */
  /* Ranges n>1, 0 <= p,q,k < n */
  /* Note: p^2-4q=0 is not tested, no proper Lucas sequence!! */
  /* Returns (non-negative) value of V_k(p,q) mod n, returns -1 if there was an error. */

  int i;
  int v0,v1,x,t,mu;

  if (n < 2) {printf("lucasvmod: n<2\n"); return -1;}
  if (p < 0 || p >= n) {printf("lucasvmod: p<0 or p>=n\n"); return -1;}
  if (q < 0 || q >= n) {printf("lucasvmod: q<0 or q>=n\n"); return -1;}
  if (k < 0 || k >= n) {printf("lucasvmod: k<0 or k>=n\n"); return -1;}
  if ((p*p - 4*q) == 0) {printf("lucasvmod: p^2-4q = 0 **error**\n"); return -1;}

  /* {v[k+2] = p*v[k+1] - q*v[k], v[0]=2, v[1]=p} */

  /* {easy out for k<=0, k<0 is handled as k=0} */
  if (k < 1) return 2;

  /* {v0 = p} */
  v0 = p%n;

  /* {start with x=q} */
  x = q%n;

  v1 = (p*p - 2*q)%n;

  /* {loop is never executed for k=1} */
  /* set i to the number of bits in k, ie k=5 gives i=3 */
  for (i = log_base_2(k)-2; i >= 0; i--)
  {
    if (k&(1<<i))
    {
      v0 = v0*v1 - x*p;
      v0 = v0%n; if (v0 < 0) v0 += n;
      v1 = v1*v1 - 2*x*q;
      v1 = v1%n; if (v1 < 0) v1 += n;
      if (i != 0) x = (x*x*q)%n;
    }
    else
    {
      v1 = v0*v1 - x*p;
      v1 = v1%n; if (v1 < 0) v1 += n;
      v0 = v0*v0 - 2*x;
      v0 = v0%n; if (v0 < 0) v0 += n;
      if (i != 0) x = (x*x)%n;
    }
  }
  return v0;
}/* method lucasvmod */

[only_human]

This short paper looks useful. It also computes the value based on the binary expansion of the index.

http://joye.site88.net/papers/JQ96lucas.pdf

Published in Electronics Letters 32(6):537–538, 1996.
Efficient computation of full Lucas sequences
M. Joye and J.-J. Quisquater
Indexing terms: Number theory, Cryptography
Recently, Yen and Laih (1995) proposed an algorithm to compute
LUC digital signatures quickly. This signature is based on a special
type of the Lucas sequence, V_k. The authors generalise their
method to any type of Lucas sequence, and extend it to the ‘sister’
Lucas sequence, Uk. As an application, the order of an elliptic
curve over GF(2m) is computed quickly.

I also notice that your example assigns an error value to negative indices. In skimming around I wonder if it is perhaps useful to extend the sequence to the left.On Fibonacci Numbers & Mathematical Phyllotaxis

A number of wonderful mathematical properties of the " extended" series F_n
and L_n follow from Table 2.1. For example:
n= 0 1 2 3 4 5 6 7
F_n 0 1 1 2 3 5 8 13
F_-n 0 1 -1 2 -3 5 -8 13
L_n 2 1 3 4 7 11 18 29
L_-n 2 -1 3 -4 7 -11 18 -29

Theorem 2.2.1 for odd integers n = 2k + 1 the terms of the sequences F_n and
F_-n coincide, i.e., F_2k+1 = F_-2k-1, while for even integers n = 2k they have op-
posite signs, i.e., F_2k = -F_-2k, whereas the reverse occurs for the Lucas numbers
L_n, with L_2k = L_-2k and L_2k+1 = -L_-2k-1

[science_man_88]

I couldn't find this thread earlier when i wanted it lol.