Liczba osiągalnych orientacji węża

11

To wyzwanie nie dotyczy gry Snake.

Wyobraź sobie węża 2d utworzonego przez narysowanie poziomej linii długości n. W punktach całkowitych wzdłuż ciała, wąż ten może obracać ciało o 90 stopni. Jeśli na początku zdefiniujemy przód węża, który będzie po lewej stronie, obrót spowoduje przesunięcie tylnej części węża, a przednia część pozostanie na swoim miejscu. Powtarzając obroty, można uzyskać wiele różnych kształtów ciała węża.

Zasady

  1. Jedna część ciała węża nie może pokrywać się z inną.
  2. Musi być możliwe osiągnięcie ostatecznej orientacji bez nakładania się na siebie części ciała węża. Dwa dotykające się punkty są liczone jako pokrywające się w tym problemie.
  3. Uważam, że wąż i jego odwrotność mają ten sam kształt.

Zadanie

Jeśli chodzi o rotację, translację i symetrię lustrzaną, jaka jest łączna liczba różnych kształtów ciała węża, które można wykonać?

Przykład obrotu części ciała węży. Wyobraź sobie, n=10że wąż jest w orientacji początkowej linii prostej. Teraz obróć w punkcie o 490 stopni w lewo. Dostajemy węża od 4do 10(ogona węża) leżącego pionowo, a węża od 0do 4leżącego poziomo. Wąż ma teraz jeden kąt prosty w ciele.

Oto kilka przykładów dzięki Martinowi Büttnerowi.

Zaczynamy od poziomego węża.

wprowadź opis zdjęcia tutaj

Teraz obracamy od pozycji 4.

wprowadź opis zdjęcia tutaj

W tej orientacji kończymy po rotacji.

wprowadź opis zdjęcia tutaj

Rozważmy teraz orientację innego węża.

wprowadź opis zdjęcia tutaj

Teraz możemy zobaczyć nielegalny ruch, w którym podczas rotacji wystąpiłoby nakładanie się.

przykład kolizji

Wynik

Twój wynik jest najwyższy, ndla którego Twój kod może rozwiązać problem w mniej niż minutę na moim komputerze.

Kiedy nastąpi obrót, przesunie on wraz z nim o połowę węża. Musimy się martwić, czy którakolwiek z obracanych części może zachodzić na część węża podczas obrotu. Dla uproszczenia możemy założyć, że wąż ma szerokość zero. Można obracać tylko w określonym punkcie węża o maksymalnie 90 stopni w kierunku zgodnym z ruchem wskazówek zegara. Ponieważ nigdy nie można całkowicie złożyć węża na pół, ponieważ wymagałoby to dwóch obrotów w tym samym punkcie w tym samym kierunku.

Kształty, których nie można wykonać

Prostym przykładem kształtu, którego nie można wykonać, jest stolica T. Bardziej wyrafinowana wersja to

wprowadź opis zdjęcia tutaj

(Dziękuję Haraldowi Hanche-Olsenowi za ten przykład)

W tym przykładzie wszystkie sąsiednie linie poziome są od siebie oddalone o 1, podobnie jak pionowe. W związku z tym nie ma legalnego przesunięcia z tej pozycji, a ponieważ problem jest odwracalny, nie ma zatem możliwości, aby dostać się z pozycji początkowej.

Języki i biblioteki

Możesz używać dowolnego języka, który ma swobodnie dostępny kompilator / tłumacz / etc. dla systemu Linux i dowolnych bibliotek, które są również bezpłatnie dostępne dla systemu Linux.

Moja maszyna Czasy zostaną uruchomione na moim komputerze. Jest to standardowa instalacja ubuntu na ośmiordzeniowym procesorze AMD FX-8350. Oznacza to również, że muszę być w stanie uruchomić Twój kod. W związku z tym korzystaj tylko z łatwo dostępnego bezpłatnego oprogramowania i dołącz pełne instrukcje, jak skompilować i uruchomić kod.


źródło
1
@TApicella Dzięki za pytania. Kiedy mówię „Kiedy nastąpi obrót, poruszy się nim o połowę węża”, nie miałem na myśli 50 procent. Miałem tylko na myśli część przed punktem obrotu i część po nim. Jeśli obracasz się od 0 wzdłuż węża, obracasz całą rzecz!
2
@TApicella O twoim drugim pytaniu. Chodzi o to, że nie ma legalnej rotacji w stosunku do stanowiska, które podałem. Jeśli był osiągalny, musi istnieć możliwość powrotu do poziomej orientacji początkowej poprzez sekwencję obrotów (odwrotność tych, które podjąłbyś, aby dostać się do orientacji końcowej). Czy możesz opisać jeden legalny obrót, który Twoim zdaniem możesz wykonać z tej pozycji? Żeby było jasne, wąż nie rośnie. Zawsze ma tę samą długość przez cały czas.
3
@TApicella Wygląda na to, że spodziewasz się wzrostu węża. Jednak jego rozmiar jest ustalony. Zaczynasz z jednym długim wężem i wszystko, co możesz zrobić, to złożyć jego części o 90 stopni. Z obecnej pozycji nie można w ogóle zastosować żadnego zagięcia, które prowadziłoby do poprzedniego etapu węża.
Martin Ender
1
Czy możesz spasować w punkcie więcej niż jeden raz (tam iz powrotem)? Jeśli możesz, to sprawia, że ​​jest to dość skomplikowane.
randomra
1
@ randomra Rzeczywiście możesz tak długo, jak nigdy nie pójdziesz dalej niż dziewięćdziesiąt stopni.

Odpowiedzi:

5

Python 3 - wynik tymczasowy: n = 11 (n = 13 z PyPy *)

Ponieważ w pierwszym tygodniu nie było odpowiedzi, oto przykład w Pythonie zachęcający do rywalizacji. Starałem się, aby był on wystarczająco czytelny, aby łatwo można było zidentyfikować nieefektywności, dając pomysły na inne odpowiedzi.

Podejście

  • Zacznij od prostego węża i znajdź wszystkie pozycje, które można legalnie osiągnąć jednym ruchem.
  • Znajdź wszystkie pozycje, do których można legalnie dotrzeć z tych pozycji, które nie zostały jeszcze zidentyfikowane.
  • Powtarzaj, aż nie będzie już więcej, i zwróć liczbę znalezionych pozycji.

Kod

(teraz z kilkoma doktrynami i stwierdzeniami po mojej niepoprawnej pierwszej próbie)

'''
Snake combinations

A snake is represented by a tuple giving the relative orientation at each joint.
A length n snake has n-1 joints.
Each relative orientation is one of the following:

0: Clockwise 90 degrees
1: Straight
2: Anticlockwise 90 degrees

So a straight snake of length 4 has 3 joints all set to 1:

(1, 1, 1)

x increases to the right
y increases upwards

'''


import turtle


def all_coords(state):
    '''Return list of coords starting from (0,0) heading right.'''
    current = (1, 0)
    heading = 0
    coords = [(0,0), (1,0)]
    for item in state:
        heading += item + 3
        heading %= 4
        offset = ((1,0), (0,1), (-1,0), (0,-1))[heading]
        current = tuple(current[i]+offset[i] for i in (0,1))
        coords.append(current)
    return coords


def line_segments(coords, pivot):
    '''Return list of line segments joining consecutive coords up to pivot-1.'''
    return [(coords[i], coords[i+1]) for i in range(pivot+1)]


def rotation_direction(coords, pivot, coords_in_final_after_pivot):
    '''Return -1 if turning clockwise, 1 if turning anticlockwise.'''
    pivot_coord = coords[pivot + 1]
    initial_coord = coords[pivot + 2]
    final_coord = coords_in_final_after_pivot[0]
    initial_direction = tuple(initial_coord[i] - pivot_coord[i] for i in (0,1))
    final_direction = tuple(final_coord[i] - pivot_coord[i] for i in (0,1))
    return (initial_direction[0] * final_direction[1] -
            initial_direction[1] * final_direction[0]
            )


def intersects(arc, line):
    '''Return True if the arc intersects the line segment.'''
    if line_segment_cuts_circle(arc, line):
        cut_points = points_cutting_circle(arc, line)
        if cut_points and cut_point_is_on_arc(arc, cut_points):
            return True


def line_segment_cuts_circle(arc, line):
    '''Return True if the line endpoints are not both inside or outside.'''
    centre, point, direction = arc
    start, finish = line
    point_distance_squared = distance_squared(centre, point)
    start_distance_squared = distance_squared(centre, start)
    finish_distance_squared = distance_squared(centre, finish)
    start_sign = start_distance_squared - point_distance_squared
    finish_sign = finish_distance_squared - point_distance_squared
    if start_sign * finish_sign <= 0:
        return True


def distance_squared(centre, point):
    '''Return the square of the distance between centre and point.'''
    return sum((point[i] - centre[i]) ** 2 for i in (0,1))


def cut_point_is_on_arc(arc, cut_points):
    '''Return True if any intersection point with circle is on arc.'''
    centre, arc_start, direction = arc
    relative_start = tuple(arc_start[i] - centre[i] for i in (0,1))
    relative_midpoint = ((relative_start[0] - direction*relative_start[1])/2,
                         (relative_start[1] + direction*relative_start[0])/2
                         )
    span_squared = distance_squared(relative_start, relative_midpoint)
    for cut_point in cut_points:
        relative_cut_point = tuple(cut_point[i] - centre[i] for i in (0,1))
        spacing_squared = distance_squared(relative_cut_point,
                                           relative_midpoint
                                           )
        if spacing_squared <= span_squared:
            return True


def points_cutting_circle(arc, line):
    '''Return list of points where line segment cuts circle.'''
    points = []
    start, finish = line
    centre, arc_start, direction = arc
    radius_squared = distance_squared(centre, arc_start)
    length_squared = distance_squared(start, finish)
    relative_start = tuple(start[i] - centre[i] for i in (0,1))
    relative_finish = tuple(finish[i] - centre[i] for i in (0,1))
    relative_midpoint = tuple((relative_start[i] +
                               relative_finish[i]
                               )*0.5 for i in (0,1))
    determinant = (relative_start[0]*relative_finish[1] -
                   relative_finish[0]*relative_start[1]
                   )
    determinant_squared = determinant ** 2
    discriminant = radius_squared * length_squared - determinant_squared
    offset = tuple(finish[i] - start[i] for i in (0,1))
    sgn = (1, -1)[offset[1] < 0]
    root_discriminant = discriminant ** 0.5
    one_over_length_squared = 1 / length_squared
    for sign in (-1, 1):
        x = (determinant * offset[1] +
             sign * sgn * offset[0] * root_discriminant
             ) * one_over_length_squared
        y = (-determinant * offset[0] +
             sign * abs(offset[1]) * root_discriminant
             ) * one_over_length_squared
        check = distance_squared(relative_midpoint, (x,y))
        if check <= length_squared * 0.25:
            points.append((centre[0] + x, centre[1] + y))
    return points


def potential_neighbours(candidate):
    '''Return list of states one turn away from candidate.'''
    states = []
    for i in range(len(candidate)):
        for orientation in range(3):
            if abs(candidate[i] - orientation) == 1:
                state = list(candidate)
                state[i] = orientation
                states.append(tuple(state))
    return states


def reachable(initial, final):
    '''
    Return True if final state can be reached in one legal move.

    >>> reachable((1,0,0), (0,0,0))
    False

    >>> reachable((0,1,0), (0,0,0))
    False

    >>> reachable((0,0,1), (0,0,0))
    False

    >>> reachable((1,2,2), (2,2,2))
    False

    >>> reachable((2,1,2), (2,2,2))
    False

    >>> reachable((2,2,1), (2,2,2))
    False

    >>> reachable((1,2,1,2,1,1,2,2,1), (1,2,1,2,1,1,2,1,1))
    False

    '''
    pivot = -1
    for i in range(len(initial)):
        if initial[i] != final[i]:
            pivot = i
            break

    assert pivot > -1, '''
        No pivot between {} and {}'''.format(initial, final)
    assert initial[pivot + 1:] == final[pivot + 1:], '''
        More than one pivot between {} and {}'''.format(initial, final)

    coords_in_initial = all_coords(initial)
    coords_in_final_after_pivot = all_coords(final)[pivot+2:]
    coords_in_initial_after_pivot = coords_in_initial[pivot+2:]
    line_segments_up_to_pivot = line_segments(coords_in_initial, pivot)

    direction = rotation_direction(coords_in_initial,
                                   pivot,
                                   coords_in_final_after_pivot
                                   )

    pivot_point = coords_in_initial[pivot + 1]

    for point in coords_in_initial_after_pivot:
        arc = (pivot_point, point, direction)
        if any(intersects(arc, line) for line in line_segments_up_to_pivot):
            return False
    return True


def display(snake):
    '''Display a line diagram of the snake.

    Accepts a snake as either a tuple:

    (1, 1, 2, 0)

    or a string:

    "1120"

    '''
    snake = tuple(int(s) for s in snake)
    coords = all_coords(snake)

    turtle.clearscreen()
    t = turtle.Turtle()
    t.hideturtle()
    s = t.screen
    s.tracer(0)

    width, height = s.window_width(), s.window_height()

    x_min = min(coord[0] for coord in coords)
    x_max = max(coord[0] for coord in coords)
    y_min = min(coord[1] for coord in coords)
    y_max = max(coord[1] for coord in coords)
    unit_length = min(width // (x_max - x_min + 1),
                      height // (y_max - y_min + 1)
                      )

    origin_x = -(x_min + x_max) * unit_length // 2
    origin_y = -(y_min + y_max) * unit_length // 2

    pen_width = max(1, unit_length // 20)
    t.pensize(pen_width)
    dot_size = max(4, pen_width * 3)

    t.penup()
    t.setpos(origin_x, origin_y)
    t.pendown()

    t.forward(unit_length)
    for joint in snake:
        t.dot(dot_size)
        t.left((joint - 1) * 90)
        t.forward(unit_length)
    s.update()


def neighbours(origin, excluded=()):
    '''Return list of states reachable in one step.'''
    states = []
    for candidate in potential_neighbours(origin):
        if candidate not in excluded and reachable(origin, candidate):
            states.append(candidate)
    return states


def mirrored_or_backwards(candidates):
    '''Return set of states that are equivalent to a state in candidates.'''
    states = set()
    for candidate in candidates:
        mirrored = tuple(2 - joint for joint in candidate)
        backwards = candidate[::-1]
        mirrored_backwards = mirrored[::-1]
        states |= set((mirrored, backwards, mirrored_backwards))
    return states


def possible_snakes(snake):
    '''
    Return the set of possible arrangements of the given snake.

    >>> possible_snakes((1, 1, 1, 1, 2, 1, 2, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 2, 2, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1))
    {(1, 1, 1, 1, 2, 1, 2, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 2, 2, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1)}

    '''
    reached = set()
    candidates = set((snake,))

    while candidates:
        candidate = candidates.pop()
        reached.add(candidate)
        new_candidates = neighbours(candidate, reached)
        reached |= mirrored_or_backwards(new_candidates)
        set_of_new_candidates = set(new_candidates)
        reached |= set_of_new_candidates
        candidates |= set_of_new_candidates

    excluded = set()
    final_answers = set()
    while reached:
        candidate = reached.pop()
        if candidate not in excluded:
            final_answers.add(candidate)
            excluded |= mirrored_or_backwards([candidate])

    return final_answers


def straight_derived_snakes(length):
    '''Return the set of possible arrangements of a snake of this length.'''
    straight_line = (1,) * max(length-1, 0)
    return possible_snakes(straight_line)


if __name__ == '__main__':
    import doctest
    doctest.testmod()
    import sys
    arguments = sys.argv[1:]
    if arguments:
        length = int(arguments[0])
    else:
        length = int(input('Enter the length of the snake:'))
    print(len(straight_derived_snakes(length)))

Wyniki

Na mojej maszynie najdłuższy wąż, który można obliczyć w czasie krótszym niż 1 minuta, ma długość 11 (lub długość 13 z PyPy *). Jest to oczywiście tylko wynik tymczasowy, dopóki nie dowiemy się, jaki jest oficjalny wynik z maszyny Lembika.

Dla porównania oto wyniki dla pierwszych kilku wartości n:

 n | reachable orientations
-----------------------------
 0 | 1
 1 | 1
 2 | 2
 3 | 4
 4 | 9
 5 | 22
 6 | 56
 7 | 147
 8 | 388
 9 | 1047
10 | 2806
11 | 7600
12 | 20437
13 | 55313
14 | 148752
15 | 401629
16 | 1078746
17 | MemoryError (on my machine)

Daj mi znać, jeśli którekolwiek z nich okażą się nieprawidłowe.

Jeśli masz przykład aranżacji, która nie powinna być możliwa do rozwinięcia, możesz użyć funkcji, neighbours(snake)aby znaleźć aranżacje osiągalne w jednym kroku, jako test kodu. snakejest krotką kierunków połączeń (0 dla wskazówek zegara, 1 dla linii prostych, 2 dla wskazówek zegara). Na przykład (1,1,1) jest prostym wężem o długości 4 (z 3 stawami).

Wyobrażanie sobie

Aby zwizualizować węża, o którym myślisz, lub któregoś z węży powracających neighbours, możesz użyć tej funkcji display(snake). To akceptuje krotkę jak inne funkcje, ale ponieważ nie jest używana przez program główny i dlatego nie musi być szybka, akceptuje również ciąg, dla twojej wygody.

display((1,1,2,0)) jest równa display("1120")

Jak wspomina Lembik w komentarzach, moje wyniki są identyczne jak OEIS A037245, który nie uwzględnia skrzyżowań podczas rotacji. Wynika to z faktu, że dla pierwszych kilku wartości n nie ma różnicy - wszystkie kształty, które nie przecinają się, można osiągnąć, składając prostego węża. Poprawność kodu można sprawdzić, wywołując neighbours()węża, którego nie można rozwinąć bez skrzyżowania. Ponieważ nie ma sąsiadów, przyczyni się tylko do sekwencji OEIS, a nie do tej. Najmniejszy przykład, jaki znam, to ten wąż o długości 31, o którym wspomniał Lembik, dzięki Davidowi K :

(1,1,1,1,2,1,2,1,1,1,1,1,1,2,1,1,1,2,1,1,2,2,1,0,1,0,1,1,1,1)

display('111121211111121112112210101111') daje następujący wynik:

Obraz najkrótszego węża bez sąsiadów

Wskazówka: jeśli zmienisz rozmiar okna wyświetlacza, a następnie ponownie wywołasz ekran, wąż zostanie dopasowany do nowego rozmiaru okna.

Chciałbym usłyszeć od każdego, kto ma krótszy przykład bez sąsiadów. Podejrzewam, że najkrótszy taki przykład wyznaczy najmniejsze n, dla których dwie sekwencje różnią się.


* Zauważ, że każdy przyrost n zajmuje około 3 razy tyle, więc zwiększenie z n = 11 do n = 13 wymaga prawie 10 razy więcej czasu. Właśnie dlatego PyPy pozwala tylko zwiększać n o 2, nawet jeśli działa znacznie szybciej niż standardowy interpreter Pythona.

trichopaks
źródło
6
Jeśli ten komentarz otrzyma 5 głosów pozytywnych, przyjrzę się dodaniu opcji włączenia wizualizacji możliwych ustaleń, na wypadek, gdyby pomogło to w analizie.
trichoplax
To okazuje się niepoprawne : P
Geobity
@Geobits Myślę, że tym razem mam to dobrze ...
trichoplax
1
@Jakube Można to otworzyć na wiele sposobów, np. W kolejności # 1 # 3 # 2 # 4 # 5 # 6.
randomra
1

C ++ 11 - prawie działa :)

Po przeczytaniu tego artykułu zebrałem odrobinę mądrości od tego faceta, który najwyraźniej pracował przez 25 lat nad mniej skomplikowanym problemem liczenia ścieżek unikania na kwadratowej sieci.

#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;
}

Budowanie pliku wykonywalnego

Kompiluj z Używam MinGW pod Win7 z g ++ 4.8 dla kompilacji „linux”, więc przenośność nie jest gwarantowana w 100%.g++ -O3 -std=c++11

Działa również (w pewnym sensie) ze standardowym projektem MSVC2013

Niezdefiniowanie NDEBUGpozwala uzyskać ślady wykonania algorytmu i podsumowanie znalezionych konfiguracji.

Występy

z tablicami skrótów lub bez nich kompilator Microsoft działa marnie: kompilacja g ++ jest 3 razy szybsza .

Algorytm praktycznie nie wykorzystuje pamięci.

Ponieważ kontrola kolizji jest w przybliżeniu w O (n), czas obliczeń powinien być w O (nk n ), przy k nieco mniejszym niż 3.
W moim i3-2100@3,1 GHz, n = 17 zajmuje około 1:30 (około 2 milionów węże / minutę).

Nie skończyłem optymalizacji, ale nie spodziewałbym się więcej niż x3, więc w zasadzie mogę mieć nadzieję, że osiągnę n = 20 poniżej godziny lub n = 24 poniżej dnia.

Osiągnięcie pierwszego znanego nieugiętego kształtu (n = 31) zajęłoby od kilku lat do dekady, przy założeniu braku przerw w dostawie prądu.

Liczenie kształtów

N wielkość wąż N1 stawów.
Każde złącze może być lewe lub wygięte w lewo lub w prawo (3 możliwości).
Liczba możliwych składania wynosi zatem 3 N-1 .
Zderzenia nieco zmniejszą tę liczbę, więc faktyczna liczba jest bliska 2,7 N-1

Jednak wiele takich fałd prowadzi do identycznych kształtów.

dwa kształty są identyczne, jeśli istnieje obrót lub symetria, które mogą przekształcić jeden w drugi.

Zdefiniujmy segment jako dowolną prostą część złożonego korpusu.
Na przykład wąż wielkości 5 złożony na 2. stawie miałby 2 segmenty (jeden 2 jednostki długości i drugi 3 jednostki długości).
Pierwszy segment zostanie nazwany głową , a ostatni ogon .

Konwencjonalnie orientujemy głowę węży poziomo z ciałem skierowanym w prawo (jak na pierwszej figurze OP).

Daną liczbę oznaczamy za pomocą listy podpisanych długości segmentów, przy czym długości dodatnie wskazują na prawe fałdowanie, a ujemne na lewe.
Początkowa długość jest zgodna z konwencją.

Rozdzielanie segmentów i zagięć

Jeśli weźmiemy pod uwagę tylko różne sposoby podziału węża o długości N na segmenty, otrzymamy podział identyczny do składu N.

Przy użyciu tego samego algorytmu, jak przedstawiono na stronie wiki jest łatwy do wygenerowania wszystkich 2 n-1 możliwe przegrody węża.

Każda przegroda z kolei wygeneruje wszystkie możliwe zagięcia, stosując lewe lub prawe zagięcie do wszystkich swoich połączeń. Jedno takie składanie będzie nazywane konfiguracją .

Wszystkie możliwe partycje mogą być reprezentowane przez liczbę całkowitą N-1 bitów, przy czym każdy bit reprezentuje obecność połączenia. Nazwamy tę liczbę całkowitą generatorem .

Przycinanie partycji

Zauważając, że wygięcie danej przegrody od głowy w dół jest równoważne zgięciu symetrycznej przegrody od ogona w górę, możemy znaleźć wszystkie pary symetrycznych przegród i wyeliminować jedną z dwóch.
Generator symetrycznej partycji jest generatorem partycji zapisanym w odwrotnej kolejności bitów, który jest banalnie łatwy i tani do wykrycia.

To wyeliminuje prawie połowę możliwych partycji, wyjątkami są partycje z generatorami „palindromowymi”, które pozostają niezmienione przez odwrócenie bitów (na przykład 00100100).

Dbanie o poziome symetrie

Dzięki naszym konwencjom (wąż zaczyna wskazywać w prawo), pierwsze zagięcie zastosowane po prawej stronie stworzy rodzinę zagięć, które będą poziomymi symetriami od tych, które różnią się tylko pierwszym zgięciem.

Jeśli zdecydujemy, że pierwsze zakręt zawsze będzie po prawej stronie, eliminujemy wszystkie poziome symetrie za jednym razem.

Mycie palindromów

Te dwa cięcia są wydajne, ale niewystarczające, aby zająć się tymi nieznośnymi palindromami.
Najbardziej szczegółowa kontrola w ogólnym przypadku wygląda następująco:

rozważ konfigurację C z partycją palindromową.

  • jeśli odwrócimy każde zagięcie w C, otrzymamy poziomą symetrię C.
  • jeśli odwrócimy C (stosując zakręty od ogona do góry), otrzymamy tę samą figurę obróconą w prawo
  • jeśli zarówno cofniemy, jak i odwrócimy C, otrzymamy tę samą liczbę obróconą w lewo.

Możemy sprawdzić każdą nową konfigurację z 3 innymi. Ponieważ jednak generujemy już tylko konfiguracje zaczynające się od skrętu w prawo, mamy tylko jedną możliwą symetrię do sprawdzenia:

  • odwrócone C rozpocznie się od skrętu w lewo, którego z założenia nie można powielić
  • z konfiguracji odwróconej i odwróconej-odwróconej tylko jedna rozpocznie się od skrętu w prawo.
    To jedyna konfiguracja, którą możemy powielić.

Eliminowanie duplikatów bez przechowywania

Moje początkowe podejście polegało na przechowywaniu wszystkich konfiguracji w ogromnej tabeli mieszającej, aby wyeliminować duplikaty, sprawdzając obecność wcześniej obliczonej konfiguracji symetrycznej.

Dzięki wspomnianemu artykułowi stało się jasne, że ponieważ partycje i foldery są przechowywane jako pola bitowe, można je porównać jak każdą wartość liczbową.
Aby więc wyeliminować jednego członka pary symetrycznej, możesz po prostu porównać oba elementy i systematycznie zachować najmniejszy (lub największy, jak chcesz).

Zatem testowanie konfiguracji duplikacji sprowadza się do obliczenia symetrycznej partycji, a jeśli obie są identyczne, zwijanie. W ogóle nie jest potrzebna pamięć.

Kolejność generacji

Oczywiste jest, że kontrola kolizji będzie najbardziej czasochłonną częścią, więc zmniejszenie tych obliczeń jest znaczną oszczędnością czasu.

Możliwym rozwiązaniem jest posiadanie „węża ragdoll”, który rozpocznie się w płaskiej konfiguracji i będzie stopniowo wyginany, aby uniknąć ponownego obliczenia całej geometrii węża dla każdej możliwej konfiguracji.

Wybierając kolejność testowania konfiguracji, aby maksymalnie przechowywać ragdoll dla każdej łącznej liczby połączeń, możemy ograniczyć liczbę instancji do N-1.

Używam rekurencyjnego skanu sake od ogona w dół, dodając pojedynczy staw na każdym poziomie. Zatem nowa instancja ragdoll jest zbudowana na konfiguracji nadrzędnej, z jednym dodatkowym zgięciem.

Oznacza to, że zakręty są stosowane w kolejności sekwencyjnej, co wydaje się wystarczające, aby uniknąć kolizji we wszystkich przypadkach.

Po wykryciu kolizji zgięcia, które prowadzą do wykroczenia, są stosowane we wszystkich możliwych rozkazach, aż do znalezienia prawidłowego złożenia lub wyczerpania wszystkich kombinacji.

Kontrola statyczna

Zanim nawet pomyślałem o ruchomych częściach, bardziej efektywne było przetestowanie statycznego ostatecznego kształtu węża pod kątem skrzyżowań.

Odbywa się to poprzez narysowanie węża na siatce. Każdy możliwy punkt jest wykreślany od głowy w dół. Jeśli istnieje skrzyżowanie, co najmniej para punktów spadnie na to samo miejsce. Wymaga to dokładnie N wykresów dla dowolnej konfiguracji węża, dla stałego czasu O (N).

Główną zaletą tego podejścia jest to, że sam test statyczny po prostu wybierze prawidłowe sam unikające się ścieżki na kwadratowej sieci, co pozwala przetestować cały algorytm przez hamowanie dynamicznego wykrywania kolizji i upewnienie się, że znajdziemy prawidłową liczbę takich ścieżek.

Kontrola dynamiczna

Kiedy wąż składa się wokół jednego stawu, każdy obrócony segment zamiata obszar, którego kształt nie jest trywialny.
Oczywiście można sprawdzić kolizje, testując indywidualnie włączenie we wszystkich takich obszarach. Globalna kontrola byłaby bardziej wydajna, ale biorąc pod uwagę złożoność obszarów, o której nie mogę myśleć (oprócz użycia GPU do narysowania wszystkich obszarów i przeprowadzenia globalnej kontroli trafień).

Ponieważ test statyczny zajmuje się początkową i końcową pozycją każdego segmentu, musimy tylko sprawdzić przecięcia z łukami przeciągniętymi przez każdy obracający się segment.

Po interesującej dyskusji z trichoplaxem i odrobiną JavaScript, aby uzyskać orientację , wymyśliłem tę metodę:

Spróbuj zadzwonić w kilku słowach, jeśli zadzwonisz

  • C środek obrotu,
  • S obrotowy segment o dowolnej długości i kierunku, który nie zawiera C ,
  • L linia przedłużająca S
  • H linia prostopadła do L przechodząca przez C ,
  • I przecięcie L i H ,

matematyka
(źródło: free.fr )

Dla każdego segmentu, który nie zawiera I , obszar przemiatania jest ograniczony przez 2 łuki (i 2 segmenty, które zostały już zajęte przez kontrolę statyczną).

Jeśli ja wchodzi w segmencie łuk zmieciony przez Muszę być również brane pod uwagę.

Oznacza to, że możemy sprawdzić każdy nieruchliwy segment w stosunku do każdego segmentu obrotowego za pomocą 2 lub 3 odcinków z łukiem

Użyłem geometrii wektorowej, aby całkowicie uniknąć funkcji trygonometrycznych.
Operacje wektorowe dają zwarty i (względnie) czytelny kod.

Przecięcie segmentu do łuku wymaga wektora zmiennoprzecinkowego, ale logika powinna być odporna na błędy zaokrąglania.
Znalazłem to eleganckie i wydajne rozwiązanie w niejasnym wpisie na forum. Zastanawiam się, dlaczego nie jest szerzej rozpowszechniany.

Czy to działa?

Hamowanie dynamicznego wykrywania kolizji powoduje powstanie prawidłowej liczby unikających ścieżek do n = 19, więc jestem pewien, że globalny układ działa.

Dynamiczne wykrywanie kolizji daje spójne wyniki, chociaż brakuje kontroli zgięć w innej kolejności (na razie).
W rezultacie program zlicza węże, które można zgiąć od głowy w dół (tj. Ze stawami złożonymi w kolejności rosnącej od głowy).

Glorfindel
źródło