読者です 読者をやめる 読者になる 読者になる

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)}