numintsynth.py

# Synthesize a procedure to numerically integrate the 3rd order poly spiral
# Edited for Rust by Fredrick R. Brennan 
# (c) 2007–2022 Raph Levien, Fredrick R. Brennan. Apache 2.0 licensed.

from __future__ import division
from __future__ import print_function

tex = False

if tex:
    mulsym = ' '
else:
    mulsym = ' * '


class Poly:
    def __init__(self, p0, coeffs):
        self.p0 = p0
        self.coeffs = coeffs

    def eval(self, x):  # TODO: method was broken. Investigate remove possibility.
        y = x ** self.p0
        z = 0
        for c in self.coeffs:
            z += y * c
            y *= x
        return z


def add(poly0, poly1, nmax):
    lp0 = len(poly0.coeffs)
    lp1 = len(poly1.coeffs)
    p0 = min(poly0.p0, poly1.p0)
    n = min(max(poly0.p0 + lp0, poly1.p1 + lp1), nmax) - p0
    if n <= 0:
        return Poly(0, [])
    coeffs = []
    for i in range(n):
        c = 0
        if poly0.p0 - p0 <= i < lp0 + poly0.p0 - p0:
            c += poly0.coeffs[i + p0 - poly0.p0]
        if poly1.p0 - p0 <= i < lp1 + poly1.p0 - p0:
            c += poly1.coeffs[i + p0 - poly1.p0]
        coeffs.append(c)
    return Poly(p0, coeffs)


def pr(string):
    if tex:
        print(string, '\\\\')
    else:
        print('\t\t' + string + ';')


def prd(string):
    if tex:
        print(string, '\\\\')
    else:
        print('\t\t' + string + ';')


def polymul(p0, p1, degree, basename, suppress_odd=False):
    result = []
    for i in range(min(degree, len(p0) + len(p1) - 1)):
        terms = []
        for j in range(i + 1):
            if j < len(p0) and i - j < len(p1):
                t0 = p0[j]
                t1 = p1[i - j]
                if t0 is not None and t1 is not None:
                    terms.append(t0 + mulsym + t1)
        if not terms:
            result.append(None)
        else:
            var = basename % i  # type: str
            if (j % 2 == 0) or not suppress_odd:
                prd('let ' + var + ': f64 = ' + ' + '.join(terms))
            result.append(var)
    return result


def polysquare(p0, degree, basename):
    result = []
    for i in range(min(degree, 2 * len(p0) - 1)):
        terms = []
        for j in range((i + 1) // 2):
            if i - j < len(p0):
                t0 = p0[j]
                t1 = p0[i - j]
                if t0 is not None and t1 is not None:
                    terms.append(t0 + mulsym + t1)
        if len(terms) >= 1:
            if tex and len(terms) == 1:
                terms = ['2. ' + terms[0]]
            else:
                terms = ['2.' + mulsym + '(' + ' + '.join(terms) + ')']
        if (i % 2) == 0:
            t = p0[i // 2]
            if t is not None:
                if tex:
                    terms.append(t + '^2')
                else:
                    terms.append(t + mulsym + t)
        if not terms:
            result.append(None)
        else:
            var = basename % i  # type: str
            prd('let ' + var + ': f64 = ' + ' + '.join(terms))
            result.append(var)
    return result


def mkspiro(degree):
    if tex:
        us = ['u = 1.']
        vs = ['v =']
    else:
        us = ['u = 1.']
        vs = ['v = 0.']
    if tex:
        tp = [None, 't_{11}', 't_{12}', 't_{13}', 't_{14}']
    else:
        tp = [None, 't1_1', 't1_2', 't1_3', 't1_4']
    if tex:
        prd(tp[1] + ' = k_0')
        prd(tp[2] + ' = \\frac{k_1}{2}')
        prd(tp[3] + ' = \\frac{k_2}{6}')
        prd(tp[4] + ' = \\frac{k_3}{24}')
    else:
        prd('let ' + tp[1] + ': f64 = km0')
        prd('let ' + tp[2] + ': f64 = 0.5 * km1')
        prd('let ' + tp[3] + ': f64 = (1./6.) * km2')
        prd('let ' + tp[4] + ': f64 = (1./24.) * km3')
    tlast = tp
    coef = 1.
    for i in range(1, degree - 1):
        tmp = []
        tcoef = coef
        # print(tlast)
        for j in range(len(tlast)):
            c = tcoef / (j + 1)
            if (j % 2) == 0 and tlast[j] is not None:
                if tex:
                    tmp.append('\\frac{%s}{%d}.' % (tlast[j], 1. / c))
                else:
                    if c < 1e-9:
                        cstr = '%.16e' % c
                    else:
                        cstr = '(1./%d.)' % int(.5 + (1. / c))
                    tmp.append(cstr + ' * ' + tlast[j])
            tcoef *= .5
        if tmp:
            sign = ('+', '-')[(i // 2) % 2]
            var = ('u', 'v')[i % 2]
            if tex:
                if i == 1:
                    pref = ''
                else:
                    pref = sign + ' '
                string = pref + (' ' + sign + ' ').join(tmp)
            else:
                string = var + ' ' + sign + '= ' + ' + '.join(tmp)
            if var == 'u':
                us.append(string)
            else:
                vs.append(string)
        if i < degree - 1:
            if tex:
                basename = 't_{%d%%d}' % (i + 1)
            else:
                basename = 't%d_%%d' % (i + 1)
            if i == 1:
                tnext = polysquare(tp, degree - 1, basename)
                t2 = tnext
            elif i == 3:
                tnext = polysquare(t2l, degree - 1, basename)
            elif (i % 2) == 0:
                tnext = polymul(tlast, tp, degree - 1, basename, True)
            else:
                tnext = polymul(t2l, t2, degree - 1, basename)
            t2l = tlast
            tlast = tnext
        coef /= (i + 1)
    if tex:
        pr(' '.join(us))
        pr(' '.join(vs))
    else:
        for u in us:
            pr(u)
        for v in vs:
            pr(v)


if __name__ == '__main__':
    print("""
pub fn integrate_spiro(ks: [f64; 4]) -> [f64; 2] {
    let th1: f64 = ks[0];
    let th2: f64 = 0.5 * ks[1];
    let th3: f64 = 1.0 / 6. * ks[2];
    let th4: f64 = 1.0 / 24. * ks[3];
    let mut x: f64 = 0.;
    let mut y: f64 = 0.;
    let ds: f64 = 1.0 / 4 as f64;
    let ds2: f64 = ds * ds;
    let ds3: f64 = ds2 * ds;
    let k0: f64 = ks[0] * ds;
    let k1: f64 = ks[1] * ds;
    let k2: f64 = ks[2] * ds;
    let k3: f64 = ks[3] * ds;
    let mut i: isize = 0;
    let mut s: f64 = 0.5 * ds - 0.5;
    while i < 4 {
        let mut u: f64;
        let mut v: f64;
        let km0: f64 = 0.;
        let km1: f64 = 0.;
        let km2: f64 = 0.;
        let km3: f64 = 0.;""")
    mkspiro(128)
    print("""
//
        let th: f64 = (((th4 * s + th3) * s + th2) * s + th1) * s;
        let cth: f64 = (th).cos();
        let sth: f64 = (th).sin();
        x += cth * u - sth * v;
        y += cth * v + sth * u;
        s += ds;

        i += 1
    }
    let mut xy: [f64; 2] = [0.; 2];
    xy[0] = x * ds;
    xy[1] = y * ds;
    xy
}""")
    print("""fn main() { eprintln!("{:?}", integrate_spiro([1.0, 2.0, 3.0, 4.0])) }""")