Numeric evaluation for iterated integrals by SageMath

The following is a source code for a arbitrary iterated integrals import sys sys.setrecursionlimit(1000000) from functools import lru_cache C = ComplexField(300) @lru_cache(maxsize = None) def ZetaN(N, ks, zs): if len(ks)==0: return C(1) if N<len(ks): return C(0) return ZetaN(N-1, ks, zs) + C(zs[-1]^N) / (N^ks[-1]) * ZetaN(N-1, ks[:-1], zs[:-1])

def Zeta(ks, zs): for N in (len(ks)..100000000): if ZetaN(N,ks,zs)==ZetaN(N+1,ks,zs)==ZetaN(N+2,ks,zs)==ZetaN(N+3,ks,zs): return ZetaN(N,ks,zs)

def ks_and_zs_from_bs(bs): zs,ks = [1/bs[0]], [1] for b in bs[1:]: if b==0: ks[-1] += 1 else: ks.append(1) zs[-1] *= b           zs.append(1/b) return tuple(ks), tuple(zs)

def iterated_integral_convergent_series_case(p,q,bs): assert all(b==p or abs( (b-p)/(q-p) ) > 1 for b in bs ) p, q, bs = 0, 1, [(b-p)/(q-p) for b in bs] ks, zs = ks_and_zs_from_bs(bs) return (-1)^len(ks) * Zeta(ks, zs)

@lru_cache(maxsize = None) def iterated_integral(p,q,bs): if len(bs)==0: return C(1) if p==q: return C(0) D = 1.2 if all(b==p or abs( (b-p)/(q-p) ) > D for b in bs ): return iterated_integral_convergent_series_case(p, q, bs) if all(b==q or abs( (b-q)/(p-q) ) > D for b in bs ): return (-1)^len(bs) * iterated_integral_convergent_series_case(q, p, bs[::-1]) return sum( iterated_integral(p, (p+q)/2, bs[:i] ) * iterated_integral( (p+q)/2, q, bs[i:] )  for i in (0..len(bs)) )

iterated_integral(p = 0, q = 1, bs = (C(99,20)/101, C(99,20)/101, 0))