# 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 = 120;
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