diff options
Diffstat (limited to 'third_party/spiro/curves/numintsynth.py')
-rw-r--r-- | third_party/spiro/curves/numintsynth.py | 176 |
1 files changed, 176 insertions, 0 deletions
diff --git a/third_party/spiro/curves/numintsynth.py b/third_party/spiro/curves/numintsynth.py new file mode 100644 index 0000000..9a51b66 --- /dev/null +++ b/third_party/spiro/curves/numintsynth.py @@ -0,0 +1,176 @@ +# 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) |