# Synthesize a procedure to numerically integrate the 3rd order poly spiral tex = False if tex: mulsym = ' ' else: mulsym = ' * ' class Poly: def __init__(self, p0, coeffs): self.p0 = p0 self.coeffs = coeffs def eval(self, x): y = x ** self.p0 z = 0 for c in 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 i >= poly0.p0 - p0 and i < lp0 + poly0.p0 - p0: c += poly0.coeffs[i + p0 - poly0.p0] if i >= poly1.p0 - p0 and i < lp1 + poly1.p0 - p0: c += poly1.coeffs[i + p0 - poly1.p0] coeffs.append(c) return Poly(p0, coeffs) def pr(str): if tex: print str, '\\\\' else: print '\t' + str + ';' def prd(str): if tex: print str, '\\\\' else: print '\tdouble ' + str + ';' 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 != None and t1 != None: terms.append(t0 + mulsym + t1) if terms == []: result.append(None) else: var = basename % i if (j % 2 == 0) or not suppress_odd: prd(var + ' = ' + ' + '.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 != None and t1 != 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 != None: if tex: terms.append(t + '^2') else: terms.append(t + mulsym + t) if terms == []: result.append(None) else: var = basename % i prd(var + ' = ' + ' + '.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(tp[1] + ' = km0') prd(tp[2] + ' = .5 * km1') prd(tp[3] + ' = (1./6) * km2') prd(tp[4] + ' = (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] != None: if tex: tmp.append('\\frac{%s}{%.0f}' % (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 + ' ' str = pref + (' ' + sign + ' ').join(tmp) else: str = var + ' ' + sign + '= ' + ' + '.join(tmp) if var == 'u': us.append(str) else: vs.append(str) 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__': mkspiro(12)