/* 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(ITER_TRIVIAL_SECOND, ITER_TRIVIAL_FIRST, ITER_CONTFRAC, PN, QN); global(THETA, ONEOVERTHETA, LOGALPHA, LOGONEOVERALPHA); /***************************************** MAIN FUNCTION; the one to be called ****************************************/ {npq_alpha(boundprime, alpha = 2) = 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, lastelement, timebound, partelaptime, partialtime, exp_vec, found, logten, digits_n_start, digits_n_end, oneminute, prev_elapsed ); 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("---------"); if (alpha <= 1, error("alpha must be greater than 1")); if (alpha > 2, error("alpha must be less or equal to 2")); print("Searching for n_{p,q} (alpha); alpha = ", alpha ); print("with the continued fractions algorithm; left neighbor" ); print("---------"); name="npq_data"; namefile = Strprintf("%s-%s-(%1.3f).csv", name, boundprime, alpha); seqfile = fileopen(namefile, "w"); string = "p;q;alpha;u_start;v_start;u_end;v_end;digits_n_start; digits_n_end; ITER_CONTFRAC;ITER_TRIVIAL_FIRST;ITER_TRIVIAL_SECOND;iter_total;npqalpha"; filewrite(seqfile, string); /*** SET THE GLOBAL VARIABLE LOGALPHA, LOGONEOVERALPHA ********/ LOGALPHA = log(alpha); LOGONEOVERALPHA = - LOGALPHA; logten = log(10); 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; logq = log(q); \\ 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("--------------"); 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 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) while ( ok == 0, \\ 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_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]; if (found == 1 && uleft <> u && vleft <> v, \\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 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(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(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 u = exp_vec[1]; v = exp_vec[2]; 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 u = exp_vec[1]; v = exp_vec[2]; ok = 0; msg_iterations(p,q,u,v,ITER_TRIVIAL_FIRST,1000, "first trivial"); ); ); ); \\ lastelement is n0 of the paper lastelement = p^uend * q^vend ; digits_n_end = floor(uend *logp/logten + vend*logq/logten)+1; npqalpha = ceil(lastelement / alpha); output_data(p,q,npqalpha, alpha, lastelement, uend, vend, digits_n_end); print_data(seqfile, p,q, alpha, ustart, vstart, uend, vend, digits_n_start, digits_n_end, npqalpha); ); ); print("*************************************"); fileflush(seqfile); fileclose(seqfile); time = getwalltime() - starttime; elaptime = time; seconds = floor(elaptime /1000)%60; minutes = floor(elaptime /oneminute); millisec = elaptime - minutes*oneminute - seconds*1000; print("npqalpha 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(exp_vec, ubad, vbad, logp, logq, p, q, alpha) = local(u1, v1, logn, lbound, ubound, z, y, astart, bstart, aend, bend); \\ 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) ); 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 , \\print("*** second trivial search found one element !!!"); exp_vec[1] = a; exp_vec[2] = 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 , \\print("*** second trivial search found one element !!!"); exp_vec[1] = a; exp_vec[2] = b; exp_vec[3] = 1; ITER_TRIVIAL_SECOND+=1; return(exp_vec); ); ); ); \\endfor a );\\endfor b return(exp_vec); } /***************************************** FIRST TRIVIAL SEARCH (using also data from failed iteration of the continued fraction search) ****************************************/ {first_trivial_search(exp_vec, ubad, vbad, logp, logq, p, q) = local(u1,v1, logn, lbound, ubound, z, y, astart, bstart, aend, bend); \\ 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; 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 , \\print("*** first trivial search found one element !!!"); exp_vec[1] = a; exp_vec[2] = b; exp_vec[3] = 1; ITER_TRIVIAL_FIRST+=1; return(exp_vec); ); ); ); \\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 , \\print("*** first trivial search found one element !!!"); exp_vec[1] = a; exp_vec[2] = b; exp_vec[3] = 1; ITER_TRIVIAL_FIRST+=1; return(exp_vec); ); ); ); \\endfor a );\\endfor b return(exp_vec); } /***************************************** CONTINUED FRACTION SEARCH [uses the new algortithm] ****************************************/ {cont_frac_search(exp_vec, logp, logq) = local( u_local, v_local, u1, v1, u2, v2, index, x1, x2, xmax, xmin, usummand_min, vsummand_min, usummand_max, vsummand_max ); \\ 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 at least one between xmin, xmax is ok if (xmin >= LOGONEOVERALPHA, \\choosing xmin exp_vec[1] = u_local + usummand_min; \\ uleft exp_vec[2] = v_local + vsummand_min; \\ vleft exp_vec[3] = 1; \\ found a good point return(exp_vec), \\ else [xmin< -log(alpha)] \\choosing 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("---"); ); } {output_data(p,q,npqalpha, alpha, lastelement, uend, vend, digits_n_end) = local(); print("Searches terminated with [",p,",",q,"]-exponents: u = ", uend, ", v = ", vend); print("--------------"); print("number of left neigbors generated (cont_frac) = ", ITER_CONTFRAC); print("number of left neigbors generated (trivial first) = ", ITER_TRIVIAL_FIRST); print("number of left neigbors generated (trivial second) = ", ITER_TRIVIAL_SECOND); print("TOTAL number of left neigbors generated = ", ITER_CONTFRAC + ITER_TRIVIAL_FIRST + ITER_TRIVIAL_SECOND); print("--------------"); print("Last generated element of the sequence = ", lastelement); print ("with [",p,",",q,"]-exponents: u = " , uend,"; v = ", vend); print ("n_end has ", digits_n_end , " digits"); print("--------------"); print ("n_{", p,",", q,"}(", alpha, ") = ", npqalpha); print("--------------"); } {print_data(seqfile, p,q, alpha, ustart, vstart, uend, vend, digits_n_start, digits_n_end, npqalpha) = local(); filewrite1(seqfile, p); filewrite1(seqfile,";"); filewrite1(seqfile, q); filewrite1(seqfile,";"); filewrite1(seqfile, alpha); filewrite1(seqfile,";"); filewrite1(seqfile, ustart); filewrite1(seqfile,";"); filewrite1(seqfile, vstart); filewrite1(seqfile,";"); filewrite1(seqfile, uend); filewrite1(seqfile,";"); filewrite1(seqfile, vend); filewrite1(seqfile,";"); filewrite1(seqfile, digits_n_start); filewrite1(seqfile,";"); filewrite1(seqfile, digits_n_end); filewrite1(seqfile,";"); filewrite1(seqfile, ITER_CONTFRAC); filewrite1(seqfile,";"); filewrite1(seqfile, ITER_TRIVIAL_FIRST); filewrite1(seqfile,";"); filewrite1(seqfile, ITER_TRIVIAL_SECOND); filewrite1(seqfile,";"); filewrite1(seqfile, ITER_CONTFRAC + ITER_TRIVIAL_FIRST + ITER_TRIVIAL_SECOND); filewrite1(seqfile,";"); filewrite(seqfile, npqalpha); } /****************** Last login: Fri Aug 11 22:11:38 from fe80::1452:dc77:8b14:3e83%en0 languasc@languasc-pro ~ % cd /Users/languasc/Documents/LANGUASCO/matematica/lavori/\[74\]Moree-Tijdeman80/programs/npqalpha/source_codes_npq_alpha languasc@languasc-pro source_codes_npq_alpha % gp2c-run -pmy_ -g -W npq_alpha_warp9.gp Reading GPRC: /Users/languasc/.gprc GPRC Done. GP/PARI CALCULATOR Version 2.15.4 (released) arm64 running darwin (aarch64/GMP-6.2.1 kernel) 64-bit version compiled: Jul 10 2023, Apple clang version 14.0.3 (clang-1403.0.22.14.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. parisizemax = 2048000000, primelimit = 500000, nbthreads = 8 ? npq_alpha(100,2) ................... .................. -------------- ************************************* p = 89; q = 97; n = p^u*q^v; u,v >= 0 -------------- Searches starting with [89,97]-exponents: u = 1, v = 52 n_start has 106 digits Searches terminated with [89,97]-exponents: u = 45, v = 0 -------------- number of left neigbors generated (cont_frac) = 203 number of left neigbors generated (trivial first) = 35 number of left neigbors generated (trivial second) = 0 TOTAL number of left neigbors generated = 238 -------------- Last generated element of the sequence = 5278983432532964757288775688474865115065412211545081783600102880716531908000481333645049 with [89,97]-exponents: u = 45; v = 0 n_end has 88 digits -------------- n_{89,97}(2) = 2639491716266482378644387844237432557532706105772540891800051440358265954000240666822525 -------------- ************************************* npqalpha computation time (output time included): 0 min, 40 sec, 462 millisec ****** END PROGRAM ******** ? ******************/