----------------------------------------------------- ----------------------------------------------------- ----------------------------------------------------- -- PART I. Header ----------------------------------------------------- ----------------------------------------------------- ----------------------------------------------------- -* This file contains code we wrote for the preprint "The worst 1-PS for rational curves with a unibranch singu- larity." Perhaps someday we will make it a package in Macaulay2. For now, we gives some examples and tests for each function. *- ----------------------------------------------------- ----------------------------------------------------- ----------------------------------------------------- -- Functions we wrote for this project ----------------------------------------------------- ----------------------------------------------------- ----------------------------------------------------- ----------------------------------------------------- ----------------------------------------------------- -- Functions for setting up and solving the KKT equations ----------------------------------------------------- ----------------------------------------------------- ----------------------------------------------------- -- smallElementsViaGap -- Uses a crude interface between GAP and Macaulay2 -- to get the list of elements in a semigroup S up to the conductor -- S: a list of generators of the semigroup -- gcd(S) should be 1 smallElementsViaGap = memoize ((S) -> ( L := sort unique S; if gcd(L) != 1 then error "The elements of L must be coprime"; gapInputFN:=temporaryFileName(); gapOutputFN:=temporaryFileName(); s1:=///LoadPackage("numericalsgps"); S:=NumericalSemigroup(///; s2:=toString(L); s2=substring(s2,1,#s2-2); s3:=///); L:=SmallElementsOfNumericalSemigroup(S); PrintTo("///; s4:=gapOutputFN; s5:=///",L); quit; ///; gapInput:=concatenate(s1,s2,s3,s4,s5); gapInputFile:=openOut gapInputFN; gapInputFile << gapInput << endl; close gapInputFile; run concatenate(///gap < ///,gapInputFN,/// > ///,gapOutputFN); s := get gapOutputFN; i0:=0; while s_i0!="[" do i0=i0+1; i1:=i0+1; while s_i1!="]" do i1=i1+1; t:=substring(s,i0+2,i1-2-(i0+2)+1); E := value concatenate("{",t,"}"); run concatenate("rm ",gapInputFN); run concatenate("rm ",gapOutputFN); return E )); -* -- Examples smallElementsViaGap({2,3}) smallElementsViaGap({2,5}) smallElementsViaGap({4,9}) -- Tests assert(smallElementsViaGap({2,3})=={0,2}) assert(smallElementsViaGap({2,5})=={0,2,4}) assert(smallElementsViaGap({4,9})=={0, 4, 8, 9, 12, 13, 16, 17, 18, 20, 21, 22, 24}) assert(smallElementsViaGap({5,7})=={0, 5, 7, 10, 12, 14, 15, 17, 19, 20, 21, 22, 24}) assert(smallElementsViaGap({3,5,7})=={0,3,5}) *- ----------------------------------------------------- ----------------------------------------------------- -- gammaVectorOfSemigroup -- Computes the list of elements of the semigroup -- starting at gamma_0 = 0 and going up to gamma_N -- S: a list of generators of the semigroup -- N+1 the desired length of c gammaVectorOfSemigroup = (N,S) -> ( se:=smallElementsViaGap(S); if #se>N+1 then error "N+1 is smaller than the conductor"; join(se,apply(N+1-#se, i -> last(se)+i+1)) ); -* -- Examples gammaVectorOfSemigroup(9,{2,5}) gammaVectorOfSemigroup(10,{2,5}) gammaVectorOfSemigroup(15,{4,9}) -- Tests assert(gammaVectorOfSemigroup(9,{2,5})=={0, 2, 4, 5, 6, 7, 8, 9, 10, 11}) assert(gammaVectorOfSemigroup(10,{2,5})=={0, 2, 4, 5, 6, 7, 8, 9, 10, 11, 12}) assert(gammaVectorOfSemigroup(15,{4,9})=={0, 4, 8, 9, 12, 13, 16, 17, 18, 20, 21, 22, 24, 25, 26, 27}) *- ----------------------------------------------------- ----------------------------------------------------- -- avector -- Create the vector a -- The option Simp=>true gives the vector a for the Simplified Problem -- The option Simp=>false (default) gives the vector a for the -- Unsimplified Problem avector = {Simp => false } >> o -> (N,S) -> ( c:=gammaVectorOfSemigroup(N,S); a:=apply(N-1,i -> c_(i+2)-c_i); a = prepend(c_1,a); if not o.Simp then ( append(a,c_N - c_(N-1)) ) else ( append(a,2) ) ); -* -- Examples avector(3,{2,3}) avector(15,{4,9},Simp=>true) avector(15,{4,9}) -- Tests assert(avector(3,{2,3})=={2,3,2,1}) assert(avector(15,{4,9},Simp=>true)=={4, 8, 5, 4, 4, 4, 4, 2, 3, 3, 2, 3, 3, 2, 2, 2}) assert(avector(15,{4,9})=={4, 8, 5, 4, 4, 4, 4, 2, 3, 3, 2, 3, 3, 2, 2, 1}) *- ----------------------------------------------------- ----------------------------------------------------- -- conductorIndex -- Computes the index of the conductor conductorIndex = (S) -> ( seS:=smallElementsViaGap(S); #seS-1 ); -* -- Examples conductorIndex({2,3}) conductorIndex({2,5}) conductorIndex({4,9}) -- Tests assert(conductorIndex({2,3})==1) assert(conductorIndex({2,5})==2) assert(conductorIndex({4,9})==12) assert(conductorIndex({8,13})==42) *- ----------------------------------------------------- ----------------------------------------------------- -- Winequalities -- Internal function -- Given the vector c = gamma, -- compute the coefficient vectors of the inequalities defining W -- Return this as a list of lists, not a matrix -- Assume the inequalities are in the direction b_0 x_0 + ... + b_N x_N <= 0 -- to match the convex optimization literature Winequalities = (c) -> ( N := #c-1; M := for i from 0 to N-2 list ( for j from 0 to N list ( if j==i then -(c_(i+2)-c_(i+1)) else if j==i+1 then (c_(i+2)-c_i) else if j==i+2 then -(c_(i+1)-c_i) else 0 ) ); M ); -* -- Test assert(Winequalities({0,2,3,4,5,6})=={{-1, 3, -2, 0, 0, 0}, {0, -1, 2, -1, 0, 0}, {0, 0, -1, 2, -1, 0}, {0, 0, 0, -1, 2, -1}}) *- ----------------------------------------------------- ----------------------------------------------------- -- The piecewise linear functions Fk Fk = (k,x) -> (max {-x+k,0}) -* -- Examples Fk(3,2) Fk(3,5) -- Tests assert(Fk(3,2)==1) assert(Fk(3,5)==0) *- ----------------------------------------------------- ----------------------------------------------------- -- The linear function Ld Ld = (d,x) -> (d-x) -* -- Examples Ld(3,2) Ld(3,5) -- Tests assert(Ld(3,2)==1) assert(Ld(3,5)==-2) *- ----------------------------------------------------- ----------------------------------------------------- -- The linear function L1 -- This might be the silliest function I've ever written in Macaulay2, -- but it allows us to be consistent with the notes L1 = (x) -> 1 -* -- Example L1(2) -- Test assert(L1(2)==1) *- ----------------------------------------------------- ----------------------------------------------------- -- KKTmatrix -- Computes the KKT matrix A, given N, S, and I KKTmatrix = (N,S,I) -> ( c := gammaVectorOfSemigroup(N,S); M := {}; for i from 0 to N-2 do ( if member(i+1,I) then ( M = append(M,apply(N+1, j -> 2/1*Fk(c_(i+1),c_j))) ) else ( M = append(M,apply(N+1, j -> if j==i then -(c_(i+2)-c_(i+1)) else if j==i+1 then (c_(i+2)-c_i) else if j==i+2 then -(c_(i+1)-c_i) else 0)) ); ); d:=max c; M = append(M, apply(N+1, j -> 2*Ld(d,c_j))); M = append(M, apply(N+1, j -> 2*L1(c_j))); transpose matrix((1/1)*M) ); -* -- Examples KKTmatrix(4,{2,3},{}) KKTmatrix(15,{4,9},{7,8,10}) -- Tests assert(KKTmatrix(4,{2,3},{})==matrix {{-1/1, 0, 0, 10, 2}, {3, -1, 0, 6, 2}, {-2, 2, -1, 4, 2}, {0, -1, 2, 2, 2}, {0, 0, -1, 0, 2}}) assert(KKTmatrix(15,{4,9},{7,8,10})==matrix {{-4/1, 0, 0, 0, 0, 0, 34, 36, 0, 42, 0, 0, 0, 0, 54, 2}, {8, -1, 0, 0, 0, 0, 26, 28, 0, 34, 0, 0, 0, 0, 46, 2}, {-4, 5, -3, 0, 0, 0, 18, 20, 0, 26, 0, 0, 0, 0, 38, 2}, {0, -4, 4, -1, 0, 0, 16, 18, 0, 24, 0, 0, 0, 0, 36, 2}, {0, 0, -1, 4, -3, 0, 10, 12, 0, 18, 0, 0, 0, 0, 30, 2}, {0, 0, 0, -3, 4, -1, 8, 10, 0, 16, 0, 0, 0, 0, 28, 2}, {0, 0, 0, 0, -1, 4, 2, 4, 0, 10, 0, 0, 0, 0, 22, 2}, {0, 0, 0, 0, 0, -3, 0, 2, 0, 8, 0, 0, 0, 0, 20, 2}, {0, 0, 0, 0, 0, 0, 0, 0, -1, 6, 0, 0, 0, 0, 18, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 3, 2, 0, 0, 0, 0, 14, 2}, {0, 0, 0, 0, 0, 0, 0, 0, -2, 0, -2, 0, 0, 0, 12, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, -1, 0, 0, 10, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 3, -1, 0, 6, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 2, -1, 4, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 2, 2, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 2}}) *- ----------------------------------------------------- ----------------------------------------------------- -- KKTsolutionVector -- Computes the solution x to the KKT equation KKTsolutionVector = {Simp=>false} >> o -> (N,S,I) -> ( A := KKTmatrix(N,S,I); a := avector(N,S,Simp=>o.Simp); b := transpose matrix({2/1*a}); x := solve(A,b,MaximalRank=>true); flatten entries x ); -* -- Examples KKTsolutionVector(3,{2,3},{}) KKTsolutionVector(15,{4,9},{7,8,10},Simp=>true) KKTsolutionVector(15,{4,9},{7,8,10}) -- Tests assert(KKTsolutionVector(3,{2,3},{})=={36/35, 6/5, 8/35, 8/5}) assert(KKTsolutionVector(15,{4,9},{7,8,10},Simp=>true)=={3847/3822, 5258/1911, 17396/5733, 12452/5733, 1210/637, 1520/1911, 5345/143913, 24/623623, 126/979, 26/979, 538/979, 784/979, 584/979, 128/979, 100/979, 2022/979}) assert(KKTsolutionVector(15,{4,9},{7,8,10})=={3847/3822, 5258/1911, 17396/5733, 12452/5733, 1210/637, 1520/1911, 1376/143913, 51621/623623, 72/979, -125/979, 727/979, 1427/979, 2012/979, 1192/979, 197/979, 1575/979}) *- ----------------------------------------------------- ----------------------------------------------------- -- Internal only: -- Compute the KKT matrix equation -- indexing the rows and columns starting at 1 -- to check the printed formulas newKKTmatrixEntry = (c,I,i,j) -> ( N:=#c-1; if j<=N-1 and not member(j,I) then ( if i==j then return c_j-c_(j+1); if i==j+1 then return c_(j+1)-c_(j-1); if i==j+2 then return c_(j-1)-c_j; return 0 ); if j<=N-1 and member(j,I) then ( if i<=j then return 2*(c_j-c_(i-1)) else return 0 ); if j==N then return 2*(c_N-c_(i-1)); if j==N+1 then return 2/1 ); newKKTmatrix = (c,I) -> ( N:=#c-1; matrix apply(toList(1..(N+1)), i -> apply(toList(1..(N+1)), j -> newKKTmatrixEntry(c,I,i,j))) ); newavec = (c) -> ( N:=#c-1; apply(toList(1..N+1), i -> if i==1 then c_1 else if i==N+1 then 2/1 else c_(i)-c_(i-2)) ); ----------------------------------------------------- ----------------------------------------------------- -- w -- Compute the stationary point on face_{gamma(I)} w = {Simp=>false} >> o -> (N,S,I) -> ( A := KKTmatrix(N,S,I); x := KKTsolutionVector(N,S,I,Simp=>o.Simp); w := sum apply(join(I,{#x-1,#x}), k -> (x_(k-1))*(1/2)*A_{k-1}); flatten entries w ); -* -- Examples w(3,{2,3},{}) w(15,{4,9},{7,8,10},Simp=>true) w(15,{4,9},{7,8,10}) -- Tests assert(w(3,{2,3},{})=={88/35, 72/35, 64/35, 8/5}) assert(w(15,{4,9},{7,8,10},Simp=>true)=={11491/1911, 10223/1911, 2985/637, 1234/273, 7687/1911, 7370/1911, 131/39, 2034/637, 3000/979, 2748/979, 2622/979, 2522/979, 2322/979, 202/89, 2122/979, 2022/979}) assert(w(15,{4,9},{7,8,10})=={11491/1911, 10223/1911, 2985/637, 1234/273, 7687/1911, 7370/1911, 131/39, 2034/637, 2973/979, 2829/979, 2757/979, 2560/979, 2166/979, 179/89, 1772/979, 1575/979}) *- ----------------------------------------------------- ----------------------------------------------------- -- approx -- Compute the decimal approximations to four places -- of a list of rational numbers approx = (w) -> ( apply(w, i -> round(4,toRR i)) ); -* -- Example approx {88/35,72/35,64/35,56/35} approx w(3,{2,3},{}) -- Test assert(approx {88/35,72/35,64/35,56/35}=={2.5143, 2.0571, 1.8286, 1.6}) *- ----------------------------------------------------- ----------------------------------------------------- -- Internal functions to interface with Octave vectToString = (v) -> ( s:=concatenate("[",toString(v_0)); for i from 1 to #v-1 do ( s = concatenate(s," ",toString(v_i)) ); s=concatenate(s,"]") ); matRowToString = (r) -> ( s:=toString(r_0); for i from 1 to #r-1 do ( s = concatenate(s," ",toString(r_i)) ); s ); matToString = (M) -> ( E:=entries M; s:=concatenate("[",matRowToString(E_0),";"); for i from 1 to #E-2 do ( s = concatenate(s," ",matRowToString(E_i),";") ); s = concatenate(s," ",matRowToString(last E),"]"); s ); computeI = (w,M) -> ( y:=flatten entries(M*(transpose matrix {w})); select(toList(1..(#w-2)), j -> abs(y_(j-1)) > 10^(-12)) ); parseOctaveNumberString = s -> ( for i from 0 to #s-2 do ( if s_i=="e" and s_(i+1)=="+" then ( t1=substring(s,0,i+1); t2=substring(s,i+2); return concatenate(t1,t2) ) ); return s ); ----------------------------------------------------- ----------------------------------------------------- -- proximumaW -- Returns two things: -- 1. prox(a,W) as a vector -- 2. the set of corners I -- This function first calls Octave to get a numerical approximation to w and -- an educated guess for I -- Then computes x and w using exact linear algebra over QQ in M2 -- Then returns w and I if they are correct proximumaW = {Simp => false} >> o -> (N,S) -> ( -- First, call Octave to get the candidate w and I a:=avector(N,S,Simp=>o.Simp); M:=matrix(Winequalities(gammaVectorOfSemigroup(N,S))); s0:=concatenate("N=",toString(N),";"); s1:="H = diag(ones(1,N+1));"; s2:=concatenate("f=",vectToString(-a),";"); s3:=concatenate("A=",matToString(M),";"); s4:="b = (0)*ones(1,N-1);"; octaveInput:=concatenate("warning('off','all')","\n","pkg load optim;","\n",s0,"\n",s1,"\n",s2,"\n",s3,"\n",s4,"\n","format long","\n","quadprog(H,f,A,b)"); octaveInputFN:=temporaryFileName(); octaveInputFile:=openOut octaveInputFN; octaveInputFile << octaveInput << endl; close octaveInputFile; s:=get concatenate(///!octave-cli < ///,octaveInputFN); t:= lines s; t=drop(t,{0,1}); t=drop(t,{#t-1,#t-1}); wOctave:=apply(t, i -> value parseOctaveNumberString i); I:=computeI(wOctave,M); run concatenate("rm ",octaveInputFN); -- Check that I is correct before returning these values x:=KKTsolutionVector(N,S,I,Simp=>o.Simp); if not all(#x-2, i -> x_i>=0) then error concatenate("The stationary point does not belong in W: I,x=",toString({I,x}))<< endl; if not all(I, i -> x_(i-1)>0) then error concatenate("This is not the smallest face containing the proximum: I,x=",toString({I,x}))<< endl; -- Now compute w with exact linear algebra -- and check that it approximately matches the Octave output wM2:=w(N,S,I,Simp=>o.Simp); if any(#wOctave, i -> abs(wOctave_i-wM2_i) > 10^-12) then error concatenate("wOctave does not approximate wM2 within the tolerance allowed: ",toString {wOctave,wM2}) << endl; return {wM2,I} ); -* -- Examples proximumaW(20,{2,3},Simp=>true) proximumaW(20,{2,3}) proximumaW(15,{4,9},Simp=>true) proximumaW(15,{4,9}) -- Tests assert(proximumaW(20,{2,3},Simp=>true) == {{33/14, 157/70, 153/70, 149/70, 29/14, 141/70, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2},{5,6}}) assert(proximumaW(20,{2,3}) == {{2777/1187, 2668/1187, 5227/2374, 2559/1187, 5009/2374, 42159/20179, 83483/40358, 41324/20179, 81813/40358, 40489/20179, 80143/40358, 39654/20179, 78473/40358, 38819/20179, 76803/40358, 32/17, 75133/40358, 37149/20179, 73463/40358, 36314/20179, 71793/40358},{4}}) assert(proximumaW(15,{4,9},Simp=>true) == {{11491/1911, 10223/1911, 2985/637, 1234/273, 7687/1911, 7370/1911, 131/39, 2034/637, 3000/979, 2748/979, 2622/979, 2522/979, 2322/979, 202/89, 2122/979, 2022/979},{7,8,10}}) assert(proximumaW(15,{4,9}) == {{3035219/507345, 2709481/507345, 794581/169115, 4604617/1014690, 411601/101469, 3953141/1014690, 1732267/507345, 220111/67646, 3519/1135, 472438/169115, 892983/338230, 84109/33823, 368652/169115, 685411/338230, 316759/169115, 116325/67646},{7}}) *- ----------------------------------------------------- ----------------------------------------------------- -- isPersistentCornerSet -- Tests whether a corner set I is the persistent solution isPersistentCornerSet = {Simp => false} >> o -> (S,I) -> ( if I=={} then return false; l:=last I; if l<=conductorIndex(S) then return false; N:=l+1; x:=KKTsolutionVector(N,S,I,Simp=>o.Simp); if o.Simp then ( if any(#x-2, k -> x_k < 0) then return false; if any(I, k -> x_(k-1)<=0) then return false; if x_(-2) != 0 then return false; if x_(-1) != 2 then return false; return true ); if not o.Simp then ( return (isPersistentCornerSet(S,I,Simp=>true) or isPersistentCornerSet(S,delete(l-1,I),Simp=>true)) ); ); -* -- Examples isPersistentCornerSet({2,3},{4,5},Simp=>true) isPersistentCornerSet({2,3},{5,6},Simp=>true) isPersistentCornerSet({2,3},{5,6}) isPersistentCornerSet({5,6,7,9},{5,6}) isPersistentCornerSet({5,6,7,9},{5,6},Simp=>true) isPersistentCornerSet({5,6,7,9},{6},Simp=>true) -- Tests assert(not isPersistentCornerSet({2,3},{4,5},Simp=>true)) assert(not isPersistentCornerSet({2,3},{4,5})) assert(isPersistentCornerSet({2,3},{5,6},Simp=>true)) assert(isPersistentCornerSet({2,3},{5,6})) assert(isPersistentCornerSet({4,9},{7,8,10,15,16},Simp=>true)) assert(isPersistentCornerSet({4,9},{7,8,10,15,16})) assert(not isPersistentCornerSet({4,9},{7,10},Simp=>true)) assert(not isPersistentCornerSet({4,9},{7,10})) assert(isPersistentCornerSet({5,6,7,9},{5,6})) assert(not isPersistentCornerSet({5,6,7,9},{5,6},Simp=>true)) assert(isPersistentCornerSet({5,6,7,9},{6},Simp=>true)) *- ----------------------------------------------------- ----------------------------------------------------- -- Internal function used in persistentCornerSet persistentCornerSetSimp = S -> ( ci := conductorIndex(S); N:=max(ci,4); Isimp:={}; while Ntrue); if isPersistentCornerSet(S,Isimp,Simp=>true) then return Isimp else N=N+1; ); error "No Isimp found with N false} >> o -> (S) -> ( Isimp:=persistentCornerSetSimp(S); if o.Simp then return Isimp; l:=max Isimp; if member(l-1,Isimp) then return Isimp else return insert(-2,l-1,Isimp) ); -* -- Examples persistentCornerSet({2,3}) persistentCornerSet({2,5}) persistentCornerSet({4,9}) persistentCornerSet({4,9},Simp=>true) persistentCornerSet({5,7}) persistentCornerSet({5,6,7,9}) persistentCornerSet({5,6,7,9},Simp=>true) persistentCornerSet({8,13}) -- Tests assert(persistentCornerSet({2,3})=={5, 6}) assert(persistentCornerSet({2,5})=={6, 7}) assert(persistentCornerSet({4,9})=={7, 8, 10, 15, 16}) assert(persistentCornerSet({4,9},Simp=>true)=={7, 8, 10, 15, 16}) assert(persistentCornerSet({5,7})=={9, 16, 17}) assert(persistentCornerSet({5,6,7,9})=={5, 6}) assert(persistentCornerSet({5,6,7,9},Simp=>true)=={6}) assert(persistentCornerSet({8,13})=={6, 11, 26, 38, 49, 50}) -- Test this code against the output from my old code load "March2024output.m2" time for L in March2024output do ( print toString({L_0,persistentCornerSet(L_0,Simp=>true)==L_3})< ( Isimp:=persistentCornerSet(S,Simp=>true); l:=max Isimp; x:=KKTsolutionVector(l+2,S,Isimp,Simp=>true); N0:=floor(l+2/x_(l-1)); Nstart:=max(conductorIndex(S)+1,4); I:=last proximumaW(Nstart,S); Ns = {Nstart}; Is = {I}; print concatenate("N=",toString(Nstart), ", I = ",toString(I)) << endl; for N from Nstart to N0 do ( I = last proximumaW(N,S); if I!=last Is then ( print concatenate("N=",toString(N), ", I = ",toString(I)) << endl; Ns = append(Ns,N); Is = append(Is,I) ) ); returnValue:=apply(#Ns-1, i -> {Ns_i,Ns_(i+1)-1,Is_i}); returnValue = append(returnValue,{last Ns,infinity,last Is}); returnValue ); -* -- Examples worst1PSAsNVaries({2,3}) -- Tests assert(worst1PSAsNVaries({2,3})=={{4, 15, {}}, {16, 28, {4}}, {29, 74, {4, 5}}, {75, 145, {5}}, {146, infty, {5, 6}}}) *- ----------------------------------------------------- ----------------------------------------------------- ----------------------------------------------------- -- Functions for analyzing solutions to the KKT equation -- as rational functions ----------------------------------------------------- ----------------------------------------------------- ----------------------------------------------------- -- Internal function subCol = (j,b,M) -> ( if numrows(b)!=numrows(M) then error "b and M have different numbers of rows"; if ring(b) =!= ring(M) then error "b and M have different base rings"; E1:=entries transpose M; b2:=flatten entries b; E2:=apply(#E1, i -> if i+1==j then b2 else E1_i); transpose matrix E2 ); ----------------------------------------------------- ----------------------------------------------------- -- Internal function KKTA = {Simp => false} >> o -> (N,j,S,I) -> ( A := KKTmatrix(N,S,I); a := avector(N,S,Simp=>o.Simp); b := transpose matrix({apply(a, i -> 2/1*i)}); subCol(j,b,A) ) ----------------------------------------------------- ----------------------------------------------------- -- Aprime -- Replace column j in the KKT matrix A with the RHS -- as arises in Cramer's Rule -- Here we number the columns starting at j=1 Aprime = {Simp => false} >> o -> (N,j,S,I) -> ( A := KKTmatrix(N,S,I); a := avector(N,S,Simp=>o.Simp); b := transpose matrix({apply(a, i -> 2*i-4/1)}); subCol(j,b,A) ) -* -- Example Aprime(15,12,{4,9},{7,8,10},Simp=>true) -- Test assert(Aprime(15,12,{4,9},{7,8,10},Simp=>true)==matrix {{-4, 0, 0, 0, 0, 0, 34, 36, 0, 42, 0, 4, 0, 0, 54, 2}, {8, -1, 0, 0, 0, 0, 26, 28, 0, 34, 0, 12, 0, 0, 46, 2}, {-4, 5, -3, 0, 0, 0, 18, 20, 0, 26, 0, 6, 0, 0, 38, 2}, {0, -4, 4, -1, 0, 0, 16, 18, 0, 24, 0, 4, 0, 0, 36, 2}, {0, 0, -1, 4, -3, 0, 10, 12, 0, 18, 0, 4, 0, 0, 30, 2}, {0, 0, 0, -3, 4, -1, 8, 10, 0, 16, 0, 4, 0, 0, 28, 2}, {0, 0, 0, 0, -1, 4, 2, 4, 0, 10, 0, 4, 0, 0, 22, 2}, {0, 0, 0, 0, 0, -3, 0, 2, 0, 8, 0, 0, 0, 0, 20, 2}, {0, 0, 0, 0, 0, 0, 0, 0, -1, 6, 0, 2, 0, 0, 18, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 3, 2, 0, 2, 0, 0, 14, 2}, {0, 0, 0, 0, 0, 0, 0, 0, -2, 0, -2, 0, 0, 0, 12, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 2, 0, 0, 10, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 2, -1, 0, 6, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, -1, 4, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 2, 2, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 2/1}}) *- ----------------------------------------------------- ----------------------------------------------------- --polynomialInterpolation -- Given a list of x values and a list of y values, -- compute the polynomial interpolating the points (x_i,y_i) -- Also enter the interpolation ring IR polynomialInterpolation = (x,y,IR) -> ( vdM:=matrix apply(x, k-> apply(#x, l -> k^(#x-l-1)/1)); N:=vdM^(-1); O:=N*(transpose matrix {y}); P:=matrix {apply(#x, l -> (IR_0)^(#x-l-1))}; (P*O)_(0,0) ); -* -- Example xvalues = {-5,-4,-3,-2,-1,0,1,2,3,4,5}; yvalues = {-345/2, -93, -83/2, -12, 3/2, 5, 9/2, 6, 31/2, 39, 165/2}; polynomialInterpolation(xvalues,yvalues,QQ[x]) -- Test xvalues = {-5,-4,-3,-2,-1,0,1,2,3,4,5}; yvalues = {-345/2, -93, -83/2, -12, 3/2, 5, 9/2, 6, 31/2, 39, 165/2}; f=polynomialInterpolation(xvalues,yvalues,QQ[x]); R = ring(f); assert(f==(R_0)^3-2*(R_0)^2+(1/2)*(R_0)+5) *- ----------------------------------------------------- ----------------------------------------------------- -- Internal function matrixRecurrencePhiEntries = (M,i,j) -> ( n:=numRows(M); -- Keep the upper left nx(n-2) block if i < n and j ( n:=numRows(M); matrix apply(n+1, i -> apply(n+1, j -> matrixRecurrencePhiEntries(M,i,j))) ); -* -- Examples S={4,9} I={7,8,10} N=15 matrixRecurrencePhi(KKTmatrix(N,S,I))==KKTmatrix(N+1,S,I) -- Tests -- Here, M1 = KKTmatrix(15,{4,9},{7,8,10}), and M2 = KKTmatrix(16,{4,9},{7,8,10}) M1 = matrix {{-4, 0, 0, 0, 0, 0, 34, 36, 0, 42, 0, 0, 0, 0, 54, 2}, {8, -1, 0, 0, 0, 0, 26, 28, 0, 34, 0, 0, 0, 0, 46, 2}, {-4, 5, -3, 0, 0, 0, 18, 20, 0, 26, 0, 0, 0, 0, 38, 2}, {0, -4, 4, -1, 0, 0, 16, 18, 0, 24, 0, 0, 0, 0, 36, 2}, {0, 0, -1, 4, -3, 0, 10, 12, 0, 18, 0, 0, 0, 0, 30, 2}, {0, 0, 0, -3, 4, -1, 8, 10, 0, 16, 0, 0, 0, 0, 28, 2}, {0, 0, 0, 0, -1, 4, 2, 4, 0, 10, 0, 0, 0, 0, 22, 2}, {0, 0, 0, 0, 0, -3, 0, 2, 0, 8, 0, 0, 0, 0, 20, 2}, {0, 0, 0, 0, 0, 0, 0, 0, -1, 6, 0, 0, 0, 0, 18, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 3, 2, 0, 0, 0, 0, 14, 2}, {0, 0, 0, 0, 0, 0, 0, 0, -2, 0, -2, 0, 0, 0, 12, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, -1, 0, 0, 10, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 3, -1, 0, 6, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 2, -1, 4, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 2, 2, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 2/1}}; M2 = matrix {{-4, 0, 0, 0, 0, 0, 34, 36, 0, 42, 0, 0, 0, 0, 0, 56, 2}, {8, -1, 0, 0, 0, 0, 26, 28, 0, 34, 0, 0, 0, 0, 0, 48, 2}, {-4, 5, -3, 0, 0, 0, 18, 20, 0, 26, 0, 0, 0, 0, 0, 40, 2}, {0, -4, 4, -1, 0, 0, 16, 18, 0, 24, 0, 0, 0, 0, 0, 38, 2}, {0, 0, -1, 4, -3, 0, 10, 12, 0, 18, 0, 0, 0, 0, 0, 32, 2}, {0, 0, 0, -3, 4, -1, 8, 10, 0, 16, 0, 0, 0, 0, 0, 30, 2}, {0, 0, 0, 0, -1, 4, 2, 4, 0, 10, 0, 0, 0, 0, 0, 24, 2}, {0, 0, 0, 0, 0, -3, 0, 2, 0, 8, 0, 0, 0, 0, 0, 22, 2}, {0, 0, 0, 0, 0, 0, 0, 0, -1, 6, 0, 0, 0, 0, 0, 20, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 3, 2, 0, 0, 0, 0, 0, 16, 2}, {0, 0, 0, 0, 0, 0, 0, 0, -2, 0, -2, 0, 0, 0, 0, 14, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, -1, 0, 0, 0, 12, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 3, -1, 0, 0, 8, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 2, -1, 0, 6, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 2, -1, 4, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 2, 2, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 2/1}}; assert(matrixRecurrencePhi(M1)==M2) *- ----------------------------------------------------- ----------------------------------------------------- -- Internal function -- To get A'(N+1,N,S,I) from A'(N,N-1,S,I): -- Insert -1 2 -1 as column N-2 -- Column N-1 is the same, insert 0 at bottom -- Column N increases by 2, insert 0 at bottom -- Column N+1 is the same, insert 2 at bottom matrixRecurrenceChiEntries = (M,i,j) -> ( n:=numRows(M); -- Keep the upper left nx(n-3) block if i < n and j ( n:=numRows(M); matrix apply(n+1, i -> apply(n+1, j -> matrixRecurrenceChiEntries(M,i,j))) ); -* -- Examples S = {4,9}; I = {7,8,10}; N = 15; matrixRecurrenceChi(Aprime(N,N-1,S,I,Simp=>true))==Aprime(N+1,N,S,I,Simp=>true) -- Tests M1=matrix {{-4, 0, 0, 0, 0, 0, 34, 36, 0, 42, 0, 0, 0, 4, 54, 2}, {8, -1, 0, 0, 0, 0, 26, 28, 0, 34, 0, 0, 0, 12, 46, 2}, {-4, 5, -3, 0, 0, 0, 18, 20, 0, 26, 0, 0, 0, 6, 38, 2}, {0, -4, 4, -1, 0, 0, 16, 18, 0, 24, 0, 0, 0, 4, 36, 2}, {0, 0, -1, 4, -3, 0, 10, 12, 0, 18, 0, 0, 0, 4, 30, 2}, {0, 0, 0, -3, 4, -1, 8, 10, 0, 16, 0, 0, 0, 4, 28, 2}, {0, 0, 0, 0, -1, 4, 2, 4, 0, 10, 0, 0, 0, 4, 22, 2}, {0, 0, 0, 0, 0, -3, 0, 2, 0, 8, 0, 0, 0, 0, 20, 2}, {0, 0, 0, 0, 0, 0, 0, 0, -1, 6, 0, 0, 0, 2, 18, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 3, 2, 0, 0, 0, 2, 14, 2}, {0, 0, 0, 0, 0, 0, 0, 0, -2, 0, -2, 0, 0, 0, 12, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, -1, 0, 2, 10, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 3, -1, 2, 6, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 2, 0, 4, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 2, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2/1}}; M2=matrix {{-4, 0, 0, 0, 0, 0, 34, 36, 0, 42, 0, 0, 0, 0, 4, 56, 2}, {8, -1, 0, 0, 0, 0, 26, 28, 0, 34, 0, 0, 0, 0, 12, 48, 2}, {-4, 5, -3, 0, 0, 0, 18, 20, 0, 26, 0, 0, 0, 0, 6, 40, 2}, {0, -4, 4, -1, 0, 0, 16, 18, 0, 24, 0, 0, 0, 0, 4, 38, 2}, {0, 0, -1, 4, -3, 0, 10, 12, 0, 18, 0, 0, 0, 0, 4, 32, 2}, {0, 0, 0, -3, 4, -1, 8, 10, 0, 16, 0, 0, 0, 0, 4, 30, 2}, {0, 0, 0, 0, -1, 4, 2, 4, 0, 10, 0, 0, 0, 0, 4, 24, 2}, {0, 0, 0, 0, 0, -3, 0, 2, 0, 8, 0, 0, 0, 0, 0, 22, 2}, {0, 0, 0, 0, 0, 0, 0, 0, -1, 6, 0, 0, 0, 0, 2, 20, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 3, 2, 0, 0, 0, 0, 2, 16, 2}, {0, 0, 0, 0, 0, 0, 0, 0, -2, 0, -2, 0, 0, 0, 0, 14, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, -1, 0, 0, 2, 12, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 3, -1, 0, 2, 8, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 2, -1, 0, 6, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 2, 0, 4, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 2, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2/1}}; assert(matrixRecurrenceChi(M1)==M2) *- ----------------------------------------------------- ----------------------------------------------------- -- Internal function matrixRecurrencePsiEntries = (M,i,j) -> ( n:=numRows(M); -- Keep the upper left nx(n-2) block if i < n and j ( n:=numRows(M); matrix apply(n+1, i -> apply(n+1, j -> matrixRecurrencePsiEntries(M,i,j))) ); -* -- Examples S = {4,9}; I = {7,8,10}; N = 15; matrixRecurrencePsi(Aprime(N,N,S,I,Simp=>true))==Aprime(N+1,N+1,S,I,Simp=>true) -- Tests M1=matrix {{-4, 0, 0, 0, 0, 0, 34, 36, 0, 42, 0, 0, 0, 0, 4, 2}, {8, -1, 0, 0, 0, 0, 26, 28, 0, 34, 0, 0, 0, 0, 12, 2}, {-4, 5, -3, 0, 0, 0, 18, 20, 0, 26, 0, 0, 0, 0, 6, 2}, {0, -4, 4, -1, 0, 0, 16, 18, 0, 24, 0, 0, 0, 0, 4, 2}, {0, 0, -1, 4, -3, 0, 10, 12, 0, 18, 0, 0, 0, 0, 4, 2}, {0, 0, 0, -3, 4, -1, 8, 10, 0, 16, 0, 0, 0, 0, 4, 2}, {0, 0, 0, 0, -1, 4, 2, 4, 0, 10, 0, 0, 0, 0, 4, 2}, {0, 0, 0, 0, 0, -3, 0, 2, 0, 8, 0, 0, 0, 0, 0, 2}, {0, 0, 0, 0, 0, 0, 0, 0, -1, 6, 0, 0, 0, 0, 2, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 3, 2, 0, 0, 0, 0, 2, 2}, {0, 0, 0, 0, 0, 0, 0, 0, -2, 0, -2, 0, 0, 0, 0, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, -1, 0, 0, 2, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 3, -1, 0, 2, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 2, -1, 0, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 2, 0, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 2/1}}; M2=matrix {{-4, 0, 0, 0, 0, 0, 34, 36, 0, 42, 0, 0, 0, 0, 0, 4, 2}, {8, -1, 0, 0, 0, 0, 26, 28, 0, 34, 0, 0, 0, 0, 0, 12, 2}, {-4, 5, -3, 0, 0, 0, 18, 20, 0, 26, 0, 0, 0, 0, 0, 6, 2}, {0, -4, 4, -1, 0, 0, 16, 18, 0, 24, 0, 0, 0, 0, 0, 4, 2}, {0, 0, -1, 4, -3, 0, 10, 12, 0, 18, 0, 0, 0, 0, 0, 4, 2}, {0, 0, 0, -3, 4, -1, 8, 10, 0, 16, 0, 0, 0, 0, 0, 4, 2}, {0, 0, 0, 0, -1, 4, 2, 4, 0, 10, 0, 0, 0, 0, 0, 4, 2}, {0, 0, 0, 0, 0, -3, 0, 2, 0, 8, 0, 0, 0, 0, 0, 0, 2}, {0, 0, 0, 0, 0, 0, 0, 0, -1, 6, 0, 0, 0, 0, 0, 2, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 3, 2, 0, 0, 0, 0, 0, 2, 2}, {0, 0, 0, 0, 0, 0, 0, 0, -2, 0, -2, 0, 0, 0, 0, 0, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, -1, 0, 0, 0, 2, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 3, -1, 0, 0, 2, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 2, -1, 0, 0, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 2, -1, 0, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 2, 0, 2}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 2/1}}; assert(matrixRecurrencePsi(M1)==M2) *- ----------------------------------------------------- ----------------------------------------------------- -- taylor -- Compute the Taylor expansion of a polynomial f at the center c -- Returns the list of coefficients {a_n,...,a_0} in order of descending degree -- That is, f = a_n*(x-c)^n + a_(n-1)*(x-c)^(n-1) + ... taylor = (f,c) -> ( g:=f; j:=(ring g)_0; t:={substitute(g,{j=>c})}; for k from 1 to first(degree(g)) do ( g = diff(j,g); t = prepend(1/(k!)*substitute(g,{j=>c}),t) ); return t ); -* -- Examples R=QQ[x]; f = 4*x^4-30*x^3+77*x^2-85*x+48; taylor(f,2) taylor(f,3) -- Tests R=QQ[x]; f = 4*x^4-30*x^3+77*x^2-85*x+48; assert(taylor(f,2)=={4, 2, -7, -9, 10}) assert(taylor(f,3)=={4, 18, 23, -1, 0}) *- ----------------------------------------------------- ----------------------------------------------------- -- printTaylor -- Prints the Taylor expansion of f at c printTaylor = (f,c,var) -> ( t:=taylor(f,c); n:=#t-1; s:=concatenate(toString(t_0),"*(",var,"-",toString(c),")^",toString(n)); for k from 1 to n do ( if t_k==0 then continue; if t_k>0 then s=concatenate(s,"+"); s=concatenate(s,toString(t_k),"*(",var,"-",toString(c),")^",toString(n-k)); ); return s ); -* -- Examples R=QQ[x]; f = 4*x^4-30*x^3+77*x^2-85*x+48; printTaylor(f,2,"x") printTaylor(f,3,"x") -- Tests R=QQ[x]; f = 4*x^4-30*x^3+77*x^2-85*x+48; assert(printTaylor(f,2,"x")=="4*(x-2)^4+2*(x-2)^3-7*(x-2)^2-9*(x-2)^1+10*(x-2)^0") assert(printTaylor(f,3,"x")=="4*(x-3)^4+18*(x-3)^3+23*(x-3)^2-1*(x-3)^1") *- ----------------------------------------------------- ----------------------------------------------------- -- Del -- Delete a subset of the rows I and a subset of the columns J Del = (I,J,M) -> ( n:=numRows(M); s:=apply(n, i -> i); rList:=select(s,i -> not member(i+1,I)); cList:=select(s,j -> not member(j+1,J)); M_cList^rList ); -* -- Examples M = matrix apply(9, i-> apply(9, j -> 10*(i+1)+j+1)) Del({2,5},{1,3},M) -- Tests M = matrix apply(9, i-> apply(9, j -> 10*(i+1)+j+1)); assert(Del({2,5},{1,3},M)==matrix {{12, 14, 15, 16, 17, 18, 19}, {32, 34, 35, 36, 37, 38, 39}, {42, 44, 45, 46, 47, 48, 49}, {62, 64, 65, 66, 67, 68, 69}, {72, 74, 75, 76, 77, 78, 79}, {82, 84, 85, 86, 87, 88, 89}, {92, 94, 95, 96, 97, 98, 99}}) *- ----------------------------------------------------- ----------------------------------------------------- -- del -- Compute the determinant of Del(I,J,M) del = (I,J,M) -> ( det Del(I,J,M) ); -* -- Examples M = matrix apply(9, i-> apply(9, j -> 10*(i+1)+j+1)) del({2,5},{1,3},M) -- Tests M = matrix apply(9, i-> apply(9, j -> 10*(i+1)+j+1)); assert(del({2,5},{1,3},M)==0) *- ----------------------------------------------------- ----------------------------------------------------- -- Internal functions h4 = (M) -> ( n:=numrows(M); return 1/3*(-1)^n*del({n-1,n},{n-1,n},M) ); h3 = (M) -> ( n:=numrows(M); s:=0; s = s + del({n-1},{n-1},M); s = s + del({n},{n-1},M); return 2/3*(-1)^n*s ); h2 = (M) -> ( n:=numrows(M); s:=0; s = s + 3*del({n},{n},M); s = s + 3*del({n-1},{n},M); s = s + 3*del({n-1},{n-1},M); s = s + 6*del({n},{n-1},M); s = s - del({n-1,n},{n-1,n},M); return 1/3*(-1)^n*s ); h1 = (M) -> ( n:=numrows(M); s:=0; s = s + 9*del({n},{n},M); s = s + 3*del({n-1},{n},M); s = s + del({n-1},{n-1},M); s = s + 4*del({n},{n-1},M); return 1/3*(-1)^n*s ); h0 = M -> ( n:=numrows(M); return (-1)^n*det(M) ); ----------------------------------------------------- ----------------------------------------------------- -- hCoefficients -- Compute the coefficients of the polynomial -- h(N) = h_4*(x-(n+1))^4 + h_3*(x-(n+1))^3 + h_2*(x-(n+1))^2 + h_1*(x-(n+1)) + h_0 -- giving h(N) = (-1)^N det(M(N)) where M(N) is given -- by the recurrence Phi starting at M(n), where M(n) is nxn hCoefficients = (M) -> ( {h4(M),h3(M),h2(M),h1(M),h0(M)} ); -* -- Examples M = KKTmatrix(14,{4,9},{7,8,10}) hCoefficients(M) hCoefficients(KKTmatrix(16,{4,9},{7,8,10,15})) -- Tests assert(hCoefficients(KKTmatrix(14,{4,9},{7,8,10}))=={-3669120, -79252992, -731622528, -3121687296, -4684732416}) assert(hCoefficients(KKTmatrix(16,{4,9},{7,8,10,15}))=={656038656, 8371464192, 29143086336, 38669589504, 17241928704}) *- ----------------------------------------------------- ----------------------------------------------------- -- Internal functions chi3 = (M) -> ( n:=numrows(M); s:=0; s = 2*del({n-1,n},{n-1,n},M)+2*del({n-2,n},{n-1,n},M)+2*del({n-2,n-1},{n-1,n},M); return 1/3*(-1)^n*s ); chi2 = (M) -> ( n:=numrows(M); s:=0; s = 4*del({n-1,n},{n-1,n},M)+2*del({n-2,n},{n-1,n},M); return (-1)^n*s ); chi1 = (M) -> ( n:=numrows(M); s:=0; s = 8*del({n-2,n},{n-1,n},M)+14*del({n-1,n},{n-1,n},M)+6*del({n+1},{n},matrixRecurrenceChi M); return 1/3*(-1)^(n+1)*s ); chi0 = M -> ( n:=numrows(M); return (-1)^n*det(M) ); ----------------------------------------------------- ----------------------------------------------------- -- Compute the coefficients of the polynomial -- chi(N) = chi_3*(x-c)^3 + chi_2*(x-c)^2 + chi_1*(x-c) + chi_0 -- giving chi(N) = (-1)^N det(M(N)) where M(N) is given -- by the recurrence Chi chiCoefficients = (M) -> ( {chi3(M),chi2(M),chi1(M),chi0(M)} ); -* -- Examples M = Aprime(14,13,{4,9},{7,8,10},Simp=>true) chiCoefficients(M) chiCoefficients(Aprime(17,16,{4,9},{7,8,10,15},Simp=>true)) -- Tests assert(chiCoefficients(Aprime(14,13,{4,9},{7,8,10},Simp=>true))=={35223552, 352235520, 457906176, -1972518912}) assert(chiCoefficients(Aprime(17,16,{4,9},{7,8,10,15},Simp=>true))=={-375717888, -2254307328, -4132896768, -2254307328}) *- ----------------------------------------------------- ----------------------------------------------------- -- Internal functions psi2 = (M) -> ( n:=numrows(M); s:=0; s = s + del({n},{n},M); s = s + del({n-1},{n},M); return (-1)^n*s ); psi1 = (M) -> ( n:=numrows(M); s:=0; s = s + 3*del({n},{n},M); s = s + del({n-1},{n},M); return (-1)^n*s ); psi0 = M -> ( n:=numrows(M); return (-1)^n*det(M) ); ----------------------------------------------------- ----------------------------------------------------- -- Compute the coefficients of the polynomial -- psi(N) = psi_2*(x-c)^2 + psi_1*(x-c) + psi_0 -- giving chi(N) = (-1)^N det(M(N)) where M(N) is given -- by the recurrence Psi psiCoefficients = (M) -> ( {psi2(M),psi1(M),psi0(M)} ); -* -- Examples M = Aprime(14,14,{4,9},{7,8,10},Simp=>true) psiCoefficients(M) psiCoefficients(Aprime(17,17,{4,9},{7,8,10,15},Simp=>true)) -- Tests assert(psiCoefficients(Aprime(14,14,{4,9},{7,8,10},Simp=>true))=={-52835328, -405070848, -422682624}) assert(psiCoefficients(Aprime(17,17,{4,9},{7,8,10,15},Simp=>true))=={563576832, 2817884160, 3381460992}) *- ----------------------------------------------------- ----------------------------------------------------- Q = (S,I,IR) -> ( ci := conductorIndex(S); N0 := max({ci+2,max(I)+2}); xs:=apply(10, i -> N0+5+i); -- Denom ys:=apply(xs, N -> (-1)^(N+1)*det(KKTmatrix(N,S,I))); f:=polynomialInterpolation(xs,ys,IR); print concatenate("Q = ",printTaylor(f,N0,toString(IR_0)))< ( ci := conductorIndex(S); N0 := max({ci+2,max(I)+2,j+1}); xs:=apply(10, i -> N0+5+i); -- Denom ys:=apply(xs, N -> (-1)^(N+1)*det(Aprime(N,j,S,I,Simp=>true))); f:=polynomialInterpolation(xs,ys,IR); print concatenate("P_",toString(j)," = ",printTaylor(f,N0,toString(IR_0)))<true) substitute(p3,{N=>20})/substitute(q,{N=>20}) -- Tests *- ----------------------------------------------------- ----------------------------------------------------- chi(List,List,Ring) := (S,I,IR) -> ( ci := conductorIndex(S); N0 := max({ci+2,max(I)+2}); xs:=apply(10, i -> N0+5+i); ys:=apply(xs, N -> (-1)^(N+1)*det(Aprime(N,N-1,S,I,Simp=>true))); f:=polynomialInterpolation(xs,ys,IR); print concatenate("chi = ",printTaylor(f,N0,toString(IR_0)))<true) substitute(f,{N=>20})/substitute(q,{N=>20}) -- Tests *- ----------------------------------------------------- ----------------------------------------------------- psi = (S,I,IR) -> ( ci := conductorIndex(S); N0 := max({ci+2,max(I)+2}); xs:=apply(10, i -> N0+5+i); ys:=apply(xs, N -> (-1)^(N+1)*det(Aprime(N,N,S,I,Simp=>true))); f:=polynomialInterpolation(xs,ys,IR); print concatenate("psi = ",printTaylor(f,N0,toString(IR_0)))<true) substitute(f,{N=>20})/substitute(q,{N=>20}) -- Tests *- ----------------------------------------------------- ----------------------------------------------------- omega = (S,I,IR) -> ( ci := conductorIndex(S); N0 := max({ci+1,max(I)+1}); xs:=apply(10, i -> N0+5+i); ys:=apply(xs, N -> (-1)^(N+1)*det(KKTA(N,N+1,S,I,Simp=>true))); f:=polynomialInterpolation(xs,ys,IR); print concatenate("omega = ",printTaylor(f,N0,toString(IR_0)))<true) substitute(f,{N=>20})/substitute(q,{N=>20}) -- Tests *- ----------------------------------------------------- interpolateXs = (S,I,IR) -> ( ci := conductorIndex(S); N0 := max({ci+1,max(I)+1}); xs:=apply(10, i -> N0+5+i); -- Denom ys:=apply(xs, N -> (-1)^(N+1)*det(KKTmatrix(N,S,I))); f:=polynomialInterpolation(xs,ys,IR); print concatenate("Denominator: ",printTaylor(f,N0,toString(IR_0)))< (-1)^(N+1)*det(Aprime(N,N-1,S,I,Simp=>true))); f=polynomialInterpolation(xs,ys,IR); print concatenate("xNm1 num: ",printTaylor(f,N0,toString(IR_0)))< (-1)^(N+1)*det(Aprime(N,N,S,I,Simp=>true))); f=polynomialInterpolation(xs,ys,IR); print concatenate("xN num: ",printTaylor(f,N0,toString(IR_0)))< (-1)^(N+1)*det(KKTA(N,N+1,S,I,Simp=>true))); f=polynomialInterpolation(xs,ys,IR); print concatenate("xNp1 num: ",printTaylor(f,N0,toString(IR_0)))<true) KKTsolutionVector(3,{2,3},{},Simp=>true) -- Example 2c KKTmatrix(14,{4,9},{7,10}) avector(14,{4,9},Simp=>true) KKTsolutionVector(14,{4,9},{7,10},Simp=>true) -- Example 3 load "Worst1PS.m2" proximumaW(3,{2,3}) -- Example 4 load "Worst1PS.m2" ci = conductorIndex({4,9}) x = KKTsolutionVector(16,{4,9},{7,8,10},Simp=>true) R=QQ[N] q=Q({4,9},{7,8,10},R) p3=Pj(3,{4,9},{7,8,10},R) substitute(p3,{N=>16})/substitute(q,{N=>16}) f1=chi({4,9},{7,8,10},R) substitute(f1,{N=>16})/substitute(q,{N=>16}) f2=psi({4,9},{7,8,10},R) substitute(f2,{N=>16})/substitute(q,{N=>16}) f3=omega({4,9},{7,8,10},R) substitute(f3,{N=>16})/substitute(q,{N=>16}) -- Check the coefficients of chi k=14 M = Aprime(k,k-1,{4,9},{7,8,10},Simp=>true) 1/3*(-1)^(k+1)*(2*del({k,k+1},{k,k+1},M)+2*del({k-1,k+1},{k,k+1},M)) (-1)^(k+1)*(4*del({k,k+1},{k,k+1},M)+2*del({k-1,k+1},{k,k+1},M)) 1/3*(-1)^(k+1)*(-14*del({k,k+1},{k,k+1},M)-8*del({k-1,k+1},{k,k+1},M)-6*del({k+2},{k+1},Aprime(k+1,k,{4,9},{7,8,10},Simp=>true))) (-1)^(k+1)*det(M) -- Check the coefficients of psi -52835328*(N-14)^2-405070848*(N-14)^1-422682624*(N-14)^0 k=14 M = Aprime(k,k,{4,9},{7,8,10},Simp=>true) (-1)^(k+1)*(del({k+1},{k+1},M)+del({k},{k+1},M)) (-1)^(k+1)*(3*del({k+1},{k+1},M)+del({k},{k+1},M)) (-1)^(k+1)*det(M) -- Check Lemma 3.12 (1/2)*leadCoefficient(f1)==-(1/3)*leadCoefficient(f2) -- Check Lemma 3.13 p14=Pj(14,{4,9},{7,8,10},R) p14 == (N-14)*(N-14+1)*(1/2*f1+1/3*(N-14-1)*f2) -- Do a different example to check Proposition 5.5 -- Check f1=chi({4,9},{7,8,10,20},R) factor(f1) f2=psi({4,9},{7,8,10,20},R) factor(f2)