Return to Answer

maths http://petiteleve.free.fr/SO/sweep.png

_{(source: free.fr)}

On my modest PC it seems to be able to twist over 4 million snakes per minute, but just a few reptiles seem to manage to crawl under the radar.

#include <cmath>
#include <cassert>
#include <ctime>
#include <sstream>
#include <vector>
#include <map>
#include <unordered_set>
#include <algorithm> // sort

using namespace std;

// theroretical max snake lenght (but the Sun would burn out before you'll reach that)
#define MAX_LENGTH ((int)(1+8*sizeof(unsigned)))
    
#ifndef QT_DEBUG // using Qt IDE for g++ builds
#ifndef NDEBUG
#define NDEBUG
#endif
#endif

#ifdef NDEBUG
inline void tprintf(const char *, ...){}
#else
#define tprintf printf
#endif

void panic(const char * msg)
{
    printf("PANIC: %s\n", msg);
    exit(-1);
}

// ============================================================================
// fast bit reversal
// ============================================================================
unsigned bit_reverse(register unsigned x, unsigned len)
{
    x = (((x & 0xaaaaaaaa) >> 1) | ((x & 0x55555555) << 1));
    x = (((x & 0xcccccccc) >> 2) | ((x & 0x33333333) << 2));
    x = (((x & 0xf0f0f0f0) >> 4) | ((x & 0x0f0f0f0f) << 4));
    x = (((x & 0xff00ff00) >> 8) | ((x & 0x00ff00ff) << 8));
    return((x >> 16) | (x << 16)) >> (32-len);
}

// ============================================================================
// 2D geometry (restricted to integer coordinates and right angle rotations)
// ============================================================================

// points using integer- or float-valued coordinates
template<typename T>struct tTypedPoint;

typedef int    tCoord;
typedef double tFloatCoord;

typedef tTypedPoint<tCoord> tPoint;
typedef tTypedPoint<tFloatCoord>  tFloatPoint;

template <typename T>
struct tTypedPoint {
    T x, y;

    template<typename U> tTypedPoint(const tTypedPoint<U>& from) : x((T)from.x), y((T)from.y) {} // conversion constructor

    tTypedPoint() {}
    tTypedPoint(T x, T y) : x(x), y(y) {}
    tTypedPoint(const tTypedPoint & p) { *this = p; }
    tTypedPoint operator+ (const tTypedPoint & p) const { return{ x + p.x, y + p.y }; }
    tTypedPoint operator- (const tTypedPoint & p) const { return{ x - p.x, y - p.y }; }
    tTypedPoint operator* (T scalar) const { return{ x * scalar, y * scalar }; }
    tTypedPoint operator/ (T scalar) const { return{ x / scalar, y / scalar }; }
    bool operator== (const tTypedPoint & p) const { return x == p.x && y == p.y; }
    bool operator!= (const tTypedPoint & p) const { return !operator==(p); }
    T dot(const tTypedPoint &p) const { return x*p.x + y * p.y; } // dot product  
    int cross(const tTypedPoint &p) const { return x*p.y - y * p.x; } // z component of cross product
    T norm2(void) const { return dot(*this); }

    // works only with direction = 1 (90° right) or -1 (90° left)
    tTypedPoint rotate(int direction) const { return{ direction * y, -direction * x }; }
    tTypedPoint rotate(int direction, const tTypedPoint & center) const { return (*this - center).rotate(direction) + center; }

    // used to compute the length of a ragdoll snake segment
    unsigned manhattan_distance(const tPoint & p) const { return abs(x-p.x) + abs(y-p.y); }
};


struct tArc {
    tPoint c;                        // circle center
    tFloatPoint middle_vector;       // vector splitting the arc in half
    tCoord      middle_vector_norm2; // precomputed for speed
    tFloatCoord dp_limit;            // dot product of point with the middle vector 
                                     // will be >= dp_limit if the point is inside the arc

    tArc() {}
    tArc(tPoint c, tPoint p, int direction) : c(c)
    {
        tPoint r = p - c;
        tPoint end = r.rotate(direction);
        middle_vector = ((tFloatPoint)(r+end)) / sqrt(2);
        middle_vector_norm2 = r.norm2();
        dp_limit = r.dot(middle_vector);
    }

    bool contains(tFloatPoint p) // p must be a point on the circle
    {
        if ((p-c).dot(middle_vector) >= dp_limit)
        {
            return true;
        }
        else return false;
    }
};

// returns the point of line (p1 p2) that is closest to c
// handles degenerate case p1 = p2
tPoint line_closest_point(tPoint p1, tPoint p2, tPoint c)
{
    if (p1 == p2) return{ p1.x, p1.y };
    tPoint p1p2 = p2 - p1;
    tPoint p1c =  c  - p1;
    return p1 + p1p2 * (p1c.dot(p1p2) / p1p2.norm2());
}

// variant of closest point computation that checks if the projection falls within the segment
bool closest_point_within(tPoint p1, tPoint p2, tPoint c, tPoint & res)
{
    tPoint p1p2 = p2 - p1;
    tPoint p1c = c - p1;
    tCoord nk = p1c.dot(p1p2);
    if (nk <= 0) return false;
    tCoord n = p1p2.norm2();
    if (nk >= n) return false;
    res = p1 + p1p2 * (nk / n);
    return true;
}

// tests intersection of line (p1 p2) with an arc
bool inter_seg_arc(tPoint p1, tPoint p2, tArc arc)
{
    tPoint m = line_closest_point(p1, p2, arc.c);
    tCoord r2 = arc.middle_vector_norm2;
    tPoint cm = m - arc.c;
    tCoord h2 = cm.norm2();
    if (r2 < h2) return false; // no circle intersection

    tPoint p1p2 = p2 - p1;
    tCoord n2p1p2 = p1p2.norm2();

    // works because by construction p is on (p1 p2)
    auto in_segment = [&](const tFloatPoint & p) -> bool
    {
        tFloatCoord nk = p1p2.dot(p - p1);
        return nk >= 0 && nk <= n2p1p2;
    };

    if (r2 == h2) return arc.contains(m) && in_segment(m); // tangent intersection

    //if (p1 == p2) return false; // degenerate segment located inside circle
    assert(p1 != p2);

    tFloatPoint u = (tFloatPoint)p1p2 * sqrt((r2-h2)/n2p1p2); // displacement on (p1 p2) from m to one intersection point

    tFloatPoint i1 = m + u;
    if    (arc.contains(i1) && in_segment(i1)) return true;
    tFloatPoint i2 = m - u;
    return arc.contains(i2) && in_segment(i2);
}

// test intersection of segments [a1 b1] and [a2 b2]
// segments may degenerate to points
bool inter_seg_seg(tPoint a1, tPoint b1, tPoint a2, tPoint b2)
{
    tPoint u1 = b1 - a1;
    tPoint u2 = b2 - a2;
    tCoord det = u1.cross(u2);
    tPoint link = a2 - a1;
    if (det == 0) // parallel
    {
        if (link.cross(u1) != 0) return false; // disjoint

        // 1st range bounds are not ordered while 2nd range's are
        auto ranges_overlap = [](tCoord a1, tCoord a2, tCoord bmin, tCoord bmax)
        {
            return (a1 <= a2)
                ? (a1 <= bmax && a2 >= bmin)
                : (a2 <= bmax && a1 >= bmin);
        };

        tCoord k1 = u1.dot(link);
        tCoord k2 = u1.dot(b2 - a1);
        return ranges_overlap(k1, k2, 0, u1.norm2());
    }
    else // secant
    {
        tCoord det2 = det*det;
        tCoord k1 = link.cross(u1) * det;
        if (k1 < 0 || k1 > det2) return false;
        tCoord k2 = link.cross(u2) * det;
        return k2 >= 0 && k2 <= det2;
    }
}

// ============================================================================
// partitions of a snake into segments
// ============================================================================
class sPartition {
    vector<int>segment;

public:
    sPartition() {} // only used in vector inits

    // split a flat snake into straight segments
    sPartition(unsigned generator, unsigned num_joints) 
    {
        int segment_len = 1;
        for (size_t i = 0; i != num_joints; i++)
        {
            if (generator & 1)
            {
                segment.push_back(segment_len);
                segment_len = 1;
            }
            else
            {
                segment_len++;
            }
            generator >>= 1;
        }
        segment.push_back(segment_len);
    }

    // direct access to underlying array of segment lengths
    size_t size(void) const { return segment.size(); }
    int & operator[] (size_t i) { return segment[i]; }
    int operator[] (size_t i) const { return segment[i]; }
};

// ============================================================================
// compact storage of a configuration (64 bits)
// ============================================================================
struct sConfiguration {
    unsigned partition;
    unsigned folding;

    sConfiguration(unsigned partition, unsigned folding) : partition(partition), folding(folding) {}

    // textual representation
    string text(unsigned length)
    {
        ostringstream res;

        unsigned f = folding;
        unsigned p = partition;

        int segment_len = 1;
        int direction = 1;
        for (size_t i = 1; i != length; i++)
        {
            if (p & 1)
            {
                res << segment_len * direction << ',';
                direction = (f & 1) ? -1 : 1;
                segment_len = 1;
            }
            else segment_len++;

            p >>= 1;
            f >>= 1;
        }
        res << segment_len * direction;
        return res.str();
    }

    // for final sorting
    bool operator< (const sConfiguration& c) const
    {
        return (partition == c.partition) ? folding < c.folding : partition < c.partition;
    }

    // hash table setup
    struct Hasher {
        size_t operator() (const sConfiguration& c) const
        {
            return c.partition ^ c.folding;
        }
    };
    struct Equalizer {
        bool operator() (const sConfiguration& c1, const sConfiguration& c2) const
        {
            return c1.partition == c2.partition && c1.folding == c2.folding;
        }
    };
    typedef unordered_set<sConfiguration, Hasher, Equalizer> set;

};

// ============================================================================
// snake ragdoll 
// ============================================================================
class tSnakeDoll {
    vector<tPoint>point; // snake geometry. Head at (0,0) pointing right

    // allows to check for collision with the area swept by a rotating segment
    struct rotatedSegment {
        struct segment { tPoint a, b; };
        segment org, end;
        tArc    arc[3];
        bool extra_arc; // see if third arc is needed

        // empty constructor to avoid wasting time in vector initializations
        rotatedSegment(){}
        // copy constructor is mandatory for vectors *but* shall never be used, since we carefully pre-allocate vector memory
        rotatedSegment(const rotatedSegment &){ assert(!"rotatedSegment should never have been copy-constructed"); }

        // rotate a segment
        rotatedSegment(tPoint pivot, int rotation, tPoint o1, tPoint o2)
        {
            arc[0] = tArc(pivot, o1, rotation);
            arc[1] = tArc(pivot, o2, rotation);
            tPoint middle;
            extra_arc = closest_point_within(o1, o2, pivot, middle);
            if (extra_arc) arc[2] = tArc(pivot, middle, rotation);
            org = { o1, o2 };
            end = { o1.rotate(rotation, pivot), o2.rotate(rotation, pivot) };
        }

        // check if a segment intersects the area swept during rotation
        bool intersects(tPoint p1, tPoint p2) const
        {
            auto print_seg = [&](const segment & s) { tprintf("(%d,%d)(%d,%d) -> (%d,%d)(%d,%d)", p1.x, p1.y, p2.x, p2.y, s.a.x, s.a.y, s.b.x, s.b.y); };
            auto print_arc = [&](int a) { tprintf("(%d,%d)(%d,%d) -> %d (%d,%d)[%f,%f]", p1.x, p1.y, p2.x, p2.y, a, arc[a].c.x, arc[a].c.y, arc[a].middle_vector.x, arc[a].middle_vector.y); };

            if (p1 == org.a) return false; // pivot is the only point allowed to intersect
            if (inter_seg_seg(p1, p2, org.a, org.b)) { print_seg(org);  return true; }
            if (inter_seg_seg(p1, p2, end.a, end.b)) { print_seg(end);  return true; }
            if (inter_seg_arc(p1, p2, arc[0])) { print_arc(0);  return true; }
            if (inter_seg_arc(p1, p2, arc[1])) { print_arc(1);  return true; }
            if (extra_arc && inter_seg_arc(p1, p2, arc[2])) { print_arc(2);  return true; }
            return false;
        }
    };

public:
    unsigned folding; // compressed folding list
    int level; // ;!;

    // holds results of a folding attempt
    class snakeFolding {
        friend class tSnakeDoll;
        vector<rotatedSegment>segment; // rotated segments
        unsigned joint;
        int direction;
        size_t i_rotate;

        // pre-allocate rotated segments
        void reserve(size_t length)
        {
            segment.clear(); // this supposedly does not release vector storage memory
            segment.reserve(length);
        }

        // handle one segment rotation
        void rotate(tPoint pivot, int rotation, tPoint o1, tPoint o2)
        {
            segment.emplace_back(pivot, rotation, o1, o2);
        }
    public:
        // nothing done during construction
        snakeFolding() {};
    };

    // empty default constructor to avoid wasting time in array/vector inits
    tSnakeDoll() {}

    // constructs ragdoll from compressed configuration
    tSnakeDoll(unsigned size, unsigned generator, unsigned folding) : point(size), folding(folding)
    {
        tPoint direction(1, 0);
        tPoint current = { 0, 0 };
        size_t p = 0;
        point[p++] = current;
        for (size_t i = 1; i != size; i++)
        {
            current = current + direction;
            if (generator & 1)
            {
                direction.rotate((folding & 1) ? -1 : 1);
                point[p++] = current;
            }
            folding >>= 1;
            generator >>= 1;
        }
        point[p++] = current;
        point.resize(p);
    }

    // constructs the initial flat snake
    tSnakeDoll(int size) : point(2), folding(0)
    {
        point[0] = { 0, 0 };
        point[1] = { size, 0 };
        level = 0;
    }

    // constructs a new folding with one added rotation
    tSnakeDoll(const tSnakeDoll & start, const snakeFolding & fold)
    {
        point.resize(start.point.size()+1);
        size_t i;

        // copy unmoved points
        for (i = 0; i != fold.i_rotate; i++) point[i] = start.point[i];

        // copy rotated points
        for (; i != start.point.size(); i++) point[i] = fold.segment[i - fold.i_rotate].end.a;
        point[i] = fold.segment[i - 1 - fold.i_rotate].end.b;

        // update folding bitfield
        folding = start.folding;
        if (fold.direction == -1) folding |= (1 << fold.joint);
        level = start.level + 1;
//tprintf("snake doll level %d fold %X\n", level, folding);
    }

    // prepare one folding at a specified joint
    // rotation is 1 for right, -1 for left
    bool fold_at(unsigned joint, int rotation, snakeFolding & res) const
    {
        // first bend is only on the right
        if (level == 0 && rotation == -1)
        {
            tprintf("1st left ");
            return false;
        }

        res.joint = joint;
        res.direction = rotation;

        // locate folding point
        joint++;
        tPoint pivot;
        size_t i_rotate = 0;
        for (size_t i = 1 ; i != point.size(); i++)
        {
            unsigned len = point[i].manhattan_distance(point[i-1]);
            if (len > joint)
            {
                pivot = point[i-1] + ((point[i] - point[i - 1]) / len) * joint;
                i_rotate = i;
                break;
            }
            else joint -= len;
        }
        res.i_rotate = i_rotate;

        // rotate around joint
        res.reserve(point.size() - i_rotate); // number of rotated segments
        res.rotate(pivot, rotation, pivot, point[i_rotate]);
        for (size_t i = i_rotate + 1; i != point.size(); i++) res.rotate (pivot, rotation, point[i-1], point[i]);

        // check collisions
        for (const rotatedSegment & fold : res.segment)
        {
            for (size_t i = 1; i != i_rotate; i++) if (fold.intersects(point[i], point[i - 1])) return false;
            if (fold.intersects(pivot, point[i_rotate - 1])) return false;
        }
        return true;
    }
};

// ============================================================================
// snake twisting engine
// ============================================================================
class cSnakeFolder {
    int length;
    unsigned num_joints;
    unsigned reverse_start;

    map<unsigned,sPartition>partitions;
    sConfiguration::set config_hash;

    // recursive folding
    void fold(unsigned generator, tSnakeDoll snake, unsigned first_joint)
    {
        // filter redundant partitions
        if (partitions.count(generator) == 0)
        {
            tprintf("%s %s P cut\n", string(snake.level, '|').c_str(), sConfiguration(generator, snake.folding).text(length).c_str());
            return;
        }

        // filter redundant foldings
        unsigned reverse = bit_reverse(snake.folding, num_joints);
        if (reverse & 1) reverse = ~reverse & ((1<<num_joints)-1); // invert only significant bits
        if (config_hash.count({ generator, reverse }))
        {
            tprintf("%s %s F cut\n", string(snake.level, '|').c_str(), sConfiguration(generator, snake.folding).text(length).c_str());
            return;
        }

        // store this configuration
        config_hash.emplace(generator, snake.folding);

        // try to bend remaining joints
        for (size_t joint = first_joint; joint != num_joints; joint++)
        {
            unsigned next_partition = generator | (1 << joint);
            tSnakeDoll::snakeFolding bend;

            tprintf("%s %s -> %s ", string(snake.level, '|').c_str(), sConfiguration(generator, snake.folding).text(length).c_str(), sConfiguration(next_partition, snake.folding).text(length).c_str());
            if (snake.fold_at(joint, 1, bend))
            {
                tprintf("ok\n");
                fold(next_partition, tSnakeDoll(snake, bend), joint + 1);
            }
            else tprintf("failed\n");

            tprintf("%s %s -> %s ", string(snake.level, '|').c_str(), sConfiguration(generator, snake.folding).text(length).c_str(), sConfiguration(next_partition, snake.folding | (1 << joint)).text(length).c_str());
            if (snake.fold_at(joint, -1, bend))
            {
                tprintf("ok\n");
                fold(next_partition, tSnakeDoll(snake, bend), joint + 1);
            }
            else tprintf("failed\n");
        }
    }

public:
    // all found configurations, in compressed format
    vector<sConfiguration>configurations;

    // constructor does all the job
    cSnakeFolder(int n) : length(n)
    {
        num_joints = length - 1;
        reverse_start = 1 << num_joints;

        // generate snake partitions
        unsigned num_part = 1 << num_joints;
        for (unsigned partition = 0; partition != num_part; partition++)
        {
            // filter out symetric partitions
            if (partitions.count(bit_reverse(partition, num_joints))) continue;

            // generate and store partition
            partitions[partition] = sPartition (partition, num_joints);
        }

        // launch recursive folding
        fold(0, tSnakeDoll(length), 0);

        // collect configurations
        for (auto c : config_hash) configurations.emplace_back(c);
        sort(configurations.begin(), configurations.end());
        config_hash.clear(); // that's a huge lump of memory taken care of
    }

    // get a ragdoll representing a particular configuration
    /*
       For a snake of length N, the ragdoll is a list of N+1 points defining the N unitary segments.
       By convention, the snake head starts at coordinates (0,0) pointing to the right.
    */
    tSnakeDoll doll(size_t conf)
    {
        if (conf >= configurations.size()) conf = 0; // out of bounds queries will get a flat snake
        sConfiguration & c = configurations[conf];
        return tSnakeDoll(length, c.partition, c.folding);
    }

    // get text description of a particular configuration
    /* 
       Absolute values are segment lengths
       Sign indicate a preceding turn (by convention, positive for right, negative for left)
       1st segment is always positive by convention.
       Snake length is the sum of segment lengths absolute values.

           ___   For instance 2,-1,3 represents a snake of length 6 bent as
        __|      2 segments, left turn, 1 segment, right turn, 3 segments
    */
    string text(size_t conf)
    {
        if (conf >= configurations.size()) return "**out of bounds**";
        return configurations[conf].text(length);
    }
};

// ============================================================================
// here we go
// ============================================================================
int main(int argc, char * argv[])
{
#ifdef NDEBUG
    if (argc != 2) panic("give me a snake length or else");
    int length = atoi(argv[1]);
#else
    (void)argc; (void)argv;
    int length = 5;
#endif // NDEBUG

    if (length <= 0 || length > MAX_LENGTH) panic("a snake of that length is hardly foldable");

    time_t start = time(NULL);
    cSnakeFolder snakes(length);
    time_t duration = time(NULL) - start;
#ifndef NDEBUG
    for (size_t i = 0; i != snakes.configurations.size(); i++) printf("%3d: %s\n", i, snakes.text(i).c_str());
#endif // NDEBUG
    printf ("Found %d configuration%c of length %d in %lds\n", snakes.configurations.size(), (snakes.configurations.size() == 1) ? '\0' : 's', length, duration);
    return 0;
}

It is my second computation-intensive challenge using STL new-ish hash tables (unordered_set), and again Microsoft implementation performs poorly, to say the least.

g++ build is 2.5 times faster and uses about 15% less memory.

The g++ build uses about 32 bytes/configuration.
A snake configuration is 8 bytes, the remaining 24 are mostly used by hash table internals.
For convenience I transfer the results from the hash table into a sorted vector in the end, which increases memory consumption significantly (from 32 to 40 bytes/configuration).
Even so, the program can compute up to n=19 before hitting the Win32 2Gb allocation limit, and compiling for x64 would allow to go a couple of steps beyond.

Since collision check is roughly in O(n), computation time should be in O(nkⁿ), with k slightly lower than 3.
On my [email protected], n=19 takes about 4 minutes.

This means memory should be less an issue than speed.

Clearly collision check will be the most time-consuming part, so the main avenue for optimization is to reduce these computations.

Collision check

When a snake folds around one joint, each rotated segment will sweep an area whose shape is anything but trivial.
Clearly you can check collisions by testing inclusion within all such swept areas individually, though a more global check would be more efficient.

After an interesting discussion with trichoplax and a bit of JavaScript to get my bearings, I devised a method that should enable to get an analytic description of the swept area under certain conditions.

For any segment that does not contain I, the swept area is bound by 2 arcs and 2 segments.
If I falls within the segment, a third arc must be taken into account.

This means we can check each unmoving segment against each rotating segment with

• 2 segment-with-segment intersections
• 2 or 3 segment-with-arc intersections

I used vector geometry to avoid trigonometric functions altogether.
Vector operations produce compact and (relatively) readable code.
Better than most of the matrix crunching stuff I found on MathWorld or many gaming sites, anyway.

This allows to compute segment-with-segment intersections with integer values.

Why does it not work?

I am not done tweaking the collision code yet, but it does not seem to explain why a very few configurations are not detected.
For n=5 I miss 2 snakes though collision detection works as planned.

Most likely the cuts are just a tad bit too selective, but for now I am not able to pinpoint the logic error.

After reading this article, I collected bits of wisdom from that guy who apparently worked for 25 years on the less complicated problem of counting self-avoiding paths on a square lattice.

#include <cassert>
#include <ctime>
#include <sstream>
#include <vector>
#include <algorithm> // sort

using namespace std;

// theroretical max snake lenght (the code would need a few decades to process that value)
#define MAX_LENGTH ((int)(1+8*sizeof(unsigned)))

#ifndef _MSC_VER
#ifndef QT_DEBUG // using Qt IDE for g++ builds
#define NDEBUG
#endif
#endif

#ifdef NDEBUG
inline void tprintf(const char *, ...){}
#else
#define tprintf printf
#endif

void panic(const char * msg)
{
    printf("PANIC: %s\n", msg);
    exit(-1);
}

// ============================================================================
// fast bit reversal
// ============================================================================
unsigned bit_reverse(register unsigned x, unsigned len)
{
    x = (((x & 0xaaaaaaaa) >> 1) | ((x & 0x55555555) << 1));
    x = (((x & 0xcccccccc) >> 2) | ((x & 0x33333333) << 2));
    x = (((x & 0xf0f0f0f0) >> 4) | ((x & 0x0f0f0f0f) << 4));
    x = (((x & 0xff00ff00) >> 8) | ((x & 0x00ff00ff) << 8));
    return((x >> 16) | (x << 16)) >> (32-len);
}

// ============================================================================
// 2D geometry (restricted to integer coordinates and right angle rotations)
// ============================================================================

// points using integer- or float-valued coordinates
template<typename T>struct tTypedPoint;

typedef int    tCoord;
typedef double tFloatCoord;

typedef tTypedPoint<tCoord> tPoint;
typedef tTypedPoint<tFloatCoord>  tFloatPoint;

template <typename T>
struct tTypedPoint {
    T x, y;

    template<typename U> tTypedPoint(const tTypedPoint<U>& from) : x((T)from.x), y((T)from.y) {} // conversion constructor

    tTypedPoint() {}
    tTypedPoint(T x, T y) : x(x), y(y) {}
    tTypedPoint(const tTypedPoint& p) { *this = p; }
    tTypedPoint operator+ (const tTypedPoint & p) const { return{ x + p.x, y + p.y }; }
    tTypedPoint operator- (const tTypedPoint & p) const { return{ x - p.x, y - p.y }; }
    tTypedPoint operator* (T scalar) const { return{ x * scalar, y * scalar }; }
    tTypedPoint operator/ (T scalar) const { return{ x / scalar, y / scalar }; }
    bool operator== (const tTypedPoint & p) const { return x == p.x && y == p.y; }
    bool operator!= (const tTypedPoint & p) const { return !operator==(p); }
    T dot(const tTypedPoint &p) const { return x*p.x + y * p.y; } // dot product  
    int cross(const tTypedPoint &p) const { return x*p.y - y * p.x; } // z component of cross product
    T norm2(void) const { return dot(*this); }

    // works only with direction = 1 (90° right) or -1 (90° left)
    tTypedPoint rotate(int direction) const { return{ direction * y, -direction * x }; }
    tTypedPoint rotate(int direction, const tTypedPoint & center) const { return (*this - center).rotate(direction) + center; }

    // used to compute length of a ragdoll snake segment
    unsigned manhattan_distance(const tPoint & p) const { return abs(x-p.x) + abs(y-p.y); }
};


struct tArc {
    tPoint c;                        // circle center
    tFloatPoint middle_vector;       // vector splitting the arc in half
    tCoord      middle_vector_norm2; // precomputed for speed
    tFloatCoord dp_limit;

    tArc() {}
    tArc(tPoint c, tPoint p, int direction) : c(c)
    {
        tPoint r = p - c;
        tPoint end = r.rotate(direction);
        middle_vector = ((tFloatPoint)(r+end)) / sqrt(2); // works only for +-90° rotations. The vector should be normalized to circle radius in the general case
        middle_vector_norm2 = r.norm2();
        dp_limit = ((tFloatPoint)r).dot(middle_vector);
        assert (middle_vector == tPoint(0, 0) || dp_limit != 0);
    }

    bool contains(tFloatPoint p) // p must be a point on the circle
    {
        if ((p-c).dot(middle_vector) >= dp_limit)
        {
            return true;
        }
        else return false;
    }
};

// returns the point of line (p1 p2) that is closest to c
// handles degenerate case p1 = p2
tPoint line_closest_point(tPoint p1, tPoint p2, tPoint c)
{
    if (p1 == p2) return{ p1.x, p1.y };
    tPoint p1p2 = p2 - p1;
    tPoint p1c =  c  - p1;
    tPoint disp = (p1p2 * p1c.dot(p1p2)) / p1p2.norm2();
    return p1 + disp;
}

// variant of closest point computation that checks if the projection falls within the segment
bool closest_point_within(tPoint p1, tPoint p2, tPoint c, tPoint & res)
{
    tPoint p1p2 = p2 - p1;
    tPoint p1c = c - p1;
    tCoord nk = p1c.dot(p1p2);
    if (nk <= 0) return false;
    tCoord n = p1p2.norm2();
    if (nk >= n) return false;
    res = p1 + p1p2 * (nk / n);
    return true;
}

// tests intersection of line (p1 p2) with an arc
bool inter_seg_arc(tPoint p1, tPoint p2, tArc arc)
{
    tPoint m = line_closest_point(p1, p2, arc.c);
    tCoord r2 = arc.middle_vector_norm2;
    tPoint cm = m - arc.c;
    tCoord h2 = cm.norm2();
    if (r2 < h2) return false; // no circle intersection

    tPoint p1p2 = p2 - p1;
    tCoord n2p1p2 = p1p2.norm2();

    // works because by construction p is on (p1 p2)
    auto in_segment = [&](const tFloatPoint & p) -> bool
    {
        tFloatCoord nk = p1p2.dot(p - p1);
        return nk >= 0 && nk <= n2p1p2;
    };

    if (r2 == h2) return arc.contains(m) && in_segment(m); // tangent intersection

    //if (p1 == p2) return false; // degenerate segment located inside circle
    assert(p1 != p2);

    tFloatPoint u = (tFloatPoint)p1p2 * sqrt((r2-h2)/n2p1p2); // displacement on (p1 p2) from m to one intersection point

    tFloatPoint i1 = m + u;
    if    (arc.contains(i1) && in_segment(i1)) return true;
    tFloatPoint i2 = m - u;
    return arc.contains(i2) && in_segment(i2);
}

// ============================================================================
// compact storage of a configuration (64 bits)
// ============================================================================
struct sConfiguration {
    unsigned partition;
    unsigned folding;

    explicit sConfiguration() {}
    sConfiguration(unsigned partition, unsigned folding) : partition(partition), folding(folding) {}

    // add a bend
    sConfiguration bend(unsigned joint, int rotation) const
    {
        sConfiguration res;
        unsigned joint_mask = 1 << joint;
        res.partition = partition | joint_mask;
        res.folding = folding;
        if (rotation == -1) res.folding |= joint_mask;
        return res;
    }

    // textual representation
    string text(unsigned length) const
    {
        ostringstream res;

        unsigned f = folding;
        unsigned p = partition;

        int segment_len = 1;
        int direction = 1;
        for (size_t i = 1; i != length; i++)
        {
            if (p & 1)
            {
                res << segment_len * direction << ',';
                direction = (f & 1) ? -1 : 1;
                segment_len = 1;
            }
            else segment_len++;

            p >>= 1;
            f >>= 1;
        }
        res << segment_len * direction;
        return res.str();
    }

    // for final sorting
    bool operator< (const sConfiguration& c) const
    {
        return (partition == c.partition) ? folding < c.folding : partition < c.partition;
    }
};

// ============================================================================
// static snake geometry checking grid
// ============================================================================
typedef unsigned tConfId;

class tGrid {
    vector<tConfId>point;
    tConfId current;
    size_t snake_len;
    int min_x, max_x, min_y, max_y;
    size_t x_size, y_size;

    size_t raw_index(const tPoint& p) { bound_check(p);  return (p.x - min_x) + (p.y - min_y) * x_size; }
    void bound_check(const tPoint& p) const { assert(p.x >= min_x && p.x <= max_x && p.y >= min_y && p.y <= max_y); }

    void set(const tPoint& p)
    {
        point[raw_index(p)] = current;
    }
    bool check(const tPoint& p)
    {
        if (point[raw_index(p)] == current) return false;
        set(p);
        return true;
    } 

public:
    tGrid(int len) : current(-1), snake_len(len)
    {
        min_x = -max(len - 3, 0);
        max_x = max(len - 0, 0);
        min_y = -max(len - 1, 0);
        max_y = max(len - 4, 0);
        x_size = max_x - min_x + 1;
        y_size = max_y - min_y + 1;
        point.assign(x_size * y_size, current);
    }

    bool check(sConfiguration c)
    {
        current++;
        tPoint d(1, 0);
        tPoint p(0, 0);
        set(p);
        for (size_t i = 1; i != snake_len; i++)
        {
            p = p + d;
            if (!check(p)) return false;
            if (c.partition & 1) d = d.rotate((c.folding & 1) ? -1 : 1);
            c.folding >>= 1;
            c.partition >>= 1;
        }
        return check(p + d);
    }

};

// ============================================================================
// snake ragdoll 
// ============================================================================
class tSnakeDoll {
    vector<tPoint>point; // snake geometry. Head at (0,0) pointing right

    // allows to check for collision with the area swept by a rotating segment
    struct rotatedSegment {
        struct segment { tPoint a, b; };
        tPoint  org;
        segment end;
        tArc    arc[3];
        bool extra_arc; // see if third arc is needed

        // empty constructor to avoid wasting time in vector initializations
        rotatedSegment(){}
        // copy constructor is mandatory for vectors *but* shall never be used, since we carefully pre-allocate vector memory
        rotatedSegment(const rotatedSegment &){ assert(!"rotatedSegment should never have been copy-constructed"); }

        // rotate a segment
        rotatedSegment(tPoint pivot, int rotation, tPoint o1, tPoint o2)
        {
            arc[0] = tArc(pivot, o1, rotation);
            arc[1] = tArc(pivot, o2, rotation);
            tPoint middle;
            extra_arc = closest_point_within(o1, o2, pivot, middle);
            if (extra_arc) arc[2] = tArc(pivot, middle, rotation);
            org = o1;
            end = { o1.rotate(rotation, pivot), o2.rotate(rotation, pivot) };
        }

        // check if a segment intersects the area swept during rotation
        bool intersects(tPoint p1, tPoint p2) const
        {
            auto print_arc = [&](int a) { tprintf("(%d,%d)(%d,%d) -> %d (%d,%d)[%f,%f]", p1.x, p1.y, p2.x, p2.y, a, arc[a].c.x, arc[a].c.y, arc[a].middle_vector.x, arc[a].middle_vector.y); };

            if (p1 == org) return false; // pivot is the only point allowed to intersect
            if (inter_seg_arc(p1, p2, arc[0]))  
            {  
                print_arc(0);   
                return true;
            }
            if (inter_seg_arc(p1, p2, arc[1]))
            { 
                print_arc(1); 
                return true;
            }
            if (extra_arc && inter_seg_arc(p1, p2, arc[2])) 
            { 
                print_arc(2);
                return true;
            }
            return false;
        }
    };

public:
    sConfiguration configuration;
    bool valid;

    // holds results of a folding attempt
    class snakeFolding {
        friend class tSnakeDoll;
        vector<rotatedSegment>segment; // rotated segments
        unsigned joint;
        int direction;
        size_t i_rotate;

        // pre-allocate rotated segments
        void reserve(size_t length)
        {
            segment.clear(); // this supposedly does not release vector storage memory
            segment.reserve(length);
        }

        // handle one segment rotation
        void rotate(tPoint pivot, int rotation, tPoint o1, tPoint o2)
        {
            segment.emplace_back(pivot, rotation, o1, o2);
        }
    public:
        // nothing done during construction
        snakeFolding(unsigned size)
        {
            segment.reserve (size);
        }
    };

    // empty default constructor to avoid wasting time in array/vector inits
    tSnakeDoll() {}

    // constructs ragdoll from compressed configuration
    tSnakeDoll(unsigned size, unsigned generator, unsigned folding) : point(size), configuration(generator,folding)
    {
        tPoint direction(1, 0);
        tPoint current = { 0, 0 };
        size_t p = 0;
        point[p++] = current;
        for (size_t i = 1; i != size; i++)
        {
            current = current + direction;
            if (generator & 1)
            {
                direction.rotate((folding & 1) ? -1 : 1);
                point[p++] = current;
            }
            folding >>= 1;
            generator >>= 1;
        }
        point[p++] = current;
        point.resize(p);
    }

    // constructs the initial flat snake
    tSnakeDoll(int size) : point(2), configuration(0,0), valid(true)
    {
        point[0] = { 0, 0 };
        point[1] = { size, 0 };
    }

    // constructs a new folding with one added rotation
    tSnakeDoll(const tSnakeDoll & parent, unsigned joint, int rotation, tGrid& grid)
    {
        // update configuration
        configuration = parent.configuration.bend(joint, rotation);

        // locate folding point
        unsigned p_joint = joint+1;
        tPoint pivot;
        size_t i_rotate = 0;
        for (size_t i = 1; i != parent.point.size(); i++)
        {
            unsigned len = parent.point[i].manhattan_distance(parent.point[i - 1]);
            if (len > p_joint)
            {
                pivot = parent.point[i - 1] + ((parent.point[i] - parent.point[i - 1]) / len) * p_joint;
                i_rotate = i;
                break;
            }
            else p_joint -= len;
        }

        // rotate around joint
        snakeFolding fold (parent.point.size() - i_rotate);
        fold.rotate(pivot, rotation, pivot, parent.point[i_rotate]);
        for (size_t i = i_rotate + 1; i != parent.point.size(); i++) fold.rotate(pivot, rotation, parent.point[i - 1], parent.point[i]);

        // copy unmoved points
        point.resize(parent.point.size()+1);
        size_t i;
        for (i = 0; i != i_rotate; i++) point[i] = parent.point[i];

        // copy rotated points
        for (; i != parent.point.size(); i++) point[i] = fold.segment[i - i_rotate].end.a;
        point[i] = fold.segment[i - 1 - i_rotate].end.b;

        // static configuration check
        valid = grid.check (configuration);

        // check collisions with swept arcs
        if (valid && parent.valid) // ;!; parent.valid test is temporary
        {
            for (const rotatedSegment & s : fold.segment)
            for (size_t i = 0; i != i_rotate; i++)
            {
                if (s.intersects(point[i+1], point[i])) 
                {
                    //printf("! %s => %s\n", parent.trace().c_str(), trace().c_str());//;!;
                    valid = false;
                    break;
                }
            }
        }
    }

    // trace
    string trace(void) const
    {
        size_t len = 0;
        for (size_t i = 1; i != point.size(); i++) len += point[i - 1].manhattan_distance(point[i]);
        return configuration.text(len);
    }
};

// ============================================================================
// snake twisting engine
// ============================================================================
class cSnakeFolder {
    int length;
    unsigned num_joints;
    tGrid grid;

    // filter redundant configurations
    bool is_unique (sConfiguration c)
    {
        unsigned reverse_p = bit_reverse(c.partition, num_joints);
        if (reverse_p < c.partition)
        {
            tprintf("P cut %s\n", c.text(length).c_str());
            return false;
        }
        else if (reverse_p == c.partition) // filter redundant foldings
        {
            unsigned first_joint_mask = c.partition & (-c.partition); // insulates leftmost bit
            unsigned reverse_f = bit_reverse(c.folding, num_joints);
            if (reverse_f & first_joint_mask) reverse_f = ~reverse_f & c.partition;

            if (reverse_f > c.folding)
            {
                tprintf("F cut %s\n", c.text(length).c_str());
                return false;
            }
        }
        return true;
    }

    // recursive folding
    void fold(tSnakeDoll snake, unsigned first_joint)
    {
        // count unique configurations
        if (snake.valid && is_unique(snake.configuration)) num_configurations++;

        // try to bend remaining joints
        for (size_t joint = first_joint; joint != num_joints; joint++)
        {
            // right bend
            tprintf("%s -> %s\n", snake.configuration.text(length).c_str(), snake.configuration.bend(joint,1).text(length).c_str());
            fold(tSnakeDoll(snake, joint, 1, grid), joint + 1);

            // left bend, except for the first joint
            if (snake.configuration.partition != 0)
            {
                tprintf("%s -> %s\n", snake.configuration.text(length).c_str(), snake.configuration.bend(joint, -1).text(length).c_str());
                fold(tSnakeDoll(snake, joint, -1, grid), joint + 1);
            }
        }
    }

public:
    // count of found configurations
    unsigned num_configurations;

    // constructor does all the work :)
    cSnakeFolder(int n) : length(n), grid(n), num_configurations(0)
    {
        num_joints = length - 1;

        // launch recursive folding
        fold(tSnakeDoll(length), 0);
    }
};

// ============================================================================
// here we go
// ============================================================================
int main(int argc, char * argv[])
{
#ifdef NDEBUG
    if (argc != 2) panic("give me a snake length or else");
    int length = atoi(argv[1]);
#else
    (void)argc; (void)argv;
    int length = 12;
#endif // NDEBUG

    if (length <= 0 || length >= MAX_LENGTH) panic("a snake of that length is hardly foldable");

    time_t start = time(NULL);
    cSnakeFolder snakes(length);
    time_t duration = time(NULL) - start;

    printf ("Found %d configuration%c of length %d in %lds\n", snakes.num_configurations, (snakes.num_configurations == 1) ? '\0' : 's', length, duration);
    return 0;
}

with or without hash tables, Microsoft compiler performs miserably: g++ build is 3 times faster.

The algorithm uses practically no memory.

Since collision check is roughly in O(n), computation time should be in O(nkⁿ), with k slightly lower than 3.
On my [email protected], n=17 takes about 1:30 (about 2 million snakes/minute).

I am not done optimizing, but I would not expect more than a x3 gain, so basically I can hope to reach maybe n=20 under an hour, or n=24 under a day.

Reaching the first known unbendable shape (n=31) would take between a few years and a decade, assuming no power outages.

Eliminating duplicates without any storage

My initial approach was to store all configurations in a huge hash table, to eliminate duplicates by checking the presence of a previously computed symetric configuration.

Thanks to the aforementioned article, it became clear that, since partitions and foldings are stored as bitfields, they can be compared like any numerical value.
So to eliminate one member of a symetric pair, you can simply compare both elements and systematically keep the smallest one (or the greatest one, as you like).

Thus, testing a configuration for duplication amounts to computing the symetric partition, and if both are identical, the folding. No memory is needed at all.

Clearly collision check will be the most time-consuming part, so reducing these computations is a major time-saver.

I use a recursive scan of the sake from the tail down, adding a single joint at each level. Thus a new ragdoll instance is built on top of the parent configuration, with a single aditional bend.

This means bends are applied in a sequential order, which seems to be enough to avoid self-collisions in almost all cases.

When self-collision is detected, the bends that lead to the offending move are applied in all possible orders until legit folding is found or all combinations are exhausted.

Static check

Before even thinking about moving parts, I found it more efficient to test the static final shape of a snake for self-intersections.

This is done by drawing the snake on a grid. Each possible point is plotted from the head down. If there is a self-intersection, at least a pair of points will fall on the same location. This requires exactly N plots for any snake configuration, for a constant O(N) time.

The main advantage of this approach is that the static test alone will simply select valid self-avoiding paths on a square lattice, which allows to test the whole algorithm by inhibiting dynamic collision detection and making sure we find the correct count of such paths.

Dynamic check

When a snake folds around one joint, each rotated segment will sweep an area whose shape is anything but trivial.
Clearly you can check collisions by testing inclusion within all such swept areas individually. A global check would be more efficient, but given the areas complexity I can't think of any (except maybe using a GPU to draw all areas and perform a global hit check).

Since the static test takes care of the starting and ending positions of each segment, we just need to check intersections with the arcs swept by each rotating segments.

After an interesting discussion with trichoplax and a bit of JavaScript to get my bearings, I came up with this method:

For any segment that does not contain I, the swept area is bound by 2 arcs (and 2 segments already taken care of by the static check).

If I falls within the segment, the arc swept by I must also be taken into account.

This means we can check each unmoving segment against each rotating segment with 2 or 3 segment-with-arc intersections

I used vector geometry to avoid trigonometric functions altogether.
Vector operations produce compact and (relatively) readable code.

Does it work?

Inhibiting dynamic collision detection produces the correct self-avoiding paths count up to n=19, so I'm pretty confident the global layout works.

Dynamic collision detection produces consistent results, though the check of bends in different order is missing (for now).
As a consequence, the program counts snakes that can be bent from the head down (i.e. with joints folded in order of increasing distance from the head).

Compile with g++ -O3 -std=c++11
I use MinGW under Win7 with g++4.8 for "linux" builds, so portabiity is not 100% guaranteed.

g++ build is 2.5 times faster and uses about 15% less memory.
On the other hand, I find Microsoft IDE a lot more convenient than Qt Creator, with better on-the-fly error checking and more powerful debug features.

The g++ build uses about 32 bytes/configuration.
A snake configuration is 8 bytes, the remaining 24 are mostly used by hash table internals.
For convenience I transfer the results from the hash table into a sorted vector in the end, which increases memory consumption significantly.
Even so, the program can compute up to n=19 before hitting the Win32 2Gb allocation limit, and compiling for x64 would allow to go a couple of steps beyond.

A N size snake has N-1 joints.
Each joint can be left straight or bent to the left or right (3 possibilities).
The number of possible foldings is thus 3^N-1

For any segment that does not contain I, the swept area is bound by 2 arcs and 2 segments.
If I falls within the segment, a third arc must be taken into account.

Compile with g++ -O3 -std=c++11
I use MinGW under Win7 with g++4.8 for "linux" builds, so portability is not 100% guaranteed.

g++ build is 2.5 times faster and uses about 15% less memory.

The g++ build uses about 32 bytes/configuration.
A snake configuration is 8 bytes, the remaining 24 are mostly used by hash table internals.
For convenience I transfer the results from the hash table into a sorted vector in the end, which increases memory consumption significantly (from 32 to 40 bytes/configuration).
Even so, the program can compute up to n=19 before hitting the Win32 2Gb allocation limit, and compiling for x64 would allow to go a couple of steps beyond.

A N size snake has N-1 joints.
Each joint can be left straight or bent to the left or right (3 possibilities).
The number of possible foldings is thus 3^N-1.
Collisions will reduce that number somewhat, so the actual number is close to 2.7^N-1

For any segment that does not contain I, the swept area is bound by 2 arcs and 2 segments.
If I falls within the segment, a third arc must be taken into account.