| Bill Allombert on Wed, 14 Jun 2023 21:34:48 +0200 |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
| Re: Give up with MT, need help!!! |
On Wed, Jun 14, 2023 at 05:56:38PM +0200, Jean-Luc ARNAUD wrote:
> Hi all,
>
> As a newby to PARI/GP, I'm trying to understand how it works in MT
> environment.
> And I'm getting a lot of difficulties.
>
> As an example, how could I replace the for ... loop by a parfor ... one in
> the chudnovsky_parallel(n) function ?
>
> After many and many tries, getting always either "Please export ..." or
> "Impossible div ..." error message, I give up ...
>
> Would somebody be so kind as to modify the below script in order to use
> parfor loop instead of for?
>
> |chudnovsky_parallel(n) = {||
> || my(k, res, Result);||
> ||
> || for(k = 0, ceil(n/14),||
> || Result=pareval([||
> || ()-> (-1)^k *(6*k)! * (13591409 + 545140134*k),||
> || ()-> (3*k)! *k!^3 * (640320^(3*k + 3 / 2))]); /* +3/2 -> 583 ms,
> +1.5 -> 18,536 s !!! */||
> || res += Result[1] / Result[2];||
> || );||
> ||
> || res = 1 / (res * 12);||
> ||
> ||};||
You can do this (using parfor)
chudnovsky_parsum(n) =
{
my(res);
parfor(k = 0, ceil(n/14),
my(d,n);
d = (-1)^k *(6*k)! * (13591409 + 545140134*k);
n = (3*k)! *k!^3 * (640320^(3*k + 3 / 2));
d/n,
dn,
res += dn
);
res = 1 / (res * 12);
}
You could also use parsum:
chudnovsky_parsum(n) =
{
my(res);
res = parsum(k = 0, ceil(n/14),
my(d,n);
d = (-1)^k *(6*k)! * (13591409 + 545140134*k);
n = (3*k)! *k!^3 * (640320^(3*k + 3 / 2));
d/n
);
res = 1 / (res * 12);
}
or parapply
chudnovsky_parapply(n) =
{
my(res);
res = vecsum(parapply(k->
my(d,n);
d = (-1)^k *(6*k)! * (13591409 + 545140134*k);
n = (3*k)! *k!^3 * (640320^(3*k + 3 / 2));
d/n,
[0..ceil(n/14)]
));
res = 1 / (res * 12);
}
But this is not the right way to compute this. You should use
the recurrence formula for the factorial
k! = k*(k-1)! and not recompute k!^3, (3*k)!, (6*k)! for all k from scratch.
(3*k)!= (3*k)*(3*k-1)*(3*k-2)*(3*(k-1))!
etc.
Also you should factor out 640320^(3/2) to avoid recomputing it.
Cheers,
Bill