Thank you for clarification about tilde. Lists, vectors and matrices are not the only types that could be bulky, as I can see.
But I don't usually use series and polynoms so I don't care :)
Anyway, I beleive there will be documentation for 2.12 when it's goes stable, so we users can read it.
Thanx for simplification.
I have looked through the code and discovered that the case n==0 && k==0 actualy never happens in main cycle.
So that first "if" could be moved out of main cycle, thus making further simplification.
I also noticed that numbpart() is called quite often, so using pre-populated vector witn numbpart()values saves a bunch of time.
I didn't catch, however, what values are actually used. So I populated all k values, though I belive only half of them are used.
For example, if we call p(1000,200), then numbpart(200) is calulated, but numbpart from 167 to 199 - are not used in recursion.
If we call p(1000,499), than numbpart 1 to 333 are used, and 334-499 are not.
But even with that overcalculation (all numbers of partitions from 1 to k), total running time of p() compared to p2() is about 6 (six) times faster:
p2(n,k) =
{
if(n==0 && k==0,return(0));
local(M = Map(),NP=vector(k,i,numbpart(i))); /* pre-calculate all numbparts */
my(memo_p = (N,K)-> my(res);
if(K<2,return(K));
if(mapisdefined(M,[N,K],&res),return(res));
res=if(K>=N,NP[N],self()(N,K-1)+self()(N-K,K));
mapput(M,[N,K],res);return(res));
memo_p(n,k);
print(#M)
}
I went farther and instead of using map, created zero-ed array to use as a map, thus not using mapisdefined() and mapput():
p3(n,k) =
{
if(n==0 && k==0,return(0)); /* that case can't happen later in rucursion, so returning instantly */
local(NP=vector(k,i,numbpart(i)),M=vector(n,i,vector(k))); /* pre-calculate all numbparts and create zeroed p(n,k) index array */
my(memo_p = (N,K)-> my(res);
if(K<2,return(K)); /* instead of if(K==0,return(0),if(K==1,return(1))) */
if(res=M[N][K],return(res)); /* instead of mapisdefined() Note here is = , not == */
res=if(K>=N,NP[N],self()(N,K-1)+self()(N-K,K)); /* use pre-calculated NP[N] instead of numbpart(N) */
M[N][K]=res; /* instead of mapput() */
return(res));
memo_p(n,k);
}
That gave about two-fold speed-up of p3() compared to p2(), making total speed-up of p3() to p() about 10 to 12 (twelve!) times...
By the way, I also wanted to ask about partitions()function.
Before going here into mailing list, I went through pari/gp source code of partitions(), and in there, you first calculate how many partitions there will be by generating all the partition without storing them, and after that create and populate output vector with partitions, generating them all again, for the second time.
Actually, I wanted to just steal calculation of partitions number algorithm (with limitation on parts number/part size) from pari/gp source code, but it seems not very effective...
That doesn't affect overall performance of partitions() as anyways calculation of partitions quantity is a fraction of time of generating them, but it would be good to have numbpart() function in pari/gp with limitation on parts size/quantity, like in partitions() generating function.
Regards,
Dmitry.
* Дмитрий Рыбас [2019-07-16 09:44]:
> So, for system functions pointers are marked with ampersand (&) and for
> user functions pointers are marked with tilde (~) in 2.12 ?
Not quite: we don't have true user function pointers (hence the tilde
instead of the customary ampersand). It's really a mutable *container*
(we're allowing to modify in place an object of list/vector/matrix-type),
this is not a side effect changing a variable value. Compare:
? f(~x) = x = 1;
? x = 0; f(~x); x \\ x didn't change
%2 = 0
? g(~x) = x[1] = 1;
? x = [0]; g(~x); x \\ but the *content* of a container type x would change
%4 = [1]
This can also be used to make sure that a given bulky argument is
*never* copied by user functions. GP's copy optimizer is usually clever
enough but it sometimes errs on the paranoid side to survive problematic
constructs (usually involving global variables). E.g.,
x = [1]; /* global; no lexical scoping */
f(y) = /*...complicated code...*/;
f(x)
Must must we make a deep copy of 'x' in case some user subroutine called by
f messes with it, or is a shallow copy safe enough ? [ In fact it's
rarely safe because GP is not compiled: global subroutines called by f may be
redefined at any time... ]
With
f(~y) = ...
f(~x)
there's no problem: no copy.
> That's great! But a little bit strange if that difference really exist.
> Will be waiting for 2.12 stable.
>
> As for "cumbersome" way, here is what I has written:
>
> p(n,k) =
> {
> local(M = Map()); /* Do not use my() here, use local() !! */
> my(memo_p = (N,K)-> my(res);
> if(N==0 && K==0,return(1));
> if(k<1,return(0));
^--- you mean K here
> if(mapisdefined(M,[N,K],&res),return(res));
> if(K>=N,res=numbpart(N);mapput(M,[N,K],res);return(res));
> res=self()(N,K-1)+self()(N-K,K);
> mapput(M,[N,K],res);
> return(res));
> memo_p(n,k);
> }
You can simplify this a little:
p(n,k) =
{
local(M = Map()); /* yes, local is needed: sorry... */
my(memo_p = (N,K)-> my(res);
if(N==0 && K==0,return(1));
if(K<1,return(0));
if(mapisdefined(M,[N,K],&res),return(res));
res = if(K>=N, numbpart(N), self()(N,K-1)+self()(N-K,K));
mapput(M,[N,K],res); return(res));
memo_p(n,k);
}
Cheers,
K.B.
--
Karim Belabas, IMB (UMR 5251) Tel: (+33) (0)5 40 00 26 17
Universite de Bordeaux Fax: (+33) (0)5 40 00 21 23
351, cours de la Liberation http://www.math.u-bordeaux.fr/~kbelabas/
F-33405 Talence (France) http://pari.math.u-bordeaux.fr/ [PARI/GP]
`