Algorithms for Trigonometric Curves (Simplifcation,Implicitization,Parameterization)
H.HONG,J.SCHICHO.Jurnal of Symbolic Computation 26(1998),279-300.
(* 次数 *)
TrigDegree[ft_, t_] := Exponent[TrigToExp@ft, Exp[I t]];
(* 主係数 *)
TrigLeadingCoefficient[ft_, t_] :=
Coefficient[TrigReduce@ft, Through[{Cos, Sin}[TrigDegree[ft, t] t]]];
(* 擬除法（主項のみ） *)
TrigQuotientRemainder0[ft_, gt_, t_] :=
Block[{X = TrigLeadingCoefficient[ft, t][[1]],
Y = TrigLeadingCoefficient[ft, t][[2]],
x = TrigLeadingCoefficient[gt, t][[1]],
y = TrigLeadingCoefficient[gt, t][[2]],
mn = TrigDegree[ft, t] - TrigDegree[gt, t],
q = {(2 (x X + y Y))/(x^2 + y^2), -((2 (X y - x Y))/(x^2 + y^2))}
.{Cos[mn t], Sin[mn t]}}, TrigReduce /@ {q, ft - gt*q}];
(* 擬除法 *)
TrigQuotientRemainder1[ft_, gt_, t_] :=
Block[{p, r = ft, n = 0, q, c}, p[0] = 0;
While[TrigDegree[r, t] > TrigDegree[gt, t],
n++; {p[n], r} = TrigQuotientRemainder0[r, gt, t]];
q = Sum[p[k], {k, 0, n}] + c /.
If[TrigDegree[r, t] == TrigDegree[gt, t] &&
Det[{TrigLeadingCoefficient[r, t],
TrigLeadingCoefficient[gt, t]}] == 0,
Solve[TrigLeadingCoefficient[r, t] ==
c TrigLeadingCoefficient[gt, t], c][[1]], c -> 0];
TrigReduce /@ {q, ft - gt*q}];
(* 除法 *)
TrigQuotientRemainder[ft_, gt_, t_] :=
Block[{qr = TrigQuotientRemainder1[TrigExpand@ft, TrigExpand@gt, t]},
If[TrigDegree[qr[[2]], t] ==
TrigDegree[TrigExpand@gt, t], {"nonQR"}, qr]];
(* 三角多項式から（代数）多項式 *)
TrigToPoly2[ft_, gt_, t_, par_: s] :=
Block[{S = {{1, TrigExpand@ft}, {2, TrigExpand@gt}}, i, j, M = {}, n = 2, tqr},
While[1 < Length@S,
If[TrigDegree[S[[1, 2]], t] <
TrigDegree[S[[2, 2]], t], {S[[1]], S[[2]]} = {S[[2]],
S[[1]]}]; {i, j} = {S[[1, 1]], S[[2, 1]]};
If[TrigDegree[S[[2, 2]], t] <= 0, AppendTo[M, u[j] == S[[2, 2]]];
S = Delete[S, 2];,
tqr = TrigQuotientRemainder[S[[1, 2]], S[[2, 2]], t];
If[tqr == {"nonQR"}, Return["nonP"],
S = Join[S, {{n + 1, tqr[[1]]}, {n + 2, tqr[[2]]}}];
AppendTo[M, u[i] == u[j]*u[n + 1] + u[n + 2]]; S = Delete[S, 1];
n = n + 2];(*Print[{S,M}]*)]];
Return[{u[1], u[2],
MinValue[S[[1, 2]], t] <= par <= MaxValue[S[[1, 2]], t]} /.
Solve[AppendTo[M, u[S[[1, 1]]] == Glaisher]][[1]] /. Glaisher -> par]];
(* （代数）多項式から陰関数 *)
PolyToImp[Xs_, Ys_, s_, Rs_, X_: x, Y_: y] :=
Block[{}, {XY, sXY} =
Take[If[Exponent[Xs, s] >= Exponent[Ys, s],
SubresultantPolynomials[Xs - X, Ys - Y, s],
SubresultantPolynomials[Ys - Y, Xs - X, s]], 2];
Return[{XY, Rs /. Solve[sXY == 0, s][[1]]}]];
(* 三角多項式から陰関数 *)
TrigToImp[ft_, gt_, t_, X_: x, Y_: y] :=
Block[{s, ttp = TrigToPoly2[ft, gt, t, s], z},
If[ttp === "nonP",
Factor[Resultant @@
AppendTo[
Numerator[
Factor[TrigToExp[{ft - X, gt - Y}] /. t -> -I Log[z]]], z]],
PolyToImp[ttp[[1]], ttp[[2]], s, ttp[[3]], X, Y]]];
(* 実行例 *)
In[5]:= TrigToImp[Cos[t], Sin[t], t]
Out[5]= 4 (-1 + x^2 + y^2)
In[6]:= TrigToImp[Cos[3 t], Cos[2 t], t]
Out[6]= {-4 + 8 x^2 + 12 y - 16 y^3, -1 <= x/(-1 + 2 y) <= 1}
In[7]:= TrigToImp[-4 Cos[t] - 6 Cos[3 t] + 6 Cos[5 t] + 3 Cos[7 t] +
Cos[9 t], -56 Cos[2 t] - 33 Cos[4 t] - 12 Cos[6 t] + 10 Cos[8 t] +
4 Cos[10 t] + Cos[12 t] + 10 Sin[t], t] // AbsoluteTiming
Out[7]= {7.566019,
262144 (757789974567443369603663333556224 -
41469986362222101588417046397084 x^2 +
973732507123909398819975575344 x^4 -
12800507764943231018833032000 x^6 +
103171447469438771599603360 x^8 - 526200579837205609181824 x^10 +
1703918904823696794624 x^12 - 3483525971359964864 x^14 +
4522213285333760 x^16 - 3715664486400 x^18 + 1872341504 x^20 -
528384 x^22 + 64 x^24 - 5828663013423184303908390785024 y -
246087078955285245633237210624 x^2 y +
18034878313025376995333650560 x^4 y -
380597934592624427607223040 x^6 y +
4028008142215886344969840 x^8 y -
24011438511810778040256 x^10 y + 82420815154523057664 x^12 y -
161432650926010560 x^14 y + 183623608143360 x^16 y -
120319344640 x^18 y + 42199296 x^20 y - 6144 x^22 y -
2055210476445866013709115366400 y^2 +
95705079672674210708853790080 x^2 y^2 -
1703741545381291262689908320 x^4 y^2 +
13948322966725935765598400 x^6 y^2 -
41348895479757190234880 x^8 y^2 -
111498097826565338880 x^10 y^2 + 1037637196794232320 x^12 y^2 -
2361575848160400 x^14 y^2 + 2419860055680 x^16 y^2 -
1175838720 x^18 y^2 + 220800 x^20 y^2 -
4682739992223613035276828672 y^3 +
1447430440681944563202784960 x^2 y^3 -
48700190263096619154696480 x^4 y^3 +
658190578375288641832320 x^6 y^3 -
4226070069531954275840 x^8 y^3 + 12262179782564638592 x^10 y^3 -
11192575145396480 x^12 y^3 - 2042422640640 x^14 y^3 +
6753246720 x^16 y^3 - 1925120 x^18 y^3 - 192 x^20 y^3 +
2373329561181171131151249920 y^4 -
80292526116290952845368960 x^2 y^4 +
752705909069344935896820 x^4 y^4 +
1419513402065161964560 x^6 y^4 - 56612895009215650880 x^8 y^4 +
274284187557180000 x^10 y^4 - 402384531820800 x^12 y^4 +
252017172480 x^14 y^4 - 81318720 x^16 y^4 + 15360 x^18 y^4 +
27165079892249234553595136 y^5 -
2062365201225957172020416 x^2 y^5 +
41317273789301336389376 x^4 y^5 - 307483552751047797760 x^6 y^5 +
607695931383793920 x^8 y^5 + 1612077252246096 x^10 y^5 -
3417094851456 x^12 y^5 + 2037313536 x^14 y^5 - 429120 x^16 y^5 -
1318886866729352346975936 y^6 + 20650460373791821605504 x^2 y^6 +
224193950107843583360 x^4 y^6 - 5238613665598379520 x^6 y^6 +
24297798913309760 x^8 y^6 - 15123826876416 x^10 y^6 -
3375097856 x^12 y^6 + 2549760 x^14 y^6 + 240 x^16 y^6 -
28497642159377727554560 y^7 + 1172900780403546648640 x^2 y^7 -
11777614983731396160 x^4 y^7 + 6112547155994240 x^6 y^7 +
253519508878080 x^8 y^7 - 261107712000 x^10 y^7 +
88104960 x^12 y^7 - 15360 x^14 y^7 + 250953353308865668032 y^8 +
7119252912927910560 x^2 y^8 - 220306445322132480 x^4 y^8 +
1131223929041400 x^6 y^8 + 755024495760 x^8 y^8 -
1183684608 x^10 y^8 + 306240 x^12 y^8 +
12551823115420924480 y^9 - 240137471011187440 x^2 y^9 -
396733299880320 x^4 y^9 + 14655086722560 x^6 y^9 -
5472388480 x^8 y^9 - 276480 x^10 y^9 - 160 x^12 y^9 +
69023763862465520 y^10 - 4482397934459320 x^2 y^10 +
30656372687520 x^4 y^10 + 82150922240 x^6 y^10 -
47199840 x^8 y^10 + 7680 x^10 y^10 - 2072074443150080 y^11 -
11940349153980 x^2 y^11 + 438503739600 x^4 y^11 +
142955520 x^6 y^11 - 91680 x^8 y^11 - 35582159985280 y^12 +
467718920960 x^2 y^12 + 2802564960 x^4 y^12 - 508160 x^6 y^12 +
60 x^8 y^12 - 90010893520 y^13 + 6762393600 x^2 y^13 +
8593200 x^4 y^13 - 1920 x^6 y^13 + 3125601660 y^14 +
43269120 x^2 y^14 + 7560 x^4 y^14 + 42331008 y^15 +
141312 x^2 y^15 - 12 x^4 y^15 + 253932 y^16 + 192 x^2 y^16 +
780 y^17 + y^18)}
