libs/evoral/test/CurveTest.cpp

   1 #include "CurveTest.hpp"
   2 #include "evoral/ControlList.hpp"
   3 #include "evoral/Curve.hpp"
   4 #include <stdlib.h>
   5
   6 CPPUNIT_TEST_SUITE_REGISTRATION (CurveTest);
   7
   8 using namespace Evoral;
   9
  10 // linear y = Y0 + YS * x ;  with x = i * (X1 - X0) + X0; and i = [0..1023]
  11 #define VEC1024LINCMP(X0, X1, Y0, YS)                                        \
  12     cl->curve ().get_vector ((X0), (X1), vec, 1024);                         \
  13     for (int i = 0; i < 1024; ++i) {                                         \
  14         char msg[64];                                                        \
  15         snprintf (msg, 64, "at i=%d (x0=%.1f, x1=%.1f, y0=%.1f, ys=%.3f)",   \
  16             i, X0, X1, Y0, YS);                                              \
  17         CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE (                               \
  18             msg,                                                             \
  19             (Y0) + i * (YS), vec[i],                                         \
  20             1e-24                                                            \
  21             );                                                               \
  22     }
  23
  24 void
  25 CurveTest::twoPointLinear ()
  26 {
  27         float vec[1024];
  28
  29         boost::shared_ptr<Evoral::ControlList> cl = TestCtrlList();
  30
  31         cl->create_curve ();
  32         cl->set_interpolation (ControlList::Linear);
  33
  34         // add two points to curve
  35         cl->fast_simple_add (   0.0 , 2048.0);
  36         cl->fast_simple_add (8192.0 , 4096.0);
  37
  38         cl->curve ().get_vector (1024.0, 2047.0, vec, 1024);
  39
  40         VEC1024LINCMP (1024.0, 2047.0, 2304.f,  .25f);
  41         VEC1024LINCMP (2048.0, 2559.5, 2560.f,  .125f);
  42         VEC1024LINCMP (   0.0, 4092.0, 2048.f, 1.f);
  43
  44         // greetings to tartina
  45         cl->curve ().get_vector (2048.0, 2048.0, vec, 1);
  46         CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=1 @ 2048..2048", 2560.f, vec[0]);
  47
  48         /* XXX WHAT DO WE EXPECT WITH veclen=1 AND  x1 > x0 ? */
  49 #if 0
  50         /* .. interpolated value at (x1+x0)/2 */
  51         cl->curve ().get_vector (2048.0, 2049.0, vec, 1);
  52         CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=1 @ 2048-2049", 2560.125f, vec[0]);
  53
  54         cl->curve ().get_vector (2048.0, 2056.0, vec, 1);
  55         CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=1 @ 2048-2049", 2561.f, vec[0]);
  56 #else
  57         /* .. value at x0 */
  58         cl->curve ().get_vector (2048.0, 2049.0, vec, 1);
  59         CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=1 , 2048..2049", 2560.f, vec[0]);
  60
  61         cl->curve ().get_vector (2048.0, 2056.0, vec, 1);
  62         CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=1 , 2048..2049", 2560.f, vec[0]);
  63 #endif
  64
  65         cl->curve ().get_vector (2048.0, 2048.0, vec, 2);
  66         CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=2 , 2048..2048 @ 0", 2560.f, vec[0]);
  67         CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=2 , 2048..2048 @ 1", 2560.f, vec[1]);
  68
  69         cl->curve ().get_vector (2048.0, 2056.0, vec, 2);
  70         CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=2 , 2048..2056 @ 0", 2560.f, vec[0]);
  71         CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=2 , 2048..2056 @ 0", 2562.f, vec[1]);
  72
  73         cl->curve ().get_vector (2048.0, 2056.0, vec, 3);
  74         CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=3 , 2048..2056 @ 0", 2560.f, vec[0]);
  75         CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=3 , 2048..2056 @ 1", 2561.f, vec[1]);
  76         CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=3 , 2048..2056 @ 2", 2562.f, vec[2]);
  77
  78         /* check out-of range..
  79          * we expect the first and last value - no interpolation
  80          */
  81         cl->curve ().get_vector (-1, -1, vec, 1);
  82         CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=1 @ -1", 2048.f, vec[0]);
  83
  84         cl->curve ().get_vector (9999.0, 9999.0, vec, 1);
  85         CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=1 @ 9999", 4096.f, vec[0]);
  86
  87         cl->curve ().get_vector (-999.0, 0, vec, 13);
  88         for (int i = 0; i < 13; ++i) {
  89                 CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=13 @ -999..0", 2048.f, vec[i]);
  90         }
  91
  92         cl->curve ().get_vector (9998.0, 9999.0, vec, 8);
  93         for (int i = 0; i < 8; ++i) {
  94                 CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=8 @ 9998..9999", 4096.f, vec[i]);
  95         }
  96 }
  97
  98 void
  99 CurveTest::threePointLinear ()
 100 {
 101         float vec[4];
 102
 103         boost::shared_ptr<Evoral::ControlList> cl = TestCtrlList();
 104
 105         cl->create_curve ();
 106         cl->set_interpolation (ControlList::Linear);
 107
 108         // add 3 points to curve
 109         cl->fast_simple_add (   0.0 , 2.0);
 110         cl->fast_simple_add ( 100.0 , 4.0);
 111         cl->fast_simple_add ( 200.0 , 0.0);
 112
 113         cl->curve ().get_vector (50.0, 60.0, vec, 1);
 114         CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=1 @ 50", 3.f, vec[0]);
 115
 116         cl->curve ().get_vector (100.0, 100.0, vec, 1);
 117         CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=1 @ 100", 4.f, vec[0]);
 118
 119         cl->curve ().get_vector (150.0, 150.0, vec, 1);
 120         CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=1 @ 150", 2.f, vec[0]);
 121
 122         cl->curve ().get_vector (130.0, 150.0, vec, 3);
 123         CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=3 130..150 @ 0", 2.8f, vec[0]);
 124         CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=3 130..150 @ 2", 2.4f, vec[1]);
 125         CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=3 130..150 @ 3", 2.0f, vec[2]);
 126
 127         cl->curve ().get_vector (80.0, 160.0, vec, 3);
 128         CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=3 80..160 @ 0", 3.6f, vec[0]);
 129         CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=3 80..160 @ 2", 3.2f, vec[1]);
 130         CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=3 80..160 @ 3", 1.6f, vec[2]);
 131 }
 132
 133 void
 134 CurveTest::threePointDiscete ()
 135 {
 136         boost::shared_ptr<Evoral::ControlList> cl = TestCtrlList();
 137         cl->set_interpolation (ControlList::Discrete);
 138
 139         // add 3 points to curve
 140         cl->fast_simple_add (   0.0 , 2.0);
 141         cl->fast_simple_add ( 100.0 , 4.0);
 142         cl->fast_simple_add ( 200.0 , 0.0);
 143
 144         CPPUNIT_ASSERT_EQUAL(2.0, cl->unlocked_eval(80.));
 145         CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(120.));
 146         CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(160.));
 147
 148         cl->set_interpolation (ControlList::Linear);
 149
 150         CPPUNIT_ASSERT_EQUAL(3.6, cl->unlocked_eval(80.));
 151         CPPUNIT_ASSERT_EQUAL(3.2, cl->unlocked_eval(120.));
 152         CPPUNIT_ASSERT_EQUAL(1.6, cl->unlocked_eval(160.));
 153 }
 154
 155 void
 156 CurveTest::ctrlListEval ()
 157 {
 158         boost::shared_ptr<Evoral::ControlList> cl = TestCtrlList();
 159
 160         cl->fast_simple_add (   0.0 , 2.0);
 161
 162         cl->set_interpolation (ControlList::Discrete);
 163         CPPUNIT_ASSERT_EQUAL(2.0, cl->unlocked_eval(80.));
 164         CPPUNIT_ASSERT_EQUAL(2.0, cl->unlocked_eval(120.));
 165         CPPUNIT_ASSERT_EQUAL(2.0, cl->unlocked_eval(160.));
 166
 167         cl->set_interpolation (ControlList::Linear);
 168         CPPUNIT_ASSERT_EQUAL(2.0, cl->unlocked_eval(80.));
 169         CPPUNIT_ASSERT_EQUAL(2.0, cl->unlocked_eval(120.));
 170         CPPUNIT_ASSERT_EQUAL(2.0, cl->unlocked_eval(160.));
 171
 172         cl->fast_simple_add ( 100.0 , 4.0);
 173
 174         cl->set_interpolation (ControlList::Discrete);
 175         CPPUNIT_ASSERT_EQUAL(2.0, cl->unlocked_eval(80.));
 176         CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(120.));
 177         CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(160.));
 178
 179         cl->set_interpolation (ControlList::Linear);
 180         CPPUNIT_ASSERT_EQUAL(3.6, cl->unlocked_eval(80.));
 181         CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(120.));
 182         CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(160.));
 183
 184         cl->fast_simple_add ( 200.0 , 0.0);
 185
 186         cl->set_interpolation (ControlList::Discrete);
 187         CPPUNIT_ASSERT_EQUAL(2.0, cl->unlocked_eval(80.));
 188         CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(120.));
 189         CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(160.));
 190
 191         cl->set_interpolation (ControlList::Linear);
 192         CPPUNIT_ASSERT_EQUAL(3.6, cl->unlocked_eval(80.));
 193         CPPUNIT_ASSERT_EQUAL(3.2, cl->unlocked_eval(120.));
 194         CPPUNIT_ASSERT_EQUAL(1.6, cl->unlocked_eval(160.));
 195
 196         cl->fast_simple_add ( 300.0 , 8.0);
 197
 198         cl->set_interpolation (ControlList::Discrete);
 199         CPPUNIT_ASSERT_EQUAL(2.0, cl->unlocked_eval(80.));
 200         CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(120.));
 201         CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(160.));
 202         CPPUNIT_ASSERT_EQUAL(0.0, cl->unlocked_eval(250.));
 203         CPPUNIT_ASSERT_EQUAL(8.0, cl->unlocked_eval(999.));
 204
 205         cl->set_interpolation (ControlList::Linear);
 206         CPPUNIT_ASSERT_EQUAL(3.6, cl->unlocked_eval(80.));
 207         CPPUNIT_ASSERT_EQUAL(3.2, cl->unlocked_eval(120.));
 208         CPPUNIT_ASSERT_EQUAL(1.6, cl->unlocked_eval(160.));
 209         CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(250.));
 210         CPPUNIT_ASSERT_EQUAL(8.0, cl->unlocked_eval(999.));
 211
 212         cl->fast_simple_add ( 400.0 , 9.0);
 213
 214         cl->set_interpolation (ControlList::Discrete);
 215         CPPUNIT_ASSERT_EQUAL(2.0, cl->unlocked_eval(80.));
 216         CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(120.));
 217         CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(160.));
 218         CPPUNIT_ASSERT_EQUAL(0.0, cl->unlocked_eval(250.));
 219         CPPUNIT_ASSERT_EQUAL(8.0, cl->unlocked_eval(350.));
 220         CPPUNIT_ASSERT_EQUAL(9.0, cl->unlocked_eval(999.));
 221
 222         cl->set_interpolation (ControlList::Linear);
 223         CPPUNIT_ASSERT_EQUAL(3.6, cl->unlocked_eval(80.));
 224         CPPUNIT_ASSERT_EQUAL(3.2, cl->unlocked_eval(120.));
 225         CPPUNIT_ASSERT_EQUAL(1.6, cl->unlocked_eval(160.));
 226         CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(250.));
 227         CPPUNIT_ASSERT_EQUAL(8.5, cl->unlocked_eval(350.));
 228         CPPUNIT_ASSERT_EQUAL(9.0, cl->unlocked_eval(999.));
 229 }
 230
 231 void
 232 CurveTest::constrainedCubic ()
 233 {
 234
 235         struct point {
 236                 int x, y;
 237         };
 238
 239         static const struct point data[] = {
 240                 /* values from worked example in www.korf.co.uk/spline.pdf */
 241                 {   0,  30 },
 242                 {  10, 130 },
 243                 {  30, 150 },
 244                 {  50, 150 },
 245                 {  70, 170 },
 246                 {  90, 220 },
 247                 { 100, 320 },
 248         };
 249
 250         int32_t type = 0;
 251         Evoral::Parameter p(type);
 252         Evoral::ParameterDescriptor pd;
 253         Evoral::ControlList l(p,pd);
 254
 255         size_t i;
 256         l.set_interpolation(Evoral::ControlList::Curved);
 257
 258         for (i=0; i<sizeof(data)/sizeof(data[0]); i++) {
 259                 l.add (data[i].x, data[i].y);
 260         }
 261
 262         Evoral::Curve curve(l);
 263
 264         float f[121];
 265         curve.get_vector(-10, 110, f, 121);
 266
 267         const float *g = &f[10]; /* so g starts at x==0 */
 268
 269         /* given points - should be exactly equal */
 270         CPPUNIT_ASSERT_EQUAL( 30.0f, g[-10]);
 271         CPPUNIT_ASSERT_EQUAL( 30.0f, g[  0]);
 272         CPPUNIT_ASSERT_EQUAL(130.0f, g[ 10]);
 273         CPPUNIT_ASSERT_EQUAL(150.0f, g[ 30]);
 274         CPPUNIT_ASSERT_EQUAL(150.0f, g[ 40]);
 275         CPPUNIT_ASSERT_EQUAL(150.0f, g[ 50]);
 276         CPPUNIT_ASSERT_EQUAL(320.0f, g[100]);
 277         CPPUNIT_ASSERT_EQUAL(320.0f, g[110]);
 278
 279         /*
 280            First segment, i=1, for 0 <= x <= 10
 281            f'1(x1) = 2/((x2 – x1)/(y2 – y1) + (x1 – x0)/(y1 – y0))
 282                    = 2/((30 – 10)/(150 – 130) + (10 – 0)/(130 – 30))
 283                    = 1.8181
 284            f'1(x0) = 3/2*(y1 – y0)/(x1 – x0) - f'1(x1)/2
 285                    = 3/2*(130 – 30)/(10 – 0) – 1.818/2
 286                    = 14.0909
 287            f"1(x0) = -2*(f'1(x1) + 2* f'1(x0))/(x1 – x0) + 6*(y1 – y0)/ (x1 – x0)^2
 288                    = -2*(1.8181 + 2*14.0909)/(10 – 0) + 6*(130 – 30)/(10 – 0)^2
 289                    = 0
 290            f"1(x1) = 2*(2*f'1(x1) + f'1(x0))/(x1 – x0) - 6*(y1 – y0)/ (x1 – x0)^2
 291                    = 2*(2*1.818 + 14.0909)/(10 – 0) – 6*(130 – 30)/(10 – 0)^2
 292                    = -2.4545
 293            d1 = 1/6 * (f"1(x1) - f"1(x0))/(x1 – x0)
 294               = 1/6 * (-2.4545 – 0)/(10 – 0)
 295               = -0.0409
 296            c1 = 1/2 * (x1*f"1(x0) – x0*f"1(x1))/(x1 – x0)
 297               = 1/2 * (10*0 – 0*1.8181)/(10 – 0)
 298               = 0
 299            b1 = ((y1 – y0) – c1*(x21 – x20) – d1*( x31 – x30))/(x1 – x0)
 300               = ((130 – 30) – 0*(102 – 02) + 0.0409*(103 – 03))/(10 – 0)
 301               = 14.09
 302            a1 = y0 – b1*x0 – c1*x20 – d1*x30
 303               = 30
 304            y1 = 30 + 14.09x - 0.0409x3 for 0 <= x <= 10
 305          */
 306         /*
 307            Second segment, i=2, for 10 <= x <= 30
 308            f'2(x2) = 2/((x3 – x2)/(y3 – y2) + (x2 – x1)/(y2 – y1))
 309                    = 2/((50 – 30)/(150 – 150) + (30 – 10)/(150 – 130))
 310                    = 0
 311            f'2(x1) = 2/((x2 – x1)/(y2 – y1) + (x1 – x0)/(y1 – y0))
 312                    = 1.8181
 313
 314            f"2(x1) = -2*(f'2(x2) + 2* f'2(x1))/(x2 – x1) + 6*(y2 – y1)/ (x2 – x1)^2
 315                    = -2*(0 + 2*1.8181)/(30 – 10) + 6*(150 – 130)/(30 – 10)2
 316                    = -0.063636
 317            f"2(x2) = 2*(2*f'2(x2) + f'2(x1))/(x2 – x1) - 6*(y2 – y1)/ (x2 – x1)^2
 318                    = 2*(2*0 + 1.8181)/(30 – 10) – 6*(150 – 130)/(30 – 10)^2
 319                    = -0.11818
 320
 321            d2 = 1/6 * (f"2(x2) - f"2(x1))/(x2 – x1)
 322               = 1/6 * (-0.11818 + 0.063636)/(30 – 10)
 323               = -0.0004545
 324            c2 = 1/2 * (x2*f"2(x1) – x1*f"2(x2))/(x2 – x1)
 325               = 1/2 * (-30*0.063636 + 10*0.11818)/(30 – 10)
 326               = -0.01818
 327            b2 = ((y2 – y1) – c2*(x2^2 – x1^2) – d2*( x2^3 – x1^3))/(x2 – x1)
 328               = ((150 – 130) + 0.01818*(302 – 102) + 0.0004545*(303 – 103))/(30 – 10)
 329               = 2.31818
 330            a2 = y1 – b2*x1 – c2*x1^2 – d2*x1^3
 331               = 130 – 2.31818*10 + 0.01818*102 + 0.0004545*103
 332               = 109.09
 333            y2 = 109.09 + 2.31818x - 0.01818x^2 - 0.0004545x^3 for 10 <= x <= 30
 334          */
 335
 336
 337         int x;
 338         long double a1, b1, c1, d1, a2, b2, c2, d2, fdx0, fddx0, fdx1, fdx2, fddx1, fddx2;
 339         double x0 = data[0].x;
 340         double y0 = data[0].y;
 341         double x1 = data[1].x;
 342         double y1 = data[1].y;
 343         double x2 = data[2].x;
 344         double y2 = data[2].y;
 345         double x3 = data[3].x;
 346         double y3 = data[3].y;
 347
 348         double dx0 = x1 - x0;
 349         double dy0 = y1 - y0;
 350         double dx1 = x2 - x1;
 351         double dy1 = y2 - y1;
 352         double dx2 = x3 - x2;
 353         double dy2 = y3 - y2;
 354
 355         // First (leftmost) segment
 356         fdx1 = 2.0 / ( dx1 / dy1 + dx0 / dy0 );
 357         fdx0 = 3.0 / 2.0 * dy0 / dx0 - fdx1 / 2.0;
 358
 359         fddx0 = -2.0 * (fdx1 + 2.0 * fdx0) / dx0 + 6.0 * dy0 / (dx0*dx0);
 360         fddx1 =  2.0 * (2.0 * fdx1 + fdx0) / dx0 - 6.0 * dy0 / (dx0*dx0);
 361         d1 = 1.0 / 6.0 * (fddx1 - fddx0) / dx0;
 362         c1 = 1.0 / 2.0 * (x1 * fddx0 - x0 * fddx1) / dx0;
 363         b1 = (dy0 - c1 * (x1* x1 - x0*x0) - d1 * (x1*x1*x1 - x0*x0*x0)) / dx0;
 364         a1 = y0 - b1*x0 - c1*x0*x0 - d1*x0*x0*x0;
 365
 366         // printf("dx0=%f, dy0=%f, dx1=%f, dy1=%f\n", dx0, dy0, dx1, dy1);
 367         // printf("fdx0=%Lf, fdx1=%Lf, fddx0=%Lf, fddx1=%Lf\n", fdx0, fdx1, fddx0, fddx1);
 368         // printf("a1=%Lf, b1=%Lf, c1=%Lf, d1=%Lf\n", a1, b1, c1, d1);
 369
 370         // values from worked example: deltas rather arbitrary, I'm afraid
 371         CPPUNIT_ASSERT_DOUBLES_EQUAL(30.0, a1, 0);
 372         CPPUNIT_ASSERT_DOUBLES_EQUAL(14.09, b1, 0.001);
 373         CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, c1, 0);
 374         CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.0409, d1, 0.00001);
 375
 376         for (x = 0; x <= 10; x++) {
 377                 double v = a1 + b1*x + c1*x*x  + d1*x*x*x;
 378                 char msg[64];
 379                 snprintf(msg, 64, "interpolating %d: v=%f, x=%f...\n", x, v, g[x]);
 380                 CPPUNIT_ASSERT_DOUBLES_EQUAL(v, g[x], 0.000004);
 381         }
 382
 383         // Second segment
 384         fdx2 = 2.0 / ( dx2 / dy2 + dx1 / dy1 );
 385
 386         fddx1 = -2.0 * (fdx2 + 2.0 * fdx1) / dx1 + 6.0 * dy1 / (dx1*dx1);
 387         fddx2 =  2.0 * (2.0 * fdx2 + fdx1) / dx1 - 6.0 * dy1 / (dx1*dx1);
 388         d2 = 1.0 / 6.0 * (fddx2 - fddx1) / dx1;
 389         c2 = 1.0 / 2.0 * (x2 * fddx1 - x1 * fddx2) / dx1;
 390         b2 = (dy1 - c2 * (x2*x2 - x1*x1) - d2 * (x2*x2*x2 - x1*x1*x1)) / dx1;
 391         a2 = y1 - b2*x1 - c2*x1*x1 - d2*x1*x1*x1;
 392
 393         // printf("dx0=%f, dy0=%f, dx1=%f, dy1=%f dx2=%f, dy2=%f\n", dx0, dy0, dx1, dy1, dx2, dy2);
 394         // printf("fdx1=%Lf, fdx2=%Lf, fddx1=%Lf, fddx2=%Lf\n", fdx1, fdx2, fddx1, fddx2);
 395         // printf("a2=%Lf, b2=%Lf, c2=%Lf, d2=%Lf\n", a2, b2, c2, d2);
 396
 397         // values from worked example: deltas rather arbitrary, I'm afraid
 398         CPPUNIT_ASSERT_DOUBLES_EQUAL(109.09, a2, 0.001);
 399         CPPUNIT_ASSERT_DOUBLES_EQUAL(2.31818, b2, 0.00001);
 400         CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.01818, c2, 0.00001);
 401         CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.0004545, d2, 0.0000001);
 402
 403         for (x = 10; x <= 30; x++) {
 404                 double v = a2 + b2*x + c2*x*x  + d2*x*x*x;
 405                 char msg[64];
 406                 snprintf(msg, 64, "interpolating %d: v=%f, x=%f...\n", x, v, g[x]);
 407                 CPPUNIT_ASSERT_DOUBLES_EQUAL(v, g[x], 0.000008);
 408         }
 409 }