x02_numeric_adaptor_example.cpp

The Numeric Adaptor is recently introduced to the Boost mailing list It is a proxy to high precision arithmetic libraries such as GMP or CLN. However, it might be that the same effect can be used using the Boost.Math bindings.

This example shows how the GGL can be combined with non-numeric value types such as from the Numeric Adaptor

// Generic Geometry Library
//
// Copyright Barend Gehrels 1995-2009, Geodan Holding B.V. Amsterdam, the Netherlands.
// Copyright Bruno Lalande 2008, 2009
// Use, modification and distribution is subject to the Boost Software License,
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)

// Numeric Adaptor Example

// The Numeric Adaptor is introduced to the Boost mailing list
// It is a proxy to high precision arithmetic libraries such as gmp
// However, it might be that the same effect can be used using the
// Boost.Math bindings.

// To build and run this example:
// 1) download gmp, it should be somewhere in the include path
// 2) download numeric_adaptor from the Boost.Sandbox

#include <iostream>
#include <iomanip>
#include <typeinfo>

#include <boost/numeric_adaptor/numeric_adaptor.hpp>
#include <boost/numeric_adaptor/gmp_value_type.hpp>

#include <ggl/ggl.hpp>
#include <ggl/geometries/cartesian2d.hpp>

template <typename Type, typename AssignType>
void calculate(AssignType const& x1,
               AssignType const& y1,
               AssignType const& x2,
               AssignType const& y2,
               AssignType const& x3,
               AssignType const& y3
               )
{
    // gmp can be used instead of "double" for any point type
    typedef ggl::point_xy<Type> point_type;

    point_type a, b, c;
    ggl::assign(a, boost::to<Type>(x1), boost::to<Type>(y1));
    ggl::assign(b, boost::to<Type>(x2), boost::to<Type>(y2));
    ggl::assign(c, boost::to<Type>(x3), boost::to<Type>(y3));

    ggl::linear_ring<point_type> r;
    r.push_back(a);
    r.push_back(b);
    r.push_back(c);
    r.push_back(a);

    // What also is possible is define point coordinates using IEEE double,
    // but doing calculations using the gmp type.
    // To do that, specify the strategy explicitly
    Type ab = ggl::distance(a, b);
    Type bc = ggl::distance(b, c);
    Type ca = ggl::distance(c, a);

    std::cout << std::endl << typeid(Type).name() << std::endl;

    std::cout << "a-b: " << ab << std::endl;
    std::cout << "b-c: " << bc << std::endl;
    std::cout << "c-a: " << ca << std::endl;

    std::cout << "area:  " << ggl::area(r,
            ggl::strategy::area::by_triangles<point_type, Type>())
        << std::endl;

    // Heron formula is "famous" for its imprecision. It should give
    // same result as area, but is sensible for rounding errors.
    Type s = ab + bc + ca;
    s /= 2.0;
    Type ar = boost::sqrt(s * (s - ab) * (s - bc) * (s - ca));
    std::cout << "heron: " << ar << std::endl;

    // Area's given:
    // float:       740.74078369140625
    // double:      740.74073407406990554591
    // long double: 740.74073407406991376156
    // GMP:        0.74074073407407e3 (right!)

    // SQL Server:  740.740724252642
    // Postgis:     740.74073407407 (might be rounded from (long) double)
    // MySQL:       740.74073407406991000000
}

int main(int argc, char *argv[])
{
    typedef boost::numeric_adaptor::gmp_value_type type;
    //typedef long double type;

    // gmp can be used instead of "double" for any point type
    typedef ggl::point_xy<type> point_type;

    // Points, polygons or other geometries are gmp now
    point_type p;

    // They can be used normally
    ggl::set<0>(p, 123456.78900001);
    std::cout << "x coordinate: " << ggl::get<0>(p) << std::endl;

    // But the value above cannot be expressed with that precision in IEEE 64 bits.
    // Points can therefore also be assigned by string with ggl::set
    // (boost::to is a converser included in the Numeric Adaptor sources)
    ggl::set<0>(p, boost::to<type>(std::string("123456.78900001")));

    // and streamed (as a string representation)
    std::cout << "x coordinate: " << ggl::get<0>(p) << std::endl;


    // The ggl::assign function also supports custom numeric types
    point_type p1, p2;
    ggl::assign(p1,
        boost::to<type>(std::string("123456.78900001")),
        boost::to<type>(std::string("234567.89100001")));
    ggl::assign(p2,
        boost::to<type>(std::string("987654.32100001")),
        boost::to<type>(std::string("876543.21900001")));

    type d = ggl::distance(p1, p2);
    std::cout << "Exact distance: " << d << std::endl;
    // It gives:          0.1076554 54858339556783e7
    // my calculator gives: 1076554.5485833955678294387789057
    // CLN gives          : 1076554.5485833955

    // All algorithms will automatically use the gmp-type
    // We show and compare that in the calculate function, with type and
    // assigning type as template parameters

    std::cout << std::fixed << std::setprecision(20);
    calculate<float>(0.0, 0.0, 0.0, 0.0012, 1234567.89012345, 0.0);
    calculate<double>(0.0, 0.0, 0.0, 0.0012, 1234567.89012345, 0.0);
    calculate<long double>(0.0, 0.0, 0.0, 0.0012, 1234567.89012345, 0.0);
    calculate<type, std::string>("0.0", "0.0", "0.0", "0.0012", "1234567.89012345", "0.0");

    return 0;
}