PARI's commands

[CRGreathouse]

Would you paste in your startup message, like

          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.

? Also try this not-actually-parallel version:

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
};

[science_man_88]

Originally Posted by CRGreathouse

Would you paste in your startup message, like

          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.

? Also try this not-actually-parallel version:

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
};

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.

(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]

it appears the parallel loops are the problem.

[CRGreathouse]

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]

Originally Posted by CRGreathouse

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.

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;

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]

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

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)

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]

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]

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

[CRGreathouse]

Originally Posted by science_man_88

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 ?

If you want to loop over partitions, use

for(i=0,n!-1, p = numtopart(n, i); print(p))

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

for(i=0,2^#v-1, p = vecextract(v,i); print(p))

where v is the set and print(p) is replaced by useful code.

[science_man_88]

Originally Posted by CRGreathouse

If you want to loop over subsets of a set, use

for(i=0,2^#v-1, p = vecextract(v,i); print(p))

where v is the set and print(p) is replaced by useful code.

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]

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:

a(n)=if(n,4*a(n-1)+1,0)

// calculates A002450

b(n)=p=prime(n);(3*a((p-1)/2))/p

// calculates the Fermat quotients assuming p is calculable ( yes I know that form of finding it has a limit)

[CRGreathouse]

Originally Posted by science_man_88

I was hoping someone could suggest an edit that wouldn't crash it:

a(n)=if(n,4*a(n-1)+1,0)

// calculates A002450

b(n)=p=prime(n);(3*a((p-1)/2))/p

// calculates the Fermat quotients assuming p is calculable ( yes I know that form of finding it has a limit)

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