mersenneforum.org

mersenneforum.org (https://www.mersenneforum.org/index.php)
-   PARI/GP (https://www.mersenneforum.org/forumdisplay.php?f=155)
-   -   PARI's commands (https://www.mersenneforum.org/showthread.php?t=13636)

CRGreathouse 2011-05-24 18:49

Faster version:
[code]pick(n,k)=binomial(if(n<k,n+k-1,n),k)[/code]

Basically, your function is binomial(n+k-1, k) if n<k and binomial(n,k) otherwise. I could have written
[code]pick(n,k)=if(n<k,binomial(n+k-1,k), binomial(n,k))[/code]
but I chose to combine the similar parts and use just
[code]pick(n,k)=binomial([COLOR="Red"]if(n<k,n+k-1,n)[/COLOR],k)[/code]

When PARI can compute the binomial directly instead of indirectly through the factorial it can save a lot of steps. When you calculate [TEX]{9\choose2}=\frac{9!}{7!2!}=\frac{9\cdot8\cdot7\cdot6\cdot5\cdot4\cdot3\cdot2\cdot1}{(7\cdot6\cdot5\cdot4\cdot3\cdot2\cdot1)(2\cdot1)}[/TEX]
you probably just cross out the 7 * ... * 1 on top and bottom. When you give PARI binomial(9,2) it can do the same; when you give it 9!/(7!*2!) it's forced to compute the large numbers and then divide.

Xitami 2011-05-24 20:42

[code]binom(n,k)={ my(t=i=1);
if(2*k>n, k=n-k);
while(i<=k,
t=(t*n)/i;
i++; n--
);
t
}[/code]-----------------------------

science_man_88 2011-05-24 22:17

[QUOTE=CRGreathouse;262186]Faster version:
[code]pick(n,k)=binomial(if(n<k,n+k-1,n),k)[/code]

Basically, your function is binomial(n+k-1, k) if n<k and binomial(n,k) otherwise. I could have written
[code]pick(n,k)=if(n<k,binomial(n+k-1,k), binomial(n,k))[/code]
but I chose to combine the similar parts and use just
[code]pick(n,k)=binomial([COLOR="Red"]if(n<k,n+k-1,n)[/COLOR],k)[/code]

When PARI can compute the binomial directly instead of indirectly through the factorial it can save a lot of steps. When you calculate [TEX]{9\choose2}=\frac{9!}{7!2!}=\frac{9\cdot8\cdot7\cdot6\cdot5\cdot4\cdot3\cdot2\cdot1}{(7\cdot6\cdot5\cdot4\cdot3\cdot2\cdot1)(2\cdot1)}[/TEX]
you probably just cross out the 7 * ... * 1 on top and bottom. When you give PARI binomial(9,2) it can do the same; when you give it 9!/(7!*2!) it's forced to compute the large numbers and then divide.[/QUOTE]

thanks:[CODE](19:14)>allocatemem(1)
*** allocatemem: Warning: new stack size = 1024 (0.001 Mbytes).
(19:15)>pick(11,1000000000)
%269 = 275573207396384843474431809413623724923088088349707046408843005955895496034675000001
(19:15)>allocatemem(931987423)
*** allocatemem: Warning: new stack size = 931987420 (888.812 Mbytes).
(19:15)>pick1(11,1000000000)
*** _!: the PARI stack overflows !
current stack size: 931987420 (888.812 Mbytes)
[hint] you can increase GP stack with allocatemem()[/CODE]

I changed my original to pick1 for this test.

science_man_88 2011-05-24 23:27

found something I find interesting now CRG:

[CODE](20:20)>for(a=1,38,print1(sum(x=1,a,pick(Me(x)%6,Me(x+1)%6)%4)%4","))
[COLOR="Red"]0,1,2,3,[/COLOR][COLOR="DarkOrange"]0,1,2,3,[/COLOR][COLOR="Lime"]0,1,2,3,[/COLOR][COLOR="Red"]0,1,2,3,[/COLOR][COLOR="DarkOrange"]0,1,2,3,[/COLOR][COLOR="Lime"]0,1,2,3,[/COLOR][COLOR="Red"]0,1,2,3,[/COLOR][COLOR="DarkOrange"]0,1,2,3,[/COLOR][COLOR="Lime"]0,1,2,3,[/COLOR][COLOR="Red"]0,1,[/COLOR][/CODE]

okay it works out this way for all primes it looks like.

science_man_88 2011-05-25 00:48

[QUOTE=science_man_88;262207]found something I find interesting now CRG:

[CODE](20:20)>for(a=1,38,print1(sum(x=1,a,pick(Me(x)%6,Me(x+1)%6)%4)%4","))
[COLOR="Red"]0,1,2,3,[/COLOR][COLOR="DarkOrange"]0,1,2,3,[/COLOR][COLOR="Lime"]0,1,2,3,[/COLOR][COLOR="Red"]0,1,2,3,[/COLOR][COLOR="DarkOrange"]0,1,2,3,[/COLOR][COLOR="Lime"]0,1,2,3,[/COLOR][COLOR="Red"]0,1,2,3,[/COLOR][COLOR="DarkOrange"]0,1,2,3,[/COLOR][COLOR="Lime"]0,1,2,3,[/COLOR][COLOR="Red"]0,1,[/COLOR][/CODE]

okay it works out this way for all primes it looks like.[/QUOTE]

I made an extended version of MeVec to test this and so far the exponents for x=40 to 47( or guessing) fit this pattern.

CRGreathouse 2011-05-25 05:07

What happens if you reverse the order of the last 8?

science_man_88 2011-05-25 11:37

[QUOTE=CRGreathouse;262228]What happens if you reverse the order of the last 8?[/QUOTE]

nothing I bet because late last night I realized why the pattern is set up it's because pick(a,b) with a = 1 or 5, b=1,5 will be 1 mod 4. so the sum is mostly a sum of ones ( which will only cycle mod anything) the 2 that start it 2 and 3 and 3 and 5 are 0 and 1 mod 4 so 0+1 is the start and +1 for each step after. so really all this states is that the first 3 exponents are used to set it up.

science_man_88 2011-05-25 22:57

[QUOTE=science_man_88;262207]found something I find interesting now CRG:

[CODE](20:20)>for(a=1,38,print1(sum(x=1,a,pick(Me(x)%6,Me(x+1)%6)%4)%4","))
[COLOR="Red"]0,1,2,3,[/COLOR][COLOR="DarkOrange"]0,1,2,3,[/COLOR][COLOR="Lime"]0,1,2,3,[/COLOR][COLOR="Red"]0,1,2,3,[/COLOR][COLOR="DarkOrange"]0,1,2,3,[/COLOR][COLOR="Lime"]0,1,2,3,[/COLOR][COLOR="Red"]0,1,2,3,[/COLOR][COLOR="DarkOrange"]0,1,2,3,[/COLOR][COLOR="Lime"]0,1,2,3,[/COLOR][COLOR="Red"]0,1,[/COLOR][/CODE]

okay it works out this way for all primes it looks like.[/QUOTE]

I've got another pattern to figure out before revealing.

science_man_88 2011-05-25 23:11

[QUOTE=science_man_88;262323]I've got another pattern to figure out before revealing.[/QUOTE]

okay I gave in:

[CODE](19:34)>for(x=1,47,print1(pick(prime(x),x^2-10)%x","))
0,0,0,3,0,0,0,0,0,4,0,6,0,0,0,0,0,0,0,0,0,0,0,0,0,13,0,20,0,0,0,0,0,0,30,0,0,0,0,0,0,0,0,0,0,0,0,
(19:34)>for(x=1,47,print1(pick(prime(x),x^2-5)%x","))
0,0,2,0,0,2,0,0,0,4,0,9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,16,0,0,25,0,0,0,13,0,0,28,0,0,0,0,0,
(19:35)>for(x=1,47,print1(pick(MPE[x],x^2-10)%x","))
0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,15,0,0,0,0,0,0,0,0,0,0,
(19:35)>for(x=1,47,print1(pick(MPE[x],x^2-5)%x","))
0,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
(19:35)>for(x=1,47,print1(MPE[x]<x^2-5","))
0,0,0,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
(19:47)>for(x=1,47,print1(MPE[x]<x^2-10","))
0,0,0,0,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,[/CODE]

[QUOTE]MPE[x]!/(k!*(MPE[x]-k)!)
(MPE[x]+k-1)!/((k!)*(MPE[x]-1)!)[/QUOTE]

is what I consider important

science_man_88 2011-05-26 17:45

[QUOTE=science_man_88;262324]okay I gave in:

[CODE][COLOR="Red"](19:34)>for(x=1,47,print1(pick(prime(x),x^2-10)%x","))
0,0,0,3,0,0,0,0,0,4,0,6,0,0,0,0,0,0,0,0,0,0,0,0,0,13,0,20,0,0,0,0,0,0,30,0,0,0,0,0,0,0,0,0,0,0,0,
(19:34)>for(x=1,47,print1(pick(prime(x),x^2-5)%x","))
0,0,2,0,0,2,0,0,0,4,0,9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,16,0,0,25,0,0,0,13,0,0,28,0,0,0,0,0,[/COLOR]
[COLOR="Lime"](19:35)>for(x=1,47,print1(pick(MPE[x],x^2-10)%x","))
0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,15,0,0,0,0,0,0,0,0,0,0,
(19:35)>for(x=1,47,print1(pick(MPE[x],x^2-5)%x","))
0,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,[/COLOR]
[COLOR="DarkOrange"](19:35)>for(x=1,47,print1(MPE[x]<x^2-5","))
0,0,0,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
(19:47)>for(x=1,47,print1(MPE[x]<x^2-10","))
0,0,0,0,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,[/COLOR][/CODE]



is what I consider important[/QUOTE]


the red is for proving not all primes follow the patterns 1 after the other.

the green shows the patterns I think I have ( pick(Me(x),x^2-5) is 0 mod 2 and pick(Me(x),x^2-10) is 0 mod 3).

the orange is to show which path of pick each Me(x) brings us through.

science_man_88 2011-06-30 15:11

sieve of atkins
 
not sure how to create this , I know I ignore 2,3,and 5, and I check mod 60 but I have no idea ( the pseudo-code in the Wikipedia article didn't help me as far as I can tell) otherwise.


All times are UTC. The time now is 22:58.

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