Re: handling several substacks ?

[Bill Allombert <allomber@math.u-bordeaux.fr>]

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