Bill Allombert on Wed, 05 Feb 2014 11:49:21 +0100

 Re: handling several substacks ?

• To: pari-dev@pari.math.u-bordeaux.fr
• Subject: Re: handling several substacks ?
• From: Bill Allombert <allomber@math.u-bordeaux.fr>
• Date: Wed, 5 Feb 2014 11:49:19 +0100
• Delivery-date: Wed, 05 Feb 2014 11:49:21 +0100
• References: <CAAuMLAkkKw564_wW=e4LssfMQF5NQ+AZUNLoqUM2gT2EVK7c1A@mail.gmail.com>
• User-agent: Mutt/1.5.21 (2010-09-15)

```On Wed, Feb 05, 2014 at 10:52:13AM +0100, Pascal Molin wrote:
> I have to deal with the following kind of loops:
>   {
>     /* main loop */
>     pari_sp av1 = avma, lim1 = stack_lim(avma, 1);
>     long n;
>     for (n = 1; n <= N; ++n)
>     {
>         long i, j, k;
>         i = random_Fl(L)+1;
>         j = random_Fl(L)+1;
>           for (k = 1; k <= K; ++k)
>           {
>             gel(gel(a, i), k) = gmul( gel(gel(a, j), k), gel(gel(m, i), k));
>             if (avma > lim1)
>               a = gerepilecopy(av1, a);

Hello Pascal,

There is two issues: First
if (avma > lim1)
is incorrect: the inequality is in the wrong direction.
The recommended idiom is
if (low_stack(lim, stack_lim(av,1)))
which does
if (avma < lim1)

Secondly, you should only need to gerepile after the
for (k = 1; k <= K; ++k)
is done.

Thirdly, it is not fair to compare GP and C with the same amount of stack,
because GP allocates a lot of memory outside the stack.

Try the new C file in attachment.

? gploop()
time = 2,280 ms.
%4 = 0
? cloop()
time = 1,028 ms.
%5 = 0

Cheers,
Bill.
```
```/*-*- compile-command: "/usr/local/bin/gcc-4.7 -c -o stack.o -O3 -Wall -fno-strict-aliasing -fomit-frame-pointer -fPIC -I"/Users/pascal/devel/install/include" stack.gp.c && /usr/local/bin/gcc-4.7 -o stack.gp.so -bundle -flat_namespace -undefined suppress -O3 -Wall -fno-strict-aliasing -fomit-frame-pointer -fPIC stack.o "; -*-*/
#include <pari/pari.h>
/*
GP;install("cloop","D50000,L,D1000,L,D8,L,p","cloop","./stack.so");
*/
GEN
cloop(long N, long K, long L, long prec)
{
pari_sp ltop = avma;
GEN a, m;
{
/* init */
long l;
a = cgetg(L+1, t_VEC);
m = cgetg(L+1, t_VEC);
for (l = 1; l <= L; ++l)
{
long k;
GEN p1, p2;
p1 = cgetg(K+1, t_VEC);
p2 = cgetg(K+1, t_VEC);
for (k = 1; k <= K; ++k)
{
gel(p1, k) = gexp(gmul(gen_I(), genrand(real_1(prec))), prec);
gel(p2, k) = gexp(gmul(gen_I(), genrand(real_1(prec))), prec);
}
gel(a, l) = p1;
gel(m, l) = p2;
}
}
{
/* main loop */
pari_sp av1 = avma, lim1 = stack_lim(avma, 1);
long n;
for (n = 1; n <= N; ++n)
{
long i, j, k;
i = random_Fl(L)+1;
j = random_Fl(L)+1;
for (k = 1; k <= K; ++k)
gel(gel(a, i), k) = gmul( gel(gel(a, j), k), gel(gel(m, i), k));
if (low_stack(lim1, stack_lim(avma, 1)))
{
if (DEBUGLEVEL>=2) pari_warn(warnmem,"cloop: %ld\n",n);
a = gerepilecopy(av1, a);
}
}
}
/* use A = a[i] here */
avma = ltop;
return gen_0;
}
```