1 files changed, 225 insertions, 0 deletions
diff --git a/src/3rdparty/2geom/tests/sbasis-text-test.cpp b/src/3rdparty/2geom/tests/sbasis-text-test.cpp
new file mode 100644
index 0000000..ef96407
--- /dev/null
+++ b/src/3rdparty/2geom/tests/sbasis-text-test.cpp
@@ -0,0 +1,225 @@
+#include <iostream>
+#include <math.h>
+#include <cassert>
+#include <2geom/sbasis.h>
+#include <2geom/sbasis-poly.h>
+#include <iterator>
+#include <2geom/point.h>
+#include <2geom/sbasis-to-bezier.h>
+#include <2geom/solver.h>
+
+using namespace Geom;
+
+Poly roots_to_poly(double *a, unsigned n) {
+    Poly r;
+    r.push_back(1);
+    
+    for(unsigned i = 0; i < n; i++) {
+        Poly p;
+        p.push_back(-a[i]);
+        p.push_back(1);
+        r = r*p;
+    }
+    return r;
+}
+
+unsigned small_rand() {
+    return (rand() & 0xff) + 1;
+}
+
+double uniform() {
+    return double(rand()) / RAND_MAX;
+}
+
+int main() {
+    SBasis P0(Linear(0.5, -1)), P1(Linear(3, 1));
+    Linear one(1,1);
+
+    std::cout << "round tripping of poly conversions\n";
+    std::cout << P0 
+              << "=>" << sbasis_to_poly(P0) 
+              << "=>" << poly_to_sbasis(sbasis_to_poly(P0)) 
+              << std::endl;
+    
+    std::cout << "derivatives and integrals\n";
+    
+    Poly test;
+    for(int i = 0; i < 4; i++)
+        test.push_back(1);
+    
+    SBasis test_sb = poly_to_sbasis(test);
+    std::cout << test << "(" << test.size() << ")"
+              << "   ==   "
+              << test_sb << "(" << test_sb.size() << ")"
+              << std::endl;
+    std::cout << "derivative\n";
+    std::cout << derivative(test)
+              << "   ==   "
+              << sbasis_to_poly(derivative(test_sb))
+              << std::endl;
+    
+    std::cout << "integral\n";
+    std::cout << integral(test)
+              << "   ==   "
+              << sbasis_to_poly(integral(test_sb))
+              << std::endl;
+    
+    std::cout << "evaluation\n";
+    std::cout << integral(test)(0.3) - integral(test)(0.)
+              << "   ==   "
+              << integral(test_sb)(0.3) - integral(test_sb)(0.)
+              << std::endl;
+    
+    std::cout << "multiplication\n";
+    std::cout << (test*test)
+              << "\n   ==   \n"
+              << sbasis_to_poly(multiply(test_sb,test_sb))
+              << std::endl;
+    std::cout << poly_to_sbasis(test*test)
+              << "\n   ==   \n"
+              << multiply(test_sb,test_sb)
+              << std::endl;
+    
+    std::cout << "sqrt\n";
+    std::cout << test
+              << "\n   ==   \n"
+              << sbasis_to_poly(sqrt(multiply(test_sb,test_sb),10))
+              << std::endl;
+    SBasis radicand = sqrt(test_sb,10);
+    std::cout << sbasis_to_poly(truncate(multiply(radicand, radicand),5))
+              << "\n   ==   \n"
+              << test
+              << std::endl;
+    
+    std::cout << "division\n";
+    std::cout << test
+              << "\n   ==   \n"
+              << sbasis_to_poly(divide(multiply(test_sb,test_sb),test_sb, 20))
+              << std::endl;
+    std::cout << divide(test_sb, radicand,5)
+              << "\n   ==   \n"
+              << truncate(radicand,6)
+              << std::endl;
+    
+    std::cout << "composition\n";
+    std::cout << (compose(test,sbasis_to_poly(Linear(0.5,1))))
+              << "\n   ==   \n"
+              << sbasis_to_poly(compose(test_sb,Linear(0.5,1)))
+              << std::endl;
+    std::cout << poly_to_sbasis(compose(test,test))
+              << "\n   ==   \n"
+              << compose(test_sb,test_sb)
+              << std::endl
+              << std::endl;
+    std::cout << (compose(sbasis_to_poly(Linear(1,2)),sbasis_to_poly(Linear(-1,0))))
+              << std::endl;
+    std::cout << (compose(SBasis(Linear(1,2)),SBasis(Linear(-1,0))))
+              << std::endl;
+
+    std::cout << "inverse of x - 1\n";
+    std::cout << sbasis_to_poly(inverse(Linear(-1,0),2))
+              << "   ==   y + 1\n";
+    std::cout << "f^-1(f(x)) = " 
+              << sbasis_to_poly(compose(inverse(Linear(-1,0),2), 
+                                        Linear(-1,0))) 
+              << std::endl
+              << std::endl;
+    
+    std::cout << "inverse of 3x - 2\n";
+    std::cout << sbasis_to_poly(inverse(Linear(-2,1),2))
+              << "   ==   (y + 2)/3\n";
+    std::cout << "f^-1(f(x)) = " 
+              << sbasis_to_poly(compose(inverse(Linear(-2,1),2), 
+                                        Linear(-2,1))) 
+              << std::endl
+              << std::endl;
+    
+    
+    std::cout << "inverse of sqrt(" << sbasis_to_poly(Linear(1,4)) << ") - 1\n";
+    SBasis A = sqrt(Linear(1,4), 5) - one;
+    Poly P;
+    P.push_back(0);
+    P.push_back(2./3);
+    P.push_back(1./3);
+    
+    std::cout << "2 term approximation\n";
+    std::cout << sbasis_to_poly(inverse(A,2))
+              << "\n   ==   \n"
+              << P
+              << std::endl;
+    std::cout << "general approximation\n";
+    std::cout << sbasis_to_poly(inverse(A,5))
+              << "\n   ==   \n"
+              << P
+              << std::endl;
+    
+    {
+    std::cout << "inverse of (x^2+2x)/3\n";
+    SBasis A = poly_to_sbasis(P);
+    SBasis I = inverse(A,10);
+    std::cout << sbasis_to_poly(truncate(compose(A, I), 10))
+              << "   ==   x\n"
+              << std::endl;
+    std::cout << sbasis_to_poly(truncate(compose(I, A), 10))
+              << "   ==   x\n"
+              << std::endl;
+    std::cout << sbasis_to_poly(truncate(I - (sqrt(Linear(1,4), 10) - one), 10))
+              << std::endl;
+    }
+#ifdef HAVE_GSL
+    for(int i = 0 ; i < 10; i++) {
+        Poly P;
+        P.push_back(0);
+        P.push_back(1);
+        for(int j = 0 ; j < 2; j++) {
+            P.push_back((uniform()-0.5)/10);
+        }
+        std::vector<std::complex<double> > prod_root = solve(derivative(P));
+        copy(prod_root.begin(), prod_root.end(), 
+             std::ostream_iterator<std::complex<double> >(std::cout, ",\t"));
+        std::cout << std::endl;
+        std::cout << "inverse of (" << P << " )\n";
+        for(int k = 1; k < 30; k++) {
+            SBasis A = poly_to_sbasis(P);
+            SBasis I = inverse(A,k);
+            SBasis err = compose(A, I) - Linear(0,1); // ideally is 0
+            std::cout << truncate(err, k).tailError(0)
+                      << std::endl;
+            /*std::cout << sbasis_to_poly(err)
+                      << "   ==   x\n"
+                      << std::endl;*/
+        }
+    }
+#endif
+    /*double roots[] = {0.1,0.2,0.6};
+    Poly prod = roots_to_poly(roots, sizeof(roots)/sizeof(double));
+    std::cout << "real solve\n";
+    std::cout << prod 
+              << " solves as ";
+    std::vector<std::complex<double> > prod_root = solve(prod);
+    copy(prod_root.begin(), prod_root.end(), 
+         std::ostream_iterator<std::complex<double> >(std::cout, ",\t"));
+    std::cout << std::endl;
+    
+    SBasis prod_sb = poly_to_sbasis(prod);
+    std::vector<double> bez = sbasis_to_bezier(prod_sb, prod_sb.size());
+    
+    copy(bez.begin(), bez.end(), 
+         std::ostream_iterator<double>(std::cout, ",\t"));
+         std::cout << std::endl;*/
+    
+    /*std::cout << "crossing count = " 
+              << crossing_count(bez, bez.size())
+              << std::endl;*/
+}
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :