/usr/share/perl5/Math/PlanePath/GosperSide.pm

# Copyright 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/>.


# arms begin at 0,0 or at 1 in ?


# math-image --path=GosperSide --lines --scale=10
# math-image --path=GosperSide --output=numbers


package Math::PlanePath::GosperSide;
use 5.004;
use strict;
use List::Util 'min','max';
use POSIX 'ceil';
use Math::PlanePath::GosperIslands;
use Math::PlanePath::SacksSpiral;

use vars '$VERSION', '@ISA', '@_xend','@_yend';
$VERSION = 117;
use Math::PlanePath;
@ISA = ('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 Devel::Comments;

use constant n_start => 0;

# secret experimental as yet ...
#
# use constant parameter_info_array => [ { name      => 'arms',
#                                          share_key => 'arms_6',
#                                          type      => 'integer',
#                                          minimum   => 1,
#                                          maximum   => 6,
#                                          default   => 1,
#                                          width     => 1,
#                                          description => 'Arms',
#                                        } ];

use constant x_negative_at_n => 113;
use constant y_negative_at_n => 11357;

use constant dx_minimum => -2;
use constant dx_maximum => 2;
use constant dy_minimum => -1;
use constant dy_maximum => 1;

*_UNDOCUMENTED__dxdy_list = \&Math::PlanePath::_UNDOCUMENTED__dxdy_list_six;
# 2,0,   # E  N=0
# 1,1,   # NE  N=1
# -1,1,   # NW  N=4
# -2,0,   # W  N=13
# -1,-1,   # SW  N=40
# 1,-1,   # SE  N=121
use constant _UNDOCUMENTED__dxdy_list_at_n => 121;

use constant absdx_minimum => 1;
use constant dsumxy_minimum => -2; # diagonals
use constant dsumxy_maximum => 2;
use constant ddiffxy_minimum => -2;
use constant ddiffxy_maximum => 2;
use constant dir_maximum_dxdy => (1,-1); # South-East


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

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

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

  my $x;
  my $y = my $yend = ($n * 0);  # inherit bignum 0
  my $xend = $y + 2;            # inherit bignum 2
  {
    my $int = int($n);
    $x = 2 * ($n - $int);
    $n = $int;
  }


  if ((my $arms = $self->{'arms'}) > 1) {
    my $rot = _divrem_mutate ($n, $arms);
    if ($rot >= 3) {
      $rot -= 3;
      $x = -$x;    # rotate 180, knowing y=0,yend=0
      $xend = -2;
    }
    if ($rot == 1) {
      $x = $y = $x/2;   # rotate +60, knowing y=0,yend=0
      $xend = $yend = $xend/2;
    } elsif ($rot == 2) {
      $y = $x/2;   # rotate +120, knowing y=0,yend=0
      $x = -$y;
      $yend = $xend/2;
      $xend = -$yend;
    }
  }

  foreach my $digit (digit_split_lowtohigh($n,3)) {
    my $xend_offset = 3*($xend-$yend)/2;   # end and end +60
    my $yend_offset = ($xend+3*$yend)/2;

    ### at: "$x,$y"
    ### $digit
    ### $xend
    ### $yend
    ### $xend_offset
    ### $yend_offset

    if ($digit == 1) {
      ($x,$y) = (($x-3*$y)/2  + $xend,   # rotate +60
                 ($x+$y)/2    + $yend);
    } elsif ($digit == 2) {
      $x += $xend_offset;   # offset and offset +60
      $y += $yend_offset;
    }
    $xend += $xend_offset;   # offset and offset +60
    $yend += $yend_offset;
  }

  ### final: "$x,$y"
  return ($x, $y);
}

# level = (log(hypot) + log(2*.99)) * 1/log(sqrt(7))
#       = (log(hypot^2)/2 + log(2*.99)) * 1/log(sqrt(7))
#       = (log(hypot^2) + 2*log(2*.99)) * 1/(2*log(sqrt(7)))
#
sub xy_to_n {
  my ($self, $x, $y) = @_;
  $x = round_nearest ($x);
  $y = round_nearest ($y);
  ### GosperSide xy_to_n(): "$x, $y"

  if (($x ^ $y) & 1) {
    return undef;
  }

  my $h2 = $x*$x + $y*$y*3 + 1;
  my $level = max (0,
                   ceil ((log($h2) + 2*log(2*.99)) * (1/2*log(sqrt(7)))));
  if (is_infinite($level)) {
    return $level;
  }
  return Math::PlanePath::GosperIslands::_xy_to_n_in_level($x,$y,$level);
}


# Points beyond N=3^level only go a small distance back before that N
# hypotenuse.
#     hypot = .99 * 2 * sqrt(7)^level
#     sqrt(7)^level = hypot / (2*.99)
#     sqrt(7)^level = hypot / (2*.99)
#     level = log(hypot / (2*.99)) / log(sqrt(7))
#           = (log(hypot) + log(2*.99)) * 1/log(sqrt(7))
#
# not exact
sub rect_to_n_range {
  my ($self, $x1,$y1, $x2,$y2) = @_;
  $y1 *= sqrt(3);
  $y2 *= sqrt(3);
  my ($r_lo, $r_hi) = Math::PlanePath::SacksSpiral::_rect_to_radius_range
    ($x1,$y1, $x2,$y2);
  my $level = max (0,
                   ceil ((log($r_hi+.1) + log(2*.99)) * (1/log(sqrt(7)))));
  return (0,
          $self->{'arms'} * 3 ** $level - 1);
}

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

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

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

=for stopwords eg Ryde Math-PlanePath Gosper

=head1 NAME

Math::PlanePath::GosperSide -- one side of the Gosper island

=head1 SYNOPSIS

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

=head1 DESCRIPTION

X<Gosper, William>This path is a single side of the Gosper island, in
integers (L<Math::PlanePath/Triangular Lattice>).

                                        20-...        14
                                       /
                               18----19               13
                              /
                            17                        12
                              \
                               16                     11
                              /
                            15                        10
                              \
                               14----13                9
                                       \
                                        12             8
                                       /
                                     11                7
                                       \
                                        10             6
                                       /
                                8---- 9                5
                              /
                       6---- 7                         4
                     /
                    5                                  3
                     \
                       4                               2
                     /
              2---- 3                                  1
            /
     0---- 1                                       <- Y=0

     ^
    X=0 1  2  3  4  5  6  7  8  9 10 11 12 13 ...

The path slowly spirals around counter clockwise, with a lot of wiggling in
between.  The N=3^level point is at

   N = 3^level
   angle = level * atan(sqrt(3)/5)
         = level * 19.106 degrees
   radius = sqrt(7) ^ level

A full revolution for example takes roughly level=19 which is about
N=1,162,000,000.

Both ends of such levels are in fact sub-spirals, like an "S" shape.

The path is both the sides and the radial spokes of the C<GosperIslands>
path, as described in L<Math::PlanePath::GosperIslands/Side and Radial
Lines>.  Each N=3^level point is the start of a C<GosperIslands> ring.

The path is the same as the C<TerdragonCurve> except the turns here are by
60 degrees each, whereas C<TerdragonCurve> is by 120 degrees.  See
L<Math::PlanePath::TerdragonCurve> for the turn sequence and total direction
formulas etc.

=head1 FUNCTIONS

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

=over 4

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

Create and return a new path object.

=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 C<$n> gives a point on the straight line between integer N.

=back

=head2 Level Methods

=over

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

Return C<(0, 3**$level)>.

=back

=head1 FORMULAS

=head2 Level Endpoint

The endpoint of each level N=3^k is at

    X + Y*i*sqrt(3) = b^k
    where b = 2 + w = 5/2 + sqrt(3)/2*i
          where w=1/2 + sqrt(3)/2*i sixth root of unity

    X(k) = ( 5*X(k-1) - 3*Y(k-1) )/2        for k>=1
    Y(k) = (   X(k-1) + 5*Y(k-1) )/2
           starting X(0)=2 Y(0)=0

    X(k) = 5*X(k-1) - 7*X(k-2)        for k>=2
           starting X(0)=2 X(1)=5
         = 2, 5, 11, 20, 23, -25, -286, -1255, -4273, -12580, -32989,..

    Y(k) = 5*Y(k-1) - Y*X(k-2)        for k>=2
           starting Y(0)=0 Y(1)=1
         = 0, 1,  5, 18, 55, 149,  360,   757,  1265, 1026, -3725, ...
                                                            (A099450)

=for Test-Pari-DEFINE  X(k) = if(k==0,2, k==1,5, 5*X(k-1)-7*X(k-2))

=for Test-Pari-DEFINE  Y(k) = if(k==0,0, k==1,1, 5*Y(k-1)-7*Y(k-2))

=for Test-Pari-DEFINE  X_samples = [ 2, 5, 11, 20, 23, -25, -286, -1255, -4273, -12580, -32989 ]

=for Test-Pari-DEFINE  Y_samples = [ 0, 1,  5, 18, 55, 149,  360,   757,  1265, 1026, -3725 ]

=for Test-Pari  vector(length(X_samples),k,my(k=k-1); X(k)) == X_samples

=for Test-Pari  vector(length(Y_samples),k,my(k=k-1); Y(k)) == Y_samples

=for Test-Pari  vector(20,k, X(k)) == vector(20,k, (5*X(k-1)-3*Y(k-1))/2) /* k>=1 */

=for Test-Pari  vector(20,k, Y(k)) == vector(20,k, (X(k-1)+5*Y(k-1))/2)   /* k>=1 */

The curve base figure is XY(k)=XY(k-1)+rot60(XY(k-1))+XY(k-1) giving XY(k) =
(2+w)^k = b^k where w is the sixth root of unity giving the rotation by +60
degrees.

The mutual recurrences are similar with the rotation done by (X-3Y)/2,
(Y+X)/2 per L<Math::PlanePath/Triangular Lattice>.  The separate recurrences
are found by using the first to get Y(k-1) = -2/3*X(k) + 5/3*X(k-1) and
substitute into the other to get X(k+1).  Similar the other way around for
Y(k+1).

=cut

# w=1/2 + sqrt(3)/2*I; b = 2+w
# nearly(x) = if(abs(x-round(x))<1e-10,round(x),x)
# Y(k) = nearly(2*imag(b^k/sqrt(3)))
# vector(20,k,my(k=k-1); Y(k))
# 0, 1, 5, 18, 55, 149, 360, 757, 1265, 1026, -3725, -25807, -102960, -334151,
# A099450 Expansion of 1/(1-5x+7x^2).

# X(k) = nearly(2*real(b^k))
# vector(20,k,my(k=k-1); X(k))
# 2, 5, 11, 20, 23, -25, -286, -1255, -4273, -12580, -32989, -76885, -153502,

#       *---*  2+w
#      /
# *---*
# X+I*Y = a+b*w
#       = a+b/2 + sqrt(3)/2*b*i
# X = a+b/2  Y=sqrt(3)/2*b
# b=2*Y/sqrt(3)
# a = X - b = X - 2*Y/sqrt(3)
#
# a(k) = X(k) - 2*Y(k)
# vector(20,k,my(k=k-1); a(k))
# 2, 3, 1, -16, -87, -323, -1006, -2769, -6803, -14632, -25539, -25271, 52418,

# XY(k) = 2*XY(k-1) + rot60 XY(k-1)
# X(k) = 2*X(k-1) + (X(k-1) - 3*Y(k-1))/2
# Y(k) = 2*Y(k-1) + (X(k-1) +   Y(k-1))/2
#
# X(k) = (5*X(k-1) - 3*Y(k-1))/2
# Y(k) = (  X(k-1) + 5*Y(k-1))/2
# 2 5 11
# 0 1  5
#
# sep
# Y(k-1) = -2/3*X(k) + 5/3*X(k-1)
# X(k-1) =   2 *Y(k) -  5 *Y(k-1)
#
# subst Y
# Y(k) = (  X(k-1) + 5*Y(k-1))/2
# -2/3*X(k+1) + 5/3*X(k) = (  X(k-1) + 5*(-2/3*X(k) + 5/3*X(k-1)))/2
# -2/3*X(k+1) + 5/3*X(k) = 1/2*X(k-1) + -5/2*2/3*X(k) + 5/2*5/3*X(k-1)
# -2/3*X(k+1) + 10/3*X(k) =  14/3*X(k-1)
# X(k+1) = 5*X(k) - 7*X(k-1)
#
# subst X
# X(k) = (5*X(k-1) - 3*Y(k-1))/2
# 2 *Y(k+1) -  5 *Y(k) = (5*(2 *Y(k) -  5 *Y(k-1)) - 3*Y(k-1))/2
# 2 *Y(k+1) -  5 *Y(k) = 5*2/2 *Y(k) -  5*5/2 *Y(k-1) - 3/2*Y(k-1)
# 2 *Y(k+1)  = 10*Y(k) - 14*Y(k-1)
# Y(k+1)  = 5*Y(k) - 7*Y(k-1)
# 5*5-7=18

=pod

=head1 OEIS

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

=over

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

=back

    A099450   Y at N=3^k (for k>=1)

The turn sequence is the same as the terdragon curve, see
L<Math::PlanePath::TerdragonCurve/OEIS> for the numerous turn forms, N
positions of turns, etc.

=head1 SEE ALSO

L<Math::PlanePath>,
L<Math::PlanePath::GosperIslands>,
L<Math::PlanePath::TerdragonCurve>,
L<Math::PlanePath::KochCurve>

L<Math::Fractal::Curve>

=head1 HOME PAGE

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

=head1 LICENSE

Copyright 2011, 2012, 2013, 2014 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
libmath-planepath-perl 117-1 / usr / share / perl5 / Math / PlanePath / GosperSide.pm
