src/live_effects/spiro.cpp

   1 /*
   2 ppedit - A pattern plate editor for Spiro splines.
   3 Copyright (C) 2007 Raph Levien
   4
   5 This program is free software; you can redistribute it and/or
   6 modify it under the terms of the GNU General Public License
   7 as published by the Free Software Foundation; either version 2
   8 of the License, or (at your option) any later version.
   9
  10 This program is distributed in the hope that it will be useful,
  11 but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 GNU General Public License for more details.
  14
  15 You should have received a copy of the GNU General Public License
  16 along with this program; if not, write to the Free Software
  17 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
  18 02110-1301, USA.
  19
  20 */
  21 /* C implementation of third-order polynomial spirals. */
  22
  23 #include <math.h>
  24 #include <stdlib.h>
  25 #include <string.h>
  26
  27 #include "bezctx_intf.h"
  28 #include "spiro.h"
  29
  30 struct spiro_seg_s {
  31     double x;
  32     double y;
  33     char ty;
  34     double bend_th;
  35     double ks[4];
  36     double seg_ch;
  37     double seg_th;
  38     double l;
  39 };
  40
  41 typedef struct {
  42     double a[11]; /* band-diagonal matrix */
  43     double al[5]; /* lower part of band-diagonal decomposition */
  44 } bandmat;
  45
  46 #ifndef M_PI
  47 #define M_PI            3.14159265358979323846  /* pi */
  48 #endif
  49
  50 int n = 4;
  51
  52 #ifndef ORDER
  53 #define ORDER 12
  54 #endif
  55
  56 /* Integrate polynomial spiral curve over range -.5 .. .5. */
  57 void
  58 integrate_spiro(const double ks[4], double xy[2])
  59 {
  60 #if 0
  61     int n = 1024;
  62 #endif
  63     double th1 = ks[0];
  64     double th2 = .5 * ks[1];
  65     double th3 = (1./6) * ks[2];
  66     double th4 = (1./24) * ks[3];
  67     double x, y;
  68     double ds = 1. / n;
  69     double ds2 = ds * ds;
  70     double ds3 = ds2 * ds;
  71     double k0 = ks[0] * ds;
  72     double k1 = ks[1] * ds;
  73     double k2 = ks[2] * ds;
  74     double k3 = ks[3] * ds;
  75     int i;
  76     double s = .5 * ds - .5;
  77
  78     x = 0;
  79     y = 0;
  80
  81     for (i = 0; i < n; i++) {
  82
  83 #if ORDER > 2
  84         double u, v;
  85         double km0, km1, km2, km3;
  86
  87         if (n == 1) {
  88             km0 = k0;
  89             km1 = k1 * ds;
  90             km2 = k2 * ds2;
  91         } else {
  92             km0 = (((1./6) * k3 * s + .5 * k2) * s + k1) * s + k0;
  93             km1 = ((.5 * k3 * s + k2) * s + k1) * ds;
  94             km2 = (k3 * s + k2) * ds2;
  95         }
  96         km3 = k3 * ds3;
  97 #endif
  98
  99         {
 100
 101 #if ORDER == 4
 102         double km0_2 = km0 * km0;
 103         u = 24 - km0_2;
 104         v = km1;
 105 #endif
 106
 107 #if ORDER == 6
 108         double km0_2 = km0 * km0;
 109         double km0_4 = km0_2 * km0_2;
 110         u = 24 - km0_2 + (km0_4 - 4 * km0 * km2 - 3 * km1 * km1) * (1./80);
 111         v = km1 + (km3 - 6 * km0_2 * km1) * (1./80);
 112 #endif
 113
 114 #if ORDER == 8
 115         double t1_1 = km0;
 116         double t1_2 = .5 * km1;
 117         double t1_3 = (1./6) * km2;
 118         double t1_4 = (1./24) * km3;
 119         double t2_2 = t1_1 * t1_1;
 120         double t2_3 = 2 * (t1_1 * t1_2);
 121         double t2_4 = 2 * (t1_1 * t1_3) + t1_2 * t1_2;
 122         double t2_5 = 2 * (t1_1 * t1_4 + t1_2 * t1_3);
 123         double t2_6 = 2 * (t1_2 * t1_4) + t1_3 * t1_3;
 124         double t3_4 = t2_2 * t1_2 + t2_3 * t1_1;
 125         double t3_6 = t2_2 * t1_4 + t2_3 * t1_3 + t2_4 * t1_2 + t2_5 * t1_1;
 126         double t4_4 = t2_2 * t2_2;
 127         double t4_5 = 2 * (t2_2 * t2_3);
 128         double t4_6 = 2 * (t2_2 * t2_4) + t2_3 * t2_3;
 129         double t5_6 = t4_4 * t1_2 + t4_5 * t1_1;
 130         double t6_6 = t4_4 * t2_2;
 131         u = 1;
 132         v = 0;
 133         v += (1./12) * t1_2 + (1./80) * t1_4;
 134         u -= (1./24) * t2_2 + (1./160) * t2_4 + (1./896) * t2_6;
 135         v -= (1./480) * t3_4 + (1./2688) * t3_6;
 136         u += (1./1920) * t4_4 + (1./10752) * t4_6;
 137         v += (1./53760) * t5_6;
 138         u -= (1./322560) * t6_6;
 139 #endif
 140
 141 #if ORDER == 10
 142         double t1_1 = km0;
 143         double t1_2 = .5 * km1;
 144         double t1_3 = (1./6) * km2;
 145         double t1_4 = (1./24) * km3;
 146         double t2_2 = t1_1 * t1_1;
 147         double t2_3 = 2 * (t1_1 * t1_2);
 148         double t2_4 = 2 * (t1_1 * t1_3) + t1_2 * t1_2;
 149         double t2_5 = 2 * (t1_1 * t1_4 + t1_2 * t1_3);
 150         double t2_6 = 2 * (t1_2 * t1_4) + t1_3 * t1_3;
 151         double t2_7 = 2 * (t1_3 * t1_4);
 152         double t2_8 = t1_4 * t1_4;
 153         double t3_4 = t2_2 * t1_2 + t2_3 * t1_1;
 154         double t3_6 = t2_2 * t1_4 + t2_3 * t1_3 + t2_4 * t1_2 + t2_5 * t1_1;
 155         double t3_8 = t2_4 * t1_4 + t2_5 * t1_3 + t2_6 * t1_2 + t2_7 * t1_1;
 156         double t4_4 = t2_2 * t2_2;
 157         double t4_5 = 2 * (t2_2 * t2_3);
 158         double t4_6 = 2 * (t2_2 * t2_4) + t2_3 * t2_3;
 159         double t4_7 = 2 * (t2_2 * t2_5 + t2_3 * t2_4);
 160         double t4_8 = 2 * (t2_2 * t2_6 + t2_3 * t2_5) + t2_4 * t2_4;
 161         double t5_6 = t4_4 * t1_2 + t4_5 * t1_1;
 162         double t5_8 = t4_4 * t1_4 + t4_5 * t1_3 + t4_6 * t1_2 + t4_7 * t1_1;
 163         double t6_6 = t4_4 * t2_2;
 164         double t6_7 = t4_4 * t2_3 + t4_5 * t2_2;
 165         double t6_8 = t4_4 * t2_4 + t4_5 * t2_3 + t4_6 * t2_2;
 166         double t7_8 = t6_6 * t1_2 + t6_7 * t1_1;
 167         double t8_8 = t6_6 * t2_2;
 168         u = 1;
 169         v = 0;
 170         v += (1./12) * t1_2 + (1./80) * t1_4;
 171         u -= (1./24) * t2_2 + (1./160) * t2_4 + (1./896) * t2_6 + (1./4608) * t2_8;
 172         v -= (1./480) * t3_4 + (1./2688) * t3_6 + (1./13824) * t3_8;
 173         u += (1./1920) * t4_4 + (1./10752) * t4_6 + (1./55296) * t4_8;
 174         v += (1./53760) * t5_6 + (1./276480) * t5_8;
 175         u -= (1./322560) * t6_6 + (1./1.65888e+06) * t6_8;
 176         v -= (1./1.16122e+07) * t7_8;
 177         u += (1./9.28973e+07) * t8_8;
 178 #endif
 179
 180 #if ORDER == 12
 181         double t1_1 = km0;
 182         double t1_2 = .5 * km1;
 183         double t1_3 = (1./6) * km2;
 184         double t1_4 = (1./24) * km3;
 185         double t2_2 = t1_1 * t1_1;
 186         double t2_3 = 2 * (t1_1 * t1_2);
 187         double t2_4 = 2 * (t1_1 * t1_3) + t1_2 * t1_2;
 188         double t2_5 = 2 * (t1_1 * t1_4 + t1_2 * t1_3);
 189         double t2_6 = 2 * (t1_2 * t1_4) + t1_3 * t1_3;
 190         double t2_7 = 2 * (t1_3 * t1_4);
 191         double t2_8 = t1_4 * t1_4;
 192         double t3_4 = t2_2 * t1_2 + t2_3 * t1_1;
 193         double t3_6 = t2_2 * t1_4 + t2_3 * t1_3 + t2_4 * t1_2 + t2_5 * t1_1;
 194         double t3_8 = t2_4 * t1_4 + t2_5 * t1_3 + t2_6 * t1_2 + t2_7 * t1_1;
 195         double t3_10 = t2_6 * t1_4 + t2_7 * t1_3 + t2_8 * t1_2;
 196         double t4_4 = t2_2 * t2_2;
 197         double t4_5 = 2 * (t2_2 * t2_3);
 198         double t4_6 = 2 * (t2_2 * t2_4) + t2_3 * t2_3;
 199         double t4_7 = 2 * (t2_2 * t2_5 + t2_3 * t2_4);
 200         double t4_8 = 2 * (t2_2 * t2_6 + t2_3 * t2_5) + t2_4 * t2_4;
 201         double t4_9 = 2 * (t2_2 * t2_7 + t2_3 * t2_6 + t2_4 * t2_5);
 202         double t4_10 = 2 * (t2_2 * t2_8 + t2_3 * t2_7 + t2_4 * t2_6) + t2_5 * t2_5;
 203         double t5_6 = t4_4 * t1_2 + t4_5 * t1_1;
 204         double t5_8 = t4_4 * t1_4 + t4_5 * t1_3 + t4_6 * t1_2 + t4_7 * t1_1;
 205         double t5_10 = t4_6 * t1_4 + t4_7 * t1_3 + t4_8 * t1_2 + t4_9 * t1_1;
 206         double t6_6 = t4_4 * t2_2;
 207         double t6_7 = t4_4 * t2_3 + t4_5 * t2_2;
 208         double t6_8 = t4_4 * t2_4 + t4_5 * t2_3 + t4_6 * t2_2;
 209         double t6_9 = t4_4 * t2_5 + t4_5 * t2_4 + t4_6 * t2_3 + t4_7 * t2_2;
 210         double t6_10 = t4_4 * t2_6 + t4_5 * t2_5 + t4_6 * t2_4 + t4_7 * t2_3 + t4_8 * t2_2;
 211         double t7_8 = t6_6 * t1_2 + t6_7 * t1_1;
 212         double t7_10 = t6_6 * t1_4 + t6_7 * t1_3 + t6_8 * t1_2 + t6_9 * t1_1;
 213         double t8_8 = t6_6 * t2_2;
 214         double t8_9 = t6_6 * t2_3 + t6_7 * t2_2;
 215         double t8_10 = t6_6 * t2_4 + t6_7 * t2_3 + t6_8 * t2_2;
 216         double t9_10 = t8_8 * t1_2 + t8_9 * t1_1;
 217         double t10_10 = t8_8 * t2_2;
 218         u = 1;
 219         v = 0;
 220         v += (1./12) * t1_2 + (1./80) * t1_4;
 221         u -= (1./24) * t2_2 + (1./160) * t2_4 + (1./896) * t2_6 + (1./4608) * t2_8;
 222         v -= (1./480) * t3_4 + (1./2688) * t3_6 + (1./13824) * t3_8 + (1./67584) * t3_10;
 223         u += (1./1920) * t4_4 + (1./10752) * t4_6 + (1./55296) * t4_8 + (1./270336) * t4_10;
 224         v += (1./53760) * t5_6 + (1./276480) * t5_8 + (1./1.35168e+06) * t5_10;
 225         u -= (1./322560) * t6_6 + (1./1.65888e+06) * t6_8 + (1./8.11008e+06) * t6_10;
 226         v -= (1./1.16122e+07) * t7_8 + (1./5.67706e+07) * t7_10;
 227         u += (1./9.28973e+07) * t8_8 + (1./4.54164e+08) * t8_10;
 228         v += (1./4.08748e+09) * t9_10;
 229         u -= (1./4.08748e+10) * t10_10;
 230 #endif
 231
 232 #if ORDER == 14
 233         double t1_1 = km0;
 234         double t1_2 = .5 * km1;
 235         double t1_3 = (1./6) * km2;
 236         double t1_4 = (1./24) * km3;
 237         double t2_2 = t1_1 * t1_1;
 238         double t2_3 = 2 * (t1_1 * t1_2);
 239         double t2_4 = 2 * (t1_1 * t1_3) + t1_2 * t1_2;
 240         double t2_5 = 2 * (t1_1 * t1_4 + t1_2 * t1_3);
 241         double t2_6 = 2 * (t1_2 * t1_4) + t1_3 * t1_3;
 242         double t2_7 = 2 * (t1_3 * t1_4);
 243         double t2_8 = t1_4 * t1_4;
 244         double t3_4 = t2_2 * t1_2 + t2_3 * t1_1;
 245         double t3_6 = t2_2 * t1_4 + t2_3 * t1_3 + t2_4 * t1_2 + t2_5 * t1_1;
 246         double t3_8 = t2_4 * t1_4 + t2_5 * t1_3 + t2_6 * t1_2 + t2_7 * t1_1;
 247         double t3_10 = t2_6 * t1_4 + t2_7 * t1_3 + t2_8 * t1_2;
 248         double t3_12 = t2_8 * t1_4;
 249         double t4_4 = t2_2 * t2_2;
 250         double t4_5 = 2 * (t2_2 * t2_3);
 251         double t4_6 = 2 * (t2_2 * t2_4) + t2_3 * t2_3;
 252         double t4_7 = 2 * (t2_2 * t2_5 + t2_3 * t2_4);
 253         double t4_8 = 2 * (t2_2 * t2_6 + t2_3 * t2_5) + t2_4 * t2_4;
 254         double t4_9 = 2 * (t2_2 * t2_7 + t2_3 * t2_6 + t2_4 * t2_5);
 255         double t4_10 = 2 * (t2_2 * t2_8 + t2_3 * t2_7 + t2_4 * t2_6) + t2_5 * t2_5;
 256         double t4_11 = 2 * (t2_3 * t2_8 + t2_4 * t2_7 + t2_5 * t2_6);
 257         double t4_12 = 2 * (t2_4 * t2_8 + t2_5 * t2_7) + t2_6 * t2_6;
 258         double t5_6 = t4_4 * t1_2 + t4_5 * t1_1;
 259         double t5_8 = t4_4 * t1_4 + t4_5 * t1_3 + t4_6 * t1_2 + t4_7 * t1_1;
 260         double t5_10 = t4_6 * t1_4 + t4_7 * t1_3 + t4_8 * t1_2 + t4_9 * t1_1;
 261         double t5_12 = t4_8 * t1_4 + t4_9 * t1_3 + t4_10 * t1_2 + t4_11 * t1_1;
 262         double t6_6 = t4_4 * t2_2;
 263         double t6_7 = t4_4 * t2_3 + t4_5 * t2_2;
 264         double t6_8 = t4_4 * t2_4 + t4_5 * t2_3 + t4_6 * t2_2;
 265         double t6_9 = t4_4 * t2_5 + t4_5 * t2_4 + t4_6 * t2_3 + t4_7 * t2_2;
 266         double t6_10 = t4_4 * t2_6 + t4_5 * t2_5 + t4_6 * t2_4 + t4_7 * t2_3 + t4_8 * t2_2;
 267         double t6_11 = t4_4 * t2_7 + t4_5 * t2_6 + t4_6 * t2_5 + t4_7 * t2_4 + t4_8 * t2_3 + t4_9 * t2_2;
 268         double t6_12 = t4_4 * t2_8 + t4_5 * t2_7 + t4_6 * t2_6 + t4_7 * t2_5 + t4_8 * t2_4 + t4_9 * t2_3 + t4_10 * t2_2;
 269         double t7_8 = t6_6 * t1_2 + t6_7 * t1_1;
 270         double t7_10 = t6_6 * t1_4 + t6_7 * t1_3 + t6_8 * t1_2 + t6_9 * t1_1;
 271         double t7_12 = t6_8 * t1_4 + t6_9 * t1_3 + t6_10 * t1_2 + t6_11 * t1_1;
 272         double t8_8 = t6_6 * t2_2;
 273         double t8_9 = t6_6 * t2_3 + t6_7 * t2_2;
 274         double t8_10 = t6_6 * t2_4 + t6_7 * t2_3 + t6_8 * t2_2;
 275         double t8_11 = t6_6 * t2_5 + t6_7 * t2_4 + t6_8 * t2_3 + t6_9 * t2_2;
 276         double t8_12 = t6_6 * t2_6 + t6_7 * t2_5 + t6_8 * t2_4 + t6_9 * t2_3 + t6_10 * t2_2;
 277         double t9_10 = t8_8 * t1_2 + t8_9 * t1_1;
 278         double t9_12 = t8_8 * t1_4 + t8_9 * t1_3 + t8_10 * t1_2 + t8_11 * t1_1;
 279         double t10_10 = t8_8 * t2_2;
 280         double t10_11 = t8_8 * t2_3 + t8_9 * t2_2;
 281         double t10_12 = t8_8 * t2_4 + t8_9 * t2_3 + t8_10 * t2_2;
 282         double t11_12 = t10_10 * t1_2 + t10_11 * t1_1;
 283         double t12_12 = t10_10 * t2_2;
 284         u = 1;
 285         v = 0;
 286         v += (1./12) * t1_2 + (1./80) * t1_4;
 287         u -= (1./24) * t2_2 + (1./160) * t2_4 + (1./896) * t2_6 + (1./4608) * t2_8;
 288         v -= (1./480) * t3_4 + (1./2688) * t3_6 + (1./13824) * t3_8 + (1./67584) * t3_10 + (1./319488) * t3_12;
 289         u += (1./1920) * t4_4 + (1./10752) * t4_6 + (1./55296) * t4_8 + (1./270336) * t4_10 + (1./1.27795e+06) * t4_12;
 290         v += (1./53760) * t5_6 + (1./276480) * t5_8 + (1./1.35168e+06) * t5_10 + (1./6.38976e+06) * t5_12;
 291         u -= (1./322560) * t6_6 + (1./1.65888e+06) * t6_8 + (1./8.11008e+06) * t6_10 + (1./3.83386e+07) * t6_12;
 292         v -= (1./1.16122e+07) * t7_8 + (1./5.67706e+07) * t7_10 + (1./2.6837e+08) * t7_12;
 293         u += (1./9.28973e+07) * t8_8 + (1./4.54164e+08) * t8_10 + (1./2.14696e+09) * t8_12;
 294         v += (1./4.08748e+09) * t9_10 + (1./1.93226e+10) * t9_12;
 295         u -= (1./4.08748e+10) * t10_10 + (1./1.93226e+11) * t10_12;
 296         v -= (1./2.12549e+12) * t11_12;
 297         u += (1./2.55059e+13) * t12_12;
 298 #endif
 299
 300 #if ORDER == 16
 301         double t1_1 = km0;
 302         double t1_2 = .5 * km1;
 303         double t1_3 = (1./6) * km2;
 304         double t1_4 = (1./24) * km3;
 305         double t2_2 = t1_1 * t1_1;
 306         double t2_3 = 2 * (t1_1 * t1_2);
 307         double t2_4 = 2 * (t1_1 * t1_3) + t1_2 * t1_2;
 308         double t2_5 = 2 * (t1_1 * t1_4 + t1_2 * t1_3);
 309         double t2_6 = 2 * (t1_2 * t1_4) + t1_3 * t1_3;
 310         double t2_7 = 2 * (t1_3 * t1_4);
 311         double t2_8 = t1_4 * t1_4;
 312         double t3_4 = t2_2 * t1_2 + t2_3 * t1_1;
 313         double t3_6 = t2_2 * t1_4 + t2_3 * t1_3 + t2_4 * t1_2 + t2_5 * t1_1;
 314         double t3_8 = t2_4 * t1_4 + t2_5 * t1_3 + t2_6 * t1_2 + t2_7 * t1_1;
 315         double t3_10 = t2_6 * t1_4 + t2_7 * t1_3 + t2_8 * t1_2;
 316         double t3_12 = t2_8 * t1_4;
 317         double t4_4 = t2_2 * t2_2;
 318         double t4_5 = 2 * (t2_2 * t2_3);
 319         double t4_6 = 2 * (t2_2 * t2_4) + t2_3 * t2_3;
 320         double t4_7 = 2 * (t2_2 * t2_5 + t2_3 * t2_4);
 321         double t4_8 = 2 * (t2_2 * t2_6 + t2_3 * t2_5) + t2_4 * t2_4;
 322         double t4_9 = 2 * (t2_2 * t2_7 + t2_3 * t2_6 + t2_4 * t2_5);
 323         double t4_10 = 2 * (t2_2 * t2_8 + t2_3 * t2_7 + t2_4 * t2_6) + t2_5 * t2_5;
 324         double t4_11 = 2 * (t2_3 * t2_8 + t2_4 * t2_7 + t2_5 * t2_6);
 325         double t4_12 = 2 * (t2_4 * t2_8 + t2_5 * t2_7) + t2_6 * t2_6;
 326         double t4_13 = 2 * (t2_5 * t2_8 + t2_6 * t2_7);
 327         double t4_14 = 2 * (t2_6 * t2_8) + t2_7 * t2_7;
 328         double t5_6 = t4_4 * t1_2 + t4_5 * t1_1;
 329         double t5_8 = t4_4 * t1_4 + t4_5 * t1_3 + t4_6 * t1_2 + t4_7 * t1_1;
 330         double t5_10 = t4_6 * t1_4 + t4_7 * t1_3 + t4_8 * t1_2 + t4_9 * t1_1;
 331         double t5_12 = t4_8 * t1_4 + t4_9 * t1_3 + t4_10 * t1_2 + t4_11 * t1_1;
 332         double t5_14 = t4_10 * t1_4 + t4_11 * t1_3 + t4_12 * t1_2 + t4_13 * t1_1;
 333         double t6_6 = t4_4 * t2_2;
 334         double t6_7 = t4_4 * t2_3 + t4_5 * t2_2;
 335         double t6_8 = t4_4 * t2_4 + t4_5 * t2_3 + t4_6 * t2_2;
 336         double t6_9 = t4_4 * t2_5 + t4_5 * t2_4 + t4_6 * t2_3 + t4_7 * t2_2;
 337         double t6_10 = t4_4 * t2_6 + t4_5 * t2_5 + t4_6 * t2_4 + t4_7 * t2_3 + t4_8 * t2_2;
 338         double t6_11 = t4_4 * t2_7 + t4_5 * t2_6 + t4_6 * t2_5 + t4_7 * t2_4 + t4_8 * t2_3 + t4_9 * t2_2;
 339         double t6_12 = t4_4 * t2_8 + t4_5 * t2_7 + t4_6 * t2_6 + t4_7 * t2_5 + t4_8 * t2_4 + t4_9 * t2_3 + t4_10 * t2_2;
 340         double t6_13 = t4_5 * t2_8 + t4_6 * t2_7 + t4_7 * t2_6 + t4_8 * t2_5 + t4_9 * t2_4 + t4_10 * t2_3 + t4_11 * t2_2;
 341         double t6_14 = t4_6 * t2_8 + t4_7 * t2_7 + t4_8 * t2_6 + t4_9 * t2_5 + t4_10 * t2_4 + t4_11 * t2_3 + t4_12 * t2_2;
 342         double t7_8 = t6_6 * t1_2 + t6_7 * t1_1;
 343         double t7_10 = t6_6 * t1_4 + t6_7 * t1_3 + t6_8 * t1_2 + t6_9 * t1_1;
 344         double t7_12 = t6_8 * t1_4 + t6_9 * t1_3 + t6_10 * t1_2 + t6_11 * t1_1;
 345         double t7_14 = t6_10 * t1_4 + t6_11 * t1_3 + t6_12 * t1_2 + t6_13 * t1_1;
 346         double t8_8 = t6_6 * t2_2;
 347         double t8_9 = t6_6 * t2_3 + t6_7 * t2_2;
 348         double t8_10 = t6_6 * t2_4 + t6_7 * t2_3 + t6_8 * t2_2;
 349         double t8_11 = t6_6 * t2_5 + t6_7 * t2_4 + t6_8 * t2_3 + t6_9 * t2_2;
 350         double t8_12 = t6_6 * t2_6 + t6_7 * t2_5 + t6_8 * t2_4 + t6_9 * t2_3 + t6_10 * t2_2;
 351         double t8_13 = t6_6 * t2_7 + t6_7 * t2_6 + t6_8 * t2_5 + t6_9 * t2_4 + t6_10 * t2_3 + t6_11 * t2_2;
 352         double t8_14 = t6_6 * t2_8 + t6_7 * t2_7 + t6_8 * t2_6 + t6_9 * t2_5 + t6_10 * t2_4 + t6_11 * t2_3 + t6_12 * t2_2;
 353         double t9_10 = t8_8 * t1_2 + t8_9 * t1_1;
 354         double t9_12 = t8_8 * t1_4 + t8_9 * t1_3 + t8_10 * t1_2 + t8_11 * t1_1;
 355         double t9_14 = t8_10 * t1_4 + t8_11 * t1_3 + t8_12 * t1_2 + t8_13 * t1_1;
 356         double t10_10 = t8_8 * t2_2;
 357         double t10_11 = t8_8 * t2_3 + t8_9 * t2_2;
 358         double t10_12 = t8_8 * t2_4 + t8_9 * t2_3 + t8_10 * t2_2;
 359         double t10_13 = t8_8 * t2_5 + t8_9 * t2_4 + t8_10 * t2_3 + t8_11 * t2_2;
 360         double t10_14 = t8_8 * t2_6 + t8_9 * t2_5 + t8_10 * t2_4 + t8_11 * t2_3 + t8_12 * t2_2;
 361         double t11_12 = t10_10 * t1_2 + t10_11 * t1_1;
 362         double t11_14 = t10_10 * t1_4 + t10_11 * t1_3 + t10_12 * t1_2 + t10_13 * t1_1;
 363         double t12_12 = t10_10 * t2_2;
 364         double t12_13 = t10_10 * t2_3 + t10_11 * t2_2;
 365         double t12_14 = t10_10 * t2_4 + t10_11 * t2_3 + t10_12 * t2_2;
 366         double t13_14 = t12_12 * t1_2 + t12_13 * t1_1;
 367         double t14_14 = t12_12 * t2_2;
 368         u = 1;
 369         u -= 1./24 * t2_2 + 1./160 * t2_4 + 1./896 * t2_6 + 1./4608 * t2_8;
 370         u += 1./1920 * t4_4 + 1./10752 * t4_6 + 1./55296 * t4_8 + 1./270336 * t4_10 + 1./1277952 * t4_12 + 1./5898240 * t4_14;
 371         u -= 1./322560 * t6_6 + 1./1658880 * t6_8 + 1./8110080 * t6_10 + 1./38338560 * t6_12 + 1./176947200 * t6_14;
 372         u += 1./92897280 * t8_8 + 1./454164480 * t8_10 + 4.6577500191e-10 * t8_12 + 1.0091791708e-10 * t8_14;
 373         u -= 2.4464949595e-11 * t10_10 + 5.1752777990e-12 * t10_12 + 1.1213101898e-12 * t10_14;
 374         u += 3.9206649992e-14 * t12_12 + 8.4947741650e-15 * t12_14;
 375         u -= 4.6674583324e-17 * t14_14;
 376         v = 0;
 377         v += 1./12 * t1_2 + 1./80 * t1_4;
 378         v -= 1./480 * t3_4 + 1./2688 * t3_6 + 1./13824 * t3_8 + 1./67584 * t3_10 + 1./319488 * t3_12;
 379         v += 1./53760 * t5_6 + 1./276480 * t5_8 + 1./1351680 * t5_10 + 1./6389760 * t5_12 + 1./29491200 * t5_14;
 380         v -= 1./11612160 * t7_8 + 1./56770560 * t7_10 + 1./268369920 * t7_12 + 8.0734333664e-10 * t7_14;
 381         v += 2.4464949595e-10 * t9_10 + 5.1752777990e-11 * t9_12 + 1.1213101898e-11 * t9_14;
 382         v -= 4.7047979991e-13 * t11_12 + 1.0193728998e-13 * t11_14;
 383         v += 6.5344416654e-16 * t13_14;
 384 #endif
 385
 386         }
 387
 388         if (n == 1) {
 389 #if ORDER == 2
 390             x = 1;
 391             y = 0;
 392 #else
 393             x = u;
 394             y = v;
 395 #endif
 396         } else {
 397             double th = (((th4 * s + th3) * s + th2) * s + th1) * s;
 398             double cth = cos(th);
 399             double sth = sin(th);
 400
 401 #if ORDER == 2
 402             x += cth;
 403             y += sth;
 404 #else
 405             x += cth * u - sth * v;
 406             y += cth * v + sth * u;
 407 #endif
 408             s += ds;
 409         }
 410     }
 411
 412 #if ORDER == 4 || ORDER == 6
 413     xy[0] = x * (1./24 * ds);
 414     xy[1] = y * (1./24 * ds);
 415 #else
 416     xy[0] = x * ds;
 417     xy[1] = y * ds;
 418 #endif
 419 }
 420
 421 static double
 422 compute_ends(const double ks[4], double ends[2][4], double seg_ch)
 423 {
 424     double xy[2];
 425     double ch, th;
 426     double l, l2, l3;
 427     double th_even, th_odd;
 428     double k0_even, k0_odd;
 429     double k1_even, k1_odd;
 430     double k2_even, k2_odd;
 431
 432     integrate_spiro(ks, xy);
 433     ch = hypot(xy[0], xy[1]);
 434     th = atan2(xy[1], xy[0]);
 435     l = ch / seg_ch;
 436
 437     th_even = .5 * ks[0] + (1./48) * ks[2];
 438     th_odd = .125 * ks[1] + (1./384) * ks[3] - th;
 439     ends[0][0] = th_even - th_odd;
 440     ends[1][0] = th_even + th_odd;
 441     k0_even = l * (ks[0] + .125 * ks[2]);
 442     k0_odd = l * (.5 * ks[1] + (1./48) * ks[3]);
 443     ends[0][1] = k0_even - k0_odd;
 444     ends[1][1] = k0_even + k0_odd;
 445     l2 = l * l;
 446     k1_even = l2 * (ks[1] + .125 * ks[3]);
 447     k1_odd = l2 * .5 * ks[2];
 448     ends[0][2] = k1_even - k1_odd;
 449     ends[1][2] = k1_even + k1_odd;
 450     l3 = l2 * l;
 451     k2_even = l3 * ks[2];
 452     k2_odd = l3 * .5 * ks[3];
 453     ends[0][3] = k2_even - k2_odd;
 454     ends[1][3] = k2_even + k2_odd;
 455
 456     return l;
 457 }
 458
 459 static void
 460 compute_pderivs(const spiro_seg *s, double ends[2][4], double derivs[4][2][4],
 461                 int jinc)
 462 {
 463     double recip_d = 2e6;
 464     double delta = 1./ recip_d;
 465     double try_ks[4];
 466     double try_ends[2][4];
 467     int i, j, k;
 468
 469     compute_ends(s->ks, ends, s->seg_ch);
 470     for (i = 0; i < jinc; i++) {
 471         for (j = 0; j < 4; j++)
 472             try_ks[j] = s->ks[j];
 473         try_ks[i] += delta;
 474         compute_ends(try_ks, try_ends, s->seg_ch);
 475         for (k = 0; k < 2; k++)
 476             for (j = 0; j < 4; j++)
 477                 derivs[j][k][i] = recip_d * (try_ends[k][j] - ends[k][j]);
 478     }
 479 }
 480
 481 static double
 482 mod_2pi(double th)
 483 {
 484     double u = th / (2 * M_PI);
 485     return 2 * M_PI * (u - floor(u + 0.5));
 486 }
 487
 488 static spiro_seg *
 489 setup_path(const spiro_cp *src, int n)
 490 {
 491     int n_seg = src[0].ty == '{' ? n - 1 : n;
 492     spiro_seg *r = (spiro_seg *)malloc((n_seg + 1) * sizeof(spiro_seg));
 493     int i;
 494     int ilast;
 495
 496     for (i = 0; i < n_seg; i++) {
 497         r[i].x = src[i].x;
 498         r[i].y = src[i].y;
 499         r[i].ty = src[i].ty;
 500         r[i].ks[0] = 0.;
 501         r[i].ks[1] = 0.;
 502         r[i].ks[2] = 0.;
 503         r[i].ks[3] = 0.;
 504     }
 505     r[n_seg].x = src[n_seg % n].x;
 506     r[n_seg].y = src[n_seg % n].y;
 507     r[n_seg].ty = src[n_seg % n].ty;
 508
 509     for (i = 0; i < n_seg; i++) {
 510         double dx = r[i + 1].x - r[i].x;
 511         double dy = r[i + 1].y - r[i].y;
 512         r[i].seg_ch = hypot(dx, dy);
 513         r[i].seg_th = atan2(dy, dx);
 514     }
 515
 516     ilast = n_seg - 1;
 517     for (i = 0; i < n_seg; i++) {
 518         if (r[i].ty == '{' || r[i].ty == '}' || r[i].ty == 'v')
 519             r[i].bend_th = 0.;
 520         else
 521             r[i].bend_th = mod_2pi(r[i].seg_th - r[ilast].seg_th);
 522         ilast = i;
 523     }
 524     return r;
 525 }
 526
 527 static void
 528 bandec11(bandmat *m, int *perm, int n)
 529 {
 530     int i, j, k;
 531     int l;
 532
 533     /* pack top triangle to the left. */
 534     for (i = 0; i < 5; i++) {
 535         for (j = 0; j < i + 6; j++)
 536             m[i].a[j] = m[i].a[j + 5 - i];
 537         for (; j < 11; j++)
 538             m[i].a[j] = 0.;
 539     }
 540     l = 5;
 541     for (k = 0; k < n; k++) {
 542         int pivot = k;
 543         double pivot_val = m[k].a[0];
 544         double pivot_scale;
 545
 546         l = l < n ? l + 1 : n;
 547
 548         for (j = k + 1; j < l; j++)
 549             if (fabs(m[j].a[0]) > fabs(pivot_val)) {
 550                 pivot_val = m[j].a[0];
 551                 pivot = j;
 552             }
 553
 554         perm[k] = pivot;
 555         if (pivot != k) {
 556             for (j = 0; j < 11; j++) {
 557                 double tmp = m[k].a[j];
 558                 m[k].a[j] = m[pivot].a[j];
 559                 m[pivot].a[j] = tmp;
 560             }
 561         }
 562
 563         if (fabs(pivot_val) < 1e-12) pivot_val = 1e-12;
 564         pivot_scale = 1. / pivot_val;
 565         for (i = k + 1; i < l; i++) {
 566             double x = m[i].a[0] * pivot_scale;
 567             m[k].al[i - k - 1] = x;
 568             for (j = 1; j < 11; j++)
 569                 m[i].a[j - 1] = m[i].a[j] - x * m[k].a[j];
 570             m[i].a[10] = 0.;
 571         }
 572     }
 573 }
 574
 575 static void
 576 banbks11(const bandmat *m, const int *perm, double *v, int n)
 577 {
 578     int i, k, l;
 579
 580     /* forward substitution */
 581     l = 5;
 582     for (k = 0; k < n; k++) {
 583         i = perm[k];
 584         if (i != k) {
 585             double tmp = v[k];
 586             v[k] = v[i];
 587             v[i] = tmp;
 588         }
 589         if (l < n) l++;
 590         for (i = k + 1; i < l; i++)
 591             v[i] -= m[k].al[i - k - 1] * v[k];
 592     }
 593
 594     /* back substitution */
 595     l = 1;
 596     for (i = n - 1; i >= 0; i--) {
 597         double x = v[i];
 598         for (k = 1; k < l; k++)
 599             x -= m[i].a[k] * v[k + i];
 600         v[i] = x / m[i].a[0];
 601         if (l < 11) l++;
 602     }
 603 }
 604
 605 int compute_jinc(char ty0, char ty1)
 606 {
 607     if (ty0 == 'o' || ty1 == 'o' ||
 608         ty0 == ']' || ty1 == '[')
 609         return 4;
 610     else if (ty0 == 'c' && ty1 == 'c')
 611         return 2;
 612     else if (((ty0 == '{' || ty0 == 'v' || ty0 == '[') && ty1 == 'c') ||
 613              (ty0 == 'c' && (ty1 == '}' || ty1 == 'v' || ty1 == ']')))
 614         return 1;
 615     else
 616         return 0;
 617 }
 618
 619 int count_vec(const spiro_seg *s, int nseg)
 620 {
 621     int i;
 622     int n = 0;
 623
 624     for (i = 0; i < nseg; i++)
 625         n += compute_jinc(s[i].ty, s[i + 1].ty);
 626     return n;
 627 }
 628
 629 static void
 630 add_mat_line(bandmat *m, double *v,
 631              double derivs[4], double x, double y, int j, int jj, int jinc,
 632              int nmat)
 633 {
 634     int k;
 635
 636     if (jj >= 0) {
 637         int joff =  (j + 5 - jj + nmat) % nmat;
 638         if (nmat < 6) {
 639             joff = j + 5 - jj;
 640         } else if (nmat == 6) {
 641             joff = 2 + (j + 3 - jj + nmat) % nmat;
 642         }
 643 #ifdef VERBOSE
 644         printf("add_mat_line j=%d jj=%d jinc=%d nmat=%d joff=%d\n", j, jj, jinc, nmat, joff);
 645 #endif
 646         v[jj] += x;
 647         for (k = 0; k < jinc; k++)
 648             m[jj].a[joff + k] += y * derivs[k];
 649     }
 650 }
 651
 652 static double
 653 spiro_iter(spiro_seg *s, bandmat *m, int *perm, double *v, int n)
 654 {
 655     int cyclic = s[0].ty != '{' && s[0].ty != 'v';
 656     int i, j, jj;
 657     int nmat = count_vec(s, n);
 658     double norm;
 659     int n_invert;
 660
 661     for (i = 0; i < nmat; i++) {
 662         v[i] = 0.;
 663         for (j = 0; j < 11; j++)
 664             m[i].a[j] = 0.;
 665         for (j = 0; j < 5; j++)
 666             m[i].al[j] = 0.;
 667     }
 668
 669     j = 0;
 670     if (s[0].ty == 'o')
 671         jj = nmat - 2;
 672     else if (s[0].ty == 'c')
 673         jj = nmat - 1;
 674     else
 675         jj = 0;
 676     for (i = 0; i < n; i++) {
 677         char ty0 = s[i].ty;
 678         char ty1 = s[i + 1].ty;
 679         int jinc = compute_jinc(ty0, ty1);
 680         double th = s[i].bend_th;
 681         double ends[2][4];
 682         double derivs[4][2][4];
 683         int jthl = -1, jk0l = -1, jk1l = -1, jk2l = -1;
 684         int jthr = -1, jk0r = -1, jk1r = -1, jk2r = -1;
 685
 686         compute_pderivs(&s[i], ends, derivs, jinc);
 687
 688         /* constraints crossing left */
 689         if (ty0 == 'o' || ty0 == 'c' || ty0 == '[' || ty0 == ']') {
 690             jthl = jj++;
 691             jj %= nmat;
 692             jk0l = jj++;
 693         }
 694         if (ty0 == 'o') {
 695             jj %= nmat;
 696             jk1l = jj++;
 697             jk2l = jj++;
 698         }
 699
 700         /* constraints on left */
 701         if ((ty0 == '[' || ty0 == 'v' || ty0 == '{' || ty0 == 'c') &&
 702             jinc == 4) {
 703             if (ty0 != 'c')
 704                 jk1l = jj++;
 705             jk2l = jj++;
 706         }
 707
 708         /* constraints on right */
 709         if ((ty1 == ']' || ty1 == 'v' || ty1 == '}' || ty1 == 'c') &&
 710             jinc == 4) {
 711             if (ty1 != 'c')
 712                 jk1r = jj++;
 713             jk2r = jj++;
 714         }
 715
 716         /* constraints crossing right */
 717         if (ty1 == 'o' || ty1 == 'c' || ty1 == '[' || ty1 == ']') {
 718             jthr = jj;
 719             jk0r = (jj + 1) % nmat;
 720         }
 721         if (ty1 == 'o') {
 722             jk1r = (jj + 2) % nmat;
 723             jk2r = (jj + 3) % nmat;
 724         }
 725
 726         add_mat_line(m, v, derivs[0][0], th - ends[0][0], 1, j, jthl, jinc, nmat);
 727         add_mat_line(m, v, derivs[1][0], ends[0][1], -1, j, jk0l, jinc, nmat);
 728         add_mat_line(m, v, derivs[2][0], ends[0][2], -1, j, jk1l, jinc, nmat);
 729         add_mat_line(m, v, derivs[3][0], ends[0][3], -1, j, jk2l, jinc, nmat);
 730         add_mat_line(m, v, derivs[0][1], -ends[1][0], 1, j, jthr, jinc, nmat);
 731         add_mat_line(m, v, derivs[1][1], -ends[1][1], 1, j, jk0r, jinc, nmat);
 732         add_mat_line(m, v, derivs[2][1], -ends[1][2], 1, j, jk1r, jinc, nmat);
 733         add_mat_line(m, v, derivs[3][1], -ends[1][3], 1, j, jk2r, jinc, nmat);
 734         if (jthl >= 0)
 735                 v[jthl] = mod_2pi(v[jthl]);
 736         if (jthr >= 0)
 737                 v[jthr] = mod_2pi(v[jthr]);
 738         j += jinc;
 739     }
 740     if (cyclic) {
 741         memcpy(m + nmat, m, sizeof(bandmat) * nmat);
 742         memcpy(m + 2 * nmat, m, sizeof(bandmat) * nmat);
 743         memcpy(v + nmat, v, sizeof(double) * nmat);
 744         memcpy(v + 2 * nmat, v, sizeof(double) * nmat);
 745         n_invert = 3 * nmat;
 746         j = nmat;
 747     } else {
 748         n_invert = nmat;
 749         j = 0;
 750     }
 751 #ifdef VERBOSE
 752     for (i = 0; i < n; i++) {
 753         int k;
 754         for (k = 0; k < 11; k++)
 755             printf(" %2.4f", m[i].a[k]);
 756         printf(": %2.4f\n", v[i]);
 757     }
 758     printf("---\n");
 759 #endif
 760     bandec11(m, perm, n_invert);
 761     banbks11(m, perm, v, n_invert);
 762     norm = 0.;
 763     for (i = 0; i < n; i++) {
 764         char ty0 = s[i].ty;
 765         char ty1 = s[i + 1].ty;
 766         int jinc = compute_jinc(ty0, ty1);
 767         int k;
 768
 769         for (k = 0; k < jinc; k++) {
 770             double dk = v[j++];
 771
 772 #ifdef VERBOSE
 773             printf("s[%d].ks[%d] += %f\n", i, k, dk);
 774 #endif
 775             s[i].ks[k] += dk;
 776             norm += dk * dk;
 777         }
 778         s[i].ks[0] = 2.0*mod_2pi(s[i].ks[0]/2.0);
 779     }
 780     return norm;
 781 }
 782
 783 int
 784 solve_spiro(spiro_seg *s, int nseg)
 785 {
 786     bandmat *m;
 787     double *v;
 788     int *perm;
 789     int nmat = count_vec(s, nseg);
 790     int n_alloc = nmat;
 791     double norm;
 792     int i;
 793
 794     if (nmat == 0)
 795         return 0;
 796     if (s[0].ty != '{' && s[0].ty != 'v')
 797         n_alloc *= 3;
 798     if (n_alloc < 5)
 799         n_alloc = 5;
 800     m = (bandmat *)malloc(sizeof(bandmat) * n_alloc);
 801     v = (double *)malloc(sizeof(double) * n_alloc);
 802     perm = (int *)malloc(sizeof(int) * n_alloc);
 803
 804     for (i = 0; i < 10; i++) {
 805         norm = spiro_iter(s, m, perm, v, nseg);
 806 #ifdef VERBOSE
 807         printf("%% norm = %g\n", norm);
 808 #endif
 809         if (norm < 1e-12) break;
 810     }
 811
 812     free(m);
 813     free(v);
 814     free(perm);
 815     return 0;
 816 }
 817
 818 static void
 819 spiro_seg_to_bpath(const double ks[4],
 820                    double x0, double y0, double x1, double y1,
 821                    bezctx *bc, int depth)
 822 {
 823     double bend = fabs(ks[0]) + fabs(.5 * ks[1]) + fabs(.125 * ks[2]) +
 824         fabs((1./48) * ks[3]);
 825
 826     if (!bend > 1e-8) {
 827         bezctx_lineto(bc, x1, y1);
 828     } else {
 829         double seg_ch = hypot(x1 - x0, y1 - y0);
 830         double seg_th = atan2(y1 - y0, x1 - x0);
 831         double xy[2];
 832         double ch, th;
 833         double scale, rot;
 834         double th_even, th_odd;
 835         double ul, vl;
 836         double ur, vr;
 837
 838         integrate_spiro(ks, xy);
 839         ch = hypot(xy[0], xy[1]);
 840         th = atan2(xy[1], xy[0]);
 841         scale = seg_ch / ch;
 842         rot = seg_th - th;
 843         if (depth > 5 || bend < 1.) {
 844             th_even = (1./384) * ks[3] + (1./8) * ks[1] + rot;
 845             th_odd = (1./48) * ks[2] + .5 * ks[0];
 846             ul = (scale * (1./3)) * cos(th_even - th_odd);
 847             vl = (scale * (1./3)) * sin(th_even - th_odd);
 848             ur = (scale * (1./3)) * cos(th_even + th_odd);
 849             vr = (scale * (1./3)) * sin(th_even + th_odd);
 850             bezctx_curveto(bc, x0 + ul, y0 + vl, x1 - ur, y1 - vr, x1, y1);
 851         } else {
 852             /* subdivide */
 853             double ksub[4];
 854             double thsub;
 855             double xysub[2];
 856             double xmid, ymid;
 857             double cth, sth;
 858
 859             ksub[0] = .5 * ks[0] - .125 * ks[1] + (1./64) * ks[2] - (1./768) * ks[3];
 860             ksub[1] = .25 * ks[1] - (1./16) * ks[2] + (1./128) * ks[3];
 861             ksub[2] = .125 * ks[2] - (1./32) * ks[3];
 862             ksub[3] = (1./16) * ks[3];
 863             thsub = rot - .25 * ks[0] + (1./32) * ks[1] - (1./384) * ks[2] + (1./6144) * ks[3];
 864             cth = .5 * scale * cos(thsub);
 865             sth = .5 * scale * sin(thsub);
 866             integrate_spiro(ksub, xysub);
 867             xmid = x0 + cth * xysub[0] - sth * xysub[1];
 868             ymid = y0 + cth * xysub[1] + sth * xysub[0];
 869             spiro_seg_to_bpath(ksub, x0, y0, xmid, ymid, bc, depth + 1);
 870             ksub[0] += .25 * ks[1] + (1./384) * ks[3];
 871             ksub[1] += .125 * ks[2];
 872             ksub[2] += (1./16) * ks[3];
 873             spiro_seg_to_bpath(ksub, xmid, ymid, x1, y1, bc, depth + 1);
 874         }
 875     }
 876 }
 877
 878 spiro_seg *
 879 run_spiro(const spiro_cp *src, int n)
 880 {
 881     int nseg = src[0].ty == '{' ? n - 1 : n;
 882     spiro_seg *s = setup_path(src, n);
 883     if (nseg > 1)
 884         solve_spiro(s, nseg);
 885     return s;
 886 }
 887
 888 void
 889 free_spiro(spiro_seg *s)
 890 {
 891     free(s);
 892 }
 893
 894 void
 895 spiro_to_bpath(const spiro_seg *s, int n, bezctx *bc)
 896 {
 897     int i;
 898     int nsegs = s[n - 1].ty == '}' ? n - 1 : n;
 899
 900     for (i = 0; i < nsegs; i++) {
 901         double x0 = s[i].x;
 902         double y0 = s[i].y;
 903         double x1 = s[i + 1].x;
 904         double y1 = s[i + 1].y;
 905
 906         if (i == 0)
 907             bezctx_moveto(bc, x0, y0, s[0].ty == '{');
 908         bezctx_mark_knot(bc, i);
 909         spiro_seg_to_bpath(s[i].ks, x0, y0, x1, y1, bc, 0);
 910     }
 911 }
 912
 913 double
 914 get_knot_th(const spiro_seg *s, int i)
 915 {
 916     double ends[2][4];
 917
 918     if (i == 0) {
 919         compute_ends(s[i].ks, ends, s[i].seg_ch);
 920         return s[i].seg_th - ends[0][0];
 921     } else {
 922         compute_ends(s[i - 1].ks, ends, s[i - 1].seg_ch);
 923         return s[i - 1].seg_th + ends[1][0];
 924     }
 925 }
 926
 927 #ifdef UNIT_TEST
 928 #include <stdio.h>
 929 #include <sys/time.h> /* for gettimeofday */
 930
 931 static double
 932 get_time (void)
 933 {
 934   struct timeval tv;
 935   struct timezone tz;
 936
 937   gettimeofday (&tv, &tz);
 938
 939   return tv.tv_sec + 1e-6 * tv.tv_usec;
 940 }
 941
 942 int
 943 test_integ(void) {
 944     double ks[] = {1, 2, 3, 4};
 945     double xy[2];
 946     double xynom[2];
 947     double ch, th;
 948     int i, j;
 949     int nsubdiv;
 950
 951     n = ORDER < 6 ? 4096 : 1024;
 952     integrate_spiro(ks, xynom);
 953     nsubdiv = ORDER < 12 ? 8 : 7;
 954     for (i = 0; i < nsubdiv; i++) {
 955         double st, en;
 956         double err;
 957         int n_iter = (1 << (20 - i));
 958
 959         n = 1 << i;
 960         st = get_time();
 961         for (j = 0; j < n_iter; j++)
 962             integrate_spiro(ks, xy);
 963         en = get_time();
 964         err = hypot(xy[0] - xynom[0], xy[1] - xynom[1]);
 965         printf("%d %d %g %g\n", ORDER, n, (en - st) / n_iter, err);
 966         ch = hypot(xy[0], xy[1]);
 967         th = atan2(xy[1], xy[0]);
 968 #if 0
 969         printf("n = %d: integ(%g %g %g %g) = %g %g, ch = %g, th = %g\n", n,
 970                ks[0], ks[1], ks[2], ks[3], xy[0], xy[1], ch, th);
 971         printf("%d: %g %g\n", n, xy[0] - xynom[0], xy[1] - xynom[1]);
 972 #endif
 973     }
 974     return 0;
 975 }
 976
 977 void
 978 print_seg(const double ks[4], double x0, double y0, double x1, double y1)
 979 {
 980     double bend = fabs(ks[0]) + fabs(.5 * ks[1]) + fabs(.125 * ks[2]) +
 981         fabs((1./48) * ks[3]);
 982
 983     if (bend < 1e-8) {
 984         printf("%g %g lineto\n", x1, y1);
 985     } else {
 986         double seg_ch = hypot(x1 - x0, y1 - y0);
 987         double seg_th = atan2(y1 - y0, x1 - x0);
 988         double xy[2];
 989         double ch, th;
 990         double scale, rot;
 991         double th_even, th_odd;
 992         double ul, vl;
 993         double ur, vr;
 994
 995         integrate_spiro(ks, xy);
 996         ch = hypot(xy[0], xy[1]);
 997         th = atan2(xy[1], xy[0]);
 998         scale = seg_ch / ch;
 999         rot = seg_th - th;
1000         if (bend < 1.) {
1001             th_even = (1./384) * ks[3] + (1./8) * ks[1] + rot;
1002             th_odd = (1./48) * ks[2] + .5 * ks[0];
1003             ul = (scale * (1./3)) * cos(th_even - th_odd);
1004             vl = (scale * (1./3)) * sin(th_even - th_odd);
1005             ur = (scale * (1./3)) * cos(th_even + th_odd);
1006             vr = (scale * (1./3)) * sin(th_even + th_odd);
1007             printf("%g %g %g %g %g %g curveto\n",
1008                    x0 + ul, y0 + vl, x1 - ur, y1 - vr, x1, y1);
1009
1010         } else {
1011             /* subdivide */
1012             double ksub[4];
1013             double thsub;
1014             double xysub[2];
1015             double xmid, ymid;
1016             double cth, sth;
1017
1018             ksub[0] = .5 * ks[0] - .125 * ks[1] + (1./64) * ks[2] - (1./768) * ks[3];
1019             ksub[1] = .25 * ks[1] - (1./16) * ks[2] + (1./128) * ks[3];
1020             ksub[2] = .125 * ks[2] - (1./32) * ks[3];
1021             ksub[3] = (1./16) * ks[3];
1022             thsub = rot - .25 * ks[0] + (1./32) * ks[1] - (1./384) * ks[2] + (1./6144) * ks[3];
1023             cth = .5 * scale * cos(thsub);
1024             sth = .5 * scale * sin(thsub);
1025             integrate_spiro(ksub, xysub);
1026             xmid = x0 + cth * xysub[0] - sth * xysub[1];
1027             ymid = y0 + cth * xysub[1] + sth * xysub[0];
1028             print_seg(ksub, x0, y0, xmid, ymid);
1029             ksub[0] += .25 * ks[1] + (1./384) * ks[3];
1030             ksub[1] += .125 * ks[2];
1031             ksub[2] += (1./16) * ks[3];
1032             print_seg(ksub, xmid, ymid, x1, y1);
1033         }
1034     }
1035 }
1036
1037 void
1038 print_segs(const spiro_seg *segs, int nsegs)
1039 {
1040     int i;
1041
1042     for (i = 0; i < nsegs; i++) {
1043         double x0 = segs[i].x;
1044         double y0 = segs[i].y;
1045         double x1 = segs[i + 1].x;
1046         double y1 = segs[i + 1].y;
1047
1048         if (i == 0)
1049             printf("%g %g moveto\n", x0, y0);
1050         printf("%% ks = [ %g %g %g %g ]\n",
1051                segs[i].ks[0], segs[i].ks[1], segs[i].ks[2], segs[i].ks[3]);
1052         print_seg(segs[i].ks, x0, y0, x1, y1);
1053     }
1054     printf("stroke\n");
1055 }
1056
1057 int
1058 test_curve(void)
1059 {
1060     spiro_cp path[] = {
1061         {334, 117, 'v'},
1062         {305, 176, 'v'},
1063         {212, 142, 'c'},
1064         {159, 171, 'c'},
1065         {224, 237, 'c'},
1066         {347, 335, 'c'},
1067         {202, 467, 'c'},
1068         {81, 429, 'v'},
1069         {114, 368, 'v'},
1070         {201, 402, 'c'},
1071         {276, 369, 'c'},
1072         {218, 308, 'c'},
1073         {91, 211, 'c'},
1074         {124, 111, 'c'},
1075         {229, 82, 'c'}
1076     };
1077     spiro_seg *segs;
1078     int i;
1079
1080     n = 1;
1081     for (i = 0; i < 1000; i++) {
1082         segs = setup_path(path, 15);
1083         solve_spiro(segs, 15);
1084     }
1085     printf("100 800 translate 1 -1 scale 1 setlinewidth\n");
1086     print_segs(segs, 15);
1087     printf("showpage\n");
1088     return 0;
1089 }
1090
1091 int main(int argc, char **argv)
1092 {
1093     return test_curve();
1094 }
1095 #endif