semigroup.gp

print("\n\nType '?semigroup' for help.\n\n");
parigp_version = version();
semigroup_library = strprintf("./libsemigroup-%d-%d-%d.so", parigp_version[1], parigp_version[2], parigp_version[3]);

/*semigroup.c*/
  addhelp(semigroup,"Basic (semi)group methods:\n\tLRword, semigroup_growth, semigroup_mats, semigroup_mgens.\n\nSemigroup orbits:\n\tsemigroup_missing, semigroup_missing_parabolic, semigroup_missinglist, semigroup_orbit.\n\nPsi methods:\n\tkron_all, kron_seq_both, psi_mats, psi_rep.\n\nContinued fractions:\n\tcontfraceven, contfractoQ, contfractoword, contfrac_tail_missing.\n\nLinear regression:\n\tOLS, OLS_nointercept, OLS_single, rsquared.\n\nPaper methods:\n\nTesting:\n\trunalltests, test_evencontfrac, test_gamma14geq0_kronequal, test_kronaction, test_kronaction_many, test_psioogens, test_psisemigroup, test_table1_psi, test_table1_psi_many.\n\nSupplementary computations:\n\ttest_psi_23orbit, test_psi1_23orbit, test_psi1_hdim, test_table1_psi1, test_table1_psi1_many.\n\nSupporting methods:\n\tgamma14geq_random, psi_random, sl2zgeq0_random, oddpart, psi_isreciprocity, psi_missingsquares, psi1_missingsquares, table1_line.");
  
/*SECTION 1: BASIC (SEMI)GROUP METHODS*/
  install(LRword,"G",,semigroup_library);
  addhelp(LRword,"LRword(M): for M a SL(2, Z)^{>=0} hyperbolic matrix with positive entries, returns M as a word in L and R, where L=[1,1;0,1] and R=[1,0;1,1]. This is a Vecsmall with 0 representing L and 1 representing R. M must itself have entries that fit into the sice of a C long.");
  install(semigroup_growth,"GD100,L,D500,L,DGp");
  addhelp(semigroup_growth,"semigroup_growth(mats, {binsize=100}, {Nbins=500}, {start=[1, 1]}): estimates the growth rate of the orbits of the semigroup generated by mats, whose elements all have infinite order and nonnegative entries. We assume the growth rate is of the form c*N^nu, and we estimate this by computing the orbit of the column vector in [binsize*Nbins, 2*binsize*Nbins], putting the sizes of the vectors into bins, and then running a linear regression. The return value is [c, nu, R^2], where R^2 is the R^2 value of the regression. The larger this value, the more confidence we have in (c, nu).");
  install(semigroup_mats,"GL");
  addhelp(semigroup_mats,"semigroup_mats(mats, N): assume the matrices in mats all have positive entries and infinite order. This returns the matrices in the semigroup they generate that have all entries at most N. If there are relations, the corresponding matrices will get counted multiple times.");
  install(semigroup_mgens,"G");
  addhelp(semigroup_mgens,"semigroup_mgens(mats): returns a minimal set of generators for the semigroup generated by mats. We assume that these matrices have all nonnegative entries, determinant +/-1, and remove the identity element if it occurs.");

/*SECTION 2: SEMIGROUP ORBITS*/
  install(semigroup_missing,"GGGGLD1,L,");
  addhelp(semigroup_missing,"semigroup_missing(mats, B, start, congs, entry, {load=1}): prints to a file missing entries in a semigroup orbit, and al loads them if load=1. mats is the vector of matrices that generate the semigroup, necessarily all containing nonnegative entries. B is the bounds we wish to search, either a positive integer (for 1 to B), or a vector of two integers [B1, B2] for between B1 and B2. Note that we stop once any entry in the orbit is > B, which means that for some choices of matrices, some entries near B are declared missing but are not actually, as the first time they appear is with a larger entry. start is the starting vector for the orbit, which should be primitive and nonnegative. congs=[[r1, r2, ..., rk], modulus] is supplied as the congruence restictions: the entryth entry of a vector is the element we track, and the user must supply the possible congruence classes it can fall into, i.e. x mod(modulus) is one of r1, r2, ..., rk. If there are no congruence obstructions, you should supply [[0], 1].");
  install(semigroup_missing_parabolic,"GGGGLD1,L,");
  addhelp(semigroup_missing_parabolic,"semigroup_missing_parabolic(mats, B, start, congs, entry, {load=1}): same as semigroup_missing, except up to 80% slower but uses less memory. If one of the matrices is parabolic, the memory savings is extreme, and it is suggested to use this method when B gets large (10^9 or so). This method will work when none of the matrices are parabolic, though it is unlikely that the difference in memory will matter, and using semigroup_missing is suggested.");
  install(semigroup_missinglist,"iGGGL");
  addhelp(semigroup_missinglist,"semigroup_missinglist(mats, miss, start, entry): this returns 1 if the sorted vector of positive integers miss is entirely missed by the semigroup orbit of <mats>*start, where we look at the entryth entry of the orbit vectors. Calling this on a vector of square numbers gives evidence towards this being missing from an orbit. This is not optimized for parabolic generators, so should not be used for vectors miss which go past 10^9.");
  install(semigroup_orbit,"GLG");
  addhelp(semigroup_orbit,"semigroup_orbit(mats, B, start): for the semigroup given by mats (with all nonnegative entries and infinite order), this returns the orbit of the semigroup on the (column) vector start, where the entries are bounded by B. The output can have repeats, especially if there are relations among the matrices in mats. This is mainly useful for some initial testing, as the other semigroup_missing methods are significantly more efficient when searching for missing entries.");

/*SECTION 3: PSI METHODS*/
  install(kron_all,"U");
  addhelp(kron_all,"kron_all(n): computes the sequence of kronecker values kron(x, n), returning a Vecsmall v where v[i]=kron(i-1, n). The length of the sequence, M, is the smallest guaranteed period. In other words, M is the smallest positive integer such that n/M is a rational square, M has same prime factors as n, and if n is even then 4|M (thus if n=2*odd prime, then M=4n. If n is not 2 mod 4, then M<=n).");
  install(kron_seq_both,"lU");
  addhelp(kron_seq_both,"kron_seq_both(n): returns the smallest m such that in any sequence of m consecutive integers, there is at least one for which kronecker(r+, n)=1 and kronecker(r-, n)=-1. If n is a square, returns 0. Should not be used on n larger than 10^9 or so, as it requires memory proportional to n.");
  install(psi_mats,"L");
  addhelp(psi_mats,"psi_mats(N): Returns all the matrices in Psi with maximum entry N (Psi as defined in the paper: elements of Gamma_1(4) with nonnegative entries and Kronecker symbol of a row/column 1.");
  install(psi_rep,"iLLLD1,L,");
  addhelp(psi_rep,"psi_rep(x, y, n, {entry=1}): returns 1 if n is a numerator (entry=1) or denominator (entry=2) in the orbit Psi*[x,y]~. This is based on Lemma 4.1. All inputs must fit as C long's on your operating system, and (x, y) need to be nonnegative and coprime.");

/*SECTION 4: CONTINUED FRACTIONS*/
  install(contfraceven,"G");
  addhelp(contfraceven,"contfraceven(q): returns the (unique) even length continued fraction representation of the rational number q.");
  install(contfractoQ,"G");
  addhelp(contfractoQ,"contfractoQ(v): returns the rational number corresponding to the continued fraction v.");
  install(contfractoword,"G");
  addhelp(contfractoword,"contfractoword(v): given an even length continued fraction [a0;a1, a2, ..., a(2n-1)] with a0>=0, this returns the corresponding word W in SL(2, Z), i.e. W=L^a0*R^a1*...*L^a(2n-2)*R^a(2n-1), where L=[1,1;0,1] and R=[1,0;1,1]. In particular, W[1,0]~=[x,y], where contfraceven(x/y)=v. If v has odd length, we convert it to having even length.");
  install(contfrac_tail_missing,"GD1,L,");
  addhelp(contfrac_tail_missing,"contfrac_tail_missing(bounds, load): searches for the continued fractions of the form [4a_0;4a_1,...,4a_n,a_{n+1},1,2] with denoninator to to bounds (or between B1 and B2, if bounds=[B1, B2]), saving to a file those that are missing, and loading them if load=1.");

/*SECTION 5: LINEAR REGRESSION*/
  install(OLS,"GGD1,L,");
  addhelp(OLS,"OLS(X, y, {retrsqr=1}): performs ordinary least squares regression on the data, with the inputs being the n columns of the matrix X, and the outputs being the entries of the column vector y. We include a constant term, so the first row of X must be all 1's. If retrsqr=1, returns [params, R^2], and otherwise returns params, where params is the column vector of best fit parameters.");
  install(OLS_nointercept,"GGD1,L,");
  addhelp(OLS_nointercept,"OLS_nointercept(x, y, {retrsqr=1}): performs ordinary least squares regression assuming that y[i]=c*X[i] (where X and y are both (column) vectors), i.e. the y-intercept is 0. Returns c if retrsqr=0, or [c, R^2] otherwise.");
  install(OLS_single,"GGD1,L,");
  addhelp(OLS_single,"OLS_single(x, y, {retrsqr=1}): performs linear regression for a single variable (essentially a macro for OLS with y=mx+b). Requires y to be a column vector and x is either a vector or column vector.");
  install(rsquared,"GGG");
  addhelp(rsquared,"rsquared(X, y, fit): returns the R^2 value for the proposed linear regression, where the input is X, output is y, and fit is the proposed parameters.");

/*paper.gp*/

/*SECTION 1: TESTING*/
  addhelp(runalltests,"runalltests(): runs basic tests for all of the testing and supplementary methods.");
  addhelp(test_evencontfrac,"test_evencontfrac(n, B): tests that the even continued fraction does correspond to orbits of [1,0]~ as described at the start of section 2. We pick random positive rational numbers x/y with 1<=x,y<=B, and do this test n times.");
  addhelp(test_gamma14geq0_kronequal,"test_gamma14geq0_kronequal(n, B): generates n elements of Gamma_1(4)^{>=0} whose coefficients are at most B, and checks Lemma 3.1: the Kronecker symbols of their rows and columns are all equal.");
  addhelp(test_kronaction,"test_kronaction(M, xy): tests Proposition 3.2: M=[a, b;c, d] in SL(2, Z)^{>=0} a, b, c, d >= 0, xy=[x, y] with x, y>=0 coprime, gcd(x, d)=1, this proposition gives a formula for kron(ax+by/cx+dy). This function returns 1 if and only if the formula is correct for the given inputs. If the method returns anything other than 1, then the formula has failed.");
  addhelp(test_kronaction_many,"test_kronaction_many(n, B, xymin, xymax): tests Proposition 3.2 by calling test_kronaction on n random matrices in SL(2, Z)^{>=0} with entries bounded by B. The values of x, y tried are all valid pairs with xymin<=x, y<=xymax. If the formula fails, we print the failing inputs, and raise an error. If no error occurs, all tests passed successfully.");
  addhelp(test_psioogens,"test_psioogens(n): tests that the matrices M_k=[12k+1,4k;12k+4,4k+1] are all in a minimal generating set for Psi up to k=n.");
  addhelp(test_psisemigroup,"test_psisemigroup(n, B): tests that Psi is a semigroup by generating n pairs of elements in Psi with coefficients bounded by B, and checking that their product lies in Psi.");
  addhelp(test_table1_psi,"test_table1_psi(xy, B, entry): returns 1 if there was a reciprocity obstruction for Psi*xy and we did not find any squares up to B^2, or if the squares are ruled out by a congruence obstruction, or if squares are not ruled out and there is no reciprocity obstruction and we find squares. If there is a reciprocity obstruction and squares are found, raises an error. If there is no reciprocity obstruction, squares are not ruled out by congruence, and we find no squares, returns 0 (presumably B is not large enough).");
  addhelp(test_table1_psi_many,"test_table1_psi_many(xymax = 200, B = 300): tests that the reciprocity obstructions listed in Table 1 are both correct and complete for Psi, for all 1<=x, y<=xymax, and computing squares up to B^2. If squares are found when not allowed, raises an error. If no squares are found when they should be, prints the inputs to the terminal. If B is too small, then this may happen, and you can try to increase B to check that eventually squares do appear. Returns the vector of counts of how many pairs of each line of Table 1 (1 to 9) were tested.");

/*SECTION 2: SUPPLEMENTARY COMPUTATIONS*/
  addhelp(test_psi_23orbit,"test_psi_23orbit(n = 109120000): tests the non-square integers between 1 and n, printing any which are not numerators of Psi*[2, 3]~. Based on the proof of Theorem 2.7, any non-square i>=default value of n is a numerator in the orbit, hence this gives the complete list.");
  addhelp(test_psi1_23orbit,"test_psi_23orbit(n = 10^7): finds all missing numerators in Psi_1*[2, 3]~ up to n, printing any non-squares.");
  addhelp(test_psi1_hdim,"test_psi1_hdim(): estimates the Hausdorff dimension of Psi_1 using a least squares regression to estimate the critical exponent (the method semigroup_growth).");
  addhelp(test_table1_psi1,"test_table1_psi1(xy, B, entry): returns 1 if there was a reciprocity obstruction for Psi_1*xy and we did not find any squares up to B^2, or if the squares are ruled out by a congruence obstruction, or if squares are not ruled out and there is no reciprocity obstruction and we find squares. If there is a reciprocity obstruction and squares are found, raises an error. If there is no reciprocity obstruction, squares are not ruled out by congruence, and we find no squares, returns 0 (presumably B is not large enough).");
  addhelp(test_table1_psi1_many,"test_table1_psi1_many(xymax = 200, B = 300): tests that the reciprocity obstructions listed in Table 1 are both correct and complete for Psi_1, for all 1<=x, y<=xymax, and computing squares up to B^2. If squares are found when not allowed, raises an error. If no squares are found when they should be, prints the inputs to the terminal. If B is too small, then this may happen, and you can try to increase B to check that eventually squares do appear. Returns the vector of counts of how many pairs of each line of Table 1 (1 to 9) were tested.");

/*SECTION 3: SUPPORTING METHODS*/
  addhelp(gamma14geq0_random,"gamma14geq0_random(n): returns a random element of Gamma_1(4)^{>=0} where the size of the coefficients are bounded by n (not necessarily uniformly, but should be reasonably close).");
  addhelp(psi_random,"psi_random(n): returns a random element of Psi where the size of the coefficients are bounded by n (not necessarily uniformly, but should be reasonably close).");
  addhelp(sl2zgeq0_random,"sl2zgeq0_random(n): returns a random element of SL(2, Z)^{>=0} where the size of the coefficients are bounded by n (not necessarily uniformly, but should be reasonably close).");
  addhelp(oddpart,"oddpart(n): returns the odd part of the integer n.");
  addhelp(psi_isreciprocity,"psi_isreciprocity(xy, entry: returns 1 if Table 1 predicts there to be no squares in Psi*xy (due to a reciprocity obstruction), 0 if squares are possible and there is no reciprocity obstruction listed, and -1 if squares are not possible due to a congruence obstruction. entry = 1 corresponds to the numerator, 2 the denominator.");
  addhelp(psi_missingsquares,"psi_missingsquares(xy, B, entry): tests if all the squares in the orbit of Psi*xy are missing up to B^2. Returns 1 if they are missing, 0 else.");
  addhelp(psi1_missingsquares,"psi1_missingsquares(xy, B, entry): tests if all the squares in the orbit of Psi_1*xy are missing up to B^2. Returns 1 if they are missing, 0 else.");
  addhelp(table1_line,"table1_line(xy): given the pair [x, y]~, returns the line (1 - 9) in Table 1 that the pair corresponds to (i.e. what the congruence and reciprocity obstructions are).");

read("paper.gp");
default(parisize, "4096M");/*Must come at the end*/