--- /dev/null
+#!/usr/bin/env python3
+
+from pypol import *
+
+def affine_derivative_closure(T, x0s):
+
+ xs = [Symbol("{}'".format(x0.name)) for x0 in x0s]
+ dxs = [Symbol('d{}'.format(x0.name)) for x0 in x0s]
+ k = Symbol('k')
+
+ for x in T.symbols:
+ x = Symbol(x)
+ assert x in x0s + xs
+ for dx in dxs:
+ assert dx.name not in T.symbols
+ assert k.name not in T.symbols
+
+ T0 = T
+
+ T1 = T0
+ for i, x0 in enumerate(x0s):
+ x, dx = xs[i], dxs[i]
+ T1 &= Eq(dx, x - x0)
+
+ T2 = T1.project_out(T0.symbols)
+
+ T3_eqs = []
+ T3_ins = []
+ for T2_eq in T2.equalities:
+ c = T2_eq.constant
+ T3_eq = T2_eq + (k - 1) * c
+ T3_eqs.append(T3_eq)
+ for T2_in in T2.inequalities:
+ c = T2_in.constant
+ T3_in = T2_in + (k - 1) * c
+ T3_ins.append(T3_in)
+ T3 = Polyhedron(T3_eqs, T3_ins)
+ T3 &= Ge(k, 0)
+
+ T4 = T3.project_out([k])
+ for i, dx in enumerate(dxs):
+ x0, x = x0s[i], xs[i]
+ T4 &= Eq(dx, x - x0)
+ T4 = T4.project_out(dxs)
+
+ return T4
+
+i0, j0, i, j = symbols(['i', 'j', "i'", "j'"])
+T = Eq(i, i0 + 2) & Eq(j, j0 + 1)
+
+print('T =', T)
+Tstar = affine_derivative_closure(T, [i0, j0])
+print('T* =', Tstar)