(* Last looked at: 12/18/2000 *) (* Programs to polarize parametric rational curves, and surfaces, rectangular and triangular (11/25/2002) *) (* gbicontnet[{cpoly__},p_,q_,r1_,s1_,r2_,s2_] Computes a rectangular control net of bidegree (p, q) from some polar forms, *) (* w.r.t. affine frames ((r1, s1), (r2, s2) *) (* The input is a list of polynomials, one for each coordinate, including the denominator (x, y, z, w); Each polynomial is a list of monomials; Each monomial c U^i V^k is a triple {c, i, k}. Example: Moebius strip, rectangular, degree (6, 1): poly = {{{2, 0, 0}, {-10, 2, 0}, {-10, 4, 0}, {2, 6, 0}, {2, 1, 1}, {-12, 3, 1}, {2, 5, 1}}, {{8, 1, 0}, {-8, 5, 0}, {8, 2, 1}, {-8, 4, 1}}, {{1, 0, 1}, {1, 2, 1}, {-1, 4, 1}, {-1, 6, 1}}, {{1, 0, 0}, {3, 2, 0}, {3, 4, 0}, {1, 6, 0}}}; *) (* fcontpoly[{cpoly__},m_,r1_,s1_] computes the control polygon of a rational curve of degree m *) (* specified by polynomials in cpoly, w.r.t. frame [r1, s1] *) (* fast method using a recurrence relation *) (* The input is a list of polynomials, one for each coordinate, including the denominator (x, y, z, w); Each polynomial is a list of monomials; Each monomial c t^i is a triple {c, i}. Example: seven-leaved rose (degree 8): poly = {{{7, 1}, {-35, 3}, {21, 5}, {-1, 7}}, {{7, 2}, {-35, 4}, {21, 6}, {-1, 8}}, {{1, 0}, {4, 2}, {6, 4}, {4, 6}, {1, 8}}} *) (* fgcontnet[{cpoly__},m_,reftrig_] Computes a triangular control net of degree m *) (* w.r.t. reference triangle *) (* (r, s, t) = ((r1, r2, r3), (s1, s2, s3), (t1, t2, t3)) *) (* fast method using a recurrence relation *) (* The input is a list of polynomials, one for each coordinate, including the denominator (x, y, z, w); Each polynomial is a list of monomials; Each monomial c U^i V^k is a triple {c, i, k}. Example, sphere, total degree 2: sphpoly = {{{2, 1, 0}, {-2, 2, 0}, {-2, 1, 1}}, {{2, 0, 1}, {-2, 0, 2}, {-2, 1, 1}}, {{2, 1, 0}, {2, 0, 1}, {-2, 1, 1}, {-1, 0, 0}}, {{2, 2, 0}, {2, 0, 2}, {-2, 1, 0}, {-2, 0, 1}, {2, 1, 1}, {1, 0, 0}}}; *) (*************************************************************************) (* This program computes the polar form of a polynomial in *) (* two variables of total degree d with respect to m >= d *) (* Returns the relative complement A - B *) relcomp[alist_,blist_] := Block[ {A = alist, B = blist, l1, l2, i, j, result, found}, l1 = Length[A]; l2 = Length[B]; result = {}; If[A =!= {}, If[B =!= {}, Do[ found = 0; j = 1; While[found === 0 && j <= l2, If[A[[i]] === B[[j]], found = 1]; j = j + 1 ]; If[found === 0, result = Append[result, A[[i]]]], {i, 1, l1} ], result = A ] ]; result ]; (* Adds an element at the end of every list in a list of lists *) addright[{alist__},a_] := Block[ {S = {alist}, i, row, l, result}, l = Length[S]; result = {}; Do[ row = S[[i]]; row = Append[row, a]; result = Join[result, {row}], {i, 1, l} ]; result ]; (* Adds an element in the front of every list in a list of lists *) addleft[{alist__},a_] := Block[ {S = {alist}, i, row, l, result}, l = Length[S]; result = {}; Do[ row = S[[i]]; row = Prepend[row, a]; result = Join[result, {row}], {i, 1, l} ]; result ]; (* Computes the set of all subsets of size n of a set of size m *) (* Simple-minded method, first compute all subsets, and then *) (* weed out those not of size n *) subsetslow[aset_,n_] := Block[ {S = aset, m, i, temp, l1, l, result}, m = Length[S]; temp = {{}}; result = {}; If[m =!= 0 && n =!= 0, Do[ temp = Join[temp, addright[temp, S[[i]]]], {i, 1, m} ] ]; l1 = Length[temp]; Do[ l = Length[temp[[i]]]; If[l === n, result = Join[result, {temp[[i]]}]], {i, 1, l1} ]; result ]; (* Computes the set of all subsets of size n of a set of size m *) (* Uses a stack to keep track of the recursion tree *) (* Entries corresponding to left branches are coded as *) (* {0, R1, p1, S1, done1}, where R1 is a subset of the original set *) (* such that if card{R1} = n1, then R1 is obtained from S by *) (* dropping the first m - n1 elements. *) (* Entries corresponding to right branches are coded as *) (* {1, R2, p2, S2, done2}, where R2 is a subset of the original set *) (* such that if card{R2} = n2, then R2 is obtained from S by *) (* dropping the first m - n2 elements. *) subsets[aset_,n_] := Block[ {S = aset, WS, temp, a, b, c, d, m, stack, S1, S2, S3, p1, p2, R1, R2, R3, n1, n2, done1, done2, nn, i, result, tag1, tag2, top, tag}, stack = {}; m = Length[S]; If[n =!= 0 && n < m, WS = Drop[S, 1]; a = {0, WS, n - 1, {}, 0}; b = {1, WS, n, {}, 0}; c = {0, S, n, {}, 0}; stack = {a, b, c}; While[Length[stack] =!= 1, a = stack[[1]]; stack = Drop[stack, 1]; b = stack[[1]]; stack = Drop[stack, 1]; tag1 = a[[1]]; R1 = a[[2]]; p1 = a[[3]]; S1 = a[[4]]; n1 = Length[R1]; done1 = a[[5]]; tag2 = b[[1]]; R2 = b[[2]]; p2 = b[[3]]; S2 = b[[4]]; n2 = Length[R2]; done2 = b[[5]]; If[done1 === 0, (* In this case, S1 is not computed yet *) If[p1 <= n1 && p1 =!= 0, stack = Prepend[stack, b]; (* puts back b on stack *) (* and recursive calls on a *) R1 = Drop[R1, 1]; stack = Prepend[stack, a]; (* puts back a on stack *) b = {1, R1, p1, {}, 0}; stack = Prepend[stack, b]; a = {0, R1, p1 - 1, {}, 0}; stack = Prepend[stack, a], (* In this case, we can compute S1 *) If[p1 === 0, S1 = {{}}]; If[p1 === n1, S1 = {R1}]; If[p1 > n1, S1 = {}]; a = {tag1, R1, p1, S1, 1}; stack = Prepend[stack, a]; stack = Prepend[stack, b] (* puts back b on stack *) ], (* In this case, S1 HAS been computed *) If[done2 === 0, (* In this case, S2 is not computed yet *) stack = Prepend[stack, a]; (* puts back a on stack *) (* and recursive calls on b *) If[p2 <= n2 && p2 =!= 0, R2 = Drop[R2, 1]; stack = Prepend[stack, b]; (* puts back b on stack *) b = {1, R2, p2, {}, 0}; stack = Prepend[stack, b]; a = {0, R2, p2 - 1, {}, 0}; stack = Prepend[stack, a], (* In this case, we can compute S2 *) If[p2 === 0, S2 = {{}}]; If[p2 === n2, S = {R}]; If[p2 > n2, S2 = {}]; b = {tag2, R2, p2, S2, 1}; stack = Prepend[stack, b] ], (* In this case, both S1 AND S2 HAVE been computed *) d = S[[m - n1]]; If[tag1 === 1, S3 = S2; S2 = S1; S1 = S3]; S3 = addleft[S1, d]; top = stack[[1]]; stack = Drop[stack, 1]; tag = top[[1]]; R3 = top[[2]]; nn = top[[3]]; S3 = Join[S3, S2]; c = {tag, R3, nn, S3, 1}; stack = Prepend[stack, c] ] ] ]; (* end While *) result = stack[[1]][[4]] ]; If[m < n, result = {}]; If[m === n, result = {S}]; If[n === 0, result = {{}}]; result ]; (* Computes the polar form of a monomial of total degree d = h + k, *) (* w.r.t degree m, as a formal polynomial in U[i], V[i] *) polarmon[{amon__}, m_] := Block[ {mon = {amon}, Ilist, Jlist, III, JJJ, setm, rest, a, h, k, d, uu, vv, i, j, l1, l2, I1, J1, ii, jj, result, b, p, q}, a = mon[[1]]; h = mon[[2]]; k = mon[[3]]; d = h + k; setm = Table[i, {i, m}]; result = 0; Ilist = subsets[setm, h]; l1 = Length[Ilist]; Do[ III = Ilist[[I1]]; rest = relcomp[setm, III]; Jlist = subsets[rest, k]; l2 = Length[Jlist]; Do[ JJJ = Jlist[[J1]]; b = 1; Do[ ii = III[[i]]; b = b * u[ii], {i, 1, Length[III]} ]; Do[ jj = JJJ[[j]]; b = b * v[jj], {j, 1, Length[JJJ]} ]; result = result + b, {J1, 1, l2} ], {I1, 1, l1} ]; result = (a * result)/Multinomial[h, k, m - d]; result ]; (* Computes the polar form of a monomial of total degree d = h + k, *) (* w.r.t degree m, as a list {coeff, {mononial list}, coeff} *) (* where {monomial list} is a list of pairs of *) (* lists {i1, ... ih} and {j1, ... jk} indicating *) (* the indices of variables Ui and Vj of power 1 in a monomial *) polarlismon[{amon__}, m_] := Block[ {mon = {amon}, Ilist, Jlist, III, JJJ, setm, rest, a, h, k, d, i, j, l1, l2, I1, J1, ii, jj, result, Umon, Vmon, pmon, p, q}, a = mon[[1]]; h = mon[[2]]; k = mon[[3]]; d = h + k; setm = Table[i, {i, m}]; result = {}; Ilist = subsets[setm, h]; l1 = Length[Ilist]; Do[ III = Ilist[[I1]]; rest = relcomp[setm, III]; Jlist = subsets[rest, k]; l2 = Length[Jlist]; Do[ JJJ = Jlist[[J1]]; Umon = {}; Vmon = {}; Do[ ii = III[[i]]; Umon = Append[Umon, ii], {i, 1, Length[III]} ]; Do[ jj = JJJ[[j]]; Vmon = Append[Vmon, jj], {j, 1, Length[JJJ]} ]; pmon = Join[{Umon}, {Vmon}]; result = Join[result, {pmon}], {J1, 1, l2} ], {I1, 1, l1} ]; result = {a , result, Multinomial[h, k, m - d]}; result ]; (* Computes the polar form of a monomial of bidegree (h, k), *) (* w.r.t bidegree (p, q), as a formal polynomial in U[i], V[i] *) bipolarmon[{amon__}, p_, q_] := Block[ {mon = {amon}, Ilist, Jlist, III, JJJ, setp, setq, a, h, k, uu, vv, i, j, l1, l2, I1, J1, ii, jj, result, b}, a = mon[[1]]; h = mon[[2]]; k = mon[[3]]; setp = Table[i, {i, p}]; setq = Table[i, {i, q}]; result = 0; Ilist = subsets[setp, h]; Jlist = subsets[setq, k]; l1 = Length[Ilist]; l2 = Length[Jlist]; Do[ Do[ III = Ilist[[I1]]; JJJ = Jlist[[J1]]; b = 1; Do[ ii = III[[i]]; b = b * u[ii], {i, 1, Length[III]} ]; Do[ jj = JJJ[[j]]; b = b * v[jj], {j, 1, Length[JJJ]} ]; result = result + b, {J1, 1, l2} ], {I1, 1, l1} ]; result = a*result/(Binomial[p, h]*Binomial[q, k]); result ]; (* Computes the polar form of a monomial of bidegree (h, k), *) (* w.r.t bidegree (p, q), as a list {coeff, {mononial list}, coeff} *) (* where {monomial list} is a list of pairs of *) (* lists {i1, ... ih} and {j1, ... jk} indicating *) (* the indices of variables Ui and Vj of power 1 in a monomial *) bipolarlismon[{amon__}, p_, q_] := Block[ {mon = {amon}, Ilist, Jlist, III, JJJ, setp, setq, a, h, k, i, j, l1, l2, I1, J1, ii, jj, result, Umon, Vmon, pmon}, a = mon[[1]]; h = mon[[2]]; k = mon[[3]]; setp = Table[i, {i, p}]; setq = Table[i, {i, q}]; result = {}; Ilist = subsets[setp, h]; Jlist = subsets[setq, k]; l1 = Length[Ilist]; l2 = Length[Jlist]; Do[ Do[ III = Ilist[[I1]]; JJJ = Jlist[[J1]]; Umon = {}; Vmon = {}; b = 1; Do[ ii = III[[i]]; Umon = Append[Umon, ii], {i, 1, Length[III]} ]; Do[ jj = JJJ[[j]]; Vmon = Append[Vmon, jj], {j, 1, Length[JJJ]} ]; pmon = Join[{Umon}, {Vmon}]; result = Join[result, {pmon}], {J1, 1, l2} ], {I1, 1, l1} ]; result = {a, result, Binomial[p, h]*Binomial[q, k]}; result ]; (* Computes the polar form of a monomial of degree h, *) (* w.r.t degree m, as a formal polynomial in U[i] *) cpolarmon[{amon__}, m_] := Block[ {mon = {amon}, Ilist, III, setm, a, h, uu, i, l1, I1, ii, result, b}, a = mon[[1]]; h = mon[[2]]; setm = Table[i, {i, m}]; result = 0; Ilist = subsets[setm, h]; l1 = Length[Ilist]; Do[ III = Ilist[[I1]]; b = 1; Do[ ii = III[[i]]; b = b * u[ii], {i, 1, Length[III]} ]; result = result + b, {I1, 1, l1} ]; result = a*result/Binomial[m, h]; result ]; (* Computes the polar form of a monomial of degree h, *) (* w.r.t degree m, as a list {coeff, {monomial list}, coeff} *) (* where {monomial list} is a list of *) (* lists {i1, ... ih} indicating *) (* the indices of variables Ui of power 1 in a monomial *) cpolarlismon[{amon__}, m_] := Block[ {mon = {amon}, Ilist, III, setm, a, h, i, l1, I1, ii, result, Umon}, a = mon[[1]]; h = mon[[2]]; setm = Table[i, {i, m}]; result = {}; Ilist = subsets[setm, h]; l1 = Length[Ilist]; Do[ III = Ilist[[I1]]; Umon = {}; b = 1; Do[ ii = III[[i]]; Umon = Append[Umon, ii], {i, 1, Length[III]} ]; result = Join[result, {Umon}], {I1, 1, l1} ]; result = {a, result, Binomial[m, h]}; result ]; (* Computes the polar form of a polynomial of one variable *) (* w.r.t degree m, symbolically *) cpolarize[{cpoly__},m_] := Block[ {poly = {cpoly}, mon1, mon2, i, l, result}, result = 0; l = Length[poly]; Do[ mon1 = poly[[i]]; mon2 = cpolarmon[mon1, m]; result = result + mon2, {i, 1, l} ]; result ]; (* Computes the polar form of a polynomial of one variable *) (* w.r.t degree m, as a list of monomials *) cpolarizelis[{cpoly__},m_] := Block[ {poly = {cpoly}, mon1, mon2, i, l, result}, result = {}; l = Length[poly]; Do[ mon1 = poly[[i]]; mon2 = cpolarlismon[mon1, m]; result = Join[result, {mon2}], {i, 1, l} ]; result ]; (* Computes the polar form of a polynomial of two variables *) (* w.r.t total degree m, symbolically *) polarize[{cpoly__},m_] := Block[ {poly = {cpoly}, mon1, mon2, i, l, result}, result = 0; l = Length[poly]; Do[ mon1 = poly[[i]]; mon2 = polarmon[mon1, m]; result = result + mon2, {i, 1, l} ]; result ]; (* Computes the polar forms of a list of polynomials of two variables *) (* w.r.t total degree m, symbolically *) polarizesurf[{polyl__},m_] := Block[ {lpoly = {polyl}, poly1, polar, i, l, result}, result = {}; l = Length[lpoly]; Do[ poly1 = lpoly[[i]]; polar = polarize[poly1, m]; result = Append[result, polar], {i, 1, l} ]; result ]; (* Computes the polar form of a polynomial of two variables *) (* w.r.t bidegree (p, q), symbolically *) bipolarize[{cpoly__},p_, q_] := Block[ {poly = {cpoly}, mon1, mon2, i, l, result}, result = 0; l = Length[poly]; Do[ mon1 = poly[[i]]; mon2 = bipolarmon[mon1, p, q]; result = result + mon2, {i, 1, l} ]; result ]; (* Computes the polar form of a polynomial of two variables *) (* w.r.t total degree m, as a list of monomials *) polarizelis[{cpoly__},m_] := Block[ {poly = {cpoly}, mon1, mon2, i, l, result}, result = {}; l = Length[poly]; Do[ mon1 = poly[[i]]; mon2 = polarlismon[mon1, m]; result = Join[result, {mon2}], {i, 1, l} ]; result ]; (* Computes the polar form of a polynomial of two variables *) (* w.r.t bidegree (p, q), as a list of monomials *) bipolarizelis[{cpoly__}, p_, q_] := Block[ {poly = {cpoly}, mon1, mon2, i, l, result}, result = {}; l = Length[poly]; Do[ mon1 = poly[[i]]; mon2 = bipolarlismon[mon1, p, q]; result = Join[result, {mon2}], {i, 1, l} ]; result ]; (* Computes the polar forms of list of polynomials of two variables *) (* w.r.t bidegree (p, q), symbolically *) bipolarizesurf[{polyl__}, p_, q_] := Block[ {lpoly = {polyl}, poly1, polar, i, l, result}, result = {}; l = Length[lpoly]; Do[ poly1 = lpoly[[i]]; polar = bipolarize[poly1, p, q]; result = Append[result, polar], {i, 1, l} ]; result ]; (* To project down onto affine space, keeping the weights *) piomeg[{net__}] := Block[ {snet = {net}, newnet = {}, pt, newpt, j, l2, h, w}, l2 = Length[snet]; Do[ pt = snet[[j]]; w = Last[pt]; h = Drop[pt, -1]; If[w =!= 0 && w =!= 0., h = h/w]; newpt = Append[h, w]; newnet = Append[newnet,newpt], {j, 1, l2} ]; newnet ]; zpiomeg[{net__}] := Block[ {snet = {net}, newnet = {}, pt, newpt, j, l2, h, w}, l2 = Length[snet]; Do[ pt = snet[[j]]; w = Last[pt]; h = Drop[pt, -1]; If[Abs[w] > 10^(-20), h = h/w, w = 0]; newpt = Append[h, w]; newnet = Append[newnet,newpt], {j, 1, l2} ]; newnet ]; (* Computes the value of a polynomial in the variables *) (* U1, ... , Um, V1, ..., Vm, *) (* assuming that Ui, Vj take the values 0 or 1 *) evalpoly[{cpoly__}, {uu__}, {vv__}] := Block[ {poly = {cpoly}, U = {uu}, V = {vv}, mon1, mon2, monlis, result, i, j, ih, jk, l1, l2, ulis, vlis, a, c, val, monval, ii, jj, debug}, l1 = Length[poly]; result = 0; debug = {}; Do[ mon1 = poly[[i]]; a = mon1[[1]]; c = mon1[[3]]; monlis = mon1[[2]]; l2 = Length[monlis]; monval = 0; Do[ mon2 = monlis[[j]]; ulis = mon2[[1]]; vlis = mon2[[2]]; val = 1; Do[ih = ulis[[ii]]; If[U[[ih]] === 0, val = 0], {ii, 1, Length[ulis]} ]; Do[jk = vlis[[jj]]; If[V[[jk]] === 0, val = 0], {jj, 1, Length[vlis]} ]; monval = monval + val, {j, 1, l2} ]; monval = a * (monval/c); result = result + monval, {i, 1, l1} ]; result ]; (* Computes the value of a polynomial in the variables *) (* U1, ... , Um, V1, ..., Vm *) (* for arbitrary values of Ui, Vj *) gevalpoly[{cpoly__}, {uu__}, {vv__}] := Block[ {poly = {cpoly}, U = {uu}, V = {vv}, mon1, mon2, monlis, result, i, j, ih, jk, l1, l2, ulis, vlis, a, c, val, monval, ii, jj, debug}, l1 = Length[poly]; result = 0; debug = {}; Do[ mon1 = poly[[i]]; a = mon1[[1]]; c = mon1[[3]]; monlis = mon1[[2]]; l2 = Length[monlis]; monval = 0; Do[ mon2 = monlis[[j]]; ulis = mon2[[1]]; vlis = mon2[[2]]; val = 1; Do[ih = ulis[[ii]]; val = U[[ih]]*val, {ii, 1, Length[ulis]} ]; Do[jk = vlis[[jj]]; val = V[[jk]]*val, {jj, 1, Length[vlis]} ]; monval = monval + val, {j, 1, l2} ]; monval = a * (monval/c); result = result + monval, {i, 1, l1} ]; result ]; (* Computes the value of a polynomial in the variables *) (* U1, ... , Um, *) (* for arbitrary values of Ui *) gevalpoly1[{cpoly__}, {uu__}] := Block[ {poly = {cpoly}, U = {uu}, mon1, mon2, monlis, result, i, j, ih, l1, l2, ulis, a, c, val, monval, ii, debug}, l1 = Length[poly]; result = 0; debug = {}; Do[ mon1 = poly[[i]]; a = mon1[[1]]; c = mon1[[3]]; monlis = mon1[[2]]; l2 = Length[monlis]; monval = 0; Do[ ulis = monlis[[j]]; val = 1; Do[ih = ulis[[ii]]; val = U[[ih]]*val, {ii, 1, Length[ulis]} ]; monval = monval + val, {j, 1, l2} ]; monval = a * (monval/c); result = result + monval, {i, 1, l1} ]; result ]; (* computes the control polygon of a rational curve of degree m from some polar forms *) gcontpoly[{cpoly__},m_,r1_,s1_] := Block[ {poly = {cpoly}, poly1, poly2, poly3, poly4, i, r, s, l1, affmap1, affmap2, affmap3, affmap4, p1, p2, p3, p4, result, U, zo}, result = {}; l1 = Length[poly]; poly1 = poly[[1]]; poly2 = poly[[2]]; poly3 = poly[[3]]; If[l1 === 4, poly4 = poly[[4]]; affmap4 = cpolarizelis[poly4, m] ]; affmap1 = cpolarizelis[poly1, m]; affmap2 = cpolarizelis[poly2, m]; affmap3 = cpolarizelis[poly3, m]; Do[ U = {}; Do[ U = Append[U, r1], {r, 1, m - i} ]; Do[ U = Append[U, s1], {s, m + 1 - i, m} ]; (* Print["U = ", U]; *) p1 = gevalpoly1[affmap1, U]; p2 = gevalpoly1[affmap2, U]; p3 = gevalpoly1[affmap3, U]; If[l1 === 4, p4 = gevalpoly1[affmap4, U]; zo = {p1, p2, p3, p4}, zo = {p1, p2, p3} ]; result = Join[result, {zo}], {i, 0, m} ]; result = piomeg[result]; result ]; (* Computes a triangular control net of degree m from some polar forms, *) (* w.r.t. reference triangle (r, s, t) = ((1, 0, 0), (0, 1, 0), (0, 0, 1)) *) contnet[{cpoly__},m_] := Block[ {poly = {cpoly}, poly1, poly2, poly3, poly4, i, j, r, s, t, l1, affmap1, affmap2, affmap3, affmap4, p1, p2, p3, p4, result, U, V, zo, count}, result = {}; count = 0; l1 = Length[poly]; poly1 = poly[[1]]; poly2 = poly[[2]]; poly3 = poly[[3]]; If[l1 === 3, affmap4 = {{1, {{{}, {}}}, 1}}, poly4 = poly[[4]]; affmap4 = polarizelis[poly4, m] ]; Print[" poly4 polarized!"]; affmap1 = polarizelis[poly1, m]; Print[" poly1 polarized!"]; affmap2 = polarizelis[poly2, m]; Print[" poly2 polarized!"]; affmap3 = polarizelis[poly3, m]; Print[" poly3 polarized!"]; Do[ Do[ U = {}; V = {}; Do[ U = Append[U, 1]; V = Append[V, 0], {r, 1, i} ]; Do[ U = Append[U, 0]; V = Append[V, 1], {s, i + 1, i + j} ]; Do[ U = Append[U, 0]; V = Append[V, 0], {t, i + j + 1, m} ]; p1 = evalpoly[affmap1, U, V]; p2 = evalpoly[affmap2, U, V]; p3 = evalpoly[affmap3, U, V]; If[l1 === 3, p4 = 1, p4 = evalpoly[affmap4, U, V]]; zo = {p1, p2, p3, p4}; count = count + 1; Print[" count = ", count]; Print[" zo = ", zo]; result = Join[result, {zo}], {j, 0, m - i} ], {i, 0, m} ]; result = piomeg[result]; result ]; (* Computes a triangular control net of degree m from some polar forms, *) (* w.r.t. reference triangle (r, s, t) = ((1, 0, 0), (0, 1, 0), (0, 0, 1)) *) (* version for 4D surfaces *) contnet4[{cpoly__},m_] := Block[ {poly = {cpoly}, poly1, poly2, poly3, poly4, poly5, i, j, r, s, t, l1, affmap1, affmap2, affmap3, affmap4, affmap5, p1, p2, p3, p4, p5, result, U, V, zo, count}, result = {}; count = 0; l1 = Length[poly]; poly1 = poly[[1]]; poly2 = poly[[2]]; poly3 = poly[[3]]; poly4 = poly[[4]]; poly5 = poly[[5]]; affmap1 = polarizelis[poly1, m]; Print[" poly1 polarized!"]; affmap2 = polarizelis[poly2, m]; Print[" poly2 polarized!"]; affmap3 = polarizelis[poly3, m]; Print[" poly3 polarized!"]; affmap4 = polarizelis[poly4, m]; Print[" poly4 polarized!"]; affmap5 = polarizelis[poly5, m]; Print[" poly5 polarized!"]; Do[ Do[ U = {}; V = {}; Do[ U = Append[U, 1]; V = Append[V, 0], {r, 1, i} ]; Do[ U = Append[U, 0]; V = Append[V, 1], {s, i + 1, i + j} ]; Do[ U = Append[U, 0]; V = Append[V, 0], {t, i + j + 1, m} ]; p1 = evalpoly[affmap1, U, V]; p2 = evalpoly[affmap2, U, V]; p3 = evalpoly[affmap3, U, V]; p4 = evalpoly[affmap4, U, V]; p5 = evalpoly[affmap5, U, V]; zo = {p1, p2, p3, p4, p5}; count = count + 1; Print[" count = ", count]; Print[" zo = ", zo]; result = Join[result, {zo}], {j, 0, m - i} ], {i, 0, m} ]; result = piomeg[result]; result ]; (* Computes a triangular control net of degree m from some polar forms, *) (* w.r.t. reference triangle *) (* (r, s, t) = ((r1, r2, r3), (s1, s2, s3), (t1, t2, t3)) *) gcontnet[{cpoly__},m_,reftrig_] := Block[ {poly = {cpoly}, poly1, poly2, poly3, poly4, i, j, r, s, t, l1, affmap1, affmap2, affmap3, affmap4, p1, p2, p3, p4, result, U, V, zo, r1, r2, s1, s2, t1, t2}, result = {}; r1 = reftrig[[1,1]]; r2 = reftrig[[1,2]]; s1 = reftrig[[2,1]]; s2 = reftrig[[2,2]]; t1 = reftrig[[3,1]]; t2 = reftrig[[3,2]]; l1 = Length[poly]; poly1 = poly[[1]]; poly2 = poly[[2]]; poly3 = poly[[3]]; If[l1 === 3, affmap4 = {{1, {{{}, {}}}, 1}}, poly4 = poly[[4]]; affmap4 = polarizelis[poly4, m] ]; affmap1 = polarizelis[poly1, m]; affmap2 = polarizelis[poly2, m]; affmap3 = polarizelis[poly3, m]; Do[ Do[ U = {}; V = {}; Do[ U = Append[U, r1]; V = Append[V, r2], {r, 1, i} ]; Do[ U = Append[U, s1]; V = Append[V, s2], {s, i + 1, i + j} ]; Do[ U = Append[U, t1]; V = Append[V, t2], {t, i + j + 1, m} ]; p1 = gevalpoly[affmap1, U, V]; p2 = gevalpoly[affmap2, U, V]; p3 = gevalpoly[affmap3, U, V]; If[l1 === 3, p4 = 1, p4 = gevalpoly[affmap4, U, V]]; zo = {p1, p2, p3, p4}; result = Join[result, {zo}], {j, 0, m - i} ], {i, 0, m} ]; result = piomeg[result]; result ]; (* Computes a rectangular control net of bidegree (p, q) from some polar forms, *) (* w.r.t. affine frames (0, 1), (0, 1) *) bicontnet[{cpoly__},p_,q_] := Block[ {poly = {cpoly}, poly1, poly2, poly3, poly4, i, j, r, s, l1, affmap1, affmap2, affmap3, affmap4, p1, p2, p3, p4, result, U, V, zo}, result = {}; l1 = Length[poly]; poly1 = poly[[1]]; poly2 = poly[[2]]; poly3 = poly[[3]]; If[l1 === 3, affmap4 = {{1, {{{}, {}}}, 1}}, poly4 = poly[[4]]; affmap4 = bipolarizelis[poly4, p, q] ]; affmap1 = bipolarizelis[poly1, p, q]; affmap2 = bipolarizelis[poly2, p, q]; affmap3 = bipolarizelis[poly3, p, q]; Do[ Do[ U = {}; V = {}; Do[ U = Append[U, 0], {r, 1, p - i} ]; Do[ U = Append[U, 1], {s, p + 1 - i, p} ]; Do[ V = Append[V, 0], {r, 1, q - j} ]; Do[ V = Append[V, 1], {s, q + 1 - j, q} ]; p1 = evalpoly[affmap1, U, V]; p2 = evalpoly[affmap2, U, V]; p3 = evalpoly[affmap3, U, V]; If[l1 === 3, p4 = 1, p4 = evalpoly[affmap4, U, V]]; zo = {p1, p2, p3, p4}; result = Join[result, {zo}], {j, 0, q} ], {i, 0, p} ]; result = piomeg[result]; result ]; (* Computes a rectangular control net of bidegree (p, q) from some polar forms, *) (* w.r.t. affine frames ((r1, s1), (r2, s2) *) gbicontnet[{cpoly__},p_,q_,r1_,s1_,r2_,s2_] := Block[ {poly = {cpoly}, poly1, poly2, poly3, poly4, i, j, r, s, l1, affmap1, affmap2, affmap3, affmap4, p1, p2, p3, p4, result, U, V, zo}, result = {}; l1 = Length[poly]; poly1 = poly[[1]]; poly2 = poly[[2]]; poly3 = poly[[3]]; If[l1 === 3, affmap4 = {{1, {{{}, {}}}, 1}}, poly4 = poly[[4]]; affmap4 = bipolarizelis[poly4, p, q] ]; affmap1 = bipolarizelis[poly1, p, q]; affmap2 = bipolarizelis[poly2, p, q]; affmap3 = bipolarizelis[poly3, p, q]; Do[ Do[ U = {}; V = {}; Do[ U = Append[U, r1], {r, 1, p - i} ]; Do[ U = Append[U, s1], {s, p + 1 - i, p} ]; Do[ V = Append[V, r2], {r, 1, q - j} ]; Do[ V = Append[V, s2], {s, q + 1 - j, q} ]; (* Print["U = ", U]; *) (* Print["V = ", V]; *) p1 = gevalpoly[affmap1, U, V]; p2 = gevalpoly[affmap2, U, V]; p3 = gevalpoly[affmap3, U, V]; If[l1 === 3, p4 = 1, p4 = gevalpoly[affmap4, U, V]]; zo = {p1, p2, p3, p4}; result = Join[result, {zo}], {j, 0, q} ], {i, 0, p} ]; result = piomeg[result]; result ]; (* Computes a rectangular control net of bidegree (p, q) from some polar forms, *) (* w.r.t. affine frames ((r1, s1), (r2, s2), IN THE HAT SPACE *) (* does not divide by the weight *) ghbicontnet[{cpoly__},p_,q_,r1_,s1_,r2_,s2_] := Block[ {poly = {cpoly}, poly1, poly2, poly3, poly4, i, j, r, s, l1, affmap1, affmap2, affmap3, affmap4, p1, p2, p3, p4, result, U, V, zo}, result = {}; l1 = Length[poly]; poly1 = poly[[1]]; poly2 = poly[[2]]; poly3 = poly[[3]]; If[l1 === 3, affmap4 = {{1, {{{}, {}}}, 1}}, poly4 = poly[[4]]; affmap4 = bipolarizelis[poly4, p, q] ]; affmap1 = bipolarizelis[poly1, p, q]; affmap2 = bipolarizelis[poly2, p, q]; affmap3 = bipolarizelis[poly3, p, q]; Do[ Do[ U = {}; V = {}; Do[ U = Append[U, r1], {r, 1, p - i} ]; Do[ U = Append[U, s1], {s, p + 1 - i, p} ]; Do[ V = Append[V, r2], {r, 1, q - j} ]; Do[ V = Append[V, s2], {s, q + 1 - j, q} ]; (* Print["U = ", U]; *) (* Print["V = ", V]; *) p1 = gevalpoly[affmap1, U, V]; p2 = gevalpoly[affmap2, U, V]; p3 = gevalpoly[affmap3, U, V]; If[l1 === 3, p4 = 1, p4 = gevalpoly[affmap4, U, V]]; zo = {p1, p2, p3, p4}; result = Join[result, {zo}], {j, 0, q} ], {i, 0, p} ]; result ]; (* Computes the polar value of all monomials of degree k < = m *) (* times Binomial[m, k]] *) (* wrt val = (r, ..., r, s, ..., s) *) fastpol[{val__}, m_] := Block[ {v = {val}, s, ns, i, j, x}, s = {{1, 0}}; Do[ Do[ If[j === 1, ns = {1} , x = s[[i, j]] + v[[i]]*s[[i, j - 1]]; ns = Append[ns, x] ], {j, 1, i + 1} ]; ns = Append[ns, 0]; s = Append[s, ns], {i, 1, m} ]; s ]; (* Computes the polar value of all monomials of degree k < = m *) (* wrt val = (r, ..., r, s, ..., s) *) yfastpol[{val__}, m_] := Block[ {v = {val}, s, ns, i, j, x}, s = {{1, 0}}; Do[ Do[ If[j === 1, ns = {1} , x = ((i + 1 - j) * s[[i, j]] + (j - 1)*v[[i]]*s[[i, j - 1]])/i; ns = Append[ns, x] ], {j, 1, i + 1} ]; ns = Append[ns, 0]; s = Append[s, ns], {i, 1, m} ]; s ]; (* Computes the polar value of a monomial of degree k < = m *) (* wrt val = (r, ..., r, s, ..., s) *) zfastpol[{val__}, m_, k_] := Block[ {v = {val}, res, s, ns, i, j, x}, s = {1, 0}; Do[ Do[ If[j === 1, ns = {1} , x = s[[j]] + v[[i]]*s[[j - 1]]; ns = Append[ns, x] ], {j, 1, i + 1} ]; ns = Append[ns, 0]; s = ns, {i, 1, m} ]; If[0 <= k && k <= m, res = s[[k + 1]]/Binomial[m, k], res = 0]; res ]; (* Computes the polar value of a polynomial of one variable *) (* w.r.t degree m, for val = (r, ..., r, s, ..., s) *) yfastpolarval[{inpoly__}, {val__},m_] := Block[ {poly = {inpoly}, v = {val}, mon, l, s, a, h, x, i, result}, l = Length[poly]; result = 0; s = yfastpol[v, m]; Do[ mon = poly[[i]]; a = mon[[1]]; h = mon[[2]]; If[a === 0, x = 0, x = a*s[[m + 1, h + 1]]]; result = result + x, {i, 1, l} ]; result ]; (* Computes the polar value of a polynomial of one variable *) (* w.r.t degree m, for val = (r, ..., r, s, ..., s) *) fastpolarval[{inpoly__}, {val__},m_] := Block[ {poly = {inpoly}, v = {val}, mon, l, s, a, h, x, i, result}, l = Length[poly]; result = 0; s = fastpol[v, m]; Do[ mon = poly[[i]]; a = mon[[1]]; h = mon[[2]]; If[a === 0, x = 0, x = a*s[[m + 1, h + 1]]/Binomial[m, h]]; result = result + x, {i, 1, l} ]; result ]; (* Computes the polar value of a polynomial of one variable *) (* w.r.t degree m, for val = (r, ..., r, s, ..., s) *) zfastpolarval[{inpoly__}, {val__},m_] := Block[ {poly = {inpoly}, v = {val}, mon, l, a, h, x, i, result}, l = Length[poly]; result = 0; Do[ mon = poly[[i]]; a = mon[[1]]; h = mon[[2]]; If[a === 0, x = 0, x = a*zfastpol[v, m, h]]; result = result + x, {i, 1, l} ]; result ]; (* computes the control polygon of a rational curve of degree m *) (* fast method using a recurrence relation *) fcontpoly[{cpoly__},m_,r1_,s1_] := Block[ {poly = {cpoly}, poly1, poly2, poly3, poly4, i, r, s, l1, p1, p2, p3, p4, result, U, zo}, result = {}; l1 = Length[poly]; poly1 = poly[[1]]; poly2 = poly[[2]]; poly3 = poly[[3]]; If[l1 === 4, poly4 = poly[[4]] ]; Do[ U = {}; Do[ U = Append[U, r1], {r, 1, m - i} ]; Do[ U = Append[U, s1], {s, m + 1 - i, m} ]; (* Print["U = ", U]; *) p1 = fastpolarval[poly1, U, m]; p2 = fastpolarval[poly2, U, m]; p3 = fastpolarval[poly3, U, m]; If[l1 === 4, p4 = fastpolarval[poly4, U, m]; zo = {p1, p2, p3, p4}, zo = {p1, p2, p3} ]; Print[" point = ", zo]; result = Join[result, {zo}], {i, 0, m} ]; result = piomeg[result]; result ]; (* computes the control polygon of a rational curve of degree m *) (* fast method using a recurrence relation *) yfcontpoly[{cpoly__},m_,r1_,s1_] := Block[ {poly = {cpoly}, poly1, poly2, poly3, poly4, i, r, s, l1, p1, p2, p3, p4, result, U, zo}, result = {}; l1 = Length[poly]; poly1 = poly[[1]]; poly2 = poly[[2]]; poly3 = poly[[3]]; If[l1 === 4, poly4 = poly[[4]] ]; Do[ U = {}; Do[ U = Append[U, r1], {r, 1, m - i} ]; Do[ U = Append[U, s1], {s, m + 1 - i, m} ]; (* Print["U = ", U]; *) p1 = yfastpolarval[poly1, U, m]; p2 = yfastpolarval[poly2, U, m]; p3 = yfastpolarval[poly3, U, m]; If[l1 === 4, p4 = yfastpolarval[poly4, U, m]; zo = {p1, p2, p3, p4}, zo = {p1, p2, p3} ]; Print[" point = ", zo]; result = Join[result, {zo}], {i, 0, m} ]; result = piomeg[result]; result ]; (* Computes all the polar values of monomials U^hV^k where h + k < = m *) (* wrt u = (r1,..., r1, s1, ..., s1, t1, ..., t1), *) (* v = (r2,..., r2, s2, ..., s2, t2, ..., t2) *) fastpol2[{uu__}, {vv__}, m_] := Block[ {u = {uu}, v = {vv}, res, s, nl, ns, i, j, k, x}, s = {{{0, 0, 0}, {0, 0, 0}, {0, 0}}, {{0, 0, 0}, {0, 1, 0}, {0, 0}}}; Do[ ns = {}; Do[ nl = {0}; Do[ If[j === 1, nl = Append[nl, 0] , x = s[[i + 1, j, k + 1]] + u[[i]]*s[[i + 1, j - 1, k + 1]] + v[[i]]*s[[i + 1, j, k]] + u[[i]]*v[[i]]*s[[i, j - 1, k]]; nl = Append[nl, x] ], {k, 1, i + 3 - j} ]; If[j =!= 1, nl = Append[nl, 0]]; ns = Append[ns, nl], {j, 1, i + 2} ]; ns = Append[ns, {0, 0}]; s = Append[s, ns], {i, 1, m} ]; s ]; (* Computes the polar value of a monomial U^hV^k where h + k < = m *) (* wrt U = (r1,..., r1, s1, ..., s1, t1, ..., t1), *) (* V = (r2,..., r2, s2, ..., s2, t2, ..., t2) *) zfastpol2[{uu__}, {vv__}, m_, hh_, kk_] := Block[ {u = {uu}, v = {vv}, res, s, nl, ns, i, j, k, x}, s = {{{0, 0, 0}, {0, 0, 0}, {0, 0}}, {{0, 0, 0}, {0, 1, 0}, {0, 0}}}; Do[ ns = {}; Do[ nl = {0}; Do[ If[j === 1, nl = Append[nl, 0] , x = s[[i + 1, j, k + 1]] + u[[i]]*s[[i + 1, j - 1, k + 1]] + v[[i]]*s[[i + 1, j, k]] + u[[i]]*v[[i]]*s[[i, j - 1, k]]; nl = Append[nl, x] ], {k, 1, i + 3 - j} ]; If[j =!= 1, nl = Append[nl, 0]]; ns = Append[ns, nl], {j, 1, i + 2} ]; ns = Append[ns, {0, 0}]; s = Append[s, ns]; Print[" s = ", s], {i, 1, m} ]; If[0 <= hh && 0 <= kk && hh + kk <= m, res = s[[m + 2, hh + 2, kk + 2]]/Multinomial[hh, kk, m - hh - kk], res = 0 ]; res ]; (* Computes the value of the polar form of a polynomial in 2 variables *) (* w.r.t. total degree m *) (* For the values U1, ... , Um, V1, ..., Vm *) (* for arbitrary values of Ui, Vj *) fastpolartval[{inpoly__}, {uu__}, {vv__}, m_] := Block[ {poly = {inpoly}, U = {uu}, V = {vv}, s, l1, mon, result, i, a, h, k, d, monval}, l1 = Length[poly]; result = 0; s = fastpol2[U, V, m]; Do[ mon = poly[[i]]; a = mon[[1]]; h = mon[[2]]; k = mon[[3]]; d = h + k; If[a === 0, monval = 0, monval = a*s[[m + 2, h + 2, k + 2]]/Multinomial[h, k, m - d] ]; result = result + monval, {i, 1, l1} ]; result ]; (* Computes a triangular control net of degree m *) (* w.r.t. reference triangle *) (* (r, s, t) = ((r1, r2, r3), (s1, s2, s3), (t1, t2, t3)) *) (* fast method using a recurrence relation *) fgcontnet[{cpoly__},m_,reftrig_] := Block[ {poly = {cpoly}, poly1, poly2, poly3, poly4, i, j, r, s, t, l1, affmap1, affmap2, affmap3, affmap4, p1, p2, p3, p4, result, U, V, zo, r1, r2, s1, s2, t1, t2}, result = {}; r1 = reftrig[[1,1]]; r2 = reftrig[[1,2]]; s1 = reftrig[[2,1]]; s2 = reftrig[[2,2]]; t1 = reftrig[[3,1]]; t2 = reftrig[[3,2]]; l1 = Length[poly]; poly1 = poly[[1]]; poly2 = poly[[2]]; poly3 = poly[[3]]; If[l1 === 3, poly4 = {{1, {{{}, {}}}, 1}}, poly4 = poly[[4]]; ]; Do[ Do[ U = {}; V = {}; Do[ U = Append[U, r1]; V = Append[V, r2], {r, 1, i} ]; Do[ U = Append[U, s1]; V = Append[V, s2], {s, i + 1, i + j} ]; Do[ U = Append[U, t1]; V = Append[V, t2], {t, i + j + 1, m} ]; p1 = fastpolartval[poly1, U, V, m]; p2 = fastpolartval[poly2, U, V, m]; p3 = fastpolartval[poly3, U, V, m]; If[l1 === 3, p4 = 1, p4 = fastpolartval[poly4, U, V, m]]; zo = {p1, p2, p3, p4}; Print[" point = ", zo]; result = Join[result, {zo}], {j, 0, m - i} ], {i, 0, m} ]; result = piomeg[result]; result ]; reftrig = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}; (* cnet = fgcontnet[enpoly, 9, reftrig] *)