forked from sympy/sympy
-
Notifications
You must be signed in to change notification settings - Fork 1
/
risch.py
1707 lines (1412 loc) · 63 KB
/
risch.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
"""
The Risch Algorithm for transcendental function integration.
The core algorithms for the Risch algorithm are here. The subproblem
algorithms are in the rde.py and prde.py files for the Risch
Differential Equation solver and the parametric problems solvers,
respectively. All important information concerning the differential extension
for an integrand is stored in a DifferentialExtension object, which in the code
is usually called DE. Throughout the code and Inside the DifferentialExtension
object, the conventions/attribute names are that the base domain is QQ and each
differential extension is x, t0, t1, ..., tn-1 = DE.t. DE.x is the variable of
integration (Dx == 1), DE.D is a list of the derivatives of
x, t1, t2, ..., tn-1 = t, DE.T is the list [x, t1, t2, ..., tn-1], DE.t is the
outer-most variable of the differential extension at the given level (the level
can be adjusted using DE.increment_level() and DE.decrement_level()),
k is the field C(x, t0, ..., tn-2), where C is the constant field. The
numerator of a fraction is denoted by a and the denominator by
d. If the fraction is named f, fa == numer(f) and fd == denom(f).
Fractions are returned as tuples (fa, fd). DE.d and DE.t are used to
represent the topmost derivation and extension variable, respectively.
The docstring of a function signifies whether an argument is in k[t], in
which case it will just return a Poly in t, or in k(t), in which case it
will return the fraction (fa, fd). Other variable names probably come
from the names used in Bronstein's book.
"""
from __future__ import print_function, division
from sympy import real_roots
from sympy.abc import z
from sympy.core.function import Lambda
from sympy.core.numbers import ilcm, oo
from sympy.core.mul import Mul
from sympy.core.power import Pow
from sympy.core.relational import Eq, Ne
from sympy.core.singleton import S
from sympy.core.symbol import Symbol, Dummy
from sympy.core.compatibility import reduce, ordered, xrange
from sympy.integrals.heurisch import _symbols
from sympy.functions import (acos, acot, asin, atan, cos, cot, exp, log,
Piecewise, sin, tan)
from sympy.functions import sinh, cosh, tanh, coth, asinh, acosh , atanh , acoth
from sympy.integrals import Integral, integrate
from sympy.polys import gcd, cancel, PolynomialError, Poly, reduced, RootSum, DomainError
from sympy.utilities.iterables import numbered_symbols
def integer_powers(exprs):
"""
Rewrites a list of expressions as integer multiples of each other.
For example, if you have [x, x/2, x**2 + 1, 2*x/3], then you can rewrite
this as [(x/6) * 6, (x/6) * 3, (x**2 + 1) * 1, (x/6) * 4]. This is useful
in the Risch integration algorithm, where we must write exp(x) + exp(x/2)
as (exp(x/2))**2 + exp(x/2), but not as exp(x) + sqrt(exp(x)) (this is
because only the transcendental case is implemented and we therefore cannot
integrate algebraic extensions). The integer multiples returned by this
function for each term are the smallest possible (their content equals 1).
Returns a list of tuples where the first element is the base term and the
second element is a list of `(item, factor)` terms, where `factor` is the
integer multiplicative factor that must multiply the base term to obtain
the original item.
The easiest way to understand this is to look at an example:
>>> from sympy.abc import x
>>> from sympy.integrals.risch import integer_powers
>>> integer_powers([x, x/2, x**2 + 1, 2*x/3])
[(x/6, [(x, 6), (x/2, 3), (2*x/3, 4)]), (x**2 + 1, [(x**2 + 1, 1)])]
We can see how this relates to the example at the beginning of the
docstring. It chose x/6 as the first base term. Then, x can be written as
(x/2) * 2, so we get (0, 2), and so on. Now only element (x**2 + 1)
remains, and there are no other terms that can be written as a rational
multiple of that, so we get that it can be written as (x**2 + 1) * 1.
"""
# Here is the strategy:
# First, go through each term and determine if it can be rewritten as a
# rational multiple of any of the terms gathered so far.
# cancel(a/b).is_Rational is sufficient for this. If it is a multiple, we
# add its multiple to the dictionary.
terms = {}
for term in exprs:
for j in terms:
a = cancel(term/j)
if a.is_Rational:
terms[j].append((term, a))
break
else:
terms[term] = [(term, S(1))]
# After we have done this, we have all the like terms together, so we just
# need to find a common denominator so that we can get the base term and
# integer multiples such that each term can be written as an integer
# multiple of the base term, and the content of the integers is 1.
newterms = {}
for term in terms:
common_denom = reduce(ilcm, [i.as_numer_denom()[1] for _, i in
terms[term]])
newterm = term/common_denom
newmults = [(i, j*common_denom) for i, j in terms[term]]
newterms[newterm] = newmults
return sorted(iter(newterms.items()), key=lambda item: item[0].sort_key())
class DifferentialExtension(object):
"""
A container for all the information relating to a differential extension.
The attributes of this object are (see also the docstring of __init__):
- f: The original (Expr) integrand.
- x: The variable of integration.
- T: List of variables in the extension.
- D: List of derivations in the extension; corresponds to the elements of T.
- fa: Poly of the numerator of the integrand.
- fd: Poly of the denominator of the integrand.
- Tfuncs: Lambda() representations of each element of T (except for x).
For back-substitution after integration.
- backsubs: A (possibly empty) list of further substitutions to be made on
the final integral to make it look more like the integrand.
- E_K: List of the positions of the exponential extensions in T.
- E_args: The arguments of each of the exponentials in E_K.
- L_K: List of the positions of the logarithmic extensions in T.
- L_args: The arguments of each of the logarithms in L_K.
(See the docstrings of is_deriv_k() and is_log_deriv_k_t_radical() for
more information on E_K, E_args, L_K, and L_args)
- cases: List of string representations of the cases of T.
- t: The top level extension variable, as defined by the current level
(see level below).
- d: The top level extension derivation, as defined by the current
derivation (see level below).
- case: The string representation of the case of self.d.
(Note that self.T and self.D will always contain the complete extension,
regardless of the level. Therefore, you should ALWAYS use DE.t and DE.d
instead of DE.T[-1] and DE.D[-1]. If you want to have a list of the
derivations or variables only up to the current level, use
DE.D[:len(DE.D) + DE.level + 1] and DE.T[:len(DE.T) + DE.level + 1]. Note
that, in particular, the derivation() function does this.)
The following are also attributes, but will probably not be useful other
than in internal use:
- newf: Expr form of fa/fd.
- level: The number (between -1 and -len(self.T)) such that
self.T[self.level] == self.t and self.D[self.level] == self.d.
Use the methods self.increment_level() and self.decrement_level() to change
the current level.
"""
# __slots__ is defined mainly so we can iterate over all the attributes
# of the class easily (the memory use doesn't matter too much, since we
# only create one DifferentialExtension per integration). Also, it's nice
# to have a safeguard when debugging.
__slots__ = ('f', 'x', 'T', 'D', 'fa', 'fd', 'Tfuncs', 'backsubs', 'E_K',
'E_args', 'L_K', 'L_args', 'cases', 'case', 't', 'd', 'newf', 'level',
'ts')
def __init__(self, f=None, x=None, handle_first='log', dummy=True, extension=None, rewrite_complex=False):
"""
Tries to build a transcendental extension tower from f with respect to x.
If it is successful, creates a DifferentialExtension object with, among
others, the attributes fa, fd, D, T, Tfuncs, and backsubs such that
fa and fd are Polys in T[-1] with rational coefficients in T[:-1],
fa/fd == f, and D[i] is a Poly in T[i] with rational coefficients in
T[:i] representing the derivative of T[i] for each i from 1 to len(T).
Tfuncs is a list of Lambda objects for back replacing the functions
after integrating. Lambda() is only used (instead of lambda) to make
them easier to test and debug. Note that Tfuncs corresponds to the
elements of T, except for T[0] == x, but they should be back-substituted
in reverse order. backsubs is a (possibly empty) back-substitution list
that should be applied on the completed integral to make it look more
like the original integrand.
If it is unsuccessful, it raises NotImplementedError.
You can also create an object by manually setting the attributes as a
dictionary to the extension keyword argument. You must include at least
D. Warning, any attribute that is not given will be set to None. The
attributes T, t, d, cases, case, x, and level are set automatically and
do not need to be given. The functions in the Risch Algorithm will NOT
check to see if an attribute is None before using it. This also does not
check to see if the extension is valid (non-algebraic) or even if it is
self-consistent. Therefore, this should only be used for
testing/debugging purposes.
"""
# XXX: If you need to debug this function, set the break point here
if extension:
if 'D' not in extension:
raise ValueError("At least the key D must be included with "
"the extension flag to DifferentialExtension.")
for attr in extension:
setattr(self, attr, extension[attr])
self._auto_attrs()
return
elif f is None or x is None:
raise ValueError("Either both f and x or a manual extension must "
"be given.")
from sympy.integrals.prde import is_deriv_k
if handle_first not in ['log', 'exp']:
raise ValueError("handle_first must be 'log' or 'exp', not %s." %
str(handle_first))
# f will be the original function, self.f might change if we reset
# (e.g., we pull out a constant from an exponential)
self.f = f
self.x = x
self.reset(dummy=dummy)
exp_new_extension, log_new_extension = True, True
if rewrite_complex:
rewritables = {
(sin, cos, cot, tan, sinh, cosh, coth, tanh): exp,
(asin, acos, acot, atan): log,
}
#rewrite the trigonometric components
for candidates, rule in rewritables.items():
self.newf = self.newf.rewrite(candidates, rule)
else:
if any(i.has(x) for i in self.f.atoms(sin, cos, tan, atan, asin, acos)):
raise NotImplementedError("Trigonometric extensions are not "
"supported (yet!)")
def update(seq, atoms, func):
s = set(seq)
new = atoms - s
s = atoms.intersection(s)
s.update(list(filter(func, new)))
return list(s)
exps = set()
pows = set()
numpows = set()
sympows = set()
logs = set()
symlogs = set()
while True:
restart = False
if self.newf.is_rational_function(*self.T):
break
if not exp_new_extension and not log_new_extension:
# We couldn't find a new extension on the last pass, so I guess
# we can't do it.
raise NotImplementedError("Couldn't find an elementary "
"transcendental extension for %s. Try using a " % str(f) +
"manual extension with the extension flag.")
# Pre-preparsing.
#################
# Get all exp arguments, so we can avoid ahead of time doing
# something like t1 = exp(x), t2 = exp(x/2) == sqrt(t1).
# Things like sqrt(exp(x)) do not automatically simplify to
# exp(x/2), so they will be viewed as algebraic. The easiest way
# to handle this is to convert all instances of (a**b)**Rational
# to a**(Rational*b) before doing anything else. Note that the
# _exp_part code can generate terms of this form, so we do need to
# do this at each pass (or else modify it to not do that).
ratpows = [i for i in self.newf.atoms(Pow).union(self.newf.atoms(exp))
if (i.base.is_Pow or i.base.func is exp and i.exp.is_Rational)]
ratpows_repl = [
(i, i.base.base**(i.exp*i.base.exp)) for i in ratpows]
self.backsubs += [(j, i) for i, j in ratpows_repl]
self.newf = self.newf.xreplace(dict(ratpows_repl))
# To make the process deterministic, the args are sorted
# so that functions with smaller op-counts are processed first.
# Ties are broken with the default_sort_key.
# XXX Although the method is deterministic no additional work
# has been done to guarantee that the simplest solution is
# returned and that it would be affected be using different
# variables. Though it is possible that this is the case
# one should know that it has not been done intentionally so
# further improvements may possible.
# TODO: This probably doesn't need to be completely recomputed at
# each pass.
exps = update(exps, self.newf.atoms(exp),
lambda i: i.exp.is_rational_function(*self.T) and
i.exp.has(*self.T))
pows = update(pows, self.newf.atoms(Pow),
lambda i: i.exp.is_rational_function(*self.T) and
i.exp.has(*self.T))
numpows = update(numpows, set(pows),
lambda i: not i.base.has(*self.T))
sympows = update(sympows, set(pows) - set(numpows),
lambda i: i.base.is_rational_function(*self.T) and
not i.exp.is_Integer)
# The easiest way to deal with non-base E powers is to convert them
# into base E, integrate, and then convert back.
for i in ordered(pows):
old = i
new = exp(i.exp*log(i.base))
# If exp is ever changed to automatically reduce exp(x*log(2))
# to 2**x, then this will break. The solution is to not change
# exp to do that :)
if i in sympows:
if i.exp.is_Rational:
raise NotImplementedError("Algebraic extensions are "
"not supported (%s)." % str(i))
# We can add a**b only if log(a) in the extension, because
# a**b == exp(b*log(a)).
basea, based = frac_in(i.base, self.t)
A = is_deriv_k(basea, based, self)
if A is None:
# Nonelementary monomial (so far)
# TODO: Would there ever be any benefit from just
# adding log(base) as a new monomial?
# ANSWER: Yes, otherwise we can't integrate x**x (or
# rather prove that it has no elementary integral)
# without first manually rewriting it as exp(x*log(x))
continue
ans, u, const = A
newterm = exp(i.exp*(log(const) + u))
# Under the current implementation, exp kills terms
# only if they are of the form a*log(x), where a is a
# Number. This case should have already been killed by the
# above tests. Again, if this changes to kill more than
# that, this will break, which maybe is a sign that you
# shouldn't be changing that. Actually, if anything, this
# auto-simplification should be removed. See
# http://groups.google.com/group/sympy/browse_thread/thread/a61d48235f16867f
self.newf = self.newf.xreplace({i: newterm})
elif i not in numpows:
continue
else:
# i in numpows
newterm = new
# TODO: Just put it in self.Tfuncs
self.backsubs.append((new, old))
self.newf = self.newf.xreplace({old: newterm})
exps.append(newterm)
atoms = self.newf.atoms(log)
logs = update(logs, atoms,
lambda i: i.args[0].is_rational_function(*self.T) and
i.args[0].has(*self.T))
symlogs = update(symlogs, atoms,
lambda i: i.has(*self.T) and i.args[0].is_Pow and
i.args[0].base.is_rational_function(*self.T) and
not i.args[0].exp.is_Integer)
# We can handle things like log(x**y) by converting it to y*log(x)
# This will fix not only symbolic exponents of the argument, but any
# non-Integer exponent, like log(sqrt(x)). The exponent can also
# depend on x, like log(x**x).
for i in ordered(symlogs):
# Unlike in the exponential case above, we do not ever
# potentially add new monomials (above we had to add log(a)).
# Therefore, there is no need to run any is_deriv functions
# here. Just convert log(a**b) to b*log(a) and let
# log_new_extension() handle it from there.
lbase = log(i.args[0].base)
logs.append(lbase)
new = i.args[0].exp*lbase
self.newf = self.newf.xreplace({i: new})
self.backsubs.append((new, i))
# remove any duplicates
logs = list(set(logs))
if handle_first == 'exp' or not log_new_extension:
exp_new_extension = self._exp_part(exps, dummy=dummy)
if exp_new_extension is None:
# reset and restart
self.f = self.newf
self.reset(dummy=dummy)
exp_new_extension = True
continue
if handle_first == 'log' or not exp_new_extension:
log_new_extension = self._log_part(logs, dummy=dummy)
self.fa, self.fd = frac_in(self.newf, self.t)
self._auto_attrs()
return
def __getattr__(self, attr):
# Avoid AttributeErrors when debugging
if attr not in self.__slots__:
raise AttributeError("%s has no attribute %s" % (repr(self), repr(attr)))
return None
def _auto_attrs(self):
"""
Set attributes that are generated automatically.
"""
if not self.T:
# i.e., when using the extension flag and T isn't given
self.T = [i.gen for i in self.D]
if not self.x:
self.x = self.T[0]
self.cases = [get_case(d, t) for d, t in zip(self.D, self.T)]
self.level = -1
self.t = self.T[self.level]
self.d = self.D[self.level]
self.case = self.cases[self.level]
def _exp_part(self, exps, dummy=True):
"""
Try to build an exponential extension.
Returns True if there was a new extension, False if there was no new
extension but it was able to rewrite the given exponentials in terms
of the existing extension, and None if the entire extension building
process should be restarted. If the process fails because there is no
way around an algebraic extension (e.g., exp(log(x)/2)), it will raise
NotImplementedError.
"""
from sympy.integrals.prde import is_log_deriv_k_t_radical
new_extension = False
restart = False
expargs = [i.exp for i in exps]
ip = integer_powers(expargs)
for arg, others in ip:
# Minimize potential problems with algebraic substitution
others.sort(key=lambda i: i[1])
arga, argd = frac_in(arg, self.t)
A = is_log_deriv_k_t_radical(arga, argd, self)
if A is not None:
ans, u, n, const = A
# if n is 1 or -1, it's algebraic, but we can handle it
if n == -1:
# This probably will never happen, because
# Rational.as_numer_denom() returns the negative term in
# the numerator. But in case that changes, reduce it to
# n == 1.
n = 1
u **= -1
const *= -1
ans = [(i, -j) for i, j in ans]
if n == 1:
# Example: exp(x + x**2) over QQ(x, exp(x), exp(x**2))
self.newf = self.newf.xreplace({exp(arg): exp(const)*Mul(*[
u**power for u, power in ans])})
self.newf = self.newf.xreplace(dict([(exp(p*exparg),
exp(const*p) * Mul(*[u**power for u, power in ans]))
for exparg, p in others]))
# TODO: Add something to backsubs to put exp(const*p)
# back together.
continue
else:
# Bad news: we have an algebraic radical. But maybe we
# could still avoid it by choosing a different extension.
# For example, integer_powers() won't handle exp(x/2 + 1)
# over QQ(x, exp(x)), but if we pull out the exp(1), it
# will. Or maybe we have exp(x + x**2/2), over
# QQ(x, exp(x), exp(x**2)), which is exp(x)*sqrt(exp(x**2)),
# but if we use QQ(x, exp(x), exp(x**2/2)), then they will
# all work.
#
# So here is what we do: If there is a non-zero const, pull
# it out and retry. Also, if len(ans) > 1, then rewrite
# exp(arg) as the product of exponentials from ans, and
# retry that. If const == 0 and len(ans) == 1, then we
# assume that it would have been handled by either
# integer_powers() or n == 1 above if it could be handled,
# so we give up at that point. For example, you can never
# handle exp(log(x)/2) because it equals sqrt(x).
if const or len(ans) > 1:
rad = Mul(*[term**(power/n) for term, power in ans])
self.newf = self.newf.xreplace(dict((exp(p*exparg),
exp(const*p)*rad) for exparg, p in others))
self.newf = self.newf.xreplace(dict(list(zip(reversed(self.T),
reversed([f(self.x) for f in self.Tfuncs])))))
restart = True
break
else:
# TODO: give algebraic dependence in error string
raise NotImplementedError("Cannot integrate over "
"algebraic extensions.")
else:
arga, argd = frac_in(arg, self.t)
darga = (argd*derivation(Poly(arga, self.t), self) -
arga*derivation(Poly(argd, self.t), self))
dargd = argd**2
darga, dargd = darga.cancel(dargd, include=True)
darg = darga.as_expr()/dargd.as_expr()
self.t = next(self.ts)
self.T.append(self.t)
self.E_args.append(arg)
self.E_K.append(len(self.T) - 1)
self.D.append(darg.as_poly(self.t, expand=False)*Poly(self.t,
self.t, expand=False))
if dummy:
i = Dummy("i")
else:
i = Symbol('i')
self.Tfuncs = self.Tfuncs + [Lambda(i, exp(arg.subs(self.x, i)))]
self.newf = self.newf.xreplace(
dict((exp(exparg), self.t**p) for exparg, p in others))
new_extension = True
if restart:
return None
return new_extension
def _log_part(self, logs, dummy=True):
"""
Try to build a logarithmic extension.
Returns True if there was a new extension and False if there was no new
extension but it was able to rewrite the given logarithms in terms
of the existing extension. Unlike with exponential extensions, there
is no way that a logarithm is not transcendental over and cannot be
rewritten in terms of an already existing extension in a non-algebraic
way, so this function does not ever return None or raise
NotImplementedError.
"""
from sympy.integrals.prde import is_deriv_k
new_extension = False
logargs = [i.args[0] for i in logs]
for arg in ordered(logargs):
# The log case is easier, because whenever a logarithm is algebraic
# over the base field, it is of the form a1*t1 + ... an*tn + c,
# which is a polynomial, so we can just replace it with that.
# In other words, we don't have to worry about radicals.
arga, argd = frac_in(arg, self.t)
A = is_deriv_k(arga, argd, self)
if A is not None:
ans, u, const = A
newterm = log(const) + u
self.newf = self.newf.xreplace({log(arg): newterm})
continue
else:
arga, argd = frac_in(arg, self.t)
darga = (argd*derivation(Poly(arga, self.t), self) -
arga*derivation(Poly(argd, self.t), self))
dargd = argd**2
darg = darga.as_expr()/dargd.as_expr()
self.t = next(self.ts)
self.T.append(self.t)
self.L_args.append(arg)
self.L_K.append(len(self.T) - 1)
self.D.append(cancel(darg.as_expr()/arg).as_poly(self.t,
expand=False))
if dummy:
i = Dummy("i")
else:
i = Symbol('i')
self.Tfuncs = self.Tfuncs + [Lambda(i, log(arg.subs(self.x, i)))]
self.newf = self.newf.xreplace({log(arg): self.t})
new_extension = True
return new_extension
@property
def _important_attrs(self):
"""
Returns some of the more important attributes of self.
Used for testing and debugging purposes.
The attributes are (fa, fd, D, T, Tfuncs, backsubs, E_K, E_args,
L_K, L_args).
"""
# XXX: This might be easier to read as a dict or something
# Maybe a named tuple.
return (self.fa, self.fd, self.D, self.T, self.Tfuncs,
self.backsubs, self.E_K, self.E_args, self.L_K, self.L_args)
# TODO: Implement __repr__
def __str__(self):
return str(self._important_attrs)
def reset(self, dummy=True):
"""
Reset self to an initial state. Used by __init__.
"""
self.t = self.x
self.T = [self.x]
self.D = [Poly(1, self.x)]
self.level = -1
self.L_K, self.E_K, self.L_args, self.E_args = [], [], [], []
if dummy:
self.ts = numbered_symbols('t', cls=Dummy)
else:
# For testing
self.ts = numbered_symbols('t')
# For various things that we change to make things work that we need to
# change back when we are done.
self.backsubs = []
self.Tfuncs = []
self.newf = self.f
def increment_level(self):
"""
Increment the level of self.
This makes the working differential extension larger. self.level is
given relative to the end of the list (-1, -2, etc.), so we don't need
do worry about it when building the extension.
"""
if self.level >= -1:
raise ValueError("The level of the differential extension cannot "
"be incremented any further.")
self.level += 1
self.t = self.T[self.level]
self.d = self.D[self.level]
self.case = self.cases[self.level]
return None
def decrement_level(self):
"""
Decrease the level of self.
This makes the working differential extension smaller. self.level is
given relative to the end of the list (-1, -2, etc.), so we don't need
do worry about it when building the extension.
"""
if self.level <= -len(self.T):
raise ValueError("The level of the differential extension cannot "
"be decremented any further.")
self.level -= 1
self.t = self.T[self.level]
self.d = self.D[self.level]
self.case = self.cases[self.level]
return None
class DecrementLevel(object):
"""
A context manager for decrementing the level of a DifferentialExtension.
"""
__slots__ = ('DE',)
def __init__(self, DE):
self.DE = DE
return
def __enter__(self):
self.DE.decrement_level()
def __exit__(self, exc_type, exc_value, traceback):
self.DE.increment_level()
class NonElementaryIntegralException(Exception):
"""
Exception used by subroutines within the Risch algorithm to indicate to one
another that the function being integrated does not have an elementary
integral in the given differential field.
"""
# TODO: Rewrite algorithms below to use this (?)
# TODO: Pass through information about why the integral was nonelementary,
# and store that in the resulting NonElementaryIntegral somehow.
pass
def gcdex_diophantine(a, b, c):
"""
Extended Euclidean Algorithm, Diophantine version.
Given a, b in K[x] and c in (a, b), the ideal generated by a and b,
return (s, t) such that s*a + t*b == c and either s == 0 or s.degree()
< b.degree().
"""
# Extended Euclidean Algorithm (Diophantine Version) pg. 13
# TODO: This should go in densetools.py.
# XXX: Bettter name?
s, g = a.half_gcdex(b)
q = c.exquo(g) # Inexact division means c is not in (a, b)
s = q*s
if not s.is_zero and b.degree() >= b.degree():
q, s = s.div(b)
t = (c - s*a).exquo(b)
return (s, t)
def frac_in(f, t, **kwargs):
"""
Returns the tuple (fa, fd), where fa and fd are Polys in t.
This is a common idiom in the Risch Algorithm functions, so we abstract
it out here. f should be a basic expresion, a Poly, or a tuple (fa, fd),
where fa and fd are either basic expressions or Polys, and f == fa/fd.
**kwargs are applied to Poly.
"""
cancel = kwargs.pop('cancel', False)
if type(f) is tuple:
fa, fd = f
f = fa.as_expr()/fd.as_expr()
fa, fd = f.as_expr().as_numer_denom()
fa, fd = fa.as_poly(t, **kwargs), fd.as_poly(t, **kwargs)
if cancel:
fa, fd = fa.cancel(fd, include=True)
if fa is None or fd is None:
raise ValueError("Could not turn %s into a fraction in %s." % (f, t))
return (fa, fd)
def as_poly_1t(p, t, z):
"""
(Hackish) way to convert an element p of K[t, 1/t] to K[t, z].
In other words, z == 1/t will be a dummy variable that Poly can handle
better.
See issue 5131.
Examples
========
>>> from sympy import Symbol, random_poly
>>> from sympy.integrals.risch import as_poly_1t
>>> from sympy.abc import x, z
>>> p1 = random_poly(x, 10, -10, 10)
>>> p2 = random_poly(x, 10, -10, 10)
>>> p = p1 + p2.subs(x, 1/x)
>>> as_poly_1t(p, x, z).as_expr().subs(z, 1/x) == p
True
"""
# TODO: Use this on the final result. That way, we can avoid answers like
# (...)*exp(-x).
pa, pd = frac_in(p, t, cancel=True)
if not pd.is_monomial:
# XXX: Is there a better Poly exception that we could raise here?
# Either way, if you see this (from the Risch Algorithm) it indicates
# a bug.
raise PolynomialError("%s is not an element of K[%s, 1/%s]." % (p, t, t))
d = pd.degree(t)
one_t_part = pa.slice(0, d + 1)
r = pd.degree() - pa.degree()
t_part = pa - one_t_part
try:
t_part = t_part.to_field().exquo(pd)
except DomainError as e:
# issue 4950
raise NotImplementedError(e)
# Compute the negative degree parts.
one_t_part = Poly.from_list(reversed(one_t_part.rep.rep), *one_t_part.gens,
domain=one_t_part.domain)
if 0 < r < oo:
one_t_part *= Poly(t**r, t)
one_t_part = one_t_part.replace(t, z) # z will be 1/t
if pd.nth(d):
one_t_part *= Poly(1/pd.nth(d), z, expand=False)
ans = t_part.as_poly(t, z, expand=False) + one_t_part.as_poly(t, z,
expand=False)
return ans
def derivation(p, DE, coefficientD=False, basic=False):
"""
Computes Dp.
Given the derivation D with D = d/dx and p is a polynomial in t over
K(x), return Dp.
If coefficientD is True, it computes the derivation kD
(kappaD), which is defined as kD(sum(ai*Xi**i, (i, 0, n))) ==
sum(Dai*Xi**i, (i, 1, n)) (Definition 3.2.2, page 80). X in this case is
T[-1], so coefficientD computes the derivative just with respect to T[:-1],
with T[-1] treated as a constant.
If basic=True, the returns a Basic expression. Elements of D can still be
instances of Poly.
"""
if basic:
r = 0
else:
r = Poly(0, DE.t)
t = DE.t
if coefficientD:
if DE.level <= -len(DE.T):
# 'base' case, the answer is 0.
return r
DE.decrement_level()
D = DE.D[:len(DE.D) + DE.level + 1]
T = DE.T[:len(DE.T) + DE.level + 1]
for d, v in zip(D, T):
pv = p.as_poly(v)
if pv is None or basic:
pv = p.as_expr()
if basic:
r += d.as_expr()*pv.diff(v)
else:
r += (d*pv.diff(v)).as_poly(t)
if basic:
r = cancel(r)
if coefficientD:
DE.increment_level()
return r
def get_case(d, t):
"""
Returns the type of the derivation d.
Returns one of {'exp', 'tan', 'base', 'primitive', 'other_linear',
'other_nonlinear'}.
"""
if not d.has(t):
if d.is_one:
return 'base'
return 'primitive'
if d.rem(Poly(t, t)).is_zero:
return 'exp'
if d.rem(Poly(1 + t**2, t)).is_zero:
return 'tan'
if d.degree(t) > 1:
return 'other_nonlinear'
return 'other_linear'
def splitfactor(p, DE, coefficientD=False, z=None):
"""
Splitting factorization.
Given a derivation D on k[t] and p in k[t], return (p_n, p_s) in
k[t] x k[t] such that p = p_n*p_s, p_s is special, and each square
factor of p_n is normal.
Page. 100
"""
kinv = [1/x for x in DE.T[:DE.level]]
if z:
kinv.append(z)
One = Poly(1, DE.t, domain=p.get_domain())
Dp = derivation(p, DE, coefficientD=coefficientD)
# XXX: Is this right?
if p.is_zero:
return (p, One)
if not p.has(DE.t):
s = p.as_poly(*kinv).gcd(Dp.as_poly(*kinv)).as_poly(DE.t)
n = p.exquo(s)
return (n, s)
if not Dp.is_zero:
h = p.gcd(Dp).to_field()
g = p.gcd(p.diff(DE.t)).to_field()
s = h.exquo(g)
if s.degree(DE.t) == 0:
return (p, One)
q_split = splitfactor(p.exquo(s), DE, coefficientD=coefficientD)
return (q_split[0], q_split[1]*s)
else:
return (p, One)
def splitfactor_sqf(p, DE, coefficientD=False, z=None, basic=False):
"""
Splitting Square-free Factorization
Given a derivation D on k[t] and p in k[t], returns (N1, ..., Nm)
and (S1, ..., Sm) in k[t]^m such that p =
(N1*N2**2*...*Nm**m)*(S1*S2**2*...*Sm**m) is a splitting
factorization of p and the Ni and Si are square-free and coprime.
"""
# TODO: This algorithm appears to be faster in every case
# TODO: Verify this and splitfactor() for multiple extensions
kkinv = [1/x for x in DE.T[:DE.level]] + DE.T[:DE.level]
if z:
kkinv = [z]
S = []
N = []
p_sqf = p.sqf_list_include()
if p.is_zero:
return (((p, 1),), ())
for pi, i in p_sqf:
Si = pi.as_poly(*kkinv).gcd(derivation(pi, DE,
coefficientD=coefficientD,basic=basic).as_poly(*kkinv)).as_poly(DE.t)
pi = Poly(pi, DE.t)
Si = Poly(Si, DE.t)
Ni = pi.exquo(Si)
if not Si.is_one:
S.append((Si, i))
if not Ni.is_one:
N.append((Ni, i))
return (tuple(N), tuple(S))
def canonical_representation(a, d, DE):
"""
Canonical Representation.
Given a derivation D on k[t] and f = a/d in k(t), return (f_p, f_s,
f_n) in k[t] x k(t) x k(t) such that f = f_p + f_s + f_n is the
canonical representation of f (f_p is a polynomial, f_s is reduced
(has a special denominator), and f_n is simple (has a normal
denominator).
"""
# Make d monic
l = Poly(1/d.LC(), DE.t)
a, d = a.mul(l), d.mul(l)
q, r = a.div(d)
dn, ds = splitfactor(d, DE)
b, c = gcdex_diophantine(dn.as_poly(DE.t), ds.as_poly(DE.t), r.as_poly(DE.t))
b, c = b.as_poly(DE.t), c.as_poly(DE.t)
return (q, (b, ds), (c, dn))
def hermite_reduce(a, d, DE):
"""
Hermite Reduction - Mack's Linear Version.
Given a derivation D on k(t) and f = a/d in k(t), returns g, h, r in
k(t) such that f = Dg + h + r, h is simple, and r is reduced.
"""
# Make d monic
l = Poly(1/d.LC(), DE.t)
a, d = a.mul(l), d.mul(l)
fp, fs, fn = canonical_representation(a, d, DE)
a, d = fn
l = Poly(1/d.LC(), DE.t)
a, d = a.mul(l), d.mul(l)
ga = Poly(0, DE.t)
gd = Poly(1, DE.t)
dd = derivation(d, DE)
dm = gcd(d, dd).as_poly(DE.t)
ds, r = d.div(dm)
while dm.degree(DE.t)>0:
ddm = derivation(dm, DE)
dm2 = gcd(dm, ddm)
dms, r = dm.div(dm2)
ds_ddm = ds.mul(ddm)
ds_ddm_dm, r = ds_ddm.div(dm)
b, c = gcdex_diophantine(-ds_ddm_dm.as_poly(DE.t), dms.as_poly(DE.t), a.as_poly(DE.t))
b, c = b.as_poly(DE.t), c.as_poly(DE.t)
db = derivation(b, DE).as_poly(DE.t)
ds_dms, r = ds.div(dms)
a = c.as_poly(DE.t) - db.mul(ds_dms).as_poly(DE.t)
ga = ga*dm + b*gd
gd = gd*dm
ga, gd = ga.cancel(gd, include=True)
dm = dm2
d = ds
q, r = a.div(d)
ga, gd = ga.cancel(gd, include=True)
r, d = r.cancel(d, include=True)