diff options
Diffstat (limited to '')
-rw-r--r-- | tests/sbasis-text-test.cpp | 225 |
1 files changed, 225 insertions, 0 deletions
diff --git a/tests/sbasis-text-test.cpp b/tests/sbasis-text-test.cpp new file mode 100644 index 0000000..ef96407 --- /dev/null +++ b/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 : |