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 2014-05-23 16:46

Would you paste in your startup message, like
[code] GP/PARI CALCULATOR Version 2.8.0 (development 16370-23a2ae3)
i686 running mingw (ix86/GMP-5.1.3 kernel) 32-bit version
compiled: Apr 14 2014, gcc version 4.8.2 (GCC)
threading engine: single
(readline v6.2 enabled, extended help enabled)

Copyright (C) 2000-2014 The PARI Group

PARI/GP is free software, covered by the GNU General Public License, and comes
WITHOUT ANY WARRANTY WHATSOEVER.[/code]
? Also try this not-actually-parallel version:
[code]parmatrix(m,n,expr=(a,b)->0)={
my(a=matrix(m,n));
for(x=1,m,
for(y=1,n,
a[x,y]=expr(x,y)
)
);
a
};[/code]

science_man_88 2014-05-23 17:07

[QUOTE=CRGreathouse;374095]Would you paste in your startup message, like
[code] GP/PARI CALCULATOR Version 2.8.0 (development 16370-23a2ae3)
i686 running mingw (ix86/GMP-5.1.3 kernel) 32-bit version
compiled: Apr 14 2014, gcc version 4.8.2 (GCC)
threading engine: single
(readline v6.2 enabled, extended help enabled)

Copyright (C) 2000-2014 The PARI Group

PARI/GP is free software, covered by the GNU General Public License, and comes
WITHOUT ANY WARRANTY WHATSOEVER.[/code]
? Also try this not-actually-parallel version:
[code]parmatrix(m,n,expr=(a,b)->0)={
my(a=matrix(m,n));
for(x=1,m,
for(y=1,n,
a[x,y]=expr(x,y)
)
);
a
};[/code][/QUOTE]

[QUOTE]Reading GPRC: gprc.txt ...Done.

GP/PARI CALCULATOR Version 2.7.0 (released)
i686 running mingw (ix86/GMP-5.1.3 kernel) 32-bit version
compiled: Mar 23 2014, gcc version 4.6.3 (GCC)
threading engine: single
(readline v6.2 enabled, extended help enabled)

Copyright (C) 2000-2014 The PARI Group

PARI/GP is free software, covered by the GNU General Public License, and comes WITHOUT ANY WARRANTY WHATSOEVER.[/QUOTE]

[CODE](12:44) gp > parmatrix(m,n,expr=(a,b)->0)={
my(a=matrix(m,n));
for(x=1,m,
for(y=1,n,
a[x,y]=expr(x,y)
)
);
a
};
(14:03) gp > parmatrix(10,10,f(x,y)=x^2+sigma(y))
%7 =
[ 2 4 5 8 7 13 9 16 14 19]

[ 5 7 8 11 10 16 12 19 17 22]

[ 10 12 13 16 15 21 17 24 22 27]

[ 17 19 20 23 22 28 24 31 29 34]

[ 26 28 29 32 31 37 33 40 38 43]

[ 37 39 40 43 42 48 44 51 49 54]

[ 50 52 53 56 55 61 57 64 62 67]

[ 65 67 68 71 70 76 72 79 77 82]

[ 82 84 85 88 87 93 89 96 94 99]

[101 103 104 107 106 112 108 115 113 118]

(14:04) gp > parmatrix(m,n,expr=(a,b)->0)=my(a=matrix(m,n));parfor(x=1,m,for(y=1,n,a[x,y]=expr(x,y)));a;
(14:04) gp > parmatrix(10,10,f(x,y)=x^2+sigma(y))
%9 =
[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 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 0 0]

[0 0 0 0 0 0 0 0 0 0]

(14:04) gp > parmatrix(m,n,expr=(a,b)->0)=my(a=matrix(m,n));for(x=1,m,parvector(n,y,a[x,y]=expr(x,y)));a;
(14:05) gp > parmatrix(10,10,f(x,y)=x^2+sigma(y))
%11 =
[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 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 0 0]

[0 0 0 0 0 0 0 0 0 0]
[/CODE]

it appears the parallel loops are the problem.

CRGreathouse 2014-05-23 17:08

Actually, on second thought, this might be a problem with the threading model -- the other threads can't change the global matrix a. Hmm. I don't know of an easy solution, maybe there is something that can be done with inline().

The 'right' solution would require working in PARI. You'd generate each column in parallel, then combine them with something like mkmat.

science_man_88 2014-05-23 17:50

[QUOTE=CRGreathouse;374097]Actually, on second thought, this might be a problem with the threading model -- the other threads can't change the global matrix a. Hmm. I don't know of an easy solution, maybe there is something that can be done with inline().

The 'right' solution would require working in PARI. You'd generate each column in parallel, then combine them with something like mkmat.[/QUOTE]

[CODE]parmatrix(m,n,expr=(a,b)->0)={
my(a=matrix(m,n));
for(x=1,m,
a[x,]=parvector(n,y,eval(expr(x,y)))
);a;[/CODE]

one step closer it appears to fall flat with both parallel though. I even tried to turn it into a vector and do two parvector(); edit: if you want to alter the reading to column's I could make a flag that if 1 could go by column and if not go by row.

science_man_88 2014-05-24 12:00

I tried the vector one again but with only one parallel function ( though I wonder if breaking parvector into pareval and vector could work):

[CODE]parmatrix1(m,n,expr(a,b)->0)
=my(a=Vec(matrix(m,n)));
for(x=1,n,
a[x]=parvector(m,y,eval(expr(x,y)))
);
a=Mat(a)[/CODE]

I'm thinking of making a check to see which is greater and working it that way the problem is this rotates the values of m and n if n is greater than m because doing the greater of the two in parvector got me down to around 8200 ms for one of the things I tried in expr ( admittedly in a 1000 by 1100 matrix) but I'd also have to change the order in which m and n are declared in a to not mess things up and this rotates the matrix.

science_man_88 2014-05-24 21:35

I've been reading up on set's etc. ( partly brought on my the fog comment by RDS) and I'm wondering why forpart() doesn't allow the variable to go over the partitions of a set (aka the elements of the power set of the set x=[1,2,3,4,5,6,7,8,9,10,11,12] for example). I'm guessing there's probably already a Set based function like this ?

CRGreathouse 2014-05-27 14:19

I notice also that you have "threading model: single" which means that you can't execute code in parallel anyway.

CRGreathouse 2014-05-27 14:23

[QUOTE=science_man_88;374200]I've been reading up on set's etc. ( partly brought on my the fog comment by RDS) and I'm wondering why forpart() doesn't allow the variable to go over the partitions of a set (aka the elements of the power set of the set x=[1,2,3,4,5,6,7,8,9,10,11,12] for example). I'm guessing there's probably already a Set based function like this ?[/QUOTE]

If you want to loop over partitions, use
[code]for(i=0,n!-1, p = numtopart(n, i); print(p))[/code]
where you replace prime(p) with useful code and n with the number of items in the partition.

If you want to loop over subsets of a set, use
[code]for(i=0,2^#v-1, p = vecextract(v,i); print(p))[/code]
where v is the set and print(p) is replaced by useful code.

science_man_88 2014-05-28 01:38

[QUOTE=CRGreathouse;374375]
If you want to loop over subsets of a set, use
[code]for(i=0,2^#v-1, p = vecextract(v,i); print(p))[/code]
where v is the set and print(p) is replaced by useful code.[/QUOTE]

once I saw them in order I thought it would cut a check out of aligentry0 but for n=73 it increases the number of vectors to checks to >19 * numbpart(73)^3 when the whole before took numbpart(73)+14.

science_man_88 2014-10-03 22:36

maybe make a version that doesn't crash gp.exe
 
I made a recent edit on the OEIS but I know a further formula from it that leads to the Fermat quotients however my codes ( based on quick ways of setting things up that I've read in the OEIS) crashes at b(1152), I was hoping someone could suggest an edit that wouldn't crash it:

[CODE]a(n)=if(n,4*a(n-1)+1,0)[/CODE]// calculates A002450
[CODE]b(n)=p=prime(n);(3*a((p-1)/2))/p[/CODE] // calculates the Fermat quotients assuming p is calculable ( yes I know that form of finding it has a limit)

CRGreathouse 2014-10-04 04:38

[QUOTE=science_man_88;384315]I was hoping someone could suggest an edit that wouldn't crash it:

[CODE]a(n)=if(n,4*a(n-1)+1,0)[/CODE]// calculates A002450
[CODE]b(n)=p=prime(n);(3*a((p-1)/2))/p[/CODE] // calculates the Fermat quotients assuming p is calculable ( yes I know that form of finding it has a limit)[/QUOTE]

You want to know why a(4655) fails. Does (4^4655 - 1)/3 fail also?


All times are UTC. The time now is 21:15.

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