C++20 Concepts

This week (2021) I'm at the C++Now conference. The conference is normally in Aspen, Colorado (I was there more than 10 years ago), but this year it's online, for obvious reasons.

Yesterday there were two excellent interesting talks from Jeff Garland about concepts in C++20, how it works, what you can do with it, how you can use them and how you can write them.

Some compilers (notably clang) don't yet support concepts, but you can already play with toy projects online using https://godbolt.org, where your code is compiled and run on the fly.

So here is my toy project, which I wrote during the talks (there were two consecutive talks):

 

#include <array>

#include <type_traits>

#include <iostream>

template <typename C> concept ValidCoordinateType = std::is_arithmetic_v<C>;

template <int D> concept ValidDimension = D >= 2 and D <= 3;

template <int D> concept HasZ = D >= 3;

template <typename C, int D>

requires(ValidDimension<D> and ValidCoordinateType<C>)

struct mypoint

{

  mypoint() = default;

 

  mypoint(C x, C y) 

    : coors{x, y} {}

 

  // Only available for 3D

  mypoint(C x, C y, C z) requires HasZ<D>

    : coors{x, y, z} {}

 

  auto x() const { return coors[0]; }

  auto y() const { return coors[1]; }

  // Avoids compilation for Dim < 3

  auto z() const requires HasZ<D> { return coors[2]; }

 

private :

   std::array<C, D> coors;

};

struct mytype {};

int main()

{

  mypoint<double, 2> two(1, 2);

  mypoint<float, 3> three(3, 4, 5);

  std::cout << "Hi " << two.x() << " " << two.y() << " " << three.z() << "\n";

  // These declarations will not compile

  //mypoint<mytype, 2> p1;

  //mypoint<double, 4> p2;

 

  // This line will not compile

  //std::cout << two.z(); // Fails because .z() is not available

  return 0;

}

So this code is using C++20 concepts. You can see it clearly at and around the keywords concept and requires. It is really cool, especially compared with what we had to do with C++03 (Boost.Geometry was written in C++03 and only a few releases ago went to C++14) to achieve this.

Some things could earlier be done with static_assert too (such as checking the number of dimensions), but now the compiler neatly warns that that is better written using a concept.

Other things needed SFINAE earlier (such as the enabling  / disabling the functions for the Z coordinate).

If the last declarations in main are uncommented, you get a neat compiler error report (for instance for p2):

<source>: In function 'int main()':
<source>:40:20: error: template constraint failure for 'template<class C, int D>  requires (ValidDimension<D>) && (ValidCoordinateType<C>) struct mypoint'
   40 |   mypoint<double, 4> p2; // Fails because 4 coordinates are not allowed by the concept
      |                    ^
<source>:40:20: note: constraints not satisfied
<source>:6:26:   required for the satisfaction of 'ValidDimension<D>' [with D = 4]
<source>:6:56: note: the expression 'D <= 3 [with D = 4]' evaluated to 'false'
    6 | template <int D> concept ValidDimension = D >= 2 and D <= 3;
      |                                                      ~~^~~~

Wow, it's only a few lines! And so clear! In the past we got hundreds of hard to interpret error messages about templates...

And if you uncomment the line streaming z for a 2d coordinate, you get:

<source>: In function 'int main()':
<source>:43:22: error: no matching function for call to 'mypoint<double, 2>::z()'
   43 |   std::cout << two.z();
      |                      ^
<source>:23:8: note: candidate: 'auto mypoint<C, D>::z() const requires  HasZ<D> [with C = double; int D = 2]'
   23 |   auto z() const requires HasZ<D>
      |        ^
<source>:23:8: note: constraints not satisfied
<source>: In instantiation of 'auto mypoint<C, D>::z() const requires  HasZ<D> [with C = double; int D = 2]':
<source>:43:22:   required from here
<source>:7:26:   required for the satisfaction of 'HasZ<D>' [with D = 2]
<source>:7:35: note: the expression 'D >= 3 [with D = 2]' evaluated to 'false'
    7 | template <int D> concept HasZ = D >= 3;
      |                                 ~~^~~~

Union, Tangencies and Colors

I like the picture below very much.

What you see here are, more or less, pink, green and blue. The colors are chosen as such that they always form distinct colors, regardless of how you combine them. There are of course other combinations with these properties, but... this is at least one of them.

When the pink layer is slightly moved to the right, you see better how the image is constructed:

The pink (transparent) layer is the union of the two underlying (multi-)polygons, the green one and the blue one. It has one hole, one kind of indentation and for the rest it is square. The underlying multi-polygons are better seen if the pink one is removed:

So this is clear: there are two multi-polygons here, complete with holes and self-tangencies. These multi-polygons have overlap in some places, no overlap in other places. They nearly fill the area, there are only three "boxes" free. The blue and green colors nicely mix into blue-green. The colors and transparencies are selected for that. SVG definitions are:

Green: style="fill-opacity:0.5;fill:rgb(153,204,0);stroke:rgb(153,204,0);stroke-width:3"

Blue: style="fill-opacity:0.3;fill:rgb(51,51,153);stroke:rgb(51,51,153);stroke-width:3"

I like this combination.

And with the pink style on top of it...

Pink: style="fill-opacity:0.2;stroke-opacity:0.4;fill:rgb(255,0,0);stroke:rgb(255,0,255);stroke-width:8"

...all colors change, however, they are still clearly distinct. Yellow-orange like for green+pink; lila for blue+pink; and brown for the blue+green+pink. So we have seven unique distinct colors now, which happens to be the ideal number of colors for a map.

By the way, this image is created using Boost.Geometry. Boost.Geometry can easily create SVG images. And Boost.Geometry can overlay multi-polygons (already possible for a long time) having all kinds of self-tangencies (new!). A self-tangency in a multi-polygon is a place where two polygon-members touch each other (orange dots in the map below), or where a polygon touches itself (the red dot in the map).

Actually... this could also be a polygon, having five interior rings (holes) which might be mutually tangent and tangent to they exterior ring.

OK, one more image...

Tag Dispatching and Inheritance

In the previous blog we described Tag Dispatching by Type.

Though very convenient, when many classes (tagged by as many tags) share the same implementation, we get as many dispatch structs with the same implementation (or, better, forwarding to the same struct). This can be avoided by using similarities between classes, or between tags. In Boost.Geometry, for example, some geometries are single while others are multi. That multi-ness can be recognized by the metafunction is_multi which results in true for a multi-polygon and in false for a single polygon. We can include an IsMulti boolean template parameter in the dispatch structure, so that it is selected not only on tag dispatching, but also in is_multi.

This similarity can also be implemented using a tag hierarchy and inheritance.

Let's go back to the fruit example. All citrus fruit species have pulp vesicles. There are many citrus fruits species. And, of course, we don't want to repeat the dispatch structure has_vesicles so many times. So we try to define our tag hierarchy:

struct citrus_tag {};
struct pome_tag {};

struct apple_tag : pome_tag {};
struct pear_tag : pome_tag {};

struct banana_tag {};

struct orange_tag : citrus_tag {};
struct lemon_tag : citrus_tag {};
struct lime_tag : citrus_tag {};

And then we couple the tags, like we did in previous blog:

template <typename T> struct tag {};
template <> struct tag<apple> { typedef apple_tag type; };
template <> struct tag<pear> { typedef pear_tag type; };
template <> struct tag<orange> { typedef orange_tag type; };
template <> struct tag<lime> { typedef lime_tag type; };
template <> struct tag<banana> { typedef banana_tag type; };

The new thing, besides the inheritance in the tags above, is that we now define a tag_cast metafunction:

template
<
    typename Tag, typename BaseTag,
    typename BT2 = void, typename BT3 = void, typename BT4 = void,
    typename BT5 = void, typename BT6 = void, typename BT7 = void
>
struct tag_cast
{
    typedef typename boost::mpl::if_
        <
          typename boost::is_base_of<BaseTag, Tag>::type,
          BaseTag,
          // Try next one in line:
          typename tag_cast<Tag, BT2, BT3, BT4, BT5, BT6, BT7, void>::type
        >::type type;
};

template <typename Tag>
struct tag_cast<Tag, void, void, void, void, void, void, void>
{
    // If not found, take specified tag, so do not cast
    typedef Tag type;
};

We can now check if a tag is a citrus tag, by calling defining a new tag as tag_cast<tag, citrus_tag>::type. This new tag is either a citrus tag, or it is the tag it already was.

And we can dispatch by this tag as well. So it looks like:

namespace dispatch
{
    template <typename T> struct has_vesicles : boost::false_type {};
    template <> struct has_vesicles<citrus_tag> : boost::true_type {};
}

template <typename Fruit>
std::string has_vesicles(Fruit const& fruit)
{
    // (Potentially) go up in hierachy: take tag corresponding to Fruit,
    // downcast to citrus_tag if possible
    typedef typename tag_cast<typename tag<Fruit>::type, citrus_tag>::type tag;

    return std::string("has vesicles: ")
        + (dispatch::has_vesicles<tag>::value ? "true" : "false");
}

This tag_cast, as proposed above, can not only cast to one base-class, but also to a range of base-classes. The definition

 typedef typename tag_cast<typename tag<Fruit>::type, citrus_tag, pome_tag>::type tag;

will result in a citrus_tag if it is a citrus, a pome_tag if it is a pome, and otherwise it results in the input tag. Besides this, it also walks through a tag hierarchy, so if citrus_tag and pome_tag were both derived from, e.g., a rosid_tag, and rosid_tag was specified in this call, it would result in a rosid_tag.

The sample program here shows two base classes, not listed here.

Our main program now displays correctly if the specific fruits have vesicles or not

int main()
{
    using namespace fruit;

    apple a("my apple");
    pear p("my pear");
    orange o("my orange");

    std::cout << has_vesicles(a) << std::endl;
    std::cout << has_vesicles(p) << std::endl;
    std::cout << has_vesicles(o) << std::endl;

    return 0;
}

Back to Boost.Geometry. Until now it is not done, but we might depricate the is_multi meta-function and replace it by this:

struct multi_tag;
struct multi_point_tag : multi_tag {};
struct multi_linestring_tag : multi_tag {};
struct multi_polygon_tag : multi_tag {};

And instead of calling is_multi, we call tag_cast<tag, multi_tag>. We then know if it is a multi and can dispatch on it. The same for linear. This will probably be very convenient, it has to be worked out a little more.


Casting and (multiple) inheritance
The tag_cast as implemented above casts to specified tags. We can not just call the parent-tag. One of the reasons for this is multiple inheritance.
Tag inheritance can be multiple. This is convenient because, of course, there can be different tag classifications for the same set of tags. For example: if a geometry contains segments, or is linear, or is areal. So we might get:

struct multi_tag;
struct multi_point_tag : multi_tag {};
struct multi_linestring_tag : multi_tag, linear_tag {};
struct multi_polygon_tag : multi_tag, areal_tag {};

When calling tag_cast, the results depends on the order of the specified base tags. And so does the dispatching.

Tag Dispatching by Type

Tag Dispatching (TD) is a Generic Programming (GP) technique for C++. It is well explained in this blog. Google Hit Number one (GH#1) is this site, where the example of the STL function advance is given. Tag Dispatching is not described a lot; GH#7 is from ggl (= Boost.Geometry), as is GH#8..

Nearly all TD descriptions I know of use the example of STL-advance. But in Boost.Geometry TD is done in another way: not by instance, but by type. The difference is explained in this blog.

Suppose you make a generic program doing something with fruit. An apple is normally being eaten differently than a banana or an orange. So if we program a generic function eat, giving it any type of fruit, it should redirect this to an implementation specific to that kind of fruit, so for apple differently than for a banana.

We first implement the tags. Tags are completely empty structures. But (of course) not completely useless.

struct apple_tag {};
struct banana_tag {};
struct orange_tag {};

We then create some sample implementations. Most properties are not used, it is just to show there are (sometimes completely) different classes which can be handled genericly.

struct apple
{
    double radius;

    std::string name;
    apple(std::string const& n) : name(n) {}
};

struct banana
{
    double length;

    std::string name;
    banana(std::string const& n) : name(n) {}
};

All this is quite simple C++ code, though you have to understand templates (generics) and specialization.

OK, now we implement Tag Dispatching by Instance:

namespace dispatch
{
    void eat(apple const& a, apple_tag)
    {
        std::cout << "bite" << std::endl;
    }

    void eat(banana const& b, banana_tag)
    {
        std::cout << "peel" << std::endl;
    }
}

template <typename T>
void eat(T const& fruit)
{
    typename tag<T>::type the_tag;
    dispatch::eat(fruit, the_tag);
}

The function eat at the bottom first declares an instance of the tag. If an apple is entered in this function, it will be an instance of the apple_tag. If a banana is entered, it is a banana_tag. Then it forwards its call to the eat function in the namespace dispatch. There are two versions (overloads) there, one specified with an apple_tag, one with a banana_tag. The compiler selects, based on the tag the right function. Quite easy, and quite powerful, this tag dispatching system.

As said, within Boost.Geometry it is not applied like this. One of the reasons is that the instance of the tag is not necessary at all. If the dispatch::eat function would be implemented within a struct, as a static method, the tag dispatching can be done by type and not by instance. This is still tag dispatching!

So the piece above can be replaced by this Tag Dispatching by Type:

namespace dispatch
{
    template <typename Tag> struct eat {};

    template <> struct eat<apple_tag>
    {
        static void apply(apple const& a)
        {
          std::cout << "bite" << std::endl;
        }
    };

    template <> struct eat<banana_tag>
    {
        static void apply(banana const& b)
        {
          std::cout << "peel" << std::endl;
        }
    };
}

template <typename T>
void eat(T const& fruit)
{
    dispatch::eat<typename tag<T>::type>::apply(fruit);
}

So no instance at all, type only. Boost.Geometry has this structure everywhere, and all static dispatching methods are called apply.

For completeness we list the main program as well here:

int main()
{
    apple a("my apple");
    banana b("my banana");
    eat(a);
    eat(b);

    return 0;
}

It is the same for the instance-version and the type-version.

More advantages of the Tag Dispatching by Type approach are that arguments can be reversed, so you can call distance(point, polygon) and distance(polygon, point) but there is only one implementation with point_tag, polygon_tag necessary. With instances (and function overloads) this is not possible, or at least not that easy. Furthermore, tag dispatching by type can be used to define types, or define constants, as well. Suppose you need to define a constant value dependant on a type. We approach geometry now and want to define if a specific fruit-type is spherical or not. So we define this structure, specialized by tag:

    template <typename Tag> struct spherical {};

    template <> struct spherical<apple_tag>
    {
        static const bool value = true;
    };

    template <> struct spherical<banana_tag>
    {
        static const bool value = false;
    };

We can create a generic structure, which calls this dispatch, just like the free function eat, but now with a struct (or meta-function):

template <typename T>
struct spherical
{
    static const bool value = dispatch::spherical<typename tag<T>::type>::value;
};

Next blog will handle Tag Dispatching and Inheritance.

Efficient distances

Square root

Every programmer knows that sqrt (square root) is an expensive function. There are sometimes rumours that it is nowadays efficient enough. But that is not always the case. See e.g. this link.

It can easily be tested and my personal conclusions are that sqrt is still slow. See e.g. this test-program I created. You can run it even remotely in codepad, showing sqrt is about 3 slower than addition. On my computer the factor is more than 8.  (Normally I use Boost.Timer for timings, but for codepad the clock() function will do).

So sqrt is relatively slow. If you calculate many (thousands) of them, and you can avoid them, avoid them, it will increase performance. In distances sqrt is often easy to avoid.

Cartesian distance

A distance in a Cartesian system is normally calculated using the Pythagorean theorem, by the well-known function:  d = √(Δx² + Δy²). So that function contains an expensive square root.

If you compare distances the square root can easily be avoided, because √(Δx²+Δy²) < d is equivalent to Δx²+Δy² < d² [Paul Hsieh's Square Root page].
Stated differently, the function sqrt is strictly increasing. Therefore, if (and only if) √a < √b then a²< b² (a and b are positive).

Avoiding the square root is a well-known performance trick in geometry. For example, in calculating the distance from a point to a linestring or a polygon, many distances are compared, sometimes 10000 or 100000 times. The smallest distance is selected. There is only one square root calculation necessary: the final result. The closest distance of a point to all polygon segments is returned as the distance, and there a (only) square root is applied.

Comparable distance

The Boost.Geometry library has a comparable_distance function. It returns a distance-measure, so not a real distance, but something related to distance. In the case of a Cartesian coordinate system, it returns the square distance.

So library users can use the distance function (returning the real distance) and the comparable_distance function (returning the distance-measure which is comparable to other distance measures). The comparable_distance is (of course) more efficient. Internally Boost.Geometry uses comparable distance measures where possible. For example in the distance of a point to a polygon (see above), but also in simplification (removing less important points from a geometry).

Great circle distance

Besides the distance between two points in metric space, often the distance on a sphere is necessary. This cannot be calculated using the Pythagorean theorem. Distances on a sphere are called great circle distances, because the shortest path between two points on a sphere is called the great circle. For a perfect sphere, the following formula is applicable:

 d = 2 * asin(sqrt((sin((lat1-lat2)/2))^2 
                    + cos(lat1) * cos(lat2) 
                      * (sin((lon1-lon2)/2))^2))

Substituting sin(a/2)² with the haversine function hav makes it somewhat simpler:

 d = 2 * asin(sqrt(hav(lat1-lat2) 
                   + cos(lat1) * cos(lat2) 
                     * hav(lon1-lon2) ))

This distance is usually multiplied by the radius of the earth (R ≈ 6378.137 km), to get the distance in kilometers (use other units to get it in e.g. miles).

So let's try, for fun:

Amsterdam: 52° 22′ 23″ N, 4° 53′ 32″ E => 52.373056, 4.892222 (degrees)
Barcelona: 41° 23′ 0″ N, 2° 11′ 0″ E => 41.383333, 2.183333 (degrees)

 d = R * 2 * asin(sqrt(hav(52.373056 - 41.383333) 
                       + cos(52.373056) * cos(41.383333) 
                         * hav(2.183333 - 4.892222) ))
 = R * 2 * asin(sqrt(0.0091693 
                     + 0.61051768 * 0.7503034 
                       * 5.58722e-4))
 = R * 2 * asin(sqrt(0.00942524))
 = R * 2 * 0.194473661
 = 1240.38 kilometers

Which is about right. Note that we still assume a perfect sphere, there are more precise calculations for the ellipsoidal Earth.

Comparable great circle distance

What is less widely known, but becomes clear now, is that this function also has a much faster comparable equivalent. Of course we can save multiplying by R * 2, because if a > b, then R * 2 * a > R * 2 * b. But the asin (arc sinus) function is also strictly monotone increasing. So if a > b then asin(a) > asin(b). Therefore, we can save that calculation too, and we then of course can also save the square root again for this case. What's left are four trigoniometric functions. Which are still expensive, but cannot be avoided here.

So the comparable variant of this function is just

 cd = hav(lat1-lat2) + cos(lat1) * cos(lat2) * hav(lon1-lon2)

Comparing the performance of the great dircle distance and its comparable variant gives, on my computer, a factor 37. So to compare great circle distances (still assuming a perfect sphere) it is enough to compare the distance using the formula above, and it is 37 times faster. This will save a lot when calculating the distance of a GPS-point to a polygon (city, community, country) specified in lat/long coordinates.

Boost.Geometry implements this comparable haversine now.

Boost.Build and bjam

Just to document the essentials.
I had some difficulties in the past using bjam in combination with Boost, and had them again with each new non-Linux platform. The Boost SVN sandbox shows that many Boost authors encounter the same problems, resulting in a variety of creative solutions.

The canonical way to build a project using boost header-only libraries, making with bjam (Boost.Build) is:

• add a clause "<dependency>/boost//headers" into the jamfile (either Jamfile.v2, or Jamfile.jam, or Jamroot.jam)
• add a clause "use-project boost : c:/software/boost_1_43_0 ;" (adapt path if necessary) into the file "user-config.jam" which lives in c:\software\boost_1_43_0\tools\build\v2\user-config.jam
  
• bjam.exe should be in your path
• how does it find Boost.Build? There are two ways to do this:
  • either add a file named 'boost-build.jam' in your directory tree, either at the level of 'Jamroot', or a level higher
  • or set the environment variable BOOST_BUILD_PATH e.g. using export BOOST_BUILD_PATH='c:\software\boost_1_43_0\tools\build\v2' (note that 1) setting environment variables might be inconvenient and 2) if both are set, the boost-build.jam file setting is preferred)

Then your project should compile without errors. For example, this is your source-file "hello.cpp":

#include <iostream>
#include <boost/algorithm/string.hpp>

int main()
{
    std::cout
        << boost::to_upper_copy(std::string("Hello boost!"))
        << std::endl;
    return 0;
}

And this is your jamfile "jamroot.jam":

exe hello
    : hello.cpp
    : <dependency>/boost//headers
    ;

And this is your file "boost-build.jam":

boost-build c:/software/boost_1_43_0/tools/build/v2 ;

And this is at the bottom of your "user-config.jam":

use-project boost : c:/software/boost_1_43_0 ;

This works in at least MSVC (command prompt) and MinGW using GCC

• Thanks to the Boost-List
• See also this link for a well written Boost.Build tutorial
• See also this link with the Boost.Build documentation
• See also this link with a Wiki about Boost.Build