# John Miller # Final project - Math 640 - Rutgers - Experimental Mathematics - Spring 2012 # This file contains Maple procedures that were used in the developments # described in the related paper: # "C-finite sequences mod m, toral autmorphisms, and Euclideanity of Q(sqrt(2+sqrt(2))" # *** Run Help() to get a summary of the procedures contained herein. *** Help := proc(): print(`Please find below a summary of the Maple procedures that were used in the development of the related paper, "C-finite sequences mod m, toral autmorphisms, and Euclideanity of Q(sqrt(2+sqrt(2))".`): print(`*********************************************************************`): Help_Torus(): print(`*********************************************************************`): Help_Cfinite_mod_m(): print(`*********************************************************************`): print(`Run "HelpOld()" for a list of older procedures of tertiary interest that are primarily related specifically to the study of the Euclideanity of Q(sqrt(2+sqrt(2)). These procedures are less well documented and maintained, but are included for completeness.`): end: ######################################################################### Help_Torus := proc(): print(`Help_Torus(): Below are utilities to calculate the intersection of the image of polygons under an automorpishm of the torus. Polygons are stored as lists of vertices.`): print(``): print(`The procedures of primary interest are:`): print(`RecurMatrix(L): Inputs a 2nd order C-finite recurrence. Outputs a 2x2 companion matrix.`): print(`Box(Left,Right,Bottom,Top): Inputs 4 rational numbers and outputs a rectangle.`): print(`SelfIntersect(T,P,n): Inputs a toral automorphism T (2x2 matrix) and a polygon P, and outputs a set of polygons which is the intersection (mod 1) of (T^k)(P), k=0,1,...,n.`): print(`PlotPolygons(Polygons): Inputs a set of polygons, and plots them.`): print(``): print(`Procedures of secondary interest are:`): print(`TransformPolygonR2(T,P), BreakPolygon(P), MovePolygonToTorus(P), PolygonMod1(P), TransformMod1(T,P,n), IntersectPolygon(P,Q), InsideClip(E,clipEdge), CalcIntersect(P1,P2,P3,P4), TransformManyMod1(T,S,n), IntersectManyPolygons(S1,S2)`); end: ######################################################################### Help_Cfinite_mod_m := proc(): print(`Help_Cfinite_mod_m(): Below are procedures that were used to study the distribution of C-finite recurrences modulo m.`): print(``): print(`The procedure of primary interest is:`): print(`DistRecur(L,M): Inputs a C-finite recurrence L, and outputs the "minimum distance spectrum" up to maximum modulus M.`): print(``): print(`Procedures of secondary interest are:`): print(`DistRecurMod(L,m), MinMaxRecurMod(L,m), MinMaxSeqMod(C,m), RecToSeqMod(C,n,m), RecToSeq(C,n), AllWords(A,n)`): end: # DistRecur(L,M): Inputs list L, giving the recurrence relation for a # C-finite sequence, and a maximum modulus M, and outputs # the set of all "distances" from modulus 1 to M, representing the # subset of the "minimum distance spectrum" up to modulo M. ######################################################################### # Let T^2 be the 2-torus given by R^2/Z^2. # We would like to define a polygon in T^2 by a list of coordinates # [x_i,y_i], but on the torus it is ambiguous as to which direction # represents the interior of the polygon. So we orient our list # of points by writing it in the counterclockwise direction going # around the interior of the polygon, and the last element of this # list should equal the first element. # The toral automorphism is given by a 2x2 matrix. with(LinearAlgebra): with(ListTools): with(plots): with(geometry): interface(rtablesize=50): ######################################################################### # RecurMatrix(L): Inputs a list L indicating the second order C-finite # recurrence, a[n] + L[1]*a[n-1] + L[2]*a[n-2] = 0, # i.e., a[n] = -L[1]*a[n-1] + -L[2]*a[n-2], # and outputs the 2x2 matrix T = | 0 1 | # | -L[2] -L[1] | # so that T( [ a[n-2],a[n-1] ] ) = [ a[n-1],a[n] ] # Note that this is a toral automorphism iff |det(T)| = 1. RecurMatrix := proc(L::list): if nops(L) <> 2 then return(FAIL): fi: Matrix([ [0,1],[-L[2],-L[1]] ]): end: # end of RecurMatrix(L) ######################################################################### # Box(Left,Right,Bottom,Top): Inputs 4 numbers and returns # the corresponding "polygon" in the form of a list of 5 pairs of # numbers [x,y] in R^2, oriented counterclockwise. Box := proc(Left::rational,Right::rational,Bottom::rational,Top::rational): if Left > Right or Bottom > Top then return(FAIL): fi: [ [Left,Bottom], [Right,Bottom], [Right,Top], [Left,Top], [Left,Bottom] ]: end: # end of Box(Left,Right,Bottom,Top) ######################################################################### # TransformPolygonR2(T,P): Inputs a 2x2 matrix T, and a list of # column vectors in R^2, representing a polygon P, and applies the # transformation given by T to it. The polygon's list should be # oriented in the counterclockwise direction (i.e. right-hand rule # point to the interior). TransformPolygonR2 := proc(T::Matrix, P::list) local pts::list, P2::list: P2 := [seq(convert(MatrixVectorMultiply(T,Vector(2,pts)),list), pts in P)]: if Determinant(T) < 0 then P2 := Reverse(P2): fi: P2: end: # end of TransformPolygonR2(T,P) ######################################################################### # BreakPolygon(P): Inputs a polygon P, and outputs a set of polygons # after breaking up P mod 1. BreakPolygon := proc(P::list) local minX, minY, maxX, maxY, i, j, pt, Q::set, P1:list: minX := floor(min(seq(pt[1],pt in P))): minY := floor(min(seq(pt[2],pt in P))): maxX := ceil(max(seq(pt[1],pt in P))): maxY := ceil(max(seq(pt[2],pt in P))): Q := {}: for i from minX to maxX-1 do: for j from minY to maxY-1 do: Q := Q union {IntersectPolygon(P,Box(i,i+1,j,j+1))}: od: od: Q: end: # end of BreakPolygon(P) ######################################################################### # MovePolygonToTorus(P): Inputs a polygon P, and moves the lower # left extremum to be in the box [0,1]x[0,1]. MovePolygonToTorus := proc(P::list) local minX, minY, pt: minX := floor(min(seq(pt[1],pt in P))): minY := floor(min(seq(pt[2],pt in P))): [seq(pt-[minX,minY], pt in P)]: end: # end of MovePolygonToTorus(P): ######################################################################### # PolygonMod1(P): Inputs a polygon, and returns a set of polygons # representing the polygon "mod 1" on the torus. PolygonMod1 := proc(P::list) local S, P1: S := BreakPolygon(P): {seq(MovePolygonToTorus(P1),P1 in S)} end: # end of PolygonMod1(P) ######################################################################### # TransformMod1(T,P,n): Inputs a 2x2 matrix T, and a list of # column vectors in the torus, representing a polygon P, and applies the # transformation given by T^n to it. The polygon's list should be # oriented in the counterclockwise direction (i.e. right-hand rule # point to the interior). TransformMod1 := proc(T::Matrix, P::list, n::nonnegint) local pts::list, P2::list, i: if n=0 then return(PolygonMod1(P)): fi: P2 := P: for i from 1 to n do: P2 := TransformPolygonR2(T,P2): od: PolygonMod1(P2): end: # end of TransformMod1(T,P) ######################################################################### # TransformManyMod1(T,S,n): Inputs a 2x2 matrix T, and a list of # column vectors in the torus, representing a polygon P, and applies the # transformation given by T^n to it. The polygon's list should be # oriented in the counterclockwise direction (i.e. right-hand rule # point to the interior). TransformManyMod1 := proc(T::Matrix, S::set, n::nonnegint) local pts::list, S2::set, i, P::list, P2::list: S2 := {}: for P in S do: P2 := P: for i from 1 to n do: P2 := TransformPolygonR2(T,P2): od: S2 := S2 union PolygonMod1(P2): od: S2: end: # end of TransformManyMod1(T,S,n) ######################################################################### # IntersectManyPolygons(S1,S2): Inputs sets of polygons S1, S2 and # outputs a set of intersected polygons. IntersectManyPolygons := proc(S1::set, S2::set) local P1::list, P2::list, S3::set: S3 := {}: for P1 in S1 do: for P2 in S2 do: S3 := S3 union {IntersectPolygon(P1,P2)}: od: od: S3: end: # end of IntersectManyPolygons(S1,S2) ######################################################################### # SelfIntersect(T,P,n): Inputs a 2x2 matrix T, and a list of # column vectors in the torus, representing a polygon P, and outputs # the intersection Intersect((T^k)(P),k=0..n): # The polygon's list should be oriented in the counterclockwise direction # (i.e. right-hand rule point to the interior). SelfIntersect := proc(T::Matrix, P::list, n::nonnegint) local k, S::set: S := {P}: if n = 0 then return(S): fi: for k from 1 to n do: S := IntersectManyPolygons(S,TransformManyMod1(T,S,1)): od: S: end: # end of SelfIntersect(T,P,n) ######################################################################### # PlotPolygons(Polygons): Inputs a set of polygons, and plots them. PlotPolygons := proc(Polygons::set) local minX, minY, maxX, maxY, minC, maxC, pt, P::list: minX := floor(min(seq(seq(pt[1],pt in P),P in Polygons))): minY := floor(min(seq(seq(pt[2],pt in P),P in Polygons))): maxX := ceil(max(seq(seq(pt[1],pt in P),P in Polygons))): maxY := ceil(max(seq(seq(pt[2],pt in P),P in Polygons))): minC := min(0,minX,minY): maxC := max(maxX,maxY): display(seq(plot(P),P in Polygons),view=[minC..maxC,minC..maxC]): end: # end of PlotPolygon(Polygons) ######################################################################### # IntersectPolygon(P,Q): Intersects two convex polygons P and Q using # the Sutherland-Hodgman algorithm. IntersectPolygon := proc(P::list,Q::list) local inputList::list, outputList::list, i, S::list, E::list, clipEdge::list, intersectPt: if nops(convert(P,set)) = 1 then return(NULL): fi: if nops(convert(Q,set)) = 1 then return(NULL): fi: outputList := P[1..nops(P)-1]: for i from 1 to nops(Q)-1 do: clipEdge := Q[i..i+1]: inputList := outputList: outputList := []: S := inputList[nops(inputList)]: for E in inputList do: if InsideClip(E,clipEdge) then if not InsideClip(S,clipEdge) then intersectPt := CalcIntersect(S,E,clipEdge[1],clipEdge[2]): outputList := [op(outputList),intersectPt]: fi: outputList := [op(outputList),E]: elif InsideClip(S,clipEdge) then intersectPt := CalcIntersect(S,E,clipEdge[1],clipEdge[2]): outputList := [op(outputList),intersectPt]: fi: S := E: od: #print(outputList): if nops(outputList)=0 then return(NULL): fi: od: if nops(outputList)=0 then return(NULL): fi: [op(MakeUnique(outputList)),outputList[1]]: end: # end of IntersectPolygon(P,Q) ######################################################################### # InsideClip(E,clipEdge): Inputs a point E (given as a list of 2 numbers) # and and a list of 2 points clipEdge, and returns "true" if E is to # the left of clipEdge, and "false" if not. InsideClip := proc(E::list,clipEdge::list) local v::list, w::list: v := clipEdge[2] - clipEdge[1]: w := E - clipEdge[1]: evalb(v[1]*w[2]-w[1]*v[2] >= 0): end: # end of InsideClip(E,clipEdge) ######################################################################### # CalcIntersect(P1,P2,P3,P4): Calculates the point of intersection # of the liens given by the line segments P1P2 and P3P4. # "Points" here are lists of two elements. CalcIntersect := proc(P1::list,P2::list,P3::list,P4::list) local x,y,x1,x2,x3,x4,y1,y2,y3,y4, detL1,detL2,detL1x,detL1y,detL2x,detL2y,Denom: x1 := P1[1]: y1 := P1[2]: x2 := P2[1]: y2 := P2[2]: x3 := P3[1]: y3 := P3[2]: x4 := P4[1]: y4 := P4[2]: detL1 := x1*y2 - x2*y1: detL2 := x3*y4 - x4*y3: detL1x := x1-x2: detL1y := y1-y2: detL2x := x3-x4: detL2y := y3-y4: Denom := detL1x*detL2y - detL2x*detL1y: x := (detL1 * detL2x - detL2 * detL1x) / Denom: y := (detL1 * detL2y - detL2 * detL1y) / Denom: [x,y]: end: # end of CalcIntersect(P1,P2,P3,P4) ######################################################################### # DistRecur(L,M): Inputs list L, giving the recurrence relation for a # C-finite sequence, and a maximum modulus M, and outputs # the set of all "distances" from modulus 1 to M, representing the # subset of the "minimum distance spectrum" up to modulo M. # See the comments for RecToSeq(C,n) for details regarding the format # for C-fintie recurrence L. # See the related paper: # "Toral Automorphisms, Linear Recurrences Modulo m, and the Euclideanity # of Q(sqrt(2+sqrt(2))" # for details of the definition of "minimum distance spectrum". DistRecur := proc(L::list,M::posint) local m: `union`(seq(DistRecurMod(L,m),m=1..M)): end: # end of DistRecur(L,M) ######################################################################### # DistRecurMod(L,m): Inputs a list L, giving the recurrence relation # for a C-finite sequence, and a modulus m. For each (MIN,MAX) pair # (see output of MinMaxRecurMod), we calculate a "minimum distance" # min(MIN,m-MAX)/m, which represents the distance from the nearest # integer point when considered as a mod 1 sequence after being rescaled # by the modulus. We then output the set of all such distances. # This is the mod m subset of the "minimum distance spectum". DistRecurMod := proc(L::list,m::posint) local r, IniSet::set, Alphabet::set, i, Ini::list, Pair, Dist, Distances::set, A::Array: Distances := {}: r := nops(L): # r is the order of the recurrence sequence A := Array((0..m-1)$r): Alphabet := {seq(i,i=0..m-1)}: IniSet := AllWords(Alphabet,r): for Ini in IniSet do: Pair := MinMaxSeqMod([L,Ini],m): Dist := evalf(min(Pair[1],m-Pair[2])/m): Distances := Distances union {Dist}: od: Distances: end: # end of DistRecurMod(L,m) ######################################################################### # MinMaxRecurMod(L,m): Inputs a list L, giving the recurrence relation # for a C-finite sequence, and a modulus m, and outputs an Array of # pairs (min,max) for each of the m^r possible initial data the will # define a C-finite sequence mod m, where r is the order of the # recurrence. MinMaxRecurMod := proc(L::list,m::posint) local r, A::Array, IniSet::set, Alphabet::set, i, Ini::list: r := nops(L): # r is the order of the recurrence sequence A := Array((0..m-1)$r): Alphabet := {seq(i,i=0..m-1)}: IniSet := AllWords(Alphabet,r): for Ini in IniSet do: A[op(Ini)] := MinMaxSeqMod([L,Ini],m): od: A: end: # end of MinMaxRecurMod(L,m) ######################################################################### # MinMaxSeqMod(C,m) # Inputs a C-finite sequence C, and a modulus m, and outputs the pair # the is the minimum and maximum of the sequence. MinMaxSeqMod := proc(C::list,m::posint) local r, maxperiod, A::Array: r := nops(C[1]): # r is the order of the C-finite sequence maxperiod := m^r: A := RecToSeqMod(C,maxperiod,m): min(A), max(A): end: # end of MinMaxSeqMod(C,m) ######################################################################### # RecToSeqMod(C,n,m) # Input: A C-finite sequence C, and positive integers n, m. # Output: An ARRAY of length n with the first n terms of the sequence, # modulo m. # e.g. RecToSeq([[-1,-1],[1,1]],16,3); should output: # [1,1,2,0,2,2,1,0,1,1,2,0,2,2,1,0]. # Using an Array is faster and less memory intensive than lists. # Also, performing the mod inside the procedure is faster than # calculating the whole array and during mod internally, for large # sequences. RecToSeqMod := proc(C::list,n::posint,m::posint) local i, j, r, recurrence::list, initialData::list, L::Array: if nops(C) <> 2 then print("Not the format of a C-finite sequence."); return(FAIL): fi: recurrence := C[1]: initialData := C[2]: r := nops(C[1]): # r is the order of the C-finite sequence if nops(initialData) <> r then print("Length of initial data should not differ from the order of the recurrence."); return(FAIL): fi: L := Array(1..n): L[1] := initialData[1] mod m: for i from 2 to n do if i <= r then L[i] := initialData[i] mod m: else L[i] := add(-L[i-j]*recurrence[j],j=1..r) mod m: fi: od: L: end: # End of RecToSeqMod(C,n,m) ######################################################################### # RecToSeq(C,n) # Input: A C-finite sequence C, and a positive integer n. # Output: An ARRAY of length n with the first n terms of the sequence. # e.g. RecToSeq([[-1,-1],[1,1]],10); should output: # [1,1,2,3,5,8,13,21,34,55]. # e.g. the sqrt(2) coeff of (a+b*sqrt(2))*(1+sqrt(2))^n is given # by RecToSeq([[-2,-1],[b,a+b]],5); which outputs: # [b, a+b, 2*a+3*b, 5*a+7*b, 12*a+17*b] # Using an Array is faster and less memory intensive than lists. RecToSeq := proc(C::list,n::posint) local i, j, r, recurrence::list, initialData::list, L::Array: if nops(C) <> 2 then print("Not the format of a C-finite sequence."); return(FAIL): fi: recurrence := C[1]: initialData := C[2]: r := nops(C[1]): # r is the order of the C-finite sequence if nops(initialData) <> r then print("Length of initial data should not differ from the order of the recurrence."); return(FAIL): fi: L := Array(1..n): L[1] := initialData[1]: for i from 2 to n do if i <= r then L[i] := initialData[i]: else L[i] := add(-L[i-j]*recurrence[j],j=1..r): fi: od: L: end: # End of RecToSeq(C,n) ######################################################################### # AllWords(A,n): Inputs an alphabet A, and a positive integer n, # and outputs a set of all nops(A)^n words of length n. AllWords := proc(A::set,n::nonnegint) option remember: local S, S1, s1, a: if n = 0 then return({[]}): fi: S := {}: S1 := AllWords(A,n-1); for s1 in S1 do: s1 := op(s1): S := S union {seq([s1,a], a in A)}: od: S: end: # end of AllWords(A,n) ########################## old stuff #########################3 # Let K be a number field, and R an order. # We study the action of the group of units U of R on the # quotient space K/R. # The orbits in K/R under the action of U are always finite. # In particular, we are interested in K = Q(sqrt(2)) and # R = Z[sqrt(2)]. Then U = {u^n : n integer}, where # u = 1 + sqrt(2). # This second version attempts to make the procedures more efficient. with(numtheory): HelpOld := proc(): print(`SolveRoots(Roots,a), bcdRootsK8(b,c,d,e,f,g,h), MaxRangeSolns(Roots,a), ScanK4Exceptions(m,a), ScanK8Exceptions(m,a), ScanK8ExceptionsNoHalf(m,a), OKrootDiff(g,m), Rmod(m), RmodExcl(m), OldOrbits(m), Orbits(m), OrbitsK4(m), OrbitsK4_v2(m), OrbitsK4_v3(m), Norm4(a,b,c,d), Norm4Mod1(a,b,c,d), NumberOrbits(OrbitList), MaxOrbitSize(OrbitList), MinOrbitSize(OrbitList), MaxD(OrbitList), MaxD2(Orbit,m), PeriodMod(x,m), ithprimeQR(i,a), ithprimeNonQR(i,a)`): end: ######################################################################## # OKrootDiff(g,m) Inputs "g", and denominator m, # and outputs a set S of points [s,t] that are the subset of R^2, of # denominator m, s.t. the MaxRangeSolutions([r1,r2,r3,r4],1/2) >= 1, # where r1 < r2 < r3 < r4, s = (r2-r1)/2, t = (r4-r3)/2, # and (r1+r2)/2 = -g, and (r3+r4)/2 = g. OKrootDiff := proc(g,m) local r1, r2, r3, r4, s, t, stmax, S: S := {}: stmax := 2: # usually we let s and t range from -3 to 3 for s from 0 to stmax by 1/m do: for t from 0 to s by 1/m do: r1 := -g - s: r2 := -g + s: r3 := g - t: r4 := g + t: if MaxRangeSolns([r1,r2,r3,r4],0.5) >= 1 then S := S union {[s,t],[-s,t],[s,-t],[-s,-t],[t,s],[-t,s],[t,-s],[-t,-s]}: fi: od: od: S: end: ######################################################################## # SolveRoots(Roots,a): SolveRoots := proc(Roots::list,a) local f,r,X: f := mul(X-r, r in Roots): [solve(abs(f)<=a, X)]: end: ######################################################################## # PlotRoots(Roots): PlotRoots := proc(Roots::list) local f,r,X: f := mul(X-r, r in Roots): plot(f,X=-2..2,view=[-2..2,-0.5..0.5]): end: ######################################################################## # MaxRangeSolns(Roots,a): Inputs a list of roots # and outputs the maximal length of all the real ranges that # satisfy abs(product(X-Roots[i]))<= a. # NOTE - We use LIST and not set of roots because of multiplicity. MaxRangeSolns := proc(Roots::list,a) local i, f, X, r, Rngs, RngLengths, maxLength, u, v: #RngLengths := {}: maxLength := 0: f := mul(X-r, r in Roots): Rngs := [solve(abs(f)<=a, X)]: for i from 1 to nops(Rngs) do u := op([1],Rngs[i]): v := op([2],Rngs[i]): if whattype(u) = function then u := op([1],u): fi: if whattype(v) = function then v := op([1],v): fi: maxLength := max(maxLength, v-u): # RngLengths := RngLengths union {v-u}: od: #RngLengths := [seq(op([1],op([2],Rngs[i]))-op([1],op([1],Rngs[i])),i=1..nops(Rngs))]: maxLength: #max(RngLengths): end: ######################################################################## # bcdRoots(b,c,d) bcdRoots := proc(b,c,d) local w, wb, s2, r1, r2, r3, r4: s2 := evalf(sqrt(2)): w := sqrt(2+s2): wb := sqrt(2-s2): r1 := b*w + c*s2 + d*wb: r2 := -d*w - c*s2 + b*wb: r3 := d*w - c*s2 - b*wb: r4 := -b*w + c*s2 - d*wb: [r1,r2,r3,r4]: end: ######################################################################## # bcdRootsK8(b,c,d,e,f,g,h) bcdRootsK8 := proc(b,c,d,e,f,g,h) local w, wb, s2, n1, n2, n3, n4, r1, r2, r3, r4, r5, r6, r7, r8: s2 := evalf(sqrt(2)): w := sqrt(2+s2): wb := sqrt(2-s2): n1 := sqrt(2+w): n2 := sqrt(2-w): n3 := sqrt(2+wb): n4 := sqrt(2-wb): r1 := b*n1 + c*w + d*n3 + e*s2 + f*n4 + g*wb + h*n2: r2 := b*n3 + c*wb - d*n2 - e*s2 - f*n1 - g*w - h*n4: r3 := b*n4 - c*wb - d*n1 - e*s2 + f*n2 + g*w + h*n3: r4 := b*n2 - c*w - d*n4 + e*s2 + f*n3 - g*wb - h*n1: r5 := -b*n2 - c*w + d*n4 + e*s2 - f*n3 - g*wb + h*n1: r6 := -b*n4 - c*wb + d*n1 - e*s2 - f*n2 + g*w - h*n3: r7 := -b*n3 + c*wb + d*n2 - e*s2 + f*n1 - g*w + h*n4: r8 := -b*n1 + c*w - d*n3 + e*s2 - f*n4 + g*wb - h*n2: [r1,r2,r3,r4,r5,r6,r7,r8]: end: ######################################################################## # ScanK4Exceptions(m,a) ScanK4Exceptions := proc(m,a) local b,c,d, MinRng, Rng, Roots: MinRng := 9999: for b from 0 to m-1 do: for c from 0 to m-1 do: for d from 0 to m-1 do: if igcd(b,c,d) = 1 then Roots := bcdRoots(b/m - 0.5, c/m - 0.5, d/m - 0.5): Rng := MaxRangeSolns(Roots,a): if Rng < 1 then print(b/m - 0.5, c/m - 0.5, d/m - 0.5): print(Roots): print(Rng): fi: MinRng := min(MinRng,Rng): fi: od: od: od: print(MinRng): end: ######################################################################## # ScanK8Exceptions(m,a) ScanK8Exceptions := proc(m,a) local b,c,d,e,f,g,h, MinRng, Rng, Roots, SumAbsRoots, r, maxGood, minGood, avgGood, totGood, ctGood, maxBad, minBad, avgBad, totBad, ctBad: maxGood := 0: maxBad := 0: minGood := 9999: minBad := 9999: totGood := 0: totBad := 0: ctGood := 0: ctBad := 0: MinRng := 9999: for b from 0 to m-1 do: for c from 0 to m-1 do: for d from 0 to m-1 do: print(b,c,d): for e from 0 to m-1 do: for f from 0 to m-1 do: for g from 0 to m-1 do: for h from 0 to m-1 do: if igcd(b,c,d,e,f,g,h) = 1 then Roots := bcdRootsK8(b/m-0.5,c/m-0.5,d/m-0.5,e/m-0.5,f/m-0.5,g/m-0.5,h/m-0.5): SumAbsRoots := add(abs(r), r in Roots): Rng := MaxRangeSolns(Roots,a): if Rng < 1 then print(b/m-0.5,c/m-0.5,d/m-0.5,e/m-0.5,f/m-0.5,g/m-0.5,h/m-0.5, Rng, SumAbsRoots): # print(Roots): # print(Rng): totBad := totBad + SumAbsRoots: ctBad := ctBad + 1: maxBad := max(maxBad,SumAbsRoots): minBad := min(minBad,SumAbsRoots): else totGood := totGood + SumAbsRoots: ctGood := ctGood + 1: maxGood := max(maxGood,SumAbsRoots): minGood := min(minGood,SumAbsRoots): fi: MinRng := min(MinRng,Rng): fi: od: od: od: od: od: od: od: avgGood := evalf(totGood/ctGood): avgBad := evalf(totBad /ctBad ): print("min range ", MinRng): print("Good max/min/avg ", maxGood, minGood, avgGood): print("Good max/min/avg ", maxBad, minBad, avgBad ): end: ######################################################################## # ScanK8ExceptionsNoHalf(m,a) ScanK8ExceptionsNoHalf := proc(m,a) local b,c,d,e,f,g,h, MinRng, Rng, Roots: MinRng := 9999: for b from 1 to m-1 do: for c from 1 to m-1 do: for d from 1 to m-1 do: print(b,c,d): for e from 1 to m-1 do: for f from 1 to m-1 do: for g from 1 to m-1 do: for h from 1 to m-1 do: if igcd(b,c,d,e,f,g,h) = 1 then Roots := bcdRootsK8(b/m-0.5,c/m-0.5,d/m-0.5,e/m-0.5,f/m-0.5,g/m-0.5,h/m-0.5): Rng := MaxRangeSolns(Roots,a): if Rng < 1 then print(b/m-0.5,c/m-0.5,d/m-0.5,e/m-0.5,f/m-0.5,g/m-0.5,h/m-0.5): print(Roots): print(Rng): fi: MinRng := min(MinRng,Rng): fi: od: od: od: od: od: od: od: print(MinRng): end: ######################################################################## # Rmod(m): Input a positive integer m. # Returns a list of elements of (Z/mZ)[sqrt(2)], excluding zero. Rmod := proc(m::posint) local i, j, L: L := [seq(seq(i + j*sqrt(2),i=0..m-1),j=0..m-1)]: L[2..nops(L)]: end: ######################################################################## # RmodExcl(m): Input a positive integer m. # Returns a list of (Z/mZ)[sqrt(2)], excluding zero and also # excluding elements that would have been generated by divisors of m. RmodExcl := proc(m::posint) local L, L1, factors, i: L := Rmod(m): factors := convert(factorset(m),list): for i from 1 to nops(factors) do if factors[i] <> m then L1 := Rmod(m/factors[i])* factors[i]: L := convert(convert(L,set) minus convert(L1,set), list): fi: od: L: end: ######################################################################## # OldOrbits(m): Input a positive integer m # Return the elements of (Z/mZ)[sqrt(2)] partitioned into orbits # under the action of multiplication by 1+sqrt(2), as a list of lists. Orbits := proc(m::posint) local L, Orbits, orbit, x, u: #L := RmodExcl(m): L := convert(RmodExcl(m),set): u := 1 + sqrt(2): Orbits := []: #while L <> [] do while L <> {} do x := L[1]: orbit := [x]: if (m=2) and (x = sqrt(2)) then x := sqrt(2): # avoids Maple error that sqrt(2) mod 2 = 0 else x := expand(u * x) mod m: fi: if not is(x in L) then print("Oops! Not in L!"): return(x, L, Orbits): fi: L := L[2..nops(L)]: while x <> orbit[1] do: orbit := [op(orbit), x]: # L := convert(convert(L,set) minus {x}, list): L := L minus {x}: x := expand(u * x) mod m: od: #print(nops(orbit)): Orbits := [op(Orbits),orbit]: od: Orbits: end: ######################################################################## # Orbits(m): Input a positive integer m # The elements of (Z/mZ)[sqrt(2)] are partitioned into orbits # under the action of multiplication by 1+sqrt(2). # Returns a list L of orbit representatives (as pairs of integers), # and also returns an m-by-m array A, with A[i,j]=k means i+j*sqrt(2) # has orbit L[k]. # Note - we want to exclude orbits which are isomorphic to those # which would have been contained in a smaller m. We'll leave these # entries as 0 in A. We test these using igcd(). # For size m=100, this is approx 1000x faster than Orbit(m) and similarly # 1/1000 of the memory usage! This is likely do to the enormous # inefficiency of using Maple lists, which copy themselves when changed. Orbits2 := proc(m::posint) local A::Array, i, j, nextPair, ct, L::list: L := []: A := Array(0..m-1,0..m-1): ct := 0: for i from 0 to m-1 do for j from 0 to m-1 do if A[i,j] = 0 and igcd(i,j,m) = 1 then ct := ct + 1: A[i,j] := ct: L := [op(L),[i,j]]: nextPair := i+2*j mod m, i+j mod m: while nextPair[1] <> i or nextPair[2] <> j do A[nextPair] := ct: nextPair := nextPair[1]+2*nextPair[2] mod m, nextPair[1]+nextPair[2] mod m: od: fi: od: od: L, A: end: ######################################################################## # HalfSimplex(L,A): Inputs a ######################################################################## # OrbitsK4(m): Input a positive integer m. # Let w := sqrt(2+sqrt(2)) and wb := sqrt(2-sqrt(2)). # The elements of Z[w]/mZ[w] are partitioned into orbits # under the action of multiplication by the unit 1+w. # Returns a list L of orbit representatives (as 4-tuples of integers), # and also returns an m^4 array A, with A[i,j,k,l]=c means # i+j*w+k*sqrt(2)+l*wb has orbit L[c]. # We also print out an analysis of the orbits. # RMK: (i+j*w+k*sqrt(2)+l*wb)*(1+w) # = (i+j*w+k*sqrt(2)+l*wb) + (i*w+j*w^2+k*sqrt(2)*w+l*wb*w) # = (i+j*w+k*sqrt(2)+l*wb) + (i*w+j*(2+sqrt(2))+k*(w+wb)+l*sqrt(2)) # = (i+j*w+k*sqrt(2)+l*wb) + ((i+k)*w+j*2+(j+l)*sqrt(2)+k*wb) # = (i+j*2)+(i+j+k)*w+(j+k+l)*sqrt(2)+(k+l)*wb # Note - we want to exclude orbits which are isomorphic to those # which would have been contained in a smaller m. We'll leave these # entries as 0 in A. We test these using igcd(). OrbitsK4 := proc(m::posint) local A::Array, i, j, k, l, nextPair, ct, L::list, MaxNorm, AbsNorm, MinNorm, MaxMinNorm, MaxMinSumCoord, MinSumCoord, SumCoord, w, wb, sq2, MaxGal, MinMaxGal, MaxMinMaxGal, Gal1, Gal2, Gal3, Gal4, a, b, c, d: sq2 := evalf(2^0.5): w := evalf((2+sq2)^0.5): wb := evalf((2-sq2)^0.5): L := []: A := Array(0..m-1,0..m-1,0..m-1,0..m-1): ct := 0: MaxNorm := 0: MaxMinNorm := 0: MaxMinSumCoord := 0: MaxMinMaxGal := 0: for i from 0 to m-1 do for j from 0 to m-1 do for k from 0 to m-1 do for l from 0 to m-1 do if A[i,j,k,l] = 0 and igcd(i,j,k,l,m) = 1 then AbsNorm := abs(Norm4Mod1(i/m,j/m,k/m,l/m)): MaxNorm := max(MaxNorm,AbsNorm): MinNorm := AbsNorm: a := i/m: a := a - floor(a+1/2): b := j/m: b := b - floor(b+1/2): c := k/m: c := c - floor(c+1/2): d := l/m: d := d - floor(c+1/2): SumCoord := abs(a)+abs(b)+abs(c)+abs(d): MinSumCoord := SumCoord: Gal1 := a + b*w + c*sq2 + d*wb: Gal2 := a - d*w - c*sq2 + b*wb: Gal3 := a + d*w - c*sq2 - b*wb: Gal4 := a - b*w + c*sq2 - d*wb: MaxGal := max(abs(Gal1),abs(Gal2),abs(Gal3),abs(Gal4)): MinMaxGal := MaxGal: ct := ct + 1: A[i,j,k,l] := ct: L := [op(L),[i,j,k,l]]: nextPair := i+2*j mod m, i+j+k mod m, j+k+l mod m, k+l mod m: while nextPair[1] <> i or nextPair[2] <> j or nextPair[3] <> k or nextPair[4] <> l do AbsNorm := abs(Norm4Mod1(nextPair/m)): MaxNorm := max(MaxNorm,AbsNorm): MinNorm := min(MinNorm,AbsNorm): a := nextPair[1]/m: a := a - floor(a+1/2): b := nextPair[2]/m: b := b - floor(b+1/2): c := nextPair[3]/m: c := c - floor(c+1/2): d := nextPair[4]/m: d := d - floor(d+1/2): SumCoord := abs(a)+abs(b)+abs(c)+abs(d): MinSumCoord := min(MinSumCoord,SumCoord): Gal1 := a + b*w + c*sq2 + d*wb: Gal2 := a - d*w - c*sq2 + b*wb: Gal3 := a + d*w - c*sq2 - b*wb: Gal4 := a - b*w + c*sq2 - d*wb: MaxGal := max(abs(Gal1),abs(Gal2),abs(Gal3),abs(Gal4)): MinMaxGal := min(MinMaxGal,MaxGal): A[nextPair] := ct: nextPair := nextPair[1]+2*nextPair[2] mod m, nextPair[1]+nextPair[2]+nextPair[3] mod m, nextPair[2]+nextPair[3]+nextPair[4] mod m, nextPair[3]+nextPair[4] mod m: od: MaxMinNorm := max(MaxMinNorm,MinNorm): MaxMinSumCoord := max(MaxMinSumCoord,MinSumCoord): MaxMinMaxGal := max(MaxMinMaxGal,MinMaxGal): fi: od: od: od: od: printf(`m=%d: orbits %d, max %f, maxmin %f, maxsumcoord %f maxgal %f\n`, m, nops(L),MaxNorm,MaxMinNorm,evalf(MaxMinSumCoord),MaxMinMaxGal): L, A: end: ######################################################################## # OrbitsK4_v2(m): Input a positive integer m. # RMK: (i+j*w+k*sqrt(2)+l*wb)*(1+w) # = (i+j*w+k*sqrt(2)+l*wb) + (i*w+j*w^2+k*sqrt(2)*w+l*wb*w) # = (i+j*w+k*sqrt(2)+l*wb) + (i*w+j*(2+sqrt(2))+k*(w+wb)+l*sqrt(2)) # = (i+j*w+k*sqrt(2)+l*wb) + ((i+k)*w+j*2+(j+l)*sqrt(2)+k*wb) # = (i+j*2)+(i+j+k)*w+(j+k+l)*sqrt(2)+(k+l)*wb # Note - we want to exclude orbits which are isomorphic to those # which would have been contained in a smaller m. We'll leave these # entries as 0 in A. We test these using igcd(). OrbitsK4_v2 := proc(m::posint) local A::Array, i, j, k, l, nextPair, ct, L::list, MaxNorm, AbsNorm, MinNorm, MaxMinNorm, MaxNorm1, MaxNorm2, MaxNorm3, nextPair2, nextPair3: MaxNorm := 0: MaxMinNorm := 0: MaxNorm1 := 0: MaxNorm2 := 0: MaxNorm3 := 0: for i from 0 to m-1 do for j from 0 to m-1 do for k from 0 to m-1 do for l from 0 to m-1 do if igcd(i,j,k,l,m) = 1 then AbsNorm := abs(Norm4Mod1(i/m,j/m,k/m,l/m)): MaxNorm := max(MaxNorm,AbsNorm): # MinNorm := AbsNorm: nextPair := i+2*j mod m, i+j+k mod m, j+k+l mod m, k+l mod m: nextPair2 := nextPair[1]+2*nextPair[2] mod m, nextPair[1]+nextPair[2]+nextPair[3] mod m, nextPair[2]+nextPair[3]+nextPair[4] mod m, nextPair[3]+nextPair[4] mod m: nextPair3 := nextPair2[1]+2*nextPair2[2] mod m, nextPair2[1]+nextPair2[2]+nextPair2[3] mod m, nextPair2[2]+nextPair2[3]+nextPair2[4] mod m, nextPair2[3]+nextPair2[4] mod m: MaxNorm1 := max(MaxNorm1,min(AbsNorm,abs(Norm4Mod1(nextPair/m)))): MaxNorm2 := max(MaxNorm2,min(AbsNorm,abs(Norm4Mod1(nextPair/m)),abs(Norm4Mod1(nextPair2/m)))): MaxNorm3 := max(MaxNorm3,min(AbsNorm,abs(Norm4Mod1(nextPair/m)),abs(Norm4Mod1(nextPair2/m)),abs(Norm4Mod1(nextPair3/m)))): # while AbsNorm > 1/2 do # AbsNorm := abs(Norm4Mod1(nextPair/m)): # MaxNorm := max(MaxNorm,AbsNorm): # MinNorm := min(MinNorm,AbsNorm): # A[nextPair] := ct: # nextPair := nextPair[1]+2*nextPair[2] mod m, nextPair[1]+nextPair[2]+nextPair[3] mod m, nextPair[2]+nextPair[3]+nextPair[4] mod m, nextPair[3]+nextPair[4] mod m: # od: # MaxMinNorm := max(MaxMinNorm, MinNorm): fi: od: od: od: od: printf(`m=%d: orbits %d, max %f, max1 %f, max2 %f, max3 %f\n`, m, nops(L),MaxNorm,MaxNorm1,MaxNorm2,MaxNorm3): end: ######################################################################## # OrbitsK4_v3(m): Input a positive integer m. # RMK: (i+j*w+k*sqrt(2)+l*wb)*(1+w) # = (i+j*w+k*sqrt(2)+l*wb) + (i*w+j*w^2+k*sqrt(2)*w+l*wb*w) # = (i+j*w+k*sqrt(2)+l*wb) + (i*w+j*(2+sqrt(2))+k*(w+wb)+l*sqrt(2)) # = (i+j*w+k*sqrt(2)+l*wb) + ((i+k)*w+j*2+(j+l)*sqrt(2)+k*wb) # = (i+j*2)+(i+j+k)*w+(j+k+l)*sqrt(2)+(k+l)*wb # Note - we want to exclude orbits which are isomorphic to those # which would have been contained in a smaller m. We'll leave these # entries as 0 in A. We test these using igcd(). OrbitsK4_v3 := proc(m::posint) local A::Array, i, j, k, l, nextPair, ct, L::list, MaxNorm, AbsNorm, MinNorm, MaxMinNorm, MaxNorm1, MaxNorm2, MaxNorm3, nextPair2, nextPair3, OKnormCt0, OKnormCt1, OKnormCt2, OKnormCt3, TotnormCt, AbsNorm1, AbsNorm2, AbsNorm3: MaxNorm := 0: MaxMinNorm := 0: OKnormCt0 := 0: OKnormCt1 := 0: OKnormCt2 := 0: OKnormCt3 := 0: TotnormCt := 0: for i from 0 to m-1 do for j from 0 to m-1 do for k from 0 to m-1 do for l from 0 to m-1 do if igcd(i,j,k,l,m) = 1 then TotnormCt := TotnormCt + 1: AbsNorm := abs(Norm4Mod1(i/m,j/m,k/m,l/m)): if AbsNorm <= 1/2 then OKnormCt0 := OKnormCt0 + 1: fi: # MaxNorm := max(MaxNorm,AbsNorm): # MinNorm := AbsNorm: nextPair := i+2*j mod m, i+j+k mod m, j+k+l mod m, k+l mod m: nextPair2 := nextPair[1]+2*nextPair[2] mod m, nextPair[1]+nextPair[2]+nextPair[3] mod m, nextPair[2]+nextPair[3]+nextPair[4] mod m, nextPair[3]+nextPair[4] mod m: nextPair3 := nextPair2[1]+2*nextPair2[2] mod m, nextPair2[1]+nextPair2[2]+nextPair2[3] mod m, nextPair2[2]+nextPair2[3]+nextPair2[4] mod m, nextPair2[3]+nextPair2[4] mod m: AbsNorm1 := abs(Norm4Mod1(nextPair/m)): AbsNorm2 := abs(Norm4Mod1(nextPair2/m)): AbsNorm3 := abs(Norm4Mod1(nextPair3/m)): #print(AbsNorm, AbsNorm1, AbsNorm2, AbsNorm3): if AbsNorm1 <= 1/2 then OKnormCt1 := OKnormCt1 + 1: fi: if AbsNorm2 <= 1/2 then OKnormCt2 := OKnormCt2 + 1: fi: if AbsNorm3 <= 1/2 then OKnormCt3 := OKnormCt3 + 1: fi: #if min(AbsNorm,AbsNorm1) <= 1/2 then OKnormCt1 := OKnormCt1 + 1: fi: #if min(AbsNorm,AbsNorm1,AbsNorm2) <= 1/2 then OKnormCt2 := OKnormCt2 + 1: fi: #if min(AbsNorm,AbsNorm1,AbsNorm2,AbsNorm3) <= 1/2 then OKnormCt3 := OKnormCt3 + 1: fi: # while AbsNorm > 1/2 do # AbsNorm := abs(Norm4Mod1(nextPair/m)): # MaxNorm := max(MaxNorm,AbsNorm): # MinNorm := min(MinNorm,AbsNorm): # A[nextPair] := ct: # nextPair := nextPair[1]+2*nextPair[2] mod m, nextPair[1]+nextPair[2]+nextPair[3] mod m, nextPair[2]+nextPair[3]+nextPair[4] mod m, nextPair[3]+nextPair[4] mod m: # od: # MaxMinNorm := max(MaxMinNorm, MinNorm): fi: od: od: od: od: #print(OKnormCt0, OKnormCt1, OKnormCt2, OKnormCt3, TotnormCt): printf(`m=%d: norm0 %f, norm1 %f, norm2 %f, norm3 %f\n`, m, evalf(OKnormCt0/TotnormCt),evalf(OKnormCt1/TotnormCt),evalf(OKnormCt2/TotnormCt),evalf(OKnormCt3/TotnormCt)): end: ######################################################################## # Norm4(a,b,c,d): Norm of Q(w), with basis 1,w,sqrt(2),wb, with # w = sqrt(2+sqrt(2)), wb = sqrt(2-sqrt(2)) Norm4 := proc(a,b,c,d): -8*b^3*d+8*b*d^3-8*c*d^2*a+8*c*b^2*a-8*c^2*d^2-8*b^2*c^2+4*b^2*d^2-4*a^2*c^2-4*a^2*d^2-4*a^2*b^2+16*a*b*c*d+a^4+2*b^4+4*c^4+2*d^4: end: ######################################################################## # Norm4Mod1(a,b,c,d): Norm of Q(w), with basis 1,w,sqrt(2),wb, with # w = sqrt(2+sqrt(2)), wb = sqrt(2-sqrt(2)) Norm4Mod1 := proc(a,b,c,d) Norm4(a-floor(a+1/2),b-floor(b+1/2),c-floor(c+1/2),d-floor(d+1/2)): #max(abs(a-floor(a+1/2)),abs(b-floor(b+1/2)),abs(c-floor(c+1/2)),abs(d-floor(d+1/2))): end: ######################################################################## NumberOrbits := proc(OrbitList::list) local i: [seq(nops(OrbitList[i]),i=1..nops(OrbitList))]: end: ######################################################################## MaxOrbitSize := proc(OrbitList::list) local i, j: [seq(max(seq(nops(OrbitList[i][j]),j=1..nops(OrbitList[i]))),i=1..nops(OrbitList))]: end: ######################################################################## MinOrbitSize := proc(OrbitList::list) local i, j: [seq(min(seq(nops(OrbitList[i][j]),j=1..nops(OrbitList[i]))),i=1..nops(OrbitList))]: end: ######################################################################## MaxD := proc(OrbitList::list) local i, j, k: [seq(max(seq( min( seq( (OrbitList[i][j][k] - subs(sqrt(2)=0,OrbitList[i][j][k]))/sqrt(2) ,k=1..nops(OrbitList[i][j])),seq( i-(OrbitList[i][j][k] - subs(sqrt(2)=0,OrbitList[i][j][k]))/sqrt(2) ,k=1..nops(OrbitList[i][j])) ) ,j=1..nops(OrbitList[i]))),i=1..nops(OrbitList))]: end: ######################################################################## MaxD2 := proc(Orbit,m) local j, k: max(seq( min( seq( (Orbit[j][k] - subs(sqrt(2)=0,Orbit[j][k]))/sqrt(2) ,k=1..nops(Orbit[j])),seq( m-(Orbit[j][k] - subs(sqrt(2)=0,Orbit[j][k]))/sqrt(2) ,k=1..nops(Orbit[j])) ) ,j=1..nops(Orbit))): end: ######################################################################## # PeriodMod(x,m) - Input x and modulus m. Outputs smallest positive # integer f such that x^f = 1 mod m. PeriodMod := proc(x,m) local f, x1: f := 1: x1 := x: if m=1 then return 1: fi: for f from 1 while x1 <> 1 do x1 := expand(x1 * x) mod m: od: f: end: ######################################################################## # ithprimeQR(i,a) - Returns the ith prime p for which "a" is a quadratic # residue mod p. ithprimeQR := proc(i,a) local j, k: k := 0: j := 0: while j < i do: k := k+1: if legendre(a,ithprime(k)) = 1 then j := j + 1: fi: od: ithprime(k): end: ######################################################################## # ithprimeNonQR(i,a) - Returns the ith prime p for which "a" is a quadratic # nonresidue mod p. ithprimeNonQR := proc(i,a) local j, k: k := 0: j := 0: while j < i do: k := k+1: if legendre(a,ithprime(k)) = -1 then j := j + 1: fi: od: ithprime(k): end: