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

Algorithms for Trigonometric Curves

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]]];