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