/*** A. LANGUASCO & A. ZACCAGNINI ***/ /*** Implementation of Pintz-Ruzsa's method paper ***/ /*** (as described in their paper paper on Acta Arith. 109,(2003)) ***/ /*** Improved by K. BELABAS using the gexpo function ***/ /*** Improved by K. BELABAS using a more clever matrices generation ***/ /*** this version is about 15% faster for small precisions ***/ /*** and 5% faster for large precisions ***/ global(L, L1, majpsiU, minpsiU, majpsiL, minpsiL, momentmajU, momentminU, momentmajL, momentminL); global(SIN, COS); \\ Goal function and its first derivative : gam e gamprime \\ gamma is internally used for the Euler Gamam function gam(i,lambda) = exp(lambda*COS[i+1]); gamprime(i,lambda) = { local(g = gam(i,lambda)); [g, -lambda*SIN[i+1]*g]; } \\ gexpo is a libpari function; returns the binary exponent of x or \\ the maximal binary exponent of the coefficients of x. \\ lg is a libpari function; returns the length of x install(gexpo, lG); \\ main function: it realizes the approximation via the matrix method \\ c is the level \\ k is the degree of the polynomials \\ lambda is the point in which we evaluate psi \\ numdigits is the precision for the final result {PintzRuzsa(c, k, lambda, numdigits) = \\ local variables local(log2 = log(2)); local(c1 = c*log2); local(eU, eL, expU, expL, eps2, K2, g1, g2); local(sigmaupper, sigmalower, sigma0upper,sigma0lower, rhoupper, delta, rholower, Uupper, Ulower, Bupper, Blower, K, j, l, eps1, i, i1,i2, m, j1, j2,alphaupp, alphalow, betalow,betaupp); \\ Building the coefficients of the polynomials K = k+1; K2 = 2*K; sigma0upper = 0; sigmaupper=vector(k); rhoupper=vector(k); sigma0lower = 0; sigmalower=vector(k); rholower=vector(k); for(i=1,K, g1 = gamprime(2*i-2, lambda); sigma0upper += g1[1]; g2 = gamprime(2*i-1, lambda); sigma0lower += g2[1]; for (j=1, k, sigmaupper[j] += g1[1] * COS[(j*(2*i-2)) % K2 + 1]; sigmalower[j] += g2[1] * COS[(j*(2*i-1)) % K2 + 1]; rhoupper[j] += g1[2] * SIN[(j*(2*i-2)) % K2 + 1]; rholower[j] += g2[2] * SIN[(j*(2*i-1)) % K2 + 1]; ); ); sigma0upper /= K; sigmaupper /= K; rhoupper /= K; sigma0lower /= K; sigmalower /= K; rholower /= K; \\ print(sigma0upper); \\ print(sigma0lower); \\ print(sigmaupper); \\ print(sigmalower); \\ print(rhoupper); \\ print(rholower); Bupper=vector(k+1); Blower=vector(k+1); Bupper[1]=sigma0upper; Blower[1]=sigma0lower; for (i=1, k, Bupper[i+1]= (2/K) * ((K-i)*sigmaupper[i]-rhoupper[i]); Blower[i+1]= (2/K) * ((K-i)*sigmalower[i]-rholower[i]); ); \\ printp("Bupper =", Bupper); \\ printp("Blower =", Blower); \\ Building the matricies U m=k-1; Uupper = matrix(m+1,m+1); Ulower = matrix(m+1,m+1); \\ print("k= ", k); \\ This definition follows from line -6 and -5 at p.178 pf Pintz-Ruzsa paper. for (j=1,m+1, forstep (n = (j-1)%2, m, 2, j1=(j-1+n)\2; j2= abs(j-1-n)\2; Uupper[n+1,j1+1] += Bupper[j]/2; Uupper[n+1,j2+1] += Bupper[j]/2; Ulower[n+1,j1+1] += Blower[j]/2; Ulower[n+1,j2+1] += Blower[j]/2; ) ); \\ printp("Uupper =", Uupper); \\ printp("Ulower =", Ulower); \\ eps1 is the desired precision eps1= 10^(-(numdigits)); \\ print(eps1); L = 0; \\ expU and expL are used to keep track of the exponents used expU = 0; expL = 0; eps2 = eps1^3; until (delta>= eU; expU += eU); eL = gexpo(Ulower); if (eL, Ulower >>= eL; expL += eL); \\ second normalization step: the entry in the lower right corner \\ decreases too quickly comparing to the others if (abs(Uupper[m+1,m+1]) < eps2, Uupper[m+1,m+1] = 0); if (abs(Ulower[m+1,m+1]) < eps2, Ulower[m+1,m+1] = 0); \\ squaring out the matrices and keeping track of the exponents Uupper = Uupper^2; expU *= 2; Ulower = Ulower^2; expL *= 2; \\ keeping track of the number of squares L++; \\ approximations, see Lemma 7, p. 179, of Pintz-Ruzsa paper betaupp = sum(i=1,m+1, Uupper[1,i]); alphaupp = Uupper[1,1]; \\ computing the precision delta = log(betaupp / alphaupp) / (lambda*2^L); \\ print("delta= ", delta); ); alphalow = Ulower[1,1]; betalow = sum(i=1,m+1, Ulower[1,i]); \\ computing the moments and the approximations L1 = 2^L; momentmajU = (log(betaupp) + expU * log2)/L1; momentminU = (log(alphaupp)+ expU * log2)/L1; momentmajL = (log(betalow) + expL * log2)/L1; momentminL = (log(alphalow)+ expL * log2)/L1; \\ computing the approximations for the function psi of Pintz-Ruzsa paper majpsiU= (momentmajU + c1)/lambda; minpsiU= (momentminU + c1)/lambda; majpsiL= (momentmajL + c1)/lambda; minpsiL= (momentminL + c1)/lambda; \\ if the degree of the polynomials is too small we are not able to reach \\ the fixed precision if (abs(majpsiU-minpsiL) > eps1, error("increase the degree of the polynomials") ); \\ return the approximated value of psi(lambda) return(majpsiU); } \\ Auxiliary functions needed to generate the sin and cos values used \\ to build the approximation matrices /* COS[i] = cos((i-1)*Pi/K), SIN[i] = sin((i-1)*Pi/K) */ initSINCOS(K) = { COS = vector(2*K); SIN = vector(2*K); COS[1] = 1; COS[2] = cos(Pi/K); SIN[1] = 0; SIN[2] = sin(Pi/K); for (i = 3, K, /* could use 3M/Karatsuba to gain 1 mul*/ COS[i] = COS[i-1]*COS[2] - SIN[i-1]*SIN[2]; SIN[i] = COS[i-1]*SIN[2] + SIN[i-1]*COS[2]; ); for (i = K+1, 2*K, COS[i] = -COS[i-K]; SIN[i] = -SIN[i-K]; ); } \\ This function realizes a dyadic approximation using the previous function \\ it starts using the approximated values of psi in 1-T, 1, 1+T \\ and it moves to a new set of three points to look for a saddle point \\ which is approximated using a precision 10^(-numdigits) \\ c is the level \\ k is the degree of the polynomials \\ numdigits is the precision for the final result {PintzRuzsa_psiapprox(c,k,numdigits) = local(iteration = 0,eps0,r0,r1,r2,l0,T); numdigits += 2; default(realprecision,numdigits+10); \\ set working precision; print("The expected number of correct decimal digits is = ", numdigits-2); eps0=10^(-numdigits); print("Approximation for the minimal lambda is = 10^(-", numdigits-2,")"); r0=0; r1=0; r2=0; T=1/10; l0=1; initSINCOS(k+1); until (T