# Copyright 2011, 2012, 2013, 2014, 2015, 2016 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::KochSquareflakes;
use 5.004;
use strict;

use vars '$VERSION', '@ISA';
$VERSION = 123;
use Math::PlanePath;
@ISA = ('Math::PlanePath');
*_divrem = \&Math::PlanePath::_divrem;

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

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


use constant n_frac_discontinuity => 0;

use constant parameter_info_array =>
  [ { name        => 'inward',
      display     => 'Inward',
      type        => 'boolean',
      default     => 0,
      description => 'Whether to direct the sides of the square inward, rather than outward.',
    } ];

use constant x_negative_at_n => 1;
use constant y_negative_at_n => 1;
use constant sumabsxy_minimum => 1;
use constant rsquared_minimum => 0.5; # minimum X=0.5,Y=0.5

# jump across rings is South-West, so
use constant dx_maximum => 1;
use constant dy_maximum => 1;
use constant dsumxy_maximum => 2; # diagonal NE
use constant ddiffxy_maximum => 2;
use constant ddiffxy_minimum => -2;
use constant dir_maximum_dxdy => (1,-1); # South-East

# N=1,2,3,4  gcd(1/2,1/2) = 1/2
use constant gcdxy_minimum => 1/2;

use constant turn_any_straight => 0; # never straight


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

# level 0 inner square
# sidelen = 4^level
# ring points 4*4^level
# Nend = 4 * [ 1 + ... + 4^level ]
#      = 4 * (4^(level+1) - 1) / 3
#      = (4^(level+2) - 4) / 3
# Nstart = Nend(level-1) + 1
#        = (4^(level+1) - 4) / 3 + 1
#        = (4^(level+1) - 4 + 3) / 3
#        = (4^(level+1) - 1) / 3
#
# level    Nstart             Nend
#    0     (4-1)/3=1          (16-4)/3=12/3=4
#    1     (16-1)/3=15/3=5    (64-4)/3=60/3=20
#    2     (64-1)/3=63/3=21   (256-4)/3=252/3=84
#    3     (256-1)/3=255/3=85
#

sub n_to_xy {
  my ($self, $n) = @_;
  ### KochSquareflakes n_to_xy(): $n
  if ($n < 1) { return; }

  my $frac;
  {
    my $int = int($n);
    $frac = $n - $int;
    $n = $int;       # BigFloat int() gives BigInt, use that
  }

  # (4^(level+1) - 1) / 3 = N
  # 4^(level+1) - 1 = 3*N
  # 4^(level+1) = 3*N+1
  #
  my ($pow,$level) = round_down_pow (3*$n + 1, 4);
  ### $level
  ### $pow
  if (is_infinite($level)) { return ($level,$level); }

  # Nstart = (4^(level+1)-1)/3 with $power=4^(level+1) here
  #
  $n -= ($pow-1)/3;

  ### base: ($pow-1)/3
  ### next base would be: (4*$pow-1)/3
  ### n remainder from base: $n

  my $sidelen = $pow/4;
  (my $rot, $n) = _divrem ($n, $sidelen);  # high part is rot
  ### $sidelen
  ### n remainder: $n
  ### $rot

  ### assert: $n>=0
  ### assert: $n < 4 ** $level

  my @horiz = (1);
  my @diag = (1);
  my $i = 0;
  while (--$level > 0) {
    $horiz[$i+1] = 2*$horiz[$i] + 2*$diag[$i];
    $diag[$i+1]  = $horiz[$i] + 2*$diag[$i];
    ### horiz: $horiz[$i+1]
    ### diag: $diag[$i+1]
    $i++;
  }

  ### horiz: join(', ',@horiz)
  ### $i
  my $x
    = my $y
      = ($n * 0) + $horiz[$i]/-2;  # inherit bignum
  if ($rot & 1) {
    ($x,$y) = (-$y,$x);
  }
  if ($rot & 2) {
    $x = -$x;
    $y = -$y;
  }
  $rot *= 2;

  my $inward = $self->{'inward'};
  my @digits = digit_split_lowtohigh($n,4);

  while ($i > 0) {
    $i--;
    my $digit = $digits[$i] || 0;

    my ($dx, $dy, $drot);
    if ($digit == 0) {
      $dx = 0;
      $dy = 0;
      $drot = 0;
    } elsif ($digit == 1) {
      if ($rot & 1) {
        $dx = $diag[$i];
        $dy = $diag[$i];
      } else {
        $dx = $horiz[$i];
        $dy = 0;
      }
      $drot = ($inward ? 1 : -1);
    } elsif ($digit == 2) {
      if ($rot & 1) {
        if ($inward) {
          $dx = $diag[$i];
          $dy = $diag[$i] + $horiz[$i];
        } else {
          $dx = $diag[$i] + $horiz[$i];
          $dy = $diag[$i];
        }
      } else {
        $dx = $horiz[$i] + $diag[$i];
        $dy = $diag[$i];
        unless ($inward) { $dy = -$dy; }
      }
      $drot = ($inward ? -1 : 1);
    } elsif ($digit == 3) {
      if ($rot & 1) {
        $dx = $dy = $diag[$i] + $horiz[$i];
      } else {
        $dx = $horiz[$i] + 2*$diag[$i];
        $dy = 0;
      }
      $drot = 0;
    }
    ### delta: "$dx,$dy   rot=$rot   drot=$drot"

    if ($rot & 2) {
      ($dx,$dy) = (-$dy,$dx);
    }
    if ($rot & 4) {
      $dx = -$dx;
      $dy = -$dy;
    }
    ### delta with rot: "$dx,$dy"

    $x += $dx;
    $y += $dy;
    $rot += $drot;
  }

  {
    my $dx = $frac;
    my $dy = ($rot & 1 ? $frac : 0);
    if ($rot & 2) {
      ($dx,$dy) = (-$dy,$dx);
    }
    if ($rot & 4) {
      $dx = -$dx;
      $dy = -$dy;
    }
    $x = $dx + $x;
    $y = $dy + $y;
  }

  return ($x,$y);
}

my @inner_to_n = (1,2,4,3);

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

  # +/- 0.75
  if (4*$x < 3 && 4*$y < 3 && 4*$x >= -3 && 4*$y >= -3) {
    return $inner_to_n[($x >= 0) + 2*($y >= 0)];
  }

  $x = round_nearest($x);
  $y = round_nearest($y);

  # quarter curve segment and high digit
  my $n;
  {
    my $negx = -$x;
    if (($y > 0 ? $x > $y : $x >= $y)) {
      ### below leading diagonal ...
      if ($negx > $y) {
        ### bottom quarter ...
        $n = 1;
      } else {
        ### right quarter ...
        $n = 2;
        ($x,$y) = ($y, $negx);   # rotate -90
      }
    } else {
      ### above leading diagonal
      if ($y > $negx) {
        ### top quarter ...
        $n = 3;
        $x = $negx;   # rotate 180
        $y = -$y;
      } else {
        ### right quarter ...
        $n = 4;
        ($x,$y) = (-$y, $x);   # rotate +90
      }
    }
  }
  $y = -$y;
  ### rotate to: "$x,$y   n=$n"

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

  my @horiz;
  my @diag;
  my $horiz = 1;
  my $diag = 1;
  for (;;) {
    push @horiz, $horiz;
    push @diag, $diag;
    my $offset = $horiz+$diag;
    my $nextdiag = $offset + $diag;  # horiz + 2*diag
    ### $horiz
    ### $diag
    ### $offset
    ### $nextdiag

    if ($y <= $nextdiag) {
      ### found level at: "top=$nextdiag vs y=$y"
      $y -= $offset;
      $x += $offset;
      last;
    }
    $horiz = 2*$offset;  # 2*horiz+2*diag
    $diag = $nextdiag;
  }
  ### base subtract to: "$x,$y"


  if ($self->{'inward'}) {
    $y = -$y;
    ### inward invert to: "$x,$y"
  }

  ### origin based side: "$x,$y   horiz=$horiz diag=$diag  with levels ".scalar(@horiz)

  # loop 4*1, 4*4, 4*4^2 etc, extra +1 on the digits to include that in the sum
  #
  my $slope;
  while (@horiz) {
    ### at: "$x,$y slope=".($slope||0)." n=$n"
    $horiz = pop @horiz;
    $diag = pop @diag;
    $n *= 4;

    if ($slope) {
      if ($y < $diag) {
        ### slope digit 0 ...
        $n += 1;
      } else {
        $x -= $diag;
        $y -= $diag;
        ### slope not digit 0, move to: "$x,$y"

        if ($y < $horiz) {
          ### digit 1 ...
          $n += 2;
          ($x,$y) = ($y, -$x);   # rotate -90
          $slope = 0;
        } else {
          $y -= $horiz;
          ### slope not digit 1, move to: "$x,$y"

          if ($x < $horiz) {
            ### digit 2 ...
            $n += 3;
            $slope = 0;

          } else {
            ### digit 3 ...
            $n += 4;
            $x -= $horiz;
          }
        }
      }

    } else {
      if ($x < $horiz) {
        ### digit 0 ...
        $n += 1;
      } else {
        $x -= $horiz;
        ### not digit 0, move to: "$x,$y"

        if ($x < $diag) {
          ### digit 1 ...
          $n += 2;
          $slope = 1;
        } else {
          $x -= $diag;
          ### not digit 1, move to: "$x,$y"

          if ($x < $diag) {
            ### digit 2 ...
            $n += 3;
            $slope = 1;
            ($x,$y) = ($diag-$y, $x);   # offset and rotate +90

          } else {
            ### digit 3 ...
            $n += 4;
            $x -= $diag;
          }
        }
      }
    }
  }
  ### final: "$x,$y n=$n"

  if ($x == 0 && $y == 0) {
    return $n;
  } else {
    return undef;
  }
}

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

  foreach ($x1,$y1, $x2,$y2) {
    if (is_infinite($_)) {
      return (0, $_);
    }
    $_ = abs(round_nearest($_));
  }
  if ($x1 > $x2) { ($x1,$x2) = ($x2,$x1); }
  if ($y1 > $y2) { ($y1,$y2) = ($y2,$y1); }
  my $max = ($x2 > $y2 ? $x2 : $y2);

  # Nend = 4 * [ 1 + ... + 4^level ]
  #      = 4 + 16 + ... + 4^(level+1)
  #
  my $horiz = 4;
  my $diag = 3;
  my $nhi = 4;
  for (;;) {
    $nhi += 1;
    $nhi *= 4;
    my $nextdiag = $horiz + 2*$diag;
    if (($self->{'inward'} ? $horiz : $nextdiag) >= 2*$max) {
      return (1, $nhi);
    }
    $horiz = $nextdiag + $horiz;   # 2*$horiz + 2*$diag;
    $diag = $nextdiag;
  }
}


#------------------------------------------------------------------------------
# Nstart = (4^(k+1) - 1)/3
# Nend = Nstart(k+1) - 1
#      = (4*4^(k+1) - 1)/3 - 1
#      = (4*4^(k+1) - 1 - 3)/3
#      = (4*4^(k+1) - 4)/3
#      = 4*(4^(k+1) - 1)/3
#      = 4*Nstart(k)

sub level_to_n_range {
  my ($self, $level) = @_;
  my $n_lo = (4**($level+1) - 1)/3;
  return ($n_lo, 4*$n_lo);
}
sub n_to_level {
  my ($self, $n) = @_;
  if ($n < 1) { return undef; }
  if (is_infinite($n)) { return $n; }
  my ($pow,$exp) = round_down_pow (3*$n + 1, 4);
  return $exp-1;
}

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

    #                             15                     3
    #                            /  \
    #                     17--16      14--13             2
    #                      |               |
    #                     18              12             1
    #                   /       4 -- 3      \
    #                 19             |        11     <- Y=0
    #                   \       1 -- 2      /
    #                     20              10            -1
    #                                      |
    #                      5-- 6       8-- 9            -2
    #                            \   /
    #                              7                    -3
    #
    #                                                   -4
    #
    #                                                   -5
    #
    # ...                                               -6
    #
    # 21--22      24--25                      33--...   -7
    #       \   /       \                   /
    #         23          26              32            -8
    #                      |               |
    #                     27--28      30--31            -9
    #                           \   /
    #                             29                   -10




=for stopwords eg Ryde ie Math-PlanePath Koch Nstart Xstart,Ystart OEIS Xstart

=head1 NAME

Math::PlanePath::KochSquareflakes -- four-sided Koch snowflakes

=head1 SYNOPSIS

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

=head1 DESCRIPTION

This is the Koch curve shape arranged as four-sided concentric snowflakes.

=cut

# math-image --path=KochSquareflakes --all --output=numbers_dash --size=132x50

=pod

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

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

The innermost square N=1 to N=4 is the initial figure.  Its sides expand as
the Koch curve pattern in subsequent rings.  The initial figure is on
X=+/-0.5,Y=+/-0.5 fractions.  The points after that are integer X,Y.

=head1 Inward

The C<inward=E<gt>1> option can direct the sides inward.  The shape and
lengths etc are the same.  The angles and sizes mean there's no overlaps.

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

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

=head2 Level Ranges

Counting the innermost N=1 to N=4 square as level 0, a given level has

    looplen = 4*4^level

many points.  The start of a level is N=1 plus the preceding loop lengths so

    Nstart = 1 + 4*[ 1 + 4 + 4^2 + ... + 4^(level-1) ]
           = 1 + 4*(4^level - 1)/3
           = (4^(level+1) - 1)/3

and the end of a level similarly the total loop lengths, or simply one less
than the next Nstart,

    Nend = 4 * [ 1 + ... + 4^level ]
         = (4^(level+2) - 4) / 3

         = Nstart(level+1) - 1

For example,

    level  Nstart   Nend                       (A002450,A080674)
      0       1       4
      1       5      20
      2      21      84
      3      85     340                       

X<Lucas Sequence>The Xstart,Ystart position of the Nstart corner is a Lucas
sequence,

    Xstart(0) = -0.5
    Xstart(1) = -2
    Xstart(2) = 4*Xstart(1) - 2*Xstart(0) = -7
    Xstart(3) = 4*Xstart(2) - 2*Xstart(1) = -24
    ...
    Xstart(level+1) = 4*Xstart(level) - 2*Xstart(level-1)

    0.5, 2, 7, 24, 82, 280, 956, 3264, ...             (A003480)

This recurrence occurs because the replications are 4 wide when horizontal
but 3 wide when diagonal.

=head1 FUNCTIONS

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

=over 4

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

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

Create and return a new path object.

=back

=head2 Level Methods

=over

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

Return per L</Level Ranges> above,

    (  (4**$level - 1)/3,
     4*(4**$level - 1)/3 )

=back

=head1 OEIS

Entries in Sloane's Online Encyclopedia of Integer Sequences related to this
path include

=over

L<http://oeis.org/A003480> (etc)

=back

    A003480    -X and -Y coordinate first point of each ring
               likewise A020727
    A007052    X,Y coordinate of axis crossing,
               and also maximum height of a side
    A072261    N on Y negative axis (half way along first side)
    A206374    N on South-East diagonal (end of first side)

=head1 SEE ALSO

L<Math::PlanePath>,
L<Math::PlanePath::KochSnowflakes>

=head1 HOME PAGE

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

=head1 LICENSE

Copyright 2011, 2012, 2013, 2014, 2015, 2016 Kevin Ryde

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
