/* Copyright (C) 2023 Alessandro Languasco */ \\ default(parisizemax,"32G"); /***************************************** GLOBAL VARIABLES : in uppercase. USED TO TRACK THE FLOW OR TO SHARE INFORMATION ABOUT THE PROBLEM ****************************************/ \\global variables global(TOT_ITER, ITER_TRIVIAL_SECOND, ITER_TRIVIAL_FIRST, ITER_CONTFRAC, PN, QN, EXTRA); global(THETA, ONEOVERTHETA, LOGALPHA, LOGRATIO, LOGONEOVERALPHA, MAXVALUE, MINVALUE); global(C1_GLOBAL, C2_GLOBAL, D1_GLOBAL, D2_GLOBAL); global(UMAX_GLOBAL_C1, VMAX_GLOBAL_C1, PMAX_GLOBAL_C1, QMAX_GLOBAL_C1); global(UMIN_GLOBAL_C2, VMIN_GLOBAL_C2, PMIN_GLOBAL_C2, QMIN_GLOBAL_C2); global(UMAX_GLOBAL_D1, VMAX_GLOBAL_D1, PMAX_GLOBAL_D1, QMAX_GLOBAL_D1); global(UMIN_GLOBAL_D2, VMIN_GLOBAL_D2, PMIN_GLOBAL_D2, QMIN_GLOBAL_D2); /***************************************** MAIN FUNCTION; the one to be called ****************************************/ {constants_pq_left_right(boundprime, LOOP = 10^4) = local(millisec, minutes, seconds, npqalpha, ok, name, namefile, seqfile, logp, logq, u, v , uleft, vleft, string, ustart, vstart, time, ubad, vbad, starttime, elaptime, uend, vend, timebound, partelaptime, partialtime, exp_vec, found, logten, digits_n_start, oneminute, prev_elapsed, prec, MAX_constant,min_constant, Sums, logn, lognleft, num, denom, value, pmin, qmin, pmax, qmax, umin, vmin, umax, vmax, aux, alpha, fine, uright, vright, aux1, logtwo, non_negative_rho, N1, N, logN ); \\accuracy = 15; prec = 19; default(realprecision,prec); default(format , f); \\ print numbers with floating point notation (not E notation) \\eps = 1.*10^(-accuracy); timebound = 7*86400000; \\ one week in milliseconds \\timebound = 86400000; \\ one day in milliseconds \\timebound = 3600000; \\ one hour in milliseconds \\timebound = 60000; \\ one minute in milliseconds oneminute = 60000; seconds = floor(timebound /1000)%60; minutes = floor(timebound /oneminute); millisec = timebound - minutes*oneminute - seconds*1000; print("**********"); print("TIMEBOUND for this computation is : ", minutes, " min, ", seconds, " sec, ", millisec, " millisec"); print("**********"); starttime = getwalltime(); print("---------"); print("Searching for C1, C2, D1, D2 "); print("with the continued fractions algorithm; left neighbors and right neighbors" ); print("---------"); name="constants_pq_data"; namefile = Strprintf("%s-%s-%s.csv", name, boundprime, LOOP); seqfile = fileopen(namefile, "w"); string = "p;q;TOT_ITER;non_negative_rho;LOOP;D1;D2;C1;C2;C3;C4;logp*logq;Sums;mu"; filewrite(seqfile, string); /*** SET THE GLOBAL VARIABLE LOGALPHA, LOGONEOVERALPHA ********/ alpha = 2; LOGALPHA = log(alpha); LOGONEOVERALPHA = - LOGALPHA; logten = log(10); logtwo = log(2); N1 = 3; \\ extra number of neighbors EXTRA = 10; MAXVALUE = 10^1000000; MINVALUE = - MAXVALUE; C1_GLOBAL = MINVALUE; C2_GLOBAL = MAXVALUE; D1_GLOBAL = MINVALUE; D2_GLOBAL = MAXVALUE; exp_vec = vector(3); \\ meaning: exp_vec[1] = u ; exp_vec[2] = v ; exp_vec[3] = found ; prev_elapsed = 0; forprime(p=2, boundprime, logp = log(p); forprime(q=p+1,boundprime, \\ starting point ITER_CONTFRAC = 0; ITER_TRIVIAL_FIRST = 0; ITER_TRIVIAL_SECOND = 0; TOT_ITER = 0; MAX_constant = MINVALUE; min_constant = MAXVALUE; Sums = 0; logq = log(q); N = max(N1,exp(logp*logq)); logN =log(N); \\ getting the starting point n_start exp_vec = get_starting_point(exp_vec, logp, logq); \\ n_start is p^exp_vec[1] * q^exp_vec[2] ustart = exp_vec[1]; vstart = exp_vec[2]; u = ustart; v = vstart; print("*************************************"); print("p = ", p, "; q = ", q,"; n = p^u*q^v; u,v >= 0"); print("n_i >= N = max ( N1, exp((log p) * (log q)); N1 = ", N1, "; N = ", N ); print("--------------"); print("Searches starting with [",p,",",q,"]-exponents: u = ", ustart, ", v = ", vstart); digits_n_start = floor(ustart *logp/logten + vstart*logq/logten)+1; print ("n_start = ", p^ustart*q^vstart); print ("n_start has ", digits_n_start, " digits"); npqalpha = 0; \\ initialization ITER_CONTFRAC += 1; \\ updating the iterations counter ok = 0 ; \\ when ok = 1 we got a (possibly tight) upper bound for npq(alpha) non_negative_rho = 0; while ( ok == 0 && TOT_ITER < LOOP, \\ with ok = 1: we did not find any left neighbor \\ with ok = 0: we found a left neighbor; looking for another one. /******** checking running time **********/ partialtime = getwalltime; partelaptime = partialtime - starttime; if (partelaptime >= (prev_elapsed + oneminute), print("****************************************************"); minutes = floor(partelaptime/oneminute); print("Computing .... elapsed time = ", minutes, " minutes"); print("****************************************************"); prev_elapsed = partelaptime; ); if (partelaptime > timebound , seconds = floor(timebound /1000)%60; minutes = floor(timebound /oneminute); millisec = timebound - minutes*oneminute - seconds*1000; print("TIMEBOUND for this computation was : ", minutes, " min, ", seconds, " sec, ", millisec, " millisec"); error (" The computation is taking too much time - aborted") ); \\ searching the left neighbor with the continued fractions exp_vec[1] = u; exp_vec[2] = v; found = 0 ; exp_vec = cont_frac_left_search(exp_vec, logp, logq); \\ if found == 1; the cont frac search was fine uleft = exp_vec[1]; vleft = exp_vec[2]; found = exp_vec[3]; lognleft = uleft * logp + vleft * logq; if (found == 1 && uleft <> u && vleft <> v && lognleft >= logN, \\found a left neighbor in the cont frac search \\ that satisfies Bertrand's postulate condition \\ we can generate another left element \\ setting data for the next iteration \\ we compute the tho_i if (uleft <> 0 || vleft <> 0, logn = u * logp + v * logq; LOGRATIO = logn - lognleft; if (LOGRATIO <= logtwo, non_negative_rho +=1; num = - log( exp(LOGRATIO) - 1 ); denom = log( lognleft ); value = num / denom ; if (value < min_constant, min_constant = value; pmin = p; qmin = q; umin=uleft; vmin =vleft ); if (value > MAX_constant, MAX_constant = value; pmax = p; qmax = q; umax=uleft; vmax =vleft ); Sums += value; ); ); ok = 0; \\ will stay in the external while ITER_CONTFRAC +=1 ; u = uleft; v = vleft; msg_iterations(p,q,u,v,ITER_CONTFRAC,5000, "continued fraction") , \\ else \\ else: not found a good left neighbor in the cont frac search ubad = exp_vec[1]; vbad = exp_vec[2]; \\ initialise the trivial search exp_vec[1] = u; exp_vec[2] = v; found = 0 ; exp_vec = first_trivial_search_closest(exp_vec, ubad, vbad, logp, logq, p, q); found = exp_vec[3]; if (found == 0, \\ not found a good left neighbor in the first trivial search \\ starting the second trivial search \\ ubad and vbad comes from before exp_vec[1]=u; exp_vec[2]=v; found = 0 ; exp_vec = second_trivial_search_closest(exp_vec, ubad, vbad, logp, logq, p, q, alpha); found = exp_vec[3]; if (found == 1 && exp_vec[1] == 0 && exp_vec[2] == 0, uend = exp_vec[1]; vend = exp_vec[2]; \\ will leave in the external while ok = 1, \\else if (found == 1 && (exp_vec[1] <> 0 || exp_vec[2] <> 0), \\ found a good left neighbor in the second trivial search \\ looking for another one. will stay in the external while \\ we compute the tho_i uleft = exp_vec[1]; vleft = exp_vec[2]; \\print("lognleft second trivial"); lognleft = uleft * logp + vleft * logq; if ((uleft <> 0 || vleft <> 0) && lognleft >= logN, logn = u * logp + v * logq; LOGRATIO = logn - lognleft; if (LOGRATIO <= logtwo, non_negative_rho +=1; num = - log( exp(LOGRATIO) - 1 ); denom = log( lognleft ); value = num / denom ; if (value < min_constant, min_constant = value; pmin = p; qmin = q; umin=uleft; vmin =vleft ); if (value > MAX_constant, MAX_constant = value; pmax = p; qmax = q; umax=uleft; vmax =vleft ); Sums += value; ); ); u = uleft; v = vleft; ok = 0; msg_iterations(p,q,u,v,ITER_TRIVIAL_SECOND,10, "second trivial") , \\ else: not found a good left neighbor in the second trivial search \\ we have determined n_0 = p^u*q^v uend = u; vend = v; \\ will leave in the external while ok = 1 ) ) , \\else: found a left neighbor in the first trivial search \\ looking for another one. will stay in the external while uleft = exp_vec[1]; vleft = exp_vec[2]; lognleft = uleft * logp + vleft * logq; if ((uleft <> 0 || vleft <> 0) && lognleft >= logN, logn = u * logp + v * logq; LOGRATIO = logn - lognleft; if (LOGRATIO <= logtwo, non_negative_rho +=1; num = - log( exp(LOGRATIO) - 1 ); denom = log( lognleft ); value = num / denom ; if (value < min_constant, min_constant = value; pmin = p; qmin = q; umin=uleft; vmin =vleft ); if (value > MAX_constant, MAX_constant = value; pmax = p; qmax = q; umax=uleft; vmax =vleft ); Sums += value; ); ); u = uleft; v = vleft; ok = 0; msg_iterations(p,q,u,v,ITER_TRIVIAL_FIRST,1000, "first trivial"); ); ); TOT_ITER = ITER_TRIVIAL_FIRST+ITER_TRIVIAL_SECOND+ITER_CONTFRAC; ); if (TOT_ITER >= LOOP, fine = 1, fine = 0 ); \\print("fine = ", fine); if (fine == 1, \\print("Number of non negative rho = ", non_negative_rho); constants_comp(seqfile, p,q, Sums, TOT_ITER, non_negative_rho, LOOP, min_constant, MAX_constant, pmin, qmin, pmax, qmax, umin, vmin, umax, vmax, logp, logq ); , \\else \\fine == 0; we need to generate rightneighbors print("starting right neighbors generation"); u = ustart; v = vstart; \\print("ustart = ", ustart); \\print("vstart = ", vstart); for ( jj = 1, LOOP - TOT_ITER + EXTRA, \\ print(jj); \\ right neighbors generation \\print("u = ", u); \\print("v = ", v); exp_vec[1] = u; exp_vec[2] = v; exp_vec = cont_frac_right_search(exp_vec, logp, logq); \\print("uright = ", uright); \\print("vright = ", vright); found = exp_vec[3]; if (found == 0, error ("wring right cfrac neighbor generation")); uright = exp_vec[1]; vright = exp_vec[2]; aux = u * logp + v * logq ; aux1= uright * logp + vright * logq ; \\ delta = element * (ratio -1); log( delta ) = log(element) + log(ratio-1) = aux + log(ratio-1) \\num = aux - log( delta ); if (LOGRATIO <= logtwo, non_negative_rho +=1; num = - log( exp(LOGRATIO) - 1 ); denom = log( aux ); value = num / denom ; if (value < min_constant, min_constant = value; pmin = p; qmin = q; umin=u; vmin =v ); if (value > MAX_constant, MAX_constant = value; pmax = p; qmax = q; umax=u; vmax =v ); Sums += value; ); \\ right neighbor exponents \\ start again u = uright ; v = vright ; );\\ end jj loop ); \\end if fine \\print("Number of non negative rho = ", non_negative_rho); constants_comp(seqfile, p,q, Sums, TOT_ITER, non_negative_rho, LOOP, min_constant, MAX_constant, pmin, qmin, pmax, qmax, umin, vmin, umax, vmax, logp, logq ); ); ); print("*************************************"); fileflush(seqfile); fileclose(seqfile); final_output(); time = getwalltime() - starttime; elaptime = time; seconds = floor(elaptime /1000)%60; minutes = floor(elaptime /oneminute); millisec = elaptime - minutes*oneminute - seconds*1000; print("Constants_pq computation time (output time included): ", minutes, " min, ", seconds, " sec, ", millisec, " millisec"); print("****** END PROGRAM ********"); } \\ Searches function /***************************************** SECOND TRIVIAL SEARCH (using also data from failed iteration of the continued fraction search) ****************************************/ {second_trivial_search_closest(exp_vec, ubad, vbad, logp, logq, p, q, alpha) = local(u1, v1, logn, lbound, ubound, z, y, astart, bstart, aend, bend, amax, bmax, ymax); \\ logN = ubad* logp + vbad * logq; \\ we know that logN < log(n/alpha); otherwise we do not enter this routine \\ we look for (a,b) such that at least one of them is > than the corresponding element \\ of the pair (ubad , vbad); so that for sure a * logp + b * logq > logN u1 = exp_vec[1]; v1 = exp_vec[2]; if (u1 == 0 && v1 == 0, exp_vec[3] = 1; return(exp_vec); ); exp_vec[3] = 0; logn = u1 * logp + v1 * logq; lbound = log( floor (p^u1*q^v1/alpha) ); ymax = lbound-1; ubound = logn; \\ b < = logn /logq = u1*THETA + v1 \\bstart = vbad + 1; bstart = vbad; bend = floor(u1*THETA) + v1 ; astart = 0; for ( b = bstart, bend, z = b * logq ; if (z >= ubound, \\exp_vec[3] = 0; break(1) ); \\ a <= (logn - b*logq)/logp = logn/logp - b/THETA = u1 + v1/THETA - b/THETA = u1 + (v1-b)/THETA aend = floor((v1-b)*ONEOVERTHETA) + u1; for (a = astart, aend , y = z + a * logp; if (y >= ubound, \\exp_vec[3] = 0; break(1), \\ else y= lbound && y > ymax, \\print("*** second trivial search found one element !!!"); ymax = y; amax = a; bmax = b ; exp_vec[3] = 1; \\ITER_TRIVIAL_SECOND+=1; \\return(exp_vec); ); ); ); \\endfor a );\\endfor b \\ exp_vec[3] == 0: for sure it is \\ not found a left neighbor in the second trivial search (upper part) \\ having b>vbad, and a "free" \\ now we test a > ubad and b <= vbad bstart = 0; bend = vbad; astart = ubad+1; \\astart = ubad; for (b = bstart , bend, z = b * logq ; if (z >= ubound, \\exp_vec[3] = 0; break(1) ); \\ a <= (logn - b*logq)/logp = logn/logp - b/THETA = u1 + v1/THETA - b/THETA = u1 + (v1-b)/THETA aend = floor((v1-b)*ONEOVERTHETA) + u1; for (a = astart, aend, y = z + a * logp; if (y >= ubound, \\exp_vec[3] = 0; break(1), \\ else y < ubound if (y >= lbound && y > ymax, \\print("*** second trivial search found one element !!!"); ymax = y; amax = a; bmax = b ; exp_vec[3] = 1; ); ); ); \\endfor a );\\endfor b if (exp_vec[3] == 1, ITER_TRIVIAL_SECOND+=1; exp_vec[1] = amax; exp_vec[2] = bmax , \\else exp_vec[1] = 0; exp_vec[2] = 0; ); return(exp_vec); } /***************************************** FIRST TRIVIAL SEARCH (using also data from failed iteration of the continued fraction search) ****************************************/ {first_trivial_search_closest(exp_vec, ubad, vbad, logp, logq, p, q) = local(u1,v1, logn, lbound, ubound, z, y, astart, bstart, aend, bend, ymax, amax, bmax); \\ we know that logN < log(n/alpha); otherwise we do not enter this routine \\ we look for (a,b) such that at least one of them is > than the corresponding element \\ of the pair (ubad, vbad); so that for sure a * logp + b * logq > logN u1 = exp_vec[1]; v1 = exp_vec[2]; exp_vec[3] = 0; logn = u1 * logp + v1 * logq; lbound = logn + LOGONEOVERALPHA; ymax = lbound-1; ubound = logn; \\ THETA = logp/logq; \\bstart = vbad + 1; bstart = vbad; bend = floor(u1*THETA) + v1 ; astart = 0; for ( b = bstart, bend, z = b * logq ; if (z >= ubound, \\exp_vec[3] = 0; break(1) ); \\ a <= (logn - b*logq)/logp = logn/logp - b/THETA = u1 + v1/THETA - b/THETA = u1 + (v1-b)/THETA aend = floor((v1-b)*ONEOVERTHETA) + u1; for (a = astart, aend , y = z + a * logp; if (y >= ubound, \\exp_vec[3] = 0; break(1), \\ else y= lbound && y > ymax, \\print("*** first trivial search found one element !!!"); ymax = y; amax = a; bmax = b ; exp_vec[3] = 1; ); ); ); \\endfor a );\\endfor b \\ exp_vec[3] == 0, \\ it surely is \\ not found a left neighbor in the first trivial search (upper part) \\ having b>vbad, and a "free" \\ now we test a > ubad and b <= vbad bstart = 0; bend = vbad; astart = ubad+1; for (b = bstart , bend, z = b * logq ; if (z >= ubound, \\exp_vec[3] = 0; break(1) ); \\ a <= (logn - b*logq)/logp = logn/logp - b/THETA = u1 + v1/THETA - b/THETA = u1 + (v1-b)/THETA aend = floor((v1-b)*ONEOVERTHETA) + u1; for (a = astart, aend, y = z + a * logp; if (y >= ubound, \\exp_vec[3] = 0; break(1), \\ else y < ubound if (y >= lbound && y > ymax, \\print("*** first trivial search found one element !!!"); ymax = y; amax = a; bmax = b ; exp_vec[3] = 1; ); ); ); \\endfor a );\\endfor b if (exp_vec[3] == 1, ITER_TRIVIAL_FIRST+=1; exp_vec[1] = amax; exp_vec[2] = bmax , \\else exp_vec[1] = 0; exp_vec[2] = 0; ); return(exp_vec); } /***************************************** CONTINUED FRACTION SEARCHES ****************************************/ {cont_frac_right_search(exp_vec, logp, logq) = local(u_local, v_local, u1, v1, u2, v2, index, x, x1, x2, uright, vright ); \\ GOING FOR THE RIGHT NEIGHBOR IN THE n_i SEQUENCE \\ initalizations u1=0; v1=0; u2=0; v2=0; x1=0; x2=0; x=0; u_local = exp_vec[1]; v_local = exp_vec[2]; \\ Searching for a denominator <= u_local index = setsearch(QN, u_local); \\ if u_local is not in the list of denominators, setsearch (with two parameters) returns 0 if (index == 0, \\ u_local was not in the list of denominators; setsearch (with three parameters) \\ returns the position in which u should be inserted in such a list \\ so this value minus 1 is the largest denominator < u_local index = setsearch(QN, u_local , 1 ) - 1 ); \\ upper convergents; p_n/q_n ; n odd; stored with an even index by pari/gp if ( index <= 1, \\ code for not found is the 0 pair \\print("index upper"); v1 = 0; u1 = 0; exp_vec[3]= 0 , if ( index % 2 == 0 , \\ index is even; it stores an odd convergent v1 = QN[index]; u1 = PN[index], \\ index is odd; it stores an even convergent \\ so the largest odd convergent is in position (index - 1) v1 = QN[index - 1]; u1 = PN[index - 1]; ); ); \\ Searching for a numerator <= v_local index = setsearch(PN, v_local); \\ if v_local is not in the list of numerators, setsearch (with two parameters) returns 0 if (index == 0, \\ v was not in the list of numerators; setsearch (with three parameters) \\ returns the position in which v_local should be inserted in such a list \\ so this value minus 1 is the largest numerator < v_local index = setsearch(PN, v_local , 1 ) - 1 ); \\ lower convergents; p_n/q_n ; n even; stored with an odd index by pari/gp if ( index <= 0, \\ code for not found is the 0 couple \\print("index lower"); v2 = 0; u2 = 0; exp_vec[3]= 0 , if ( index % 2 == 1 , \\ index is odd; it stores an even convergent u2 = PN[index]; v2 = QN[index], \\ index is even; it stores an odd convergent; \\ so the largest even convergent is in position (index - 1) u2 = PN[index - 1]; v2 = QN[index - 1]; ); ); x1 = u1 * logq - v1 * logp ; x2 = v2 * logp - u2 * logq ; if (x1 < 0 || x2 < 0, error("wrong cfrac approx - right neighbors "); ); x = abs(x1) - abs(x2); if (x1 == 0 && x2 == 0, \\ [both pairs are (0,0); skipping] exp_vec[3] = 0; return(exp_vec), \\ else [x <> 0] \\ x = abs (v1*logp - u1*logq) - abs(v2*logp - u2*logq); if (x1 == 0, uright = u_local + v2; vright = v_local - u2; LOGRATIO = x2; exp_vec[3] = 1 , \\else if (x2 == 0, uright = u_local - v1; vright = v_local + u1; LOGRATIO = x1; exp_vec[3] = 1 ,\\else if ( x < 0 , uright = u_local - v1; vright = v_local + u1; LOGRATIO = x1; exp_vec[3] = 1 \\rightelement = element * ratio \\else , uright = u_local + v2; vright = v_local - u2; LOGRATIO = x2; exp_vec[3] = 1 \\rightelement = element * ratio; ); ); ); ); \\ end right neighbor generation exp_vec[1] = uright; exp_vec[2] = vright; return(exp_vec); } {cont_frac_left_search(exp_vec, logp, logq) = local( u_local, v_local, u1, v1, u2, v2, index, x1, x2, xmax, xmin, usummand_max, vsummand_max \\ usummand_min, vsummand_min, ); \\ GOING FOR THE LEFT NEIGHBOR IN THE n_i SEQUENCE \\ initalizations u1=0; v1=0; u2=0; v2=0; x1=0; x2=0; \\ this is n_i = p^u_local*q^v_local u_local = exp_vec[1]; v_local = exp_vec[2]; \\ SEARCHING FOR AN LOWER CONVERGENT \\ looking for a denominator <= u_local index = setsearch(QN, u_local); \\ if u_local is not in the list of denominators, setsearch (with two parameters) returns 0 if (index == 0, \\ u was not in the list of denominators; setsearch (with three parameters) \\ returns the position in which u should be inserted in such a list \\ so this value minus 1 is the index of the largest denominator < u index = setsearch(QN, u_local , 1 ) - 1; ); \\ lower convergents; p_n/q_n ; n even; stored with an odd index by pari/gp if ( index <= 0, \\print("index - error; lower convergent"); exp_vec[1] = 0; exp_vec[2] = 0; exp_vec[3] = 0; x1 = 0 , \\else exp_vec[3] = 1; if ( index % 2 == 1 , \\ index is odd; it stores an even convergent; good u1 = QN[index]; v1 = PN[index], \\else \\ index is even; it stores an odd convergent; \\ so the largest even convergent is in position (index - 1) u1 = QN[index - 1]; v1 = PN[index - 1]; ); x1 = -u1*logp + v1*logq; ); \\ SEARCHING FOR AN UPPER CONVERGENT \\ looking for a numerator <= v_local index = setsearch(PN, v_local); \\ if v_local is not in the list of numerators, setsearch (with two parameters) returns 0 if (index == 0, \\ v_local was not in the list of numerators; setsearch (with three parameters) \\ returns the position in which v should be inserted in such a list \\ so this value minus 1 is the index of the largest numerator < v index = setsearch(PN, v_local, 1 ) - 1 ); \\ upper convergents; p_n/q_n ; n odd; stored with an even index by pari/gp if ( index <= 1, \\print("index - error; upper convergent"); exp_vec[1] = 0; exp_vec[2] = 0; exp_vec[3] = 0; x2 = 0 , \\else exp_vec[3] = 1; if ( index % 2 == 0 , \\ index is even; it stores an odd convergent; good u2 = QN[index]; v2 = PN[index], \\ else \\ index is odd; it stores an even convergent \\ so the largest odd convergent is in position (index - 1) u2 = QN[index - 1]; v2 = PN[index - 1]; ); x2 = u2*logp - v2*logq; ); if ( exp_vec[3] == 1, xmax = max(x1,x2); xmin = min(x1,x2); if (xmin == x1, \\usummand_min = -u1; \\vsummand_min = v1; usummand_max = u2; vsummand_max = -v2, \\else \\usummand_min = u2; \\vsummand_min = -v2; usummand_max = -u1; vsummand_max = v1 ); if( xmax < LOGONEOVERALPHA, exp_vec[1] = u_local + usummand_max; \\ ubad; needed for the first trivial search exp_vec[2] = v_local + vsummand_max; \\ vbad; needed for the first trivial search exp_vec[3] = 0; \\ not found a good point return(exp_vec) , \\else we always choose xmax exp_vec[1] = u_local + usummand_max; \\ uleft exp_vec[2] = v_local + vsummand_max; \\ vleft exp_vec[3] = 1; \\ found a good point return(exp_vec) ); ); return(exp_vec); } /***************************************** DETECT THE STARTING POINT OF THE SEQUENCE (using the continued fraction expansion of log p/ log q) ****************************************/ {get_starting_point(exp_vec, logp, logq) = local( cfrac, y1, index, a2, b2, x2, a3, b3, x3, mat, ustart, vstart ); \\ generate all the convergents of logp/logq that approximate it up to defaulprecision; /*** SET THE GLOBAL VARIABLE THETA, ONEOVERTHETA, PN, QN ********/ THETA = logp/logq; ONEOVERTHETA = logq/logp; cfrac = contfrac(THETA); mat = contfracpnqn(cfrac, #cfrac); \\ mat = matrix of all the convergents PN = mat[1,]; \\ numerators of the convergents QN = mat[2,]; \\ denominators of the convergents y1 = logq / LOGALPHA; \\ q is prime >= 3, alpha <=2, hence y1>0 (since q>=3>2>=alpha) cannot be an integer; so setsearch with three \\ parameters will give me the position of the minimal denominator > y1; so \\ s_{m+1} >= ceil(y1) > y1; in other words, index is M+1 in the paper index = setsearch(QN, ceil(y1)); \\ search if ceil(y1) is in the list if (index == 0, \\ ceil(y1) is not in the list \\ setsearch with three parameters returns the index where ceil(y1) should be inserted; so QN[index] is minimal and > ceil(y1) index = setsearch(QN, ceil(y1), 1); ); a2 = QN[index]; b2 = PN[index-1]; x2 = a2 * logp + b2 * logq; a3 = QN[index-1]; b3 = PN[index]; x3 = a3 * logp + b3 * logq; if( x2 < x3 , ustart = a2; vstart = b2, \\else ustart = a3; vstart = b3 ); \\ n_start is p^ustart * q^vstart exp_vec[1] = ustart; exp_vec[2] = vstart; exp_vec[3] = 0; return(exp_vec); } {msg_iterations(p,q, u,v,iterations,param, text) = local(printcontrol); printcontrol = iterations/param; if ( iterations == param, print("---")); if ( printcontrol == floor(printcontrol) , print("performed ", iterations, " steps in the ", text," search"); print("[",p,",",q,"]-exponents: u = ", u, ", v = ", v); print("---"); ); } {constants_comp(seqfile, p,q, Sums, TOT_ITER, non_negative_rho, LOOP, min_constant, MAX_constant, pmin, qmin, pmax, qmax, umin, vmin, umax, vmax, logp, logq ) = local(mu, D1, D2, C1, C2, C3, C4, prodlog, logNmin ,logNmax, loglogNmax, loglogp); \\mu = Sums/TOT_ITER; \\print(non_negative_rho); mu = Sums/non_negative_rho; D2 = min_constant; D1 = MAX_constant; C4 = 8*q; loglogp = log(logp); logNmax = umax * log(pmax) + vmax * log(qmax); loglogNmax = log(logNmax); logNmin = umin * log(pmin) + vmin * log(qmin); C2 = D2 + log(C4)/log(logNmin); \\if (loglogNmax <> loglogp, \\ C1 = D1 * loglogNmax/(loglogNmax-loglogp); \\ C3 = logp^C1, \\ \\else \\ C1 = "NA: Nmax = p"; \\ C3 = "NA: Nmax = p" \\ ); \\ corrected on Nov 30th, 2023, after Florian email C1 = D1 * loglogNmax/(loglogNmax-loglogp + log(2)); C3 = (0.5*logp)^C1; prodlog = logp*logq; output_data(p,q, TOT_ITER, non_negative_rho, LOOP, D1, D2, C1, C2, C3, C4, prodlog, Sums, mu); save_data(seqfile, p,q, TOT_ITER, non_negative_rho, LOOP, D1, D2, C1, C2, C3, C4, prodlog, Sums, mu); \\ BUILDING DATA for the global estimates if ( C1 > C1_GLOBAL, C1_GLOBAL = C1; UMAX_GLOBAL_C1 = umax; VMAX_GLOBAL_C1 = vmax; PMAX_GLOBAL_C1 = pmax; QMAX_GLOBAL_C1 = qmax ); if ( C2 < C2_GLOBAL, C2_GLOBAL = C2; UMIN_GLOBAL_C2 = umin; VMIN_GLOBAL_C2 = vmin; PMIN_GLOBAL_C2 = pmin; QMIN_GLOBAL_C2 = qmin ); if ( D1 > D1_GLOBAL, D1_GLOBAL = D1; UMAX_GLOBAL_D1 = umax; VMAX_GLOBAL_D1 = vmax; PMAX_GLOBAL_D1 = pmax; QMAX_GLOBAL_D1 = qmax ); if ( D2 < D2_GLOBAL, D2_GLOBAL = D2; UMIN_GLOBAL_D2 = umin; VMIN_GLOBAL_D2 = vmin; PMIN_GLOBAL_D2 = pmin; QMIN_GLOBAL_D2 = qmin ); } {output_data(p,q, TOT_ITER, non_negative_rho, LOOP, D1, D2, C1, C2, C3, C4, prodlog, Sums, mu) = local(); print ("We generated ", TOT_ITER , " left neighbors"); print ("We generated ", LOOP - TOT_ITER + EXTRA, " right neighbors"); print ("We generated ", non_negative_rho, " non-negative rho_i"); print ("D1(", p,",", q,") = ", D1 ); print ("D2(", p,",", q,") = ", D2 ); print ("C1(", p,",", q,") = ", C1 ); print ("C2(", p,",", q,") = ", C2 ); print ("C3(", p,",", q,") = ", C3 ); print ("C4(", p,",", q,") = ", C4 ); print ("Sum = ", Sums); print ("Averaged sum = ", mu); print ("logp*logq = ", prodlog); print("--------------"); } {final_output() = local(); print("-----------------------------"); print(" GLOBAL RESULTS "); print("-----------------------------"); print ("D1_GLOBAL = ", D1_GLOBAL, " attained at p = ", PMAX_GLOBAL_D1, ", q = ", QMAX_GLOBAL_D1); print ("with exponents: u = ", UMAX_GLOBAL_D1, ", v = ", VMAX_GLOBAL_D1); print("-----"); print ("D2_GLOBAL = ", D2_GLOBAL, " attained at p = ", PMIN_GLOBAL_D2, ", q = ", QMIN_GLOBAL_D2); print ("with exponents: u = ", UMIN_GLOBAL_D2, ", v = ", VMIN_GLOBAL_D2); print("-----"); print ("C1_GLOBAL = ", C1_GLOBAL, " attained at p = ", PMAX_GLOBAL_C1, ", q = ", QMAX_GLOBAL_C1); print ("with exponents: u = ", UMAX_GLOBAL_C1, ", v = ", VMAX_GLOBAL_C1); print("-----"); print ("C2_GLOBAL = ", C2_GLOBAL, " attained at p = ", PMIN_GLOBAL_C2, ", q = ", QMIN_GLOBAL_C2); print ("with exponents: u = ", UMIN_GLOBAL_C2, ", v = ", VMIN_GLOBAL_C2); print("-----------------------------"); print("-----------------------------"); } {save_data(seqfile, p,q, TOT_ITER, non_negative_rho, LOOP, D1, D2, C1, C2, C3, C4, prodlog, Sums, mu) = local(); filewrite1(seqfile, p); filewrite1(seqfile,";"); filewrite1(seqfile, q); filewrite1(seqfile,";"); filewrite1(seqfile, TOT_ITER); filewrite1(seqfile,";"); filewrite1(seqfile, non_negative_rho); filewrite1(seqfile,";"); filewrite1(seqfile, LOOP); filewrite1(seqfile,";"); filewrite1(seqfile, D1); filewrite1(seqfile,";"); filewrite1(seqfile, D2); filewrite1(seqfile,";"); filewrite1(seqfile, C1); filewrite1(seqfile,";"); filewrite1(seqfile, C2); filewrite1(seqfile,";"); filewrite1(seqfile, C3); filewrite1(seqfile,";"); filewrite1(seqfile, C4); filewrite1(seqfile,";"); filewrite1(seqfile, prodlog); filewrite1(seqfile,";"); filewrite1(seqfile, Sums); filewrite1(seqfile,";"); filewrite(seqfile, mu); fileflush(seqfile); } /****************** ----------------- OPTIPLEX ----------------- languasc@languasco1:~/Desktop/EK/Moree-Tijdeman80$ gp2c-run -pmy_ -g -W constants_pq_left_right.gp Warning:constants_pq_left_right.gp:43: variable undeclared f GP/PARI CALCULATOR Version 2.15.4 (released) amd64 running linux (x86-64/GMP-6.2.1 kernel) 64-bit version compiled: Jul 10 2023, gcc version 11.3.0 (Ubuntu 11.3.0-1ubuntu1~22.04.1) threading engine: pthread (readline v8.1 enabled, extended help enabled) Copyright (C) 2000-2022 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 = 500000, nbthreads = 4 --------- constants_pq(1000, 10^4) --------- gp2c-run -pmy_ -g -W constants_pq_left_right.gp report-constants-1000.txt ******************/