/*
    Copyright (C) 2022-2025 Alessandro Languasco
*/


\\ Global variables:
global(defaultprecision);

/*************  GENERAL CASE  q = +-1 mod 8, prime *********/

{Aqdivisors(q,defaultprecision)=local(minutes, millisec, seconds, ghcommon,
elaptimefinalcomp, qminusone, a, g, G, h, H, R, A, E, rho, rhoalphan, rhobetan,
alphan, nq, betan, prime1, prime2, qclass8, kroneckerchar, w, b1, b2,
rhodivisorsfile, logtwo, logten, accuracy, zetafactor, trunckroneckerchar,
startAproduct, startalphaproduct, startbetaproduct, elapbetaproducttime, elapalphaproducttime,
elapAproducttime, start, elappreptime, nqprime
);

\\ DO NOT DEFINE the variables x,p ; they are formal variables needed in the polynomials 

\\ minutes,millisec,seconds: used just to compute the elapsed computation time; local variables

\\ defaultprecision: used to fix the precision used in the computations; global variable

print("************ A. LANGUASCO *************");
print("********* COMPUTATION of Aq(1) - d(n) mod q **********"); 
print("********* q = -1,1 mod 8   **********"); 

\\ precision setting 

default(realprecision,defaultprecision);       
default(format,f);   

if (isprime(q) == 0,  print("q has to prime "); return() );
qminusone = q-1; 
qclass8 = q%8; 
if (qclass8 <> 1 && qclass8 <> 7, print("q has to be = -1,1 mod 8"); return() );

rhodivisorsfile = fileopen("rhodivfile-prodeulerrat.csv", "a");
\\filewrite(rhodivisorsfile,"q;nq;nqprime;Aq(1);rho(q);rhonorm;rhoalpha;rhoalphanorm;rhobeta;rhobetanorm;accuracy;millesecs");  


logten = log(10);
logtwo = log(2);
accuracy = defaultprecision*logten;

\\ zeta(x) = 1 + ERR(x), with |ERR(x)| < 2^(1-x) for x>=3
zetafactor = 1.;
\\ means |ERR(q-1)| < 2^(2-q) < 10^(-defaultprecision)
if (q <= 2 + accuracy / logtwo, zetafactor = zeta(qminusone));
 
print("Starting computation of Aq - divisors"); 
print("---------");
  
\\ initializations 

start=gettime();
\\ quadratic Kronecker character (chi_1) values
kroneckerchar = vector(qminusone, j, kronecker(j,q) );

\\nq determination (2|q) = 1 since q =+-1 mod 8 
forprime(j=3, qminusone, if (kroneckerchar[j] == -1 , nq = j; break));

\\nqprime determination 
forprime(j=nq+1, qminusone, if (kroneckerchar[j] == -1 , nqprime = j; break));
 
\\ primes for Lemma 2 functions
prime1 = nextprime(nq+1); 
prime2 = nextprime(prime1+1);	 

/***** generating the polynomials *****/
\\ to truncate the polynomial for g and h, see their definition in the paper
b1 = prime1-1;
b2 = prime2-1;

\\  generating coefficients for g and h
trunckroneckerchar = vector(b2, j, kroneckerchar[j]);
trunckroneckerchar[prime1] = 0; \\ to discard x^{p1-1}; see the definition in the paper


elappreptime=gettime()-start;

startAproduct= gettime();
\\ Aq Euler product
a = (1-x)*Polrev(kroneckerchar); \\ generates the Fekete polynomial for Aq(1)
A = subst(a, x, 1/p); \\ formal substitution of x->1/p 
R = prodeulerrat(A);  \\ evaluates the Euler product \\ BY FAR, THIS IS THE MOST TIME-CONSUMING PART
R *= zetafactor;  \\ R = zetafactor * R; see A(q) definition in the paper
elapAproducttime=gettime()-startAproduct;


startalphaproduct= gettime();
\\ generates the common part of the polynomials (1-x)*g(x) and (1-x)*h(x)
ghcommon = (1-x)*Polrev(trunckroneckerchar);

\\ generate the error polynomial (1-x)*E(x)
E = Pol((1-x)*x^(b1)+x^(b2));

\\ alpha Euler product
g = ghcommon-E; \\ getting the polynomial g for the product alpha
\\print(g);
G = subst(g, x, 1/p); \\ formal substitution of x->1/p 
alphan = prodeulerrat(G); \\ evaluates the Euler product 
elapalphaproducttime=gettime()-startalphaproduct;

startbetaproduct= gettime();
\\ beta Euler product
h = ghcommon + E; \\ getting the polynomial h for the product beta
\\print(h);
H = subst(h, x, 1/p); \\ formal substitution of x->1/p 
betan = prodeulerrat(H); \\ evaluates the Euler product 
elapbetaproducttime=gettime()-startbetaproduct;

elaptimefinalcomp=elappreptime+elapAproducttime+elapalphaproducttime+elapbetaproducttime;

print("q = ", q);
print("nq(", q, ") = ", nq);
print("nq'(", q, ") = ", nqprime);
w = 2^(-nq);

print("---");
print("A(", q, ",1) = ", R );     
rho = (1 + R)/(1-R); 
print("rho = ", rho);     
print("rho_tilde = ", rho * w);     
  
print("---");
print("betan(", nq, "," , q, ") = ", betan );     
rhobetan = (1 + betan)/(1-betan);
print("rhonbeta = ", rhobetan );
print("rhonbetanorm = ", rhobetan * w );     

print("---");
print("alphan(", nq, "," , q, ") = ", alphan );     
rhoalphan = (1 + alphan)/(1-alphan); 
print("rhonalpha = ", rhoalphan );    
print("rhonalphanorm = ", rhoalphan * w );    

print("---");
 
print("q & nq & nqprime & A(q,1) & rho=(1+A(q,1))/(1-A(q,1)) & rho_tilde");
print("$", q, "$ & $", nq, "$ & $", nqprime, "$ & $", R, " $ & $", rho, " $ & $", rho * w, "$ \\\\ \\hline");

print("---"); 
print("nq & q & rhonalpha & rhonbeta ");
print("$", nq, "$ & $", q, "$ & $", rhoalphan, " $ & $", rhobetan, " $ & $", rho,  "$ \\\\ \\hline");

print("---------");

filewrite1(rhodivisorsfile,q);
filewrite1(rhodivisorsfile,";"); 
filewrite1(rhodivisorsfile,nq);
filewrite1(rhodivisorsfile,";"); 
filewrite1(rhodivisorsfile,nqprime);
filewrite1(rhodivisorsfile,";"); 
filewrite1(rhodivisorsfile,R);
filewrite1(rhodivisorsfile,";"); 
filewrite1(rhodivisorsfile,rho);
filewrite1(rhodivisorsfile,";");  
filewrite1(rhodivisorsfile,rho*w);   

filewrite1(rhodivisorsfile,";"); 
filewrite1(rhodivisorsfile,rhoalphan);
filewrite1(rhodivisorsfile,";"); 
filewrite1(rhodivisorsfile,rhoalphan * w);
  
filewrite1(rhodivisorsfile,";"); 
filewrite1(rhodivisorsfile,rhobetan);
filewrite1(rhodivisorsfile,";"); 
filewrite1(rhodivisorsfile,rhobetan * w );
filewrite1(rhodivisorsfile,";"); 
filewrite1(rhodivisorsfile, defaultprecision );
filewrite1(rhodivisorsfile,";"); 
filewrite(rhodivisorsfile, elaptimefinalcomp );
fileflush(rhodivisorsfile); 						
fileclose(rhodivisorsfile);  

seconds=floor(elaptimefinalcomp/1000)%60;
minutes=floor(elaptimefinalcomp/60000);
millisec=elaptimefinalcomp- minutes*60000 - seconds*1000;
 
print("Total computation time (output time included): ", minutes, " min, ", seconds, " sec, ", millisec, " millisec");

seconds=floor(elappreptime/1000)%60;
minutes=floor(elappreptime/60000);
millisec=elappreptime- minutes*60000 - seconds*1000;
 
print("Preparation time (output time included): ", minutes, " min, ", seconds, " sec, ", millisec, " millisec");

seconds=floor(elapAproducttime/1000)%60;
minutes=floor(elapAproducttime/60000);
millisec=elapAproducttime- minutes*60000 - seconds*1000;
 
print("A(q) product computation time (output time included): ", minutes, " min, ", seconds, " sec, ", millisec, " millisec");


seconds=floor(elapalphaproducttime/1000)%60;
minutes=floor(elapalphaproducttime/60000);
millisec=elapalphaproducttime- minutes*60000 - seconds*1000;
 
print("ALPHA product computation time (output time included): ", minutes, " min, ", seconds, " sec, ", millisec, " millisec");


seconds=floor(elapbetaproducttime/1000)%60;
minutes=floor(elapbetaproducttime/60000);
millisec=elapbetaproducttime- minutes*60000 - seconds*1000;
 
print("BETA product computation time (output time included): ", minutes, " min, ", seconds, " sec, ", millisec, " millisec");



print("******  END PROGRAM   ********"); 
}




 
/**** RESULTS ****

New run july 27th 2025  
 
---------
macbook ARM
---------


languasc@rozycki ~ % gp
                  GP/PARI CALCULATOR Version 2.17.2 (released)
         arm64 running darwin (aarch64/GMP-6.3.0 kernel) 64-bit version
     compiled: Mar  7 2025, Apple clang version 16.0.0 (clang-1600.0.26.6)
                    threading engine: pthread, nbthreads = 8
                 (readline v8.2 enabled, extended help enabled)

                     Copyright (C) 2000-2024 The PARI Group

PARI/GP is free software, covered by the GNU General Public License, and comes 
WITHOUT ANY WARRANTY WHATSOEVER.

Type ? for help, \q to quit.
Type ?18 for how to get moral (and possibly technical) support.

parisize = 8000000, primelimit = 1048576, factorlimit = 1048576

------
q <= 100
------
************ A. LANGUASCO *************
********* COMPUTATION of Aq(1) - d(n) mod q **********
********* q = -1,1 mod 8   **********
Starting computation of Aq - divisors
---------
q = 7
nq(7) = 3
nq'(7) = 5
---
A(7,1) = 0.458175374302921827638937646284263464587768628852480873516
rho = 2.69123126773155290630742717482368921730707567404537665645
rho_tilde = 0.336403908466444113288428396852961152163384459255672082056
---
betan(3,7) = 0.518129296788080596043569544939347180024716847434338906307
rhonbeta = 3.15049096504302405901097604117785432217628979120422980775
rhonbetanorm = 0.393811370630378007376372005147231790272036223900528725969
---
alphan(3,7) = 0.438618380624200524969748303109148800399703288765876559432
rhonalpha = 2.56263890902555925733978510811007661265851707999689944404
rhonalphanorm = 0.320329863628194907167473138513759576582314634999612430505
---
q & nq & A(q,1) & rho=(1+A(q,1))/(1-A(q,1)) & rho_tilde
$7$ & $3$ & $0.458175374302921827638937646284263464587768628852480873516 $ & $2.69123126773155290630742717482368921730707567404537665645 $ & $0.336403908466444113288428396852961152163384459255672082056$ \\ \hline
---
nq & q & rhonalpha & rhonbeta 
$3$ & $7$ & $2.56263890902555925733978510811007661265851707999689944404 $ & $3.15049096504302405901097604117785432217628979120422980775 $ & $2.69123126773155290630742717482368921730707567404537665645$ \\ \hline
---------
Total computation time (output time included): -1 min, 59 sec, 989 millisec
Preparation time (output time included): -1 min, 59 sec, 983 millisec
A(q) product computation time (output time included): 0 min, 0 sec, 2 millisec
ALPHA product computation time (output time included): 0 min, 0 sec, 2 millisec
BETA product computation time (output time included): 0 min, 0 sec, 2 millisec
******  END PROGRAM   ********
************ A. LANGUASCO *************
********* COMPUTATION of Aq(1) - d(n) mod q **********
********* q = -1,1 mod 8   **********
Starting computation of Aq - divisors
---------
q = 17
nq(17) = 3
nq'(17) = 5
---
A(17,1) = 0.447549513648168557638511450029638552429559551867483867732
rho = 2.62023393844256274534572641894361825899149547051879521437
rho_tilde = 0.327529242305320343168215802367952282373936933814849401796
---
betan(3,17) = 0.518129296788080596043569544939347180024716847434338906307
rhonbeta = 3.15049096504302405901097604117785432217628979120422980775
rhonbetanorm = 0.393811370630378007376372005147231790272036223900528725969
---
alphan(3,17) = 0.438618380624200524969748303109148800399703288765876559432
rhonalpha = 2.56263890902555925733978510811007661265851707999689944404
rhonalphanorm = 0.320329863628194907167473138513759576582314634999612430505
---
q & nq & A(q,1) & rho=(1+A(q,1))/(1-A(q,1)) & rho_tilde
$17$ & $3$ & $0.447549513648168557638511450029638552429559551867483867732 $ & $2.62023393844256274534572641894361825899149547051879521437 $ & $0.327529242305320343168215802367952282373936933814849401796$ \\ \hline
---
nq & q & rhonalpha & rhonbeta 
$3$ & $17$ & $2.56263890902555925733978510811007661265851707999689944404 $ & $3.15049096504302405901097604117785432217628979120422980775 $ & $2.62023393844256274534572641894361825899149547051879521437$ \\ \hline
---------
Total computation time (output time included): 0 min, 0 sec, 3 millisec
Preparation time (output time included): 0 min, 0 sec, 0 millisec
A(q) product computation time (output time included): 0 min, 0 sec, 1 millisec
ALPHA product computation time (output time included): 0 min, 0 sec, 1 millisec
BETA product computation time (output time included): 0 min, 0 sec, 1 millisec
******  END PROGRAM   ********
************ A. LANGUASCO *************
********* COMPUTATION of Aq(1) - d(n) mod q **********
********* q = -1,1 mod 8   **********
Starting computation of Aq - divisors
---------
q = 23
nq(23) = 5
nq'(23) = 7
---
A(23,1) = 0.898576458864156997644651084962786708511332965074728471471
rho = 18.7192878261001825739181174131915574081394448346863422941
rho_tilde = 0.584977744565630705434941169162236169004357651083948196691
---
betan(5,23) = 0.916850779228909835416677911291295302299925792490423457856
rhonbeta = 23.0531418268609003184857276806683038586857850885743587848
rhonbetanorm = 0.720410682089403134952678990020884495583930784017948712025
---
alphan(5,23) = 0.897813044296315189281089379777004468377906648563104632058
rhonalpha = 18.5719696924866996139207933541007149680235307403926994739
rhonalphanorm = 0.580374052890209362935024792315647342750735335637271858558
---
q & nq & A(q,1) & rho=(1+A(q,1))/(1-A(q,1)) & rho_tilde
$23$ & $5$ & $0.898576458864156997644651084962786708511332965074728471471 $ & $18.7192878261001825739181174131915574081394448346863422941 $ & $0.584977744565630705434941169162236169004357651083948196691$ \\ \hline
---
nq & q & rhonalpha & rhonbeta 
$5$ & $23$ & $18.5719696924866996139207933541007149680235307403926994739 $ & $23.0531418268609003184857276806683038586857850885743587848 $ & $18.7192878261001825739181174131915574081394448346863422941$ \\ \hline
---------
Total computation time (output time included): 0 min, 0 sec, 3 millisec
Preparation time (output time included): 0 min, 0 sec, 0 millisec
A(q) product computation time (output time included): 0 min, 0 sec, 1 millisec
ALPHA product computation time (output time included): 0 min, 0 sec, 1 millisec
BETA product computation time (output time included): 0 min, 0 sec, 1 millisec
******  END PROGRAM   ********
************ A. LANGUASCO *************
********* COMPUTATION of Aq(1) - d(n) mod q **********
********* q = -1,1 mod 8   **********
Starting computation of Aq - divisors
---------
q = 31
nq(31) = 3
nq'(31) = 11
---
A(31,1) = 0.516820785012757788705384212520040784046479120252761760392
rho = 3.13925089897090810987728952626151370494442874685073440490
rho_tilde = 0.392406362371363513734661190782689213118053593356341800612
---
betan(3,31) = 0.518129296788080596043569544939347180024716847434338906307
rhonbeta = 3.15049096504302405901097604117785432217628979120422980775
rhonbetanorm = 0.393811370630378007376372005147231790272036223900528725969
---
alphan(3,31) = 0.438618380624200524969748303109148800399703288765876559432
rhonalpha = 2.56263890902555925733978510811007661265851707999689944404
rhonalphanorm = 0.320329863628194907167473138513759576582314634999612430505
---
q & nq & A(q,1) & rho=(1+A(q,1))/(1-A(q,1)) & rho_tilde
$31$ & $3$ & $0.516820785012757788705384212520040784046479120252761760392 $ & $3.13925089897090810987728952626151370494442874685073440490 $ & $0.392406362371363513734661190782689213118053593356341800612$ \\ \hline
---
nq & q & rhonalpha & rhonbeta 
$3$ & $31$ & $2.56263890902555925733978510811007661265851707999689944404 $ & $3.15049096504302405901097604117785432217628979120422980775 $ & $3.13925089897090810987728952626151370494442874685073440490$ \\ \hline
---------
Total computation time (output time included): 0 min, 0 sec, 5 millisec
Preparation time (output time included): 0 min, 0 sec, 0 millisec
A(q) product computation time (output time included): 0 min, 0 sec, 1 millisec
ALPHA product computation time (output time included): 0 min, 0 sec, 2 millisec
BETA product computation time (output time included): 0 min, 0 sec, 2 millisec
******  END PROGRAM   ********
************ A. LANGUASCO *************
********* COMPUTATION of Aq(1) - d(n) mod q **********
********* q = -1,1 mod 8   **********
Starting computation of Aq - divisors
---------
q = 41
nq(41) = 3
nq'(41) = 7
---
A(41,1) = 0.504311182221407966673019915352424830736647370961233677443
rho = 3.03478942487125969929494971797080264176721891080174709090
rho_tilde = 0.379348678108907462411868714746350330220902363850218386363
---
betan(3,41) = 0.518129296788080596043569544939347180024716847434338906307
rhonbeta = 3.15049096504302405901097604117785432217628979120422980775
rhonbetanorm = 0.393811370630378007376372005147231790272036223900528725969
---
alphan(3,41) = 0.438618380624200524969748303109148800399703288765876559432
rhonalpha = 2.56263890902555925733978510811007661265851707999689944404
rhonalphanorm = 0.320329863628194907167473138513759576582314634999612430505
---
q & nq & A(q,1) & rho=(1+A(q,1))/(1-A(q,1)) & rho_tilde
$41$ & $3$ & $0.504311182221407966673019915352424830736647370961233677443 $ & $3.03478942487125969929494971797080264176721891080174709090 $ & $0.379348678108907462411868714746350330220902363850218386363$ \\ \hline
---
nq & q & rhonalpha & rhonbeta 
$3$ & $41$ & $2.56263890902555925733978510811007661265851707999689944404 $ & $3.15049096504302405901097604117785432217628979120422980775 $ & $3.03478942487125969929494971797080264176721891080174709090$ \\ \hline
---------
Total computation time (output time included): 0 min, 0 sec, 7 millisec
Preparation time (output time included): 0 min, 0 sec, 0 millisec
A(q) product computation time (output time included): 0 min, 0 sec, 3 millisec
ALPHA product computation time (output time included): 0 min, 0 sec, 2 millisec
BETA product computation time (output time included): 0 min, 0 sec, 2 millisec
******  END PROGRAM   ********
************ A. LANGUASCO *************
********* COMPUTATION of Aq(1) - d(n) mod q **********
********* q = -1,1 mod 8   **********
Starting computation of Aq - divisors
---------
q = 47
nq(47) = 5
nq'(47) = 11
---
A(47,1) = 0.915564527190921027656098571139780059874780359316593570499
rho = 22.6867270764539244689979834380851325291210940161697311036
rho_tilde = 0.708960221139185139656186982440160391535034188005304096987
---
betan(5,47) = 0.916850779228909835416677911291295302299925792490423457856
rhonbeta = 23.0531418268609003184857276806683038586857850885743587848
rhonbetanorm = 0.720410682089403134952678990020884495583930784017948712025
---
alphan(5,47) = 0.897813044296315189281089379777004468377906648563104632058
rhonalpha = 18.5719696924866996139207933541007149680235307403926994739
rhonalphanorm = 0.580374052890209362935024792315647342750735335637271858558
---
q & nq & A(q,1) & rho=(1+A(q,1))/(1-A(q,1)) & rho_tilde
$47$ & $5$ & $0.915564527190921027656098571139780059874780359316593570499 $ & $22.6867270764539244689979834380851325291210940161697311036 $ & $0.708960221139185139656186982440160391535034188005304096987$ \\ \hline
---
nq & q & rhonalpha & rhonbeta 
$5$ & $47$ & $18.5719696924866996139207933541007149680235307403926994739 $ & $23.0531418268609003184857276806683038586857850885743587848 $ & $22.6867270764539244689979834380851325291210940161697311036$ \\ \hline
---------
Total computation time (output time included): 0 min, 0 sec, 7 millisec
Preparation time (output time included): 0 min, 0 sec, 0 millisec
A(q) product computation time (output time included): 0 min, 0 sec, 3 millisec
ALPHA product computation time (output time included): 0 min, 0 sec, 2 millisec
BETA product computation time (output time included): 0 min, 0 sec, 2 millisec
******  END PROGRAM   ********
************ A. LANGUASCO *************
********* COMPUTATION of Aq(1) - d(n) mod q **********
********* q = -1,1 mod 8   **********
Starting computation of Aq - divisors
---------
q = 71
nq(71) = 7
nq'(71) = 11
---
A(71,1) = 0.981075357290720578572892039380146350418397461528754442659
rho = 104.682312248850504704730966602883462421955082486280965609
rho_tilde = 0.817830564444144568005710676585027050171524081924070043821
---
betan(7,71) = 0.982458136849913452663005468256860875609563613667685992346
rhonbeta = 113.012974727267353049355544740857282609681507300936023080
rhonbetanorm = 0.882913865056776195698090193287947520388136775788562680309
---
alphan(7,71) = 0.980970080424897228935791455593548999293768188213773382832
rhonalpha = 104.097659089250197303031479434237189758514457201861911005
rhonalphanorm = 0.813262961634767166429933433079978044988394196889546179727
---
q & nq & A(q,1) & rho=(1+A(q,1))/(1-A(q,1)) & rho_tilde
$71$ & $7$ & $0.981075357290720578572892039380146350418397461528754442659 $ & $104.682312248850504704730966602883462421955082486280965609 $ & $0.817830564444144568005710676585027050171524081924070043821$ \\ \hline
---
nq & q & rhonalpha & rhonbeta 
$7$ & $71$ & $104.097659089250197303031479434237189758514457201861911005 $ & $113.012974727267353049355544740857282609681507300936023080 $ & $104.682312248850504704730966602883462421955082486280965609$ \\ \hline
---------
Total computation time (output time included): 0 min, 0 sec, 5 millisec
Preparation time (output time included): 0 min, 0 sec, 0 millisec
A(q) product computation time (output time included): 0 min, 0 sec, 3 millisec
ALPHA product computation time (output time included): 0 min, 0 sec, 1 millisec
BETA product computation time (output time included): 0 min, 0 sec, 1 millisec
******  END PROGRAM   ********
************ A. LANGUASCO *************
********* COMPUTATION of Aq(1) - d(n) mod q **********
********* q = -1,1 mod 8   **********
Starting computation of Aq - divisors
---------
q = 73
nq(73) = 5
nq'(73) = 7
---
A(73,1) = 0.898339084168655588699878545658343212379173567120691795780
rho = 18.6732439762592985533629405864402714364011924734863103533
rho_tilde = 0.583538874258103079792591893326258482387537264796447198541
---
betan(5,73) = 0.916850779228909835416677911291295302299925792490423457856
rhonbeta = 23.0531418268609003184857276806683038586857850885743587848
rhonbetanorm = 0.720410682089403134952678990020884495583930784017948712025
---
alphan(5,73) = 0.897813044296315189281089379777004468377906648563104632058
rhonalpha = 18.5719696924866996139207933541007149680235307403926994739
rhonalphanorm = 0.580374052890209362935024792315647342750735335637271858558
---
q & nq & A(q,1) & rho=(1+A(q,1))/(1-A(q,1)) & rho_tilde
$73$ & $5$ & $0.898339084168655588699878545658343212379173567120691795780 $ & $18.6732439762592985533629405864402714364011924734863103533 $ & $0.583538874258103079792591893326258482387537264796447198541$ \\ \hline
---
nq & q & rhonalpha & rhonbeta 
$5$ & $73$ & $18.5719696924866996139207933541007149680235307403926994739 $ & $23.0531418268609003184857276806683038586857850885743587848 $ & $18.6732439762592985533629405864402714364011924734863103533$ \\ \hline
---------
Total computation time (output time included): 0 min, 0 sec, 6 millisec
Preparation time (output time included): 0 min, 0 sec, 0 millisec
A(q) product computation time (output time included): 0 min, 0 sec, 3 millisec
ALPHA product computation time (output time included): 0 min, 0 sec, 2 millisec
BETA product computation time (output time included): 0 min, 0 sec, 1 millisec
******  END PROGRAM   ********
************ A. LANGUASCO *************
********* COMPUTATION of Aq(1) - d(n) mod q **********
********* q = -1,1 mod 8   **********
Starting computation of Aq - divisors
---------
q = 79
nq(79) = 3
nq'(79) = 7
---
A(79,1) = 0.505207323331751963627968967643873591409107209148748964215
rho = 3.04209701216126450666780151334488115444920787348129074602
rho_tilde = 0.380262126520158063333475189168110144306150984185161343253
---
betan(3,79) = 0.518129296788080596043569544939347180024716847434338906307
rhonbeta = 3.15049096504302405901097604117785432217628979120422980775
rhonbetanorm = 0.393811370630378007376372005147231790272036223900528725969
---
alphan(3,79) = 0.438618380624200524969748303109148800399703288765876559432
rhonalpha = 2.56263890902555925733978510811007661265851707999689944404
rhonalphanorm = 0.320329863628194907167473138513759576582314634999612430505
---
q & nq & A(q,1) & rho=(1+A(q,1))/(1-A(q,1)) & rho_tilde
$79$ & $3$ & $0.505207323331751963627968967643873591409107209148748964215 $ & $3.04209701216126450666780151334488115444920787348129074602 $ & $0.380262126520158063333475189168110144306150984185161343253$ \\ \hline
---
nq & q & rhonalpha & rhonbeta 
$3$ & $79$ & $2.56263890902555925733978510811007661265851707999689944404 $ & $3.15049096504302405901097604117785432217628979120422980775 $ & $3.04209701216126450666780151334488115444920787348129074602$ \\ \hline
---------
Total computation time (output time included): 0 min, 0 sec, 5 millisec
Preparation time (output time included): 0 min, 0 sec, 0 millisec
A(q) product computation time (output time included): 0 min, 0 sec, 3 millisec
ALPHA product computation time (output time included): 0 min, 0 sec, 1 millisec
BETA product computation time (output time included): 0 min, 0 sec, 1 millisec
******  END PROGRAM   ********
************ A. LANGUASCO *************
********* COMPUTATION of Aq(1) - d(n) mod q **********
********* q = -1,1 mod 8   **********
Starting computation of Aq - divisors
---------
q = 89
nq(89) = 3
nq'(89) = 7
---
A(89,1) = 0.505038282406622222829193287216864449800431531750956410184
rho = 3.04071654212870891599553414693089724917009527969569702246
rho_tilde = 0.380089567766088614499441768366362156146261909961962127808
---
betan(3,89) = 0.518129296788080596043569544939347180024716847434338906307
rhonbeta = 3.15049096504302405901097604117785432217628979120422980775
rhonbetanorm = 0.393811370630378007376372005147231790272036223900528725969
---
alphan(3,89) = 0.438618380624200524969748303109148800399703288765876559432
rhonalpha = 2.56263890902555925733978510811007661265851707999689944404
rhonalphanorm = 0.320329863628194907167473138513759576582314634999612430505
---
q & nq & A(q,1) & rho=(1+A(q,1))/(1-A(q,1)) & rho_tilde
$89$ & $3$ & $0.505038282406622222829193287216864449800431531750956410184 $ & $3.04071654212870891599553414693089724917009527969569702246 $ & $0.380089567766088614499441768366362156146261909961962127808$ \\ \hline
---
nq & q & rhonalpha & rhonbeta 
$3$ & $89$ & $2.56263890902555925733978510811007661265851707999689944404 $ & $3.15049096504302405901097604117785432217628979120422980775 $ & $3.04071654212870891599553414693089724917009527969569702246$ \\ \hline
---------
Total computation time (output time included): 0 min, 0 sec, 5 millisec
Preparation time (output time included): 0 min, 0 sec, 0 millisec
A(q) product computation time (output time included): 0 min, 0 sec, 3 millisec
ALPHA product computation time (output time included): 0 min, 0 sec, 1 millisec
BETA product computation time (output time included): 0 min, 0 sec, 1 millisec
******  END PROGRAM   ********
************ A. LANGUASCO *************
********* COMPUTATION of Aq(1) - d(n) mod q **********
********* q = -1,1 mod 8   **********
Starting computation of Aq - divisors
---------
q = 97
nq(97) = 5
nq'(97) = 7
---
A(97,1) = 0.899311551807478380244703619674742835814872682030238729256
rho = 18.8632518019931600244048574235737731655600913482061283462
rho_tilde = 0.589476618812286250762651794486680411423752854631441510818
---
betan(5,97) = 0.916850779228909835416677911291295302299925792490423457856
rhonbeta = 23.0531418268609003184857276806683038586857850885743587848
rhonbetanorm = 0.720410682089403134952678990020884495583930784017948712025
---
alphan(5,97) = 0.897813044296315189281089379777004468377906648563104632058
rhonalpha = 18.5719696924866996139207933541007149680235307403926994739
rhonalphanorm = 0.580374052890209362935024792315647342750735335637271858558
---
q & nq & A(q,1) & rho=(1+A(q,1))/(1-A(q,1)) & rho_tilde
$97$ & $5$ & $0.899311551807478380244703619674742835814872682030238729256 $ & $18.8632518019931600244048574235737731655600913482061283462 $ & $0.589476618812286250762651794486680411423752854631441510818$ \\ \hline
---
nq & q & rhonalpha & rhonbeta 
$5$ & $97$ & $18.5719696924866996139207933541007149680235307403926994739 $ & $23.0531418268609003184857276806683038586857850885743587848 $ & $18.8632518019931600244048574235737731655600913482061283462$ \\ \hline
---------
Total computation time (output time included): 0 min, 0 sec, 6 millisec
Preparation time (output time included): 0 min, 0 sec, 0 millisec
A(q) product computation time (output time included): 0 min, 0 sec, 4 millisec
ALPHA product computation time (output time included): 0 min, 0 sec, 1 millisec
BETA product computation time (output time included): 0 min, 0 sec, 1 millisec
******  END PROGRAM   ********
  ***   last result: cpu time 61 ms, real time 80 ms.
Goodbye!

------
q <= 2*10^4
------
gp2c-run -pmy_ -g -W Adivisors-file-v6.gp <istruzione-macarmM1.txt >results_rho_div_20000.txt &

? forprime(q=2,2*10^4,qmod8=q%8;if(qmod8==1||qmod8==7,init_Adivisors_file_v6();Aqdivisors(q,57)))

Output results on the files:
- results_rho_div_20000.txt
- rhodivfile-prodeulerrat-20000.csv

---------------------
 
default(parisizemax,"12G")
  
  ***   Warning: new maximum stack size = 12000002048 (11444.094 Mbytes).

? init_Adivisors_file_v6();Aqdivisors(31391,57) 
************ A. LANGUASCO *************
********* COMPUTATION of Aq(1) - d(n) mod q **********
********* q = -1,1 mod 8   **********
Starting computation of Aq - divisors
---------
  *** Aqdivisors: Warning: increasing stack size to 16000000.
  *** Aqdivisors: Warning: increasing stack size to 32000000.
  *** Aqdivisors: Warning: increasing stack size to 64000000.
  *** Aqdivisors: Warning: increasing stack size to 128000000.
  *** Aqdivisors: Warning: increasing stack size to 256000000.
  *** Aqdivisors: Warning: increasing stack size to 512000000.
  *** Aqdivisors: Warning: increasing stack size to 1024000000.
  *** Aqdivisors: Warning: increasing stack size to 2048000000.
  *** Aqdivisors: Warning: increasing stack size to 4096000000.
  *** Aqdivisors: Warning: increasing stack size to 8192000000.
  *** Aqdivisors: Warning: increasing stack size to 12000002048.


q = 31391
nq(31391) = 31
nq'(31391) = 37
---
A(31391,1) = 0.999999999052967941217892165576092177279347796168204757764
rho = 2111860924.35454277627683908925278299227080590583044322132
rho_tilde = 0.983411876649848546961713158177576489872674413877745727690
---
betan(31,31391) = 0.999999999068670949468476627417627833810040502670308532622
rhonbeta = 2147468714.65809118242180480451320876213321353346670497867
rhonbetanorm = 0.999993046120782933394333722308840957532270594251670399061
---
alphan(31,31391) = 0.999999999052300035788774249587625699316406215870501990001
rhonalpha = 2110372559.43858515148331915906495221602703184160666055006
rhonalphanorm = 0.982718802727100044248308594834502886993359700593473645888
---
q & nq & A(q,1) & rho=(1+A(q,1))/(1-A(q,1)) & rho_tilde
$31391$ & $31$ & $0.999999999052967941217892165576092177279347796168204757764 $ & $2111860924.35454277627683908925278299227080590583044322132 $ & $0.983411876649848546961713158177576489872674413877745727690$ \\ \hline
---
nq & q & rhonalpha & rhonbeta
$31$ & $31391$ & $2110372559.43858515148331915906495221602703184160666055006 $ & $2147468714.65809118242180480451320876213321353346670497867 $ & $2111860924.35454277627683908925278299227080590583044322132$ \\ \hline
---------
Total computation time (output time included): 6 min, 44 sec, 620 millisec
Preparation time (output time included): -1 min, 59 sec, 995 millisec
A(q) product computation time (output time included): 6 min, 44 sec, 625 millisec
ALPHA product computation time (output time included): 0 min, 0 sec, 0 millisec
BETA product computation time (output time included): 0 min, 0 sec, 0 millisec
******  END PROGRAM   ********
? ##
  ***   last result: cpu time 6min, 44,636 ms, real time 19min, 444 ms.
 
**********/
 
  
/********
---------
I7 assembled machine with Ubuntu
---------
                  GP/PARI CALCULATOR Version 2.17.2 (released)
          amd64 running linux (x86-64/GMP-6.3.0 kernel) 64-bit version
    compiled: Mar  7 2025, gcc version 13.3.0 (Ubuntu 13.3.0-6ubuntu2~24.04)
                   threading engine: pthread, nbthreads = 24
                (readline v8.2 disabled, extended help enabled)

                     Copyright (C) 2000-2024 The PARI Group

PARI/GP is free software, covered by the GNU General Public License, and comes 
WITHOUT ANY WARRANTY WHATSOEVER.

Type ? for help, \q to quit.
Type ?18 for how to get moral (and possibly technical) support.

parisize = 8000000, primelimit = 1048576, factorlimit = 1048576

 
--------
q <= 10^5
--------

gp2c-run -pmy_ -g -W Adivisors-file-v6.gp <istruzione-languasco0.txt >results_rho_div_100000.txt &

default(parisizemax,"58G")
forprime(q=2,10^5,qmod8=q%8;if(qmod8==1||qmod8==7,init_Adivisors_file_v6();Aqdivisors(q,57)))

Output results on the files:
- results_rho_div_100000.txt
- rhodivfile-prodeulerrat-100000.csv