libmath-planepath-perl 117-1 / usr / share / perl5 / Math / PlanePath / Flowsnake.pm

# Copyright 2010, 2011, 2012, 2013, 2014 Kevin Ryde

# This file is part of Math-PlanePath.
#
# Math-PlanePath is free software; you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the Free
# Software Foundation; either version 3, or (at your option) any later
# version.
#
# Math-PlanePath is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
# or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
# for more details.
#
# You should have received a copy of the GNU General Public License along
# with Math-PlanePath.  If not, see <http://www.gnu.org/licenses/>.


# math-image --path=Flowsnake --lines --scale=10
# math-image --path=Flowsnake --all --output=numbers_dash
# math-image --path=Flowsnake,arms=3 --all --output=numbers_dash
#
# Martin Gardner, "In which `Monster' Curves Force Redefinition of the Word
# `Curve'", Scientific American 235, December 1976, pages 124-133.
#
# http://80386.nl/pub/gosper-level21.png
#
# http://www.mathcurve.com/fractals/gosper/gosper.shtml
#
# plain hexagonal tiling http://tilingsearch.org/HTML/data136/F666.html

# Jeffrey Ventrella
# root-7 family
# "inner-flip" which is initial state reversal



package Math::PlanePath::Flowsnake;
use 5.004;
use strict;

use vars '$VERSION', '@ISA';
$VERSION = 117;

# inherit: new(), rect_to_n_range(), arms_count(), n_start(),
#          parameter_info_array(), xy_is_visited()
use Math::PlanePath::FlowsnakeCentres 55; # v.55 inheritance switch-around
@ISA = ('Math::PlanePath::FlowsnakeCentres');
use Math::PlanePath;
*_divrem_mutate = \&Math::PlanePath::_divrem_mutate;

use Math::PlanePath::Base::Generic
  'is_infinite',
  'round_nearest';
use Math::PlanePath::Base::Digits
  'digit_split_lowtohigh',
  'round_down_pow';

# uncomment this to run the ### lines
# use Smart::Comments;


# (i,j)*(2+w) = (2i-j,2j+i+j) = (2i-j,3j+i)
# (x,y)*(2+w) = 2x + (x-3y)/2, 2y + (x+y)/2
#             = (4x + x-3y)/2, (4y + x+y)/2
#             = (5x-3y)/2, (x+5y)/2

{
  my @x_negative_at_n = (undef, 23, 1, 1);
  sub x_negative_at_n {
    my ($self) = @_;
    return $x_negative_at_n[$self->{'arms'}];
  }
}
{
  my @y_negative_at_n = (undef, 8598, 7, 2);
  sub y_negative_at_n {
    my ($self) = @_;
    return $y_negative_at_n[$self->{'arms'}];
  }
}

{
  my @_UNDOCUMENTED__dxdy_list_at_n = (undef, 11, 6, 9);
  sub _UNDOCUMENTED__dxdy_list_at_n {
    my ($self) = @_;
    return $_UNDOCUMENTED__dxdy_list_at_n[$self->{'arms'}];
  }
}

# Table generated by tools/flowsnake-table.pl.
# next_state length 84
my @next_state = (0, 21,49,28, 0, 0,77, 70, 7, 7,35,42,14, 7,  # 0,7
                  14,35,63,42,14,14, 7,  0,21,21,49,56,28,21,  # 14,21
                  28,49,77,56,28,28,21, 14,35,35,63,70,42,35,  # 28,35
                  42,63, 7,70,42,42,35, 28,49,49,77, 0,56,49,  # 42,49
                  56,77,21, 0,56,56,49, 42,63,63, 7,14,70,63,  # 56,63
                  70, 7,35,14,70,70,63, 56,77,77,21,28, 0,77);  # 70,77
my @digit_to_i = (0,  1, 1, 0,-1, 0, 1,  0, 1, 2, 3, 2, 1, 1,  # 0,7
                  0,  0,-1,-1,-2,-2,-2,  0, 1, 1, 1, 0, 0,-1,  # 14,21
                  0, -1,-2,-1,-1,-2,-3,  0, 0,-1,-2,-2,-1,-2,  # 28,35
                  0, -1,-1, 0, 1, 0,-1,  0,-1,-2,-3,-2,-1,-1,  # 42,49
                  0,  0, 1, 1, 2, 2, 2,  0,-1,-1,-1, 0, 0, 1,  # 56,63
                  0,  1, 2, 1, 1, 2, 3,  0, 0, 1, 2, 2, 1,2);  # 70,77
my @digit_to_j = (0,  0, 1, 1, 2, 2, 2,  0,-1,-1,-1, 0, 0, 1,  # 0,7
                  0,  1, 2, 1, 1, 2, 3,  0, 0, 1, 2, 2, 1, 2,  # 14,21
                  0,  1, 1, 0,-1, 0, 1,  0, 1, 2, 3, 2, 1, 1,  # 28,35
                  0,  0,-1,-1,-2,-2,-2,  0, 1, 1, 1, 0, 0,-1,  # 42,49
                  0, -1,-2,-1,-1,-2,-3,  0, 0,-1,-2,-2,-1,-2,  # 56,63
                  0, -1,-1, 0, 1, 0,-1,  0,-1,-2,-3,-2,-1,-1);  # 70,77

# state 0 to 11
my @dir6_to_di = (1, 0,-1, -1, 0, 1);
my @dir6_to_dj = (0, 1, 1,  0,-1,-1);

sub n_to_xy {
  my ($self, $n) = @_;
  ### Flowsnake n_to_xy(): $n

  if ($n < 0) { return; }
  if (is_infinite($n)) { return ($n,$n); }

  my $int = int($n);
  $n -= $int;  # fraction part
  ### $int
  ### frac: $n

  my $state;
  {
    my $arm = _divrem_mutate ($int, $self->{'arms'});
    $state = 28 * $arm;  # initial rotation

    # adjust so that for arms=2 point N=1 has $int==1
    # or for arms=3 then points N=1 and N=2 have $int==1
    if ($arm) { $int += 1; }
  }
  ### initial state: $state

  my $i = my $j = $int*0;  # bignum zero

  foreach my $digit (reverse digit_split_lowtohigh($int,7)) { # high to low
    ### at: "state=$state digit=$digit  i=$i,j=$j  di=".$digit_to_i[$state+$digit]." dj=".$digit_to_j[$state+$digit]

    # (i,j) *= (2+w), being (i,j) = 2*(i,j)+rot60(i,j)
    # then add low digit pos
    #
    $state += $digit;
    ($i, $j) = (2*$i - $j + $digit_to_i[$state],
                3*$j + $i + $digit_to_j[$state]);
    $state = $next_state[$state];
  }
  ### integer: "i=$i, j=$j"

  # fraction in final $state direction
  if ($n) {
    ### apply: "frac=$n  state=$state"
    $state = int($state/14);   # divide to direction 0 to 5
    $i = $n * $dir6_to_di[$state] + $i;
    $j = $n * $dir6_to_dj[$state] + $j;
  }

  ### ret: "$i, $j  x=".(2*$i+$j)." y=$j"
  return (2*$i+$j,
          $j);

}

# Table generated by tools/flowsnake-table.pl.
my @digit_to_next_di
  = (0, -1,-1, 1, 1, 1,undef,  1, 1,-1,-1, 0, 1,undef,  # 0,7
     -1, 0,-1, 0, 0, 1,undef,  0, 0,-1, 0,-1, 0,undef,  # 14,21
     -1, 1, 0,-1,-1, 0,undef, -1,-1, 0, 1,-1,-1,undef,  # 28,35
     0,  1, 1,-1,-1,-1,undef, -1,-1, 1, 1, 0,-1,undef,  # 42,49
     1,  0, 1, 0, 0,-1,undef,  0, 0, 1, 0, 1, 0,undef,  # 56,63
     1, -1, 0, 1, 1, 0,undef,  1, 1, 0,-1, 1, 1,undef,  # 70,77
     1, -1,-1, 1, 1, 0,undef,  1, 1, 0,-1, 0, 1,undef,  # 84,91
     0, -1,-1, 0, 0, 1,undef,  1, 1,-1, 0,-1, 1,undef,  # 98,105
     -1, 0, 0,-1,-1, 1,undef,  0, 0,-1, 1,-1, 0,undef,  # 112,119
     -1, 1, 1,-1,-1, 0,undef, -1,-1, 0, 1, 0,-1,undef,  # 126,133
     0,  1, 1, 0, 0,-1,undef, -1,-1, 1, 0, 1,-1,undef,  # 140,147
     1,  0, 0, 1, 1,-1,undef,  0, 0, 1,-1, 1,0);
my @digit_to_next_dj
  = (1,  0, 1, 0, 0,-1,undef,  0, 0, 1, 0, 1, 0,undef,  # 0,7
     1, -1, 0, 1, 1, 0,undef,  1, 1, 0,-1, 1, 1,undef,  # 14,21
     0, -1,-1, 1, 1, 1,undef,  1, 1,-1,-1, 0, 1,undef,  # 28,35
     -1, 0,-1, 0, 0, 1,undef,  0, 0,-1, 0,-1, 0,undef,  # 42,49
     -1, 1, 0,-1,-1, 0,undef, -1,-1, 0, 1,-1,-1,undef,  # 56,63
     0,  1, 1,-1,-1,-1,undef, -1,-1, 1, 1, 0,-1,undef,  # 70,77
     0,  1, 1, 0, 0,-1,undef, -1,-1, 1, 0, 1,-1,undef,  # 84,91
     1,  0, 0, 1, 1,-1,undef,  0, 0, 1,-1, 1, 0,undef,  # 98,105
     1, -1,-1, 1, 1, 0,undef,  1, 1, 0,-1, 0, 1,undef,  # 112,119
     0, -1,-1, 0, 0, 1,undef,  1, 1,-1, 0,-1, 1,undef,  # 126,133
     -1, 0, 0,-1,-1, 1,undef,  0, 0,-1, 1,-1, 0,undef,  # 140,147
     -1, 1, 1,-1,-1, 0,undef, -1,-1, 0, 1, 0,-1);

sub n_to_dxdy {
  my ($self, $n) = @_;
  ### Flowsnake n_to_dxdy(): $n

  if ($n < 0) { return; }
  if (is_infinite($n)) { return ($n,$n); }

  my $int = int($n);
  $n -= $int;  # fraction part
  ### $int
  ### frac: $n

  my $state;
  {
    my $arm = _divrem_mutate ($int, $self->{'arms'});
    $state = 28 * $arm;  # initial rotation

    # adjust so that for arms=2 point N=1 has $int==1
    # or for arms=3 then points N=1 and N=2 have $int==1
    if ($arm) { $int += 1; }
  }
  ### initial state: $state

  my $turn_state = $state;
  my $turn_notlow = 0;
  foreach my $digit (reverse digit_split_lowtohigh($int,7)) { # high to low
    ### $digit
    $state += $digit;

    if ($digit == 6) {
      $turn_notlow = 84;     # is not the least significant digit
    } else {
      $turn_state = $state;  # lowest non-6
      $turn_notlow = 0;      # and is the least significant digit
    }
    $state = $next_state[$state];
  }
  ### int digits state: $state
  ### $turn_state
  ### $turn_notlow

  $state = int($state/14);
  my $di = $dir6_to_di[$state];
  my $dj = $dir6_to_dj[$state];
  ### int direction: "di=$di, dj=$dj"

  # fraction in final $state direction
  if ($n) {
    $turn_state += $turn_notlow;
    my $next_di = $digit_to_next_di[$turn_state];
    my $next_dj = $digit_to_next_dj[$turn_state];

    ### $next_di
    ### $next_dj

    $di += $n*($next_di - $di);
    $dj += $n*($next_dj - $dj);

    ### with frac: "di=$di, dj=$dj"
  }

  ### ret: "dx=".(2*$di+$dj)." dy=$dj"
  return (2*$di+$dj,
          $dj);

}

my @attempt_dx = (0, -2, -1);
my @attempt_dy = (0, 0, -1);
sub xy_to_n {
  my ($self, $x, $y) = @_;
  ### Flowsnake xy_to_n(): "$x, $y"

  $x = round_nearest($x);
  $y = round_nearest($y);
  if (($x + $y) % 2) { return undef; }
  ### round to: "$x,$y"

  my ($n, $cx, $cy);
  foreach my $i (0, 1, 2) {
    if (defined ($n = $self->SUPER::xy_to_n($x + $attempt_dx[$i],
                                            $y + $attempt_dy[$i]))
        && (($cx,$cy) = $self->n_to_xy($n))
        && $x == $cx
        && $y == $cy) {
      return $n;
    }
  }
  return undef;
}

# 0  straight
# 1  +60 rev
# 2  180 rev
# 3  +240
# 4  straight
# 5  straight
# 6  -60 rev

# 4---- 5---- 6
#  \           \
#    3---- 2    7
#        /
# 0---- 1
#
# turn(N) = tdir6(N)-tdir6(N-1)
# N-1 changes low 0s to low 6s
# N   = aaad000
# N-1 = aaac666
# low 0s no change to direction
# low 6s state 7
# N=14=20[7] dir[2]=3,dirrev[0]=5 total 3+5=2mod6
# N-1=13=16[7] dir[1]=1,dirrev[6]=0 total 1+0=1  diff 2-1=1
# dir[2]-dir[1]=2
# dirrev[0] since digit=2 goes to rev
# N=23=32[7]

#            0  1  2  3  4  5
my @turn6 = (1, 2,-1,-2, 0,-1,  # forward
             1, 0, 2, 1,-2,-1,  # reverse
             #
             1, 1,-1,-1, 1,-1,  # 0,0,-1,0,+1,+1,0
             1,-1, 1, 1,-1,-1,  # 0,0,-1,-1,0,+1,0
            );
my @digit_to_reverse = (-1,5,5,undef,-1,-1,5);  # -1=forward,5=reverse
sub _WORKING_BUT_SECRET__n_to_turn6 {
  my ($self, $n) = @_;
  unless ($n >= 1) {
    return undef;
  }
  if (is_infinite($n)) {
    return $n;
  }

  my $lowdigit = _divrem_mutate($n,7);
  ### $lowdigit

  # skip low 0 digits
  unless ($lowdigit) {
    while ($n) {
      last if ($lowdigit = _divrem_mutate($n,7));  # stop at non-zero
    }
    # flag that some zeros were skipped
    $lowdigit += 12;
    ### $lowdigit
  }

  # Forward/reverse reverse from lowest non-3.
  # Digit parts 0,4,5 always forward, 1,2,6 always reverse,
  # 3 is unchanged so following the digit above it.
  for (;;) {
    my $digit = _divrem_mutate($n,7);
    if ($digit != 3) {
      $lowdigit += $digit_to_reverse[$digit];
      last;
    }
  }

  ### lookup: $lowdigit
  return $turn6[$lowdigit];
}

#------------------------------------------------------------------------------
# levels

sub level_to_n_range {
  my ($self, $level) = @_;
  return (0, 7**$level * $self->{'arms'});
}
sub n_to_level {
  my ($self, $n) = @_;
  if ($n < 0) { return undef; }
  if (is_infinite($n)) { return $n; }
  $n = round_nearest($n) - 1;
  _divrem_mutate ($n, $self->{'arms'});
  my ($pow, $exp) = round_down_pow ($n, 7);
  return $exp + 1;
}

#------------------------------------------------------------------------------
1;
__END__

#       4-->5-->6
#       ^       ^
#        \       \
#         3-->2
#            /
#           v
#       0-->1
#
# longest to 6 is x=4,y=2 is 4*4+3*2*2 = 28
#
#             6<---
#             ^
#            /
#       0   5<--4
#        \       \
#         v       v
#         1<--2<--3
#
# longest to 3 is x=5,y=1 is 5*5+3*1*1 = 28
#
# side len 1 len sqrt(7)
# total sqrt(7)^k + ... + 1
#     = (b^(k+1) - 1) / (b - 1)
#     < b^(k+1) / (b - 1)
# squared 7^(k+1) / (7 - 2*sqrt(7) + 1)
#     = 7^k * 7/(7-2*sqrt(7)+1)
#     = 7^k * 2.584
#
# minimum = b^k - upper(k-1)
#         = b^k - b^k / (b - 1)
#         = b^k * (1 - 1/(b-1))
#         = b^k * (b-1 - 1)/(b-1)
#         = b^k * (b-2)/(b-1)
#         = b^k * 0.392
#
# sqrt((x/2)^2 + (y*sqrt(3)/2)^2)
#    = sqrt(x^2/4 + y^2*3/4)
#    = sqrt(x^2 + 3*y^2)/2

# sqrt(x^2 + 3*y^2)/2 > b^k * (b-2)/(b-1)
# sqrt(x^2 + 3*y^2) > b^k * 2*(b-2)/(b-1)
# x^2 + 3*y^2 > 7^k * (2*(b-2)/(b-1))^2
# x^2 + 3*y^2 > 7^k * (2*(b-2)/(b-1))^2
# (x^2 + 3*y^2) / (2*(b-2)/(b-1))^2 > 7^k
# 7^k < (x^2 + 3*y^2) / (2*(b-2)/(b-1))^2
# k < log7 ((x^2 + 3*y^2) / (2*(b-2)/(b-1))^2)
# k < log7 ((x^2 + 3*y^2) * 1.62
# k < log((x^2 + 3*y^2) * 1.62/log(7)
# k < log((x^2 + 3*y^2) * 0.8345



    #                     *---E
    #                    /     \
    #               *---*       *---*
    #              /     \     /     \
    #             *       *---*       *
    #              \     /     \     /
    #               *---*       *---*
    #              /     \     /     \
    #             *       *---*       *
    #              \     /     \     /
    #               *---*       *---*
    #                    \     /
    #                     *---*
    #
    #
    #                     *     *
    #                    / \   / \
    #                   /   \ /   \
    #                  *     *     *
    #                  |     |     |
    #                  |     |     |
    #                  *     *     *
    #                 / \   / \   / \
    #                /   \ /   \ /   \
    #               *     *     *     *
    #               |     |     |     |
    #               |     |     |     |
    #               *     *     *     *
    #                \   / \   / \   /
    #                 \ /   \ /   \ /
    #                  *     *     *
    #                  |     |     |
    #                  |     |     |
    #                  *     *     *
    #                   \   / \   /
    #                    \ /   \ /
    #                     *     *
    #
    #
    #
    #
    #    B
    #   / \   / \
    #  /   \ /   \
    # .  ^  .     .
    # |   | |     |
    # |    ||     |
    # .     O-->  A
    #  \   / \   /
    #   \ / | \ /
    #    . |   .
    #    | v   |
    #    |     |
    #    C     .
    #     \   /
    #      \ /


=for stopwords eg Ryde flowsnake Gosper ie Fukuda Shimizu Nakamura Math-PlanePath Ns zdigit tdigit

=head1 NAME

Math::PlanePath::Flowsnake -- self-similar path through hexagons

=head1 SYNOPSIS

 use Math::PlanePath::Flowsnake;
 my $path = Math::PlanePath::Flowsnake->new;
 my ($x, $y) = $path->n_to_xy (123);

=head1 DESCRIPTION

X<Gosper, William>This path is an integer version of the flowsnake curve by
William Gosper,

=cut

# math-image --path=Flowsnake --all --output=numbers_dash

=pod

                         39----40----41                        8
                           \           \
          32----33----34    38----37    42                     7
            \           \        /     /
             31----30    35----36    43    47----48            6
                  /                    \     \     \
          28----29    17----16----15    44    46    49--..     5
         /              \           \     \  /
       27    23----22    18----19    14    45                  4
         \     \     \        /     /
          26    24    21----20    13    11----10               3
            \  /                    \  /     /
             25     4---- 5---- 6    12     9                  2
                     \           \         /
                       3---- 2     7---- 8                     1
                           /
                    0---- 1                                  Y=0

     X=-4 -3 -2 -1  0  1  2  3  4  5  6  7  8  9 10 11

The points are spread out on every second X coordinate to make little
triangles with integer coordinates, per L<Math::PlanePath/Triangular
Lattice>.

The base pattern is the seven points 0 to 6,

    4---- 5---- 6
     \           \
       3---- 2
           /
    0---- 1

This repeats at 7-fold increasing scale, with sub-sections rotated according
to the edge direction and the 1, 2 and 6 sections in reverse.  The next
level can be seen at the multiple-of-7 points N=0,7,14,21,28,35,42,49.

                                  42
                      -----------    ---
                   35                   ---
       -----------                         ---
    28                                        49 ---
      ---
         ----                  14
             ---   -----------  |
                21              |
                               |
                              |
                              |
                        ---- 7
                   -----
              0 ---

Notice this is the same shape as N=0 to N=6, but rotated by atan(1/sqrt(7))
= 20.68 degrees anti-clockwise.  Each level rotates further and for example
after about 18 levels it goes all the way around and back to the first
quadrant.

The rotation doesn't fill the plane though, only 1/3 of it.  The shape
fattens as it curls around, but leaves a spiral gap beneath (see L</Arms>
below).

=head2 Tiling

The base pattern corresponds to a tiling by hexagons as follows, with the
"***" lines being the base figure.

          .     .
         / \   / \
        /   \ /   \
       .     .     .
       |     |     |
       |     |     |
       4*****5*****6
      /*\   / \   /*\
     / * \ /   \ / * \
    .   * .     .   * .
    |   * |     |    *|
    |    *|     |    *|
    .     3*****2     7...
     \   / \   /*\   /
      \ /   \ / * \ /
       .     . *   .
       |     | *   |
       |     |*    |
       0*****1     .
        \   / \   /
         \ /   \ /
          .     .

In the next level the parts corresponding to 1, 2 and 6 are reversed because
they have their hexagon tile to the right of the line segment, rather than
to the left.

=head2 Arms

The optional C<arms> parameter can give up to three copies of the flowsnake,
each advancing successively.  For example C<arms=E<gt>3> is as follows.

    arms => 3                        51----48----45                 5
                                       \           \
                      ...   69----66    54----57    42              4
                        \     \     \        /     /
       28----25----22    78    72    63----60    39    33----30     3
         \           \     \  /                    \  /     /
          31----34    19    75    12----15----18    36    27        2
               /     /              \           \        /
       40----37    16     4---- 1     9---- 6    21----24           1
      /              \     \              /
    43    55----58    13     7     0---- 3    74----77---...    <- Y=0
      \     \     \     \  /                    \
       46    52    61    10     2     8----11    71----68          -1
         \  /     /              \  /     /           /
          49    64    70----73     5    14    62----65             -2
                  \  /     /           /     /
                   67    76    20----17    59    53----50          -3
                        /     /              \  /     /
                      ...   23    35----38    56    47             -4
                              \     \     \        /
                               26    32    41----44                -5
                                 \  /
                                  29                               -6

                                   ^
       -9 -8 -7 -6 -5 -4 -3 -2 -1 X=0 1  2  3  4  5  6  7  8  9

Notice the N=3*k points are the plain curve, N=3*k+1 is a copy rotated by
120 degrees (1/3 around), and N=3*k+2 is a copy rotated by 240 degrees (2/3
around).  The initial N=1 of the second arm and N=2 of the third correspond
to N=3 of the first arm, rotated around.

Essentially the flowsnake fills an ever expanding hexagon with one corner at
the origin, and wiggly sides.  In the following picture the plain curve
fills "A" and there's room for two more arms to fill B and C, rotated 120
and 240 degrees respectively.

            *---*
           /     \
      *---*   A   *
     /     \     /
    *   B   O---*
     \     /     \
      *---*   C   *
           \     /
            *---*

The sides of these "hexagons" are not straight lines but instead wild wiggly
spiralling S shapes, and the endpoints rotate around by the angle described
above at each level.  Opposing sides are symmetric, so they mesh perfectly
and with three arms fill the plane.

=head2 Fractal

The flowsnake can also be thought of as successively subdividing line
segments with suitably scaled copies of the 0 to 7 figure (or its reversal).

The code here could be used for that by taking points N=0 to N=7^level.  The
Y coordinates should be multiplied by sqrt(3) to make proper equilateral
triangles, then a rotation and scaling to make the endpoint come out at some
desired point, such as X=1,Y=0.  With such a scaling the path is confined to
a finite fractal boundary.

=head1 FUNCTIONS

See L<Math::PlanePath/FUNCTIONS> for behaviour common to all path classes.

=over 4

=item C<$path = Math::PlanePath::Flowsnake-E<gt>new ()>

=item C<$path = Math::PlanePath::Flowsnake-E<gt>new (arms =E<gt> $a)>

Create and return a new flowsnake path object.

The optional C<arms> parameter gives between 1 and 3 copies of the curve
successively advancing.

=item C<($x,$y) = $path-E<gt>n_to_xy ($n)>

Return the X,Y coordinates of point number C<$n> on the path.  Points begin
at 0 and if C<$n E<lt> 0> then the return is an empty list.

Fractional positions give an X,Y position along a straight line between the
integer positions.

=item C<($n_lo, $n_hi) = $path-E<gt>rect_to_n_range ($x1,$y1, $x2,$y2)>

In the current code the returned range is exact, meaning C<$n_lo> and
C<$n_hi> are the smallest and biggest in the rectangle, but don't rely on
that yet since finding the exact range is a touch on the slow side.  (The
advantage of which though is that it helps avoid very big ranges from a
simple over-estimate.)

=back

=head2 Level Methods

=over

=item C<($n_lo, $n_hi) = $path-E<gt>level_to_n_range($level)>

Return C<(0, 7**$level)>, or for multiple arms return C<(0, $arms *
7**$level)>.

There are 7^level + 1 points in a level, numbered starting from 0.  On the
second and third arms the origin is omitted (so as not to repeat that point)
and so just 7^level for them, giving 7^level+1 + (arms-1)*7^level =
arms*7^level + 1 many points starting from 0.

=back

=head1 FORMULAS

=head2 N to X,Y

The position of a given N can be calculated from the base-7 digits of N from
high to low.  At a given digit position the state maintained is

    direction 0 to 5, multiple of 60-degrees
    plain or reverse pattern

It's convenient to work in the "i,j" coordinates per
L<Math::PlanePath/Triangular Calculations>.  This represents a point in the
triangular grid as i+j*w where w=1/2+I*sqrt(3)/2 the a complex sixth root of
unity at +60 degrees.

    foreach base-7 digit high to low
      (i,j) = (2i-j, i+3j)   # multiply by 2+w
      (i,j) += position of digit in plain or reverse,
               and rotated by "direction"

The multiply by 2+w scales up i,j by that vector, so for instance i=1,j=0
becomes i=2,j=1.  This spreads the points as per the multiple-of-7 diagram
shown above, so what was at N scales up to 7*N.

The digit is then added as either the plain or reversed base figure,

      plain             reverse

    4-->5-->6
    ^       ^
     \       \
      3-->2   *             6<---*
         /                  ^
        v                  /
    0-->1             0   5<--4
                       \       \
                        v       v
                        1<--2<--3


The arrow shown in each part is whether the state becomes plain or reverse.
For example in plain state at digit=1 the arrow points backwards so if
digit=1 then the state changes to reverse for the next digit.  The direction
likewise follows the direction of each segment in the pattern.

Notice the endpoint "*" is at at 2+w in both patterns.  When considering the
rotation it's convenient to reckon the direction by this endpoint.

The combination of direction and plain/reverse is a total of 14 different
states, and for each there's 7 digit values (0 to 6) for a total 84 i,j.

=head2 X,Y to N

The current approach uses the C<FlowsnakeCentres> code.  The tiling in
C<Flowsnake> and C<FlowsnakeCentres> is the same so the X,Y of a given N are
no more than 1 away in the grid of the two forms.

The way the two lowest shapes are arranged in fact means that if the
Flowsnake N is at X,Y then the same N in C<FlowsnakeCentres> is at one of
three locations

    X, Y         same
    X-2, Y       left      (+180 degrees)
    X-1, Y-1     left down (+240 degrees)

This is true even when the rotated "arms" multiple paths (the same number of
arms in both paths).

Is there an easy way to know which of the three offsets is right?  The
current approach is to put each through C<FlowsnakeCentres> to make an N,
and put that N back through Flowsnake C<n_to_xy()> to see if it's the target
C<$n>.

=head2 Rectangle to N Range

The current code calculates an exact C<rect_to_n_range()> by searching for
the highest and lowest Ns which are in the rectangle.

The curve at a given level is bounded by the Gosper island shape but the
wiggly sides make it difficult to calculate, so a bounding radius
sqrt(7)^level, plus a bit, is used.  The portion of the curve comprising a
given set of high digits of N can be excluded if the N point is more than
that radius away from the rectangle.

When a part of the curve is excluded it prunes a whole branch of the digits
tree.  When the lowest digit is reached then a check for that point being
actually within the rectangle is made.  The radius calculation is a bit
rough, and it doesn't take into account the direction of the curve, so it's
a rather large over-estimate, but it works.

The same sort of search can be applied for highest and lowest N in a
non-rectangular shapes, calculating a radial distance away from the shape.
The distance calculation doesn't have to be exact either, it can go from
something bounding the shape until the lowest digit is reached and an
individual X,Y is considered as a candidate for high or low N.

=head2 N to Turn

The turn made by the curve at a point NE<gt>=1 can be calculated from the
lowest non-0 digit and the plain/reverse state per the lowest non-3 above
there.

   N digits in base 7
   strip low 0 digits, zcount many of them
   zdigit = take low digit
   strip low 3 digits
   tdigit = take low digit (0 if no digits left)
   plain if tdigit=0,4,5, reverse if tdigit=1,2,6

             if zcount=0       if zcount>=1
             ie. no low 0s     ie. some low 0s
   zdigit   plain reverse     plain   reverse
   ------   ----- -------     -----   -------
      1        1      1          1        1
      2        2      0          1       -1       turn left
      3       -1      2         -1        1       multiple of
      4       -2      1         -1        1       60 degrees
      5        0     -2          1       -1
      6        1     -1         -1       -1

For example N=9079 is base-7 "35320" so a single low 0 for zcount=1 and
strip it to "3532".  Take zdigit=2 leaving "353".  Skip low 3s leaving "35".
Take tdigit=5 which is "plain".  So table "plain" with zcount>=1 is the
third column and there zdigit=2 is turn=+1.

The turns in the zcount=0 "no low 0s" columns follow the turns of the base
pattern shown above.  For example zdigit=1 is as per N=1 turning 120 degrees
left, so +2.  For the reverse pattern the turns are negated and the zdigit
value reversed, so the "reverse" column read 6 to 1 is the same as the plain
column negated and read 1 to 6.

Low 0s are stripped because the turn at a point such as N=7 ("10" in base-7)
is determined by the pattern above it, the self-similar multiple-of-7s
shape.  But when there's low 0s in the way there's an adjustment to apply
because the last segment of the base pattern is not in the same direction as
the first, but instead at -60 degrees.  Likewise the first segment of the
reverse pattern.  At some zdigit positions those two cancel out, such as at
zdigit=1 where a plain and reverse meet, but others it's not so and hence
separate table columns for with or without low 0s.

The plain or reverse pattern is determined by the lowest non-3 digit.  This
works because the digit=0, digit=4, and digit=5 segments of the base pattern
have their sub-parts "plain" in all cases, both the plain and reverse forms.
Conversely digit=1, digit=2 and digit=6 segments are "reverse" in all cases,
both plain and reverse forms.  But the digit=3 part is plain in plain and
reverse in reverse, so it inherits the orientation of the digit above and
it's therefore necessary to skip across any 3s.

When taking digits, N is treated as having infinite 0-digits at the high
end.  This only affects the tdigit plain/reverse step.  If N has a single
non-zero such as "5000" then it's taken as zdigit=5 and above that for the
plain/reverse a tdigit=0 is then assumed.  The first turn is at N=1 so
there's always at least one non-0 for the zdigit.

=head1 SEE ALSO

L<Math::PlanePath>,
L<Math::PlanePath::FlowsnakeCentres>,
L<Math::PlanePath::GosperIslands>

L<Math::PlanePath::KochCurve>,
L<Math::PlanePath::HilbertCurve>,
L<Math::PlanePath::PeanoCurve>,
L<Math::PlanePath::ZOrderCurve>

"New Gosper Space Filling Curves" by Fukuda, Shimizu and Nakamura.  On
flowsnake variations in bigger hexagons (with wiggly sides too).

=over

L<http://kilin.clas.kitasato-u.ac.jp/museum/gosperex/343-024.pdf> or
L<http://web.archive.org/web/20070630031400/http://kilin.u-shizuoka-ken.ac.jp/museum/gosperex/343-024.pdf>

=back

=head1 HOME PAGE

L<http://user42.tuxfamily.org/math-planepath/index.html>

=head1 LICENSE

Copyright 2010, 2011, 2012, 2013, 2014 Kevin Ryde

This file is part of Math-PlanePath.

Math-PlanePath is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the Free
Software Foundation; either version 3, or (at your option) any later
version.

Math-PlanePath is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for
more details.

You should have received a copy of the GNU General Public License along with
Math-PlanePath.  If not, see <http://www.gnu.org/licenses/>.

=cut