Recent improvements to ore_algebra

Marc Mezzarobba
CNRS, Sorbonne Université

FastRelax Final Workshop, Lyon, 2019-05-23

ore_algebra

Sage implementation of Ore polynomials and D-finite functions

Source code: https://github.com/mkauers/ore_algebra

$ sage -pip install git+https://github.com/mkauers/ore_algebra/

Documentation: http://www.algebra.uni-linz.ac.at/people/mkauers/ore_algebra

History & Contributors

  • Kauers, Jaroschek, Johansson 2013: Initial implementation
  • Mezzarobba ~2015–: Numerical solutions of ODEs
  • Kauers 2015–, Mezzarobba 2018–: Multivariate operators
  • Schwaiger 2017, Hofstadler 2017–2018: “D-finite function” objects

Pre-existing features

  • Basic arithmetic (diff, shift, qdiff, qshift, custom)
  • Gcrd, lclm, D-finite closure properties
  • Conversions (diff. eq. ↔ rec., etc.)
  • Polynomial & rational solutions
  • First-order right factors
  • Formal solutions at singularities
  • Desingularization
  • Guessing
  • Numerical solutions & connection matrices

Recent improvements

  • Multivariate closure properties
  • Creative telescoping
  • D-finite functions (alpha quality)
  • Numerics: faster, better support for reg. sing. points, large operators, high precision

Differential operators: Basic arithmetic

In [1]:
from ore_algebra import OreAlgebra
In [2]:
Pol.<x> = PolynomialRing(QQ)
Dop.<Dx> = OreAlgebra(Pol)
In [3]:
Dop
Out[3]:
Univariate Ore algebra in Dx over Univariate Polynomial Ring in x over Rational Field
In [4]:
Dx*x
Out[4]:
x*Dx + 1
In [5]:
Dx(sin(x))
Out[5]:
cos(x)
In [6]:
(x*Dx)(sin(x))
Out[6]:
x*cos(x)

Recurrence operators

In [7]:
Poln.<n> = PolynomialRing(GF(17))
Rop.<Sn> = OreAlgebra(Poln)
Rop
Out[7]:
Univariate Ore algebra in Sn over Univariate Polynomial Ring in n over Finite Field of size 17
In [8]:
Sn*n
Out[8]:
(n + 1)*Sn
In [9]:
(Sn - n)**10
Out[9]:
Sn^10 + (7*n + 6)*Sn^9 + (11*n^2 + 3*n + 2)*Sn^8 + (16*n^3 + 15*n^2 + 4*n + 2)*Sn^7 + (6*n^4 + 4*n^3 + 2*n^2 + 8*n + 13)*Sn^6 + (3*n^5 + 12*n^4 + 13*n^3 + 10*n^2 + 3*n + 9)*Sn^5 + (6*n^6 + 4*n^5 + 16*n^4 + 10*n^3 + 11*n^2 + 10*n + 3)*Sn^4 + (16*n^7 + 15*n^6 + 7*n^5 + 4*n^4 + 5*n^3 + 16*n^2 + 10*n + 3)*Sn^3 + (11*n^8 + 3*n^7 + 8*n^6 + 6*n^5 + 16*n^4 + 12*n^3 + 3*n^2 + 1)*Sn^2 + (7*n^9 + 6*n^8 + 16*n^7 + 11*n^6 + 3*n^5 + 11*n^4 + 16*n^3 + 6*n^2 + 7*n + 16)*Sn + n^10

⚐ D-Finite functions as objects

C. Hofstadler, experimental

In [10]:
from ore_algebra.dfinite_function import DFiniteFunctionRing, UnivariateDFiniteFunction
In [11]:
Fun = DFiniteFunctionRing(Dop);  Fun
Out[11]:
Ring of D-finite functions over Univariate Polynomial Ring in x over Rational Field
In [12]:
my_ai = Fun(airy_ai(x))
In [13]:
my_ai[10]
Out[13]:
-1/136080*3^(2/3)/gamma(1/3)
In [14]:
my_ai(1)
Out[14]:
[0.13529241631288141552414742351546630617494414298833 +/- 1.49e-51]
In [15]:
my_ai(x^2)
Out[15]:
Univariate D-finite function defined by the annihilating operator x*Dx^2 - Dx - 4*x^5 and the coefficient sequence defined by (n^10 + 22*n^9 + 186*n^8 + 708*n^7 + 777*n^6 - 2562*n^5 - 7876*n^4 - 3928*n^3 + 6912*n^2 + 5760*n)*Sn^6 - 4*n^8 - 48*n^7 - 168*n^6 + 924*n^4 + 1008*n^3 - 752*n^2 - 960*n and {0: 1/3*3^(1/3)/gamma(2/3), 1: 0, 2: -1/3*3^(2/3)/gamma(1/3), 3: 0, 4: 0, 5: 0, 6: 1/18*3^(1/3)/gamma(2/3), 7: 0, 8: -1/36*3^(2/3)/gamma(1/3)}
In [16]:
from ore_algebra.examples.stdfun import mittag_leffler_e
ivp = mittag_leffler_e(1/3, 1)
ivp
Out[16]:
IVP(pt=0, dop=1/3*Dx^3 - x^2*Dx^2 - 4*x*Dx - 2, ini=[1, 1/gamma(4/3), 1/gamma(5/3)])
In [17]:
my_e = UnivariateDFiniteFunction(Fun, ivp.dop, ivp.ini)
In [18]:
my_ai*my_e
Out[18]:
Univariate D-finite function defined by the annihilating operator (-81*x^9 - 414*x^6 - 52*x^3 - 24)*Dx^6 + (486*x^11 + 3213*x^8 + 2796*x^5 + 300*x^2)*Dx^5 + (-729*x^13 - 3483*x^10 + 5148*x^7 - 3732*x^4 + 2064*x)*Dx^4 + (-6075*x^12 - 47574*x^9 - 33384*x^6 - 7464*x^3 + 1968)*Dx^3 + (1458*x^14 + 4293*x^11 - 53604*x^8 - 2244*x^5 - 41088*x^2)*Dx^2 + (7047*x^13 + 60075*x^10 + 71808*x^7 + 136452*x^4 - 31872*x)*Dx - 729*x^15 - 2430*x^12 + 31752*x^9 - 22028*x^6 + 46720*x^3 - 2160 and the coefficient sequence defined by (-24*n^33 - 7848*n^32 - 1217904*n^31 - 119265120*n^30 - 8265240096*n^29 - 430836655392*n^28 - 17522874893376*n^27 - 569149922442240*n^26 - 14980851481012560*n^25 - 322206414583126320*n^24 - 5678126299596310560*n^23 - 81721683804628296000*n^22 - 949652588246790733920*n^21 - 8681713167577264503840*n^20 - 58748680805012448505920*n^19 - 242135956214206048051200*n^18 + 108488260436315928944040*n^17 + 10476807873210919884061080*n^16 + 84066023678942093501383440*n^15 + 328407072430892969322928800*n^14 + 161957855017977822122938944*n^13 - 5556526268492872609705821312*n^12 - 30808184921580886798491192576*n^11 - 61848743568136626694433518080*n^10 + 96516578682202812595845129216*n^9 + 853030845899725268959565893632*n^8 + 1755693662806275221362303696896*n^7 - 191121075488552116917636956160*n^6 - 6918611018576047183218416025600*n^5 - 10305536133084558124200099840000*n^4 - 651883440675078087990312960000*n^3 + 9710692999529424209156505600000*n^2 + 5748846330021297285758976000000*n)*Sn^21 + (-52*n^33 - 15768*n^32 - 2266192*n^31 - 205228464*n^30 - 13131424864*n^29 - 630793023888*n^28 - 23590821575184*n^27 - 702719853964272*n^26 - 16907977517369640*n^25 - 331025204263890240*n^24 - 5279792730717476400*n^23 - 68200672898298404880*n^22 - 701622354982097190720*n^21 - 5530245187115752498800*n^20 - 30124828013909392092720*n^19 - 68576800943396360742480*n^18 + 571745721721581523264380*n^17 + 7492473681331898408673240*n^16 + 41649036729364135964121280*n^15 + 99768141174052190820326400*n^14 - 273482617801028411195658368*n^13 - 3270682878463161988366755072*n^12 - 11711448970574997636645705728*n^11 - 10059256312341487086991773696*n^10 + 76801979904398494328865844224*n^9 + 322846351399709844933782550528*n^8 + 422861006998310521457169874944*n^7 - 488850815277845481558737092608*n^6 - 2352232797247927527470849064960*n^5 - 2540273963213858698740695040000*n^4 + 482772293714613829168005120000*n^3 + 2719501179843522132317306880000*n^2 + 1381740258268276716050841600000*n)*Sn^18 + (-414*n^33 - 117678*n^32 - 15852576*n^31 - 1345594824*n^30 - 80700092856*n^29 - 3633815655240*n^28 - 127399183240896*n^27 - 3557731583055600*n^26 - 80245669721651460*n^25 - 1472279089838630580*n^24 - 21987171811970146080*n^23 - 265389681264853335840*n^22 - 2539134234244059547320*n^21 - 18385725342062668599720*n^20 - 88138188290507684371200*n^19 - 108534146024035898174160*n^18 + 2443632611990285662853490*n^17 + 24833637843205336319220210*n^16 + 121245809437122843644887680*n^15 + 227257412847909330770915880*n^14 - 1106099032265374878229720416*n^13 - 9933584593498812982942176672*n^12 - 31972810234198476727554597504*n^11 - 18761909942962878443753678976*n^10 + 235289950832015737087257086976*n^9 + 902805274335383262934308679680*n^8 + 1108702158744622782158950680576*n^7 - 1476049778974732201091769876480*n^6 - 6591499425554580260012924928000*n^5 - 7038297709246683060184350720000*n^4 + 1380775208407863970913648640000*n^3 + 7639985744558539462921420800000*n^2 + 3899687418093882578042880000000*n)*Sn^15 + (-81*n^33 - 18900*n^32 - 1999656*n^31 - 124129380*n^30 - 4754325036*n^29 - 94831186068*n^28 + 798409333576*n^27 + 136202820180876*n^26 + 5626859890107690*n^25 + 150244980230360340*n^24 + 2958076803935621160*n^23 + 44707003575753198420*n^22 + 522317644089199847700*n^21 + 4622690166053310278820*n^20 + 28843536943661955489720*n^19 + 94346507950084965030660*n^18 - 297866764199743660754145*n^17 - 6149478723439194914314320*n^16 - 39608446972567185637607040*n^15 - 117688111266949538549650560*n^14 + 137113268102502964083036576*n^13 + 2913079673575808379081154560*n^12 + 12078165530462456499372270336*n^11 + 15374326616664284229238287360*n^10 - 63952753134498320456998269184*n^9 - 329179698517056051282299154432*n^8 - 512742818335596794644189188096*n^7 + 357125936983452978364669722624*n^6 + 2474569756741412395763355156480*n^5 + 3018711856234964620121210880000*n^4 - 254978118045128898082897920000*n^3 - 3064821762415666751880560640000*n^2 - 1655071469079834867322060800000*n)*Sn^12 + (486*n^32 + 118503*n^31 + 13524948*n^30 + 958340016*n^29 + 47116511700*n^28 + 1699711744468*n^27 + 46308339416268*n^26 + 962659458170280*n^25 + 15123563382620160*n^24 + 171523512133221690*n^23 + 1186100656517839140*n^22 - 37805196339922440*n^21 - 119715775715274835860*n^20 - 1652851915028270167500*n^19 - 12217664497859617376700*n^18 - 42469132523979420144360*n^17 + 126657823154351913282330*n^16 + 2472822893446539301290135*n^15 + 13578493106461078407258360*n^14 + 24139238444251948125814680*n^13 - 138749005788261364047779136*n^12 - 1053511283509579824159117888*n^11 - 2542721346034360611629054208*n^10 + 2052017446773332083506041984*n^9 + 28003740422874185244995320320*n^8 + 61582361453291997219785870592*n^7 - 480533425570793441674487808*n^6 - 226204572707464956969431900160*n^5 - 341180558039329149831045120000*n^4 - 22763074540137562621992960000*n^3 + 316325128579097767572602880000*n^2 + 186360211691558827956633600000*n)*Sn^9 + (-729*n^31 - 176661*n^30 - 20099340*n^29 - 1425944763*n^28 - 70643005308*n^27 - 2592496126863*n^26 - 72915039882360*n^25 - 1602498019959135*n^24 - 27763983229969350*n^23 - 378691305544396845*n^22 - 4005949233962800800*n^21 - 31524964604749905705*n^20 - 163123050618207140940*n^19 - 252729347108177982645*n^18 + 4264442956387338816600*n^17 + 44947121910163022424555*n^16 + 212775455886957615862935*n^15 + 314133074661371762672550*n^14 - 2470847682946988384781300*n^13 - 17927143197524940083201640*n^12 - 46978390051656573808514736*n^11 + 13491956267542245956894496*n^10 + 462244206112205932904523840*n^9 + 1276871492416154316862946688*n^8 + 775173181989109336044208128*n^7 - 3421420839768170941727404032*n^6 - 8336006423977666738296176640*n^5 - 5406695301365082643599360000*n^4 + 4263545613333384992194560000*n^3 + 7555321040085118538219520000*n^2 + 2884275787535563515494400000*n)*Sn^6 + (1458*n^29 + 329265*n^28 + 34682661*n^27 + 2260935909*n^26 + 102022226865*n^25 + 3374719185825*n^24 + 84459402937845*n^23 + 1624738560888105*n^22 + 24090362505771705*n^21 + 271805463177566175*n^20 + 2238594275860798635*n^19 + 11830705538422983615*n^18 + 16440329540618410995*n^17 - 340948354047289111125*n^16 - 3395596953711720654585*n^15 - 14616943639216277213565*n^14 - 10632331877126123028735*n^13 + 230072441021389691864100*n^12 + 1274975782424092590278580*n^11 + 2318852894373189748894320*n^10 - 4894627815672088825635408*n^9 - 34165273162434107755314240*n^8 - 59650619585806830886883136*n^7 + 28654842091496143828271616*n^6 + 258699952696179776376453120*n^5 + 320170609559343027732480000*n^4 - 27066930120236510300160000*n^3 - 317194158035945496821760000*n^2 - 168348741731156341555200000*n)*Sn^3 - 729*n^27 - 157464*n^26 - 15779205*n^25 - 972340200*n^24 - 41152593105*n^23 - 1264380020280*n^22 - 29021590766025*n^21 - 503177814196200*n^20 - 6550692787196415*n^19 - 62011070347323240*n^18 - 386150232022854975*n^17 - 950333662778830200*n^16 + 8397278950113060165*n^15 + 105798884663200766040*n^14 + 499294024107256828125*n^13 + 406201588063165369800*n^12 - 8081520464744734885020*n^11 - 43218065849227900693920*n^10 - 60308274941946097981200*n^9 + 247545163236333885436800*n^8 + 1148925647418746051335104*n^7 + 1122767140648429530468864*n^6 - 2853333449859192239646720*n^5 - 7415089937678104381440000*n^4 - 2697145321907074590720000*n^3 + 6087484712019018792960000*n^2 + 4469435621181141811200000*n and {0: 1/3*3^(1/3)/gamma(2/3), 1: -1/3*3^(2/3)/gamma(1/3) + 1/3*3^(1/3)/(gamma(4/3)*gamma(2/3)), 2: -1/3*3^(2/3)/(gamma(4/3)*gamma(1/3)) + 1/3*3^(1/3)/(gamma(5/3)*gamma(2/3)), 3: 7/18*3^(1/3)/gamma(2/3) - 1/3*3^(2/3)/(gamma(5/3)*gamma(1/3)), 4: -13/36*3^(2/3)/gamma(1/3) + 11/36*3^(1/3)/(gamma(4/3)*gamma(2/3)), 5: -5/18*3^(2/3)/(gamma(4/3)*gamma(1/3)) + 23/90*3^(1/3)/(gamma(5/3)*gamma(2/3)), 6: 121/540*3^(1/3)/gamma(2/3) - 41/180*3^(2/3)/(gamma(5/3)*gamma(1/3)), 7: -295/1512*3^(2/3)/gamma(1/3) + 1139/7560*3^(1/3)/(gamma(4/3)*gamma(2/3)), 8: -389/3024*3^(2/3)/(gamma(4/3)*gamma(1/3)) + 119/1080*3^(1/3)/(gamma(5/3)*gamma(2/3)), 9: 3313/38880*3^(1/3)/gamma(2/3) - 349/3780*3^(2/3)/(gamma(5/3)*gamma(1/3)), 10: -1363/19440*3^(2/3)/gamma(1/3) + 1999/38880*3^(1/3)/(gamma(4/3)*gamma(2/3)), 11: -2263/54432*3^(2/3)/(gamma(4/3)*gamma(1/3)) + 72901/2138400*3^(1/3)/(gamma(5/3)*gamma(2/3)), 12: 24737/1026432*3^(1/3)/gamma(2/3) - 81157/2993760*3^(2/3)/(gamma(5/3)*gamma(1/3)), 13: -400297/21228480*3^(2/3)/gamma(1/3) + 1563307/116756640*3^(1/3)/(gamma(4/3)*gamma(2/3)), 14: -766349/74299680*3^(2/3)/(gamma(4/3)*gamma(1/3)) + 1477331/179625600*3^(1/3)/(gamma(5/3)*gamma(2/3)), 15: 5835271/1077753600*3^(1/3)/gamma(2/3) - 3643349/583783200*3^(2/3)/(gamma(5/3)*gamma(1/3)), 16: -20629681/5094835200*3^(2/3)/gamma(1/3) + 78854057/28021593600*3^(1/3)/(gamma(4/3)*gamma(2/3)), 17: -10572301/5094835200*3^(2/3)/(gamma(4/3)*gamma(1/3)) + 7436603/4580452800*3^(1/3)/(gamma(5/3)*gamma(2/3)), 18: 332256637/329792601600*3^(1/3)/gamma(2/3) - 1125195523/952734182400*3^(2/3)/(gamma(5/3)*gamma(1/3)), 19: -1260215623/1742433638400*3^(2/3)/gamma(1/3) + 80668819493/162917545190400*3^(1/3)/(gamma(4/3)*gamma(2/3)), 20: -8560868333/24394070937600*3^(2/3)/(gamma(4/3)*gamma(1/3)) + 127635493/471132288000*3^(1/3)/(gamma(5/3)*gamma(2/3)), 21: 22139176921/138512892672000*3^(1/3)/gamma(2/3) - 2162693257679/11404228163328000*3^(2/3)/(gamma(5/3)*gamma(1/3)), 22: -17797586453/161000868188160*3^(2/3)/gamma(1/3) + 29105422369/388780505568000*3^(1/3)/(gamma(4/3)*gamma(2/3)), 23: -3295558349/64400347275264*3^(2/3)/(gamma(4/3)*gamma(1/3)) + 17135676017/438047023075200*3^(1/3)/(gamma(5/3)*gamma(2/3)), 24: 241225777519/10922730964992000*3^(1/3)/gamma(2/3) - 83210444770261/3147566973078528000*3^(2/3)/(gamma(5/3)*gamma(1/3)), 25: -7129346024401/483002604564480000*3^(2/3)/gamma(1/3) + 187269180969109/18885401838471168000*3^(1/3)/(gamma(4/3)*gamma(2/3)), 26: -41167459351903/6279033859338240000*3^(2/3)/(gamma(4/3)*gamma(1/3)) + 271991502318703/54668268479784960000*3^(1/3)/(gamma(5/3)*gamma(2/3))}
In [19]:
my_e(1)
Out[19]:
[7.28355749378782728852876426824018405440330223126660 +/- 8.00e-51]
In [20]:
ivp.dop.numerical_solution(ivp.ini, [0,1], 1e-50)
Out[20]:
[7.28355749378782728852876426824018405440330223126660 +/- 8.00e-51]

Kreweras walks

Adapted from an exercise by B. Salvy & A. Bostan

In [21]:
@cached_function
def w(n, i, j):
    if i < 0 or j < 0:
        return 0
    elif n == 0:
        return 1 if i == j == 0 else 0
    else:
        return w(n-1,i-1,j-1)+w(n-1,i,j+1)+w(n-1,i+1,j)
In [22]:
Pol.<x,y> = ZZ[]
terms = [Pol({(i,j): w(k,i,j) for i in range(k+1) for j in range(k+1)})(y=0)
         for k in range(90)]
terms
Out[22]:
[1,
 0,
 x,
 2,
 2*x^2,
 8*x,
 5*x^3 + 16,
 30*x^2,
 14*x^4 + 96*x,
 112*x^3 + 192,
 42*x^5 + 480*x^2,
 420*x^4 + 1408*x,
 132*x^6 + 2240*x^3 + 2816,
 1584*x^5 + 8320*x^2,
 429*x^7 + 10080*x^4 + 23296*x,
 6006*x^6 + 44800*x^3 + 46592,
 1430*x^8 + 44352*x^5 + 153600*x^2,
 22880*x^7 + 228480*x^4 + 417792*x,
 4862*x^9 + 192192*x^6 + 913920*x^3 + 835584,
 87516*x^8 + 1123584*x^5 + 2976768*x^2,
 16796*x^10 + 823680*x^7 + 5107200*x^4 + 7938048*x,
 335920*x^9 + 5381376*x^6 + 19066880*x^3 + 15876096,
 58786*x^11 + 3500640*x^8 + 27320832*x^5 + 59924480*x^2,
 1293292*x^10 + 25259520*x^7 + 114250752*x^4 + 157515776*x,
 208012*x^12 + 14780480*x^9 + 141453312*x^6 + 406224896*x^3 + 315031552,
 4992288*x^11 + 116688000*x^8 + 652861440*x^5 + 1243545600*x^2,
 742900*x^13 + 62078016*x^10 + 713856000*x^7 + 2571878400*x^4 + 3233218560*x,
 19315400*x^12 + 532097280*x^9 + 3597834240*x^6 + 8817868800*x^3 + 6466437120,
 2674440*x^14 + 259598976*x^11 + 3528645120*x^8 + 15498362880*x^5 + 26453606400*x^2,
 74884320*x^13 + 2400349952*x^10 + 19262251008*x^7 + 58370457600*x^4 + 68191518720*x,
 9694845*x^15 + 1081662400*x^12 + 17145356800*x^9 + 89890504704*x^6 + 194568192000*x^3 + 136383037440,
 290845350*x^14 + 10730091008*x^11 + 100706411520*x^8 + 367464480768*x^5 + 574439424000*x^2,
 35357670*x^16 + 4493059200*x^13 + 82108522496*x^10 + 505506693120*x^7 + 1336234475520*x^4 + 1470564925440*x,
 1131445440*x^15 + 47593145600*x^12 + 516110712832*x^9 + 2224229449728*x^6 + 4355134586880*x^3 + 2941129850880,
 129644790*x^17 + 18614102400*x^14 + 388360068096*x^11 + 2770699616256*x^8 + 8725823225856*x^5 + 12692106510336*x^2,
 4407922860*x^16 + 209676096000*x^13 + 2600323934208*x^10 + 13038586429440*x^7 + 30848869990400*x^4 + 32307180208128*x,
 477638700*x^18 + 76938289920*x^15 + 1817192832000*x^12 + 14858993909760*x^9 + 54762063003648*x^6 + 98716383969280*x^3 + 64614360416256,
 17194993200*x^17 + 918295718400*x^14 + 12909337878528*x^11 + 74406691307520*x^8 + 207815008321536*x^5 + 284610873262080*x^2,
 1767263190*x^19 + 317370445920*x^16 + 8422988313600*x^13 + 78204394684416*x^10 + 332641678786560*x^7 + 717906392383488*x^4 + 721014212263936*x,
 67156001220*x^18 + 4000791075840*x^15 + 63266001111040*x^12 + 414962094243840*x^9 + 1345350789758976*x^6 + 2262492872966144*x^3 + 1442028424527872,
 6564120420*x^20 + 1306819483200*x^17 + 38717332992000*x^14 + 404902407110656*x^11 + 1965609920102400*x^8 + 4967449069879296*x^5 + 6464265351331840*x^2,
 262564816800*x^19 + 17349584376960*x^16 + 306534470860800*x^13 + 2268459448532992*x^10 + 8427713121484800*x^7 + 16831852220252160*x^4 + 16309838732591104*x,
 24466267020*x^21 + 5372480097600*x^18 + 176650313656320*x^15 + 2066269396172800*x^12 + 11342297242664960*x^9 + 33036635436220416*x^6 + 52365762463006720*x^3 + 32619677465182208,
 1027583214840*x^20 + 74924317036800*x^17 + 1470186481397760*x^14 + 12185085810573312*x^11 + 51338819098378240*x^8 + 119209118077550592*x^5 + 148465568301711360*x^2,
 91482563640*x^22 + 22055444611200*x^19 + 800619844907520*x^16 + 10409596236103680*x^13 + 64104147090407424*x^10 + 212603109677989888*x^7 + 397363726925168640*x^4 + 373284857444302848*x,
 4025232800160*x^21 + 322348805856000*x^18 + 6987227737374720*x^15 + 64440357652070400*x^12 + 305257843287654400*x^9 + 811757327861415936*x^6 + 1222657621308211200*x^3 + 746569714888605696,
 343059613650*x^23 + 90427322905920*x^20 + 3606821773632000*x^17 + 51840721922457600*x^14 + 355710774239428608*x^11 + 1330281548643041280*x^8 + 2872372083201933312*x^5 + 3443403096745574400*x^2,
 15780742227900*x^22 + 1382141195635200*x^19 + 32935435395793920*x^16 + 336070886945587200*x^13 + 1776835461659754496*x^10 + 5349581414863994880*x^7 + 9440663490244116480*x^4 + 8631463762508906496*x,
 1289904147324*x^24 + 370321417614720*x^21 + 16160420133580800*x^18 + 255499135191613440*x^15 + 1941742902352281600*x^12 + 8122676396158877696*x^9 + 19971770615492247552*x^6 + 28771545875029688320*x^3 + 17262927525017812992,
 61915399071552*x^23 + 5907918429853440*x^20 + 154091789814251520*x^17 + 1730800593233510400*x^14 + 10148842902961258496*x^11 + 34278471681397751808*x^8 + 69487580603014447104*x^5 + 80560328450083127296*x^2,
 4861946401452*x^25 + 1514951253878400*x^22 + 72047785729920000*x^19 + 1247409727067750400*x^16 + 10444486338478080000*x^13 + 48537944318510366720*x^10 + 134425379142736281600*x^7 + 225609027931865088000*x^4 + 201400821125207818240*x,
 243097320072600*x^24 + 25181856397800960*x^21 + 716044147408128000*x^18 + 8812868721102028800*x^15 + 56993369896386560000*x^12 + 214323390497318502400*x^9 + 492203695938019000320*x^6 + 681840617749636710400*x^3 + 402801642250415636480,
 18367353072152*x^26 + 6191539907155200*x^23 + 319751013795333120*x^20 + 6037993891657728000*x^17 + 55435787116609536000*x^14 + 284510902522761707520*x^11 + 879853918883728588800*x^8 + 1687555528930350858240*x^5 + 1899413149445416550400*x^2,
 955102359751904*x^25 + 107056555274073600*x^22 + 3306693410956615680*x^19 + 44410061194968268800*x^16 + 315198115253059584000*x^13 + 1311224159452727869440*x^10 + 3376091055264261734400*x^7 + 5420632911109611847680*x^4 + 4737359855087391866880*x,
 69533550916004*x^27 + 25282121287550400*x^24 + 1413146529617771520*x^21 + 28997157603773399040*x^18 + 290684036912519577600*x^15 + 1639030199315909836800*x^12 + 5619532111940262297600*x^9 + 12153927798951342243840*x^6 + 16261898733328835543040*x^3 + 9474719710174783733760,
 3754811749464216*x^26 + 454046259858048000*x^23 + 15183109225195591680*x^20 + 221677112183286988800*x^17 + 1719099143031029760000*x^14 + 7867344956716367216640*x^11 + 22523630732068257792000*x^8 + 41136371011835312209920*x^5 + 45096021697466518732800*x^2,
 263747951750360*x^28 + 103151054853205632*x^25 + 6221399815927296000*x^22 + 138252701887959859200*x^19 + 1507404362846351523840*x^16 + 9294991228526395392000*x^13 + 35118003574907842068480*x^10 + 84794845108962852864000*x^7 + 130888453219475993395200*x^4 + 112238987335916668846080*x,
 14769885298020160*x^27 + 1921441217853830400*x^24 + 69347869948202926080*x^21 + 1096906052341835366400*x^18 + 9257594470813956833280*x^15 + 46381067342343831552000*x^12 + 146646608334779999846400*x^9 + 300739050653121584824320*x^6 + 390098919399222568550400*x^3 + 224477974671833337692160,
 1002242216651368*x^29 + 420538915939992192*x^26 + 27292307910740121600*x^23 + 654772911603962511360*x^20 + 7737634585438352179200*x^17 + 51961981868439628677120*x^14 + 215208152468475378401280*x^11 + 575560372562068871577600*x^8 + 1006319131031599149219840*x^5 + 1077416063102614713139200*x^2,
 58130048565779344*x^28 + 8114549648452176384*x^25 + 315197121998930595840*x^22 + 5384195370680667340800*x^19 + 49275222153616934830080*x^16 + 269094652183894252584960*x^13 + 934248100013648967106560*x^10 + 2130701928229776528506880*x^7 + 3175017579190606941388800*x^4 + 2676528325181232340008960*x,
 3814986502092304*x^30 + 1713306694570338560*x^27 + 119331612477237888000*x^24 + 3081927415100654714880*x^21 + 39346043093435645952000*x^18 + 286692201621043984465920*x^15 + 1295640917922453808742400*x^12 + 3813257551076118233088000*x^9 + 7457456748804217849774080*x^6 + 9407459493898094641152000*x^3 + 5353056650362464680017920,
 228899190125538240*x^29 + 34203831829786031616*x^26 + 1426134536054500147200*x^23 + 26232219393647433154560*x^20 + 259471203102656692224000*x^17 + 1538553457379796162969600*x^14 + 5836364011810684233842688*x^11 + 14691076459935360771686400*x^8 + 24700716471152431973007360*x^5 + 25888196802757162696704000*x^2,
 14544636039226909*x^31 + 6975605827893521280*x^28 + 520156385939764933632*x^25 + 14423176229884519219200*x^22 + 198341171025139128729600*x^19 + 1562757988972572306309120*x^16 + 7675082764400362468147200*x^13 + 24722982956862649984352256*x^10 + 53579220030352492226150400*x^7 + 77345677838962160723558400*x^4 + 64202728070837763487825920*x,
 901767434432068358*x^30 + 143917762343908439040*x^27 + 6425461238079449180160*x^24 + 126923950822983769128960*x^21 + 1352788499812487390822400*x^18 + 8679118748674120577187840*x^15 + 35817052900535024851353600*x^12 + 98891931827450599937409024*x^9 + 185321066928513326052802560*x^6 + 227966208367467421079961600*x^3 + 128405456141675526975651840,
 55534064877048198*x^32 + 28383499575566741760*x^29 + 2260817212093398024192*x^26 + 67139513344830162862080*x^23 + 991777848291222009937920*x^20 + 8423850550183705266094080*x^17 + 44795451606059977172582400*x^14 + 157185695014919423347654656*x^11 + 374748373240865431341760512*x^8 + 608233245303838608583557120*x^5 + 625278742950767783533608960*x^2,
 3554180152131084672*x^31 + 604552505084105177600*x^28 + 28836083686323340836864*x^25 + 610173692708334914764800*x^22 + 6988136600425954541568000*x^19 + 48355090171184386072903680*x^16 + 216253904305117131177984000*x^13 + 651523315569086305469988864*x^10 + 1348575460797578022821560320*x^7 + 1891634494964091366408192000*x^4 + 1548309268259044035416555520*x,
 212336130412243110*x^33 + 115426231607304749824*x^30 + 9800114292942336563200*x^27 + 310977373087800734515200*x^24 + 4922067787847234979102720*x^21 + 44939093830431523051929600*x^18 + 257893814246316725722152960*x^15 + 981724073512119039950848000*x^12 + 2559555882592839057203527680*x^9 + 4615124910285044789211561984*x^6 + 5548794518561334674797363200*x^3 + 3096618536518088070833111040,
 14014184607208045260*x^32 + 2535592628750628930560*x^29 + 128933867315874086129664*x^26 + 2915752862537513592422400*x^23 + 35789919263571367522467840*x^20 + 266322099063884112042393600*x^17 + 1286269396116865058812723200*x^14 + 4209632827219966443309236224*x^11 + 9556731871228990411725864960*x^8 + 15022511853721360366502412288*x^5 + 15174254397698343804547891200*x^2,
 812944042149730764*x^34 + 469151780081303176704*x^31 + 42374480202510510602240*x^28 + 1433679732166448957718528*x^25 + 24256582324514208822067200*x^22 + 237435561943692974783201280*x^19 + 1466039745323095397719080960*x^16 + 6032159926617022344776908800*x^13 + 17113072580220298367365808128*x^10 + 33979491097703077019469742080*x^7 + 46433218456956932041916547072*x^4 + 37521792692490450134882058240*x,
 55280194866181691952*x^33 + 10619213307872036983808*x^30 + 574508742114037238480896*x^27 + 13854888167995094969548800*x^24 + 181834528240062069095792640*x^21 + 1451179588383130489234391040*x^18 + 7545491276627959529519185920*x^15 + 26720234341607254534641418240*x^12 + 66151372999170901252002283520*x^9 + 115172590878530429476518494208*x^6 + 135614479302858341201788010496*x^3 + 75043585384980900269764116480,
 3116285494907301262*x^35 + 1905929106580294155360*x^32 + 182789737266649816934400*x^29 + 6580736500578972004417536*x^26 + 118756184297100814024704000*x^23 + 1243240727966936007306117120*x^20 + 8236424690823173047006003200*x^17 + 36510441661103029981544448000*x^14 + 112224984234750469045493956608*x^11 + 243715584733787530928429465600*x^8 + 372096062838329079847213596672*x^5 + 369857670825977294186694574080*x^2,
 218139984643511088340*x^34 + 44413035181030034061312*x^31 + 2551620806996826936053760*x^28 + 65488083989589123019702272*x^25 + 916921271665085717741568000*x^22 + 7828832965468067097005260800*x^19 + 43698305942081625792247234560*x^16 + 166856914993730628973127270400*x^13 + 448325896098926426212433657856*x^10 + 857154502345457074720420659200*x^7 + 1143671881451141327668924907520*x^4 + 913387639257196100426271817728*x,
 11959798385860453492*x^36 + 7739227281265436873280*x^33 + 786745194635389174800384*x^30 + 30082266356173117561896960*x^27 + 577836035202256967820902400*x^24 + 6455125752522203452900638720*x^21 + 45768561951967161490184601600*x^18 + 217924019243627848106791403520*x^15 + 723046631639499392216884838400*x^12 + 1707908175614957814142604410880*x^9 + 2880039127880735771060613414912*x^6 + 3327045473312411135036872458240*x^3 + 1826775278514392200852543635456,
 861105483781952651424*x^35 + 185510433040481964455040*x^32 + 11298176893452474051231744*x^29 + 308011153184505011503890432*x^26 + 4591241422559429513025945600*x^23 + 41842316188653437180958474240*x^20 + 250061914365425572549283020800*x^17 + 1026351832566763413664243384320*x^14 + 2980653408546830435679958204416*x^11 + 6216596517335331212696349573120*x^8 + 9241444234518404891754935353344*x^5 + 9051216877707056379789752401920*x^2,
 45950804324621742364*x^37 + 31412157788665596720960*x^34 + 3379143887998933014073344*x^31 + 136982415894627171265216512*x^28 + 2795346503429186991195684864*x^25 + 33252310217855953324128337920*x^22 + 251734259996776776535847731200*x^19 + 1283991380700838245987747102720*x^16 + 4583191804048133174811017871360*x^13 + 11721023548585120650451623084032*x^10 + 21648383166250094575977876160512*x^7 + 28258961708857932313630794055680*x^4 + 22326334965010739070148055924736*x,
 3400359520022008934936*x^36 + 773922728126543687328000*x^33 + 49882600251412820683939840*x^30 + 1441920167311864960686489600*x^27 + 22837798230630612673167360000*x^24 + 221682068119039688827522252800*x^21 + 1415074242585431584077250560000*x^18 + 6225412754913155132061804134400*x^15 + 19471076291707755317824585728000*x^12 + 44063998302951581392675274752000*x^9 + 72161277220833648586592920535040*x^6 + 81910033938718644387335634944000*x^3 + 44652669930021478140296111849472,
 176733862787006701400*x^38 + 127443611599728992410752*x^35 + 14485061210010235580736000*x^32 + 621488134279897438029414400*x^29 + 13449182651472485906039439360*x^26 + 170025077439470357207580672000*x^23 + 1371335584178245516933044633600*x^20 + 7474214532265908907404165120000*x^17 + 28616816695971761494155067392000*x^14 + 78922762569055434888248987484160*x^11 + 158630393890625693013630989107200*x^8 + 230108688200700305982422180167680*x^5 + 222327234976522034765625294848000*x^2,
 13431773571812509306400*x^37 + 3224981532969667930018560*x^34 + 219638097362862895205744640*x^31 + 6720498807975840092589260800*x^28 + 112894187288251223412540702720*x^25 + 1164852658202328830230659072000*x^22 + 7924415758478416870832603136000*x^19 + 37271416467565999084922103398400*x^16 + 125147782548712206453140619264000*x^13 + 305938580759679420367720331345920*x^10 + 547430378916669058243118707507200*x^7 + 700330790176044409511719678771200*x^4 + 547814306982150293662500726505472*x,
 680425371729975800390*x^39 + 516854647043345358110272*x^36 + 61975732068373618481226240*x^33 + 2809973118642341167235399680*x^30 + 64375304371136994571118182400*x^27 + 863308491027803473154723020800*x^24 + 7403285783241467676577077657600*x^21 + 43018256974597120155948417024000*x^18 + 176192150573948359310540852428800*x^15 + 522221611129194145446438633472000*x^12 + 1136343299964523561365818373570560*x^9 + 1811496890233341247277229177569280*x^6 + 2023177838286350516367190183116800*x^3 + 1095628613964300587325001453010944,
 53073178994938112430420*x^38 + 13424060421838120533932544*x^35 + 964601394043283721674966016*x^32 + 31192675836919665416149401600*x^29 + 554798077671253371394727608320*x^26 + 6073591461924732412249925222400*x^23 + 43942937797486922631275033395200*x^20 + 220439500604962540042373509939200*x^17 + 792363181898022033332033814528000*x^14 + 2084488788844109683718837155921920*x^11 + 4049824993858527579604495556935680*x^8 + 5743475090104372917788673111490560*x^5 + 5479921687701315112903132267413504*x^2,
 2622127042276492108820*x^40 + 2095356677202751451798400*x^37 + 264699782965822095035289600*x^34 + 12663484967952852961476476928*x^31 + 306639864159549253243502592000*x^28 + 4354641515306441556607673303040*x^25 + 39629107411140098009006604288000*x^22 + 244978050438951485052404367360000*x^19 + 1070706145795532337348671333990400*x^16 + 3400179171363159683263899893760000*x^13 + 7975435366012245746402507379179520*x^10 + 13860363615344693320571535595929600*x^7 + 17404469970013251266026282156032000*x^4 + 13489038000495544893300017889017856*x,
 209770163382119368705600*x^39 + 55820301880681298675909376*x^36 + 4225989578480429273780797440*x^33 + 144208665553422284744977022976*x^30 + 2711341956779172344468865024000*x^27 + 31437251581089283964814753792000*x^24 + 241432715920176597101024850739200*x^21 + 1288961434617252429044958363648000*x^18 + 4947005935654224275236641885388800*x^15 + 13958630282438234489188641669120000*x^12 + 29297517671065392537805129148006400*x^9 + 45558412579132991958052525697925120*x^6 + 50124873513638163646155692609372160*x^3 + 26978076000991089786600035778035712,
 10113918591637898134020*x^41 + 8491708639190097988867200*x^38 + 1128640624327199956789619712*x^35 + 56893173131184286641048944640*x^32 + 1453907037956634510133784739840*x^29 + 21828767608396682075105698775040*x^26 + 210437112624434390621617127424000*x^23 + 1381219723636359136903537518182400*x^20 + 6427388775321164139426886975488000*x^17 + 21809381006647655406957238419456000*x^14 + 54941168791676890949446493609656320*x^11 + 103452478666250046377512848379084800*x^8 + 143684224288034820790781042585763840*x^5 + 135502405322582288537959344856104960*x^2,
 829341324514307646989640*x^40 + 231886138943771160665689600*x^37 + 18471555006594455630838564864*x^34 + 664212166313958089312641744896*x^31 + 13180995255599319268896496680960*x^28 + 161600276154099622360473749422080*x^25 + 1314973643208004097206766272512000*x^22 + 7456340621906849324422348878643200*x^19 + 30484186762951807061281806798028800*x^16 + 91987298474499366907021092716544000*x^13 + 207707192988824022512814280188887040*x^10 + 351368009383069211225925991425638400*x^7 + 433665113305705095477630055804305408*x^4 + 333235544941461628108166685127606272*x,
 39044429911904443959240*x^42 + 34402306794667576467718400*x^39 + 4804680798914938448993088512*x^36 + 254853918351853938558816141312*x^33 + 6863525718577566922897298030592*x^30 + 108777897688314382176998456819712*x^27 + 1109021503018330741689525731328000*x^24 + 7714512040153624036946362132070400*x^21 + 38142050104369652313391246186905600*x^18 + 137948643129317268317719691368857600*x^15 + 372037518274641883935063086098022400*x^12 + 755298883595723718228415564323225600*x^9 + 1147802163984692756671358238657085440*x^6 + 1245397248467665915217809391027748864*x^3 + 666471089882923256216333370255212544,
 3279732112599973292576160*x^41 + 962393645774877772071616000*x^38 + 80560675313313488514623840256*x^35 + 3048465846383583358880508641280*x^32 + 63759528533234228245494025420800*x^29 + 825273637007211593706484160004096*x^26 + 7103308841938444938921923051520000*x^23 + 42698927105966570251005446219366400*x^20 + 185555919426663173416497954422784000*x^17 + 597230967707570516485883213905920000*x^14 + 1445631499581465606147673705980887040*x^11 + 2644410278950703178236786758385664000*x^8 + 3602332945428881882476262779785314304*x^5 + 3360595749833384215667104705947893760*x^2,
 150853479205085351660700*x^43 + 139329342518403684694259520*x^40 + 20422743080469484150195072000*x^37 + 1138439590108326762577548165120*x^34 + 32266838497106236167842922233856*x^31 + 539038183599614119472074913218560*x^28 + 5802867460214229948200310005563392*x^25 + 42706215469161836137348218224640000*x^22 + 223909007994702746438199291150336000*x^19 + 861215092704068442968000029733683200*x^16 + 2479538224551430696031046170836992000*x^13 + 5405404737565480092552171248450273280*x^10 + 8918403293716096993269163185143808000*x^7 + 10832189975765169296956594372781015040*x^4 + 8257463842447744072782028706043396096*x,
 12973399211637340242820200*x^42 + 3990667588181438870255334400*x^39 + 350617653205500103890548996096*x^36 + 13944117215488325315546117898240*x^33 + 306961775755539749363923672891392*x^30 + 4188507225673556602979184971612160*x^27 + 38073112295523228166925110896230400*x^24 + 242191693060668901827805895353958400*x^21 + 1116505551177205550112649859039232000*x^18 + 3823959837461126866384229797095014400*x^15 + 9891914398475019813795813824397312000*x^12 + 19472886632223468656399126236653158400*x^9 + 28966973897989883034138242025347088384*x^6 + 31025531782438509591283085117101178880*x^3 + 16514927684895488145564057412086792192,
 583300119592996693088040*x^44 + 564113923367195406323099520*x^41 + 86683361788852520270356377600*x^38 + 5071948517602850817923558080512*x^35 + 151095956693201853419201217822720*x^32 + 2656980616375819469904126217814016*x^29 + 30157252024849607541450131795607552*x^26 + 234432691452317953027830420445593600*x^23 + 1301076304581732937726120042482892800*x^20 + 5310945324518599373508820951105536000*x^17 + 16282667694995766011700591394081996800*x^14 + 37984951290144076084975925085685678080*x^11 + 67642658827723627964333806927321497600*x^8 + 90500368154666616816715927747829956608*x^5 + 83578983577181291143864637458313379840*x^2,
 51330410524183708991747520*x^43 + 16533748645517237250385463040*x^40 + 1522925348415060901581014384640*x^37 + 63577946206570946872562911150080*x^34 + 1471183879187141465086786558492672*x^31 + 21133026104363605706695068623044608*x^28 + 202565692846159628014268809796911104*x^25 + 1361372253625233616589982611779092480*x^22 + 6645382559987043412202277404934144000*x^19 + 24166797822546288277049161200068198400*x^16 + 66627927579522904599602419957392998400*x^13 + 140594583414181665358001628199948124160*x^10 + 226642696872184579344026684857660735488*x^7 + 271196389419708043659519110086089768960*x^4 + 205200814851562342394591799552824573952*x]

Guessing differential operators

In [23]:
from ore_algebra.guessing import guess
%time dop = guess(terms, OreAlgebra(ZZ['x']['t'], 'Dt'))
dop

CPU times: user 564 ms, sys: 4.14 ms, total: 568 ms
Wall time: 566 ms
Out[23]:
((31104*x^6 - 3888*x^3 - 972)*t^12 + (7776*x^4 + 3888*x)*t^11 + (-6912*x^8 - 8208*x^5 - 4320*x^2)*t^10 + (3168*x^6 - 936*x^3 + 36)*t^9 + (4572*x^4 - 144*x)*t^8 + (256*x^8 - 2504*x^5 + 160*x^2)*t^7 + (272*x^6 + 40*x^3)*t^6 - 180*x^4*t^5 + 104*x^5*t^4 - 16*x^6*t^3)*Dt^4 + ((497664*x^6 - 62208*x^3 - 15552)*t^11 + (77760*x^4 + 59292*x)*t^10 + (-124416*x^8 - 93312*x^5 - 57348*x^2)*t^9 + (80928*x^6 - 28692*x^3 + 414)*t^8 + (-12096*x^7 + 78228*x^4 - 1548*x)*t^7 + (3456*x^8 - 40248*x^5 + 1404*x^2)*t^6 + (3712*x^6 + 944*x^3)*t^5 + (448*x^7 - 2146*x^4)*t^4 + 1100*x^5*t^3 - 168*x^6*t^2)*Dt^3 + ((2239488*x^6 - 279936*x^3 - 69984)*t^10 + (139968*x^4 + 253692*x)*t^9 + (-635904*x^8 - 262656*x^5 - 205632*x^2)*t^8 + (525888*x^6 - 186696*x^3 + 1170)*t^7 + (-117504*x^7 + 362700*x^4 - 4005*x)*t^6 + (12672*x^8 - 171528*x^5 + 2400*x^2)*t^5 + (13008*x^6 + 4497*x^3)*t^4 + (3360*x^7 - 6486*x^4)*t^3 + 2856*x^5*t^2 - 432*x^6*t)*Dt^2 + ((2985984*x^6 - 373248*x^3 - 93312)*t^9 + (-93312*x^4 + 320760*x)*t^8 + (-967680*x^8 - 158976*x^5 - 204984*x^2)*t^7 + (953280*x^6 - 320904*x^3 + 720)*t^6 + (-253440*x^7 + 478080*x^4 - 2160*x)*t^5 + (13056*x^8 - 198744*x^5 + 135*x^2)*t^4 + (10656*x^6 + 4473*x^3)*t^3 + (5376*x^7 - 4506*x^4)*t^2 + 1608*x^5*t - 240*x^6)*Dt + (746496*x^6 - 93312*x^3 - 23328)*t^8 + (-93312*x^4 + 75816*x)*t^7 + (-276480*x^8 + 3456*x^5 - 34128*x^2)*t^6 + (311040*x^6 - 97200*x^3)*t^5 + (-92160*x^7 + 112680*x^4)*t^4 + (2304*x^8 - 35136*x^5)*t^3 - 288*x^6*t^2 + 1344*x^7*t

Guessing algebraic equations

In [24]:
%time algeq = guess(terms, OreAlgebra(ZZ['x']['t'], 'C'))
algeq

CPU times: user 950 ms, sys: 189 ms, total: 1.14 s
Wall time: 1.13 s
Out[24]:
16*x^6*t^10*C^6 + (96*x^4*t^9 - 48*x^5*t^8)*C^5 + ((48*x^5 + 192*x^2)*t^8 - 192*x^3*t^7 + 56*x^4*t^6)*C^4 + ((192*x^3 + 128)*t^7 + (-96*x^4 - 192*x)*t^6 + 128*x^2*t^5 - 32*x^3*t^4)*C^3 + ((48*x^4 + 192*x)*t^6 - 264*x^2*t^5 + (64*x^3 + 32)*t^4 - 32*x*t^3 + 9*x^2*t^2)*C^2 + (96*x^2*t^5 + (-48*x^3 - 144)*t^4 + 104*x*t^3 - 16*x^2*t^2 + 2*t - x)*C + (16*x^3 + 108)*t^4 - 72*x*t^3 + 8*x^2*t^2 - 2*t + x

The Voigt profile

$$V(x) = \frac{1}{σ \sqrt{2π}} \frac{λ}{π} \int_{-∞}^{+∞} \frac{\exp\bigl(-(u-x)²/(2σ²)\bigr)}{u² + λ²} \mathrm du, \qquad σ = 1, λ = 1/2$$
In [25]:
P.<x, u> = PolynomialRing(QQ)
A.<Dx, Du> = OreAlgebra(P)
f = exp(-(u-x)^2/2)/(u^2 + 1/4)
In [26]:
op1 = Dx + x - u
op1(f)
Out[26]:
0
In [27]:
op2 = (4*u^2 + 1)*Du + (4*u^3 - 4*u^2*x + 9*u - x)
op2(f).simplify_full()
Out[27]:
0

Creative telescoping

$$V(x) = \frac{1}{(2π)^{3/2}} \int_{-∞}^{+∞} f(x,u) \, \mathrm du, \qquad f(x,u) = \frac{\exp\bigl(-(u-x)²/2\bigr)}{u² + 1/4}$$
In [28]:
ideal = A.ideal([op1, op2]); ideal
Out[28]:
Left Ideal (Dx + x - u, (4*u^2 + 1)*Du - 4*x*u^2 + 4*u^3 - x + 9*u) of Multivariate Ore algebra in Dx, Du over Fraction Field of Multivariate Polynomial Ring in x, u over Rational Field

Idea: Find a pair $(P,Q)$ such that $P - D_u Q \in I$ and $P$ does not depend on $u$.

Then $P\left(\int f \, \mathrm du \right) = \int P(f) \, \mathrm du = \int (D_u Q(f)) \, \mathrm du = 0$.

In [29]:
[tel], [cert] = ideal.ct(Du)
tel, cert
Out[29]:
(-4*Dx^3 - 8*x*Dx^2 + (-4*x^2 - 13)*Dx - 8*x, 4*u^2 + 1)
In [30]:
(tel-Du*cert).reduce(ideal)
Out[30]:
0
In [31]:
ini = [(1/2*exp(1/8)*erfc(1/4*2^(1/2))*2^(1/2)/pi^(1/2)), 0,
       1/2*(1/2/pi-5/8*exp(1/8)*erfc(1/4*2^(1/2))*2^(1/2)/pi^(1/2))]
tel.numerical_solution(ini, [0,2])
Out[31]:
[0.08242408278858694 +/- 2.82e-18]

Another walk model

After A. Bostan, F. Chyzak, M. van Hoeij, M. Kauers and L. Pech, Hypergeometric expressions for generating functions of walks with small steps in the quarter plane, 2017.

In [32]:
from ore_algebra.examples import ssw
Rat = Frac(PolynomialRing(ZZ,'u,v,t'))
q = Rat(ssw.rat[10](x=1,y=1))
q
Out[32]:
(-u^3*v^2 - u^2*v^2 + u^3 - u*v^2 + 2*u^2 - v^2 + 2*u + 1)/(u^4*v^3*t + u^3*v^3*t + 2*u^2*v^3*t - u^3*v^2 - u^4*t + u^3*v*t + u*v^3*t + u^3*v - u^2*v^2 - 2*u^3*t + u^2*v*t + v^3*t + u^2*v - u*v^2 - 3*u^2*t + u*v*t + u*v - 2*u*t - t)
In [33]:
A.<Du,Dv,Dt> = OreAlgebra(Rat)
In [34]:
ideal = A.ideal([q*D - D(q) for D in Du,Dv,Dt])
ideal
Out[34]:
Left Ideal ((-u^7*v^4*t - u^7*v^3*t - 2*u^6*v^4*t - 2*u^6*v^3*t - 4*u^5*v^4*t + u^6*v^3 + u^7*v*t - 4*u^5*v^3*t - 5*u^4*v^4*t + 2*u^5*v^3 + u^7*t + 3*u^6*v*t - 5*u^4*v^3*t - 5*u^3*v^4*t - u^6*v + 3*u^4*v^3 + 4*u^6*t + 6*u^5*v*t - 5*u^3*v^3*t - 4*u^2*v^4*t - 3*u^5*v + 3*u^3*v^3 + 9*u^5*t + 8*u^4*v*t - 4*u^2*v^3*t - 2*u*v^4*t - 5*u^4*v + 2*u^2*v^3 + 13*u^4*t + 8*u^3*v*t - 2*u*v^3*t - v^4*t - 5*u^3*v + u*v^3 + 13*u^3*t + 6*u^2*v*t - v^3*t - 3*u^2*v + 9*u^2*t + 3*u*v*t - u*v + 4*u*t + v*t + t)*Du - u^6*v^4*t - u^6*v^3*t - 2*u^5*v^4*t - 2*u^5*v^3*t - 2*u^4*v^4*t + u^6*v*t + 2*u^5*v^2*t - 2*u^4*v^3*t - 4*u^3*v^4*t + u^6*t + 4*u^5*v*t + 4*u^4*v^2*t - 4*u^3*v^3*t - u^2*v^4*t + 4*u^5*t + 6*u^4*v*t + 2*u^3*v^2*t - u^2*v^3*t - 2*u*v^4*t - u^4*v + 3*u^2*v^3 + 7*u^4*t + 6*u^3*v*t - 2*u^2*v^2*t - 2*u*v^3*t - 2*u^3*v + 2*u*v^3 + 8*u^3*t + 2*u^2*v*t - 4*u*v^2*t - 3*u^2*v + v^3 + 5*u^2*t - 2*v^2*t - 2*u*v + 2*u*t - v*t - v, (-u^5*v^5*t - u^4*v^5*t + u^5*v^3*t - 2*u^3*v^5*t + u^4*v^4 + u^5*v^2*t + u^4*v^3*t - 2*u^2*v^5*t - u^4*v^3 + u^3*v^4 + 2*u^4*v^2*t + 2*u^3*v^3*t - u*v^5*t - u^4*v^2 - u^3*v^3 + u^2*v^4 - u^5*t + u^4*v*t + 3*u^3*v^2*t + 2*u^2*v^3*t - v^5*t + u^4*v - 2*u^3*v^2 - u^2*v^3 + u*v^4 - 3*u^4*t + 2*u^3*v*t + 3*u^2*v^2*t + u*v^3*t + 2*u^3*v - 2*u^2*v^2 - u*v^3 - 5*u^3*t + 2*u^2*v*t + 2*u*v^2*t + v^3*t + 2*u^2*v - u*v^2 - 5*u^2*t + u*v*t + v^2*t + u*v - 3*u*t - t)*Dv - u^5*v^4*t - u^4*v^4*t + 3*u^5*v^2*t - 2*u^3*v^4*t - 2*u^5*v*t + 7*u^4*v^2*t - 2*u^2*v^4*t + u^4*v^2 - 4*u^4*v*t + 10*u^3*v^2*t - u*v^4*t - 2*u^4*v + u^3*v^2 + u^4*t - 6*u^3*v*t + 10*u^2*v^2*t - v^4*t + u^4 - 4*u^3*v + u^2*v^2 + 2*u^3*t - 6*u^2*v*t + 7*u*v^2*t + 2*u^3 - 4*u^2*v + u*v^2 + 2*u^2*t - 4*u*v*t + 3*v^2*t + 2*u^2 - 2*u*v + u*t - 2*v*t + u, (-u^5*v^4*t - u^5*v^3*t - u^4*v^4*t - u^4*v^3*t - 2*u^3*v^4*t + u^4*v^3 + u^5*v*t - 2*u^3*v^3*t - 2*u^2*v^4*t + u^3*v^3 + u^5*t + 2*u^4*v*t - 2*u^2*v^3*t - u*v^4*t - u^4*v + u^2*v^3 + 3*u^4*t + 3*u^3*v*t - u*v^3*t - v^4*t - 2*u^3*v + u*v^3 + 5*u^3*t + 3*u^2*v*t - v^3*t - 2*u^2*v + 5*u^2*t + 2*u*v*t - u*v + 3*u*t + v*t + t)*Dt - u^5*v^4 - u^5*v^3 - u^4*v^4 - u^4*v^3 - 2*u^3*v^4 + u^5*v - 2*u^3*v^3 - 2*u^2*v^4 + u^5 + 2*u^4*v - 2*u^2*v^3 - u*v^4 + 3*u^4 + 3*u^3*v - u*v^3 - v^4 + 5*u^3 + 3*u^2*v - v^3 + 5*u^2 + 2*u*v + 3*u + v + 1) of Multivariate Ore algebra in Du, Dv, Dt over Fraction Field of Multivariate Polynomial Ring in u, v, t over Integer Ring

⚐ Fast creative telescoping

In [35]:
%time tel1, _ = ideal.ct(Dv)
tel1

CPU times: user 1.87 s, sys: 52.6 ms, total: 1.93 s
Wall time: 1.96 s
Out[35]:
[(-9*u^9*t^4 - 24*u^8*t^4 - 3*u^8*t^3 - 37*u^7*t^4 - u^7*t^3 - 37*u^6*t^4 + 5*u^7*t^2 - u^6*t^3 - 15*u^5*t^4 + 6*u^6*t^2 - 3*u^5*t^3 + 15*u^4*t^4 - u^6*t + u^5*t^2 + 3*u^4*t^3 + 37*u^3*t^4 - u^5*t - u^4*t^2 + u^3*t^3 + 37*u^2*t^4 + u^4*t - 6*u^3*t^2 + u^2*t^3 + 24*u*t^4 + u^3*t - 5*u^2*t^2 + 3*u*t^3 + 9*t^4)*Dt^2 + (-2*u^8 - 4*u^7 - 6*u^6 - 6*u^5 - 4*u^4 - 2*u^3)*Du + (-36*u^9*t^3 - 96*u^8*t^3 - 12*u^8*t^2 - 148*u^7*t^3 - 6*u^7*t^2 - 148*u^6*t^3 + 16*u^7*t - 6*u^6*t^2 - 60*u^5*t^3 + 20*u^6*t - 12*u^5*t^2 + 60*u^4*t^3 - 2*u^6 + 4*u^5*t + 12*u^4*t^2 + 148*u^3*t^3 - 2*u^5 - 4*u^4*t + 6*u^3*t^2 + 148*u^2*t^3 + 2*u^4 - 20*u^3*t + 6*u^2*t^2 + 96*u*t^3 + 2*u^3 - 16*u^2*t + 12*u*t^2 + 36*t^3)*Dt - 18*u^9*t^2 - 48*u^8*t^2 - 6*u^8*t - 74*u^7*t^2 - 4*u^7*t - 74*u^6*t^2 + 4*u^7 - 4*u^6*t - 30*u^5*t^2 + 4*u^6 - 6*u^5*t + 30*u^4*t^2 + 6*u^4*t + 74*u^3*t^2 - 6*u^4 + 4*u^3*t + 74*u^2*t^2 - 8*u^3 + 4*u^2*t + 48*u*t^2 - 6*u^2 + 6*u*t + 18*t^2,
 (-9*u^10*t^3 - 33*u^9*t^3 - 3*u^9*t^2 - 79*u^8*t^3 - 4*u^8*t^2 - 131*u^7*t^3 + 5*u^8*t - 8*u^7*t^2 - 168*u^6*t^3 + 11*u^7*t - 9*u^6*t^2 - 168*u^5*t^3 - u^7 + 17*u^6*t - 9*u^5*t^2 - 131*u^4*t^3 - 2*u^6 + 17*u^5*t - 8*u^4*t^2 - 79*u^3*t^3 - 2*u^5 + 11*u^4*t - 4*u^3*t^2 - 33*u^2*t^3 - u^4 + 5*u^3*t - 3*u^2*t^2 - 9*u*t^3)*Du*Dt + (-9*u^10*t^2 - 33*u^9*t^2 - 6*u^9*t - 79*u^8*t^2 - 15*u^8*t - 131*u^7*t^2 + u^8 - 30*u^7*t - 168*u^6*t^2 + 2*u^7 - 39*u^6*t - 168*u^5*t^2 + 3*u^6 - 39*u^5*t - 131*u^4*t^2 + 3*u^5 - 30*u^4*t - 79*u^3*t^2 + 2*u^4 - 15*u^3*t - 33*u^2*t^2 + u^3 - 6*u^2*t - 9*u*t^2)*Du + (-18*u^9*t^3 - 60*u^8*t^3 - 3*u^8*t^2 - 106*u^7*t^3 - u^7*t^2 - 137*u^6*t^3 + 8*u^7*t + 4*u^6*t^2 - 109*u^5*t^3 + 17*u^6*t - 7*u^5*t^2 - 59*u^4*t^3 - u^6 + 10*u^5*t - 2*u^4*t^2 + 6*u^3*t^3 - 2*u^5 + 7*u^4*t - 12*u^3*t^2 + 27*u^2*t^3 - 6*u^3*t - 3*u^2*t^2 + 27*u*t^3 - 3*u^2*t + 9*t^3)*Dt - 18*u^9*t^2 - 60*u^8*t^2 - 6*u^8*t - 106*u^7*t^2 - 12*u^7*t - 137*u^6*t^2 + 4*u^7 - 12*u^6*t - 109*u^5*t^2 + 8*u^6 - 24*u^5*t - 59*u^4*t^2 + 4*u^5 - 15*u^4*t + 6*u^3*t^2 - u^4 - 18*u^3*t + 27*u^2*t^2 - 6*u^3 - 3*u^2*t + 27*u*t^2 - 3*u^2 + 9*t^2,
 (9*u^12*t^3 + 33*u^11*t^3 + 3*u^11*t^2 + 70*u^10*t^3 + 4*u^10*t^2 + 98*u^9*t^3 - 5*u^10*t + 5*u^9*t^2 + 89*u^8*t^3 - 11*u^9*t + 5*u^8*t^2 + 37*u^7*t^3 + u^9 - 12*u^8*t + u^7*t^2 - 37*u^6*t^3 + 2*u^8 - 6*u^7*t - u^6*t^2 - 89*u^5*t^3 + u^7 + 6*u^6*t - 5*u^5*t^2 - 98*u^4*t^3 - u^6 + 12*u^5*t - 5*u^4*t^2 - 70*u^3*t^3 - 2*u^5 + 11*u^4*t - 4*u^3*t^2 - 33*u^2*t^3 - u^4 + 5*u^3*t - 3*u^2*t^2 - 9*u*t^3)*Du^2 + (36*u^11*t^3 + 126*u^10*t^3 + 18*u^10*t^2 + 172*u^9*t^3 + 34*u^9*t^2 + 98*u^8*t^3 - 14*u^9*t + 8*u^8*t^2 - 158*u^7*t^3 - 30*u^8*t - 8*u^7*t^2 - 442*u^6*t^3 + 2*u^8 + 6*u^7*t - 64*u^6*t^2 - 590*u^5*t^3 + 4*u^7 + 36*u^6*t - 68*u^5*t^2 - 514*u^4*t^3 - 4*u^6 + 60*u^5*t - 28*u^4*t^2 - 294*u^3*t^3 - 8*u^5 + 54*u^4*t - 12*u^3*t^2 - 108*u^2*t^3 - 4*u^4 + 14*u^3*t + 18*u^2*t^2 - 6*u*t^3 - 2*u^3 + 6*u^2*t + 6*u*t^2)*Du + (-6*u^9*t^5 - 6*u^9*t^4 - 12*u^8*t^5 - 6*u^9*t^3 - 12*u^8*t^4 + 6*u^7*t^5 - 12*u^8*t^3 + 6*u^7*t^4 + 30*u^6*t^5 + 6*u^7*t^3 + 30*u^6*t^4 + 18*u^5*t^5 + 30*u^6*t^3 + 18*u^5*t^4 - 18*u^4*t^5 + 18*u^5*t^3 - 18*u^4*t^4 - 30*u^3*t^5 - 18*u^4*t^3 - 30*u^3*t^4 - 6*u^2*t^5 - 30*u^3*t^3 - 6*u^2*t^4 + 12*u*t^5 - 6*u^2*t^3 + 12*u*t^4 + 6*t^5 + 12*u*t^3 + 6*t^4 + 6*t^3)*Dt + 18*u^10*t^3 - 6*u^9*t^4 + 54*u^9*t^3 - 12*u^8*t^4 + 6*u^9*t^2 + 38*u^8*t^3 + 6*u^7*t^4 + 14*u^8*t^2 - 68*u^7*t^3 + 30*u^6*t^4 - 4*u^8*t - 2*u^7*t^2 - 148*u^6*t^3 + 18*u^5*t^4 - 8*u^7*t - 16*u^6*t^2 - 272*u^5*t^3 - 18*u^4*t^4 + 8*u^6*t + 10*u^5*t^2 - 244*u^4*t^3 - 30*u^3*t^4 + 46*u^5*t - 76*u^4*t^2 - 188*u^3*t^3 - 6*u^2*t^4 + 2*u^4*t - 2*u^3*t^2 - 30*u^2*t^3 + 12*u*t^4 - 6*u^4 + 22*u^3*t - 6*u*t^3 + 6*t^4 + 12*u*t^2 + 6*t^3 + 6*t^2]
In [36]:
%time tel2, _ = tel1[0].parent().ideal(tel1).ct(Du)
tel2   # HolonomicFunctions.m: ≥ 1 min

CPU times: user 12.9 s, sys: 90.3 ms, total: 13 s
Wall time: 13.1 s
Out[36]:
[(1892352000*t^24 + 9122611200*t^23 - 18640527360*t^22 - 58034765824*t^21 - 38532706304*t^20 + 41965621248*t^19 + 94071750144*t^18 + 56013030400*t^17 - 5294450432*t^16 - 21996570880*t^15 - 10952369952*t^14 - 2543280896*t^13 - 537441232*t^12 - 276140992*t^11 - 112877022*t^10 - 19917936*t^9 + 1121240*t^8 + 1056920*t^7 + 134035*t^6 - 5391*t^5 - 2254*t^4 - 103*t^3 + 3*t^2)*Dt^5 + (47308800000*t^23 + 227888332800*t^22 - 542727905280*t^21 - 1484019662848*t^20 - 767620100096*t^19 + 1085290090496*t^18 + 1896743070208*t^17 + 918495748096*t^16 - 215512785664*t^15 - 427218085376*t^14 - 200103936864*t^13 - 53308965120*t^12 - 16198105488*t^11 - 7684582384*t^10 - 2788409498*t^9 - 526917856*t^8 - 11674372*t^7 + 14725960*t^6 + 2406665*t^5 + 42072*t^4 - 17460*t^3 - 836*t^2 + 27*t)*Dt^4 + (378470400000*t^22 + 1821691084800*t^21 - 4964213882880*t^20 - 12009478135808*t^19 - 4381901037568*t^18 + 8936076752896*t^17 + 11847464069120*t^16 + 4115996062720*t^15 - 2219866571776*t^14 - 2632683934976*t^13 - 1121810512000*t^12 - 331856938400*t^11 - 129562335520*t^10 - 61796701736*t^9 - 21296426392*t^8 - 4208335952*t^7 - 300047680*t^6 + 45411880*t^5 + 10529090*t^4 + 535888*t^3 - 18106*t^2 - 1020*t + 48)*Dt^3 + (1135411200000*t^21 + 5460826521600*t^20 - 16785832181760*t^19 - 36057594863616*t^18 - 7287121182720*t^17 + 27418579642368*t^16 + 26356460402688*t^15 + 4411389139968*t^14 - 7504283639808*t^13 - 5923884943104*t^12 - 2207397164160*t^11 - 718441676064*t^10 - 351342324480*t^9 - 173136422904*t^8 - 57960776808*t^7 - 11637201048*t^6 - 1161271776*t^5 - 2986416*t^4 + 10037526*t^3 + 711192*t^2 + 15354*t + 792)*Dt^2 + (1135411200000*t^20 + 5456579788800*t^19 - 18705041326080*t^18 - 35704801148928*t^17 - 845976502272*t^16 + 27816503623680*t^15 + 17768620793856*t^14 - 1927605172224*t^13 - 7625766414336*t^12 - 4160593317888*t^11 - 1282076116992*t^10 - 484808936256*t^9 - 302844173952*t^8 - 153142047744*t^7 - 49288460160*t^6 - 9650682192*t^5 - 1077476208*t^4 - 60418992*t^3 - 1875024*t^2 - 158760*t - 9216)*Dt + 227082240000*t^19 + 1090466611200*t^18 - 4130053816320*t^17 - 6994180227072*t^16 + 1235817283584*t^15 + 5584717234176*t^14 + 1907260735488*t^13 - 1376741382144*t^12 - 1399425761280*t^11 - 495281498112*t^10 - 109934770176*t^9 - 64845759168*t^8 - 53271954240*t^7 - 26356354176*t^6 - 7768879584*t^5 - 1370419584*t^4 - 140485008*t^3 - 8567520*t^2 - 451296*t - 22464]

Pólya walks in dimension 15

(Thanks to B. Salvy)

Task: Compute a certain solution at $z=1$

In [37]:
from ore_algebra.examples import polya
dim = 15
polya.dop[dim]
Out[37]:
(269276305858560000*z^30 - 80950584950784000*z^28 + 3629735201193984*z^26 - 57080630763520*z^24 + 409981265408*z^22 - 1499829760*z^20 + 2871232*z^18 - 2720*z^16 + z^14)*Dz^15 + (60587168818176000000*z^29 - 16999622839664640000*z^27 + 707798364232826880*z^25 - 10274513537433600*z^23 + 67646908792320*z^21 - 224974464000*z^19 + 387616320*z^17 - 326400*z^15 + 105*z^13)*Dz^14 + (5937542544181248000000*z^28 - 1551046657578762240000*z^26 + 59789931630571929600*z^24 - 798278106947641344*z^22 + 4796169345443840*z^20 - 14416495588352*z^18 + 22181540160*z^16 - 16423904*z^14 + 4550*z^12)*Dz^13 + (334481563322210304000000*z^27 - 81131900756901396480000*z^25 + 2886264219850496409600*z^23 - 35303727569401454592*z^21 + 192591205900620800*z^19 - 519787439069184*z^17 + 707833651200*z^15 - 454991264*z^13 + 106470*z^11)*Dz^12 + (12041336279599570944000000*z^26 - 2704250458682470563840000*z^24 + 88474218684064461619200*z^22 - 987113006841956179968*z^20 + 4862067116732345856*z^18 - 11694546001056256*z^16 + 13948983786816*z^14 - 7666274304*z^12 + 1479478*z^10)*Dz^11 + (291400337966309616844800000*z^25 - 60403894924097334804480000*z^23 + 1810516178566181073715200*z^21 - 18336647935059261603840*z^19 + 81033820655292578304*z^17 - 172214449187123200*z^15 + 177735512066112*z^13 - 81994323840*z^11 + 12662650*z^9)*Dz^10 + (4856672299438493614080000000*z^24 - 926086890189285634867200000*z^22 + 25324208695139969934950400*z^20 - 231569354145921664204800*z^18 + 911597715337047246336*z^16 - 1694812596013786624*z^14 + 1491624282894720*z^12 - 564676102848*z^10 + 67128490*z^8)*Dz^9 + (56198636607788283248640000000*z^23 - 9821600336936930947891200000*z^21 + 243895079188460801654784000*z^19 - 2001474041629081060638720*z^17 + 6961148053301631762432*z^15 - 11190413197947013632*z^13 + 8252818589658240*z^11 - 2491958589120*z^9 + 216627840*z^7)*Dz^8 + (449589092862306265989120000000*z^22 - 71725429608240620136038400000*z^20 + 1609161818725523357171712000*z^18 - 11770087561727055499100160*z^16 + 35825937728946311725056*z^14 - 49112927471281826304*z^12 + 29704887544668864*z^10 - 6898246766112*z^8 + 408741333*z^6)*Dz^7 + (2447762838917000781496320000000*z^21 - 354907365563571111670579200000*z^19 + 7152475405487676025208832000*z^17 - 46268511513163153122263040*z^15 + 121891943669163953209344*z^13 - 140140992502944986112*z^11 + 67667087104034880*z^9 - 11517430544256*z^7 + 420693273*z^5)*Dz^6 + (8811946220101202813386752000000*z^20 - 1155571219953730314672537600000*z^18 + 20785112953613635862082355200*z^16 - 117852071549595952781721600*z^14 + 265216641208694799482880*z^12 - 250559836896574725120*z^10 + 93303962552021952*z^8 - 10897207743264*z^6 + 210766920*z^4)*Dz^5 + (20027150500230006394060800000000*z^19 - 2362573936582134301458432000000*z^17 + 37651520030507075820847104000*z^15 - 185180175371863397892096000*z^13 + 350420651229762451537920*z^11 - 265169341550020423680*z^9 + 72752913986864640*z^7 - 5304232959840*z^5 + 42355950*z^3)*Dz^4 + (26702867333640008525414400000000*z^18 - 2816819552897219199762432000000*z^16 + 39443088818742150876364800000*z^14 - 166217350344029675008819200*z^12 + 259441279222310968688640*z^10 - 152179348773380567040*z^8 + 28891552538695680*z^6 - 1138125841504*z^4 + 2375101*z^2)*Dz^3 + (18486600461750775132979200000000*z^17 - 1732080379790685408067584000000*z^15 + 21106053778472440493506560000*z^13 - 75098351566004361206169600*z^11 + 94387689349096767160320*z^9 - 41091696296267489280*z^7 + 4930837635294720*z^5 - 82548913344*z^3 + 16383*z)*Dz^2 + (5281885846214507180851200000000*z^16 - 436217453243899431616512000000*z^14 + 4573788149507776189562880000*z^12 - 13497729262079414265446400*z^10 + 13245706827488316948480*z^8 - 4031063063164477440*z^6 + 265858264373760*z^4 - 1204325184*z^2 + 1)*Dz + 352125723080967145390080000000*z^15 - 25412342171473281024000000000*z^13 + 226231973017884794290176000*z^11 - 541574695562332078080000*z^9 + 398335457415580876800*z^7 - 77667722566041600*z^5 + 2222757642240*z^3 - 983040*z

⚐ Numerics: Logarithms

In [38]:
polya.dop[dim].leading_coefficient().factor()
Out[38]:
(2*z - 1) * (2*z + 1) * (6*z - 1) * (6*z + 1) * (10*z - 1) * (10*z + 1) * (14*z - 1) * (14*z + 1) * (18*z - 1) * (18*z + 1) * (22*z - 1) * (22*z + 1) * (26*z - 1) * (26*z + 1) * (30*z - 1) * (30*z + 1) * z^14
In [39]:
polya.dop[dim].local_basis_monomials(0)
Out[39]:
[1/87178291200*log(z)^14,
 1/6227020800*log(z)^13,
 1/479001600*log(z)^12,
 1/39916800*log(z)^11,
 1/3628800*log(z)^10,
 1/362880*log(z)^9,
 1/40320*log(z)^8,
 1/5040*log(z)^7,
 1/720*log(z)^6,
 1/120*log(z)^5,
 1/24*log(z)^4,
 1/6*log(z)^3,
 1/2*log(z)^2,
 log(z),
 1]
In [40]:
ini = [0]*(dim - 1) + [1]
# 90 s in Nov. 2017
%time polya.dop[dim].numerical_solution(ini, [0,1/(2*dim)], 1e-100)

CPU times: user 4.3 s, sys: 43.5 ms, total: 4.34 s
Wall time: 4.31 s
Out[40]:
[1.03720412092152168870320272201726012686492785463956609556680543456793165591178031741211 +/- 7.28e-87]

Iterated integrals

After J. Ablinger, J. Blümlein, C. G. Raab, and C. Schneider, Iterated Binomial Sums and their Associated Iterated Integrals, Journal of Mathematical Physics 55(11), 2014.

$$\int_{0}^1 dx_1 \, w_1(x_1) \int_{x_1}^1 dx_2 \, w_2(x_2) \, \cdots \! \int_{x_{n-1}}^1 dx_n \, w_n(x_n)$$
In [41]:
from ore_algebra.examples import iint
i = 69; iint.word[i]
Out[41]:
[x/(x - 1), 2/(sqrt(4*x - 1)*x), 2/(sqrt(4*x - 1)*x), -1/(x - 1), -1/(x - 1)]
In [42]:
dop = iint.diffop(iint.word[i])
dop
Out[42]:
(x^9 - 13/4*x^8 + 15/4*x^7 - 7/4*x^6 + 1/4*x^5)*Dx^6 + (27/2*x^8 - 35*x^7 + 30*x^6 - 9*x^5 + 1/2*x^4)*Dx^5 + (101/2*x^7 - 397/4*x^6 + 57*x^5 - 17/2*x^4 + 1/4*x^3)*Dx^4 + (111/2*x^6 - 303/4*x^5 + 45/2*x^4 + 3/4*x^3 - 3/4*x^2)*Dx^3 + (12*x^5 - 37/4*x^4 + 1/4*x^3 - 3/2*x^2 + 3/2*x)*Dx^2 + (-1/4*x^2 + 3/2*x - 3/2)*Dx

⚐ Numerics: Advanced interface

In [43]:
iint.ini[i]
Out[43]:
[0, 0, 0, -1/3, 1/3*I*pi + 2/3, -2/3*I*pi + 1/6*pi^2 - 11/12]
In [44]:
roots = dop.leading_coefficient().roots(AA, multiplicities=False)
sing = list(reversed(sorted([s for s in roots if 0 < s < 1])))
sing
Out[44]:
[0.25000000000000000?]
In [45]:
from ore_algebra.analytic.path import Point
path = [1] + [Point(s, dop, outgoing_branch=(0,-1)) for s in sing] + [0]
%time dop.numerical_solution(iint.ini[i], path, 1e-100)

CPU times: user 745 ms, sys: 41.5 ms, total: 786 ms
Wall time: 732 ms
Out[45]:
[-0.97080469562493124056011987537954344846933233520382808407829772093922068543247104239693854364491332570 +/- 9.81e-102] + [+/- 2.77e-102]*I

Bonus examples

Transition matrices

Compacted binary trees of fixed right height

After A. Genitrini, B. Gittenberger, M. Kauers, M. Wallner, Asymptotic Enumeration of Compacted Binary Trees, 2017

Exponential generating function of compacted binary trees:

In [46]:
from ore_algebra.examples import cbt
QQ[['z']](cbt.egf)
Out[46]:
1 + z + 3/2*z^2 + 5/2*z^3 + 37/8*z^4 + 373/40*z^5 + 4829/240*z^6 + 76981/1680*z^7 + 293057/2688*z^8 + 32536277/120960*z^9 + 827662693/1209600*z^10

Annihilator of EGF of CBT of right height ≤ 5:

In [47]:
dop = cbt.dop[5]; dop
Out[47]:
(-4*z^3 + 10*z^2 - 6*z + 1)*Dz^6 + (9*z^3 - 58*z^2 + 75*z - 21)*Dz^5 + (-6*z^3 + 78*z^2 - 184*z + 95)*Dz^4 + (z^3 - 33*z^2 + 141*z - 110)*Dz^3 + (3*z^2 - 32*z + 40)*Dz^2 + (z - 3)*Dz
In [48]:
dop.leading_coefficient().roots(QQbar, multiplicities=False)
Out[48]:
[0.2928932188134525?, 0.50000000000000000?, 1.707106781186548?]

⚐ Numerics: Algebraic points and exponents

In [49]:
s = _[0]
dop.local_basis_monomials(s)
Out[49]:
[1,
 z - 0.2928932188134525?,
 (z - 0.2928932188134525?)^1.771446609406727?,
 (z - 0.2928932188134525?)^2,
 (z - 0.2928932188134525?)^3,
 (z - 0.2928932188134525?)^4]
In [50]:
1/(4*(cos(pi/8)^2)).n()
Out[50]:
0.292893218813452

Singularity analysis implies

$$\#\{\text{cbt of rh ≤ 5}\} \sim \kappa \, n! α^n n^β, \qquad α = 4 \cos²(\pi/8), \quad β=-(α+5)/8$$

Singularity analysis implies

$$\#\{\text{cbt of rh ≤ 5}\} \sim \kappa_k n! α^n n^β, \qquad α = 4 \cos²(\pi/8), \quad β=-(α+5)/8$$
In [51]:
dop.local_basis_monomials(0)
Out[51]:
[1, z, z^2, z^3, z^4, z^5]
In [52]:
ini = list(cbt.egf[:6])
%time c = (dop.numerical_transition_matrix([0,s])*vector(ini))[2]
c

CPU times: user 671 ms, sys: 16.5 ms, total: 688 ms
Wall time: 662 ms
Out[52]:
[31.42732179257817 +/- 3.99e-15] + [27.45408401646515 +/- 4.94e-15]*I
In [53]:
expo = dop.local_basis_monomials(s)[2].op[1]
CBF(-s)^expo*c/CBF(-expo).gamma()
Out[53]:
[1.6248570260793 +/- 3.26e-14] + [+/- 9.96e-15]*I

A conjecture of M. Kontsevich

(via D. van Straten and A. Bostan)

  • $L = D_x\,x\,(x-1)\,(x-t)\,D_x + x$ for $t > 0$ (small)
  • $M(t)$ = transition matrix along a simple loop around $\{0, t\}$
  • $λ(t), \barλ(t)$ = eigenvalues of $M(t)$

Then

$$t \mapsto \left(\frac{\log(λ(t))}{2πi}\right)^2$$

is analytic at $0$, and its Taylor coefficients are rationals with small denominators.

Task: Compute that series (heuristically!)

⚐ Numerics: Decent high precision evaluation at singularities

In [54]:
terms = 20
sz = ceil((terms + 1)^2 * sqrt(ZZ(terms + 1).nbits())/2)
prec = terms*(sz + 4)
hprec = prec + 100
C = ComplexField(hprec)
sz, prec
Out[54]:
(494, 9960)
In [55]:
Pol.<x> = QQ[]
Dop.<Dx> = OreAlgebra(Pol)
t0 = 2^(-sz)
L = Dx * (x*(x-1)*(x-t0)) * Dx + x
L
Out[55]:
(x^3 - 51146728248377216718956089012931236753385031969422887335676427626502090568823039920051095192592252455482604439493126109519019633529459266458258243585/51146728248377216718956089012931236753385031969422887335676427626502090568823039920051095192592252455482604439493126109519019633529459266458258243584*x^2 + 1/51146728248377216718956089012931236753385031969422887335676427626502090568823039920051095192592252455482604439493126109519019633529459266458258243584*x)*Dx^2 + (3*x^2 - 51146728248377216718956089012931236753385031969422887335676427626502090568823039920051095192592252455482604439493126109519019633529459266458258243585/25573364124188608359478044506465618376692515984711443667838213813251045284411519960025547596296126227741302219746563054759509816764729633229129121792*x + 1/51146728248377216718956089012931236753385031969422887335676427626502090568823039920051095192592252455482604439493126109519019633529459266458258243584)*Dx + x
In [56]:
m1 = L.numerical_transition_matrix([0, t0/2], 2^(-prec)).change_ring(C)
m2 = L.numerical_transition_matrix([t0, t0/2], 2^(-prec)).change_ring(C)
delta = matrix(C, [[1, 0], [2*pi*I, 1]])
mat = m1*delta*~m1*m2*delta*~m2
In [57]:
tr = mat.trace().real()
Pol.<la> = C[]
char = (la+1/la-tr).numerator()
rt = char.roots(multiplicities=False)[0]
a = (rt.log()/C(2*I*pi))^2
a
Out[57]:
9.556619453472961320858516937971887816879016410846290172126354839510701168391942496813468949482302488162252696886520244471966592218665291723177173168199392132406276973789705340073690947932553910124796231579916084856421141195367398954924504111462358438957864101423930407260331976176725065395140834421424427590908128799409817584336622391673407243198108939261965265090478623967959224212373223620909718920076884504325493512957093231290793326009149024368397362080608773817261344914045941037215573589929466261952514844246399651185700946374898667516902390725364103154432169875542859847614933659911701076729537808003318825249648022560871556004599073286958174424349446361789734983296900487540550803091366045082974239981262926261619181189568107287215173880557039860815668766673099266893087110588547841898487650663644043980236519485686759761011681600245421501121353018152348931578768245258254780107760973263563775281863663012073352722917223246017584235308522688246510732993890785886249250725937793549943881255036341563534916875647477037886097657249792008428528870297161215131300381890686100703684816317908759813352398646880833758301441006787007561756610067089326048633636601655581680294114920232193088316072162202643820964032906843841817185410145450874756553779319769681938557795157701492901383950831572956075682891800333143265314437926822965689897336891281115810393634163501579557476432032208140913208912647082900790120127010816341829130398447859634549497723253118914899666245483559245438168221086949951270207323864052877486421440537213113570849672571099832804454703454119688641196701003291562809385222820524421206480989940681005905092328271286077626602437533505092641703146035760742921680842263781843854873991144874044147403548564686612183559041672386026367636029777757668704070453267111227724964978924932012765810717992589937434808382663799797773043364171093373260596575966020347101598724092360916823714857848399257003615082945323894832711059841119689648486033169179927891587569070063188570258161818486202606231937307760766826338422121070834991233113786550649092555451055378950142673684858167685590606256915969226571204633960531410689676714870414271577763534727803754179156714701688996696004086261971352341976630738528832473896535081911494009668569113201638788353744280283060834315636625364415819665437542451448080606771694562471371975196008988383248568008014743171957640899885078174330701068772917303994305295373221848857890630191210892023261909597212215222580922194025813235619871426751892331277173873023952664284270669470966003915290659729911588342236100212184270452170452531046709945743782625626814519771418724775958330968324506098946055581698099590764388244028515969168881038955796413001853465395680656243512452699320906509054087659901425563057011885470050958906345477915933201856109902467292184386952372661834542279554976803521112299906835150188799637765061595534103641294316535688856568007252802956390958368379318710244523774844628336072323710306061629273180267967815264454788464436467987614675704690957518274321357948708194043185e-299
In [58]:
coeffs = []; den_ratios = []; cur = a
for k in range(terms):
    rat = QQ(pari.bestappr(cur, 2^(sz/2+10)))
    cur = (cur - rat)/t0
    if k >= 2:
        den_ratios.append(rat.denom()/coeffs[-1].denom())
    coeffs.append(rat)
coeffs
Out[58]:
[0,
 0,
 1/4,
 1/24,
 101/576,
 239/17280,
 19153/115200,
 -1516283/72576000,
 23167560743/121927680000,
 -5350452180523/76814438400000,
 8122785754979827/32262064128000000,
 -10922037427834714189/74525368135680000000,
 257615133666208067057417/688614401573683200000000,
 -5719273111411836892974100997/20679090479257706496000000000,
 749577901737131766662423170141529/1241986174184217852149760000000000,
 -56589264915861458106843233145366159317/111890534432256186300171878400000000000,
 166060654964554378941594941573970938041843381/161283491952031357180519752400896000000000000,
 -227464761949125949651795174312351837947239351830471/247010506429294584462681416394544250880000000000000,
 1379581218345710272912991297578804281398062616806207006247/756608001823315069884260939301472713100492800000000000000,
 -37154931571287997286276581805653700634917295203190508379056930613/22016589207616772850617000970999305581601157021696000000000000000]
In [59]:
den_ratios
Out[59]:
[4,
 6,
 24,
 30,
 20/3,
 630,
 1680,
 630,
 420,
 2310,
 9240,
 30030,
 60060,
 90090,
 1441440,
 1531530,
 3063060,
 29099070]

Desingularization

In [60]:
from ore_algebra.examples import periods
dop = periods.allODEs[2][0].numerator()
dop
Out[60]:
(t^14 - 280163/1250*t^13 - 77823473621/3750000*t^12 + 15467808452531/58593750*t^11 - 25196309686447913/19531250000*t^10 - 668120393063773277/915527343750*t^9 + 1960569648056681994201/61035156250000*t^8 - 4671428124563366603059/30517578125000*t^7 + 41197718959797815058559/95367431640625*t^6 - 164817439510375233136933/190734863281250*t^5 + 5927091704553171594142627/4768371582031250*t^4 - 8659852913093432184666047/7152557373046875*t^3 + 5049362927700459914731448/7152557373046875*t^2 - 685309305030323264047472/11920928955078125*t - 18187783435081723962509824/35762786865234375)*Dt^5 + (29/2*t^13 - 8613043/2500*t^12 - 174431849123/500000*t^11 + 1255822655957817/312500000*t^10 - 20280225824735266/1220703125*t^9 - 49190254079706055749/976562500000*t^8 + 78414082510831851698897/122070312500000*t^7 - 18568806332381461314869/7629394531250*t^6 + 4003810338808418269811661/762939453125000*t^5 - 27771814401571351451778/3814697265625*t^4 + 11867297379816696715186737/1907348632812500*t^3 - 307053942520796638296658/476837158203125*t^2 - 1871856137277813387262933/476837158203125*t + 6163568693739391139347276/2384185791015625)*Dt^4 + (475/8*t^12 - 30050227/2000*t^11 - 16988836335053/10000000*t^10 + 21958673179999589/1250000000*t^9 - 9131560294338258463/156250000000*t^8 - 401564615212972941641/976562500000*t^7 + 104294200346934559073119/30517578125000*t^6 - 637131873161718618633403/61035156250000*t^5 + 4868593843047491215960457/305175781250000*t^4 - 426971833579701247214697/38146972656250*t^3 - 775947208574144990906417/953674316406250*t^2 + 9713861736043758399614103/953674316406250*t - 3186441544992332414987113/476837158203125)*Dt^3 + (1185/16*t^11 - 80462511/4000*t^10 - 52523948963397/20000000*t^9 + 120582837906110277/5000000000*t^8 - 18344488178046465327/312500000000*t^7 - 6879971713312897388329/7812500000000*t^6 + 2687398835047244023617471/488281250000000*t^5 - 130498326095692372331277/9765625000000*t^4 + 1786945003773501791092069/152587890625000*t^3 + 2499771123849829666054977/610351562500000*t^2 - 100956317858140812799725681/7629394531250000*t + 8521916368516547564271043/953674316406250)*Dt^2 + (5145/256*t^10 - 380739891/64000*t^9 - 305932960216191/320000000*t^8 + 153851731551379239/20000000000*t^7 - 6655819554769778547/625000000000*t^6 - 6447339574295055155727/15625000000000*t^5 + 157321670664515867630799/78125000000000*t^4 - 145611308765873254822287/39062500000000*t^3 - 29971841950167023004357/122070312500000*t^2 + 3523442028617667203992833/610351562500000*t - 14063874730153571703882243/3814697265625000)*Dt + 105/512*t^9 - 8901123/128000*t^8 - 11004229733343/640000000*t^7 + 4759356296189817/40000000000*t^6 + 11832518636433327/625000000000*t^5 - 57291975299178874233/6250000000000*t^4 + 214707950483595841419/6250000000000*t^3 - 3411921799561771135773/78125000000000*t^2 - 3177202651293055513191/19531250000000*t + 111251064741292878982587/610351562500000
In [61]:
dop.leading_coefficient().roots(QQbar)
Out[61]:
[(-79.81779428070589?, 1),
 (-5.251936447658327?, 1),
 (-0.5142673948169302?, 1),
 (2.501571821316740?, 1),
 (2.754410642927720?, 1),
 (292.1291388587215?, 1),
 (0.1227946785361402? - 1.519364185445195?*I, 1),
 (0.1227946785361402? + 1.519364185445195?*I, 1),
 (0.2812140893416301? - 1.755844000088624?*I, 1),
 (0.2812140893416301? + 1.755844000088624?*I, 1),
 (1.749115902402969? - 1.266768315130795?*I, 1),
 (1.749115902402969? + 1.266768315130795?*I, 1),
 (4.011513729826854? - 4.25537902618774?*I, 1),
 (4.011513729826854? + 4.25537902618774?*I, 1)]
In [62]:
dop1 = dop.desingularize().numerator()
dop1.leading_coefficient().roots(QQbar)
Out[62]:
[(2.501571821316740?, 1),
 (2.754410642927720?, 1),
 (0.1227946785361402? - 1.519364185445195?*I, 1),
 (0.1227946785361402? + 1.519364185445195?*I, 1),
 (0.2812140893416301? - 1.755844000088624?*I, 1),
 (0.2812140893416301? + 1.755844000088624?*I, 1)]

Plots

In [63]:
from ore_algebra.analytic.function import DFiniteFunction
P.<x> = QQ[]
A.<Dx> = OreAlgebra(P)
f = DFiniteFunction(Dx^2 - x,
        [1/(gamma(2/3)*3^(2/3)), -1/(gamma(1/3)*3^(1/3))],
        name='my_Ai')
f.plot((-5,5))
Out[63]: