prde.py 51 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332
  1. """
  2. Algorithms for solving Parametric Risch Differential Equations.
  3. The methods used for solving Parametric Risch Differential Equations parallel
  4. those for solving Risch Differential Equations. See the outline in the
  5. docstring of rde.py for more information.
  6. The Parametric Risch Differential Equation problem is, given f, g1, ..., gm in
  7. K(t), to determine if there exist y in K(t) and c1, ..., cm in Const(K) such
  8. that Dy + f*y == Sum(ci*gi, (i, 1, m)), and to find such y and ci if they exist.
  9. For the algorithms here G is a list of tuples of factions of the terms on the
  10. right hand side of the equation (i.e., gi in k(t)), and Q is a list of terms on
  11. the right hand side of the equation (i.e., qi in k[t]). See the docstring of
  12. each function for more information.
  13. """
  14. import itertools
  15. from functools import reduce
  16. from sympy.core import Dummy, ilcm, Add, Mul, Pow, S
  17. from sympy.integrals.rde import (order_at, order_at_oo, weak_normalizer,
  18. bound_degree)
  19. from sympy.integrals.risch import (gcdex_diophantine, frac_in, derivation,
  20. residue_reduce, splitfactor, residue_reduce_derivation, DecrementLevel,
  21. recognize_log_derivative)
  22. from sympy.polys import Poly, lcm, cancel, sqf_list
  23. from sympy.polys.polymatrix import PolyMatrix as Matrix
  24. from sympy.solvers import solve
  25. zeros = Matrix.zeros
  26. eye = Matrix.eye
  27. def prde_normal_denom(fa, fd, G, DE):
  28. """
  29. Parametric Risch Differential Equation - Normal part of the denominator.
  30. Explanation
  31. ===========
  32. Given a derivation D on k[t] and f, g1, ..., gm in k(t) with f weakly
  33. normalized with respect to t, return the tuple (a, b, G, h) such that
  34. a, h in k[t], b in k<t>, G = [g1, ..., gm] in k(t)^m, and for any solution
  35. c1, ..., cm in Const(k) and y in k(t) of Dy + f*y == Sum(ci*gi, (i, 1, m)),
  36. q == y*h in k<t> satisfies a*Dq + b*q == Sum(ci*Gi, (i, 1, m)).
  37. """
  38. dn, ds = splitfactor(fd, DE)
  39. Gas, Gds = list(zip(*G))
  40. gd = reduce(lambda i, j: i.lcm(j), Gds, Poly(1, DE.t))
  41. en, es = splitfactor(gd, DE)
  42. p = dn.gcd(en)
  43. h = en.gcd(en.diff(DE.t)).quo(p.gcd(p.diff(DE.t)))
  44. a = dn*h
  45. c = a*h
  46. ba = a*fa - dn*derivation(h, DE)*fd
  47. ba, bd = ba.cancel(fd, include=True)
  48. G = [(c*A).cancel(D, include=True) for A, D in G]
  49. return (a, (ba, bd), G, h)
  50. def real_imag(ba, bd, gen):
  51. """
  52. Helper function, to get the real and imaginary part of a rational function
  53. evaluated at sqrt(-1) without actually evaluating it at sqrt(-1).
  54. Explanation
  55. ===========
  56. Separates the even and odd power terms by checking the degree of terms wrt
  57. mod 4. Returns a tuple (ba[0], ba[1], bd) where ba[0] is real part
  58. of the numerator ba[1] is the imaginary part and bd is the denominator
  59. of the rational function.
  60. """
  61. bd = bd.as_poly(gen).as_dict()
  62. ba = ba.as_poly(gen).as_dict()
  63. denom_real = [value if key[0] % 4 == 0 else -value if key[0] % 4 == 2 else 0 for key, value in bd.items()]
  64. denom_imag = [value if key[0] % 4 == 1 else -value if key[0] % 4 == 3 else 0 for key, value in bd.items()]
  65. bd_real = sum(r for r in denom_real)
  66. bd_imag = sum(r for r in denom_imag)
  67. num_real = [value if key[0] % 4 == 0 else -value if key[0] % 4 == 2 else 0 for key, value in ba.items()]
  68. num_imag = [value if key[0] % 4 == 1 else -value if key[0] % 4 == 3 else 0 for key, value in ba.items()]
  69. ba_real = sum(r for r in num_real)
  70. ba_imag = sum(r for r in num_imag)
  71. ba = ((ba_real*bd_real + ba_imag*bd_imag).as_poly(gen), (ba_imag*bd_real - ba_real*bd_imag).as_poly(gen))
  72. bd = (bd_real*bd_real + bd_imag*bd_imag).as_poly(gen)
  73. return (ba[0], ba[1], bd)
  74. def prde_special_denom(a, ba, bd, G, DE, case='auto'):
  75. """
  76. Parametric Risch Differential Equation - Special part of the denominator.
  77. Explanation
  78. ===========
  79. Case is one of {'exp', 'tan', 'primitive'} for the hyperexponential,
  80. hypertangent, and primitive cases, respectively. For the hyperexponential
  81. (resp. hypertangent) case, given a derivation D on k[t] and a in k[t],
  82. b in k<t>, and g1, ..., gm in k(t) with Dt/t in k (resp. Dt/(t**2 + 1) in
  83. k, sqrt(-1) not in k), a != 0, and gcd(a, t) == 1 (resp.
  84. gcd(a, t**2 + 1) == 1), return the tuple (A, B, GG, h) such that A, B, h in
  85. k[t], GG = [gg1, ..., ggm] in k(t)^m, and for any solution c1, ..., cm in
  86. Const(k) and q in k<t> of a*Dq + b*q == Sum(ci*gi, (i, 1, m)), r == q*h in
  87. k[t] satisfies A*Dr + B*r == Sum(ci*ggi, (i, 1, m)).
  88. For case == 'primitive', k<t> == k[t], so it returns (a, b, G, 1) in this
  89. case.
  90. """
  91. # TODO: Merge this with the very similar special_denom() in rde.py
  92. if case == 'auto':
  93. case = DE.case
  94. if case == 'exp':
  95. p = Poly(DE.t, DE.t)
  96. elif case == 'tan':
  97. p = Poly(DE.t**2 + 1, DE.t)
  98. elif case in ('primitive', 'base'):
  99. B = ba.quo(bd)
  100. return (a, B, G, Poly(1, DE.t))
  101. else:
  102. raise ValueError("case must be one of {'exp', 'tan', 'primitive', "
  103. "'base'}, not %s." % case)
  104. nb = order_at(ba, p, DE.t) - order_at(bd, p, DE.t)
  105. nc = min([order_at(Ga, p, DE.t) - order_at(Gd, p, DE.t) for Ga, Gd in G])
  106. n = min(0, nc - min(0, nb))
  107. if not nb:
  108. # Possible cancellation.
  109. if case == 'exp':
  110. dcoeff = DE.d.quo(Poly(DE.t, DE.t))
  111. with DecrementLevel(DE): # We are guaranteed to not have problems,
  112. # because case != 'base'.
  113. alphaa, alphad = frac_in(-ba.eval(0)/bd.eval(0)/a.eval(0), DE.t)
  114. etaa, etad = frac_in(dcoeff, DE.t)
  115. A = parametric_log_deriv(alphaa, alphad, etaa, etad, DE)
  116. if A is not None:
  117. Q, m, z = A
  118. if Q == 1:
  119. n = min(n, m)
  120. elif case == 'tan':
  121. dcoeff = DE.d.quo(Poly(DE.t**2 + 1, DE.t))
  122. with DecrementLevel(DE): # We are guaranteed to not have problems,
  123. # because case != 'base'.
  124. betaa, alphaa, alphad = real_imag(ba, bd*a, DE.t)
  125. betad = alphad
  126. etaa, etad = frac_in(dcoeff, DE.t)
  127. if recognize_log_derivative(Poly(2, DE.t)*betaa, betad, DE):
  128. A = parametric_log_deriv(alphaa, alphad, etaa, etad, DE)
  129. B = parametric_log_deriv(betaa, betad, etaa, etad, DE)
  130. if A is not None and B is not None:
  131. Q, s, z = A
  132. # TODO: Add test
  133. if Q == 1:
  134. n = min(n, s/2)
  135. N = max(0, -nb)
  136. pN = p**N
  137. pn = p**-n # This is 1/h
  138. A = a*pN
  139. B = ba*pN.quo(bd) + Poly(n, DE.t)*a*derivation(p, DE).quo(p)*pN
  140. G = [(Ga*pN*pn).cancel(Gd, include=True) for Ga, Gd in G]
  141. h = pn
  142. # (a*p**N, (b + n*a*Dp/p)*p**N, g1*p**(N - n), ..., gm*p**(N - n), p**-n)
  143. return (A, B, G, h)
  144. def prde_linear_constraints(a, b, G, DE):
  145. """
  146. Parametric Risch Differential Equation - Generate linear constraints on the constants.
  147. Explanation
  148. ===========
  149. Given a derivation D on k[t], a, b, in k[t] with gcd(a, b) == 1, and
  150. G = [g1, ..., gm] in k(t)^m, return Q = [q1, ..., qm] in k[t]^m and a
  151. matrix M with entries in k(t) such that for any solution c1, ..., cm in
  152. Const(k) and p in k[t] of a*Dp + b*p == Sum(ci*gi, (i, 1, m)),
  153. (c1, ..., cm) is a solution of Mx == 0, and p and the ci satisfy
  154. a*Dp + b*p == Sum(ci*qi, (i, 1, m)).
  155. Because M has entries in k(t), and because Matrix does not play well with
  156. Poly, M will be a Matrix of Basic expressions.
  157. """
  158. m = len(G)
  159. Gns, Gds = list(zip(*G))
  160. d = reduce(lambda i, j: i.lcm(j), Gds)
  161. d = Poly(d, field=True)
  162. Q = [(ga*(d).quo(gd)).div(d) for ga, gd in G]
  163. if not all(ri.is_zero for _, ri in Q):
  164. N = max(ri.degree(DE.t) for _, ri in Q)
  165. M = Matrix(N + 1, m, lambda i, j: Q[j][1].nth(i), DE.t)
  166. else:
  167. M = Matrix(0, m, [], DE.t) # No constraints, return the empty matrix.
  168. qs, _ = list(zip(*Q))
  169. return (qs, M)
  170. def poly_linear_constraints(p, d):
  171. """
  172. Given p = [p1, ..., pm] in k[t]^m and d in k[t], return
  173. q = [q1, ..., qm] in k[t]^m and a matrix M with entries in k such
  174. that Sum(ci*pi, (i, 1, m)), for c1, ..., cm in k, is divisible
  175. by d if and only if (c1, ..., cm) is a solution of Mx = 0, in
  176. which case the quotient is Sum(ci*qi, (i, 1, m)).
  177. """
  178. m = len(p)
  179. q, r = zip(*[pi.div(d) for pi in p])
  180. if not all(ri.is_zero for ri in r):
  181. n = max(ri.degree() for ri in r)
  182. M = Matrix(n + 1, m, lambda i, j: r[j].nth(i), d.gens)
  183. else:
  184. M = Matrix(0, m, [], d.gens) # No constraints.
  185. return q, M
  186. def constant_system(A, u, DE):
  187. """
  188. Generate a system for the constant solutions.
  189. Explanation
  190. ===========
  191. Given a differential field (K, D) with constant field C = Const(K), a Matrix
  192. A, and a vector (Matrix) u with coefficients in K, returns the tuple
  193. (B, v, s), where B is a Matrix with coefficients in C and v is a vector
  194. (Matrix) such that either v has coefficients in C, in which case s is True
  195. and the solutions in C of Ax == u are exactly all the solutions of Bx == v,
  196. or v has a non-constant coefficient, in which case s is False Ax == u has no
  197. constant solution.
  198. This algorithm is used both in solving parametric problems and in
  199. determining if an element a of K is a derivative of an element of K or the
  200. logarithmic derivative of a K-radical using the structure theorem approach.
  201. Because Poly does not play well with Matrix yet, this algorithm assumes that
  202. all matrix entries are Basic expressions.
  203. """
  204. if not A:
  205. return A, u
  206. Au = A.row_join(u)
  207. Au, _ = Au.rref()
  208. # Warning: This will NOT return correct results if cancel() cannot reduce
  209. # an identically zero expression to 0. The danger is that we might
  210. # incorrectly prove that an integral is nonelementary (such as
  211. # risch_integrate(exp((sin(x)**2 + cos(x)**2 - 1)*x**2), x).
  212. # But this is a limitation in computer algebra in general, and implicit
  213. # in the correctness of the Risch Algorithm is the computability of the
  214. # constant field (actually, this same correctness problem exists in any
  215. # algorithm that uses rref()).
  216. #
  217. # We therefore limit ourselves to constant fields that are computable
  218. # via the cancel() function, in order to prevent a speed bottleneck from
  219. # calling some more complex simplification function (rational function
  220. # coefficients will fall into this class). Furthermore, (I believe) this
  221. # problem will only crop up if the integral explicitly contains an
  222. # expression in the constant field that is identically zero, but cannot
  223. # be reduced to such by cancel(). Therefore, a careful user can avoid this
  224. # problem entirely by being careful with the sorts of expressions that
  225. # appear in his integrand in the variables other than the integration
  226. # variable (the structure theorems should be able to completely decide these
  227. # problems in the integration variable).
  228. A, u = Au[:, :-1], Au[:, -1]
  229. D = lambda x: derivation(x, DE, basic=True)
  230. for j, i in itertools.product(range(A.cols), range(A.rows)):
  231. if A[i, j].expr.has(*DE.T):
  232. # This assumes that const(F(t0, ..., tn) == const(K) == F
  233. Ri = A[i, :]
  234. # Rm+1; m = A.rows
  235. DAij = D(A[i, j])
  236. Rm1 = Ri.applyfunc(lambda x: D(x) / DAij)
  237. um1 = D(u[i]) / DAij
  238. Aj = A[:, j]
  239. A = A - Aj * Rm1
  240. u = u - Aj * um1
  241. A = A.col_join(Rm1)
  242. u = u.col_join(Matrix([um1], u.gens))
  243. return (A, u)
  244. def prde_spde(a, b, Q, n, DE):
  245. """
  246. Special Polynomial Differential Equation algorithm: Parametric Version.
  247. Explanation
  248. ===========
  249. Given a derivation D on k[t], an integer n, and a, b, q1, ..., qm in k[t]
  250. with deg(a) > 0 and gcd(a, b) == 1, return (A, B, Q, R, n1), with
  251. Qq = [q1, ..., qm] and R = [r1, ..., rm], such that for any solution
  252. c1, ..., cm in Const(k) and q in k[t] of degree at most n of
  253. a*Dq + b*q == Sum(ci*gi, (i, 1, m)), p = (q - Sum(ci*ri, (i, 1, m)))/a has
  254. degree at most n1 and satisfies A*Dp + B*p == Sum(ci*qi, (i, 1, m))
  255. """
  256. R, Z = list(zip(*[gcdex_diophantine(b, a, qi) for qi in Q]))
  257. A = a
  258. B = b + derivation(a, DE)
  259. Qq = [zi - derivation(ri, DE) for ri, zi in zip(R, Z)]
  260. R = list(R)
  261. n1 = n - a.degree(DE.t)
  262. return (A, B, Qq, R, n1)
  263. def prde_no_cancel_b_large(b, Q, n, DE):
  264. """
  265. Parametric Poly Risch Differential Equation - No cancellation: deg(b) large enough.
  266. Explanation
  267. ===========
  268. Given a derivation D on k[t], n in ZZ, and b, q1, ..., qm in k[t] with
  269. b != 0 and either D == d/dt or deg(b) > max(0, deg(D) - 1), returns
  270. h1, ..., hr in k[t] and a matrix A with coefficients in Const(k) such that
  271. if c1, ..., cm in Const(k) and q in k[t] satisfy deg(q) <= n and
  272. Dq + b*q == Sum(ci*qi, (i, 1, m)), then q = Sum(dj*hj, (j, 1, r)), where
  273. d1, ..., dr in Const(k) and A*Matrix([[c1, ..., cm, d1, ..., dr]]).T == 0.
  274. """
  275. db = b.degree(DE.t)
  276. m = len(Q)
  277. H = [Poly(0, DE.t)]*m
  278. for N, i in itertools.product(range(n, -1, -1), range(m)): # [n, ..., 0]
  279. si = Q[i].nth(N + db)/b.LC()
  280. sitn = Poly(si*DE.t**N, DE.t)
  281. H[i] = H[i] + sitn
  282. Q[i] = Q[i] - derivation(sitn, DE) - b*sitn
  283. if all(qi.is_zero for qi in Q):
  284. dc = -1
  285. else:
  286. dc = max([qi.degree(DE.t) for qi in Q])
  287. M = Matrix(dc + 1, m, lambda i, j: Q[j].nth(i), DE.t)
  288. A, u = constant_system(M, zeros(dc + 1, 1, DE.t), DE)
  289. c = eye(m, DE.t)
  290. A = A.row_join(zeros(A.rows, m, DE.t)).col_join(c.row_join(-c))
  291. return (H, A)
  292. def prde_no_cancel_b_small(b, Q, n, DE):
  293. """
  294. Parametric Poly Risch Differential Equation - No cancellation: deg(b) small enough.
  295. Explanation
  296. ===========
  297. Given a derivation D on k[t], n in ZZ, and b, q1, ..., qm in k[t] with
  298. deg(b) < deg(D) - 1 and either D == d/dt or deg(D) >= 2, returns
  299. h1, ..., hr in k[t] and a matrix A with coefficients in Const(k) such that
  300. if c1, ..., cm in Const(k) and q in k[t] satisfy deg(q) <= n and
  301. Dq + b*q == Sum(ci*qi, (i, 1, m)) then q = Sum(dj*hj, (j, 1, r)) where
  302. d1, ..., dr in Const(k) and A*Matrix([[c1, ..., cm, d1, ..., dr]]).T == 0.
  303. """
  304. m = len(Q)
  305. H = [Poly(0, DE.t)]*m
  306. for N, i in itertools.product(range(n, 0, -1), range(m)): # [n, ..., 1]
  307. si = Q[i].nth(N + DE.d.degree(DE.t) - 1)/(N*DE.d.LC())
  308. sitn = Poly(si*DE.t**N, DE.t)
  309. H[i] = H[i] + sitn
  310. Q[i] = Q[i] - derivation(sitn, DE) - b*sitn
  311. if b.degree(DE.t) > 0:
  312. for i in range(m):
  313. si = Poly(Q[i].nth(b.degree(DE.t))/b.LC(), DE.t)
  314. H[i] = H[i] + si
  315. Q[i] = Q[i] - derivation(si, DE) - b*si
  316. if all(qi.is_zero for qi in Q):
  317. dc = -1
  318. else:
  319. dc = max([qi.degree(DE.t) for qi in Q])
  320. M = Matrix(dc + 1, m, lambda i, j: Q[j].nth(i), DE.t)
  321. A, u = constant_system(M, zeros(dc + 1, 1, DE.t), DE)
  322. c = eye(m, DE.t)
  323. A = A.row_join(zeros(A.rows, m, DE.t)).col_join(c.row_join(-c))
  324. return (H, A)
  325. # else: b is in k, deg(qi) < deg(Dt)
  326. t = DE.t
  327. if DE.case != 'base':
  328. with DecrementLevel(DE):
  329. t0 = DE.t # k = k0(t0)
  330. ba, bd = frac_in(b, t0, field=True)
  331. Q0 = [frac_in(qi.TC(), t0, field=True) for qi in Q]
  332. f, B = param_rischDE(ba, bd, Q0, DE)
  333. # f = [f1, ..., fr] in k^r and B is a matrix with
  334. # m + r columns and entries in Const(k) = Const(k0)
  335. # such that Dy0 + b*y0 = Sum(ci*qi, (i, 1, m)) has
  336. # a solution y0 in k with c1, ..., cm in Const(k)
  337. # if and only y0 = Sum(dj*fj, (j, 1, r)) where
  338. # d1, ..., dr ar in Const(k) and
  339. # B*Matrix([c1, ..., cm, d1, ..., dr]) == 0.
  340. # Transform fractions (fa, fd) in f into constant
  341. # polynomials fa/fd in k[t].
  342. # (Is there a better way?)
  343. f = [Poly(fa.as_expr()/fd.as_expr(), t, field=True)
  344. for fa, fd in f]
  345. B = Matrix.from_Matrix(B.to_Matrix(), t)
  346. else:
  347. # Base case. Dy == 0 for all y in k and b == 0.
  348. # Dy + b*y = Sum(ci*qi) is solvable if and only if
  349. # Sum(ci*qi) == 0 in which case the solutions are
  350. # y = d1*f1 for f1 = 1 and any d1 in Const(k) = k.
  351. f = [Poly(1, t, field=True)] # r = 1
  352. B = Matrix([[qi.TC() for qi in Q] + [S.Zero]], DE.t)
  353. # The condition for solvability is
  354. # B*Matrix([c1, ..., cm, d1]) == 0
  355. # There are no constraints on d1.
  356. # Coefficients of t^j (j > 0) in Sum(ci*qi) must be zero.
  357. d = max([qi.degree(DE.t) for qi in Q])
  358. if d > 0:
  359. M = Matrix(d, m, lambda i, j: Q[j].nth(i + 1), DE.t)
  360. A, _ = constant_system(M, zeros(d, 1, DE.t), DE)
  361. else:
  362. # No constraints on the hj.
  363. A = Matrix(0, m, [], DE.t)
  364. # Solutions of the original equation are
  365. # y = Sum(dj*fj, (j, 1, r) + Sum(ei*hi, (i, 1, m)),
  366. # where ei == ci (i = 1, ..., m), when
  367. # A*Matrix([c1, ..., cm]) == 0 and
  368. # B*Matrix([c1, ..., cm, d1, ..., dr]) == 0
  369. # Build combined constraint matrix with m + r + m columns.
  370. r = len(f)
  371. I = eye(m, DE.t)
  372. A = A.row_join(zeros(A.rows, r + m, DE.t))
  373. B = B.row_join(zeros(B.rows, m, DE.t))
  374. C = I.row_join(zeros(m, r, DE.t)).row_join(-I)
  375. return f + H, A.col_join(B).col_join(C)
  376. def prde_cancel_liouvillian(b, Q, n, DE):
  377. """
  378. Pg, 237.
  379. """
  380. H = []
  381. # Why use DecrementLevel? Below line answers that:
  382. # Assuming that we can solve such problems over 'k' (not k[t])
  383. if DE.case == 'primitive':
  384. with DecrementLevel(DE):
  385. ba, bd = frac_in(b, DE.t, field=True)
  386. for i in range(n, -1, -1):
  387. if DE.case == 'exp': # this re-checking can be avoided
  388. with DecrementLevel(DE):
  389. ba, bd = frac_in(b + (i*(derivation(DE.t, DE)/DE.t)).as_poly(b.gens),
  390. DE.t, field=True)
  391. with DecrementLevel(DE):
  392. Qy = [frac_in(q.nth(i), DE.t, field=True) for q in Q]
  393. fi, Ai = param_rischDE(ba, bd, Qy, DE)
  394. fi = [Poly(fa.as_expr()/fd.as_expr(), DE.t, field=True)
  395. for fa, fd in fi]
  396. Ai = Ai.set_gens(DE.t)
  397. ri = len(fi)
  398. if i == n:
  399. M = Ai
  400. else:
  401. M = Ai.col_join(M.row_join(zeros(M.rows, ri, DE.t)))
  402. Fi, hi = [None]*ri, [None]*ri
  403. # from eq. on top of p.238 (unnumbered)
  404. for j in range(ri):
  405. hji = fi[j] * (DE.t**i).as_poly(fi[j].gens)
  406. hi[j] = hji
  407. # building up Sum(djn*(D(fjn*t^n) - b*fjnt^n))
  408. Fi[j] = -(derivation(hji, DE) - b*hji)
  409. H += hi
  410. # in the next loop instead of Q it has
  411. # to be Q + Fi taking its place
  412. Q = Q + Fi
  413. return (H, M)
  414. def param_poly_rischDE(a, b, q, n, DE):
  415. """Polynomial solutions of a parametric Risch differential equation.
  416. Explanation
  417. ===========
  418. Given a derivation D in k[t], a, b in k[t] relatively prime, and q
  419. = [q1, ..., qm] in k[t]^m, return h = [h1, ..., hr] in k[t]^r and
  420. a matrix A with m + r columns and entries in Const(k) such that
  421. a*Dp + b*p = Sum(ci*qi, (i, 1, m)) has a solution p of degree <= n
  422. in k[t] with c1, ..., cm in Const(k) if and only if p = Sum(dj*hj,
  423. (j, 1, r)) where d1, ..., dr are in Const(k) and (c1, ..., cm,
  424. d1, ..., dr) is a solution of Ax == 0.
  425. """
  426. m = len(q)
  427. if n < 0:
  428. # Only the trivial zero solution is possible.
  429. # Find relations between the qi.
  430. if all(qi.is_zero for qi in q):
  431. return [], zeros(1, m, DE.t) # No constraints.
  432. N = max([qi.degree(DE.t) for qi in q])
  433. M = Matrix(N + 1, m, lambda i, j: q[j].nth(i), DE.t)
  434. A, _ = constant_system(M, zeros(M.rows, 1, DE.t), DE)
  435. return [], A
  436. if a.is_ground:
  437. # Normalization: a = 1.
  438. a = a.LC()
  439. b, q = b.quo_ground(a), [qi.quo_ground(a) for qi in q]
  440. if not b.is_zero and (DE.case == 'base' or
  441. b.degree() > max(0, DE.d.degree() - 1)):
  442. return prde_no_cancel_b_large(b, q, n, DE)
  443. elif ((b.is_zero or b.degree() < DE.d.degree() - 1)
  444. and (DE.case == 'base' or DE.d.degree() >= 2)):
  445. return prde_no_cancel_b_small(b, q, n, DE)
  446. elif (DE.d.degree() >= 2 and
  447. b.degree() == DE.d.degree() - 1 and
  448. n > -b.as_poly().LC()/DE.d.as_poly().LC()):
  449. raise NotImplementedError("prde_no_cancel_b_equal() is "
  450. "not yet implemented.")
  451. else:
  452. # Liouvillian cases
  453. if DE.case in ('primitive', 'exp'):
  454. return prde_cancel_liouvillian(b, q, n, DE)
  455. else:
  456. raise NotImplementedError("non-linear and hypertangent "
  457. "cases have not yet been implemented")
  458. # else: deg(a) > 0
  459. # Iterate SPDE as long as possible cumulating coefficient
  460. # and terms for the recovery of original solutions.
  461. alpha, beta = a.one, [a.zero]*m
  462. while n >= 0: # and a, b relatively prime
  463. a, b, q, r, n = prde_spde(a, b, q, n, DE)
  464. beta = [betai + alpha*ri for betai, ri in zip(beta, r)]
  465. alpha *= a
  466. # Solutions p of a*Dp + b*p = Sum(ci*qi) correspond to
  467. # solutions alpha*p + Sum(ci*betai) of the initial equation.
  468. d = a.gcd(b)
  469. if not d.is_ground:
  470. break
  471. # a*Dp + b*p = Sum(ci*qi) may have a polynomial solution
  472. # only if the sum is divisible by d.
  473. qq, M = poly_linear_constraints(q, d)
  474. # qq = [qq1, ..., qqm] where qqi = qi.quo(d).
  475. # M is a matrix with m columns an entries in k.
  476. # Sum(fi*qi, (i, 1, m)), where f1, ..., fm are elements of k, is
  477. # divisible by d if and only if M*Matrix([f1, ..., fm]) == 0,
  478. # in which case the quotient is Sum(fi*qqi).
  479. A, _ = constant_system(M, zeros(M.rows, 1, DE.t), DE)
  480. # A is a matrix with m columns and entries in Const(k).
  481. # Sum(ci*qqi) is Sum(ci*qi).quo(d), and the remainder is zero
  482. # for c1, ..., cm in Const(k) if and only if
  483. # A*Matrix([c1, ...,cm]) == 0.
  484. V = A.nullspace()
  485. # V = [v1, ..., vu] where each vj is a column matrix with
  486. # entries aj1, ..., ajm in Const(k).
  487. # Sum(aji*qi) is divisible by d with exact quotient Sum(aji*qqi).
  488. # Sum(ci*qi) is divisible by d if and only if ci = Sum(dj*aji)
  489. # (i = 1, ..., m) for some d1, ..., du in Const(k).
  490. # In that case, solutions of
  491. # a*Dp + b*p = Sum(ci*qi) = Sum(dj*Sum(aji*qi))
  492. # are the same as those of
  493. # (a/d)*Dp + (b/d)*p = Sum(dj*rj)
  494. # where rj = Sum(aji*qqi).
  495. if not V: # No non-trivial solution.
  496. return [], eye(m, DE.t) # Could return A, but this has
  497. # the minimum number of rows.
  498. Mqq = Matrix([qq]) # A single row.
  499. r = [(Mqq*vj)[0] for vj in V] # [r1, ..., ru]
  500. # Solutions of (a/d)*Dp + (b/d)*p = Sum(dj*rj) correspond to
  501. # solutions alpha*p + Sum(Sum(dj*aji)*betai) of the initial
  502. # equation. These are equal to alpha*p + Sum(dj*fj) where
  503. # fj = Sum(aji*betai).
  504. Mbeta = Matrix([beta])
  505. f = [(Mbeta*vj)[0] for vj in V] # [f1, ..., fu]
  506. #
  507. # Solve the reduced equation recursively.
  508. #
  509. g, B = param_poly_rischDE(a.quo(d), b.quo(d), r, n, DE)
  510. # g = [g1, ..., gv] in k[t]^v and and B is a matrix with u + v
  511. # columns and entries in Const(k) such that
  512. # (a/d)*Dp + (b/d)*p = Sum(dj*rj) has a solution p of degree <= n
  513. # in k[t] if and only if p = Sum(ek*gk) where e1, ..., ev are in
  514. # Const(k) and B*Matrix([d1, ..., du, e1, ..., ev]) == 0.
  515. # The solutions of the original equation are then
  516. # Sum(dj*fj, (j, 1, u)) + alpha*Sum(ek*gk, (k, 1, v)).
  517. # Collect solution components.
  518. h = f + [alpha*gk for gk in g]
  519. # Build combined relation matrix.
  520. A = -eye(m, DE.t)
  521. for vj in V:
  522. A = A.row_join(vj)
  523. A = A.row_join(zeros(m, len(g), DE.t))
  524. A = A.col_join(zeros(B.rows, m, DE.t).row_join(B))
  525. return h, A
  526. def param_rischDE(fa, fd, G, DE):
  527. """
  528. Solve a Parametric Risch Differential Equation: Dy + f*y == Sum(ci*Gi, (i, 1, m)).
  529. Explanation
  530. ===========
  531. Given a derivation D in k(t), f in k(t), and G
  532. = [G1, ..., Gm] in k(t)^m, return h = [h1, ..., hr] in k(t)^r and
  533. a matrix A with m + r columns and entries in Const(k) such that
  534. Dy + f*y = Sum(ci*Gi, (i, 1, m)) has a solution y
  535. in k(t) with c1, ..., cm in Const(k) if and only if y = Sum(dj*hj,
  536. (j, 1, r)) where d1, ..., dr are in Const(k) and (c1, ..., cm,
  537. d1, ..., dr) is a solution of Ax == 0.
  538. Elements of k(t) are tuples (a, d) with a and d in k[t].
  539. """
  540. m = len(G)
  541. q, (fa, fd) = weak_normalizer(fa, fd, DE)
  542. # Solutions of the weakly normalized equation Dz + f*z = q*Sum(ci*Gi)
  543. # correspond to solutions y = z/q of the original equation.
  544. gamma = q
  545. G = [(q*ga).cancel(gd, include=True) for ga, gd in G]
  546. a, (ba, bd), G, hn = prde_normal_denom(fa, fd, G, DE)
  547. # Solutions q in k<t> of a*Dq + b*q = Sum(ci*Gi) correspond
  548. # to solutions z = q/hn of the weakly normalized equation.
  549. gamma *= hn
  550. A, B, G, hs = prde_special_denom(a, ba, bd, G, DE)
  551. # Solutions p in k[t] of A*Dp + B*p = Sum(ci*Gi) correspond
  552. # to solutions q = p/hs of the previous equation.
  553. gamma *= hs
  554. g = A.gcd(B)
  555. a, b, g = A.quo(g), B.quo(g), [gia.cancel(gid*g, include=True) for
  556. gia, gid in G]
  557. # a*Dp + b*p = Sum(ci*gi) may have a polynomial solution
  558. # only if the sum is in k[t].
  559. q, M = prde_linear_constraints(a, b, g, DE)
  560. # q = [q1, ..., qm] where qi in k[t] is the polynomial component
  561. # of the partial fraction expansion of gi.
  562. # M is a matrix with m columns and entries in k.
  563. # Sum(fi*gi, (i, 1, m)), where f1, ..., fm are elements of k,
  564. # is a polynomial if and only if M*Matrix([f1, ..., fm]) == 0,
  565. # in which case the sum is equal to Sum(fi*qi).
  566. M, _ = constant_system(M, zeros(M.rows, 1, DE.t), DE)
  567. # M is a matrix with m columns and entries in Const(k).
  568. # Sum(ci*gi) is in k[t] for c1, ..., cm in Const(k)
  569. # if and only if M*Matrix([c1, ..., cm]) == 0,
  570. # in which case the sum is Sum(ci*qi).
  571. ## Reduce number of constants at this point
  572. V = M.nullspace()
  573. # V = [v1, ..., vu] where each vj is a column matrix with
  574. # entries aj1, ..., ajm in Const(k).
  575. # Sum(aji*gi) is in k[t] and equal to Sum(aji*qi) (j = 1, ..., u).
  576. # Sum(ci*gi) is in k[t] if and only is ci = Sum(dj*aji)
  577. # (i = 1, ..., m) for some d1, ..., du in Const(k).
  578. # In that case,
  579. # Sum(ci*gi) = Sum(ci*qi) = Sum(dj*Sum(aji*qi)) = Sum(dj*rj)
  580. # where rj = Sum(aji*qi) (j = 1, ..., u) in k[t].
  581. if not V: # No non-trivial solution
  582. return [], eye(m, DE.t)
  583. Mq = Matrix([q]) # A single row.
  584. r = [(Mq*vj)[0] for vj in V] # [r1, ..., ru]
  585. # Solutions of a*Dp + b*p = Sum(dj*rj) correspond to solutions
  586. # y = p/gamma of the initial equation with ci = Sum(dj*aji).
  587. try:
  588. # We try n=5. At least for prde_spde, it will always
  589. # terminate no matter what n is.
  590. n = bound_degree(a, b, r, DE, parametric=True)
  591. except NotImplementedError:
  592. # A temporary bound is set. Eventually, it will be removed.
  593. # the currently added test case takes large time
  594. # even with n=5, and much longer with large n's.
  595. n = 5
  596. h, B = param_poly_rischDE(a, b, r, n, DE)
  597. # h = [h1, ..., hv] in k[t]^v and and B is a matrix with u + v
  598. # columns and entries in Const(k) such that
  599. # a*Dp + b*p = Sum(dj*rj) has a solution p of degree <= n
  600. # in k[t] if and only if p = Sum(ek*hk) where e1, ..., ev are in
  601. # Const(k) and B*Matrix([d1, ..., du, e1, ..., ev]) == 0.
  602. # The solutions of the original equation for ci = Sum(dj*aji)
  603. # (i = 1, ..., m) are then y = Sum(ek*hk, (k, 1, v))/gamma.
  604. ## Build combined relation matrix with m + u + v columns.
  605. A = -eye(m, DE.t)
  606. for vj in V:
  607. A = A.row_join(vj)
  608. A = A.row_join(zeros(m, len(h), DE.t))
  609. A = A.col_join(zeros(B.rows, m, DE.t).row_join(B))
  610. ## Eliminate d1, ..., du.
  611. W = A.nullspace()
  612. # W = [w1, ..., wt] where each wl is a column matrix with
  613. # entries blk (k = 1, ..., m + u + v) in Const(k).
  614. # The vectors (bl1, ..., blm) generate the space of those
  615. # constant families (c1, ..., cm) for which a solution of
  616. # the equation Dy + f*y == Sum(ci*Gi) exists. They generate
  617. # the space and form a basis except possibly when Dy + f*y == 0
  618. # is solvable in k(t}. The corresponding solutions are
  619. # y = Sum(blk'*hk, (k, 1, v))/gamma, where k' = k + m + u.
  620. v = len(h)
  621. shape = (len(W), m+v)
  622. elements = [wl[:m] + wl[-v:] for wl in W] # excise dj's.
  623. items = [e for row in elements for e in row]
  624. # Need to set the shape in case W is empty
  625. M = Matrix(*shape, items, DE.t)
  626. N = M.nullspace()
  627. # N = [n1, ..., ns] where the ni in Const(k)^(m + v) are column
  628. # vectors generating the space of linear relations between
  629. # c1, ..., cm, e1, ..., ev.
  630. C = Matrix([ni[:] for ni in N], DE.t) # rows n1, ..., ns.
  631. return [hk.cancel(gamma, include=True) for hk in h], C
  632. def limited_integrate_reduce(fa, fd, G, DE):
  633. """
  634. Simpler version of step 1 & 2 for the limited integration problem.
  635. Explanation
  636. ===========
  637. Given a derivation D on k(t) and f, g1, ..., gn in k(t), return
  638. (a, b, h, N, g, V) such that a, b, h in k[t], N is a non-negative integer,
  639. g in k(t), V == [v1, ..., vm] in k(t)^m, and for any solution v in k(t),
  640. c1, ..., cm in C of f == Dv + Sum(ci*wi, (i, 1, m)), p = v*h is in k<t>, and
  641. p and the ci satisfy a*Dp + b*p == g + Sum(ci*vi, (i, 1, m)). Furthermore,
  642. if S1irr == Sirr, then p is in k[t], and if t is nonlinear or Liouvillian
  643. over k, then deg(p) <= N.
  644. So that the special part is always computed, this function calls the more
  645. general prde_special_denom() automatically if it cannot determine that
  646. S1irr == Sirr. Furthermore, it will automatically call bound_degree() when
  647. t is linear and non-Liouvillian, which for the transcendental case, implies
  648. that Dt == a*t + b with for some a, b in k*.
  649. """
  650. dn, ds = splitfactor(fd, DE)
  651. E = [splitfactor(gd, DE) for _, gd in G]
  652. En, Es = list(zip(*E))
  653. c = reduce(lambda i, j: i.lcm(j), (dn,) + En) # lcm(dn, en1, ..., enm)
  654. hn = c.gcd(c.diff(DE.t))
  655. a = hn
  656. b = -derivation(hn, DE)
  657. N = 0
  658. # These are the cases where we know that S1irr = Sirr, but there could be
  659. # others, and this algorithm will need to be extended to handle them.
  660. if DE.case in ('base', 'primitive', 'exp', 'tan'):
  661. hs = reduce(lambda i, j: i.lcm(j), (ds,) + Es) # lcm(ds, es1, ..., esm)
  662. a = hn*hs
  663. b -= (hn*derivation(hs, DE)).quo(hs)
  664. mu = min(order_at_oo(fa, fd, DE.t), min([order_at_oo(ga, gd, DE.t) for
  665. ga, gd in G]))
  666. # So far, all the above are also nonlinear or Liouvillian, but if this
  667. # changes, then this will need to be updated to call bound_degree()
  668. # as per the docstring of this function (DE.case == 'other_linear').
  669. N = hn.degree(DE.t) + hs.degree(DE.t) + max(0, 1 - DE.d.degree(DE.t) - mu)
  670. else:
  671. # TODO: implement this
  672. raise NotImplementedError
  673. V = [(-a*hn*ga).cancel(gd, include=True) for ga, gd in G]
  674. return (a, b, a, N, (a*hn*fa).cancel(fd, include=True), V)
  675. def limited_integrate(fa, fd, G, DE):
  676. """
  677. Solves the limited integration problem: f = Dv + Sum(ci*wi, (i, 1, n))
  678. """
  679. fa, fd = fa*Poly(1/fd.LC(), DE.t), fd.monic()
  680. # interpreting limited integration problem as a
  681. # parametric Risch DE problem
  682. Fa = Poly(0, DE.t)
  683. Fd = Poly(1, DE.t)
  684. G = [(fa, fd)] + G
  685. h, A = param_rischDE(Fa, Fd, G, DE)
  686. V = A.nullspace()
  687. V = [v for v in V if v[0] != 0]
  688. if not V:
  689. return None
  690. else:
  691. # we can take any vector from V, we take V[0]
  692. c0 = V[0][0]
  693. # v = [-1, c1, ..., cm, d1, ..., dr]
  694. v = V[0]/(-c0)
  695. r = len(h)
  696. m = len(v) - r - 1
  697. C = list(v[1: m + 1])
  698. y = -sum([v[m + 1 + i]*h[i][0].as_expr()/h[i][1].as_expr() \
  699. for i in range(r)])
  700. y_num, y_den = y.as_numer_denom()
  701. Ya, Yd = Poly(y_num, DE.t), Poly(y_den, DE.t)
  702. Y = Ya*Poly(1/Yd.LC(), DE.t), Yd.monic()
  703. return Y, C
  704. def parametric_log_deriv_heu(fa, fd, wa, wd, DE, c1=None):
  705. """
  706. Parametric logarithmic derivative heuristic.
  707. Explanation
  708. ===========
  709. Given a derivation D on k[t], f in k(t), and a hyperexponential monomial
  710. theta over k(t), raises either NotImplementedError, in which case the
  711. heuristic failed, or returns None, in which case it has proven that no
  712. solution exists, or returns a solution (n, m, v) of the equation
  713. n*f == Dv/v + m*Dtheta/theta, with v in k(t)* and n, m in ZZ with n != 0.
  714. If this heuristic fails, the structure theorem approach will need to be
  715. used.
  716. The argument w == Dtheta/theta
  717. """
  718. # TODO: finish writing this and write tests
  719. c1 = c1 or Dummy('c1')
  720. p, a = fa.div(fd)
  721. q, b = wa.div(wd)
  722. B = max(0, derivation(DE.t, DE).degree(DE.t) - 1)
  723. C = max(p.degree(DE.t), q.degree(DE.t))
  724. if q.degree(DE.t) > B:
  725. eqs = [p.nth(i) - c1*q.nth(i) for i in range(B + 1, C + 1)]
  726. s = solve(eqs, c1)
  727. if not s or not s[c1].is_Rational:
  728. # deg(q) > B, no solution for c.
  729. return None
  730. M, N = s[c1].as_numer_denom()
  731. M_poly = M.as_poly(q.gens)
  732. N_poly = N.as_poly(q.gens)
  733. nfmwa = N_poly*fa*wd - M_poly*wa*fd
  734. nfmwd = fd*wd
  735. Qv = is_log_deriv_k_t_radical_in_field(nfmwa, nfmwd, DE, 'auto')
  736. if Qv is None:
  737. # (N*f - M*w) is not the logarithmic derivative of a k(t)-radical.
  738. return None
  739. Q, v = Qv
  740. if Q.is_zero or v.is_zero:
  741. return None
  742. return (Q*N, Q*M, v)
  743. if p.degree(DE.t) > B:
  744. return None
  745. c = lcm(fd.as_poly(DE.t).LC(), wd.as_poly(DE.t).LC())
  746. l = fd.monic().lcm(wd.monic())*Poly(c, DE.t)
  747. ln, ls = splitfactor(l, DE)
  748. z = ls*ln.gcd(ln.diff(DE.t))
  749. if not z.has(DE.t):
  750. # TODO: We treat this as 'no solution', until the structure
  751. # theorem version of parametric_log_deriv is implemented.
  752. return None
  753. u1, r1 = (fa*l.quo(fd)).div(z) # (l*f).div(z)
  754. u2, r2 = (wa*l.quo(wd)).div(z) # (l*w).div(z)
  755. eqs = [r1.nth(i) - c1*r2.nth(i) for i in range(z.degree(DE.t))]
  756. s = solve(eqs, c1)
  757. if not s or not s[c1].is_Rational:
  758. # deg(q) <= B, no solution for c.
  759. return None
  760. M, N = s[c1].as_numer_denom()
  761. nfmwa = N.as_poly(DE.t)*fa*wd - M.as_poly(DE.t)*wa*fd
  762. nfmwd = fd*wd
  763. Qv = is_log_deriv_k_t_radical_in_field(nfmwa, nfmwd, DE)
  764. if Qv is None:
  765. # (N*f - M*w) is not the logarithmic derivative of a k(t)-radical.
  766. return None
  767. Q, v = Qv
  768. if Q.is_zero or v.is_zero:
  769. return None
  770. return (Q*N, Q*M, v)
  771. def parametric_log_deriv(fa, fd, wa, wd, DE):
  772. # TODO: Write the full algorithm using the structure theorems.
  773. # try:
  774. A = parametric_log_deriv_heu(fa, fd, wa, wd, DE)
  775. # except NotImplementedError:
  776. # Heuristic failed, we have to use the full method.
  777. # TODO: This could be implemented more efficiently.
  778. # It isn't too worrisome, because the heuristic handles most difficult
  779. # cases.
  780. return A
  781. def is_deriv_k(fa, fd, DE):
  782. r"""
  783. Checks if Df/f is the derivative of an element of k(t).
  784. Explanation
  785. ===========
  786. a in k(t) is the derivative of an element of k(t) if there exists b in k(t)
  787. such that a = Db. Either returns (ans, u), such that Df/f == Du, or None,
  788. which means that Df/f is not the derivative of an element of k(t). ans is
  789. a list of tuples such that Add(*[i*j for i, j in ans]) == u. This is useful
  790. for seeing exactly which elements of k(t) produce u.
  791. This function uses the structure theorem approach, which says that for any
  792. f in K, Df/f is the derivative of a element of K if and only if there are ri
  793. in QQ such that::
  794. --- --- Dt
  795. \ r * Dt + \ r * i Df
  796. / i i / i --- = --.
  797. --- --- t f
  798. i in L i in E i
  799. K/C(x) K/C(x)
  800. Where C = Const(K), L_K/C(x) = { i in {1, ..., n} such that t_i is
  801. transcendental over C(x)(t_1, ..., t_i-1) and Dt_i = Da_i/a_i, for some a_i
  802. in C(x)(t_1, ..., t_i-1)* } (i.e., the set of all indices of logarithmic
  803. monomials of K over C(x)), and E_K/C(x) = { i in {1, ..., n} such that t_i
  804. is transcendental over C(x)(t_1, ..., t_i-1) and Dt_i/t_i = Da_i, for some
  805. a_i in C(x)(t_1, ..., t_i-1) } (i.e., the set of all indices of
  806. hyperexponential monomials of K over C(x)). If K is an elementary extension
  807. over C(x), then the cardinality of L_K/C(x) U E_K/C(x) is exactly the
  808. transcendence degree of K over C(x). Furthermore, because Const_D(K) ==
  809. Const_D(C(x)) == C, deg(Dt_i) == 1 when t_i is in E_K/C(x) and
  810. deg(Dt_i) == 0 when t_i is in L_K/C(x), implying in particular that E_K/C(x)
  811. and L_K/C(x) are disjoint.
  812. The sets L_K/C(x) and E_K/C(x) must, by their nature, be computed
  813. recursively using this same function. Therefore, it is required to pass
  814. them as indices to D (or T). E_args are the arguments of the
  815. hyperexponentials indexed by E_K (i.e., if i is in E_K, then T[i] ==
  816. exp(E_args[i])). This is needed to compute the final answer u such that
  817. Df/f == Du.
  818. log(f) will be the same as u up to a additive constant. This is because
  819. they will both behave the same as monomials. For example, both log(x) and
  820. log(2*x) == log(x) + log(2) satisfy Dt == 1/x, because log(2) is constant.
  821. Therefore, the term const is returned. const is such that
  822. log(const) + f == u. This is calculated by dividing the arguments of one
  823. logarithm from the other. Therefore, it is necessary to pass the arguments
  824. of the logarithmic terms in L_args.
  825. To handle the case where we are given Df/f, not f, use is_deriv_k_in_field().
  826. See also
  827. ========
  828. is_log_deriv_k_t_radical_in_field, is_log_deriv_k_t_radical
  829. """
  830. # Compute Df/f
  831. dfa, dfd = (fd*derivation(fa, DE) - fa*derivation(fd, DE)), fd*fa
  832. dfa, dfd = dfa.cancel(dfd, include=True)
  833. # Our assumption here is that each monomial is recursively transcendental
  834. if len(DE.exts) != len(DE.D):
  835. if [i for i in DE.cases if i == 'tan'] or \
  836. ({i for i in DE.cases if i == 'primitive'} -
  837. set(DE.indices('log'))):
  838. raise NotImplementedError("Real version of the structure "
  839. "theorems with hypertangent support is not yet implemented.")
  840. # TODO: What should really be done in this case?
  841. raise NotImplementedError("Nonelementary extensions not supported "
  842. "in the structure theorems.")
  843. E_part = [DE.D[i].quo(Poly(DE.T[i], DE.T[i])).as_expr() for i in DE.indices('exp')]
  844. L_part = [DE.D[i].as_expr() for i in DE.indices('log')]
  845. # The expression dfa/dfd might not be polynomial in any of its symbols so we
  846. # use a Dummy as the generator for PolyMatrix.
  847. dum = Dummy()
  848. lhs = Matrix([E_part + L_part], dum)
  849. rhs = Matrix([dfa.as_expr()/dfd.as_expr()], dum)
  850. A, u = constant_system(lhs, rhs, DE)
  851. u = u.to_Matrix() # Poly to Expr
  852. if not A or not all(derivation(i, DE, basic=True).is_zero for i in u):
  853. # If the elements of u are not all constant
  854. # Note: See comment in constant_system
  855. # Also note: derivation(basic=True) calls cancel()
  856. return None
  857. else:
  858. if not all(i.is_Rational for i in u):
  859. raise NotImplementedError("Cannot work with non-rational "
  860. "coefficients in this case.")
  861. else:
  862. terms = ([DE.extargs[i] for i in DE.indices('exp')] +
  863. [DE.T[i] for i in DE.indices('log')])
  864. ans = list(zip(terms, u))
  865. result = Add(*[Mul(i, j) for i, j in ans])
  866. argterms = ([DE.T[i] for i in DE.indices('exp')] +
  867. [DE.extargs[i] for i in DE.indices('log')])
  868. l = []
  869. ld = []
  870. for i, j in zip(argterms, u):
  871. # We need to get around things like sqrt(x**2) != x
  872. # and also sqrt(x**2 + 2*x + 1) != x + 1
  873. # Issue 10798: i need not be a polynomial
  874. i, d = i.as_numer_denom()
  875. icoeff, iterms = sqf_list(i)
  876. l.append(Mul(*([Pow(icoeff, j)] + [Pow(b, e*j) for b, e in iterms])))
  877. dcoeff, dterms = sqf_list(d)
  878. ld.append(Mul(*([Pow(dcoeff, j)] + [Pow(b, e*j) for b, e in dterms])))
  879. const = cancel(fa.as_expr()/fd.as_expr()/Mul(*l)*Mul(*ld))
  880. return (ans, result, const)
  881. def is_log_deriv_k_t_radical(fa, fd, DE, Df=True):
  882. r"""
  883. Checks if Df is the logarithmic derivative of a k(t)-radical.
  884. Explanation
  885. ===========
  886. b in k(t) can be written as the logarithmic derivative of a k(t) radical if
  887. there exist n in ZZ and u in k(t) with n, u != 0 such that n*b == Du/u.
  888. Either returns (ans, u, n, const) or None, which means that Df cannot be
  889. written as the logarithmic derivative of a k(t)-radical. ans is a list of
  890. tuples such that Mul(*[i**j for i, j in ans]) == u. This is useful for
  891. seeing exactly what elements of k(t) produce u.
  892. This function uses the structure theorem approach, which says that for any
  893. f in K, Df is the logarithmic derivative of a K-radical if and only if there
  894. are ri in QQ such that::
  895. --- --- Dt
  896. \ r * Dt + \ r * i
  897. / i i / i --- = Df.
  898. --- --- t
  899. i in L i in E i
  900. K/C(x) K/C(x)
  901. Where C = Const(K), L_K/C(x) = { i in {1, ..., n} such that t_i is
  902. transcendental over C(x)(t_1, ..., t_i-1) and Dt_i = Da_i/a_i, for some a_i
  903. in C(x)(t_1, ..., t_i-1)* } (i.e., the set of all indices of logarithmic
  904. monomials of K over C(x)), and E_K/C(x) = { i in {1, ..., n} such that t_i
  905. is transcendental over C(x)(t_1, ..., t_i-1) and Dt_i/t_i = Da_i, for some
  906. a_i in C(x)(t_1, ..., t_i-1) } (i.e., the set of all indices of
  907. hyperexponential monomials of K over C(x)). If K is an elementary extension
  908. over C(x), then the cardinality of L_K/C(x) U E_K/C(x) is exactly the
  909. transcendence degree of K over C(x). Furthermore, because Const_D(K) ==
  910. Const_D(C(x)) == C, deg(Dt_i) == 1 when t_i is in E_K/C(x) and
  911. deg(Dt_i) == 0 when t_i is in L_K/C(x), implying in particular that E_K/C(x)
  912. and L_K/C(x) are disjoint.
  913. The sets L_K/C(x) and E_K/C(x) must, by their nature, be computed
  914. recursively using this same function. Therefore, it is required to pass
  915. them as indices to D (or T). L_args are the arguments of the logarithms
  916. indexed by L_K (i.e., if i is in L_K, then T[i] == log(L_args[i])). This is
  917. needed to compute the final answer u such that n*f == Du/u.
  918. exp(f) will be the same as u up to a multiplicative constant. This is
  919. because they will both behave the same as monomials. For example, both
  920. exp(x) and exp(x + 1) == E*exp(x) satisfy Dt == t. Therefore, the term const
  921. is returned. const is such that exp(const)*f == u. This is calculated by
  922. subtracting the arguments of one exponential from the other. Therefore, it
  923. is necessary to pass the arguments of the exponential terms in E_args.
  924. To handle the case where we are given Df, not f, use
  925. is_log_deriv_k_t_radical_in_field().
  926. See also
  927. ========
  928. is_log_deriv_k_t_radical_in_field, is_deriv_k
  929. """
  930. if Df:
  931. dfa, dfd = (fd*derivation(fa, DE) - fa*derivation(fd, DE)).cancel(fd**2,
  932. include=True)
  933. else:
  934. dfa, dfd = fa, fd
  935. # Our assumption here is that each monomial is recursively transcendental
  936. if len(DE.exts) != len(DE.D):
  937. if [i for i in DE.cases if i == 'tan'] or \
  938. ({i for i in DE.cases if i == 'primitive'} -
  939. set(DE.indices('log'))):
  940. raise NotImplementedError("Real version of the structure "
  941. "theorems with hypertangent support is not yet implemented.")
  942. # TODO: What should really be done in this case?
  943. raise NotImplementedError("Nonelementary extensions not supported "
  944. "in the structure theorems.")
  945. E_part = [DE.D[i].quo(Poly(DE.T[i], DE.T[i])).as_expr() for i in DE.indices('exp')]
  946. L_part = [DE.D[i].as_expr() for i in DE.indices('log')]
  947. # The expression dfa/dfd might not be polynomial in any of its symbols so we
  948. # use a Dummy as the generator for PolyMatrix.
  949. dum = Dummy()
  950. lhs = Matrix([E_part + L_part], dum)
  951. rhs = Matrix([dfa.as_expr()/dfd.as_expr()], dum)
  952. A, u = constant_system(lhs, rhs, DE)
  953. u = u.to_Matrix() # Poly to Expr
  954. if not A or not all(derivation(i, DE, basic=True).is_zero for i in u):
  955. # If the elements of u are not all constant
  956. # Note: See comment in constant_system
  957. # Also note: derivation(basic=True) calls cancel()
  958. return None
  959. else:
  960. if not all(i.is_Rational for i in u):
  961. # TODO: But maybe we can tell if they're not rational, like
  962. # log(2)/log(3). Also, there should be an option to continue
  963. # anyway, even if the result might potentially be wrong.
  964. raise NotImplementedError("Cannot work with non-rational "
  965. "coefficients in this case.")
  966. else:
  967. n = reduce(ilcm, [i.as_numer_denom()[1] for i in u])
  968. u *= n
  969. terms = ([DE.T[i] for i in DE.indices('exp')] +
  970. [DE.extargs[i] for i in DE.indices('log')])
  971. ans = list(zip(terms, u))
  972. result = Mul(*[Pow(i, j) for i, j in ans])
  973. # exp(f) will be the same as result up to a multiplicative
  974. # constant. We now find the log of that constant.
  975. argterms = ([DE.extargs[i] for i in DE.indices('exp')] +
  976. [DE.T[i] for i in DE.indices('log')])
  977. const = cancel(fa.as_expr()/fd.as_expr() -
  978. Add(*[Mul(i, j/n) for i, j in zip(argterms, u)]))
  979. return (ans, result, n, const)
  980. def is_log_deriv_k_t_radical_in_field(fa, fd, DE, case='auto', z=None):
  981. """
  982. Checks if f can be written as the logarithmic derivative of a k(t)-radical.
  983. Explanation
  984. ===========
  985. It differs from is_log_deriv_k_t_radical(fa, fd, DE, Df=False)
  986. for any given fa, fd, DE in that it finds the solution in the
  987. given field not in some (possibly unspecified extension) and
  988. "in_field" with the function name is used to indicate that.
  989. f in k(t) can be written as the logarithmic derivative of a k(t) radical if
  990. there exist n in ZZ and u in k(t) with n, u != 0 such that n*f == Du/u.
  991. Either returns (n, u) or None, which means that f cannot be written as the
  992. logarithmic derivative of a k(t)-radical.
  993. case is one of {'primitive', 'exp', 'tan', 'auto'} for the primitive,
  994. hyperexponential, and hypertangent cases, respectively. If case is 'auto',
  995. it will attempt to determine the type of the derivation automatically.
  996. See also
  997. ========
  998. is_log_deriv_k_t_radical, is_deriv_k
  999. """
  1000. fa, fd = fa.cancel(fd, include=True)
  1001. # f must be simple
  1002. n, s = splitfactor(fd, DE)
  1003. if not s.is_one:
  1004. pass
  1005. z = z or Dummy('z')
  1006. H, b = residue_reduce(fa, fd, DE, z=z)
  1007. if not b:
  1008. # I will have to verify, but I believe that the answer should be
  1009. # None in this case. This should never happen for the
  1010. # functions given when solving the parametric logarithmic
  1011. # derivative problem when integration elementary functions (see
  1012. # Bronstein's book, page 255), so most likely this indicates a bug.
  1013. return None
  1014. roots = [(i, i.real_roots()) for i, _ in H]
  1015. if not all(len(j) == i.degree() and all(k.is_Rational for k in j) for
  1016. i, j in roots):
  1017. # If f is the logarithmic derivative of a k(t)-radical, then all the
  1018. # roots of the resultant must be rational numbers.
  1019. return None
  1020. # [(a, i), ...], where i*log(a) is a term in the log-part of the integral
  1021. # of f
  1022. respolys, residues = list(zip(*roots)) or [[], []]
  1023. # Note: this might be empty, but everything below should work find in that
  1024. # case (it should be the same as if it were [[1, 1]])
  1025. residueterms = [(H[j][1].subs(z, i), i) for j in range(len(H)) for
  1026. i in residues[j]]
  1027. # TODO: finish writing this and write tests
  1028. p = cancel(fa.as_expr()/fd.as_expr() - residue_reduce_derivation(H, DE, z))
  1029. p = p.as_poly(DE.t)
  1030. if p is None:
  1031. # f - Dg will be in k[t] if f is the logarithmic derivative of a k(t)-radical
  1032. return None
  1033. if p.degree(DE.t) >= max(1, DE.d.degree(DE.t)):
  1034. return None
  1035. if case == 'auto':
  1036. case = DE.case
  1037. if case == 'exp':
  1038. wa, wd = derivation(DE.t, DE).cancel(Poly(DE.t, DE.t), include=True)
  1039. with DecrementLevel(DE):
  1040. pa, pd = frac_in(p, DE.t, cancel=True)
  1041. wa, wd = frac_in((wa, wd), DE.t)
  1042. A = parametric_log_deriv(pa, pd, wa, wd, DE)
  1043. if A is None:
  1044. return None
  1045. n, e, u = A
  1046. u *= DE.t**e
  1047. elif case == 'primitive':
  1048. with DecrementLevel(DE):
  1049. pa, pd = frac_in(p, DE.t)
  1050. A = is_log_deriv_k_t_radical_in_field(pa, pd, DE, case='auto')
  1051. if A is None:
  1052. return None
  1053. n, u = A
  1054. elif case == 'base':
  1055. # TODO: we can use more efficient residue reduction from ratint()
  1056. if not fd.is_sqf or fa.degree() >= fd.degree():
  1057. # f is the logarithmic derivative in the base case if and only if
  1058. # f = fa/fd, fd is square-free, deg(fa) < deg(fd), and
  1059. # gcd(fa, fd) == 1. The last condition is handled by cancel() above.
  1060. return None
  1061. # Note: if residueterms = [], returns (1, 1)
  1062. # f had better be 0 in that case.
  1063. n = reduce(ilcm, [i.as_numer_denom()[1] for _, i in residueterms], S.One)
  1064. u = Mul(*[Pow(i, j*n) for i, j in residueterms])
  1065. return (n, u)
  1066. elif case == 'tan':
  1067. raise NotImplementedError("The hypertangent case is "
  1068. "not yet implemented for is_log_deriv_k_t_radical_in_field()")
  1069. elif case in ('other_linear', 'other_nonlinear'):
  1070. # XXX: If these are supported by the structure theorems, change to NotImplementedError.
  1071. raise ValueError("The %s case is not supported in this function." % case)
  1072. else:
  1073. raise ValueError("case must be one of {'primitive', 'exp', 'tan', "
  1074. "'base', 'auto'}, not %s" % case)
  1075. common_denom = reduce(ilcm, [i.as_numer_denom()[1] for i in [j for _, j in
  1076. residueterms]] + [n], S.One)
  1077. residueterms = [(i, j*common_denom) for i, j in residueterms]
  1078. m = common_denom//n
  1079. if common_denom != n*m: # Verify exact division
  1080. raise ValueError("Inexact division")
  1081. u = cancel(u**m*Mul(*[Pow(i, j) for i, j in residueterms]))
  1082. return (common_denom, u)