mersenneforum.org  

Go Back   mersenneforum.org > Extra Stuff > Programming

Reply
 
Thread Tools
Old 2017-12-20, 23:59   #1
carpetpool
 
carpetpool's Avatar
 
"Sam"
Nov 2016

1001111002 Posts
Post Pari/GP to PFGW

Hi,

The following code is to generate prime indices of prime Fibonacci numbers in GP in a certain index range.

Code:
v=[3, 4]; forprime(p=5, 1e5, if(ispseudoprime(fibonacci(p)), v=concat(v, p)));
I find this code very useful for when the fibonacci(p) sequence is replaced by another divisiblity sequence a(p) (defined in GP).

Code:
v=[]; forprime(p=1, 1e5, if(ispseudoprime(a(p)), v=concat(v, p)));
which will list the primes p such that a(p) is prime for (p given in a certain range).

Now the bad thing or catch about this is PARI's pseudoprime(n) function is seemingly slow if n is a (considerably large, say 1000+ digits) pseudoprime (probable prime). Since PARI/GP seems to be fast at generating the terms for a(p), (which PFGW cannot do for all sequences and high-order recurrence relations), can someone come up with a modification to fiabonacci(p) code so that PARI/GP generates the prime-index Fibonacci numbers, and PFGW 3-PRP tests them. (I am aware of PFGW's built-in fibonacci function, however I don't want that being used here because if PARI is able to generate the fibonacci terms, it will also be able to for other sequences, which PFGW can easily PRP test.

Thanks for help.


The following GP code is from here.
carpetpool is offline   Reply With Quote
Old 2017-12-21, 00:17   #2
science_man_88
 
science_man_88's Avatar
 
"Forget I exist"
Jul 2009
Dumbassville

8,369 Posts
Default

Quote:
Originally Posted by carpetpool View Post
Hi,

The following code is to generate prime indices of prime Fibonacci numbers in GP in a certain index range.

Code:
v=[3, 4]; forprime(p=5, 1e5, if(ispseudoprime(fibonacci(p)), v=concat(v, p)));
I find this code very useful for when the fibonacci(p) sequence is replaced by another divisiblity sequence a(p) (defined in GP).

Code:
v=[]; forprime(p=1, 1e5, if(ispseudoprime(a(p)), v=concat(v, p)));
which will list the primes p such that a(p) is prime for (p given in a certain range).

Now the bad thing or catch about this is PARI's pseudoprime(n) function is seemingly slow if n is a (considerably large, say 1000+ digits) pseudoprime (probable prime). Since PARI/GP seems to be fast at generating the terms for a(p), (which PFGW cannot do for all sequences and high-order recurrence relations), can someone come up with a modification to fiabonacci(p) code so that PARI/GP generates the prime-index Fibonacci numbers, and PFGW 3-PRP tests them. (I am aware of PFGW's built-in fibonacci function, however I don't want that being used here because if PARI is able to generate the fibonacci terms, it will also be able to for other sequences, which PFGW can easily PRP test.

Thanks for help.


The following GP code is from here.
technically, the new PARI 2.10.0 alpha ( or the newest one I have) has a fibonacci command. maybe test other conditions to weed things out a bit first ? programming efficiently relates to testing conditions that are faster first and only going for the big ones later. maybe try changing the flag to change the test being used behind the scenes ? edit: two edits to your first code made a small speed up ( changed to parforprime and called the fibonacci outside the ispseudoprime test part in parallel).

Last fiddled with by science_man_88 on 2017-12-21 at 01:14
science_man_88 is offline   Reply With Quote
Old 2017-12-21, 00:58   #3
CRGreathouse
 
CRGreathouse's Avatar
 
Aug 2006

72×112 Posts
Default

PFGW is definitely faster, but the comparison is not a fair one: it's just doing a 3-PRP test, while gp is doing a full BPSW test (~3x slower). A better comparison would be to use ispseudoprime(..., 1) which will perform 1 Miller-Rabin test with a random base.

I call PFGW within gp like so:
Code:
pfgw(n:int)=
{
	my(prog="~/mth/pfgw/pfgw64 -k -N -u0 -q");
	Vecsmall(externstr(Str(prog, n)));
}
and you may also find the extern and system commands useful.
CRGreathouse is offline   Reply With Quote
Old 2017-12-21, 01:50   #4
danaj
 
"Dana Jacobsen"
Feb 2011
Bangkok, TH

2·3·151 Posts
Default

Admittedly the BPSW test should exit early for the composites. Still, PFGW will be faster on the Fermat / M-R once the result hits somewhere between 2000 and 6000 digits -- at least that's where it's faster than GMP's mpz_powm.


As an aside, I have an example Fibonacci prime finder in my Perl module that has a couple parallel implementations. Starting at the beginning searching for each, my timing results are

26244 seconds for Fp36 (p = 148091) on a 3930K (12 proc)
12000 seconds for Fp36 (p = 148091) on a c4.8xlarge (36 proc)
6465 seconds for Fp36 (p = 148091) on a r4.16xlarge (64 proc)

Someday I'll get a big c5 for a day to benchmark.

[Edit: that it not the timing for one test, but the time from starting search for all Fp until we see the output of the in-order 36th Fibonacci prime]

I've tried it on a Power8 with 152 logical processors but it's a shared machine so really couldn't let it run very long, and it looks like GMP is better optimized for x86-64.

It's pretty simple, just running primality tests in order on each core, with the most interesting part just being making sure the results are shown in order.

Last fiddled with by danaj on 2017-12-21 at 01:51
danaj is offline   Reply With Quote
Old 2017-12-21, 02:15   #5
danaj
 
"Dana Jacobsen"
Feb 2011
Bangkok, TH

2×3×151 Posts
Default

Thanks for the Pari/GP function to use PFGW. That's handy.

Connecting to another thread on Lucas sequences, I saw the nice OEISWiki page. Are there some good Pari/GP functions for lucasu(P,Q,n) and lucasv(P,Q,n) possibly with n = Mod(x,y)? Good in the sense of either (or both) (1) clean, (2) fast?

It looks like maybe polcoef could do it. Not sure about efficiently handling the modulo case though.

Pari/GP's fibonacci() is pretty basic. Does Mod(fibonacci(n),p) optimize? It seems unlikely.
danaj is offline   Reply With Quote
Old 2017-12-21, 02:39   #6
science_man_88
 
science_man_88's Avatar
 
"Forget I exist"
Jul 2009
Dumbassville

100000101100012 Posts
Default

Quote:
Originally Posted by danaj View Post
Are there some good Pari/GP functions for lucasu(P,Q,n) and lucasv(P,Q,n) possibly with n = Mod(x,y)? Good in the sense of either (or both) (1) clean, (2) fast?
well a basic parity argument might be able to be used mod 2 for example. if both P and Q are the same parity that results in producing certain patterns in parity in the sequence. at least if I understand this request correctly.
science_man_88 is offline   Reply With Quote
Old 2017-12-21, 06:04   #7
CRGreathouse
 
CRGreathouse's Avatar
 
Aug 2006

10111001010012 Posts
Default

Quote:
Originally Posted by danaj View Post
Pari/GP's fibonacci() is pretty basic. Does Mod(fibonacci(n),p) optimize? It seems unlikely.
No, and it couldn't -- order of operations forces fibonacci(n) to be computed in full before reducing it mod p.

I use
Code:
fibmod(n, m)=((Mod([1, 1; 1, 0], m))^n)[1, 2]
which works decently. I have a version in C, nothing special, which is faster.

Quote:
Originally Posted by danaj View Post
Connecting to another thread on Lucas sequences, I saw the nice OEISWiki page. Are there some good Pari/GP functions for lucasu(P,Q,n) and lucasv(P,Q,n) possibly with n = Mod(x,y)? Good in the sense of either (or both) (1) clean, (2) fast?

It looks like maybe polcoef could do it. Not sure about efficiently handling the modulo case though.
I'll look into it tomorrow. I don't think that signature will work, though; just because you're working mod m doesn't mean you can reduce the index mod m, right? So I think you'll need separate arguments for those two.
CRGreathouse is offline   Reply With Quote
Reply

Thread Tools


Similar Threads
Thread Thread Starter Forum Replies Last Post
LLL in GP/Pari paul0 Programming 2 2015-11-17 13:04
PARI vs GAP skan Miscellaneous Math 0 2012-12-16 00:13
pari devarajkandadai Programming 21 2012-08-31 18:08
PFGW 3.3.6 or PFGW 3.4.2 Please update now! Joe O Sierpinski/Riesel Base 5 5 2010-09-30 14:07
64-bit Pari? CRGreathouse Software 2 2009-03-13 04:22

All times are UTC. The time now is 09:17.

Fri Oct 23 09:17:20 UTC 2020 up 43 days, 6:28, 0 users, load averages: 1.23, 1.58, 1.70

Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2020, Jelsoft Enterprises Ltd.

This forum has received and complied with 0 (zero) government requests for information.

Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.2 or any later version published by the Free Software Foundation.
A copy of the license is included in the FAQ.