I updated it with many new functions. I excluded _aprcl_step4_worker and _aprcl_step6_worker which, despite having GP help entries, are not callable functions in gp. I also added a section where obsolete/deprecated/removed functions can be added. I removed them from the main list -- you shouldn't be quizzed on functions which are gone or will soon be gone -- and also prevented them from being randomly generated, since a user could legitimately recognize them and how to use them, so they're not good camouflage for true functions.

[url]http://math.crg4.com/gp-function-quiz.html[/url]

at least I can now do my powerset pastebin ( maybe a little once i bring back the concat calls) easier:

[CODE]powerset(S)=my(a=S);forsubset(#a,s,y=apply(r->a[r],Vec(s));print(y))[/CODE]

Since I've forgotten a few times, this works on MacOS:

[CODE]brew install readline
brew info readline
# Take include and lib from the output of the above and use here:
./Configure --with-readline --with-readline-include=/usr/local/opt/readline/include --with-readline-lib=/usr/local/opt/readline/lib
make -j8 gp[/CODE]

The Configure script should say something like:

[FONT="Courier New"]...Found readline header in /usr/local/opt/readline/include/readline
Using GNU readline, version 7.0[/FONT]

and not:

[FONT="Courier New"]###
### Editline wrapper detected, building without readline support
###
### Building without GNU readline support[/FONT]

A couple cool things in recent Pari/GP:

ECPP. See Jared's presentation here: [url]http://guissmo.com/files/lfantseminar/ecpplong.pdf[/url]

It makes Primo v4 certificates if desired, which is pretty cool. It seems to run at a good speed. I haven't spent much time with it (if only I didn't have to work every day...).


forfactored and forsquarefree. Pretty neat functions that run through a range and return the number and its factorization. I added this to Perl/ntheory, with both positive and negative results.

1) A bunch of Pari/GP functions take either n or the factorization of n as input. You can see a bunch of examples of this Charles added to OEIS sequences. I have ranged euler and mobius functions, but this is a more clever way to do it. It also works on some other functions like isfundamental. My only saving grace is faster factoring which helps some (and could be added to Pari/GP).

2) The current Pari/GP implementation is limited to 2^64, and is also really slow near that limit. It's chunk size is 1024 and it's sieving to sqrt(b) for each chunk's a,b.

[code]
? s=0; forfactored(n=2,10^7,s+=n[2][,1][1]); s
time = 3,500 ms.
%5 = 3203714961609
? s=0; forfactored(n=10^12,10^12+10^5,s+=n[2][,1][1]); s
time = 345 ms.
%6 = 3614000705858566
? s=0; forfactored(n=10^16,10^16+10^5,s+=n[2][,1][1]); s
time = 34,346 ms.
%7 = 26970000029808967678
? s=0; forfactored(n=10^17,10^17+10^5,s+=n[2][,1][1]); s
time = 1min, 47,390 ms.
%8 = 251500000082904735996

? s=0; for(n=10^17,10^17+10^5,d=factorint(n);s+=d[,1][1]); s
time = 10,975 ms.
%9 = 251500000082904735996
[/code]

Making 10^5 longer just makes it worse. ntheory looks to be about 10x faster for factoring at this size (I'm surprised it is that much), and the forfactored is a lot faster.

[code]
# n is in $_ and factors in @_
time perl -Mntheory=:all -E '$s=0; forfactored { $s+=$_[0]; } 10**7; say $s'
3203714961609
user 0m2.525s
time perl -Mntheory=:all -E '$s=0; forfactored { $s+=$_[0]; } 10**12,10**12+10**5; say $s'
3614000705858566
user 0m0.078s
time perl -Mntheory=:all -E 'use Math::GMPz; $s=Math::GMPz->new(0); forfactored { $s+=$_[0]; } 10**16,10**16+10**5; say $s'
26970000029808967678
user 0m0.749s
time perl -Mntheory=:all -E 'use Math::GMPz; $s=Math::GMPz->new(0); forfactored { $s+=$_[0]; } 10**17,10**17+10**5; say $s'
251500000082904735996
user 0m0.991s

time perl -Mntheory=:all -E 'use Math::GMPz; $s=Math::GMPz->new(0); for (10**17 .. 10**17+10**5) { $s+=(factor($_))[0]; }; say $s'
251500000082904735996
user 0m1.002s
[/code]

The simple solution which I used is to find crossover points where it's faster to factor each n rather than sieve. The chunk size can be tuned some -- I use 8192, Pari uses 1024. I would have gone larger but I thought 4MB of internal data was enough.

The other idea I thought about (but have not implemented) was for large enough sizes to sieve to the cube root of b which leaves only semiprimes behind. Then a single pass fixes up the remainder.

Note that Pari/GP is far cleaner with bignums than Perl as can be seen above. It also is more efficient for the call itself. I have a lot of overhead on the block call. I thought about some different calling conventions to cut overhead but they are messier to use.

Example, OEIS 053462:
[code]
? a(n)=my(s); forsquarefree(d=1, sqrtint(n=10^n-1), s += n\d[1]^2 * moebius(d)); s
%4 = (n)->my(s);forsquarefree(d=1,sqrtint(n=10^n-1),s+=n\d[1]^2*moebius(d));s
? a(15)
time = 7,891 ms.
%5 = 607927101854103
[/code]
[code]
time perl -Mntheory=:all -E 'sub a { my($s,$n)=(0,10**$_[0]-1); forsquarefree { $s += int($n / ($_*$_)) * ((scalar(@_) & 1)?-1:1); } sqrtint($n); $s; } say a(15)'
607927101854103

real 0m8.048s
user 0m7.794s
sys 0m0.060s
[/code]

[QUOTE=danaj;479893]forfactored and forsquarefree. Pretty neat functions that run through a range and return the number and its factorization. I added this to Perl/ntheory, with both positive and negative results.[/QUOTE]

I asked for this feature* back in 2011, it's good that it was finally added. Unfortunately it's pretty slow; I guess beggars can't be choosers. But at least they did good work preparing for it before it was finally written, by allowing essentially all of the arithmetic functions to take [n, factor(n)] as an argument.

* forfactored, that is -- forsquarefree was all Karim.

[QUOTE=danaj;479893]2) The current Pari/GP implementation is limited to 2^64, and is also really slow near that limit. It's chunk size is 1024 and it's sieving to sqrt(b) for each chunk's a,b.[/QUOTE]

Ick. I don't often need forfactored in that range -- 10^19 [url=https://www.computerworld.com/article/2534312/operating-systems/the--640k--quote-won-t-go-away----but-did-gates-really-say-it-.html]ought to be enough for anybody[/url] -- but it would be great if it scaled better and even handled 128 bit numbers (or larger).

With a faster forfactored implementation I would run loops up to 10^12 probably daily (and forprime up to 10^13, but that's another issue, not quite as bad).

[QUOTE=danaj;479893]The simple solution which I used is to find crossover points where it's faster to factor each n rather than sieve. The chunk size can be tuned some -- I use 8192, Pari uses 1024. I would have gone larger but I thought 4MB of internal data was enough.[/QUOTE]

I'd be more interested in tuning the crossover point if the implementation was faster.

[QUOTE=danaj;479893]The other idea I thought about (but have not implemented) was for large enough sizes to sieve to the cube root of b which leaves only semiprimes behind. Then a single pass fixes up the remainder.[/QUOTE]

I can't really see how this would help much. If you implement it I'd love to hear how it goes -- I hope I'm wrong.

Raising the chunk size to 8192 in src/language/sumiter.c speeds up the 10^16 example from 34s to 4.5s and 10^17 from 107s to 14s. In the simple and perhaps more common case of 1 -> 10^8, it makes it go from 23s to 30s so actually slows it down some. Interesting.

For my program almost 75% of the time for the simple case is just the callback, vs. the actual work. For large values mine continues to work though it just just a convenience method for now.

My idea on sieving to the cube root is basically the same as a partial sieve. The initial sieve isn't so deep it takes crazy long, but efficiently removes all the small factors from the range. Most non-square-free inputs are caught here. The remainder is either prime (single test), a perfect square (easy check), or requires one factor call. It should be faster on a range than punting every value to the generic factor call, and faster than sieving deeply which is expensive.

It looks like raising the chunk size is an easy mitigation for that however, without all the complication.

[QUOTE=danaj;479912]My idea on sieving to the cube root is basically the same as a partial sieve. The initial sieve isn't so deep it takes crazy long, but efficiently removes all the small factors from the range. Most non-square-free inputs are caught here. The remainder is either prime (single test), a perfect square (easy check), or requires one factor call. It should be faster on a range than punting every value to the generic factor call, and faster than sieving deeply which is expensive.[/QUOTE]

It sounds interesting, but I'd have to code it up to see how fast it would be in practice. Checking if you have a square takes amortized time o(1), prime testing isn't too bad with BPSW, and then factoring... at least gp is pretty efficient at those (tiny) sizes.

I think it's worthwhile, though, since it's the same code you'd want to call for a short range in any case, so even if it's not good for normal sieving it would still have value.

The proof would be putting it in Pari/GP, inside src/basemath/ifactor1.c (vecfactoru_i and/or vecfactorsquarefreeu) then evaluating performance and complexity. Then extend to the negative cases if desired.

It was rather easy to add to my code, and I think it should be just about as easy for Pari. The main loop is identical, just up to the cube root instead of usqrt, then the "complete factorisation" loop gets the new code if we chose the cube root limit.

The perfect square test should be only needed for the square-free case, as otherwise we don't care and just want to factor it. One call to ifac_next I believe.

I found the crossover point to be around 10^14, but this is dependent on a lot of different code that would be different in Pari. As an example, forfactored from 10^17 to 10^17+10^6 was 9.0s factoring each, 52.5s sieve to sqrt, 5.9s sieve to cbrt. I'm using a x86_64 Montgomery Pollard-Rho-Brent, as discussed on [URL="http://www.mersenneforum.org/showthread.php?t=22525"]this thread[/URL], which is about 10x faster than Pari at this size, which would certainly change the results. In theory we could put that faster factoring code in Pari. If the people at Bordeaux would like to discuss it, I could find some time to live in France.... :)

Apologies for what I assume is a basic question, but I cannot for the life of me figure out the answer, even after reading up on the particular functions (and number structures PARI/GP uses).

If I do Mod(5,991), I understand, I think, that that's the equivalent of 5 % 991.

If I do sqrt(Mod(5,991)), I get Mod(63, 991) as the answer. How does that happen?

[QUOTE=wombatman;496786]If I do Mod(5,991), I understand, I think, that that's the equivalent of 5 % 991.

If I do sqrt(Mod(5,991)), I get Mod(63, 991) as the answer. How does that happen?[/QUOTE]

Because Mod(63,99)^2 = 63^2 %991 = Mod(5,991)

[QUOTE=wombatman;496786]Apologies for what I assume is a basic question, but I cannot for the life of me figure out the answer, even after reading up on the particular functions (and number structures PARI/GP uses).

If I do Mod(5,991), I understand, I think, that that's the equivalent of 5 % 991.

If I do sqrt(Mod(5,991)), I get Mod(63, 991) as the answer. How does that happen?[/QUOTE]
Most precisely, Mod(5,991) is of type t_INTMOD (integer modulo to a prime modulus) and 5%991 (which is 5) is of type t_INT (integer). So you can say that 5%991 is lift(Mod(5,991)).

Pari, however, will affirm that 5 == Mod(5,991) [comparison equal] and that Mod(5,991) - 5 is Mod(0,991). It obviously converts them to the same type to make the comparison or do the arithmetic.

Taking the square root of 5%991 gives the numerical approximation 2.236... -- here Pari converts 5 to real type; sqrt() is further defined for complex type, using the "principal branch."

According to the manual for my ancient version of Pari-GP, sqrt() is also defined for intmod a [i]prime[/i] modulus and for [i]p[/i]-adics. If you try sqrt(Mod(n,p)) where n is a quadratic non-residue modulo the prime p, you get an error message.

[QUOTE=Dr Sardonicus;496815]Pari, however, will affirm that 5 == Mod(5,991) [comparison equal] and that Mod(5,991) - 5 is Mod(0,991). It obviously converts them to the same type to make the comparison or do the arithmetic.[/QUOTE]

It's a little more subtle than that, since also 996 == Mod(5,991). :smile: In this case == is asserting that the integer is in the congruence class in question. The operator === is also available when this behavior is not wanted.

[QUOTE=Dr Sardonicus;496815]According to the manual for my ancient version of Pari-GP, sqrt() is also defined for intmod a [i]prime[/i] modulus and for [i]p[/i]-adics. If you try sqrt(Mod(n,p)) where n is a quadratic non-residue modulo the prime p, you get an error message.[/QUOTE]

I have the following terrible script for that use case:
[code]sqrts(x)={
if(type(x) != "t_INTMOD", error("must be intmod"));
if (x.mod == 1, return([0]));
my(n=lift(x), f=factor(x.mod), v, t);
v = Mod(sqrtmodpk(n, f[1,1], f[1,2]), f[1,1]^f[1,2]);
for(i=2,#f[,1],
t = sqrtmodpk(n, f[i,1], f[i,2]);
v=vecchinese(v, Mod(t, f[i,1]^f[i,2]))
);
vecsort(lift(v));
}
addhelp(sqrts, "sqrts(x): Returns a vector of the square roots of the intmod x.");


\\ Should be the same as select(m->m^2%p^k==n,[0..p^k-1])
sqrtmodpk(n,p,k)={
\\return(select(m->m^2%p^k==n,[0..p^k-1])); \\\\\\\\\\\\\\\\\\
if (!isprime(p), error("p must be a prime"));
if (k < 1 || type(k) != "t_INT", error("bad exponent"));
n %= p^k;
if(n==0,
return(vector(p^(k\2),i,i-1)*p^((k+1)\2))
);
if(n%p==0,
return(if(n%p^2, [], sqrtmodpk(n/p^2,p,k-1)*p))
);
my(r);
if (p == 2,
if (k < 3,
return(if (k == 1,
[n]
,
\\ k == 2
if(n > 1, [], if(n, [1,3], [0,2]))
))
);

if (n%8 != 1, return([]));

\\ TODO: Slow; should have modulus increase throughout calculation
\\ rather than starting at full p-adic precision.
my(m=1<<k);
r=Mod(1,m);
for(e=2,k,
my(z=3-n*Mod(r,m)^2);
r *= lift(z)/2
);
m /= 2;
r=lift(r*n)%m;
vecsort([r,m-r,m+r,2*m-r]%(2*m))
,
if (kronecker(n, p) == -1, return([]));
my(m=p,t);
r=lift(sqrt(Mod(n,p)));
for(e=2,k,
t = lift(Mod((n - r^2)/m, p)/(2*r));
r += t * m;
m *= p
);
vecsort([r,m-r])
)
};[/code]

Thanks everyone. That makes it a lot clearer. :smile: