Add Spiro version 0.01.

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)