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-02-05 02:13

I don't get that, or any, error. (The program doesn't seem to do what you want, but that's another matter,)

science_man_88 2011-02-05 12:05

[QUOTE=CRGreathouse;251356]I don't get that, or any, error. (The program doesn't seem to do what you want, but that's another matter,)[/QUOTE]

[CODE]if(#A!=0&&#B!=0,a=0;for(i=1,#A,for(j=1,#B,if(A[i]==B[j],return(0),a=a+1;if(a==#A*#B,return(1))))),return(0))[/CODE] works perfect when put into PARI as far as I can tell but gives and error when in your form in file.

science_man_88 2011-02-05 12:09

I find something weird if I had it in fancy form I was always "missing something" if not it worked fine.

CRGreathouse 2011-02-05 19:37

Maybe you didn't use braces, so it read each line individually?

I always write my functions like
[code]foo=[COLOR="red"]{[/COLOR]
\\ . . .
[COLOR="Red"]}[/COLOR];[/code]
so I know that problem won't happen to me.

science_man_88 2011-02-05 19:38

[QUOTE=CRGreathouse;251439]Maybe you didn't use braces, so it read each line individually?

I always write my functions like
[code]foo=[COLOR="red"]{[/COLOR]
\\ . . .
[COLOR="Red"]}[/COLOR];[/code]
so I know that problem won't happen to me.[/QUOTE]

I never knew they were mandatory.

CRGreathouse 2011-02-05 22:41

[QUOTE=science_man_88;251440]I never knew they were mandatory.[/QUOTE]

If you don't use them then each line will be interpreted individually. When you type a command and press enter, it does it right away instead of waiting for another line. So you 'knew', really, that it worked that way -- you just never quite associated the two.

But really, do write your functions out in full form. You often seem to get frustrated with others when they don't respond to your posts, but that's one big reason they don't -- it's too hard to figure out what's going on in your scripts. When you make it easy for other people to read they're more likely to read them and hence respond.

science_man_88 2011-02-05 23:00

[QUOTE=CRGreathouse;251454]If you don't use them then each line will be interpreted individually. When you type a command and press enter, it does it right away instead of waiting for another line. So you 'knew', really, that it worked that way -- you just never quite associated the two.

But really, do write your functions out in full form. You often seem to get frustrated with others when they don't respond to your posts, but that's one big reason they don't -- it's too hard to figure out what's going on in your scripts. When you make it easy for other people to read they're more likely to read them and hence respond.[/QUOTE]
okay here's the ones I've got coded for the homework stuff I learned:

[CODE]aredisjoint(A=[],B=[]) ={
a=0;
for(i=1,#A,
for(j=1,#B,
if(A[i]==B[j],
return(0),
a=a+1;
if(a==#A*#B,
return(1),
return(0)
)
)
)
)};
addhelp(aredisjoint,"a function to tell when 2 sets (vectors(which are actually sequences)) are considered disjoint.")
ispd(C) ={
b=0;
for(i=1,#C,
for(j=1,#C,
if(i!=j && aredisjoint(C[i],C[j]),
b=b+1
)
)
);
if(b==#C*#C-#C,
return(1),
return(0)
)
};
addhelp(ispd,"a function to tell when a collection of sets (vectors(which are actually sequences)) is considered pairwise disjoint.")
CPC(C)={
Cp(C[1],C[2]);
A=d;
for(i=3,#C,
Cp(A,C[i]);
A=d
);
printVectorAsSet(A) };
addhelp(CPC,"a function to create the cartesian product of a collection of sets.")
Cp(p,c)={
d=[];
for(i=1,#p,
for(j=1,#c,
x=[p[i],c[j]];
d=concat(d,[x])
)
);
d;};
addhelp(Cp,"a function to create the cartesian product of 2 sets.")
printVectorAsSet(v)={
if(type(v) != "t_VEC",
print1(v)
,
print1("{");
if (#v, printVectorAsSet(v[1]));
for(i=2,#v,
print1(", ");
printVectorAsSet(v[i])
);
print1("}")
)
}
intersection(x,y) ={
C=[];
for(i=1,#x,
for(j=1,#y,
if(x[i]==y[j],
C=concat(C,x[i])
)
)
);
printVectorAsSet(vecsort(C,,8)); };[/CODE]

Any errors ? if not I may continue and make a union script assuming I remember what a union is.

CRGreathouse 2011-02-05 23:34

aredisjoint([1,2,3],[4,5]) returns 0.

science_man_88 2011-02-05 23:36

[QUOTE=CRGreathouse;251457]aredisjoint([1,2,3],[4,5]) returns 0.[/QUOTE]

okay I added a new thing to that and on my end it works it's and if loop checking if #a and #b are = 0 but that shouldn't make a difference in that:

[CODE]aredisjoint(A=[],B=[]) ={
if(#A!=0 && #B!=0,
a=0;
for(i=1,#A,
for(j=1,#B,
if(A[i]==B[j],
return(0),
a=a+1;
if(a==#A*#B,
return(1),
)
)
)
)
,return(0)
)};[/CODE]

science_man_88 2011-02-06 00:30

[QUOTE=science_man_88;251456]if not I may continue and make a union script assuming I remember what a union is.[/QUOTE]

I looked over the other thread and I realized it's simple:

[CODE]union(A,B)=concat(A,B);vecsort(union,,8)[/CODE]

so simple I messed it up:

[CODE]union(A,B)= U=concat(A,B);vecsort(U,,8)[/CODE]

science_man_88 2011-02-06 01:33

Just to update I've changed one thing in Cp and fixed the output for CPC.

hopefully I can program more of this it may become very useful.

science_man_88 2011-02-06 01:43

[CODE]interCol(C) ={ IC=intersection(C[1],C[2]);for(i=1,#C,IC=intersection(IC,C[i]));IC;}[/CODE] problem is intersection makes it pretty which then makes it useless here. I think I'll have to add a control to pretty it or not.

CRGreathouse 2011-02-06 02:18

[QUOTE=science_man_88;251463]I looked over the other thread and I realized it's simple:

[CODE]union(A,B)=concat(A,B);vecsort(union,,8)[/CODE]

so simple I messed it up:

[CODE]union(A,B)= U=concat(A,B);vecsort(U,,8)[/CODE][/QUOTE]

Yes... though you should *only* use variables that have been declared with my() or local() -- and use my() unless you have a good reason to use local(). So
[code]union(A,B)=my(U=concat(A,B));vecsort(U,,8)[/code]
Of course in this case you don't even need a variable at all -- can you see how to avoid it?

CRGreathouse 2011-02-06 02:21

[QUOTE=science_man_88;251458]okay I added a new thing to that and on my end it works it's and if loop checking if #a and #b are = 0 but that shouldn't make a difference in that:

[CODE]aredisjoint(A=[],B=[]) ={
if(#A!=0 && #B!=0,
a=0;
for(i=1,#A,
for(j=1,#B,
if(A[i]==B[j],
return(0),
a=a+1;
if(a==#A*#B,
return(1),
)
)
)
)
,return(0)
)};[/CODE][/QUOTE]

Still wrong.

The check that neither are empty is not only needless but actually harmful here. Also, there's no need to count the number of times you find disjoint pairs -- if you find any identical elements they're not disjoint, otherwise they are.

CRGreathouse 2011-02-06 02:23

[QUOTE=science_man_88;251469]Just to update I've changed one thing in Cp and fixed the output for CPC.

hopefully I can program more of this it may become very useful.[/QUOTE]

When you get all of them working maybe you can post them all in one swell foop. :grin:

[QUOTE=science_man_88;251470][CODE]interCol(C) ={ IC=intersection(C[1],C[2]);for(i=1,#C,IC=intersection(IC,C[i]));IC;}[/CODE] problem is intersection makes it pretty which then makes it useless here. I think I'll have to add a control to pretty it or not.[/QUOTE]

I don't know what this is trying to do.

science_man_88 2011-02-06 12:19

[QUOTE=CRGreathouse;251474]Yes... though you should *only* use variables that have been declared with my() or local() -- and use my() unless you have a good reason to use local(). So
[code]union(A,B)=my(U=concat(A,B));vecsort(U,,8)[/code]
Of course in this case you don't even need a variable at all -- can you see how to avoid it?[/QUOTE]

[CODE]vecsort(concat(A,B),,8)[/CODE]

CRGreathouse 2011-02-06 18:24

Yep, you win.

As for the disjoint program: If you have a working intersection program you can just use that. If the sets are disjoint, what (by definition!) will their intersection be? That makes for a very short (albeit potentially slow) disjoint function.

science_man_88 2011-02-06 22:05

[QUOTE=CRGreathouse;251555]Yep, you win.

As for the disjoint program: If you have a working intersection program you can just use that. If the sets are disjoint, what (by definition!) will their intersection be? That makes for a very short (albeit potentially slow) disjoint function.[/QUOTE]

yeah well I tried this:

[CODE](17:57)>aredisjoint(v,c)
%9 = 1
(17:57)>##
*** last result computed in 0 ms.
(17:57)>aredisjoint2(v,c)
%10 = 1
(17:57)>##
*** last result computed in 328 ms.
(17:57)>?aredisjoint2
aredisjoint2(A,B)=if(intersection(A,B)==[],return(1),return(0))[/CODE]

the aredisjoint is:

[CODE]aredisjoint(A=[],B=[]) ={
a=0;
for(i=1,#A,
for(j=1,#B,
if(A[i]==B[j],
return(0),
return(1)
)
)
)};[/CODE]

the reason the aredisjoint2 is slower is because my intersection script doesn't check if they have a intersection but goes through both sets creating a new vector for their intersection. my disjoint script aredisjoint only has to check the whole of both sets if no exception is found. never mind I'm an idiot.

CRGreathouse 2011-02-06 22:37

[QUOTE=science_man_88;251586]the reason the aredisjoint2 is slower is because my intersection script doesn't check if they have a intersection but goes through both sets creating a new vector for their intersection. my disjoint script aredisjoint only has to check the whole of both sets if no exception is found. never mind I'm an idiot.[/QUOTE]

Certainly you can check whether two sets are disjoint faster than you can compute their intersection in most cases. But your current code is broken and I was suggesting a way to fix it.

The whole idea of breaking programs into different functions is that you can change the internal workings of one without affecting the others. Generally I recommend making a function that gives the correct answer, even if it runs very slowly, and using that to check if your newer, fast versions are correct.


[QUOTE=science_man_88;251586][CODE]aredisjoint(A=[],B=[]) ={
a=0;
for(i=1,#A,
for(j=1,#B,
if(A[i]==B[j],
return(0),
return(1)
)
)
)};[/CODE][/QUOTE]

This function is equivalent to
[code]aredisjoint(A=[],B=[]) ={
a=0;
if(#A & #B, A[1] != B[1])
};[/code]

That is, the function does this:
1. It clobbers the global a variable. That is,
[code]a=1;aredisjoint();print(a)[/code]
prints 0 instead of 1, a bad side-effect. (This is why you should always use my() or local().)
2. It checks if either vector is of length 0. If so, it does nothing -- internally it returns a value called gnil which is interpreted as 0 if needed but which is not displayed in gp.
3. Otherwise, it looks at A[1] and B[1]. If the values are the same it returns 0, otherwise it returns 1. It never looks at any other values.

science_man_88 2011-02-06 23:39

[QUOTE=CRGreathouse;251590]Certainly you can check whether two sets are disjoint faster than you can compute their intersection in most cases. But your current code is broken and I was suggesting a way to fix it.

The whole idea of breaking programs into different functions is that you can change the internal workings of one without affecting the others. Generally I recommend making a function that gives the correct answer, even if it runs very slowly, and using that to check if your newer, fast versions are correct.




This function is equivalent to
[code]aredisjoint(A=[],B=[]) ={
a=0;
if(#A & #B, A[1] != B[1])
};[/code]

That is, the function does this:
1. It clobbers the global a variable. That is,
[code]a=1;aredisjoint();print(a)[/code]
prints 0 instead of 1, a bad side-effect. (This is why you should always use my() or local().)
2. It checks if either vector is of length 0. If so, it does nothing -- internally it returns a value called gnil which is interpreted as 0 if needed but which is not displayed in gp.
3. Otherwise, it looks at A[1] and B[1]. If the values are the same it returns 0, otherwise it returns 1. It never looks at any other values.[/QUOTE]

I suck at this:

[CODE](19:37)>aredisjoint([],[])
%1 = 0
(19:38)>aredisjoint([],[1,2,3])
%2 = 0
(19:38)>aredisjoint([4,5,6],[1,2,3])
%3 = 1
(19:38)>aredisjoint([3,5,6],[1,2,3])
%4 = 0
(19:38)>[/CODE]

[CODE]aredisjoint(A=[],B=[]) ={
if(#A!=0 [COLOR="Red"]&&[/COLOR] #B!=0,
for(i=1,#A,
for(j=1,#B,
if(A[i]==B[j],
return(0)
)
)
);return(1)
,return(0))
}[/CODE]

I think this could be changed to || because [] is a subset of every other set and hence no set can be disjoint from it so if it appears even once we know it's 0.

CRGreathouse 2011-02-07 00:36

You should remove that entire if() test on the size of the vectors, it doesn't need to be there.

science_man_88 2011-02-07 12:11

[QUOTE=CRGreathouse;251598]You should remove that entire if() test on the size of the vectors, it doesn't need to be there.[/QUOTE]

well it will return 1 for aredisjoint([],[]) it would also do that if only one is empty as we know the empty set is not disjoint from any set.

CRGreathouse 2011-02-07 15:46

[QUOTE=science_man_88;251651]well it will return 1 for aredisjoint([],[]) it would also do that if only one is empty as we know the empty set is not disjoint from any set.[/QUOTE]

What's the problem?

science_man_88 2011-02-07 15:57

[QUOTE=CRGreathouse;251663]What's the problem?[/QUOTE]

if we take away the if statement on the outside then it tries to find the first element of a 0 element set since that's impossible it would skip to the end and return 1 regardless if they are the same or the fact that [] is a subset of every set.

CRGreathouse 2011-02-07 16:16

[QUOTE=science_man_88;251665]if we take away the if statement on the outside then it tries to find the first element of a 0 element set[/QUOTE]

No it doesn't. How many times does for(i=1,#A, ...) execute?

science_man_88 2011-02-07 16:32

[QUOTE=CRGreathouse;251667]No it doesn't. How many times does for(i=1,#A, ...) execute?[/QUOTE]

that doesn't matter because the check I have inside checks with indexes if it's not an index it won't be checked which means if i take away the if it can only return 1.

axn 2011-02-07 16:42

[QUOTE=science_man_88;251665]if we take away the if statement on the outside then it tries to find the first element of a 0 element set since that's impossible it would skip to the end and return 1 regardless if they are the same or the fact that [] is a subset of every set.[/QUOTE]

So? If one or both of the sets are empty, they _are_ disjoint. Disjoint means "no common element". If one of them is empty, how can there be a common element?

science_man_88 2011-02-07 16:49

[QUOTE=axn;251670]So? If one or both of the sets are empty, they _are_ disjoint. Disjoint means "no common element". If one of them is empty, how can there be a common element?[/QUOTE]

I was told the empty is part of every set unless that's false [],[] can only return 0 because they have something in common [],set or set,[] have the same reason to return 0 but my code can't work with 0 elements unless I have that outer statement if the statement isn't there it returns 1 (true) that all these are disjoint but everything I've been told makes that a bad result.

axn 2011-02-07 17:32

[QUOTE=science_man_88;251671]I was told the empty is part of every set [/quote]
Empty set is a _subset_ of every set. It is not an _element_ of every set. "Part of" is not a mathematically meaningful term.

[QUOTE=science_man_88;251671]unless that's false [],[] can only return 0 because they have something in common [],set or set,[] have the same reason to return 0 but my code can't work with 0 elements unless I have that outer statement if the statement isn't there it returns 1 (true) that all these are disjoint but everything I've been told makes that a bad result.[/QUOTE]

Please restate the definition of disjoint sets before proceeding any further. You've hit a conceptual "wrong turn" here.

CRGreathouse 2011-02-07 20:47

This is the empty set (as we're representing it with PARI vectors): []

This is a set containing the empty set as a member: [[]]

This is another set containing the empty set as a member: [1, 2, 3, 4, 5, [], 6]

aredisjoint([], [[]]) should return 1 because they are disjoint: the first set has no members and so there are no common members.

aredisjoint([1, 2, 3, 4, 5, [], 6], [[]]) should return 0 because they aren't disjoint: [] is a member of both.

science_man_88 2011-02-08 13:18

[QUOTE=CRGreathouse;251708]This is the empty set (as we're representing it with PARI vectors): []

This is a set containing the empty set as a member: [[]]

This is another set containing the empty set as a member: [1, 2, 3, 4, 5, [], 6]

aredisjoint([], [[]]) should return 1 because they are disjoint: the first set has no members and so there are no common members.

aredisjoint([1, 2, 3, 4, 5, [], 6], [[]]) should return 0 because they aren't disjoint: [] is a member of both.[/QUOTE]

point taken.

science_man_88 2011-02-10 15:35

ispower
 
[CODE](11:29)>?ispower
ispower(x,{k},{&n}): if k > 0 is given, return true (1) if x is a k-th power, false (0) if not. If k is omitted, return the maximal k >= 2 such that x =
n^k is a perfect power, or 0 if no such k exist. If n is present, and the function returns a non-zero result, set n to the k-th root of x.

(11:31)>ispower(4,,2)
*** expected character: '&': ispower(4,,2)
^--
(11:31)>ispower(4,,&2)
*** syntax error, unexpected KINTEGER, expecting KENTRY: ispower(4,,&2)
^--
(11:31)>ispower(4,,{&2})
*** syntax error, unexpected KINTEGER, expecting KENTRY: ispower(4,,&2)
^--
(11:31)>ispower(4,,{&2})
*** syntax error, unexpected KINTEGER, expecting KENTRY: ispower(4,,&2)
^--
(11:32)>ispower(4,,{2&})
*** syntax error, unexpected ')': ispower(4,,2&)
^-
(11:32)>ispower(4,,&2)
*** syntax error, unexpected KINTEGER, expecting KENTRY: ispower(4,,&2)
^--
(11:32)>ispower(4,2,&2)
*** syntax error, unexpected KINTEGER, expecting KENTRY: ispower(4,2,&2)
^--
(11:32)>ispower(4,,&2)
*** syntax error, unexpected KINTEGER, expecting KENTRY: ispower(4,,&2)
^--
(11:32)>ispower(4,,2)
*** expected character: '&': ispower(4,,2)
^--
(11:32)>ispower(4,,&)
*** syntax error, unexpected ')', expecting KENTRY: ispower(4,,&)
^-
(11:33)>ispower(4,,&"2")
*** syntax error, unexpected KSTRING, expecting KENTRY: ispower(4,,&"2")
^--
(11:33)>ispower(4,,&2)
*** syntax error, unexpected KINTEGER, expecting KENTRY: ispower(4,,&2)
^--
(11:33)>ispower(4,,{2})
*** expected character: '&': ispower(4,,2)
^--
(11:33)>ispower(4,,&{2})
*** syntax error, unexpected KINTEGER, expecting KENTRY: ispower(4,,&2)
^--
(11:33)>ispower(4,,&})
*** unexpected closing brace.
(11:33)>ispower(4,,&}
*** unexpected closing brace.
(11:33)>ispower(4,,&2}
*** unexpected closing brace.
(11:33)>ispower(4,,&2)
*** syntax error, unexpected KINTEGER, expecting KENTRY: ispower(4,,&2)
^--
(11:34)>ispower(4,,{2})
*** expected character: '&': ispower(4,,2)
^--
(11:34)>ispower(4,,2)
*** expected character: '&': ispower(4,,2)
^--
(11:34)>ispower(4,,2)[/CODE]

axn 2011-02-10 18:57

[CODE]? n=0
%12 = 0
? ispower(1024,,&n)
%13 = 10
? n
%14 = 2[/CODE]

You must pass a variable as the third parameter so that the function can pass back the "base".

science_man_88 2011-02-10 19:11

[QUOTE=axn;252089][CODE]? n=0
%12 = 0
? ispower(1024,,&n)
%13 = 10
? n
%14 = 2[/CODE]

You must pass a variable as the third parameter so that the function can pass back the "base".[/QUOTE]

thanks for the tip I never quite figured it out I asked because i was seeing something that I sent to prime95 as weird, now I don't as much but if you look (2^mersenne1[n]-1)%n seems to be a mersenne number a lot if mersenne1 is the mersenne prime exponents. in fact it looks like 35/40 have this property in the first 100 prime numbers if you do 2^(prime(n))-1 over there index into the primes it seems about 49% of the first 100 work out.

CRGreathouse 2011-02-10 21:37

ispower(n) checks if n is a nontrivial power. ispower(n,k) checks if n is a k-th power. ispower(n,,&b) checks if n is a nontrivial power and, if so, stores the base in b. ispower(n,k,&b) checks if n is a k-th power and, if so, stores the base in b.

science_man_88 2011-02-10 21:40

[QUOTE=CRGreathouse;252111]ispower(n) checks if n is a nontrivial power. ispower(n,k) checks if n is a k-th power. ispower(n,,&b) checks if n is a nontrivial power and, if so, stores the base in b. ispower(n,k,&b) checks if n is a k-th power and, if so, stores the base in b.[/QUOTE]

it should be able to check if n is a power of b.

axn 2011-02-10 22:19

[QUOTE=science_man_88;252113]it should be able to check if n is a power of b.[/QUOTE]

use [CODE]b^valuation(n,b)==n[/CODE]

science_man_88 2011-02-10 22:56

[QUOTE=axn;252115]use [CODE]b^valuation(n,b)==n[/CODE][/QUOTE]

or just:

[CODE](18:55)>a=0;for(n=1,#mersenne1,if(valuation(((2^mersenne1[n]-1)%n)+1,2)!=0,a=a+1));return(a)
%51 = 35[/CODE] to check my count, which is accurate, for the amount in mersenne1 . my count wasn't as accurate for the primes and the indexes withing the primes. looks like from primes 1 to 100 we have 75 that fit my criteria not 67.5%

CRGreathouse 2011-02-10 23:57

[QUOTE=science_man_88;252113]it should be able to check if n is a power of b.[/QUOTE]

I've needed that too, at least in the case of powers of two. axn's solution is good. For a faster method, I implemented
[code]long
ispow2(GEN n)
{
if (typ(n) != t_INT)
pari_err(arither1, "ispow2");
if (signe(n) < 1)
return 0;
pari_sp ltop = avma;

GEN xp = int_LSW(n);
long lx = lgefint(n);
ulong u = *xp;
long i = 3;
for (; i < lx; ++i)
{
if (u != 0)
{
avma = ltop;
return 0;
}
xp = int_nextW(xp);
u = *xp;
}
avma = ltop;
return !(u & (u-1));
}[/code]
in C and load it into gp. For small numbers it takes ~11.6 ns vs. 205 ns for axn's version or 175 ns for the shifting variant. For large numbers, especially those divisible by large powers of two, the difference is greater.

Of course, this is comparing apples and oranges. The internal form of axn's algorithm (when appropriately optimized) is
[code]long
ispow2(GEN n)
{
if (typ(n) != t_INT)
pari_err(arither1, "ispow2");
if (signe(n) < 1)
return 0;
pari_sp ltop = avma;

long ret = expi(n) == vali(n);
avma = ltop;
return ret;
}[/code]
which is nearly as fast for small inputs, maybe 25% slower.

CRGreathouse 2011-02-11 00:14

[QUOTE=science_man_88;252118]or just:

[CODE](18:55)>a=0;for(n=1,#mersenne1,if(valuation(((2^mersenne1[n]-1)%n)+1,2)!=0,a=a+1));return(a)
%51 = 35[/CODE] to check my count, which is accurate, for the amount in mersenne1 . my count wasn't as accurate for the primes and the indexes withing the primes. looks like from primes 1 to 100 we have 75 that fit my criteria not 67.5%[/QUOTE]

Your code isn't testing if the numbers are powers of two, just if they're odd. Shorter version:
[code]sum(n=1,#mersenne1,((2^mersenne1[n]-1)%n)%2)[/code]

science_man_88 2011-02-11 00:42

[QUOTE=CRGreathouse;252128]Your code isn't testing if the numbers are powers of two, just if they're odd. Shorter version:
[code]sum(n=1,#mersenne1,((2^mersenne1[n]-1)%n)%2)[/code][/QUOTE]

okay it still gave a correct answer, here's the version you seem to want:

[CODE]a=0;for(n=1,#mersenne1,if(2^valuation((((2^mersenne1[n])-1)%n)+1,2)==(((2^mersenne1[n])-1)%n)+1,a=a+1));return(a)[/CODE]

using this on the primes I get:

[CODE](20:38)>a=0;for(n=1,100,if(2^valuation((((2^prime(n))-1)%n)+1,2)==(((2^prime(n))-1)%n)+1,a=a+1));return(a)
%107 = 50
(20:39)>a=0;for(n=1,1000,if(2^valuation((((2^prime(n))-1)%n)+1,2)==(((2^prime(n))-1)%n)+1,a=a+1));return(a)
%108 = 256
(20:39)>a=0;for(n=1,10000,if(2^valuation((((2^prime(n))-1)%n)+1,2)==(((2^prime(n))-1)%n)+1,a=a+1));return(a)
%109 = 971[/CODE]

each time as you get higher fewer and fewer fit what 35/40 of the Mersenne prime exponents appear to fit.

CRGreathouse 2011-02-11 02:15

Incidentally, since you're working with Mersenne primes and exponents a lot, here's an idea for a script for your file:

[code]MeVec=[2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 521, 607, 1279, 2203, 2281, 3217, 4253, 4423, 9689, 9941, 11213, 19937, 21701, 23209, 44497, 86243, 110503, 132049, 216091, 756839, 859433, 1257787, 1398269, 2976221, 3021377, 6972593, 13466917];
Me(n)={
if(n<=#MeVec,
MeVec[n]
,
error("not known")
)
};
addhelp(Me, "Me(n): Returns the n-th Mersenne prime exponent, Sloane's A000043.");

M(n,modulus=0)={
if(n>#MeVec, error("not known"));
if(modulus,
lift(Mod(2,modulus)^MeVec[n]-1)
,
1<<MeVec[n]-1
)
};
addhelp(M, "M(n, {modulus}): Returns the n-th Mersenne prime, Sloane's A000668. If modulus is given, instead return the n-th Mersenne number mod the modulus.");[/code]

(Of course you can expand the vector to include the numbers whose primality is known but whose order is not.)

This should not only save you keystrokes but time on longer calculations since the modular calculations are efficient compared to creating a large number and dividing it.

science_man_88 2011-02-11 19:21

[QUOTE=CRGreathouse;252145]Incidentally, since you're working with Mersenne primes and exponents a lot, here's an idea for a script for your file:

[code]MeVec=[2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 521, 607, 1279, 2203, 2281, 3217, 4253, 4423, 9689, 9941, 11213, 19937, 21701, 23209, 44497, 86243, 110503, 132049, 216091, 756839, 859433, 1257787, 1398269, 2976221, 3021377, 6972593, 13466917];
Me(n)={
if(n<=#MeVec,
MeVec[n]
,
error("not known")
)
};
addhelp(Me, "Me(n): Returns the n-th Mersenne prime exponent, Sloane's A000043.");

M(n,modulus=0)={
if(n>#MeVec, error("not known"));
if(modulus,
lift(Mod(2,modulus)^MeVec[n]-1)
,
1<<MeVec[n]-1
)
};
addhelp(M, "M(n, {modulus}): Returns the n-th Mersenne prime, Sloane's A000668. If modulus is given, instead return the n-th Mersenne number mod the modulus.");[/code]

(Of course you can expand the vector to include the numbers whose primality is known but whose order is not.)

This should not only save you keystrokes but time on longer calculations since the modular calculations are efficient compared to creating a large number and dividing it.[/QUOTE]

I tried something with this can it be confirmed that all Mersenne primes so far found>31 fit (31 or 43) mod 84 I bet I'm missing something trivial. If not can we limit the prime exponents by this ( did test saves 2 prime exponent checks under 1000000 both under 1000)?

CRGreathouse 2011-02-11 21:59

Can you prove that it behaves the way you expect mod 7? Mod 3? Mod 4?

science_man_88 2011-02-11 22:03

[QUOTE=CRGreathouse;252199]Can you prove that it behaves the way you expect mod 7? Mod 3? Mod 4?[/QUOTE]

[CODE]for(n=1,100,for(x=1,#MeVec,print1(M(x,n)","));print("::"n))[/CODE]

the code I showed to prime95 was :

[CODE](11:01)>for(n=1,#mersenne1,print((2^mersenne1[n]-1)%n))
0
1
1
3
1
1
1
7
1
1
6
7
5
1
7
15
1
1
12
7
10
1
15
7
1
1
1
3
26
1
1
31
7
7
31
31
31
31
1
7[/CODE]

science_man_88 2011-02-12 01:06

[QUOTE=CRGreathouse;252199]Can you prove that it behaves the way you expect mod 7? Mod 3? Mod 4?[/QUOTE]

what the valuation code if so no not really, I can't even find a use of valuation because I have no idea what it does.

it can't act modulo( checked with second argument as 3) so I'm stumped, oh never mind I think I realize which valuation it's talking of after going through :

[url]http://en.wikipedia.org/wiki/Valuation_(mathematics[/url])

science_man_88 2011-02-12 14:20

[QUOTE=CRGreathouse;252199]Can you prove that it behaves the way you expect mod 7? Mod 3? Mod 4?[/QUOTE]

[CODE]0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,::3
3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,::4
3,0,3,1,1,3,1,1,1,3,3,1,3,1,1,1,1,1,3,1,3,3,3,3,3,1,1,3,1,1,1,3,3,1,3,3,3,3,1,::7
3,7,31,43,43,31,43,43,43,31,31,43,31,43,43,43,43,43,31,43,31,31,31,31,31,43,43,31,43,43,43,31,31,43,31,31,31,31,43,::84[/CODE]

84 is a multiple of 3,4 and 7 so it's a matter of proving 31 and 43 align with the modulo 3,4 and 7 need according to this. 31%3 == 1 that aligns 43%3== 1 that also aligns, 31%4==3 that aligns and 43%4==3 so that also aligns, 31%7 =3 that aligns, 43%7 = 1 this also can align with the modulo so it is in fact plausible that this is true. Plausible is not a proof I know that. The proof needs to be more rigorous, if it's true maybe we can eliminate more checks when combined with 24x+7. though we know what x are mersenne numbers(it doesn't look to help though).

science_man_88 2011-02-12 18:18

[CODE]minimumprime(x) = if((x-2)%2==1,print(6*((x-2)-floor((x-2)*.5))-1),print(6*floor((.5*(x-2)))+1))[/CODE]

just thought I'd put this up.

science_man_88 2011-02-12 22:15

[QUOTE=science_man_88;252271][CODE]minimumprime(x) = if((x-2)%2==1,print(6*((x-2)-floor((x-2)*.5))-1),print(6*floor((.5*(x-2)))+1))[/CODE]

just thought I'd put this up.[/QUOTE]

and one thing I just learned by playing alt + 1(numpad) gets you back to the start of text and yes I get it's the control character start of header. something tells me if i every make long scripts that act as one line in PARI I'll be able to use it.

CRGreathouse 2011-02-13 01:34

[QUOTE=science_man_88;252271][CODE]minimumprime(x) = if((x-2)%2==1,print(6*((x-2)-floor((x-2)*.5))-1),print(6*floor((.5*(x-2)))+1))[/CODE]

just thought I'd put this up.[/QUOTE]

Don't print, just return! That way you can actually use the value instead of merely looking at it.

science_man_88 2011-02-13 01:40

[QUOTE=CRGreathouse;252305]Don't print, just return! That way you can actually use the value instead of merely looking at it.[/QUOTE]

I changed it to returning it the hard part now is I was trying to make a script that could use that value to get the real answer because by x = 1000 it's off by about 5000:

[CODE]v=vector(13,n,minimumprime(n));v[1]=2;v[2]=3;c=[];a=0;for(x=1,#v,for(y=1,x-1,if(v[x]%v[y]==0,a=a+1));if(!a>1,c=concat(c,[v[x]])))[/CODE] is my best attempt so far, but anyways I'm kinda repeating the same task I know.

CRGreathouse 2011-02-13 01:52

I don't know what the function is supposed to return, nor by what measure it's off by 5000.

science_man_88 2011-02-13 02:00

[QUOTE=CRGreathouse;252307]I don't know what the function is supposed to return, nor by what measure it's off by 5000.[/QUOTE]

[CODE](21:44)>prime(1000)
%38 = 7919
(22:00)>minimumprime(1000)
%39 = 2995[/CODE]

CRGreathouse 2011-02-13 02:14

[QUOTE=science_man_88;252310][CODE](21:44)>prime(1000)
%38 = 7919
(22:00)>minimumprime(1000)
%39 = 2995[/CODE][/QUOTE]

Are you trying to have it return the n-th prime? Why not use prime() in that case? :confused:

science_man_88 2011-02-13 11:35

[QUOTE=CRGreathouse;252313]Are you trying to have it return the n-th prime? Why not use prime() in that case? :confused:[/QUOTE]

because I want to try and make one that works without that so it's not limited by primelimit.

CRGreathouse 2011-02-13 21:49

[QUOTE=science_man_88;252333]because I want to try and make one that works without that so it's not limited by primelimit.[/QUOTE]

There are very few cases when you'd actually want to do that. In almost all cases the better approach would be to calculate more primes and keep them on hand -- and you can already do this with
default(primelimit, 10^9)
changing the number as appropriate.

But if you insist, here's a really, really, really bad way:
[code]primeN(n)={
my(p=1);
while(n,
if(isprime(p++), n--)
);
p
};[/code]

But please don't do this! It takes 10.9 seconds to do primeN(10^6), compared to 1.5 milliseconds for prime(10^6) -- seven thousand times slower. primeN(10^7) takes several minutes while prime(10^7) gives an answer in ~10 ms.

science_man_88 2011-02-13 21:59

[QUOTE=CRGreathouse;252391]There are very few cases when you'd actually want to do that. In almost all cases the better approach would be to calculate more primes and keep them on hand -- and you can already do this with
default(primelimit, 10^9)
changing the number as appropriate.

But if you insist, here's a really, really, really bad way:
[code]primeN(n)={
my(p=1);
while(n,
if(isprime(p++), n--)
);
p
};[/code]

But please don't do this! It takes 10.9 seconds to do primeN(10^6), compared to 1.5 milliseconds for prime(10^6) -- seven thousand times slower. primeN(10^7) takes several minutes while prime(10^7) gives an answer in ~10 ms.[/QUOTE]

fine what ever. want what my minimumprime (x) does ? or have you figured it out ? all it does is x-2 to find out how many to go into the [TEX]6n\pm1[/TEX] figures which one it's on and does a calucation based on where it lands to figure out a n value than it returns [TEX]6n\pm1[/TEX]

CRGreathouse 2011-02-14 00:51

Yes, and it's clear that that won't give the right answer because it's linear and the correct figure is about n log n.

This may be of interest to you:
[url]http://primes.utm.edu/nthprime/[/url]

science_man_88 2011-02-14 01:23

[QUOTE=CRGreathouse;252414]Yes, and it's clear that that won't give the right answer because it's linear and the correct figure is about n log n.

This may be of interest to you:
[url]http://primes.utm.edu/nthprime/[/url][/QUOTE]

Since it's built to give the lowest possible value it could be I know it's not correct. I'm just trying to figure a easy way to correct for it. The first time it's wrong is x = 10 last I checked.

CRGreathouse 2011-02-14 02:55

[QUOTE=science_man_88;252416]Since it's built to give the lowest possible value it could be I know it's not correct. I'm just trying to figure a easy way to correct for it. The first time it's wrong is x = 10 last I checked.[/QUOTE]

I gave pretty much the easiest way to do it without using the built-in prime list. There are decently fast methods that use the prime list but avoid the prime() command itself, if you're interested. For example,
[code]primeN(n)={
my(v=primes(n));
v[#v]
};[/code]
and
[code]primeN(n)={
if(n<6, return([2,3,5,7,11][n]));
my(l=log(n), lower=ceil(n*l+n*(log(l)-1)), upper=floor(n*l+n*log(l)), t);
while(upper - lower > 1,
if(primepi(t = (lower + upper) \ 2) < n, lower = t, upper = t)
);
if(primepi(lower) == n, lower, upper)
};[/code]

These might be a hundred times slower rather than 10,000 times slower.

CRGreathouse 2011-02-14 03:01

By the way, I think that re-implementing basic features of a language is a good way to learn. Working through Prime Numbers: A Computational Perspective really helped me to understand how the various algorithms fit together. So I don't think there's anything wrong with doing this -- I just want you to code it, learn, and then use the built-in functions rather than your own. :smile:

science_man_88 2011-02-16 01:37

[QUOTE=CRGreathouse;252426]By the way, I think that re-implementing basic features of a language is a good way to learn. Working through Prime Numbers: A Computational Perspective really helped me to understand how the various algorithms fit together. So I don't think there's anything wrong with doing this -- I just want you to code it, learn, and then use the built-in functions rather than your own. :smile:[/QUOTE]

yeah bad part is I'd be lucky to program all the things my pocket calculator can do. I have an idea for one thing I haven't done before mainly because it's like the base converter.

CRGreathouse 2011-02-16 12:39

[QUOTE=science_man_88;252634]yeah bad part is I'd be lucky to program all the things my pocket calculator can do. I have an idea for one thing I haven't done before mainly because it's like the base converter.[/QUOTE]

Actually that's a good project because PARI doesn't know how to do that on its own. I mean, it inputs and displays numbers in decimal while storing and calculating in binary, and it can return a number in binary (as a vector) if requested, but it doesn't have generic base-conversion routines.

science_man_88 2011-02-16 13:59

[QUOTE=CRGreathouse;252671]Actually that's a good project because PARI doesn't know how to do that on its own. I mean, it inputs and displays numbers in decimal while storing and calculating in binary, and it can return a number in binary (as a vector) if requested, but it doesn't have generic base-conversion routines.[/QUOTE]

we made the base converter last I checked, a list of what my calculator can do:

[QUOTE]list of scientific constant (not sure what all of them are)
unit conversion(40 in total, but I think I could cut that to 20)
powers (not so hard in PARI)
base 10 and natural logarithms (not so hard)
some SI prefixes ( mega,giga,tera,micro,milli,kilo,fempto,pico,nano)
combinations and permutations
trig functions
xth roots ( pretty much done)
convert answer to fraction if possible.
insert function.
9 variable storage.
Pol( and Rec( (haven't used them but I'm pretty sure I know what they do).
integrals ?
"random" numbers
rounding ( not hard)
factorial
Pi
a+bi ( seen the equation before don't know what it is)
percentage
[TEX]Re\leftrightarrow lm[/TEX]
degrees rad and grad
base conversion( DEC,HEX,BIN,OCT)[/QUOTE]

modes are : [QUOTE]COMP,CMPLX,SD,REG,BASE,EQN,MAT,VCT,Deg,Rad,Gra,Fix,Sci,Norm,Disp[/QUOTE]

CRGreathouse 2011-02-16 20:16

Most of those are of course already in GP. If you want to implement some of the others, give it a shot and I'll take a look.

science_man_88 2011-02-16 21:22

[QUOTE=CRGreathouse;252703]Most of those are of course already in GP. If you want to implement some of the others, give it a shot and I'll take a look.[/QUOTE]

the conversion among units is basically a base conversion with each digit in the number as a different conversion base.

I'm not completely sure which others aren't.

science_man_88 2011-02-16 22:13

[CODE](18:10)>C= vector(100,n,(2*n)!/((n+1)!*n!))
%10 = [1, 2, 5, 14, 42, 132, 429, 1430, 4862,
(18:10)>Catalan(n) = return(C[n+1])
%11 = (n)->return(C[n+1])
(18:12)>Catalan(0)
%12 = 1
(18:12)>Catalan(2)
%13 = 5[/CODE]

does this cover

[url]http://rosettacode.org/wiki/Catalan_numbers[/url] ?

science_man_88 2011-02-16 22:46

I forgot [url]http://www.mersenneforum.org/showthread.php?p=248234#post248234[/url] lead to:

[url]http://rosettacode.org/wiki/Luhn_test_of_credit_card_numbers[/url]

done it should be put up.

science_man_88 2011-02-17 01:45

sorry for getting off topic again I see a few others that shouldn't be too difficult in the list but anyways.

CRGreathouse 2011-02-17 04:59

[QUOTE=science_man_88;252712][CODE](18:10)>C= vector(100,n,(2*n)!/((n+1)!*n!))
%10 = [1, 2, 5, 14, 42, 132, 429, 1430, 4862,
(18:10)>Catalan(n) = return(C[n+1])
%11 = (n)->return(C[n+1])
(18:12)>Catalan(0)
%12 = 1
(18:12)>Catalan(2)
%13 = 5[/CODE]

does this cover

[url]http://rosettacode.org/wiki/Catalan_numbers[/url] ?[/QUOTE]

Sure, that'll do -- as long as you don't want Catalan(100) or higher. (Edit: I see that the link suggests memoization, in which case you did it right. As it happens it's probably not needed here.) I would just do
[code]Catalan(n)=(2*n)!/(n+1)!/n![/code]
or, probably faster,
[code]Catalan(n)=binomial(2*n,n+1)/n[/code]

CRGreathouse 2011-02-17 07:42

[QUOTE=science_man_88;252720]I forgot [url]http://www.mersenneforum.org/showthread.php?p=248234#post248234[/url] lead to:

[url]http://rosettacode.org/wiki/Luhn_test_of_credit_card_numbers[/url]

done it should be put up.[/QUOTE]

Great, put it up!

[QUOTE=science_man_88;252736]sorry for getting off topic again I see a few others that shouldn't be too difficult in the list but anyways.[/QUOTE]

You can do them if you like, or post it here if you want feedback.

science_man_88 2011-02-17 12:13

[QUOTE=CRGreathouse;252751]Great, put it up!



You can do them if you like, or post it here if you want feedback.[/QUOTE]

well I kinda see the math of the sierpinski triangle and carpet

going up a step from the triangle shown means multiplying the width by 3 and printing 8/9 of the previous carpet and leaving the center clear. since they show examples with them hitting the smallest possible it's simply working backwards I think.

science_man_88 2011-02-18 00:57

[CODE](20:28)>sdsub(x,y) = C=y;for(a=1,#C,C[a]=x[C[a]+1]);C=vecsort(C);for(b=1,#y,x[y[b]+1]=C[b]);x;
(20:48)>sdsub([7, 6, 5, 4, 3, 2, 1, 0],[6, 1, 7])
%94 = [7, 1, 5, 4, 3, 2, 0, 6][/CODE] this is trying :
[url]http://rosettacode.org/wiki/Sort_disjoint_sublist[/url]

it seems to reverse 2 thing as it places them. oh never mind I see why trying out the simpler solution they suggest, it's because the call to a unsorted y.

CRGreathouse 2011-02-18 01:19

I think this is what PARI's indirect sorting is for. Take the elements you want to sort into a new vector, call it u, and sort with vecsort(u,,1). Rather than sorting u, this will return a vector that tells you where to put the elements. Then you can just put each element where the resulting vector tells you to put it.

science_man_88 2011-02-18 01:35

[QUOTE=CRGreathouse;252861]I think this is what PARI's indirect sorting is for. Take the elements you want to sort into a new vector, call it u, and sort with vecsort(u,,1). Rather than sorting u, this will return a vector that tells you where to put the elements. Then you can just put each element where the resulting vector tells you to put it.[/QUOTE]

[CODE](21:28)>vecsort([1,2,3],,1)
%134 = Vecsmall([1, 2, 3])
(21:30)>vecsort([1,2,3],,2)
%135 = [1, 2, 3]
(21:30)>vecsort([1,2,3],,3)
%136 = Vecsmall([1, 2, 3])
(21:33)>vecsort([1,2,3],,4)
%137 = [3, 2, 1]
(21:33)>vecsort([1,2,3],,5)
%138 = Vecsmall([3, 2, 1])
(21:33)>vecsort([1,2,3],,6)
%139 = [3, 2, 1]
(21:33)>vecsort([1,2,3],,7)
%140 = Vecsmall([3, 2, 1])
(21:33)>vecsort([1,2,3],,8)
%141 = [1, 2, 3]
(21:33)>vecsort([1,2,3],,9)
%142 = Vecsmall([1, 2, 3])
(21:34)>vecsort([1,2,3],,10)
%143 = [1, 2, 3][/CODE]

science_man_88 2011-02-18 01:44

[CODE](21:25)>sdsub1(x,y) = D=y;C=y;for(a=1,#C,C[a]=x[C[a]]);C=vecsort(C);D=vecsort(D);for(b=1,#y,x[D[b]]=C[b]);x;
(21:28)>sdsub1([7, 6, 5, 4, 3, 2, 1, 0],[7,2,8])[/CODE]

okay this works for the 1 based what if I want it to act like it's 0 based.

CRGreathouse 2011-02-18 02:16

[QUOTE=science_man_88;252865]%134 = Vecsmall([1, 2, 3])
(21:30)>vecsort([1,2,3],,2)
%135 = [1, 2, 3]
(21:30)>vecsort([1,2,3],,3)
%136 = Vecsmall([1, 2, 3])
(21:33)>vecsort([1,2,3],,4)
%137 = [3, 2, 1]
(21:33)>vecsort([1,2,3],,5)
%138 = Vecsmall([3, 2, 1])
(21:33)>vecsort([1,2,3],,6)
%139 = [3, 2, 1]
(21:33)>vecsort([1,2,3],,7)
%140 = Vecsmall([3, 2, 1])
(21:33)>vecsort([1,2,3],,8)
%141 = [1, 2, 3]
(21:33)>vecsort([1,2,3],,9)
%142 = Vecsmall([1, 2, 3])
(21:34)>vecsort([1,2,3],,10)
%143 = [1, 2, 3][/CODE][/QUOTE]

Not ,,2. Not ,,3. The flag you want to set is 1 and none of the others. (You can see what they do by reading the manual, or get an overview with ?vecsort, but they're not applicable here.)

science_man_88 2011-02-18 02:41

[QUOTE=CRGreathouse;252871]Not ,,2. Not ,,3. The flag you want to set is 1 and none of the others. (You can see what they do by reading the manual, or get an overview with ?vecsort, but they're not applicable here.)[/QUOTE]

but the value 1 only returns a permutation of order in the vector by the looks of it. so unless you go through the indexes I don't see how.

CRGreathouse 2011-02-18 17:33

[QUOTE=science_man_88;252877]go through the indexes[/QUOTE]

Look at that, he has an idea!

science_man_88 2011-02-18 20:10

[QUOTE=CRGreathouse;252931]Look at that, he has an idea![/QUOTE]

It's not having an idea that's new, It's that it may work for once I'm guessing.

CRGreathouse 2011-02-18 20:52

[QUOTE=science_man_88;252964]It's not having an idea that's new, It's that it may work for once I'm guessing.[/QUOTE]

So you have a vector you need to put in order, plus two vectors that tell you how to order them. The first is a list of the indexes that should be reordered. The second is the ordering for those indexes.

science_man_88 2011-02-18 21:45

[QUOTE=CRGreathouse;252967]So you have a vector you need to put in order, plus two vectors that tell you how to order them. The first is a list of the indexes that should be reordered. The second is the ordering for those indexes.[/QUOTE]

my code (#2210) does that but adds one thing instead of ordering y ( because it may be used again afterwards), it orders D to get the order of the indexes.

science_man_88 2011-02-18 21:56

[QUOTE=CRGreathouse;252751]Great, put it up!



You can do them if you like, or post it here if you want feedback.[/QUOTE]

I think I should put up a few more looks like I forgot a few I tried out. anyone else want to put them up ? I'm pretty sure by the looks of it I forgot to put up the code for the base conversion code.

science_man_88 2011-02-20 23:31

[QUOTE=CRGreathouse;252145]Incidentally, since you're working with Mersenne primes and exponents a lot, here's an idea for a script for your file:

[code]MeVec=[2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 521, 607, 1279, 2203, 2281, 3217, 4253, 4423, 9689, 9941, 11213, 19937, 21701, 23209, 44497, 86243, 110503, 132049, 216091, 756839, 859433, 1257787, 1398269, 2976221, 3021377, 6972593, 13466917];
Me(n)={
if(n<=#MeVec,
MeVec[n]
,
error("not known")
)
};
addhelp(Me, "Me(n): Returns the n-th Mersenne prime exponent, Sloane's A000043.");

M(n,modulus=0)={
if(n>#MeVec, error("not known"));
if(modulus,
lift(Mod(2,modulus)^MeVec[n]-1)
,
1<<MeVec[n]-1
)
};
addhelp(M, "M(n, {modulus}): Returns the n-th Mersenne prime, Sloane's A000668. If modulus is given, instead return the n-th Mersenne number mod the modulus.");[/code]

(Of course you can expand the vector to include the numbers whose primality is known but whose order is not.)

This should not only save you keystrokes but time on longer calculations since the modular calculations are efficient compared to creating a large number and dividing it.[/QUOTE]

this gives errors I should of figured. M(29,233) should be 0 it's 36 apparently. I found the error in the code and fixed it on my end.

science_man_88 2011-02-21 00:37

[QUOTE=science_man_88;253203]this gives errors I should of figured. M(29,233) should be 0 it's 36 apparently. I found the error in the code and fixed it on my end.[/QUOTE]

oh wait another error popped out of it. It won't let n be greater than #MeVec so when you say n-th mersenne number you mean with n<=#MeVec. Also even with lower ones if you use your form we get n entangled with MeVec[n] in one of the equations in the second if, it doesn't work for me in that form but works fine with it replaced with n. so with n < #MeVec it could be confusing n between the index into the mersenne prime exponents or the mersenne number index, so that part doesn't work which pretty much fails the function.

the most proper to the definition you gave that I can get is:

[CODE]M(n,modulus=0)=if(modulus,return(lift(Mod(2,modulus)^n-1)),if(n>#MeVec,error("not known"),return(2^MeVec[n]-1)))[/CODE]

and I never had to learn the math behind it, I knew one fact that told me the old one was wrong.

CRGreathouse 2011-02-21 01:11

[QUOTE=science_man_88;253203]this gives errors I should of figured. M(29,233) should be 0 it's 36 apparently. I found the error in the code and fixed it on my end.[/QUOTE]

Huh? The 29-th Mersenne prime is a prime (greater than 233) so it isn't divisible by 233.

[QUOTE=science_man_88;253208]oh wait another error popped out of it. It won't let n be greater than #MeVec so when you say n-th mersenne number you mean with n<=#MeVec.[/QUOTE]

Well sure -- I can't tell you what the millionth Mersenne number is mod 101 since I don't know what that will be.

[QUOTE=science_man_88;253208]Also even with lower ones if you use your form we get n entangled with MeVec[n] in one of the equations in the second if, it doesn't work for me in that form but works fine with it replaced with n. so with n < #MeVec it could be confusing n between the index into the mersenne prime exponents or the mersenne number index, so that part doesn't work which pretty much fails the function.[/QUOTE]

This is what the help entries are for.

science_man_88 2011-02-21 01:17

[QUOTE]Huh? The 29-th Mersenne prime is a prime (greater than 233) so it isn't divisible by 233.[/QUOTE]

ah but your definition states if modulus is present it goes to mersenne numbers not mersenne primes so since the 233 is there it checks M29 [B]not[/B] 2^MeVec[29]-1 by definition.


[QUOTE]Well sure -- I can't tell you what the millionth Mersenne number is mod 101 since I don't know what that will be.



This is what the help entries are for.[/QUOTE]

yes I know and they helped me spot a flaw in the function.

CRGreathouse 2011-02-21 01:33

I see what you're saying. Yes, the documentation is in error; M(n,m) should be equal to M(n)%m.

science_man_88 2011-02-21 01:35

[QUOTE=CRGreathouse;253213]I see what you're saying. Yes, the documentation is in error; M(n,m) should be equal to M(n)%m.[/QUOTE]

I got it working to the definition completely in post 2220 is my code, the reason I knew to pursue this argument is because I had tested it and I knew that M(x,x) should return 1 and it didn't so I knew at least something was up. Then I found the rest of the errors I listed, nice battle CRG.

CRGreathouse 2011-02-21 01:41

I'm glad you got it to do what you wanted. I specifically did not want that behavior because it makes the function do two unrelated things depending on the argument given, where my version does two very much related things. (You can see an omitted argument as using an 'infinite' m, if you like...)

science_man_88 2011-02-21 01:44

[QUOTE=CRGreathouse;253215]I'm glad you got it to do what you wanted. I specifically did not want that behavior because it makes the function do two unrelated things depending on the argument given, where my version does two very much related things. (You can see an omitted argument as using an 'infinite' m, if you like...)[/QUOTE]

okay it's a multi-tool, thanks for the game of find the error. Does the missing argument use up a lot of extra memory ?

CRGreathouse 2011-02-21 01:54

[QUOTE=science_man_88;253216]okay it's a multi-tool, thanks for the game of find the error. Does the missing argument use up a lot of extra memory ?[/QUOTE]

Any optional argument that you don't give to gp is passed to PARI as null. This value takes up wordsize -- 4 bytes on 32-bit systems and 8 bytes on 64-bit systems. So no, it takes almost exactly no memory.

science_man_88 2011-03-01 01:45

could :

[CODE]createvar(x,y) = v=v;v=concat(v,[x]);v=concat(v,[y]);[/CODE]

and

[CODE]questionvar(x) = d=0;forstep(c=1,#v-1,[2],if(v[c]==x,return(v[c+1]),d=d+1;if(d==#v/2,error("no variable found"))))[/CODE]

qualify for the create a variable on a fly script needed in the list ? I know I should also make sure they can alter them, I could make another script altervar(x) if need be. edit never mind I don't see one like I thought but I do see a use in one sense.

CRGreathouse 2011-03-01 08:44

[QUOTE=science_man_88;254030]could :

[CODE]createvar(x,y) = v=v;v=concat(v,[x]);v=concat(v,[y]);[/CODE]

and

[CODE]questionvar(x) = d=0;forstep(c=1,#v-1,[2],if(v[c]==x,return(v[c+1]),d=d+1;if(d==#v/2,error("no variable found"))))[/CODE]

qualify for the create a variable on a fly script needed in the list ? I know I should also make sure they can alter them, I could make another script altervar(x) if need be. edit never mind I don't see one like I thought but I do see a use in one sense.[/QUOTE]

I don't know what these are intended to do. The first is equivalent to
[CODE]createvar(x,y) = v=concat(v,[x,y])[/CODE]
and the second (apart from the gratuitous variable-clobbering) to
[CODE]questionvar(x) = forstep(c=1,#v-1,2,if(v[c]==x,return(v[c+1])));error("no variable found")[/CODE]

Both rely on a global variable v to do... something. Actually it looks like v is a lookup table of sorts. (In other languages you'd use a hash table to do this.)

science_man_88 2011-03-01 14:38

[QUOTE=CRGreathouse;254063]I don't know what these are intended to do. The first is equivalent to
[CODE]createvar(x,y) = v=concat(v,[x,y])[/CODE]
and the second (apart from the gratuitous variable-clobbering) to
[CODE]questionvar(x) = forstep(c=1,#v-1,2,if(v[c]==x,return(v[c+1])));error("no variable found")[/CODE]

Both rely on a global variable v to do... something. Actually it looks like v is a lookup table of sorts. (In other languages you'd use a hash table to do this.)[/QUOTE]

basically the name and value go into v together the questionvar(x) checks to see if variable x is in v ( the forstep comes from skipping values the variables should be equal to). If the variable name is found it returns the value, else it returns an error telling them they are looking up a non existent variable within the variables in v). I was thinking of making one to change a value in v (so they can be changed as they are wanted to be variable).

science_man_88 2011-04-02 12:16

I made a lucas sequence code if I haven't already ( I couldn't find it if I did).

[CODE]lucas(T,p,q,x)= if(T==1,v=[0,1],v=[2,p]);for(y=3,x,v=concat(v,p*v[y-1]-q*v[y-2]));v=vector(x,i,v[i]);
addhelp(lucas,"lucas(T,p,q,x): returns lucas sequence of type T, 1 is U, 2 is V, for p and q to x length")[/CODE]

of course this could probably be done better.

science_man_88 2011-04-28 01:01

[CODE]pieces=[N,B,K,R,Q,a,b,c,d,e,f,g,h];files=[a,b,c,d,e,f,g,h];rows=[1,2,3,4,5,6,7,8];captures=[x,""];for(piece=1,#pieces,for(capture=1,#captures,for(file=1,#files,for(row=1,#rows,print1(pieces[piece]captures[capture]files[file]rows[row]",")))))[/CODE]

once again I had it working ( i thought) but can't find my error this code is for:

[url]http://www.mersenneforum.org/showthread.php?p=259371[/url]

accidently defined c as a vector that's part of it.


All times are UTC. The time now is 23:20.

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