## ---------------------------------------------------------------------- ## Title: ChowForm_Code ## Created: Mon Jun 04 17:29:42 2018 ## Author: Clement Laroche ## Contributions of: ## C. Konaxis ## I.Z. Emiris ## ## Description: Implementation of an implicitization ## algorithm based on Chow Forms ## ## Version last modified 03/2019 ## ---------------------------------------------------------------------- ## ## ---------------------------------------------------------------------- with(LinearAlgebra): with(Statistics): read "multires.mpl": __ChowFormImplicitize_Help := NULL: addtohelp := proc() global __ChowFormImplicitize_Help; __ChowFormImplicitize_Help := __ChowFormImplicitize_Help, args: end proc: # Size of the box for point randomization # - Low value (1 .. 3): integer coefficients of the computations and outputs are much lower # - High value (4 ..): extremely low chances to obtain affine dependencies (they are checked anyway so it only may spare a few randomization rerun) ChowFormImplicitize_RandomRange := 2: #XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX# #XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX# # ALGORITHM # #XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX# #XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX# #---------------------------------------------------------------------- ## HELP ChowFormImplicitize_GetRandomPointG ## ## ChowFormImplicitize_GetRandomPointG ## - Generate random point(s) that can be used as an apex of a conical surface ## ## CALLING SEQUENCE: ## ChowFormImplicitize_GetRandomPointG(d, n, param, quotient) ## ## PARAMETERS: ## d - dimension of the variety ## n - dimension of the ambient space ## param - numerators of the parametric equations ## quotient - common denominator of the parametric equations ## ## DESCRIPTION: ## - The output should be out of the curve but relatively close to it. ## - Use the global variable ChowFormImplicitize_RandomRange to specify both ## the base point on the curve around which the output is picked and ## the size of the box in which the output is picked. ## ## EXAMPLES: ##> G := ChowFormImplicitize_GetRandomPointG(1, 3, t -> (t, t^2, t^3), t -> 1) ## ## [[ 1]] ## G := [[ 2]] ## [[10]] ## ##> ComputeConicalSurface((t -> (t, t^2, t^3)), (t -> 1)) ## ## -92*x[1]^3+12*x[3]*x[1]^2+60*x[1]^2+48*x[1]^2*x[2]+80*x[1]*x[2] ## -24*x[1]*x[2]^2-28*x[3]*x[1]+x[3]^2*x[1]-8*x[1]*x[3]*x[2]-60*x[2] ## -x[3]*x[2]^2-x[3]^2+12*x[3]-20*x[2]^2+9*x[2]^3+12*x[3]*x[2] ## ## SEE ALSO: ## ComputeConicalSurface ## ChowFormImplicitize_GetRandomPointG := proc (d, n, param, quotient) local G, matG, Pbase, PbaseT; global ChowFormImplicitize_RandomRange; PbaseT := seq((rand(-ChowFormImplicitize_RandomRange .. ChowFormImplicitize_RandomRange))(), k = 1 .. d): while quotient(PbaseT) = 0 do PbaseT := seq((rand(-ChowFormImplicitize_RandomRange .. ChowFormImplicitize_RandomRange))(), k = 1 .. d): end do: Pbase := Vector[column]([seq([param(PbaseT)][k]/quotient(PbaseT), k = 1 .. n)]): G := [seq(Pbase + Vector[column]([seq((rand(1 .. ChowFormImplicitize_RandomRange))()*(2*(rand(0 .. 1))()-1), k = 1 .. n)]), i = 1 .. n-d-1)]: matG := Matrix([seq(G[k]-G[1], k = 2 .. n-d-1)]): while LinearAlgebra:-Rank(matG) < n-d-2 do G := [seq(Pbase + Vector[column]([seq((rand(1 .. ChowFormImplicitize_RandomRange))()*(2*(rand(0 .. 1))()-1), k = 1 .. n)]), i = 1 .. n-d-1)]: matG := Matrix([seq(G[k]-G[1], k = 2 .. n-d-1)]): end do: return G: end proc: #---------------------------------------------------------------------- ## HELP ChowFormImplicitize_GetRandomPointP ## ## ChowFormImplicitize_GetRandomPointP ## - Generate random points that can be used to define random hyperplanes ## ## CALLING SEQUENCE: ## ChowFormImplicitize_GetRandomPointP(n, G, amount) ## ## PARAMETERS: ## n - dimension of the ambient space ## G - point(s) around which the outputs are placed ## amount - number of points required ## ## DESCRIPTION: ## - The output forms an affinely independant set with G. ## - Note that the input point(s) G must be affinely independant. ## - Use the global variable ChowFormImplicitize_RandomRange to specify ## the size of the box in which the output is picked. ## ## EXAMPLES: ##> P := ChowFormImplicitize_GetRandomPointP(3, Vector[column]([1, 0, 0]), 2) ## ## [[ 2] [ 0]] ## P := [[-1], [ 1]] ## [[-1] [-1]] ## ## SEE ALSO: ## ComputeConicalSurface ## ChowFormImplicitize_GetRandomPointP := proc (n, G, amount) local P, matP, listG: global ChowFormImplicitize_RandomRange; if is(G, Vector) then listG := [G]: else: listG := G: end if: Vector[column]([seq((rand(1 .. ChowFormImplicitize_RandomRange))()*(2*(rand(0 .. 1))()-1), k = 1 .. n)]): P := [seq(listG[1] + Vector[column]([seq((rand(1 .. ChowFormImplicitize_RandomRange))()*(2*(rand(0 .. 1))()-1), k = 1 .. n)]), i = 1 .. amount)]: matP := Matrix([seq(listG[k] - listG[1], k = 2 .. nops(listG)), seq(P[k] - listG[1], k = 1 .. amount)]): while LinearAlgebra:-Rank(matP) < min(n, amount + nops(listG) - 1) do P := [seq(listG[1] + Vector[column]([seq((rand(1 .. ChowFormImplicitize_RandomRange))()*(2*(rand(0 .. 1))()-1), k = 1 .. n)]), i = 1 .. amount)]: matP := Matrix([seq(listG[k] - listG[1], k = 2 .. nops(listG)), seq(P[k] - listG[1], k = 1 .. amount)]): end do: return P: end proc: #---------------------------------------------------------------------- addtohelp(ComputeConicalSurface): #---------------------------------------------------------------------- ## HELP ComputeConicalSurface ## ## ComputeConicalSurface ## - Given the parametric equations of a curve C and point(s) G, ## compute the conical (hyper)surface S of directrix C and apex G ## ## CALLING SEQUENCES: ## ComputeConicalSurface(param) # Polynomial parameterization ## ComputeConicalSurface(param, quotient) # Rational parameterization ## ComputeConicalSurface(param, quotient, apex) # Rational parameterization with specified apex ## ## PARAMETERS: ## param - numerators of the parametric equations R -> R^n ## quotient - common denominator of the parametric equations ## apex - point(s) used as the apex of the conical surface, ## randomized if it is not provided ## ## DESCRIPTION: ## - If the apex G is given, it should be a single point if n = 3 or a list ## [G[1], ..., G[n-2]] if n > 3, where n = nops([param(t)]). ## - If the apex G is omitted and the curve is (hyper)planar, ## then the equation of a hyperplane containing it is returned instead. ## - If the apex G is omitted and the curve is not (hyper)planar, ## G is randomized using ChowFormImplicitize_GetRandomPointG and appended to the output. ## In that case, the output is of the form (S, [G[1], ..., G[n-2]]). ## - This algorithm is generalized to allow a rational parameterized variety ## of any dimension 1 <= d <= n-2; see ComputeConicalSurfaceExtended. ## ## EXAMPLES: ##> ComputeConicalSurface(t -> (t, t^2, t^3), t -> 1) ## ## 4*x[1]^3-4*x[1]^2*x[3]-4*x[1]^2+4*x[1]*x[3]+x[1]*x[3]^2 [[1]] ## +4*x[2]+4*x[3]*x[2]-x[3]*x[2]^2-4*x[2]^2-4*x[3]+x[2]^3-x[3]^2, [[2]] ## [[2]] ## ##> ComputeConicalSurface(t -> (1-t^2, 2t, 1+t^2), t -> 1+t^2) ## ## 24-24*x[3] ## ##> ComputeConicalSurface(t -> (1-t^2, 2t, 1+t^2), t -> 1+t^2, Vector[column]([0, 0, 0])) ## ## 4*x[1]^2+4*x[2]^2-4*x[3]^2 ## ##> factor(ComputeConicalSurface(t -> (t, t^2, t^3, t^4), t -> 1, [Vector([0, 1, 0, 0]), Vector([0, 0, 0, 1])])) ## ## x[4]^2*x[3]*x[1]-2*x[4]*x[3]*x[1]-2*x[1]*x[3]*x[2]+2*x[4]*x[1]*x[3]*x[2]-x[1]^4 ## -x[3]^4+2*x[3]*x[1]^3-2*x[1]*x[3]^3+x[3]^2*x[1]^2+x[3]*x[1]+x[3]*x[1]*x[2]^2 ## ##> ChowFormImplicitize_RandomRange := 10: ##> ComputeConicalSurface(t -> (t, t^2, t^3), t -> 1) ## ## 152*x[1]^3+43746*x[1]^2-1268*x[1]^2*x[3]-7854*x[1]^2*x[2]+2310*x[1]*x[2]^2 [[ 5]] ## -8*x[1]*x[3]^2+154*x[1]*x[3]*x[2]+10626*x[1]*x[2]+8958*x[1]*x[3]-43746*x[2] [[ 17]] ## +40*x[3]^2-1104*x[2]^2-10778*x[3]-1042*x[3]*x[2]+8*x[3]*x[2]^2-194*x[2]^3, [[-69]] ## ## SEE ALSO: ## ComputeDefiningSurfaces ## ComputeConicalSurfaceExtended ## ComputeConicalSurface := proc (param::procedure, quotient::procedure := (t -> 1), apex := NULL) local d, n, G, P, H1, H2, L, result, plane, dgr, rnk, teval, i, X, apexgiven: description "Given n parametric equations (n >= 3) of 1 variable, returns a polynomial in [x_1,..., x_n] defining a conical surface or a hyperplane that contains the parametric subspace. Also returns the apex of the conical surface if applicable": if nops([param(t)]) < 3 then error "Invalid parameter %1: expected a polynomial R -> R^n with n >= 3", param: end if: if nops([quotient(t)]) <> 1 then error "Invalid parameter %1: expected a polynomial R -> R", quotient: end if: d := 1: n := nops([param(t)]): X := Vector[column]([seq(x[i], i = 1 .. n)]): dgr := degree(quotient(t), t): for i to n do dgr := max(degree(param(t)[i], t), dgr): end do: if dgr = 0 then error "Unexpected curve of degree 0": end if: apexgiven := apex <> NULL: if apexgiven then if n = 3 and is(apex, Vector) then G := [apex]: else G := apex: end if: else teval := [seq(2*i-n+(rand(0 .. 1))(), i = 1 .. n+1)]: rnk := LinearAlgebra:-Rank(Matrix([seq(Vector[column]([param(teval[i]), quotient(teval[i])]), i = 1 .. n+1)])): if rnk <= n then # Most likely, the curve is contained in a rnk-dimensional linear subspace P := ColumnSpace(Matrix([seq(Vector[column]([param(teval[i]), quotient(teval[i])]), i = 1 .. n+1)])): plane := varx -> Determinant(Matrix([op(P), op(ChowFormImplicitize_GetRandomPointP(n + 1, P, n - rnk)), varx])): if plane(Vector[column]([param(t), quotient(t)])) = 0 then return plane(ArrayTools:-Concatenate(1, X, 1)): end if: end if: G := ChowFormImplicitize_GetRandomPointG(d, n, param, quotient): end if: P := ChowFormImplicitize_GetRandomPointP(n, G, 2): H1 := xi -> (varx -> Determinant(Matrix(n, n, [P[1]-G[1], seq(G[i]-G[1], i = 2 .. n-d-1), xi-quotient(t)*G[1], varx-G[1]]))): H2 := xi -> (varx -> Determinant(Matrix(n, n, [P[2]-G[1], seq(G[i]-G[1], i = 2 .. n-d-1), xi-quotient(t)*G[1], varx-G[1]]))): L := Determinant(Matrix(n, n, [P[1]-G[1], P[2]-G[1], seq(G[i]-G[1], i = 2 .. n-d-1), X-G[1]])): result := evala(Resultant((H1(Vector[column]([param(t)])))(X), (H2(Vector[column]([param(t)])))(X), t)): if apexgiven then return simplify(result/L^dgr): end if: return simplify(result/L^dgr), G: end proc: #---------------------------------------------------------------------- addtohelp(ComputeDefiningSurfaces): #---------------------------------------------------------------------- ## HELP ComputeDefiningSurfaces ## ## ComputeDefiningSurfaces ## - Compute a set of planes or conical surfaces that define the curve ## ## CALLING SEQUENCES: ## ComputeDefiningSurfaces(param) # Polynomial parameterization ## ComputeDefiningSurfaces(param, quotient) # Rational parameterization ## ## PARAMETERS: ## param - numerators of the parametric equations ## quotient - common denominator of the parametric equations ## ## DESCRIPTION: ## - As for the other procedures of the package, the input must be given as ## a polynomial param: R -> R^n such that nops([param(t)]) = n. ## - If the curve is included in a hyperplane, this hyperplane's equation is returned ## along with other conical surfaces the intersection of which is the curve. ## - Else, only equations of conical surfaces are returned. ## - Therefore either 2 or 3 equations are returned for a space curve (n = 3). ## ## EXAMPLES: ##> ComputeDefiningSurfaces(t -> (1-t^2, 2t, 1+t^2), t -> 1+t^2) ## ## {1-x[3], ## 4*x[1]^2-8*x[3]*x[1]+8*x[1]+36*x[3]^2+4*x[2]^2 ## -24*x[2]+40+24*x[3]*x[2]-80*x[3]} ## ##> eq := ComputeDefiningSurfaces(t -> (1-t^2, 2t, 1+t^3), t -> 1+t^2) ## ## {-8*x[2]^3+16*x[2]*x[1]*x[3]-16*x[2]*x[1]^2-16*x[3]*x[2] ## +16*x[2]+8-24*x[1]-8*x[1]^3+24*x[1]^2, ## -32*x[1]^3+64*x[3]*x[1]^2+48*x[2]*x[1]^2-24*x[2]^2*x[1] ## -48*x[2]*x[1]*x[3]-32*x[3]^2*x[1]+32*x[3]*x[1]-32*x[3]^2 ## +48*x[2]+8*x[2]^3-48*x[3]*x[2]+32+24*x[2]^2-32*x[3], ## -36*x[1]^3-24*x[3]*x[1]^2+18*x[2]*x[1]^2-18*x[1]^2-72*x[3]*x[1] ## -48*x[2]*x[1]*x[3]-4*x[3]^2*x[1]+108*x[1]-72*x[2]^2*x[1]+36*x[2]*x[1] ## -48*x[3]*x[2]-48*x[2]^2-4*x[3]^2+8*x[2]^3+98+42*x[2]-48*x[3]} ## ##> seq(simplify(eval(eq[i], [x[1] = (1-t^2)/(1+t^2), x[2] = 2t/(1+t^2), x[3] = (1+t^3)/(1+t^2)])), i = 1 .. 3) ## ## 0, 0, 0 ## ## SEE ALSO: ## ComputeConicalSurface ## ComputeDefiningSurfaces := proc (param::procedure, quotient::procedure := (t -> 1)) local s1, s2, s3, g1, g2, g3, result, i, j, k, n: description "Compute a set of planes or conical surfaces that define the curve.": n := nops([param(t)]): result := ComputeConicalSurface(param, quotient): if nops([result]) = 1 then # Defining equations consist of: # * One linear equation of an hyperplane (result) # * n conical hypersurfaces with the apex out of that hyperplane # Notes: # * An improvement would be to determine if the curve is contained # in an affine subsapce of any dimension # * If the curve is contained in a plane (of dim 2), # we can compute a set of defining equations of size n-1 instead if n = 3 then g2 := ChowFormImplicitize_GetRandomPointG(1, n, param, quotient): while eval(result, {seq(x[i] = g2[i], i = 1 .. n)}) = 0 do g2 := ChowFormImplicitize_GetRandomPointG(1, n, param, quotient): end do: s2 := ComputeConicalSurface(param, quotient, g2): return `union`({result}, {s2}): else result := {result}: for i from 1 to n do g2 := ChowFormImplicitize_GetRandomPointG(1, n, param, quotient): j := 1: while j <= n-2 do if eval(result, {seq(x[k] = g2[j][k], k = 1 .. n)}) = 0 then g2 := ChowFormImplicitize_GetRandomPointG(1, n, param, quotient): j := 1: else: j := j + 1: end if: end do: s2 := ComputeConicalSurface(param, quotient, g2): result := `union`(result, {s2}): end do: return result: end if: end if: # Defining equations consist of m conical hypersurfaces # with m = n if n = 3 # or m = n+1 if n>3 s1, g1 := result: s2, g2 := ComputeConicalSurface(param, quotient): while g1 = g2 do s2, g2 := ComputeConicalSurface(param, quotient) end do: if n = 3 then g3 := ChowFormImplicitize_GetRandomPointG(1, n, param, quotient): while LinearAlgebra:-Rank(Matrix([g1[1]-g3[1], g2[1]-g3[1]])) < 2 do g3 := ChowFormImplicitize_GetRandomPointG(1, n, param, quotient): end do: s3 := ComputeConicalSurface(param, quotient, g3): else s3, g3 := ComputeConicalSurface(param, quotient): end if: result := {s1, s2, s3}: if n = 3 then return result: end if: for i from 3 to n do result := `union`(result, {ComputeConicalSurface(param, quotient)[1]}): end do: return result: end proc: #---------------------------------------------------------------------- addtohelp(ComputeConicalSurfaceExtended): #---------------------------------------------------------------------- ## HELP ComputeConicalSurfaceExtended ## ## ComputeConicalSurfaceExtended ## - Given the parametric equations of a variety V and point(s) G, compute a ## variety containing the conical hypersurface S of directrix V and apex G ## ## CALLING SEQUENCES: ## ComputeConicalSurfaceExtended(d, n, param) # Polynomial parameterization ## ComputeConicalSurfaceExtended(d, n, param, quotient) # Rational parameterization ## ComputeConicalSurfaceExtended(d, n, param, quotient, apex) # Rational parameterization with specified apex ## ComputeConicalSurfaceExtended(d, n, param, [quotient, apex], method = "macaulay") # Return a macaulay matrix instead of sparse matrix ## ## PARAMETERS: ## d - dimension of the variety (must satisfy 1 <= d <= n-2) ## n - dimension of the ambient space ## param - numerators of the parametric equations ## quotient - common denominator of the parametric equations ## apex - set of n-d-1 point(s) used as the apex of the conical hypersurface, ## randomized if it is not provided ## method - type of resultant matrix used: "sparse" (default) or "macaulay" ## ## DESCRIPTION: ## - Unlike with "ComputeConicalSurface" the output is a matrix, the determinant ## of which is the equation of an hypersurface containing the curve. ## - The output's determinant is an equation of the form Det = S^q * E[1] * E[2]^p ## where S is the equation of the conical hypersurface, ## E[1] is an extraneous factor coming from the sparse resultant matrix ## (E[1] can be expressed as one of the output's minors), ## and E[2] is a degree d equation that is raised in Det's formula ## to some power p such that 1 <= p <= degree(param)^d. ## - The polynomial S^q * E[2] is of degree (d+1)*degree(param)^d. ## - The number and complexity of V's base points increase the power q ## and decrease the power p. ## - The non-properness of the input parameterization increases the power q. ## - The "method" argument can be either "sparse" (default) or "macaulay", ## which returns a Macaulay matrix instead of a sparse matrix. ## The Macaulay matrix not only is bigger but also has a null determinant oftentimes. ## ## EXAMPLES: ##> conical_matrix := ComputeConicalSurfaceExtended(2, 4, (s,t) -> (t+s, t^2, t^3, t^4)) ## ## [ 12 x 12 Matrix ] ## conical_matrix := [ Data Type: anything ] ## [ Storage: rectangular ] ## [ Order: Fortran_order ] ## ##> full_equation := factor(Determinant(conical_matrix)) ## ## full_equation := -11664 ## *(9*x[4]^3-3*x[2]*x[4]^3-27*x[2]*x[4]^2-14*x[3]^2*x[2]*x[4]-30*x[3]^2*x[2] ## +12*x[3]^2*x[2]^2+9*x[2]^2*x[4]^2+27*x[2]^2*x[4]+3*x[3]^2*x[4]^2 ## +24*x[3]^2*x[4]+9*x[3]^2-9*x[2]^3*x[4]-9*x[2]^3-4*x[3]^4+3*x[2]^4) ## *(-12-x[4]+16*x[3]+5*x[2]+6*x[1]+4*x[4]^2-x[3]*x[4]-4*x[3]^2-5*x[2]*x[4] ## -3*x[3]*x[2]+x[2]^2-4*x[1]*x[4]-2*x[1]*x[3]+2*x[1]*x[2])^4 ## ##> eval(op(full_equation)[2], [x[1] = t+s, x[2] = t^2, x[3] = t^3, x[4] = t^4]) ## ## 0 ## ## SEE ALSO: ## ComputeConicalSurface ## spresultant (from "multires.mpl") ## mresultant (from "multires.mpl") ## ComputeConicalSurfaceExtended := proc (d::integer, n::integer, param::procedure, quotient::procedure := (t -> 1), apex := NULL, { method::string := "sparse" }) local i, j, vart, G, P, H, Pbase, X, L, result, dgr, apexgiven: description "Given n parametric equations of d variables, returns a polynomial whose kernel contains the parametric subspace": if method <> "sparse" and method <> "macaulay" then error "Invalid method \"%1\": expected either \"sparse\" or \"macaulay\"", method: end if: #if d = 1 then # return ComputeConicalSurface(param, quotient, apex): #end if: vart := seq(t[i], i = 1 .. d): X := Vector[column]([seq(x[i], i = 1 .. n)]): dgr := degree(param(vart)[1], t): for i from 2 to n do dgr := max(degree(param(vart)[i], [vart]), dgr): end do: apexgiven := apex <> NULL: if apexgiven then if n = d+2 and is(apex, Vector) then G := [apex]: else G := apex: end if: else G := ChowFormImplicitize_GetRandomPointG(d, n, param, quotient): end if: for i to d do P[i] := ChowFormImplicitize_GetRandomPointP(n, G, d+1): end do: if quotient(vart) = 1 then H := l -> (xi -> (x -> Determinant(Matrix(n, n, [seq(G[k]-P[1][l], k = 1 .. n-d-1), seq(P[k][l]-P[1][l], k = 2 .. d), xi-P[1][l], x-P[1][l]])))): else H := l -> (xi -> (x -> Determinant(Matrix(n, n, [seq(G[k]-P[1][l], k = 1 .. n-d-1), seq(P[k][l]-P[1][l], k = 2 .. d), xi-quotient(vart)*P[1][l], x-P[1][l]])))): end if: #if d = 1 then # L := Determinant(Matrix(n, n, [P[1][1]-G[1], P[1][2]-G[1], seq(G[i]-G[1], i = 2 .. n-d-1), X-G[1]])): # result := evala(Resultant(((H(1))(Vector[column]([param(t)])))(X), ((H(2))(Vector[column]([param(t)])))(X), t)): # return simplify(result/L^dgr): #end if: if method = "sparse" then return Matrix(eval(spresultant([seq(((H(k))(Vector[column]([param(vart)])))(X), k = 1 .. d+1)], [vart])[1])): end if: return Matrix(eval(mresultant([seq(((H(k))(Vector[column]([param(vart)])))(X), k = 1 .. d+1)], [vart]))): end proc: #----------------------------------------------------------------------# [__ChowFormImplicitize_Help];