Math-PlanePath-Toothpick
view release on metacpan or search on metacpan
lib/Math/PlanePath/ToothpickTree.pm view on Meta::CPAN
}
my %parts_to_numroots = (two_horiz => 2,
# everything else 1 root
);
sub tree_num_roots {
my ($self) = @_;
return ($parts_to_numroots{$self->{'parts'}} || 1);
}
sub tree_n_parent {
my ($self, $n) = @_;
### tree_n_parent(): $n
$n = int($n);
if ($n < ($parts_to_numroots{$self->{'parts'}} || 1)) {
return undef;
}
my ($x,$y) = $self->n_to_xy($n)
or return undef;
### parent at: "xy=$x,$y"
### parent odd list: (($x%2) ^ ($y%2)) && ($self->xy_to_n_list($x,$y-1), $self->xy_to_n_list($x,$y+1))
### parent even list: !(($x%2) ^ ($y%2)) && ($self->xy_to_n_list($x-1,$y), $self->xy_to_n_list($x+1,$y))
### parent min: min($self->xy_to_n_list($x-1,$y), $self->xy_to_n_list($x+1,$y),$self->xy_to_n_list($x,$y-1), $self->xy_to_n_list($x,$y+1))
return min((($x%2) ^ ($y%2))
?
# odd X,Y, vertical to parent
($self->xy_to_n_list($x,$y-1),
$self->xy_to_n_list($x,$y+1))
:
# even X,Y, horizontal to parent
($self->xy_to_n_list($x-1,$y),
$self->xy_to_n_list($x+1,$y)));
}
sub tree_n_to_depth {
my ($self, $n) = @_;
### tree_n_to_depth(): "$n"
if ($n < 0) {
return undef;
}
my ($depth) = _n0_to_depth_and_rem($self, int($n));
### n0 depth: $depth
return $depth;
}
# Do a binary search for the bits of depth which give Ndepth <= N.
#
# Ndepth grows as roughly depth*depth, so this is about log4(N) many
# compares. For large N wouldn't want to a table to search through to
# sqrt(N).
sub _n0_to_depth_and_rem {
my ($self, $n) = @_;
### _n0_to_depth_and_rem(): "n=$n parts=$self->{'parts'}"
# For parts=4 have depth=2^exp formula
# T[2^exp] = parts*(4^exp-1)*2/3 + 3
# parts*(4^exp-1)*2/3 + 3 = N
# 4^exp = (N-3)*3/2parts, round down
# but must be bigger ... (WHY-IS-IT-SO?)
#
my ($pow,$exp) = round_down_pow (12*$n, # /$self->{'parts'}
4);
if (is_infinite($exp)) {
return ($exp,0);
}
### $pow
### $exp
my $depth = 0;
my $n_depth = 0;
$pow = 2 ** $exp; # pow=2^exp down to 1, inclusive
while ($exp-- >= 0) {
my $try_depth = $depth + $pow;
my $try_n_depth = $self->tree_depth_to_n($try_depth);
### $depth
### $pow
### $try_depth
### $try_n_depth
if ($try_n_depth <= $n) {
### use this tried depth ...
$depth = $try_depth;
$n_depth = $try_n_depth;
}
$pow /= 2;
}
### _n0_to_depth_and_rem() final ...
### $depth
### remainder: $n - $n_depth
return ($depth, $n - $n_depth);
}
# First unsorted @pending
# depth=119 parts=4 pow=64 119
# depth=119 parts=4 pow=32 56,55
# depth=119 parts=4 pow=16 25,24,23
# depth=119 parts=4 pow=8 10,9,8,7 <- list crosses pow=8 boundary
# depth=119 parts=4 pow=4 3,2,7
# depth=119 parts=4 pow=2 3
# T(2^k+rem) = T(2^k) + T(rem) + 2T(rem-1) rem>=1
#------------------------------------------------------------------------------
# tree_depth_to_n()
# initial toothpicks not counted by the blocks crunching
my %depth_to_n_initial
= (4 => 3, # 1 origin + 1 above + 1 below
3 => 2, # 1 origin + 1 above
2 => 1, # 1 middle X=0,Y=1
lib/Math/PlanePath/ToothpickTree.pm view on Meta::CPAN
$n -= $parts*($depth-3);
}
while (--$exp >= 0) {
last unless @pending;
### @pending
### @mult
### $exp
### $pow
my @new_pending;
my @new_mult;
my $tpow;
# if (1||join(',',@pending) ne join(',',reverse sort {$a<=>$b} @pending)) {
# # print "depth=$depth parts=$parts pow=$pow ",join(',',@pending),"\n";
# print "mult ",join(',',@mult),"\n";
# }
foreach my $depth (@pending) {
my $mult = shift @mult;
### assert: $depth >= 2
### assert: $depth < 2*$pow
if ($depth <= 3) {
if ($depth eq '3') {
### depth==3 total=1 ...
$n += $mult;
} else {
### depth==2 total=0 ...
}
next;
}
if ($depth < $pow) {
# Smaller than $pow, keep unchanged. Cannot stop processing
# @pending on finding one $depth<$pow because @pending is not quite
# sorted and therefore might have a later $depth>=$pow.
push @new_pending, $depth;
push @new_mult, $mult;
next;
}
my $rem = $depth - $pow;
### $depth
### $mult
### $rem
### assert: $rem >= 0 && $rem < $pow
my $basemult = $mult; # multiple of oct(2^k) base part
if ($rem == 0) {
### rem==0, so just the oct(2^k) part ...
} elsif ($rem == 1) {
### rem==1 "A" ...
$n += $mult;
} else {
### rem >= 2, formula ...
# formula oct(pow+rem) = oct(pow) + oct(rem+1) + 2*oct(rem) - rem + 4
$n += (4-$rem)*$mult;
$rem += 1; # to give rem+1
if ($rem == $pow) {
### rem+1==pow so oct(2^k) by increasing basemult ...
$basemult += $mult;
} elsif (@new_pending && $new_pending[-1] == $rem) {
### combine rem+1 here with rem of previous ...
$new_mult[-1] += $mult;
} else {
push @new_pending, $rem;
push @new_mult, $mult;
}
if ($rem -= 1) { # to give plain rem again
push @new_pending, $rem;
push @new_mult, 2*$mult;
}
}
# oct(2^k) = (4^(k-1) - 4)/3 + 2^(k-1)
$tpow ||= ($pow*$pow - 16)/12 + $pow/2;
$n += $basemult * $tpow;
}
@pending = @new_pending;
@mult = @new_mult;
$pow /= 2;
}
### return: $n
return $n;
}
#------------------------------------------------------------------------------
# $depth numbered from origin in parts=4 style.
# Return cells added at that depth,
# ie. added = depth_to_n($depth+1) - depth_to_n($depth)
#
# @$depth_list is a list of $depth values.
# @mult_list is the multiple of T[depth] desired for that @$depth_list entry.
# $depth_list->[0], ie. the first array entry, must be the biggest $depth.
#
# @$depth_list is maintained mostly high to low and growing by one more
# value at each $exp level, but sometimes it's a bit more and some values
# not high to low and possibly duplicated.
#
# added(pow) = 1
# added(pow+1) = 2
# added(pow+rem) = 2*added(rem) + added(rem+1) - 1
#
# added(pow+pow-1) = 2*added(pow-1) + added(pow) - 1
# = 2*added(pow-1) + 1 - 1
# = 2*added(pow-1)
# repeats down to added(2^k-1) = 2^(k-1)
sub _depth_to_octant_added {
my ($depth_list, $mult_list, $zero) = @_;
### _depth_to_octant_added(): join(',',@$depth_list)
### assert: scalar(@$depth_list) >= 1
### assert: max(@$depth_list) == $depth_list->[0]
my ($pow,$exp) = round_down_pow ($depth_list->[0], 2);
if (is_infinite($exp)) {
return $exp;
}
### $pow
### $exp
my $add = $zero;
while (--$exp >= 0) { # running $pow down to 2 (inclusive)
### assert: $pow >= 2
last unless @$depth_list;
### pending: join(',',@$depth_list)
### mult : join(',',@$mult_list)
### $exp
### $pow
my @new_depth_list;
my @new_mult_list;
foreach my $depth (@$depth_list) {
my $mult = shift @$mult_list;
### assert: $depth >= 0
### assert: $depth == int($depth)
if ($depth <= 3) {
### depth==2or3 add=1 ...
$add += $mult;
next;
}
if ($depth < $pow) {
# less than 2^exp so unchanged
push @new_depth_list, $depth;
push @new_mult_list, $mult;
next;
}
my $rem = $depth - $pow;
### $depth
### $mult
### $rem
### assert: $rem >= 0 && $rem <= $pow
if ($rem == 0 || $rem == $pow-1) {
### rem==0, A of each, add=1 ...
### or depth=2*pow-1, add=1 ...
$add += $mult;
} else {
### rem >= 2, formula ...
# A(pow+rem) = A(rem+1) + 2A(rem) - 1
$add -= $mult;
$rem += 1; # to make rem+1
if (@new_depth_list && $new_depth_list[-1] == $rem) {
# add to previously pushed pending depth
# print "rem=$rem ",join(',',@new_depth_list),"\n";
$new_mult_list[-1] += $mult;
} else {
push @new_depth_list, $rem;
push @new_mult_list, $mult;
}
push @new_depth_list, $rem-1; # back to plain rem
push @new_mult_list, 2*$mult;
}
}
$depth_list = \@new_depth_list;
$mult_list = \@new_mult_list;
$pow /= 2;
}
### return: $add
return $add;
}
#------------------------------------------------------------------------------
sub tree_n_to_subheight {
my ($self, $n) = @_;
### tree_n_to_subheight(): $n
if ($n < 0) { return undef; }
if (is_infinite($n)) { return $n; }
my $zero = $n * 0;
(my $depth, $n) = _n0_to_depth_and_rem($self, int($n));
### $depth
### $n
my $parts = $self->{'parts'};
$depth += $parts_depth_adjust{$parts};
### depth adjusted to: $depth
if ($parts eq 'octant') {
my $add = _depth_to_octant_added ([$depth],[1], $zero);
$n = $add-1 - $n;
### octant mirror numbering to n: $n
} elsif ($parts eq 'octant_up') {
} elsif ($parts eq 'wedge' || $parts eq 'wedge+1') {
my $add = _depth_to_octant_added ([$depth],[1], $zero);
### wedge half width: $add
if ($parts eq 'wedge+1') {
# 0, 1 to $add, $add+1 to 2*$add, 2*$add+1
if ($n == 0 || $n == 2*$add+1) {
return 0; # first,last toothpicks don't grow
}
$n -= 1;
}
lib/Math/PlanePath/ToothpickTree.pm view on Meta::CPAN
--------+--------
/ | \
oct 2 behind / | \ oct 2 behind
/ | \
/ oct | oct \
The four octants near the Y axis begin immediately. The N=0 and N=2 points
are shared by the central octants going up and down on the right, and
likewise N=1,N=3 on the left.
The four octants near the X axis begin 3 row depth levels later. This is
not the same as the quadrants of the full pattern (their opposite octants
are 1 depth offset). For example the point N=10 at depth=3 is the start of
the lower octant in that quarter. A bigger picture than what's shown above
makes this easier to see.
=head2 Octant Vertical or Horizontal
The parts=octant pattern is half a parts=1 quadrant across the diagonal.
It's also interesting to note that an octant is half a quadrant by taking
just the vertical or horizontal toothpicks in the quadrant (taking vertical
or horizontal according to the orientation of the last row in the octant).
oct(d) = quad_verticals(d) + floor(d/2) if d odd
quad_horizontals(d) + floor(d/2) if d even
d = depth starting from 0 per tree_depth_to_n()
This works because in a quadrant the vertical toothpicks above the X=Y
diagonal can be folded down across the diagonal to become horizontal and
complete the lower octant.
quadrant verticals octant made by quadrant
numbered by row depth upper verticals fold down
to become horizontals
| | | | |
12 12 12 12 ......12
| | | | | | | |
10 10 ......10
| | | | | | | | |
12 8 8 12 ......8 -12--12
| | | | | | | |
6 ......6
| | | | | | | | |
4 4 8 12 ......4 --8---8 -12--12
| | | | | | | | | | | |
2 10 10 .....2 10--10--10
| | | | | | | | | | | |
0 4 12 0 --4---4 -12-- -12--12
| | | | | |
For example the vertical depth "4" toothpick which is above the X=Y diagonal
folds down to become the horizontal "4" in the lower octant. Similarly the
block of five 8,10,12,12,12 above the diagonal fold down to make five
horizontals. And the final 12 above becomes the horizontal 12.
However the horizontals which are on the central diagonal spine don't have
corresponding verticals above. These are marked "....." in the octant shown
above. This means 1 toothpick missing at every second row and therefore the
floor(depth/2) in the oct() formula above.
The key is that a quadrant has the upper octant running 1 growth row behind
the lower. So the horizontals in the lower correspond to the verticals in
the upper (and vice-versa).
The correspondence can be seen algebraically in the formula for a quadrant,
quad(d) = oct(d) + oct(d-1) + d
reversed to
oct(d) = quad(d) - oct(d-1)
and the oct(d-1) term repeatedly expanded
oct(d) = quad(d) - (quad(d-1) - oct(d-2) + d-1) + d
= quad(d)-quad(d-1)+1 + oct(d-2)
= ...
= quad(d)-quad(d-1)+1 + quad(d-2)-quad(d-3)+1 + ...
= quad(d)-quad(d-1) + quad(d-2)-quad(d-3) + ... + floor(d/2)
The difference quad(d)-quad(d-1) is the number of toothpicks added to make
depth d, and similarly quad(d-2)-quad(d-3) the number added to make depth
d-2. This means the number added at every second row, so if d is even then
this counts only the vertical toothpicks added. Or if d is odd then only
the horizontals.
The +d, +(d-1), +(d-2) additions from the quad(d) formula have alternating
signs and so cancel out to be +1 per pair, giving floor(d/2).
The parts=wedge pattern is two octants and therefore the wedge corresponds
to the horizontals or verticals of parts=2 which is two quadrants. But
there's an adjustment to make there though since parts=2 doesn't have a
toothpick at the origin the way the wedge does.
=head2 Quadrant and 2^k-1 Sums
In OEIS A168002 (L<http://oeis.org/A168002>) Omar Pol observed that the
quadrant(d) total cells taken mod 2 gives the number of ways d can be
expressed as a sum of terms 2^k-1.
d = (2^a - 1) + (2^b - 1) + (2^c - 1) + ...
distinct a,b,c,...
There's only ever 0 or 1 way to write d as a sum of 2^k-1 terms, ie. d
either is or is not such a sum. For example,
d ways
--- ----
8 1 8 = 7+1
9 0 no sum possible
10 1 10 = 7+3
11 1 11 = 7+3+1
12 0 no sum possible
The sum can be formed by taking the highest possible 2^k-1 from d
repeatedly. This works because smaller 2^k-1 terms are not big enough to
add up to that highest term. The result is a recurrence
ways(2^k-1) = 1
ways(2^k-1 + rem) = ways(rem) 1 <= rem < 2^k-1
ways(2*(2^k-1)) = 0)
The quadrant total cells follows a similar recurrence when taken mod 2.
Using the quadrant count Q(d) in the paper by Applegate, Pol, Sloane above
Q(d) = (T(d) - 3)/4
numbered same as T,
so Q(2)=0 then first quadrant toothpick Q(3)=1
Substituting the recurrences for T in the paper gives
Q(2^k) = (4^k-4)/6
Q(2^k+1) = Q(2^k) + 1
Q(2^k + rem) = Q(2^k) + Q(rem+1) + 2*Q(rem) + 2
for 2 <= rem < 2^k
Taking these modulo 2
Q(2^k) == 0 mod 2 since (4^k-4)/6 always even
Q(2^k+1) == 1 mod 2
Q(2^k + rem) == Q(rem+1) mod 2 2 <= rem < 2^k
Q(2^k-1 + rem) == Q(rem) mod 2 1 <= rem < 2^k-1
The last formula is the key, being the same as the ways(2^k-1 + rem)
recurrence above. Then Q(2^k)=0 corresponds to ways(2^k-2)=0, and
Q(2^k+1)=1 corresponding to ways(2^k-1)=1. And therefore
Q(d-2) == ways(d) mod 2
=head1 FUNCTIONS
See L<Math::PlanePath/FUNCTIONS> for behaviour common to all path classes.
=over 4
=item C<$path = Math::PlanePath::ToothpickTree-E<gt>new ()>
=item C<$path = Math::PlanePath::ToothpickTree-E<gt>new (parts =E<gt> $str)>
Create and return a new path object. C<parts> can be
"4" full pattern (the default)
"3" three quadrants
"2" half plane
"1" single quadrant
"octant" single eighth
"octant_up" single eighth upper
"wedge" V-shaped wedge
"two_horiz" starting two horizontal toothpicks
=back
=head2 Tree Methods
=over
=item C<@n_children = $path-E<gt>tree_n_children($n)>
Return the children of C<$n>, or an empty list if C<$n> has no children
(including when C<$n E<lt> 0>, ie. before the start of the path).
The children are the new toothpicks added at the ends of C<$n> at the next
row. This can be 0, 1 or 2 points. For example in the parts=4 default N=24
has no children, N=8 has a single child N=12, and N=2 has two children
N=4,N=5. The way points are numbered means that if there are two children
then they're consecutive N values.
=item C<$n_parent = $path-E<gt>tree_n_parent($n)>
Return the parent node of C<$n>, or C<undef> if no parent due to C<$n E<lt>=
0> (the start of the path).
=item C<$depth = $path-E<gt>tree_n_to_depth($n)>
=item C<$n = $path-E<gt>tree_depth_to_n($depth)>
Return the depth of point C<$n>, or first C<$n> at given C<$depth>,
respectively.
The first point N=0 is depth=0 in all "parts" forms. The way parts=1 and
parts=2 don't start at the origin means their depth at a given X,Y differs
by 2 or 1 respectively from the full pattern at the same point.
=back
=head2 Tree Descriptive Methods
=over
=item C<$num = $path-E<gt>tree_num_roots ()>
Return the number of root nodes in C<$path>. This is 1 except for
parts=two_horiz which is 2.
=item C<$num = $path-E<gt>tree_num_children_minimum()>
=item C<$num = $path-E<gt>tree_num_children_maximum()>
Return minimum 0 and maximum 2 since each node has 0, 1 or 2 children.
=back
=head2 Level Methods
=over
=item C<($n_lo, $n_hi) = $path-E<gt>level_to_n_range($level)>
Return C<(0, tree_depth_to_n_end($depth)> where the depth for a completed
level is
parts depth
----- -----
4 \
3 | 4*2^level - 1 = 3, 7, 15, 31, ...
wedge |
two_horiz /
2 4*2^level - 2 = 2, 6, 14, 30, ...
1 4*2^level - 3 = 1, 5, 13, 29, ...
octant \ 4*2^level - 4 = 0, 4, 12, 28, ...
octant_up /
parts=octant is one depth less than parts=1 because the lower eighth is one
row ahead of the upper, so parts=1 finishes one later.
parts=octant_up is the upper eighth of parts=1 but one depth less because
the octant starts at X=0,Y=1 which is one row later than parts=1.
In each case the depth is reckoned by the slowest eighth in the parts
pattern. For example parts=two_horiz completes levels of the eighths
nearest the X axis (the "oct 2 behind" shown in L</Two Horizontal> above).
=back
=head1 FORMULAS
=head2 Depth to N
The first N at given depth is the total count of toothpicks in the preceding
rows. The paper by Applegate, Pol and Sloane above gives formulas for
parts=4 and parts=1. A similar formula can be made for parts=octant.
The depth in all the following is per the full pattern, which means the
octant starts at depth=2. So oct(2)=0 then oct(3)=1. This reckoning keeps
the replications on 2^k boundaries and is convenient for relating an octant
to the full pattern. Note though that C<tree_depth_to_n()> always counts
from C<$depth = 0> so an adjustment +1 or +2 is applied there.
for depth >= 2
depth = pow + rem where pow=2^k is the high bit of depth
so 0 <= rem < 2^k
oct(2) = 0
oct(pow) = (pow*pow - 16)/12 + pow/2
oct(pow+1) = oct(pow) + 1
oct(pow+rem) = oct(pow) + oct(rem+1) + 2*oct(rem) - rem + 4
for 2 <= rem < pow
The other parts patterns can be expressed in terms of an octant. It's
convenient to make an octant the unit and have the others as multiples and
depth offsets from it.
quad(d) = oct(d) + oct(d-1) - d + 3
half(d) = 2*quad(d) + 1
full(d) = 4*quad(d) + 3
3corner(d) = quad(d+1) + 2*quad(d) + 2
= oct(d+1) + 3*oct(d) + 2*oct(d-1) - 3*d + 10
wedge(d) = 2*oct(d-1) + 4
In quad(d) the "-d" term adjusts for the stairstep diagonal spine being
counted twice by the oct(d)+oct(d-1).
The oct() recurrence corresponds to the sub-block breakdown shown under
L</One Octant> above.
oct(pow+rem) = oct(pow) "base"
+ oct(rem+1) "lower"
+ 2*oct(rem) "upper" and "extend"
- rem + 4 unduplicate diagonal
The stairstep diagonal between the "upper" and "lower" parts is duplicated
by those two parts, hence "-(rem-1)" to subtract one copy of it. A further
+3 is the points in-between the replications, ie. the "A", "B" and one
further toothpick not otherwise counted by the replications.
oct(rem+1) + 2*oct(rem) is the important part of the recurrence. It removes
the high bit of depth and spreads rem to an adjacent pair of smaller depths
rem+1 and rem. A list of pending depth values can be maintained and
compared to a pow=2^k for reduction.
for each pending depth
if depth == 2^k or 2^k+1 then oct(pow) or oct(pow+1)
if depth >= 2^k+2 then reduce by recurrence
repeat for 2^(k-1)
rem+1,rem are adjacent so successive reductions make a list growing by one
further value each time, like
d
d+1, d
d+2, d+1, d
( run in 0.429 second using v1.01-cache-2.11-cpan-39bf76dae61 )