Alberto Zanoni on Thu, 05 Mar 2020 20:06:32 +0100


[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]

Help request for segmentation fault


Hello,
         I ask for help, as this is my first time using gp2c. I
generated the below reported file (with the main function at the end
manually added) with gp2c by a PARI/GP script I wrote, and I compiled
it with

gcc -g -O3 -Wall df8.c -o df8.exe -lpari

When I try to execute it with

./df8.exe

I obtain


Confronto tempi per fattoriale doppio.
t1: ricorsione binaria  t2: Zanoni (ricorsione 7)       t3: Zanoni
(ricorsione 6)
-----------------------------------------------------
N       Check   t1      t2      t3      t2/t1   t3/t1   t3/t2
  ***   bug in PARI/GP (Segmentation Fault), please report.
  ***   Error in the PARI system. End of program.


I can't figure out what's going wrong nor how to correct it. Suggestions ?

Thank you very much in advance.

Alberto
------------------------------------

#include <pari/pari.h>
/*
GP;install("init_df8","vp","init_df8","./df8.gp.so");
GP;install("ss","D0,G,D0,G,","ss","./df8.gp.so");
GP;install("ff","D0,G,D0,G,D0,G,","ff","./df8.gp.so");
GP;install("tt","D0,G,D0,G,D0,G,","tt","./df8.gp.so");
GP;install("prod8","D0,G,D0,G,D0,G,D0,G,","prod8","./df8.gp.so");
GP;install("fattorialeDop","D0,G,D0,G,","fattorialeDop","./df8.gp.so");
GP;install("fattorialeDoppio","D0,G,","fattorialeDoppio","./df8.gp.so");
GP;install("computeN","D0,G,","computeN","./df8.gp.so");
GP;install("prodOdd","D0,G,","prodOdd","./df8.gp.so");
GP;install("computeProductsFromTo","D0,G,D0,G,D0,G,","computeProductsFromTo","./df8.gp.so");
GP;install("computeProducts","D0,G,","computeProducts","./df8.gp.so");
GP;install("prodOdd2","D0,G,","prodOdd2","./df8.gp.so");
GP;install("prodOdd3","D0,G,","prodOdd3","./df8.gp.so");
GP;install("tempi","D0,G,DGp","tempi","./df8.gp.so");
GP;install("df8test","lp","df8test","./df8.gp.so");
*/
void init_df8(long prec);
GEN ss(GEN n, GEN i);
GEN ff(GEN n, GEN a, GEN k);
GEN tt(GEN n, GEN a, GEN k);
GEN prod8(GEN nSqr, GEN aSqr, GEN nSqrAsqrShifted8, GEN subtrahend);
GEN fattorialeDop(GEN start, GEN end);
GEN fattorialeDoppio(GEN N);
GEN computeN(GEN N);
GEN prodOdd(GEN N);
GEN computeProductsFromTo(GEN v, GEN start, GEN end);
GEN computeProducts(GEN v);
GEN prodOdd2(GEN N);
GEN prodOdd3(GEN N);
GEN tempi(GEN N, GEN volte, long prec);
long df8test(long prec);
/*End of prototype*/

static GEN NUMERO_ESEMPI;
static GEN START;
static GEN PRODOTTI_FATTORI_DISPARI;
/*End of global vars*/

void
init_df8(long prec)	  /* void */
{
  GEN p1;	  /* vec */
  NUMERO_ESEMPI = pol_x(fetch_user_var("NUMERO_ESEMPI"));
  START = pol_x(fetch_user_var("START"));
  PRODOTTI_FATTORI_DISPARI = pol_x(fetch_user_var("PRODOTTI_FATTORI_DISPARI"));
  NUMERO_ESEMPI = stoi(60);
  START = addis(mulsi(4, powiu(stoi(10), 2)), 1);
  p1 = cgetg(9, t_VEC);
  gel(p1, 1) = gen_1;
  gel(p1, 2) = stoi(3);
  gel(p1, 3) = stoi(15);
  gel(p1, 4) = stoi(105);
  gel(p1, 5) = stoi(945);
  gel(p1, 6) = stoi(10395);
  gel(p1, 7) = stoi(135135);
  gel(p1, 8) = stoi(2027025);
  /*\\\\\\\\\\\\\\\\ */
  /* Calcolo del fattoriale doppio dispari (prodotto dei numeri dispari */
  /* da N "in giù"). */
  /* PRIMA VERSIONE: DA NON USARE: PIÙ EFFICIENTE LA prodOdd2 ! */

  PRODOTTI_FATTORI_DISPARI = p1;
  df8test(prec);
  return;
}

/* \p 5 */
/* allocatemem(); */
/* allocatemem(); */
/* allocatemem(); */
/* allocatemem(); */

GEN
ss(GEN n, GEN i)
{
  return gmul(gmul(gmul(gsub(gsqr(n), gsqr(gaddgs(gmulsg(8, i), 1))),
gsub(gsqr(n), gsqr(gaddgs(gmulsg(8, i), 3)))), gsub(gsqr(n),
gsqr(gaddgs(gmulsg(8, i), 5)))), gsub(gsqr(n), gsqr(gaddgs(gmulsg(8,
i), 7))));
}

GEN
ff(GEN n, GEN a, GEN k)
{
  return gsub(gsqr(gsub(gsub(gsqr(n), gsqr(a)), gsqr(gaddgs(gmulsg(2,
k), 1)))), gsqr(gmul(gmulsg(2, gaddgs(gmulsg(2, k), 1)), a)));
}

GEN
tt(GEN n, GEN a, GEN k)
{
  return gsub(gsqr(gaddgs(gadd(gsub(gsqr(n), gsqr(a)), gmulsg(4,
gsqr(gaddgs(k, 1)))), 1)), gmulsg(4, gadd(gmul(gmulsg(4,
gaddgs(gsqr(n), 1)), gsqr(gaddgs(k, 1))), gsqr(n))));
}

/* ff(n,a,k) = tt(n,a,k)^2 - (16*(k+1)*n*a)^2; */
/* t0(n,a) = (n^2-a^2+5)^2 - 4*(5*n^2+4); */

/* print("\n** ",ff(n,a,0)*ff(n,a,1) - (tt(n,a,0)^2 - (16*((0+1)*n*a))^2)); */

/*\\\\\\\\\\\\\\\\\\\\\\\\\\\ */
/* Calcola, con a = 8*i+4, il seguente prodotto. */
/* I primi due valori passati come argomento sono n^2 e a^2. */
/*  */
/*   (n-(8*i+7))(n-(8*i+5))(n-(8*i+3))(n-(8*i+1))(n+(8*i+1))(n+(8*i+3))(n+(8*i+5))(n+(8*i+7))
= */
/* = (n^-(8*i+7)^2)(n^2-(8*i+5)^2)(n^2-(8*i+3)^2)(n^2-(8*i+1)^2) */
/* = (n^-(a-3)^2)(n^2-(a-1)^2)(n^2-(a+1)^2)(n^2-(a+3)^2) */
/* */
/* = ((n^2-a^2+5)^2-4*(5*n^2+4))^2 - 256*n^2*a^2  <= Algoritmo di calcolo */

GEN
prod8(GEN nSqr, GEN aSqr, GEN nSqrAsqrShifted8, GEN subtrahend)
{
  GEN s, t, res;
  s = gaddgs(gsub(nSqr, aSqr), 5);
  t = gsub(gsqr(s), subtrahend);
  res = gsub(gsqr(t), nSqrAsqrShifted8);
  return res;
}

/*\\\\\\\\\\\\\\\\ */
/* Funzioni "naif" per confronto: calcolo del fattoriale doppio */
/* (per argomenti N dispari) con ricorsione binaria. */

/* Questa calcola i prodotti dei numeri dispari da start ad end,
ricorsivamente. */
GEN
fattorialeDop(GEN start, GEN end)
{
  GEN tmp = gen_0, res;
  res = gcopy(start);
  if (!gequal(end, start))
  {
    if (gequalgs(gsub(end, start), 2))
      res = gmul(res, end);
    else
    {
      /* else */
      tmp = gshift(gadd(end, start), -1);
      if (gequal0(gmodgs(tmp, 2)))
        tmp = gaddgs(tmp, 1);
      /* print("\ttmp = ",tmp); */
      res = gmul(fattorialeDop(start, tmp), fattorialeDop(gaddgs(tmp, 2), end));
    }
  }
  return res;
}

/* ...e questa usa la precedente. */
GEN
fattorialeDoppio(GEN N)
{
  GEN tmp = gen_0, res = gen_1;
  if (gcmpgs(N, 5) < 0)
    res = gcopy(N);
  else
  {
    tmp = gshift(gaddgs(N, 1), -1);
    if (gequal0(gmodgs(tmp, 2)))
      tmp = gaddgs(tmp, 1);
    res = gmul(fattorialeDop(stoi(3), tmp), fattorialeDop(gaddgs(tmp, 2), N));
  }
  return res;
}

/*\\\\\\\\\\\\\\\\ */
/* Calcola l'n giusto per la formula. */
GEN
computeN(GEN N)	  /* vec */
{
  GEN quozRem;
  GEN p1;	  /* vec */
  quozRem = divrem(gaddgs(N, 1), stoi(16), -1);
  p1 = cgetg(3, t_VEC);
  gel(p1, 1) = gadd(gmulsg(8, gel(quozRem, 1)), gel(quozRem, 2));
  gel(p1, 2) = gcopy(gel(quozRem, 1));
  return p1;
}

/* (a+8)^2         = a^2 + 16a + 64 = a^2 +(a+4)<<4 */
/* (n^2(a+8)^2)<<8 = (n^2(a^2 +(a+4)<<4))<<8 = n^2a^2<<8 + n^2(a+4)<<12 */
/* (a+8)<<4 = (a<<4) + 1<<5; */

GEN
prodOdd(GEN N)
{
  GEN n = gen_0, nSqr = gen_0, aSqr = stoi(16), aFourth = stoi(128),
subtrahend = gen_0, addend = gen_0, nSqrShifted12 = gen_0,
nSqrAsqrShifted8 = gen_0, tmp = gen_0, index, res = gen_0, p1 = gen_0,
nSqrShifted12TimesA = pol_x(fetch_user_var("nSqrShifted12TimesA")),
p2;
  /* a = 8*i+4;   i+1 => a' = 8(i+1)+4 = (8i+4) + 8 = a+8 */
  index = gbitand(gaddgs(N, 1), stoi(15));
  if (gequal0(index))
    p1 = gen_1;
  else
    p1 = gcopy(gel(PRODOTTI_FATTORI_DISPARI, itos(gshift(index, -1))));
  res = p1;
  /* Preparazione costanti di base, in modo che nel ciclo principale
ci siano solo somme. */
  tmp = computeN(N);
  n = gcopy(gel(tmp, 1));
  nSqr = gsqr(n);
  nSqrShifted12 = gshift(nSqr, 12);
  nSqrAsqrShifted8 = nSqrShifted12;
  subtrahend = gshift(gaddgs(gmulsg(5, nSqr), 4), 2);
  addend = gshift(nSqrShifted12, 3);
  nSqrShifted12TimesA = addend;
  /* Ciclo principale di prodotti. */
  p2 = gsubgs(gel(tmp, 2), 1);
  {
    GEN i;
    for (i = gen_0; gcmp(i, p2) <= 0; i = gaddgs(i, 1))
    {
      res = gmul(res, prod8(nSqr, aSqr, nSqrAsqrShifted8, subtrahend));
      aSqr = gadd(aSqr, aFourth);
      aFourth = gaddgs(aFourth, 128);
      nSqrAsqrShifted8 = gadd(nSqrAsqrShifted8, nSqrShifted12TimesA);
      nSqrShifted12TimesA = gadd(nSqrShifted12TimesA, addend);
    }
  }
  return res;
}

/*\\\\\\\\\\\\\\\\\\\\\\\\ */
/* Calcola il prodotti degli elementi del vettore contenuti */
/* nelle posizioni (indici) da start a end. */
GEN
computeProductsFromTo(GEN v, GEN start, GEN end)
{
  GEN res, LIM = stoi(4), half = gen_0;
  res = gcopy(gel(v, gtos(start)));
  if (gcmp(gsub(end, start), LIM) < 0)
  {
    GEN p1;
    p1 = gaddgs(start, 1);
    {
      GEN i;
      for (i = p1; gcmp(i, end) <= 0; i = gaddgs(i, 1))
        res = gmul(res, gel(v, gtos(i)));
    }
  }
  else
  {
    /* else */
    half = gshift(gadd(start, end), -1);
    res = gmul(computeProductsFromTo(v, start, half),
computeProductsFromTo(v, gaddgs(half, 1), end));
  }
  return res;
}

/*\\\\\\\\\\\\\\\\\\\\\\\\ */
/* Calcola il prodotti degli elementi del vettore */
GEN
computeProducts(GEN v)
{
  GEN lim, res;
  lim = stoi(glength(v)>>1);
  res = gcopy(gel(v, 1));
  if (glength(v) > 1)
    res = gmul(computeProductsFromTo(v, gen_1, lim),
computeProductsFromTo(v, gaddgs(lim, 1), stoi(glength(v))));
  return res;
}

/*\\\\\\\\\\\\\\\\\\\\\\\\ Più efficiente. */
/* Fattoriale doppio. N!! = N(N-2)(N-4)···5*3*1 */
GEN
prodOdd2(GEN N)
{
  GEN n = gen_0, nSqr = gen_0, a, aSqr = stoi(16), aFourth =
stoi(128), subtrahend = gen_0, addend = gen_0, nSqrShifted12 = gen_0,
nSqrAsqrShifted8 = gen_0, tmp = gen_0, incrA, resVec = gen_0, index,
res = gen_0, p1 = gen_0, nSqrShifted12TimesA =
pol_x(fetch_user_var("nSqrShifted12TimesA")), p2;
  GEN p3;	  /* vec */
  GEN p4;
  /* a = 8*i+4;   i+1 => a' = 8(i+1)+4 = (8i+4) + 8 = a+8 */
  a = shifti(stoi(-3), 31);
  incrA = shifti(gen_1, 31);
  index = gbitand(gaddgs(N, 1), stoi(15));
  if (gequal0(index))
    p1 = gen_1;
  else
    p1 = gcopy(gel(PRODOTTI_FATTORI_DISPARI, itos(gshift(index, -1))));
  res = p1;
  tmp = computeN(N);
  n = gcopy(gel(tmp, 1));
  nSqr = gsqr(n);
  subtrahend = gshift(gaddgs(gmulsg(5, nSqr), 4), 2);
  nSqrShifted12 = gshift(nSqr, 12);
  nSqrAsqrShifted8 = nSqrShifted12;
  addend = gshift(nSqrShifted12, 3);
  nSqrShifted12TimesA = addend;
  p2 = gcopy(gel(tmp, 2));
  {
    long l5;
    p3 = cgetg(gtos(p2)+1, t_VEC);
    for (l5 = 1; gcmpsg(l5, p2) <= 0; ++l5)
      gel(p3, l5) = gen_0;
  }
  resVec = p3;
  /* Vettore per organizzare i risultati. */

  /* Ciclo di prodotti: memorizzazione dei risultati parziali in un vettore. */
  p4 = gcopy(gel(tmp, 2));
  {
    GEN i;
    for (i = gen_1; gcmp(i, p4) <= 0; i = gaddgs(i, 1))
    {
      if (gcmpgs(i, 7) > 0)
        gel(resVec, gtos(i)) = gadd(gmulsg(7,
gadd(gadd(gsub(gsub(gel(resVec, gtos(gsubgs(i, 1))), gel(resVec,
gtos(gsubgs(i, 6)))), gmulsg(3, gsub(gel(resVec, gtos(gsubgs(i, 2))),
gel(resVec, gtos(gsubgs(i, 5)))))), gmulsg(5, gsub(gel(resVec,
gtos(gsubgs(i, 3))), gel(resVec, gtos(gsubgs(i, 4)))))), gmulsg(45,
a))), gel(resVec, gtos(gsubgs(i, 7))));
      else
      {
        gel(/*else */
        resVec, gtos(i)) = prod8(nSqr, aSqr, nSqrAsqrShifted8, subtrahend);
        aSqr = gadd(aSqr, aFourth);
        aFourth = gaddgs(aFourth, 128);
        nSqrAsqrShifted8 = gadd(nSqrAsqrShifted8, nSqrShifted12TimesA);
        nSqrShifted12TimesA = gadd(nSqrShifted12TimesA, addend);
      }
      a = gadd(a, incrA);
    }
  }
  if (!gequal0(gel(tmp, 2)))
    /* Calcolo dei prodotti nel vettore. */
    res = gmul(res, computeProducts(resVec));
  return res;
}

/*\\\\\\\\\\\\\\\\\\\\\\\\ Ancora (un po') più efficiente. */
/* Fattoriale doppio. N!! = N(N-2)(N-4)···5*3*1 */
/* addend2(aa) = -7*aa^2 - 336*aa - 4251 */
/* addend2(aa+8)-addend2(aa) = -112*aa - 3136    ;     addend(4) = -5707 */

GEN
prodOdd3(GEN N)
{
  GEN n = gen_0, nSqr = gen_0, aSqr = stoi(16), aFourth = stoi(128),
subtrahend = gen_0, addend = gen_0, nSqrAsqrShifted8 = gen_0, tmp =
gen_0, addend2 = gen_0, nSqrTimes45 = gen_0, resVec = gen_0, index,
incr, step, res = gen_0, p1 = gen_0, nSqrShifted12TimesA =
pol_x(fetch_user_var("nSqrShifted12TimesA")), p2;
  GEN p3;	  /* vec */
  GEN p4;
  /* a = 8*i+4;   i+1 => a' = 8(i+1)+4 = (8i+4) + 8 = a+8 */
  index = gbitand(gaddgs(N, 1), stoi(15));
  incr = shifti(stoi(945), 31);
  step = shifti(stoi(315), 31);
  if (gequal0(index))
    p1 = gen_1;
  else
    p1 = gcopy(gel(PRODOTTI_FATTORI_DISPARI, itos(gshift(index, -1))));
  res = p1;
  tmp = computeN(N);
  n = gcopy(gel(tmp, 1));
  nSqr = gsqr(n);
  subtrahend = gshift(gaddgs(gmulsg(5, nSqr), 4), 2);
  nSqrTimes45 = gmulgs(nSqr, 45);
  addend2 = gshift(gsubsg(256815, nSqrTimes45), 24);
  nSqrAsqrShifted8 = gshift(nSqr, 12);
  addend = gshift(nSqr, 15);
  nSqrShifted12TimesA = addend;
  p2 = gcopy(gel(tmp, 2));
  {
    long l5;
    p3 = cgetg(gtos(p2)+1, t_VEC);
    for (l5 = 1; gcmpsg(l5, p2) <= 0; ++l5)
      gel(p3, l5) = gen_0;
  }
  resVec = p3;
  /* Vettore per organizzare i risultati. */

  /* Ciclo di prodotti: memorizzazione dei risultati parziali in un vettore. */
  p4 = gcopy(gel(tmp, 2));
  {
    GEN i;
    for (i = gen_1; gcmp(i, p4) <= 0; i = gaddgs(i, 1))
    {
      if (gcmpgs(i, 6) > 0)
      {
        gel(resVec, gtos(i)) = gadd(gsub(gadd(gsub(gmulsg(6,
gadd(gel(resVec, gtos(gsubgs(i, 1))), gel(resVec, gtos(gsubgs(i,
5))))), gmulsg(15, gadd(gel(resVec, gtos(gsubgs(i, 2))), gel(resVec,
gtos(gsubgs(i, 4)))))), gmulsg(20, gel(resVec, gtos(gsubgs(i, 3))))),
gel(resVec, gtos(gsubgs(i, 6)))), addend2);
        addend2 = gadd(addend2, incr = gadd(incr, step));
      }
      else
      {
        gel(/*else */
        resVec, gtos(i)) = prod8(nSqr, aSqr, nSqrAsqrShifted8, subtrahend);
        aSqr = gadd(aSqr, aFourth);
        aFourth = gaddgs(aFourth, 128);
        nSqrAsqrShifted8 = gadd(nSqrAsqrShifted8, nSqrShifted12TimesA);
        nSqrShifted12TimesA = gadd(nSqrShifted12TimesA, addend);
      }
    }
  }
  if (gcmpgs(gel(tmp, 2), 6) > 0)
  {
    long l6;
    l6 = gtos(gel(tmp, 2));
    gel(/* Prodotto parziale nel fattore più piccolo: invece */
    /* di una moltiplicazione finale lunga ne si fa una */
    /* all'inizio, molto più piccola. */
    resVec, l6) = gmul(gel(resVec, l6), res);
    /* Calcolo dei prodotti nel vettore. */
    res = computeProducts(resVec);
  }
  return res;
}

/*\\\\\\\\\\\\ */

GEN
tempi(GEN N, GEN volte, long prec)
{
  GEN fatt = gen_1, r2 = gen_0, r3 = gen_0, t;
  GEN p1;	  /* vec */
  if (!volte)
    volte = gen_1;
  {
    long l2;
    p1 = cgetg(4, t_VEC);
    for (l2 = 1; l2 <= 3; ++l2)
      gel(p1, l2) = gen_0;
  }
  t = p1;
  gel(t, 1) = stoi(gettime());
  {
    GEN j;
    for (j = gen_1; gcmp(j, volte) <= 0; j = gaddgs(j, 1))
      fatt = fattorialeDoppio(N);
  }
  gel(t, 1) = gmaxgs(gsubsg(gettime(), gel(t, 1)), 1);
  gel(t, 2) = stoi(gettime());
  {
    GEN j;
    for (j = gen_1; gcmp(j, volte) <= 0; j = gaddgs(j, 1))
      r2 = prodOdd2(N);
  }
  gel(t, 2) = gmaxgs(gsubsg(gettime(), gel(t, 2)), 1);
  gel(t, 3) = stoi(gettime());
  {
    GEN j;
    for (j = gen_1; gcmp(j, volte) <= 0; j = gaddgs(j, 1))
      r3 = prodOdd3(N);
  }
  gel(t, 3) = gsubsg(gettime(), gel(t, 3));
  pari_printf("%Ps\t%ld\t%Ps\t%Ps\t%Ps\t%Ps %Ps %%\t%Ps %%\n", N,
gequal(fatt, r2) && gequal(fatt, r3), gel(t, 1), gel(t, 2), gel(t, 3),
gdiv(gmul(stor(100, prec), gel(t, 2)), gel(t, 1)), gdiv(gmul(stor(100,
prec), gel(t, 3)), gel(t, 1)), gdiv(gmul(stor(100, prec), gel(t, 3)),
gel(t, 2)));
  return r2;
}

/*\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */
/* tempi(N, numero di ripetizioni del test) */
long
df8test(long prec)
{
  GEN end = pol_x(fetch_user_var("end"));
  pari_printf("Confronto tempi per fattoriale doppio.\nt1: ricorsione
binaria\tt2: Zanoni (ricorsione 7)\tt3: Zanoni (ricorsione 6)\n");
  pari_printf("-----------------------------------------------------\n");
  pari_printf("N\tCheck\tt1\tt2\tt3\tt2/t1\tt3/t1\tt3/t2\n");
  end = gadd(START, NUMERO_ESEMPI);
  {
    GEN h;
    for (h = START; gcmp(h, end) <= 0; h = gaddgs(h, 2))
      tempi(h, gen_1, prec);
  }
  /* print(h,"!! = ",FZ) */ return 0;
}

///////////////////////
// main function.

int main()
{
 pari_init(1e9,2);

 df8test(0);
 pari_close();

 return 0;
}
#include <pari/pari.h>
/*
GP;install("init_df8","vp","init_df8","./df8.gp.so");
GP;install("ss","D0,G,D0,G,","ss","./df8.gp.so");
GP;install("ff","D0,G,D0,G,D0,G,","ff","./df8.gp.so");
GP;install("tt","D0,G,D0,G,D0,G,","tt","./df8.gp.so");
GP;install("prod8","D0,G,D0,G,D0,G,D0,G,","prod8","./df8.gp.so");
GP;install("fattorialeDop","D0,G,D0,G,","fattorialeDop","./df8.gp.so");
GP;install("fattorialeDoppio","D0,G,","fattorialeDoppio","./df8.gp.so");
GP;install("computeN","D0,G,","computeN","./df8.gp.so");
GP;install("prodOdd","D0,G,","prodOdd","./df8.gp.so");
GP;install("computeProductsFromTo","D0,G,D0,G,D0,G,","computeProductsFromTo","./df8.gp.so");
GP;install("computeProducts","D0,G,","computeProducts","./df8.gp.so");
GP;install("prodOdd2","D0,G,","prodOdd2","./df8.gp.so");
GP;install("prodOdd3","D0,G,","prodOdd3","./df8.gp.so");
GP;install("tempi","D0,G,DGp","tempi","./df8.gp.so");
GP;install("df8test","lp","df8test","./df8.gp.so");
*/
void init_df8(long prec);
GEN ss(GEN n, GEN i);
GEN ff(GEN n, GEN a, GEN k);
GEN tt(GEN n, GEN a, GEN k);
GEN prod8(GEN nSqr, GEN aSqr, GEN nSqrAsqrShifted8, GEN subtrahend);
GEN fattorialeDop(GEN start, GEN end);
GEN fattorialeDoppio(GEN N);
GEN computeN(GEN N);
GEN prodOdd(GEN N);
GEN computeProductsFromTo(GEN v, GEN start, GEN end);
GEN computeProducts(GEN v);
GEN prodOdd2(GEN N);
GEN prodOdd3(GEN N);
GEN tempi(GEN N, GEN volte, long prec);
long df8test(long prec);
/*End of prototype*/

static GEN NUMERO_ESEMPI;
static GEN START;
static GEN PRODOTTI_FATTORI_DISPARI;
/*End of global vars*/

void
init_df8(long prec)	  /* void */
{
  GEN p1;	  /* vec */
  NUMERO_ESEMPI = pol_x(fetch_user_var("NUMERO_ESEMPI"));
  START = pol_x(fetch_user_var("START"));
  PRODOTTI_FATTORI_DISPARI = pol_x(fetch_user_var("PRODOTTI_FATTORI_DISPARI"));
  NUMERO_ESEMPI = stoi(60);
  START = addis(mulsi(4, powiu(stoi(10), 2)), 1);
  p1 = cgetg(9, t_VEC);
  gel(p1, 1) = gen_1;
  gel(p1, 2) = stoi(3);
  gel(p1, 3) = stoi(15);
  gel(p1, 4) = stoi(105);
  gel(p1, 5) = stoi(945);
  gel(p1, 6) = stoi(10395);
  gel(p1, 7) = stoi(135135);
  gel(p1, 8) = stoi(2027025);
  /*\\\\\\\\\\\\\\\\ */
  /* Calcolo del fattoriale doppio dispari (prodotto dei numeri dispari */
  /* da N "in giù"). */
  /* PRIMA VERSIONE: DA NON USARE: PIÙ EFFICIENTE LA prodOdd2 ! */
  
  PRODOTTI_FATTORI_DISPARI = p1;
  df8test(prec);
  return;
}

/* \p 5 */
/* allocatemem(); */
/* allocatemem(); */
/* allocatemem(); */
/* allocatemem(); */

GEN
ss(GEN n, GEN i)
{
  return gmul(gmul(gmul(gsub(gsqr(n), gsqr(gaddgs(gmulsg(8, i), 1))), gsub(gsqr(n), gsqr(gaddgs(gmulsg(8, i), 3)))), gsub(gsqr(n), gsqr(gaddgs(gmulsg(8, i), 5)))), gsub(gsqr(n), gsqr(gaddgs(gmulsg(8, i), 7))));
}

GEN
ff(GEN n, GEN a, GEN k)
{
  return gsub(gsqr(gsub(gsub(gsqr(n), gsqr(a)), gsqr(gaddgs(gmulsg(2, k), 1)))), gsqr(gmul(gmulsg(2, gaddgs(gmulsg(2, k), 1)), a)));
}

GEN
tt(GEN n, GEN a, GEN k)
{
  return gsub(gsqr(gaddgs(gadd(gsub(gsqr(n), gsqr(a)), gmulsg(4, gsqr(gaddgs(k, 1)))), 1)), gmulsg(4, gadd(gmul(gmulsg(4, gaddgs(gsqr(n), 1)), gsqr(gaddgs(k, 1))), gsqr(n))));
}

/* ff(n,a,k) = tt(n,a,k)^2 - (16*(k+1)*n*a)^2; */
/* t0(n,a) = (n^2-a^2+5)^2 - 4*(5*n^2+4); */

/* print("\n** ",ff(n,a,0)*ff(n,a,1) - (tt(n,a,0)^2 - (16*((0+1)*n*a))^2)); */

/*\\\\\\\\\\\\\\\\\\\\\\\\\\\ */
/* Calcola, con a = 8*i+4, il seguente prodotto. */
/* I primi due valori passati come argomento sono n^2 e a^2. */
/*  */
/*   (n-(8*i+7))(n-(8*i+5))(n-(8*i+3))(n-(8*i+1))(n+(8*i+1))(n+(8*i+3))(n+(8*i+5))(n+(8*i+7)) = */
/* = (n^-(8*i+7)^2)(n^2-(8*i+5)^2)(n^2-(8*i+3)^2)(n^2-(8*i+1)^2) */
/* = (n^-(a-3)^2)(n^2-(a-1)^2)(n^2-(a+1)^2)(n^2-(a+3)^2) */
/* */
/* = ((n^2-a^2+5)^2-4*(5*n^2+4))^2 - 256*n^2*a^2  <= Algoritmo di calcolo */

GEN
prod8(GEN nSqr, GEN aSqr, GEN nSqrAsqrShifted8, GEN subtrahend)
{
  GEN s, t, res;
  s = gaddgs(gsub(nSqr, aSqr), 5);
  t = gsub(gsqr(s), subtrahend);
  res = gsub(gsqr(t), nSqrAsqrShifted8);
  return res;
}

/*\\\\\\\\\\\\\\\\ */
/* Funzioni "naif" per confronto: calcolo del fattoriale doppio */
/* (per argomenti N dispari) con ricorsione binaria. */

/* Questa calcola i prodotti dei numeri dispari da start ad end, ricorsivamente. */
GEN
fattorialeDop(GEN start, GEN end)
{
  GEN tmp = gen_0, res;
  res = gcopy(start);
  if (!gequal(end, start))
  {
    if (gequalgs(gsub(end, start), 2))
      res = gmul(res, end);
    else
    {
      /* else */
      tmp = gshift(gadd(end, start), -1);
      if (gequal0(gmodgs(tmp, 2)))
        tmp = gaddgs(tmp, 1);
      /* print("\ttmp = ",tmp); */
      res = gmul(fattorialeDop(start, tmp), fattorialeDop(gaddgs(tmp, 2), end));
    }
  }
  return res;
}

/* ...e questa usa la precedente. */
GEN
fattorialeDoppio(GEN N)
{
  GEN tmp = gen_0, res = gen_1;
  if (gcmpgs(N, 5) < 0)
    res = gcopy(N);
  else
  {
    tmp = gshift(gaddgs(N, 1), -1);
    if (gequal0(gmodgs(tmp, 2)))
      tmp = gaddgs(tmp, 1);
    res = gmul(fattorialeDop(stoi(3), tmp), fattorialeDop(gaddgs(tmp, 2), N));
  }
  return res;
}

/*\\\\\\\\\\\\\\\\ */
/* Calcola l'n giusto per la formula. */
GEN
computeN(GEN N)	  /* vec */
{
  GEN quozRem;
  GEN p1;	  /* vec */
  quozRem = divrem(gaddgs(N, 1), stoi(16), -1);
  p1 = cgetg(3, t_VEC);
  gel(p1, 1) = gadd(gmulsg(8, gel(quozRem, 1)), gel(quozRem, 2));
  gel(p1, 2) = gcopy(gel(quozRem, 1));
  return p1;
}

/* (a+8)^2         = a^2 + 16a + 64 = a^2 +(a+4)<<4 */
/* (n^2(a+8)^2)<<8 = (n^2(a^2 +(a+4)<<4))<<8 = n^2a^2<<8 + n^2(a+4)<<12 */
/* (a+8)<<4 = (a<<4) + 1<<5; */

GEN
prodOdd(GEN N)
{
  GEN n = gen_0, nSqr = gen_0, aSqr = stoi(16), aFourth = stoi(128), subtrahend = gen_0, addend = gen_0, nSqrShifted12 = gen_0, nSqrAsqrShifted8 = gen_0, tmp = gen_0, index, res = gen_0, p1 = gen_0, nSqrShifted12TimesA = pol_x(fetch_user_var("nSqrShifted12TimesA")), p2;
  /* a = 8*i+4;   i+1 => a' = 8(i+1)+4 = (8i+4) + 8 = a+8 */
  index = gbitand(gaddgs(N, 1), stoi(15));
  if (gequal0(index))
    p1 = gen_1;
  else
    p1 = gcopy(gel(PRODOTTI_FATTORI_DISPARI, itos(gshift(index, -1))));
  res = p1;
  /* Preparazione costanti di base, in modo che nel ciclo principale ci siano solo somme. */
  tmp = computeN(N);
  n = gcopy(gel(tmp, 1));
  nSqr = gsqr(n);
  nSqrShifted12 = gshift(nSqr, 12);
  nSqrAsqrShifted8 = nSqrShifted12;
  subtrahend = gshift(gaddgs(gmulsg(5, nSqr), 4), 2);
  addend = gshift(nSqrShifted12, 3);
  nSqrShifted12TimesA = addend;
  /* Ciclo principale di prodotti. */
  p2 = gsubgs(gel(tmp, 2), 1);
  {
    GEN i;
    for (i = gen_0; gcmp(i, p2) <= 0; i = gaddgs(i, 1))
    {
      res = gmul(res, prod8(nSqr, aSqr, nSqrAsqrShifted8, subtrahend));
      aSqr = gadd(aSqr, aFourth);
      aFourth = gaddgs(aFourth, 128);
      nSqrAsqrShifted8 = gadd(nSqrAsqrShifted8, nSqrShifted12TimesA);
      nSqrShifted12TimesA = gadd(nSqrShifted12TimesA, addend);
    }
  }
  return res;
}

/*\\\\\\\\\\\\\\\\\\\\\\\\ */
/* Calcola il prodotti degli elementi del vettore contenuti */
/* nelle posizioni (indici) da start a end. */
GEN
computeProductsFromTo(GEN v, GEN start, GEN end)
{
  GEN res, LIM = stoi(4), half = gen_0;
  res = gcopy(gel(v, gtos(start)));
  if (gcmp(gsub(end, start), LIM) < 0)
  {
    GEN p1;
    p1 = gaddgs(start, 1);
    {
      GEN i;
      for (i = p1; gcmp(i, end) <= 0; i = gaddgs(i, 1))
        res = gmul(res, gel(v, gtos(i)));
    }
  }
  else
  {
    /* else */
    half = gshift(gadd(start, end), -1);
    res = gmul(computeProductsFromTo(v, start, half), computeProductsFromTo(v, gaddgs(half, 1), end));
  }
  return res;
}

/*\\\\\\\\\\\\\\\\\\\\\\\\ */
/* Calcola il prodotti degli elementi del vettore */
GEN
computeProducts(GEN v)
{
  GEN lim, res;
  lim = stoi(glength(v)>>1);
  res = gcopy(gel(v, 1));
  if (glength(v) > 1)
    res = gmul(computeProductsFromTo(v, gen_1, lim), computeProductsFromTo(v, gaddgs(lim, 1), stoi(glength(v))));
  return res;
}

/*\\\\\\\\\\\\\\\\\\\\\\\\ Più efficiente. */
/* Fattoriale doppio. N!! = N(N-2)(N-4)···5*3*1 */
GEN
prodOdd2(GEN N)
{
  GEN n = gen_0, nSqr = gen_0, a, aSqr = stoi(16), aFourth = stoi(128), subtrahend = gen_0, addend = gen_0, nSqrShifted12 = gen_0, nSqrAsqrShifted8 = gen_0, tmp = gen_0, incrA, resVec = gen_0, index, res = gen_0, p1 = gen_0, nSqrShifted12TimesA = pol_x(fetch_user_var("nSqrShifted12TimesA")), p2;
  GEN p3;	  /* vec */
  GEN p4;
  /* a = 8*i+4;   i+1 => a' = 8(i+1)+4 = (8i+4) + 8 = a+8 */
  a = shifti(stoi(-3), 31);
  incrA = shifti(gen_1, 31);
  index = gbitand(gaddgs(N, 1), stoi(15));
  if (gequal0(index))
    p1 = gen_1;
  else
    p1 = gcopy(gel(PRODOTTI_FATTORI_DISPARI, itos(gshift(index, -1))));
  res = p1;
  tmp = computeN(N);
  n = gcopy(gel(tmp, 1));
  nSqr = gsqr(n);
  subtrahend = gshift(gaddgs(gmulsg(5, nSqr), 4), 2);
  nSqrShifted12 = gshift(nSqr, 12);
  nSqrAsqrShifted8 = nSqrShifted12;
  addend = gshift(nSqrShifted12, 3);
  nSqrShifted12TimesA = addend;
  p2 = gcopy(gel(tmp, 2));
  {
    long l5;
    p3 = cgetg(gtos(p2)+1, t_VEC);
    for (l5 = 1; gcmpsg(l5, p2) <= 0; ++l5)
      gel(p3, l5) = gen_0;
  }
  resVec = p3;
  /* Vettore per organizzare i risultati. */
  
  /* Ciclo di prodotti: memorizzazione dei risultati parziali in un vettore. */
  p4 = gcopy(gel(tmp, 2));
  {
    GEN i;
    for (i = gen_1; gcmp(i, p4) <= 0; i = gaddgs(i, 1))
    {
      if (gcmpgs(i, 7) > 0)
        gel(resVec, gtos(i)) = gadd(gmulsg(7, gadd(gadd(gsub(gsub(gel(resVec, gtos(gsubgs(i, 1))), gel(resVec, gtos(gsubgs(i, 6)))), gmulsg(3, gsub(gel(resVec, gtos(gsubgs(i, 2))), gel(resVec, gtos(gsubgs(i, 5)))))), gmulsg(5, gsub(gel(resVec, gtos(gsubgs(i, 3))), gel(resVec, gtos(gsubgs(i, 4)))))), gmulsg(45, a))), gel(resVec, gtos(gsubgs(i, 7))));
      else
      {
        gel(/*else */
        resVec, gtos(i)) = prod8(nSqr, aSqr, nSqrAsqrShifted8, subtrahend);
        aSqr = gadd(aSqr, aFourth);
        aFourth = gaddgs(aFourth, 128);
        nSqrAsqrShifted8 = gadd(nSqrAsqrShifted8, nSqrShifted12TimesA);
        nSqrShifted12TimesA = gadd(nSqrShifted12TimesA, addend);
      }
      a = gadd(a, incrA);
    }
  }
  if (!gequal0(gel(tmp, 2)))
    /* Calcolo dei prodotti nel vettore. */
    res = gmul(res, computeProducts(resVec));
  return res;
}

/*\\\\\\\\\\\\\\\\\\\\\\\\ Ancora (un po') più efficiente. */
/* Fattoriale doppio. N!! = N(N-2)(N-4)···5*3*1 */
/* addend2(aa) = -7*aa^2 - 336*aa - 4251 */
/* addend2(aa+8)-addend2(aa) = -112*aa - 3136    ;     addend(4) = -5707 */

GEN
prodOdd3(GEN N)
{
  GEN n = gen_0, nSqr = gen_0, aSqr = stoi(16), aFourth = stoi(128), subtrahend = gen_0, addend = gen_0, nSqrAsqrShifted8 = gen_0, tmp = gen_0, addend2 = gen_0, nSqrTimes45 = gen_0, resVec = gen_0, index, incr, step, res = gen_0, p1 = gen_0, nSqrShifted12TimesA = pol_x(fetch_user_var("nSqrShifted12TimesA")), p2;
  GEN p3;	  /* vec */
  GEN p4;
  /* a = 8*i+4;   i+1 => a' = 8(i+1)+4 = (8i+4) + 8 = a+8 */
  index = gbitand(gaddgs(N, 1), stoi(15));
  incr = shifti(stoi(945), 31);
  step = shifti(stoi(315), 31);
  if (gequal0(index))
    p1 = gen_1;
  else
    p1 = gcopy(gel(PRODOTTI_FATTORI_DISPARI, itos(gshift(index, -1))));
  res = p1;
  tmp = computeN(N);
  n = gcopy(gel(tmp, 1));
  nSqr = gsqr(n);
  subtrahend = gshift(gaddgs(gmulsg(5, nSqr), 4), 2);
  nSqrTimes45 = gmulgs(nSqr, 45);
  addend2 = gshift(gsubsg(256815, nSqrTimes45), 24);
  nSqrAsqrShifted8 = gshift(nSqr, 12);
  addend = gshift(nSqr, 15);
  nSqrShifted12TimesA = addend;
  p2 = gcopy(gel(tmp, 2));
  {
    long l5;
    p3 = cgetg(gtos(p2)+1, t_VEC);
    for (l5 = 1; gcmpsg(l5, p2) <= 0; ++l5)
      gel(p3, l5) = gen_0;
  }
  resVec = p3;
  /* Vettore per organizzare i risultati. */
  
  /* Ciclo di prodotti: memorizzazione dei risultati parziali in un vettore. */
  p4 = gcopy(gel(tmp, 2));
  {
    GEN i;
    for (i = gen_1; gcmp(i, p4) <= 0; i = gaddgs(i, 1))
    {
      if (gcmpgs(i, 6) > 0)
      {
        gel(resVec, gtos(i)) = gadd(gsub(gadd(gsub(gmulsg(6, gadd(gel(resVec, gtos(gsubgs(i, 1))), gel(resVec, gtos(gsubgs(i, 5))))), gmulsg(15, gadd(gel(resVec, gtos(gsubgs(i, 2))), gel(resVec, gtos(gsubgs(i, 4)))))), gmulsg(20, gel(resVec, gtos(gsubgs(i, 3))))), gel(resVec, gtos(gsubgs(i, 6)))), addend2);
        addend2 = gadd(addend2, incr = gadd(incr, step));
      }
      else
      {
        gel(/*else */
        resVec, gtos(i)) = prod8(nSqr, aSqr, nSqrAsqrShifted8, subtrahend);
        aSqr = gadd(aSqr, aFourth);
        aFourth = gaddgs(aFourth, 128);
        nSqrAsqrShifted8 = gadd(nSqrAsqrShifted8, nSqrShifted12TimesA);
        nSqrShifted12TimesA = gadd(nSqrShifted12TimesA, addend);
      }
    }
  }
  if (gcmpgs(gel(tmp, 2), 6) > 0)
  {
    long l6;
    l6 = gtos(gel(tmp, 2));
    gel(/* Prodotto parziale nel fattore più piccolo: invece */
    /* di una moltiplicazione finale lunga ne si fa una */
    /* all'inizio, molto più piccola. */
    resVec, l6) = gmul(gel(resVec, l6), res);
    /* Calcolo dei prodotti nel vettore. */
    res = computeProducts(resVec);
  }
  return res;
}

/*\\\\\\\\\\\\ */

GEN
tempi(GEN N, GEN volte, long prec)
{
  GEN fatt = gen_1, r2 = gen_0, r3 = gen_0, t;
  GEN p1;	  /* vec */
  if (!volte)
    volte = gen_1;
  {
    long l2;
    p1 = cgetg(4, t_VEC);
    for (l2 = 1; l2 <= 3; ++l2)
      gel(p1, l2) = gen_0;
  }
  t = p1;
  gel(t, 1) = stoi(gettime());
  {
    GEN j;
    for (j = gen_1; gcmp(j, volte) <= 0; j = gaddgs(j, 1))
      fatt = fattorialeDoppio(N);
  }
  gel(t, 1) = gmaxgs(gsubsg(gettime(), gel(t, 1)), 1);
  gel(t, 2) = stoi(gettime());
  {
    GEN j;
    for (j = gen_1; gcmp(j, volte) <= 0; j = gaddgs(j, 1))
      r2 = prodOdd2(N);
  }
  gel(t, 2) = gmaxgs(gsubsg(gettime(), gel(t, 2)), 1);
  gel(t, 3) = stoi(gettime());
  {
    GEN j;
    for (j = gen_1; gcmp(j, volte) <= 0; j = gaddgs(j, 1))
      r3 = prodOdd3(N);
  }
  gel(t, 3) = gsubsg(gettime(), gel(t, 3));
  pari_printf("%Ps\t%ld\t%Ps\t%Ps\t%Ps\t%Ps %Ps %%\t%Ps %%\n", N, gequal(fatt, r2) && gequal(fatt, r3), gel(t, 1), gel(t, 2), gel(t, 3), gdiv(gmul(stor(100, prec), gel(t, 2)), gel(t, 1)), gdiv(gmul(stor(100, prec), gel(t, 3)), gel(t, 1)), gdiv(gmul(stor(100, prec), gel(t, 3)), gel(t, 2)));
  return r2;
}

/*\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */
/* tempi(N, numero di ripetizioni del test) */
long
df8test(long prec)
{
  GEN end = pol_x(fetch_user_var("end"));
  pari_printf("Confronto tempi per fattoriale doppio.\nt1: ricorsione binaria\tt2: Zanoni (ricorsione 7)\tt3: Zanoni (ricorsione 6)\n");
  pari_printf("-----------------------------------------------------\n");
  pari_printf("N\tCheck\tt1\tt2\tt3\tt2/t1\tt3/t1\tt3/t2\n");
  end = gadd(START, NUMERO_ESEMPI);
  {
    GEN h;
    for (h = START; gcmp(h, end) <= 0; h = gaddgs(h, 2))
      tempi(h, gen_1, prec);
  }
  /* print(h,"!! = ",FZ) */ return 0;
}

///////////////////////
// main function.

int main()
{
 pari_init(1e9,2);

 df8test(0);
 pari_close();

 return 0;
}