Initial setup for period computation (#13)

Introduces:
 - class for Picard--Fuchs operators and their ideals
 - class for periods (their complex linear combinations in the Frobenius bases)

Implements:
 - method to obtain operator from period
 - method to get power series solution (at given indicials) to PF ideal

---------

Co-authored-by: Julian Piribauer <julian.piribauer@gmail.com>
Reviewed-on: #13

README.md

@@ -12,12 +12,21 @@ elliptic_curve_D = ToricPolytopeProjectiveSpace([1, 2, 3], model_name="elliptic_
 
 CY3_quintic = ToricPolytopeProjectiveSpace([1, 1, 1, 1, 1], model_name="quintic")
 CY3_bicubic = ToricPolytopeCICY([[3, 3]], model_name="bi-cubic")
-CICY3_two_parameter_manual_nef = ToricPolytope([[1,0,0,0,0,0],[0,1,0,0,0,0],[0,0,1,0,0,0],
-                                                [0,0,0,1,0,0],[0,0,0,0,1,0],[0,0,0,0,0,1],
-                                                [-1,-1,0,0,0,0],[0,0,-1,-1,-1,-1]],
-                                                nef_partition=[[0,1,2,3], [4,5,6,7]])
+CICY3_two_parameter_manual_nef = ToricPolytope(
+    [
+        [1, 0, 0, 0, 0, 0],
+        [0, 1, 0, 0, 0, 0],
+        [0, 0, 1, 0, 0, 0],
+        [0, 0, 0, 1, 0, 0],
+        [0, 0, 0, 0, 1, 0],
+        [0, 0, 0, 0, 0, 1],
+        [-1, -1, 0, 0, 0, 0],
+        [0, 0, -1, -1, -1, -1],
+    ],
+    nef_partition=[[0, 1, 2, 3], [4, 5, 6, 7]],
+)
 
-CICY5_two_parameter = ToricPolytopeCICY([[6, 1],[0, 2]])
+CICY5_two_parameter = ToricPolytopeCICY([[6, 1], [0, 2]])
 ```
 
 The discriminant factors and topological data, e.g. for the quintic, can then be computed with the methods below.
@@ -71,4 +80,72 @@ DEBUG:__main__:
 
 The result is saved in the folder `data/topdata` as a JSON file — giving a model name helps keeping
 track of these outputs.
-Note that for Calabi–Yau dimensions larger than four, the additional 
+Note that for Calabi–Yau dimensions larger than four, the additional
+
+## period_computation
+
+`period_computation.sage` provides classes for working with Picard–Fuchs operators and their
+period solutions: `PFOperator`, `PFIdeal` and `Period`, together with Ansatz variants of the first and
+last (`PFOperatorAnsatz`, `PeriodAnsatz`) used to search for unknown operators or periods of a given
+z- and theta-degree.
+
+A `PFOperator` is parsed from a string in the variables `z0, ..., z<n-1>` and `theta0, ..., theta<n-1>`,
+the logarithmic derivatives theta_i = z_i d/dz_i. For example, the quintic's Picard&ndash;Fuchs operator:
+
+```python
+L = PFOperator(
+    "theta0^4 - 3125*z0*theta0^4 - 6250*z0*theta0^3 - 4375*z0*theta0^2 - 1250*z0*theta0 - 120*z0",
+    no_variables=1,
+)
+
+L.simplify().operator_string
+```
+
+```term
+'-5*(5*theta0 + 4)*(5*theta0 + 3)*(5*theta0 + 2)*(5*theta0 + 1)*z0 + theta0^4'
+```
+
+An operator (or a `PFIdeal` of several) can be solved for its power series solution at given indicial
+exponents and order.
+
+```python
+ideal = PFIdeal([L])
+period = ideal.find_power_series_solution(indicials=[0], order=3)[0]
+
+period.period_string
+```
+
+```term
+'168168000*z0^3 + 113400*z0^2 + 120*z0 + 1'
+```
+
+The reverse direction is supported too: given a `Period`, `find_annihilating_operators` searches for
+`PFOperator`s of a given z- and theta-degree that annihilate it, by solving an Ansatz of unknown
+coefficients via linear algebra. Both directions extend to several moduli, e.g. for the two-parameter
+model P_{2,2,2,1,1}[8]:
+
+```python
+M1 = PFOperator(
+    "theta1*(-2*theta0 + 2*theta1 - 1) + 2*(theta0 - 2*theta1 - 1)*(theta0 - 2*theta1)*z1",
+    no_variables=2,
+)
+M2 = PFOperator(
+    "theta0^2*(2*(theta0 - 2*theta1)*z1 - theta1) - 16*(2*theta0 + 1)*(4*theta0 + 1)*(4*theta0 + 3)*z0*z1",
+    no_variables=2,
+)
+ideal = PFIdeal([M1, M2])
+
+period = ideal.find_power_series_solution(indicials=[0, 1 / 2], order=6)[0]
+recovered = period.find_annihilating_operators(z_degree=1, theta_degree=2)
+
+period.period_string
+recovered[0].simplify().operator_string
+```
+
+```term
+'-1/45045*(60886425600*z0^3*z1^3 - 2767564800*z0^2*z1^4 + 100638720*z0*z1^5 - 14192640*z1^6 + 830269440*z0^2*z1^3 - 30750720*z0*z1^4 + 4193280*z1^5 - 242161920*z0^2*z1^2 + 9884160*z0*z1^3 - 1281280*z1^4 - 3459456*z0*z1^2 + 411840*z1^3 + 1441440*z0*z1 - 144144*z1^2 + 60060*z1 - 45045)*sqrt(z1)'
+'-2*(theta0 - 2*theta1)*(theta0 - 2*theta1 - 1)*z1 + (2*theta0 - 2*theta1 + 1)*theta1'
+```
+
+`recovered[0]` is, up to scale, `M1` &mdash; recovered purely from `M1`, `M2`'s shared power series
+solution.