libmath-planepath-perl 122-1 / usr / share / perl5 / Math / PlanePath / DekkingCurve.pm

# Copyright 2011, 2012, 2013, 2014, 2015 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/>.


package Math::PlanePath::DekkingCurve;
use 5.004;
use strict;
#use List::Util 'max';
*max = \&Math::PlanePath::_max;

use vars '$VERSION', '@ISA';
$VERSION = 122;
use Math::PlanePath;
use Math::PlanePath::Base::NSEW;
@ISA = ('Math::PlanePath::Base::NSEW',
        'Math::PlanePath');
*_divrem = \&Math::PlanePath::_divrem;
*_divrem_mutate = \&Math::PlanePath::_divrem_mutate;

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

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


use constant n_start => 0;
use constant class_x_negative => 1;
use constant class_y_negative => 1;

use constant parameter_info_array => [ { name      => 'arms',
                                         share_key => 'arms_4',
                                         display   => 'Arms',
                                         type      => 'integer',
                                         minimum   => 1,
                                         maximum   => 4,
                                         default   => 1,
                                         width     => 1,
                                         description => 'Arms',
                                       } ];


#------------------------------------------------------------------------------

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

sub y_negative {
  my ($self) = @_;
  return $self->{'arms'} > 2;
}
{
  my @y_negative_at_n = (undef, undef, undef, 8, 3);
  sub y_negative_at_n {
    my ($self) = @_;
    return $y_negative_at_n[$self->{'arms'}];
  }
}

#------------------------------------------------------------------------------

use Math::PlanePath::DekkingCentres;
use vars '@_next_state','@_digit_to_x','@_digit_to_y','@_yx_to_digit';
BEGIN {
  *_next_state = \@Math::PlanePath::DekkingCentres::_next_state;
  *_digit_to_x = \@Math::PlanePath::DekkingCentres::_digit_to_x;
  *_digit_to_y = \@Math::PlanePath::DekkingCentres::_digit_to_y;
  *_yx_to_digit = \@Math::PlanePath::DekkingCentres::_yx_to_digit;
}

sub new {
  my $self = shift->SUPER::new(@_);
  $self->{'arms'} ||= 1;
  return $self;
}

# tables generated by tools/dekking-curve-table.pl
#
my @edge_dx = (0,0,0,1,1, 0,0,1,1,0, 0,0,0,1,0, 0,0,1,0,1, 0,1,0,1,1,
               1,1,1,1,1, 1,1,1,0,1, 1,1,0,1,0, 0,0,1,0,0, 0,1,1,0,0,
               1,1,1,0,0, 1,1,0,0,1, 1,1,1,0,1, 1,1,0,1,0, 1,0,1,0,0,
               0,0,0,0,0, 0,0,0,1,0, 0,0,1,0,1, 1,1,0,1,1, 1,0,0,1,1,
               1,1,1,1,1, 1,0,0,0,0, 1,1,1,1,1, 0,0,0,0,1, 1,0,0,1,1,
               1,1,1,0,0, 1,1,1,1,1, 0,0,0,1,1, 0,0,1,0,1, 0,1,0,1,1,
               0,0,0,0,0, 0,1,1,1,1, 0,0,0,0,0, 1,1,1,1,0, 0,1,1,0,0,
               0,0,0,1,1, 0,0,0,0,0, 1,1,1,0,0, 1,1,0,1,0, 1,0,1,0,0);
my @edge_dy = (0,0,0,0,0, 0,0,0,1,0, 0,0,1,0,1, 1,1,0,1,1, 1,0,0,1,1,
               0,0,0,1,1, 0,0,1,1,0, 0,0,0,1,0, 0,0,1,0,1, 0,1,0,1,1,
               1,1,1,1,1, 1,1,1,0,1, 1,1,0,1,0, 0,0,1,0,0, 0,1,1,0,0,
               1,1,1,0,0, 1,1,0,0,1, 1,1,1,0,1, 1,1,0,1,0, 1,0,1,0,0,
               0,0,0,1,1, 0,0,0,0,0, 1,1,1,0,0, 1,1,0,1,0, 1,0,1,0,0,
               1,1,1,1,1, 1,0,0,0,0, 1,1,1,1,1, 0,0,0,0,1, 1,0,0,1,1,
               1,1,1,0,0, 1,1,1,1,1, 0,0,0,1,1, 0,0,1,0,1, 0,1,0,1,1,
               0,0,0,0,0, 0,1,1,1,1, 0,0,0,0,0, 1,1,1,1,0, 0,1,1,0,0);

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

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

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

  my $arms = $self->{'arms'};
  my $arm = _divrem_mutate ($int, $arms);
  if ($arm) { $int += 1; }
  my @digits = digit_split_lowtohigh($int,25);
  my $state = 0;
  my @x;
  my @y;
  foreach my $i (reverse 0 .. $#digits) {
    $state += $digits[$i];
    $x[$i] = $_digit_to_x[$state];
    $y[$i] = $_digit_to_y[$state];
    $state = $_next_state[$state];
  }

  ### @x
  ### @y
  ### $state
  ### dx: $_digit_to_x[$state+24] - $_digit_to_x[$state]
  ### dy: $_digit_to_y[$state+24] - $_digit_to_y[$state]

  my $zero = $int * 0;
  my $x = ($n * (($_digit_to_x[$state+24] - $_digit_to_x[$state])/4)
           + digit_join_lowtohigh(\@x, 5, $zero)
           + $edge_dx[$state]);
  my $y = ($n * (($_digit_to_y[$state+24] - $_digit_to_y[$state])/4)
           + digit_join_lowtohigh(\@y, 5, $zero)
           + $edge_dy[$state]);

  if ($arm < 2) {
    if ($arm < 1) { return ($x,$y); }  # arm==0
    return (-$y,$x);                   # arm==1  rotate +90
  }
  if ($arm < 3) { return (-$x,-$y); }  # arm==2
  return ($y,-$x);                     # arm==3  rotate -90
}

sub xy_to_n {
  my ($self, $x, $y) = @_;
  ### DekkingCurve xy_to_n(): "$x, $y"

  $x = round_nearest ($x);
  $y = round_nearest ($y);
  if (is_infinite($x)) {
    return $x;
  }
  if (is_infinite($y)) {
    return $y;
  }

  my $arms = $self->{'arms'};
  if (($arms < 2 && $x < 0) || ($arms < 3 && $y < 0)) {
    ### X or Y negative, no N value ...
    return undef;
  }

  foreach my $arm (0 .. $arms-1) {
    foreach my $xoffset (0,-1) {
      foreach my $yoffset (0,-1) {

        my @x = digit_split_lowtohigh($x+$xoffset,5);
        my @y = digit_split_lowtohigh($y+$yoffset,5);
        my $state = 0;
        my @n;
        foreach my $i (reverse 0 .. max($#x,$#y)) {
          my $digit = $n[$i] = $_yx_to_digit[$state + 5*($y[$i]||0) + ($x[$i]||0)];
          $state = $_next_state[$state+$digit];
        }
        my $zero = $x*0*$y;
        my $n = digit_join_lowtohigh(\@n, 25, $zero);
        $n = $n*$arms;
        if (my ($nx,$ny) = $self->n_to_xy($n)) {
          if ($nx == $x && $ny == $y) {
            return $n + ($arm ? $arm-$arms : $arm);
          }
        }
      }
    }
    ($x,$y) = ($y,-$x); # rotate -90
    ### rotate to: "$x, $y"
  }
  return undef;
}

# not exact
sub rect_to_n_range {
  my ($self, $x1,$y1, $x2,$y2) = @_;
  ### DekkingCurve rect_to_n_range(): "$x1,$y1, $x2,$y2"

  $x1 = round_nearest ($x1);
  $x2 = round_nearest ($x2);
  $y1 = round_nearest ($y1);
  $y2 = round_nearest ($y2);

  if ($x1 > $x2) { ($x1,$x2) = ($x2,$x1); }
  if ($y1 > $y2) { ($y1,$y2) = ($y2,$y1); }

  my $arms = $self->{'arms'};
  if (($arms < 2 && $x2 < 0) || ($arms < 3 && $y2 < 0)) {
    ### rectangle all negative, no N values ...
    return (1, 0);
  }

  my ($pow) = round_down_pow (max(abs($x1),abs($y1),$x2,$y2) + 1, 5);
  ### $pow
  return (0, 25*$pow*$pow*$arms - 1);
}

#------------------------------------------------------------------------------

sub level_to_n_range {
  my ($self, $level) = @_;
  return (0,  25**$level * $self->{'arms'});
}
sub n_to_level {
  my ($self, $n) = @_;
  ### n_to_level(): $n
  if ($n < 0) { return undef; }
  if (is_infinite($n)) { return $n; }
  $n = round_nearest($n);
  $n += $self->{'arms'}-1;  # division rounding up
  _divrem_mutate ($n, $self->{'arms'});
  my ($pow, $exp) = round_up_pow ($n, 25);
  return $exp;
}

#------------------------------------------------------------------------------
# Not taking into account multiple arms ...

# Return true if X axis segment $x to $x+1 is traversed
sub _UNDOCUMENTED__xseg_is_traversed {
  my ($self, $x) = @_;
  if ($x < 0 || is_infinite($x)) { return 0; }
  if ($x == 0) { return 1; }
  my $digit = _divrem_mutate($x, 5);
  if ($digit) {
    return ($digit == 1);
  }
  # find lowest non-zero
  while ($x && ! ($digit = _divrem_mutate($x, 5))) { }
  return ($digit == 1 || $digit == 2);
}

# Return true if Y axis segment $y to $y+1 is traversed
sub _UNDOCUMENTED__yseg_is_traversed {
  my ($self, $y) = @_;
  if ($y < 0 || is_infinite($y)) { return 0; }
  my $digit = _divrem_mutate($y, 5);
  if ($digit != 4) {
    return ($digit == 3);
  }
  # find lowest non-4
  while ($y && ($digit = _divrem_mutate($y, 5)) == 4) { }
  return ($digit == 2 || $digit == 3);
}

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

=for stopwords eg Ryde ie Math-PlanePath Dekking

=head1 NAME

Math::PlanePath::DekkingCurve -- 5x5 self-similar edge curve

=head1 SYNOPSIS

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

=head1 DESCRIPTION

This is an integer version of a 5x5 self-similar curve per

=over

F. M. Dekking, "Recurrent Sets", Advances in Mathematics, volume 44, 1982,
pages 79-104, section 4.9 "Gosper-Type Curves"

=back

The base pattern is N=0 to N=25.  It repeats with rotations or reversals
which make the ends join.  For example N=75 to N=100 is the base pattern in
reverse, ie. from N=25 down to N=0.  Or N=50 to N=75 is reverse and also
rotate by -90.

=cut

# math-image --path=DekkingCurve --all --output=numbers_dash --size=78x30

=pod

    10 |             123-124-125-...      86--85
       |               |                   |   |
     9 | 115-116-117 122-121  90--89--88--87  84
       |   |       |       |   |               |
     8 | 114-113 118-119-120  91--92--93  82--83
       |       |                       |   |
     7 |     112 107-106 103-102  95--94  81  78--77
       |       |   |   |   |   |   |       |   |   |
     6 |     111 108 105-104 101  96--97  80--79  76
       |       |   |           |       |           |
     5 |     110-109  14--15 100--99--98  39--40  75          66--65
       |               |   |               |   |   |           |   |
     4 |  10--11--12--13  16  35--36--37--38  41  74  71--70  67  64
       |   |               |   |               |   |   |   |   |   |
     3 |   9---8---7  18--17  34--33--32  43--42  73--72  69--68  63
       |           |   |               |   |                       |
     2 |       5---6  19  22--23  30--31  44  47--48  55--56--57  62--61
       |       |       |   |   |   |       |   |   |   |       |       |
     1 |       4---3  20--21  24  29--28  45--46  49  54--53  58--59--60
       |           |           |       |           |       |
    Y=0|   0---1---2          25--26--27          50--51--52
       +----------------------------------------------------------------
         X=0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15

The curve segments correspond to edges of squares in a 5x5 arrangement.

     +- - -+- - -+- - 14----15  ---+
     |     |     |     |  v  |>    |
        ^     ^       <|     |
    10----11----12----13- - 16   --+
     |              v        |>    |
     |>       ^           ^  |
     9-----8-----7 -- 18----17   --+
        v  |     |     |>          |
     |        ^  |>    |        ^
     +- -  5-----6 -  19    22----23
           |          <|     |    <|
     |    <|  ^        |    <|     |
     +- -  4-----3    20----21 -- 24
                 |       v        <|
        ^     ^  |>    |     |     |
     0-----1-----2  -- + -- -+-   25

The little notch marks show which square each edge represents.  This is the
side the curve expands into at the next level.  For example N=1 to N=2 has
its notch on the left so the next level N=25 to N=50 expands on the left.

All the directions are made by rotating the base pattern.  When the
expansion is on the right the segments go in reverse.  For example N=2 to
N=3 expands on the right and is made by rotating the base pattern clockwise
90 degrees.  This means that N=2 becomes the 25 end, and following the curve
to the 0 start at N=3.

Dekking writes these directions as a sequence of 25 symbols s(i,j) where i=0
for left plain or i=1 for right reverse and j=0,1,2,3 direction j*90 degrees
anti-clockwise so E,N,W,S.

=head2 Arms

The optional C<arms> parameter can give up to four copies of the curve, each
advancing successively.  Each copy is in a successive quadrant.

=cut

# math-image --path=DekkingCurve,arms=3  --expression='i<75?i:0' --output=numbers_dash --size=78x24

=pod

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

The origin is N=0 and is on the first arm only.  The second and subsequent
arms begin 1,2,etc.  The curves interleave perfectly on the axes where the
arms meet.  The result is that arms=4 fills the plane visiting each integer
X,Y exactly once and not touching or crossing.

=head1 FUNCTIONS

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

=over 4

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

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

Create and return a new path object.

The optional C<arms> parameter gives between 1 and 4 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.

=back

=head2 Level Methods

=over

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

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

There are 25^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 25^level for them, giving 25^level+1 + (arms-1)*25^level =
arms*25^level + 1 many points starting from 0.

=back

=head1 FORMULAS

=head2 X Axis Segments

In the sample points above there are some line segments on the X axis.
A segment X to X+1 is traversed or not according to

    X digits in base 5

    traversed        if X==0
    traversed        if low digit 1
    not-traversed    if low digit 2 or 3 or 4
    when low digit == 0
      traversed      if lowest non-zero 1 or 2
      not-traversed  if lowest non-zero 3 or 4

In the samples the segments at X=1, X=6 and X=11 segments traversed are low
digit 1.  Their preceding X=5 and X=10 segments are low digit==0 and the
lowest non-zero 1 or 2 (respectively).  At X=15 however the lowest non-zero
is 3 and so not-traversed there.

In general in groups of 5 there is always X==1 mod 5 traversed but its
preceding X==0 mod 5 is traversed or not according to lowest non-zero 1,2 or
3,4.

This pattern is found by considering how the base pattern expands.  The
plain base pattern has its south edge on the X axis.  The first two
sub-parts of that south edge are the base pattern unrotated, so the south
edge again, but the other parts rotated.  In general the sides are

           0 1 2 3 4
    S  ->  S,S,E,N,W
    E  ->  S,S,E,N,N
    N  ->  W,S,E,N,N
    W  ->  W,S,E,N,W

Starting in S and taking digits high to low a segment is traversed when the
final state is S again.

Any digit 1,2,3 goes to S,E,N respectively.  If no 1,2,3 at all then S
start.  At the lowest 1,2,3 there are only digits 0,4 below.  If no such
digits then only digit 1 which is S, or no digits at all for N=0, is
traversed.  If one or more 0s below then E goes to S so a lowest non-zero 2
means traversed too.  If there is any 4 then it goes to N or W and in those
states both 0,4 stay in N or W so not-traversed.

The transitions from the lowest 1,2,3 can be drawn in a state diagram,

               +--+
               v  |4                           base 5 digits of X
               North  <---+    <-------+       high to low
             /            |            |
            /0            |4           |
           /              |            |3
    +->   v               |       2    |  
    |   West             East  <--- start lowest 1,2,3
    +--   ^               |            |
    0,4    \              |            |1
            \4            |0           |or no 1,2,3 at all
             \            |            |
               South  <---+    <-------+
               ^  |0
               +--+

The full diagram, starting from the top digit, is less clear

               +--+
               v  |3,4
        +--->  North  <---+
       3|    /  | ^  \    |3,4
        |   /0  1 |  2\   |              base 5 digits of X
        |  /    | |    \  |              high to low
    +-> | v     | |     v |   <-+
    |   West 2---------> East   |        start in South,
    +-- | ^     | |     ^ |   --+        segment traversed
    0,4 |  \    | |    /  |    2         if end in South
        |   \4  | 3  2/   |
       1|    \  v |  /    |0,1
        +--->  South  <---+
               ^  |0,1
               +--+

but allows usual DFA state machine manipulations to reverse to go low to
high.

          +---------- start ----------+
          |       1    0|   2,3,4     |         base 5 digits of X
          |             |             |         low to high
          v       1,2   v   3,4       v
    traversed <------- m1 -------> not-traversed
                      0| ^
                       +-+

In state m1 a 0 digit loops back to m1 so finds the lowest non-zero.  States
start and m1 are the same except for the behaviour of digit 2 and so in the
rules above the result for digit 2 differs according to whether there are
any low 0s.

=head2 Y Axis Segments

The Y axis can be treated similarly

    Y digits in base 5  (with a single 0 digit if Y==0)

    traversed        if lowest digit 3
    not-traversed    if lowest digit 0 or 1 or 2
    when lowest digit == 4
      traversed      if lowest non-4 is 2 or 3
      not-traversed  if lowest non-4 is 0 or 1

The Y axis goes around the base square clockwise, so the digits are reversed
0E<lt>-E<gt>4 from the X axis for the state transitions.  The initial state
is W.

           0 1 2 3 4
    S  ->  W,N,E,S,S
    E  ->  N,N,E,S,S
    N  ->  N,N,E,S,W
    W  ->  W,N,E,S,W

N and W can be merged as equivalent.  Their only difference is digit 0 going
to N or W and both of those are final result not-traversed.

Final state S is reached if the lowest digit is 3, or if state S or E are
reached by digit 2 or 3 and then only 4s below.

=head2 X,Y Axis Interleaving

For arms=2 the second copy of the curve is rotated +90 degrees, and
similarly a third or fourth copy in arms=3 or 4.  This means each axis is a
Y axis of the quadrant before and an X axis of the quadrant after.  When
this happens the segments do not overlap nor does the curve touch.

This is seen from the digit rules above.  The 1 mod 5 segment is always
traversed by X and never by Y.  The 2 mod 5 segment is never traversed by
either.  The 3 mod 5 segment is always traversed by Y and never by X.

The 0 mod 5 segment is sometimes traversed by X, and never by Y.  The 4 mod
5 segment is sometimes traversed by Y, and never by Y.

        0       1       2       3       4
    *-------*-------*-------*-------*-------*
        X       X    neither    Y       Y
      maybe                            maybe

A 4 mod 5 segment has one or more trailing 4s and at +1 for the next segment
they become 0s and increment the lowest non-4.

    +--------+-----+-------+                 
    |  ...   |  d  | 4...4 |   N   == 4 mod 5    X never
    +--------+-----+-------+                     Y maybe

    +--------+-----+-------+
    |  ...   | d+1 | 0...0 |   N+1 == 0 mod 5    X maybe
    +--------+-----+-------+                     Y never

Per the Y rule, a 4 mod 5 segment is traversed when d=2,3.  The following
segment is then d+1=3,4 as lowest non-zero which in the X rule is
not-traversed.  Conversely in the Y rule not-traversed when d=0,1 which
becomes d+1=1,2 which in the X rule is traversed.

So exactly one of two consecutive 4 mod 5 and 0 mod 5 segments are
traversed.

=head1 SEE ALSO

L<Math::PlanePath>,
L<Math::PlanePath::DekkingCentres>,
L<Math::PlanePath::CincoCurve>,
L<Math::PlanePath::PeanoCurve>

=head1 HOME PAGE

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

=head1 LICENSE

Copyright 2011, 2012, 2013, 2014, 2015 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