chineseRemainder=(a, m) -> ( M := 1; N := 0; xtmp := 0; apply(m, c-> M= M*c); apply(0..(length m - 1), i -> ( N = sub(M/m_i,ZZ); xtmp = xtmp + a_i*(gcdCoefficients(m_i,N))_2*N; )); (xtmp%M,M)); ratConvert=(c, m) -> ( a1 := m; a2 := c; u1 := 0; u2 := 1; q := 0; h := 0; ret:= 0; while ( (u2^2 < m/2) and (a2^2 >= m/2)) do ( q = a1//a2; a1 = a1 - q*a2; u1 = u1 - q*u2; h = a1; a1 = a2; a2 =h; h = u1; u1 = u2; u2 =h; ); if (u2^2 >= m/2) then (ret = {0};) else ( if (a2^2 < m/2) then (ret = {a2,u2};)); ret) solveQuadratic=(eqtn, varx) -> ( B := contract(varx,eqtn - varx^2); C := eqtn - (varx^2 + B*varx); (sub(- B/2,RR) + sqrt(sub((B/2)^2 - C,RR)), sub(- B/2,RR) - sqrt(sub((B/2)^2 - C,RR)))) debug PrimaryDecomposition; removeIsolatedPoints = method(TypicalValue => Ideal) removeIsolatedPoints(Ideal):=(J) -> ( mInd := support (mInd0 = (independentSets(J))_0); C := flatt(J,mInd0); f := (saturate(intersect(values C),t_(ring J)))_0; facsf := factors f; facs := select(facsf, g -> first degree g < 4); thisdoesit := (#facs == #facsf); Jtemp := J; scan(facs, fac -> Jtemp = saturate(Jtemp,fac)); if not thisdoesit then removeIsolatedPoints(Jtemp,4) else Jtemp) removeIsolatedPoints(Ideal,ZZ):=(J,m) -> ( Jtemp:=J; J1:=intersect(apply(m,c->( counter = 0; while (Jtemp = saturate(J+ideal random(1, ring J)); not degree Jtemp == degree J) do (); Jtemp))); degs:=tally (degrees gens J1)_1; J2:=ideal (gens J1)_{0..degs_{1}+degs_{2}-1}; J2) setupSE3=method(TypicalValue => Ring) setupSE3(Ring):=(R)->( use R; A := matrix apply(3,i->apply(3,j->a_(i,j))); O := ideal mingens ideal (transpose A *A - t^2*id_(R^3)) +ideal (det A-t^3); SE3 = saturate(O,ideal t)) setupEquation=method() setupEquation(Ring):=(RPQ) -> ( use RPQ; A := matrix apply(3,i -> apply(3,j->a_(i,j))); B := matrix{{b_0},{b_1},{b_2}}; P := matrix{{p_0},{p_1},{p_2}}; Q := matrix{{q_0},{q_1},{q_2}}; eq=t^2*scalarProduct(P-Q) - scalarProduct(A*P+B-t*Q)) scalarProduct = method () scalarProduct(Matrix,Matrix):=(X,Y) -> sum(3, i -> X^{i}*Y^{i}) scalarProduct(Matrix):=(X) -> sum(3, i -> X^{i}^2) eulerAngles = method() eulerAngles(Ring,Ring):= (S,R)-> ( r := (entries vars S)_0; a := (entries vars R)_0; study={a_0=> r_0^2 - r_1^2 - r_2^2 + r_3^2, a_1=> 2*(r_0*r_1 - r_2*r_3), a_2=> 2*(r_0*r_2 + r_1*r_3), a_3=> 2*(r_0*r_1 + r_2*r_3), a_4=> -r_0^2 + r_1^2 - r_2^2 + r_3^2, a_5=> 2*(r_1*r_2 - r_0*r_3), a_6=> 2*(r_0*r_2 - r_1*r_3), a_7=> 2*(r_1*r_2 + r_0*r_3), a_8=> - r_0^2 - r_1^2 + r_2^2 + r_3^2, a_9=> r_0^2 + r_1^2 + r_2^2 + r_3^2}| apply(dim R - 11, i-> a_(10 + i) => r_(4 + i)); map(S,R,study)) assemblySpace = method() assemblySpace(List):=(Legs) -> ( eqs = sum apply(Legs, leg -> ideal homogenize(sub(eq, vars R|transpose sub(leg_0||leg_1,R)),t_R)); g := gens gb (eqs + SE3, Algorithm => LinearAlgebra); ideal(divideByVariable(g,t_R))_0) randomLegs = method() randomLegs(Ring,ZZ,Boolean) := (R,N,different) -> ( while( Legs:= apply(N, i -> apply(2, j -> random(R^3,R^1))); print(different and #(flatten Legs) < 2*N); different and #(unique flatten Legs) < 2*N) do (); Legs) -- construction of the genus 7 mechanism -- PART1 return the surface x -- the legs defined over ZZ and the parameter for the veronese familyOfSextics = method() familyOfSextics(List,Matrix):=(Legs, h) -> ( -- setup ring an SE3 R = ZZ[a_(0,0)..a_(2,2), b_0..b_2,t]; P = ZZ[p_0..p_2, q_0..q_2]; RP = R**P; setupSE3(R); setupEquation(RP); -- setup the parametrization PP^2 --> Veronese phi = eulerAngles(ZZ[s_0..s_3], ZZ[a_(0,0)..a_(2,2),t]); angles = phi(matrix{{a_(0,0), a_(1,0), a_(2,0), a_(0,1), a_(1,1), a_(2,1), a_(0,2),a_(1,2), a_(2,2),t}}); S = ZZ[s_1..s_3,L, b_0..b_2, Degrees => {4:1,3:3}]; parametrization = (L*sub(angles_{0..8},transpose (syz h * matrix{{s_1},{s_2},{s_3}}))) | matrix{{b_0,b_1,b_2}} | L*matrix{{sub(angles_{9}, transpose (syz h * matrix{{s_1},{s_2},{s_3}}))}}; -- compute the syzygy matrix N Q = apply(3, i -> map(S^1,S^{-6},sub(eq, parametrization | sub(Legs_i_0|Legs_i_1, S)))); D = gens (ideal(Q_0 - Q_1):L) | gens (ideal(Q_0-Q_2):L); M = transpose diff(transpose matrix{{L,b_0, b_1, b_2}}, D); N = syz M; -- compute the 2x2 minors and the corresponding ideal I T = ZZ[s_1,s_2,s_3,x,y,z]; LB = (s_1*x + s_2*y + s_3*z)*sub(N_{0},T) + sub(N_{1},T); f = sub(Q_2, matrix{{s_1,s_2,s_3}}|transpose LB); parametrization = sub(parametrization, matrix{{s_1,s_2,s_3}}|transpose LB); TP = ZZ[s_1,s_2,s_3, p_0..p_2,q_0..q_2,x,y,z]; generaleq = sub(eq, sub(parametrization, TP) | sub(vars P, TP)); monomsdeg6 = gens (ideal((vars TP)_{0..2}))^6; I = minors(2,contract(monomsdeg6,sub(f,TP)) || contract(monomsdeg6,generaleq)); -- return the sextic f, the line and I (f, N_(0,1),I)) triplePointsPreparation = method() -- we divide the preparation into 3 STEPS -- STEP 1 : Computation of the coefficient Matrix D -- STEP 2 : Computation of the Ideal Itriple -- STEP 3 : Computation of the distinguished z-Value triplePointsPreparation(Ideal,ZZ):=(I,prime) -> ( R = (ZZ/prime)[p_0..p_2, q_0..q_2, x,y,z, MonomialSize => 8]; -- we can compute the Groebner basis of I only over finite fields J = ideal mingens ideal mingens sub(I,R); -- simplification of I: -- STEP 1 : Dividing out the known legs (this takes some time!) Irr = ideal mingens ((ideal contract(matrix{{x,y,z}},gens J) + J)); J1 = ideal mingens (J:Irr); -- STEP 2 : Throw away a ueseless component Irr = ideal mingens ideal((gens J1)_{1}|contract(matrix{{p_1,p_2,q_0,q_1,q_2}}, (gens J1)_{1})); J2 = ideal mingens (J1:Irr); -- STEP 3 : Reduce modulo the linear equations J3 = ideal mingens ideal( gens J2 % (gens J2)_{0,1}); -- Computation of the matrix D of coefficients of the gens J and -- all multiples with the remaining variables p_2, q_0, q_1, q_2 var = matrix{{p_2, q_0, q_1, q_2, 1}}; monoms = symmetricPower(3,var); D = sub(contract(transpose monoms, flatten (transpose var * gens J3)), {p_2 => 0, q_0 => 0, q_1 => 0, q_2 => 0}); D) triplePointsPreparation(Matrix):=(D) -> ( -- iterative computation of another presentation D0 of D D0 = D; apply(30, i -> ( Syz = syz D0^{0}; D1 = D0^{1..(rank target D0)-1} * Syz; D0 = D1)); Syz = syz D0^{0,1}; -- computing the maximal minors of D0 Itriple = ideal 0_(ring D); apply(3, i -> Itriple = ideal mingens(Itriple + ideal mingens ideal(D0^{2+i} * Syz))); Itriple = ideal mingens ideal mingens ideal mingens Itriple; Itriple) triplePointsPreparation(Ideal,ZZ,ZZ):=(Itriple,prime,deg) -> ( select(toList(0..prime-1), c -> degree ideal mingens sub(Itriple, z => c) == deg)) triplePoints = method() triplePoints(Ideal):=(I) -> ( primes = select(toList(102..120), i -> isPrime i); -- pass to all primes in the list. For each p, scan through z values and save z if deg Itriple is 14 Mods = {}; zDists = {}; time apply(primes, prime -> time ( -- compute the distinguished z-value D = triplePointsPreparation(I,prime); Itriple = triplePointsPreparation(D); zDist = triplePointsPreparation(Itriple, 101, 14); if #zDist == 1 then ( Mods = append(Mods, prime); zDists = append(zDists,zDist_0);))); -- lift the distinguished z value zDist = ratConvert(chineseRemainder(zDists,Mods)); zDist = zDist_0 / zDist_1; primes = select(toList(100..500), i -> isPrime i); -- pass again to all primes and extract the coefficients of the quadric and the septic Mods = {}; coeffs = {}; time apply(primes, prime -> ( D = triplePointsPreparation(I,prime); Itriple = triplePointsPreparation(sub(D, z=> zDist)); QS = mingens Itriple; -- extract the coefficients if rank source QS == 2 then ( coeff = entries( sub(sub(contract(symmetricPower(2, matrix{{x,y,1}}), QS_(0,0)), {x => 0, y => 0}),ZZ)| sub(sub(contract(symmetricPower(7, matrix{{x,y,1}}), QS_(0,1)), {x => 0, y => 0}),ZZ)); coeffs = append(coeffs, coeff); Mods = append(Mods, prime)))); -- lift the coefficients and return ideal of triplepoints coeffvector = matrix{apply(#coeffs_0_0,i -> ( C = ratConvert(chineseRemainder(apply(coeffs, c -> c_0_i),Mods)); C_0/C_1))}; R = QQ[x,y,z]; monoms = symmetricPower(2, matrix{{x,y,1}})|symmetricPower(7, matrix{{x,y,1}}); Conic = sum apply(toList(0..5), i -> coeffvector_(0,i)* monoms_(0,i)); Septic = sum apply(toList(6..41), i -> coeffvector_(0,i)*monoms_(0,i)); Itriple = ideal(z - zDist, Conic, Septic); Itriple) remainingLegs = method() remainingLegs(List, Ideal, Ideal, Matrix) := (Legs, I, Itriple, triplePt) -> ( R = QQ[p_0..p_2, q_0..q_2,x,y,z]; -- simplification of I over QQ: -- STEP1: Plug in zDist zDist = z - (gens sub(Itriple,R))_(0,0); J = ideal mingens ideal mingens sub(sub(I,R), z => zDist); -- STEP2: Divide out the given legs: J0 = J + ideal mingens ideal(contract(matrix{{x,y}}, gens J)); apply(3, i -> J = (J : ideal((vars R)_{0..5} - (Legs_i_0|Legs_i_1)))); -- STEP3: Throw in the Conic and the Septic Conic = (gens Itriple)_(0,1); Septic = (gens Itriple)_(0,2); J2c = J + ideal(sub(Conic,R)); Jconics = J2c : ideal(sub(Septic,R)); Jrem=J2c:Jconics; -- Compute the remaining legs: Rreal = RR[p_0..p_2, q_0..q_2,x,y,z]; conicsr = ideal mingens ideal sub(sub((gens Jconics)_{0,1,2,3,5,8} ,Rreal), (vars Rreal)_{0..5}|sub(triplePt, Rreal)); sols = solveQuadratic((gens conicsr)_(0,5), q_2); LEG4=sub((matrix{{q_2,q_1,q_0,p_2,p_1,p_0}} - mingens (ideal (gens conicsr)_{0..4}+ideal (q_2-sols_0)))_{5,4,3,2,1,0},RR); LEG5=sub((matrix{{q_2,q_1,q_0,p_2,p_1,p_0}} - mingens (ideal (gens conicsr)_{0..4}+ideal (q_2-sols_1)))_{5,4,3,2,1,0},RR); Jrem = sub(Jrem, Rreal); sub(Jrem, (vars Rreal)_{0..5}|sub(triplePt, Rreal)); Qeq = mingens ideal sub((mingens Jrem)_{0,1,3,4,5,6,7},(vars Rreal)_{0..5}|sub(triplePt, Rreal)); sols2 = solveQuadratic( ( Qeq)_(0,6), q_2); LEG6=sub((matrix{{q_2,q_1,q_0,p_2,p_1,p_0}}- (matrix{{q_2 - sols2_0}}| (sub(Qeq, q_2 => sols2_0))_{1..5}))_{5,4,3,2,1,0},RR); ) end; restart load"StewartGoughPlatformsFile.m2"; h = matrix {{1,-1,0,0}}; Legs = {(matrix {{ 0, 0, 0}}, matrix {{ 3, 0, 0}}), (matrix {{-2,-2,-1}}, matrix {{-3, 3,-3}}), (matrix {{ 2, 1,-2}}, matrix {{ 3, 1, 1}})}; out = familyOfSextics(Legs,h); I = out_2; time triplePoints(out_2); -- pass the result to maple and return with numerical solution for a -- triplePoint: triplePt = matrix{{2.2072318327235894537139604979311561842021706500872, -7.4995700484259023116832401447984825811702987469853, 90/47}}; remainingLegs(Legs, I, Itriple, triplePt) time D = triplePointsPreparation(out_2,101); time Itriple = triplePointsPreparation(D); time triplePointsPreparation(Itriple, 101, 14) transpose mingens ideal mingens ideal mingens gb sub(Itriple, z => 32) LEG4 =