# Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020 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 . # math-image --path=GosperReplicate --lines --scale=10 # math-image --path=GosperReplicate --all --output=numbers_dash # math-image --path=GosperReplicate,numbering_type=rotate --all --output=numbers_dash # package Math::PlanePath::GosperReplicate; use 5.004; use strict; use List::Util qw(max); use POSIX 'ceil'; use Math::Libm 'hypot'; use Math::PlanePath::SacksSpiral; use vars '$VERSION', '@ISA'; $VERSION = 129; use Math::PlanePath; @ISA = ('Math::PlanePath'); use Math::PlanePath::Base::Generic 'is_infinite', 'round_nearest'; use Math::PlanePath::Base::Digits 'round_up_pow', 'digit_split_lowtohigh', 'digit_join_lowtohigh'; # uncomment this to run the ### lines #use Smart::Comments; use constant parameter_info_array => [ { name => 'numbering_type', display => 'Numbering', share_key => 'numbering_type_rotate', type => 'enum', default => 'fixed', choices => ['fixed','rotate'], choices_display => ['Fixed','Rotate'], description => 'Fixed or rotating sub-part numbering.', }, ]; use constant n_start => 0; *xy_is_visited = \&Math::PlanePath::Base::Generic::xy_is_even; use constant x_negative_at_n => 3; use constant y_negative_at_n => 5; use constant absdx_minimum => 1; use constant dir_maximum_dxdy => (3,-1); #------------------------------------------------------------------------------ sub new { my $self = shift->SUPER::new (@_); $self->{'numbering_type'} ||= 'fixed'; # default return $self; } sub _digits_rotate_lowtohigh { my ($aref) = @_; my $rot = 0; foreach my $digit (reverse @$aref) { if ($digit) { $rot += $digit-1; $digit = ($rot % 6) + 1; # mutate $aref } } } sub _digits_unrotate_lowtohigh { my ($aref) = @_; my $rot = 0; foreach my $digit (reverse @$aref) { if ($digit) { $digit = ($digit-1-$rot) % 6; # mutate $aref $rot += $digit; $digit++; } } } sub n_to_xy { my ($self, $n) = @_; ### GosperReplicate n_to_xy(): $n if ($n < 0) { return; } if (is_infinite($n)) { return ($n,$n); } { my $int = int($n); ### $int ### $n if ($n != $int) { my ($x1,$y1) = $self->n_to_xy($int); my ($x2,$y2) = $self->n_to_xy($int+1); my $frac = $n - $int; # inherit possible BigFloat my $dx = $x2-$x1; my $dy = $y2-$y1; return ($frac*$dx + $x1, $frac*$dy + $y1); } $n = $int; # BigFloat int() gives BigInt, use that } my $x = my $y = $n*0; # inherit bigint from $n my $sx = $x + 2; # 2 my $sy = $x; # 0 # digit # 3 2 # \ / # 4---0---1 # / \ # 5 6 my @digits = digit_split_lowtohigh($n,7); if ($self->{'numbering_type'} eq 'rotate') { _digits_rotate_lowtohigh(\@digits); } foreach my $digit (@digits) { ### digit: "$digit $x,$y side $sx,$sy" if ($digit == 1) { ### right ... # $x = -$x; # rotate 180 # $y = -$y; $x += $sx; $y += $sy; } elsif ($digit == 2) { ### up right ... # ($x,$y) = ((3*$y-$x)/2, # rotate -120 # ($x+$y)/-2); $x += ($sx - 3*$sy)/2; # at +60 $y += ($sx + $sy)/2; } elsif ($digit == 3) { ### up left ... # ($x,$y) = (($x+3*$y)/2, # -60 # ($y-$x)/2); $x += ($sx + 3*$sy)/-2; # at +120 $y += ($sx - $sy)/2; } elsif ($digit == 4) { ### left $x -= $sx; # at -180 $y -= $sy; } elsif ($digit == 5) { ### down left # ($x,$y) = (($x-3*$y)/2, # rotate +60 # ($x+$y)/2); $x += (3*$sy - $sx)/2; # at -120 $y += ($sx + $sy)/-2; } elsif ($digit == 6) { ### down right # ($x,$y) = (($x+3*$y)/-2, # rotate +120 # ($x-$y)/2); $x += ($sx + 3*$sy)/2; # at -60 $y += ($sy - $sx)/2; } # 2*(sx,sy) + rot+60(sx,sy) ($sx,$sy) = ((5*$sx - 3*$sy) / 2, ($sx + 5*$sy) / 2); } return ($x,$y); } # modulus # 1 3 # \ / # 5---0---2 # / \ # 4 6 # 0 1 2 3 4 5 6 my @modulus_to_x = (0,-1, 2, 1,-1,-2, 1); my @modulus_to_y = (0, 1, 0, 1,-1, 0,-1); my @modulus_to_digit = (0, 3, 1, 2, 5, 4, 6); sub xy_to_n { my ($self, $x, $y) = @_; ### GosperReplicate xy_to_n(): "$x, $y" $x = round_nearest($x); $y = round_nearest($y); if (($x + $y) % 2) { return undef; } my $level = _xy_to_level_ceil($x,$y); if (is_infinite($level)) { return $level; } my $zero = ($x * 0 * $y); # inherit bignum 0 my @n; # digits low to high while ($level-- >= 0 && ($x || $y)) { ### at: "$x,$y m=".(($x + 2*$y) % 7) my $m = ($x + 2*$y) % 7; push @n, $modulus_to_digit[$m]; $x -= $modulus_to_x[$m]; $y -= $modulus_to_y[$m]; ### digit: "to $x,$y" ### assert: (3 * $y + 5 * $x) % 14 == 0 ### assert: (5 * $y - $x) % 14 == 0 # shrink ($x,$y) = ((3*$y + 5*$x) / 14, (5*$y - $x) / 14); } if ($self->{'numbering_type'} eq 'rotate') { _digits_unrotate_lowtohigh(\@n); } return digit_join_lowtohigh (\@n, 7, $zero); } # 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); $r_hi *= 2; my $level_plus_1 = ceil( log(max(1,$r_hi/4)) / log(sqrt(7)) ) + 2; return (0, 7**$level_plus_1 - 1); } sub _xy_to_level_ceil { my ($x,$y) = @_; my $r = hypot($x,$y); $r *= 2; return ceil( log(max(1,$r/4)) / log(sqrt(7)) ) + 1; } #------------------------------------------------------------------------------ # levels sub level_to_n_range { my ($self, $level) = @_; return (0, 7**$level - 1); } sub n_to_level { my ($self, $n) = @_; if ($n < 0) { return undef; } if (is_infinite($n)) { return $n; } $n = round_nearest($n); my ($pow, $exp) = round_up_pow ($n+1, 7); return $exp; } #------------------------------------------------------------------------------ 1; __END__ =for stopwords PlanePath eg Ryde Gosper Math-PlanePath =head1 NAME Math::PlanePath::GosperReplicate -- self-similar hexagon replications =head1 SYNOPSIS use Math::PlanePath::GosperReplicate; my $path = Math::PlanePath::GosperReplicate->new; my ($x, $y) = $path->n_to_xy (123); =head1 DESCRIPTION This is a self-similar hexagonal tiling of the plane. At each level the shape is the Gosper island. 17----16 4 / \ 24----23 18 14----15 3 / \ \ 25 21----22 19----20 10---- 9 2 \ / \ 26----27 3---- 2 11 7---- 8 1 / \ \ 31----30 4 0---- 1 12----13 <- Y=0 / \ \ 32 28----29 5---- 6 45----44 -1 \ / \ 33----34 38----37 46 42----43 -2 / \ \ 39 35----36 47----48 -3 \ 40----41 -4 ^ -7 -6 -5 -4 -3 -2 -1 X=0 1 2 3 4 5 6 7 Points are spread out on every second X coordinate to make a triangular lattice in integer coordinates (see L). The base pattern is the inner N=0 to N=6, then six copies of that shape are arranged around as the blocks N=7,14,21,28,35,42. Then six copies of the resulting N=0 to N=48 shape are replicated around, etc. Each point can be taken as a little hexagon, so that all points tile the plane with hexagons. The innermost N=0 to N=6 are for instance, * * / \ / \ / \ / \ * * * | 3 | 2 | * * * / \ / \ / \ / \ / \ / \ * * * * | 4 | 0 | 1 | * * * * \ / \ / \ / \ / \ / \ / * * * | 5 | 6 | * * * \ / \ / \ / \ / * * The further replications are the same arrangement, but the sides become ever wigglier and the centres rotate around. The rotation can be seen N=7 at X=5,Y=1 which is up from the X axis. The C path is this same replicating shape, but starting from a side instead of the middle and traversing in such as way as to make each N adjacent. The C curve itself is this replication too, but segments across hexagons. =head2 Complex Base The path corresponds to expressing complex integers X+i*Y in a base b = 5/2 + i*sqrt(3)/2 =cut # GP-DEFINE sqrt3 = quadgen(12); # GP-DEFINE sqrt3i = quadgen(-12); # GP-Test sqrt3^2 == 3 # GP-Test sqrt3i^2 == -3 # GP-DEFINE b = 5/2 + sqrt3i/2; =pod with coordinates scaled to put equilateral triangles on a square grid. So for integer X,Y on the triangular grid (X,Y either both odd or both even), X/2 + i*Y*sqrt(3)/2 = a[n]*b^n + ... + a[2]*b^2 + a[1]*b + a[0] where each digit a[i] is either 0 or a sixth root of unity encoded into base-7 digits of N, w6 = e^(i*pi/3) sixth root of unity, b = 2 + w6 = 1/2 + i*sqrt(3)/2 N digit a[i] complex number ------- ------------------- 0 0 1 w6^0 = 1 2 w6^1 = 1/2 + i*sqrt(3)/2 3 w6^2 = -1/2 + i*sqrt(3)/2 4 w6^3 = -1 5 w6^4 = -1/2 - i*sqrt(3)/2 6 w6^5 = 1/2 - i*sqrt(3)/2 =cut # GP-DEFINE w6 = 1/2 + sqrt3i/2; # GP-Test w6^6 == 1 # GP-Test w6^0 == 1 # GP-Test w6^1 == 1/2 + sqrt3i/2 # GP-Test w6^2 == -1/2 + sqrt3i/2 # GP-Test w6^3 == -1 # GP-Test w6^4 == -1/2 - sqrt3i/2 # GP-Test w6^5 == 1/2 - sqrt3i/2 # GP-Test (5/2)^2 + (sqrt3/2)^2 == 7 # GP-DEFINE z_digit(d) = [0, 1,w6,w6^2, -1,w6^4,w6^5][d+1]; # GP-DEFINE z_point(n) = \ # GP-DEFINE subst(Pol(apply(z_digit,digits(n,7))),'x,b); # GP-Test z_point(0) == 0 # GP-Test z_point(1) == 1 # GP-Test z_point(2) == w6 # GP-Test z_point(7) == w6+2 # GP-DEFINE nearly_equal_epsilon = 1e-15; # GP-DEFINE nearly_equal(x,y, epsilon=nearly_equal_epsilon) = \ # GP-DEFINE abs(x-y) < epsilon; # GP-DEFINE to_base7(n) = fromdigits(digits(n,7)); # GP-DEFINE from_base7(n) = fromdigits(digits(n),7); =pod 7 digits suffice because norm(b) = (5/2)^2 + (sqrt(3)/2)^2 = 7 =cut # GP-Test norm(b) == 7 # GP-Test (5/2)^2 + (sqrt3/2)^2 == 7 =pod =head2 Rotate Numbering Parameter C 'rotate'> applies a rotation in each sub-part according to its location around the preceding level. The effect can be illustrated by writing N in base-7. Part 10-16 is the same as the middle 0-6. Part 20-26 has a rotation by +60 degrees. Part 30-36 has rotation by +120 degrees, and so on. =cut # start from this, then mangled by hand # math-image --path=GosperReplicate,numbering_type=rotate --all --output=numbers_dash =pod 22----21 / / numbering_type => 'rotate' 31 36 23 20 26 N shown in base-7 / \ \ \ / 32 30 35 24----25 13----12 \ / / \ 33----34 3---- 2 14 10----11 / \ \ 46----45 4 0---- 1 15----16 \ \ 41----40 44 5---- 6 64----63 \ / / \ 42----43 55----54 65 60 62 / \ \ \ / 56 50 53 66 61 / / 51----52 Notice this means in each part the 11, 21, 31, etc, points are directed away from the middle in the same way, relative to the sub-part locations. Working through the expansions gives the following rule for when an N is on the boundary of level k, write N in k many base-7 digits (empty string if k=0) if any 0 digit then non-boundary ignore high digit and all 1 digits if any 4 or 5 digit then non-boundary if any 32, 33, 66 pair then non-boundary A 0 digit is the middle of a block, or 4 or 5 digit the inner side of a block, for kE=1, hence non-boundary. After that the 6,1,2,3 parts variously expand with rotations so that a 66 is enclosed on the clockwise side and 32 and 33 on the anti-clockwise side. =cut # in decimal # 16----15 # / / # 22 27 17 14 20 # / \ \ \ / # 23 21 26 18----19 10---- 9 # \ / / \ # 24----25 3---- 2 11 7---- 8 # / \ \ # 34----33 4 0---- 1 12----13 # \ \ # 29----28 32 5---- 6 46----45 # \ / / \ # 30----31 40----39 47 42 44 # / \ \ \ / # 41 35 38 48 43 # / / # 36----37 =pod =head1 FUNCTIONS See L for behaviour common to all path classes. =over 4 =item C<$path = Math::PlanePath::GosperReplicate-Enew ()> =item C<$path = Math::PlanePath::GosperReplicate-Enew (numbering_type =E $str)> Create and return a new path object. The C parameter can be "fixed" (default) "rotate" =item C<($x,$y) = $path-En_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 0> then the return is an empty list. =back =head2 Level Methods =over =item C<($n_lo, $n_hi) = $path-Elevel_to_n_range($level)> Return C<(0, 7**$level - 1)>. =back =head1 FORMULAS =head2 Axis Rotations In the fixed numbering, digit positions 1,2,3,4,5,6 go around +60deg each, so the N for rotation of X,Y by +60 degrees is each digit +1. N = 0, 1, 2, 3, 4, 5, 6, 10, 11, 12 rot+60(N) = 0, 2, 3, 4, 5, 6, 1, 14, 16, 17, ... decimal = 0, 2, 3, 4, 5, 6, 1, 20, 22, 23, ... base7 rot+120(N) = 0, 3, 4, 5, 6, 1, 2, 21, 24, 25, ... decimal = 0, 3, 4, 5, 6, 1, 2, 30, 33, 34, ... base7 etc =cut # rot180(N) = 0, 4, 5, 6, 1, 2, 3, 28, 32, 33, ... decimal # = 0, 4, 5, 6, 1, 2, 3, 40, 44, 45, ... base7 # # rot-120(N) = 0, 5, 6, 1, 2, 3, 4, 35, 40, 41, ... decimal # = 0, 5, 6, 1, 2, 3, 4, 50, 55, 56, ... base7 # # rot-60(N) = 0, 6, 1, 2, 3, 4, 5, 42, 48, 43, ... decimal # = 0, 6, 1, 2, 3, 4, 5, 60, 66, 61, ... base7 # GP-DEFINE digit_plus1(d) = [0,2,3,4,5,6,1][d+1]; # GP-DEFINE digit_plus2(d) = [0,3,4,5,6,1,2][d+1]; # GP-DEFINE digit_plus3(d) = [0,4,5,6,1,2,3][d+1]; # GP-DEFINE digit_minus2(d) = [0,5,6,1,2,3,4][d+1]; # GP-DEFINE digit_minus1(d) = [0,6,1,2,3,4,5][d+1]; # GP-DEFINE N_rotate_plus60(n) = fromdigits(apply(digit_plus1, digits(n,7)),7); # GP-DEFINE N_rotate_plus120(n)= fromdigits(apply(digit_plus2, digits(n,7)),7); # GP-DEFINE N_rotate_180(n) = fromdigits(apply(digit_plus3, digits(n,7)),7); # GP-DEFINE N_rotate_minus120(n)=fromdigits(apply(digit_minus2,digits(n,7)),7); # GP-DEFINE N_rotate_minus60(n)= fromdigits(apply(digit_minus1,digits(n,7)),7); # GP-Test my(v=[0, 2, 3, 4, 5, 6, 1, 14, 16, 17]); /* samples shown */ \ # GP-Test vector(#v,n,n--; N_rotate_plus60(n)) == v # GP-Test my(v=[0, 2, 3, 4, 5, 6, 1, 20, 22, 23]); /* samples shown */ \ # GP-Test vector(#v,n,n--; to_base7(N_rotate_plus60(n))) == v # GP-Test my(v=[0, 3, 4, 5, 6, 1, 2, 21, 24, 25]); /* samples shown */ \ # GP-Test vector(#v,n,n--; N_rotate_plus120(n)) == v # GP-Test my(v=[0, 3, 4, 5, 6, 1, 2, 30, 33, 34]); /* samples shown */ \ # GP-Test vector(#v,n,n--; to_base7(N_rotate_plus120(n))) == v # GP-Test my(v=[0, 4, 5, 6, 1, 2, 3, 28, 32, 33]); /* samples shown */ \ # GP-Test vector(#v,n,n--; N_rotate_180(n)) == v # GP-Test my(v=[0, 4, 5, 6, 1, 2, 3, 40, 44, 45]); /* samples shown */ \ # GP-Test vector(#v,n,n--; to_base7(N_rotate_180(n))) == v # GP-Test my(v=[0, 5, 6, 1, 2, 3, 4, 35, 40, 41]); /* samples shown */ \ # GP-Test vector(#v,n,n--; N_rotate_minus120(n)) == v # GP-Test my(v=[0, 5, 6, 1, 2, 3, 4, 50, 55, 56]); /* samples shown */ \ # GP-Test vector(#v,n,n--; to_base7(N_rotate_minus120(n))) == v # GP-Test my(v=[0, 6, 1, 2, 3, 4, 5, 42, 48, 43]); /* samples shown */ \ # GP-Test vector(#v,n,n--; N_rotate_minus60(n)) == v # GP-Test my(v=[0, 6, 1, 2, 3, 4, 5, 60, 66, 61]); /* samples shown */ \ # GP-Test vector(#v,n,n--; to_base7(N_rotate_minus60(n))) == v # GP-Test vector(500,n,n--; z_point(N_rotate_plus60(n))) == \ # GP-Test vector(500,n,n--; w6*z_point(n)) # GP-Test vector(500,n,n--; z_point(N_rotate_plus120(n))) == \ # GP-Test vector(500,n,n--; w6^2*z_point(n)) # GP-Test vector(500,n,n--; z_point(N_rotate_180(n))) == \ # GP-Test vector(500,n,n--; -z_point(n)) # GP-Test vector(500,n,n--; z_point(N_rotate_minus120(n))) == \ # GP-Test vector(500,n,n--; conj(w6)^2*z_point(n)) # GP-Test vector(500,n,n--; z_point(N_rotate_minus60(n))) == \ # GP-Test vector(500,n,n--; conj(w6)*z_point(n)) # not in OEIS: 2, 3, 4, 5, 6, 1, 14, 16, 17 # not in OEIS: 2, 3, 4, 5, 6, 1, 20, 22, 23 # not in OEIS: 3, 4, 5, 6, 1, 2, 21, 24, 25 # not in OEIS: 3, 4, 5, 6, 1, 2, 30, 33, 34 # not in OEIS: 4, 5, 6, 1, 2, 3, 28, 32, 33 # not in OEIS: 4, 5, 6, 1, 2, 3, 40, 44, 45 # not in OEIS: 5, 6, 1, 2, 3, 4, 35, 40, 41 # not in OEIS: 5, 6, 1, 2, 3, 4, 50, 55, 56 # not in OEIS: 6, 1, 2, 3, 4, 5, 42, 48, 43 # not in OEIS: 6, 1, 2, 3, 4, 5, 60, 66, 61 =pod In the rotate numbering, just adding +1 (etc) at the high digit alone is rotation. =cut # GP-DEFINE n_rotate_highdigit(n,offset) = { # GP-DEFINE my(v=digits(n)); # GP-DEFINE v[1] = ((v[1]-1+offset)%6) + 1; # GP-DEFINE fromdigits(v,7); # GP-DEFINE } # for(offset=1,6,print(vector(18,n, n_rotate_highdigit(n,offset)))) # not in OEIS: 2, 3, 4, 5, 6, 1, 2, 3, 4, 14, 15, 16, 17, 18, 19, 20, 21, 22 # not in OEIS: 3, 4, 5, 6, 1, 2, 3, 4, 5, 21, 22, 23, 24, 25, 26, 27, 28, 29 # not in OEIS: 4, 5, 6, 1, 2, 3, 4, 5, 6, 28, 29, 30, 31, 32, 33, 34, 35, 36 # not in OEIS: 5, 6, 1, 2, 3, 4, 5, 6, 1, 35, 36, 37, 38, 39, 40, 41, 42, 43 # not in OEIS: 6, 1, 2, 3, 4, 5, 6, 1, 2, 42, 43, 44, 45, 46, 47, 48, 49, 50 # not in OEIS: 1, 2, 3, 4, 5, 6, 1, 2, 3, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16 =pod =head2 X,Y Extents The maximum X in a given level N=0 to 7^k-1 can be calculated from the replications. A given high digit 1 to 6 has sub-parts located at b^k*w6^(d-1). Those sub-parts are all the same, so the one with maximum real(b^k*w6^(d-1)) contains the maximum X. N_xmax_digit(j) = d=1to6 where real(w6^(d-1) * b^j) is maximum = 1,1,6,6,6,5,5,5,4,4,4,3,3,3,3,2,2, ... k-1 N_xmax(k) = digits N_xmax_digit(j) low digit j=0 j=0 = 0, 1, 8, 302, 2360, 16766, 100801, ... decimal = 0, 1, 11, 611, 6611, 66611, 566611, ... base7 k-1 z_xmax(k) = sum w6^d[j] * b^j j=0 each d[j] with real(w6^d[j] * b^j) maximum = 0, 1, 7/2+1/2*sqrt3*i, 10-sqrt3*i, 57/2-3/2*sqrt3*i,... xmax(k) = 2*real(z_xmax(k)) = 0, 2, 7, 20, 57, 151, 387, 1070, 2833, 7106, ... =cut # GP-DEFINE N_xmax_digit(j) = \ # GP-DEFINE my(p=b^j,d); vecmax(vector(6,d,real(w6^(d-1)*p)),&d); d; # GP-Test my(v=[1,1,6,6,6,5,5,5,4,4,4,3,3,3,3,2,2]); /* samples shown */ \ # GP-Test vector(#v,j,j--; N_xmax_digit(j)) == v # GP-DEFINE N_xmax(k) = fromdigits(Vecrev(vector(k,j,j--; N_xmax_digit(j))),7); # GP-Test my(v=[0, 1, 8, 302, 2360, 16766, 100801]); /* samples shown */ \ # GP-Test vector(#v,j,j--; N_xmax(j)) == v # GP-Test my(v=[0, 1, 11, 611, 6611, 66611, 566611]); /* samples shown */ \ # GP-Test vector(#v,j,j--; to_base7(N_xmax(j))) == v # GP-Test to_base7(N_xmax(51)) \ # GP-Test == 334445556661112222333444555666111222333344455566611 # GP-DEFINE z_xmax(k) = { # GP-DEFINE sum(j=0,k-1, # GP-DEFINE my(p=b^j, v=vector(6,d,(w6^(d-1)*p)), i); # GP-DEFINE vecmax(real(v),&i); # GP-DEFINE v[i]); # GP-DEFINE } # GP-Test my(v=[0, 1, 7/2+1/2*sqrt3i, 10-sqrt3i, 57/2-3/2*sqrt3i]); \ # GP-Test vector(#v,k,k--; z_xmax(k)) == v # GP-Test z_xmax(0) == 0 # GP-Test z_xmax(1) == 1 # GP-Test z_point(7) == 5/2 + 1/2*sqrt3i # GP-Test z_point(8) == 7/2 + 1/2*sqrt3i # GP-DEFINE xmax(k) = real(z_xmax(k)); # GP-Test my(v=[0, 2, 7, 20, 57, 151, 387, 1070, 2833, 7106]); \ # GP-Test vector(#v,k,k--; 2*xmax(k)) == v # GP-Test 2*xmax(45) == 12321054172600214702 # GP-Test 2*xmax(2) == 7 /* X of N=8 shown in sample numbers */ # vector(15,k,k--; N_xmax_digit(k)) # not in OEIS: 1, 1, 6, 6, 6, 5, 5, 5, 4, 4, 4, 3, 3, 3, 3 # vector(8,k,k++; N_xmax(k)) # vector(8,k,k++; to_base7(N_xmax(k))) # not in OEIS: 8, 57, 400, 10004, 77232, 547828, 3018457, 20312860 # not in OEIS: 11, 111, 1111, 41111, 441111, 4441111, 34441111, 334441111 # vector(6,k,k--; z_xmax(k)) # vector(8,k, norm(z_xmax(k))) # vector(10,k,k++; 2*real(z_xmax(k))) # vector(10,k,k++; 2*imag(z_xmax(k))) # vector(10,k,k++; real(z_xmax(k))+imag(z_xmax(k))) # not in OEIS: 1, 13, 103, 819, 5827, 39243, 291772, 2026399 \\ norm # not in OEIS: 7, 20, 57, 151, 387, 1070, 2833, 7106, 19686, 52675 \\ real # not in OEIS: 1, -2, -3, 13, -49, -86, 163, -1102, -2128, 1597 \\ imag # not in OEIS: 4, 9, 27, 82, 169, 492, 1498, 3002, 8779, 27136 \\ real+imag # GP-DEFINE z_points(k) = vector(7^k,n,n--; z_point(n)); # GP-DEFINE N_xmax_by_points(k) = my(n); vecmax(real(z_points(k)),&n); n-1; # GP-Test vector(5,k,k--; N_xmax_by_points(k)) == \ # GP-Test vector(5,k,k--; N_xmax(k)) # GP-Test z_point(302) == 10 - sqrt3i # GP-Test z_point(57) == 9 + 3*sqrt3i # GP-Test to_base7(302) == 611 # GP-Test to_base7(57) == 111 =pod For computer calculation these maximums can be calculated from the powers. The parts resulting can also be written in terms of the angle arg(b) = atan(sqrt(3)/5) = 19.106... degrees =cut # GP-DEFINE b_angle = arg(b); # GP-DEFINE b_angle_degrees = b_angle * 180/Pi; # GP-Test nearly_equal( b_angle, atan(sqrt3/5) ) # GP-Test b_angle_degrees > 19.106 # GP-Test b_angle_degrees < 19.106+1/10^3 # not in OEIS: 0.333473172251832115336090 \\ radians # not in OEIS: 19.1066053508690943945174 \\ degrees =pod For successive k, if adding this pushes the b^k angle past +30deg then the preceding digit goes past -30deg and becomes the new maximum X. Write the angle as a fraction of 60deg (pi/3), F = atan(sqrt(3)/5) / (pi/3) = 0.318443 ... =cut # GP-DEFINE angle_F = atan(sqrt3/5) / (Pi/3); # GP-Test angle_F > 0.318443 # GP-Test angle_F < 0.318443 + 1/10^6 # not in OEIS: 0.318443422514484906575291 =pod This is irrational since b^k is never on the X or Y axes. That can be seen since 2/sqrt3*imag(b^k) mod 7 goes in a repeating pattern 1,5,4,6,2,3. Similarly 2*real(b^k) mod 7 so not on the Y axis, and also anything on the Y axis would have 3*k fall on the X axis. =cut # GP-DEFINE is_integer(x) = (x==floor(x)); # GP-Test vector(100,k,k--; is_integer(imag(2*b^k))) == vector(100,k,1) # GP-Test vector(100,k,k--; imag(2*b^k)%7) == \ # GP-Test vector(100,k,k--; if(k==0,0, [1,5,4,6,2,3][(k-1)%6+1])) # # GP-Test vector(100,k,k--; is_integer(real(2*b^k))) == vector(100,k,1) # GP-Test vector(100,k,k--; real(2*b^k)%7) == \ # GP-Test vector(100,k,k--; if(k==0,2, [5, 4, 6, 2, 3, 1][(k-1)%6+1])) =pod Digits low to high are successive steps back cyclically 6,5,4,3,2,1 so that (with mod giving 0 to 5), N_xmax_digit(j) = (-floor(F*j+1/2) mod 6) + 1 =cut # GP-DEFINE N_xmax_digit_by_floor(j) = (-floor(angle_F*j+1/2) % 6) + 1; # GP-Test vector(1000,j,j--; N_xmax_digit_by_floor(j)) == \ # GP-Test vector(1000,j,j--; N_xmax_digit(j)) =pod The +1/2 is since initial direction b^0=1 is angle 0 which is half way between -30 and +30 deg. Similarly for the location, using conj(w6) for rotation back z_xmax_exp(j) = floor(F*j+1/2) = 0,0,1,1,1,2,2,2,3,3,3,4,4,4,4,5,5,5, ... z_xmax(k) = sum(j=0,k-1, conj(w6)^z_xmax_exp(j) * b^j) =cut # GP-DEFINE z_xmax_exp(j) = floor(angle_F*j+1/2); # GP-Test my(v=[0,0,1,1,1,2,2,2,3,3,3,4,4,4,4,5,5,5]); /* samples shown */ \ # GP-Test vector(#v,j,j--; z_xmax_exp(j)) == v # GP-DEFINE z_xmax_by_floor(k) = sum(j=0,k-1, conj(w6)^z_xmax_exp(j) * b^j); # GP-Test vector(200,j,j--; z_xmax_by_floor(j)) == \ # GP-Test vector(200,j,j--; z_xmax(j)) # # # vector(35,k,k++; z_xmax_exp(k)) \\ floor(angle_F*j+1/2)) # not in OEIS: 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7, 7, 8, 8, 8, 9, 9, 9, 10, 10, 10, 11, 11, 11, 11 # not A082964 a(n) = m given by arctan(tan(n)) = n - m*Pi. # GP-DEFINE A082964(n) = round((n-atan(tan(n)))/Pi); # # atan(tan(n)) gives fractional part -pi/2 to +pi/2, so how many revolutions # angle n makes around a circle, up to -pi/2, so factor 1/Pi # 1/Pi \\ 0.318309886183790671537767 is close to F # # GP-DEFINE A082964_by_floor(n) = floor(1/Pi*n+1/2); # GP-Test vector(10000,n,A082964(n)) == \ # GP-Test vector(10000,n,A082964_by_floor(n)) # GP-Test vector(1000,n,A082964(n)) != \ # GP-Test vector(1000,j, floor(angle_F*j+1/2)) =pod By symmetry the maximum extent is the same in 60deg, 120deg, etc directions, suitably rotated. The N in those cases has the digits 1,2,3,4,5,6 cycled around for the rotation. In PlanePath triangular X,Y coordinates direction 60deg means when sum X+3*Y is a maximum, etc. =cut # GP-DEFINE w12_times_sqrt3 = 1+w6; /* w12 * sqrt(3) */ # (x/2+y*sqrt3i/2) * conj(w6) == (x/4 + 3*y/4) + (-x/4 + y*1/4)*sqrt3i # (x/2+y*sqrt3i/2) * conj(w12_times_sqrt3) == (x*3/4 + y*3/4) + (-x/4 + y*3/4)*sqrt3i # GP-DEFINE z_to_x(z) = 2*real(z); # GP-DEFINE z_to_y(z) = 2*imag(z); # GP-Test z_to_x(z_point(1)) == 2 # GP-Test z_to_x(z_point(3)) == -1 # GP-Test z_to_y(z_point(3)) == 1 # GP-DEFINE N_s3max_by_points(k) = my(n); vecmax(real(z_points(k)/w6),&n); n-1; # GP-Test to_base7(N_s3max_by_points(3)) == 122 # GP-Test to_base7(N_s3max_by_points(4)) == 1122 =pod If the +1/2 in the floor is omitted then the effect is to find the maximum point in direction +30deg. In the PlanePath coordinates this means maximum sum S = X+Y. N_smax_digit(j) = (-floor(F*j) mod 6) + 1 = 1,1,1,1,6,6,6,5,5,5,4,4,4,3,3, ... k-1 N_smax(k) = digits N_smax_digit(j) low digit j=0 j=0 = 0, 1, 8, 57, 400, 14806, 115648, ... decimal = 0, 1, 11, 111, 1111, 61111, 661111, ... base7 and also N_smax() + 1 z_smax_exp(j) = floor(F*j) = 0,0,0,0,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,6, ... z_smax(k) = sum(j=0,k-1, conj(w6)^z_smax_exp(j) * b^j) = 0, 1, 7/2+1/2*sqrt3*i, 9+3*sqrt3*i, 19+12*sqrt3*i, ... and also z_smax() + w6^2 smax(k) = 2*real(z_smax(k)) + imag(z_smax(k))*2/sqrt3 = 0, 2, 8, 24, 62, 172, 470, 1190, 3202, 8740, ... coordinate sum X+Y max In the base figure, points 1 and 2 have the same X+Y=2 and this remains so in subsequent levels, so that for kE=1 N_smax(k) and N_smax(k)+1 are equal maximums. =cut # GP-DEFINE N_smax_digit(j) = (-floor(angle_F*j) % 6) + 1; # GP-Test my(v=[1,1,1,1,6,6,6,5,5,5,4,4,4,3,3]); /* samples shown */ \ # GP-Test vector(#v,j,j--; N_smax_digit(j)) == v # GP-DEFINE N_smax(k) = fromdigits(Vecrev(vector(k,j,j--; N_smax_digit(j))),7); # GP-Test N_smax(0) == 0 # GP-Test N_smax(1) == 1 # GP-Test N_smax(6) == 115648 # GP-Test to_base7(N_smax(51)) \ # GP-Test == 444555566611122233344455566661112223334445556661111 # GP-Test my(v=[0, 1, 8, 57, 400, 14806, 115648]); /* samples shown */ \ # GP-Test vector(#v,j,j--; N_smax(j)) == v # GP-Test my(v=[0, 1, 11, 111, 1111, 61111, 661111]); /* samples shown */ \ # GP-Test vector(#v,j,j--; to_base7(N_smax(j))) == v # vector(25,k,k--; N_smax_digit(k)) # vector(8,k, N_smax(k)) # vector(8,k, to_base7(N_smax(k))) # not in OEIS: 1,1,1,1,6,6,6,5,5,5,4,4,4,3,3,3,2,2,2,1,1,1,6,6,6 \\ digits # not in OEIS: 1, 8, 57, 400, 14806, 115648 \\ decimal # not in OEIS: 1, 11, 111, 1111, 61111, 661111 \\ base7 # vector(8,k, N_smax(k)+1) # vector(8,k, to_base7(N_smax(k))+1) # not in OEIS: 2, 9, 58, 401, 14807, 115649, 821543, 4939258 \\ decimal # not in OEIS: 2, 12, 112, 1112, 61112, 661112, 6661112, 56661112 \\ base7 # GP-DEFINE z_smax_exp(j) = floor(angle_F*j); # GP-Test my(v=[0,0,0,0,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,6]); /*samples shown*/ \ # GP-Test vector(#v,j,j--; z_smax_exp(j)) == v # GP-DEFINE z_smax(k) = sum(j=0,k-1, conj(w6)^z_smax_exp(j) * b^j); # GP-Test my(v=[0, 1, 7/2+1/2*sqrt3i, 9+3*sqrt3i, 19+12*sqrt3i]); /*samples*/ \ # GP-Test vector(#v,j,j--; z_smax(j)) == v # GP-Test vector(50,j,j--; real( z_smax(j) / w12_times_sqrt3 )) == \ # GP-Test vector(50,j,j--; real( (z_smax(j)+w6^2) / w12_times_sqrt3 )) # GP-DEFINE smax(k) = my(z=z_smax(k)); z_to_x(z)+z_to_y(z); # GP-Test my(v=[0, 2, 8, 24, 62, 172, 470, 1190, 3202, 8740]); /*samples*/ \ # GP-Test vector(#v,j,j--; smax(j)) == v # vector(50,k,k++; z_smax_exp(k)) \\ floor(angle_F*j) # not in OEIS: 4,4,4,5,5,5,6,6,6,7,7,7,7,8,8,8,9,9,9,10,10,10,11,11,11,12,12,12,13,13,13,14,14,14,14,15,15,15,16 # not A032615 = floor(n/Pi) # 1/Pi \\ = 0.318309886183790671537767 is close to F # GP-DEFINE A032615(n) = floor(1/Pi*n); # # is not A062300 which is same, almost, maybe, as A032615 after initial terms # A062300 a(n) = floor cosec( pi/(n+1) ) # GP-DEFINE A062300(n) = floor(1/sin(Pi/(n+1))); # GP-Test vector(10000,n,n+=4; A062300(n)) == \ # GP-Test vector(10000,n,n+=4; A032615(n+1)) # GP-Test vector(200,n,n+=4; A062300(n)) != \ # GP-Test vector(200,n,n+=4; z_smax_exp(n+1)) # sin(x)~x when x small so floor(1/sin(Pi/(n+1))) ~ floor((n+1)/Pi) # but with sin(x)n-1, Vec(select(e->e==z,v,1))); # GP-DEFINE } # GP-Test N_smax_list_by_points(0) == [0] # GP-Test N_smax_list_by_points(1) == [1,2] # GP-Test N_smax_list_by_points(2) == [8,9] # GP-Test N_smax_list_by_points(3) == [57,58] # GP-Test N_smax(3) == 57 # tan(n) # atan(tan(n)) # n-atan(tan(n)) # (n-atan(tan(n)))/Pi # n - m*Pi =pod =head1 SEE ALSO L, L, L, L, L, L =head1 HOME PAGE L =head1 LICENSE Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020 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 . =cut