20100917, 11:23  #1 
Mar 2006
2^{2}×137 Posts 
Need help with Lucas Sequences...
Hello everyone,
I'm trying to write some code that will calculate the value of the kth 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 pseudocode 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_{0} = 2 V_{1} = P V_{k}(P,Q) = P * V_{k1}(P,Q)  Q * V_{k2}(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 ) Code:
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*p2)%n; i; while (i) { if (e[i]) // if ith 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; } Code:
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 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}  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_{5}(1,2) % 349 = 7 ) If anyone can confirm or deny this, I'd appreciate it. 
20100917, 13:53  #2  
"William"
May 2003
Near Grandkid
5^{3}·19 Posts 
Quote:
V_{0}(1,2) = 2 V_{1}(1,2) = 1 V_{2}(1,2) = 1 * V_{1}(1,2)  2 * V_{0}(1,2) = 1*12*2 = 3 V_{3}(1,2) = 1 * V_{2}(1,2)  2 * V_{1}(1,2) = 1*(3)2*1 = 5 V_{4}(1,2) = 1 * V_{3}(1,2)  2 * V_{2}(1,2) = 1*(5)2*(3) = 1 V_{5}(1,2) = 1 * V_{4}(1,2)  2 * V_{3}(1,2) = 1*(1)2*(5) = 11 

20100917, 15:01  #3 
"Forget I exist"
Jul 2009
Dartmouth NS
2×5^{2}×13^{2} Posts 
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. 
20100917, 15:55  #4  
Mar 2006
2^{2}·137 Posts 
Quote:
Quote:
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. 

20100917, 16:01  #5  
Aug 2006
2^{2}×3×499 Posts 
Quote:
http://pari.math.ubordeaux.fr/cgib...2541&root=pari It has reasonable performance: the 208,988digit F_1,000,000 is calculated in 78 milliseconds on a slow PC. If you need a speed demon, though, this isn't it. 

20100917, 16:21  #6  
"Forget I exist"
Jul 2009
Dartmouth NS
10000100000010_{2} Posts 
Quote:


20100917, 16:30  #7  
"Forget I exist"
Jul 2009
Dartmouth NS
2×5^{2}×13^{2} Posts 
Quote:


20100918, 01:32  #8  
Mar 2006
1000100100_{2} Posts 
Quote:
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); } 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 kth value of a generic Lucas V Sequence, I'd love to hear about it. 

20100918, 14:27  #9 
Mar 2006
2^{2}·137 Posts 
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. Code:
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^24q=0 is not tested, no proper Lucas sequence!! */ /* Returns (nonnegative) 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^24q = 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 */ 
20100920, 00:53  #10  
"Gang aft agley"
Sep 2002
2×1,877 Posts 
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 Quote:
Quote:


20100923, 22:04  #11 
"Forget I exist"
Jul 2009
Dartmouth NS
20402_{8} Posts 
I couldn't find this thread earlier when i wanted it lol.

Thread Tools  
Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Silly Lucaslike sequences of primes  fivemack  Math  14  20170205 05:44 
LucasLehmer sequences starting with s0 other than 4  GP2  Math  6  20160521 17:00 
Any interest in all sequences/open sequences?  Greebley  Aliquot Sequences  6  20120407 10:06 
factors of lucas sequences  jocelynl  Math  2  20100923 21:32 
LucasLehmer  Dougal  Information & Answers  9  20090206 10:25 