Math-PlanePath
view release on metacpan or search on metacpan
devel/c-curve.pl view on Meta::CPAN
sub points_str {
my ($points) = @_;
### points_str(): $points
my $count = scalar(@$points);
return "count=$count ".join(' ',map{join(',',@$_)}@$points)
}
}
{
my @sdir = (2,2,0,-2, -2,-2,0,2);
sub s0_by_formula {
my ($k) = @_;
{
my $h = int($k/2);
return 2**$k/4 + $sdir[$k%8]*2**$h/4;
}
{
# (1/4)*(2^k + (1-I)^k + (1+I)^k))
require Math::Complex;
return (2**$k / 2
+ (Math::Complex->new(1,-1)**$k * Math::Complex->new(1,-1)
+ Math::Complex->new(1,1) **$k * Math::Complex->new(1,1)) / 4);
}
}
my @s1dir = (0,2,2,2, 0,-2,-2,-2);
sub s1_by_formula {
my ($k) = @_;
my $h = int($k/2);
return 2**$k/4 + $sdir[($k-2)%8]*2**$h/4;
}
my @s2dir = (0,2,2,2, 0,-2,-2,-2);
sub s2_by_formula {
my ($k) = @_;
my $h = int($k/2);
return 2**$k/4 + $sdir[($k-4)%8]*2**$h/4;
}
sub s3_by_formula {
my ($k) = @_;
my $h = int($k/2);
return 2**$k/4 + $sdir[($k-6)%8]*2**$h/4;
}
# print " 1, 1, 1, 1, 2, 6, 16, 36, 72,136,256,496,992,2016,4096,8256,16512,32896,65536,\n"; # M0
# print " 0, 1, 2, 3, 4, 6, 12, 28, 64,136,272,528,1024,2016,4032,8128,16384,32896,65792,\n"; # M1
# print " 0, 0, 1, 3, 6, 10, 16, 28, 56,120,256,528,1056,2080,4096,8128,16256,32640,65536,\n"; # M2
print " 0, 0, 0, 1, 4, 10, 20, 36, 64,120,240,496,1024,2080,4160,8256,16384,32640,65280,\n"; # M3
foreach my $k (0 .. 17) {
printf "%3d,", s3_by_formula($k);
}
exit 0;
}
devel/c-curve.pl view on Meta::CPAN
my $want = $want[$k];
my $got = $L->($k);
my $diff = $want - $got;
print "$k $want $got $diff\n";
# print "$got,";
}
exit 0;
}
{
# right boundary formula vs recurrence
# R[k] = 2*R[k-1] + R[k-2] - 4*R[k-3] + 2*R[k-4]
#
# R[2k] = 4*R[2k-2] - 5*R[2k-4] + 2*R[2k-6]
# R[2k+1] = 4*R[2k-1] - 5*R[2k-3] + 2*R[2k-5]
my $R;
$R = sub {
my ($k) = @_;
if ($k < 4) { return R_formula($k); }
return (2*$R->($k-4)
- 4*$R->($k-3)
+ $R->($k-2)
+ 2*$R->($k-1));
};
require Memoize;
$R = Memoize::memoize($R);
my $R2;
$R2 = sub {
my ($k) = @_;
if ($k < 3) { return R_formula(2*$k); }
return (2*$R2->($k-3)
- 5*$R2->($k-2)
+ 4*$R2->($k-1));
};
require Memoize;
$R2 = Memoize::memoize($R2);
my $R2P1;
$R2P1 = sub {
my ($k) = @_;
if ($k < 3) { return R_formula(2*$k+1); }
return (2*$R2P1->($k-3)
- 5*$R2P1->($k-2)
+ 4*$R2P1->($k-1));
};
require Memoize;
$R2P1 = Memoize::memoize($R2P1);
foreach my $k (0 .. 50) {
# my $want = R_formula($k);
# print "$k $want ",$R->($k),"\n";
my $want = R_formula(2*$k);
print "$k $want ",$R2->($k),"\n";
}
exit 0;
}
{
# right outer boundary with sqrt(2)
sub S_formula {
my ($h) = @_;
return 2**$h;
};
sub Z_formula {
my ($h) = @_;
return 2*2**$h - 2;
};
my $S_cum = sub { # sum S[0] .. S[h] inclusive
my ($h) = @_;
return 2**($h+1) - 1;
};
my $Z_cum = sub { # sum Z[0] .. Z[h] inclusive
my ($h) = @_;
return 2*(2**($h+1) - 1) - 2*($h+1);
};
my $S_inR = sub {
my ($k) = @_;
my ($h, $rem) = Math::PlanePath::_divrem($k,2);
if ($rem) {
return 2*$S_cum->($h);
} else {
return 2*$S_cum->($h-1) + S_formula($h);
}
};
my $Z_inR = sub {
my ($k) = @_;
my ($h, $rem) = Math::PlanePath::_divrem($k,2);
if ($rem) {
return 2*$Z_cum->($h-1) + Z_formula($h);
} else {
return 2*$Z_cum->($h-1);
}
};
my $R_bySZ = sub {
my ($k) = @_;
return $S_inR->($k) + $Z_inR->($k);
};
{
my $total = 0;
foreach my $h (0 .. 10) {
$total == $S_cum->($h-1) or die;
$total += S_formula($h);
}
}
{
my $total = 0;
foreach my $h (0 .. 10) {
$total == $Z_cum->($h-1) or die;
$total += Z_formula($h);
}
}
# {
# print $S_cum->(-1),"\n";
# foreach my $h (0 .. 10) {
# print "+ ",S_formula($h), " = ",$S_cum->($h),"\n";
# }
# }
{
foreach my $k (0 .. 10) {
my $s = $S_inR->($k);
my $z = $S_inR->($k);
my $rby = $R_bySZ->($k);
my $rformula = R_formula($k);
# print "$k $s + $z = $rby rf=$rformula\n";
$rby == $rformula or die "$k $rby $rformula";
}
}
{
foreach my $k (0 .. 100) {
my $s = $S_inR->($k);
my $z = $S_inR->($k);
my $t = sqrt(2)**$k;
my $f = ($s + $z/sqrt(2)) / $t;
print "$k $s + $z f=$f\n";
}
print "2+ 2*sqrt(2)=",2 + 2*sqrt(2),"\n"; # k odd
print "3+3/2*sqrt(2)=",3+1.5*sqrt(2),"\n"; # k even
}
exit 0;
}
{
# right outer boundary
sub R_formula {
my ($k) = @_;
my $h = int($k/2);
return ($k & 1
? 10*2**$h - 2*$k - 6 # yes
: 7*2**$h - 2*$k - 6); # yes
if ($k & 1) {
my $j = ($k+1)/2;
return 5*2**$j - 4*$j - 4; # yes
devel/c-curve.pl view on Meta::CPAN
return 2*$Scum_recurrence->($k-1) + $k; # yes
};
$R_recurrence = sub {
my ($k) = @_;
if ($k == 0) { return 1; }
if ($k == 1) { return 2; }
return 2*$R_recurrence->($k-2) + 2*$k-2;
};
for (my $k = 0; $k < 15; $k++) {
print R_formula($k),", ";
}
print "\n";
require MyOEIS;
my $path = Math::PlanePath::CCurve->new;
foreach my $k (0 .. 17) {
my $n_end = 2**$k;
my $p = MyOEIS::path_boundary_length($path, $n_end, side => 'right');
# my $b = $B->($k);
my $srec = $Scum_recurrence->($k);
my $zrec = $Zcum_recurrence->($k);
my $rszrec = $srec + $zrec + 1;
my $rrec = $R_recurrence->($k);
# my $t = $T->($k);
# my $u = $U->($k);
# my $u2 = $U2->($k);
# my $u_lr = $U_from_LsubR->($k);
# my $v = $V->($k);
my ($s, $z) = path_S_and_Z($path, $n_end);
my $r = $s + $z + 1;
my $rformula = R_formula($k);
my $drformula = $r - $rformula;
# next unless $k & 1;
print "$k $p $s $z $r $srec $zrec $rszrec $rrec $rformula small by=$drformula\n";
}
exit 0;
sub path_S_and_Z {
my ($path, $n_end) = @_;
### path_S_and_Z(): $n_end
my $s = 0;
my $z = 0;
my $x = 1;
my $y = 0;
devel/chan-tree.pl view on Meta::CPAN
$poly *= make_poly_k4($level);
foreach my $i (0 .. 30) {
print " ",$poly->coeff($i);
}
print "\n";
}
exit 0;
}
{
# children formulas
foreach my $k (3 .. 8) {
my $half_ceil = int(($k+1) / 2);
foreach my $digit (0 .. $k-1) {
my $c1 = ($digit < $half_ceil ? $digit+1 : $k-$digit);
my $c0 = ($digit <= $half_ceil ? $digit : $k-$digit+1);
my $c2 = ($digit < $half_ceil-1 ? $digit+2 : $k-$digit-1);
print "${c1}x + ${c0}y / ${c2}x + ${c1}y\n";
}
print "\n";
}
devel/complex-minus.pl view on Meta::CPAN
$max_pos = "$x,$y";
}
}
# print "$min_hypot,";
# print " min $min_hypot at $min_x,$min_y\n";
# print " max $max_hypot at $max_x,$max_y\n";
{
my $factor = $min_hypot / $prev_min;
print " min r^2 $min_hypot 0b".sprintf('%b',$min_hypot)." at $min_pos factor $factor\n";
print " cf formula ", 2**($level-7), "\n";
}
# {
# my $factor = $max_hypot / $prev_max;
# print " max r^2 $max_hypot 0b".sprintf('%b',$max_hypot)." at $max_pos factor $factor\n";
# }
$prev_min = $min_hypot;
$prev_max = $max_hypot;
}
exit 0;
}
devel/complex-revolving.pl view on Meta::CPAN
# print " max $max_hypot at $max_x,$max_y\n";
{
my $factor = $min_hypot / $min[$level-1];
my $factor4_level = max($level-4,0);
my $factor4 = $min_hypot / $min[max($factor4_level)];
# printf " min r^2 %5d", $min_hypot;
printf " 0b%-20b", $min_hypot;
# print " at $min_pos";
# print " factor $factor";
# print " factor[$factor4_level] $factor4";
# print " cf formula ", 2**($level-7), "\n";
print "\n";
}
# {
# my $factor = $max_hypot / $prev_max;
# print " max r^2 $max_hypot 0b".sprintf('%b',$max_hypot)." at $max_pos factor $factor\n";
# }
$prev_max = $max_hypot;
}
exit 0;
}
devel/dragon.pl view on Meta::CPAN
return ($path->_UNDOCUMENTED_level_to_right_line_boundary($level+1)
- $path->_UNDOCUMENTED_level_to_left_line_boundary($level+1)) / 4;
return ($path->_UNDOCUMENTED_level_to_line_boundary($level) / 2
- $path->_UNDOCUMENTED_level_to_line_boundary($level+1) / 4);
return ($path->_UNDOCUMENTED_level_to_enclosed_area($level+1)
- 2*$path->_UNDOCUMENTED_level_to_enclosed_area($level));
}
sub level_to_join_points_by_formula {
my ($level) = @_;
{
if ($level == 0) { return 1; }
if ($level == 1) { return 1; }
if ($level == 2) { return 1; }
if ($level == 3) { return 2; }
my $j0 = 1;
my $j1 = 1;
my $j2 = 1;
my $j3 = 2;
devel/dragon.pl view on Meta::CPAN
# foreach my $n (0 .. $n_end) {
# my ($x,$y) = $path->n_to_xy($n);
# $seen{"$x,$y"}++;
# }
my $u = $path->_UNDOCUMENTED_level_to_u_left_line_boundary($k);
my $ru = $path->_UNDOCUMENTED_level_to_u_right_line_boundary($k);
my $bu = $path->_UNDOCUMENTED_level_to_u_line_boundary($k);
my $ja = level_to_join_area($k);
my $join_points = path_level_to_join_points($path,$k);
my $join_area = $join_points - 1;
my $j = level_to_join_points_by_formula($k);
my $da = level_to_denclosed($k);
my $area = $path->_UNDOCUMENTED_level_to_enclosed_area($k);
my $area_next = $path->_UNDOCUMENTED_level_to_enclosed_area($k+1);
my $darea = $area_next - $area;
my $v = $path->_UNDOCUMENTED_level_to_visited($k);
my $visited = $v; # MyOEIS::path_n_to_visited($path,$n_end);
my $dvisited = $visited - $prev_visited;
my $singles = 0 && MyOEIS::path_n_to_singles($path, $n_end-1);
my $doubles = 0 && MyOEIS::path_n_to_doubles($path, $n_end-1);
devel/dragon.pl view on Meta::CPAN
print 3*$xmin/$len+.001," / 3\n";
print 6*$xmax/$len+.001," / 6\n";
print 3*$ymin/$len+.001," / 3\n";
print 3*$ymax/$len+.001," / 3\n";
}
{
# A073089 midpoint vertical/horizontal formula
require Math::NumSeq::OEIS::File;
my $A073089 = Math::NumSeq::OEIS::File->new (anum => 'A073089');
my $A014577 = Math::NumSeq::OEIS::File->new (anum => 'A014577'); # 0=left n=0
my $A014707 = Math::NumSeq::OEIS::File->new (anum => 'A014707'); # 1=left
my $A038189 = Math::NumSeq::OEIS::File->new (anum => 'A038189');
my $A082410 = Math::NumSeq::OEIS::File->new (anum => 'A082410');
my $A000035 = Math::NumSeq::OEIS::File->new (anum => 'A000035'); # n mod 2
devel/koch-squareflakes.pl view on Meta::CPAN
my $a = $polygon->area;
my $len = $polygon->perimeter;
my $a_log = log($a);
my $len_log = log($len);
my $d_a_log = $a_log - $prev_a_log;
my $d_len_log = $len_log - $prev_len_log;
my $f = $d_a_log / $d_len_log;
my $formula = area_by_formula($level);
print "$level $a $formula\n";
# print "$level $d_len_log $d_a_log $f\n";
push @values, $a;
$prev_a_log = $a_log;
$prev_len_log = $len_log;
}
shift @values;
require Math::OEIS::Grep;
Math::OEIS::Grep->search(array => \@values);
exit 0;
sub area_by_formula {
my ($n) = @_;
return (9**$n - 4**$n)/5;
# return (4 * (9**$n - 4**$n)/5 + 16**$n);
}
}
{
# max extents of a single side
# horiz: 1, 4, 14, 48, 164, 560, 1912, 6528, 22288, 76096, 259808, 887040
# A007070 a(n+1) = 4*a(n) - 2*a(n-1), starting 1,4
devel/lib/Math/PlanePath/SquaRecurve.pm view on Meta::CPAN
# vector(25,n, (-1)^valuation(n,3))
# not in OEIS: 1,1,-1,1,1,-1,1,1,1,1,1,-1,1,1,-1,1,1,1,1,1,-1,1,1,-1,1,1,-1,1
# vector(100,n, valuation(n,3)%2)
# A182581 num ternary low 0s mod 2
=pod
The power of -1 means left or right flip for each low ternary 0 of N, and
flip again if N is odd. Odd N is an odd number of ternary 1 digits.
This formula follows from the turns in a new low base-9 digit. The start
and end of the base figure are in the same directions so the turns at 9*N
are unchanged. Then 9*N+r goes as r in the base figure, but flipped
LE<lt>-E<gt>R when N odd since blocks are mirrored alternately.
turn(9N) = turn(N)
turn(9N+r) = turn(r)*(-1)^N for 1 <= r <= 8
=cut
# GP-Test vector(900,n, turn(9*n)) == \
devel/quadric.pl view on Meta::CPAN
}
{
my $factor = $max_width / $prev_max;
print " max width $max_width oct ".sprintf('%o',$max_width)." at $max_pos factor $factor\n";
}
{
my $factor = $min_width / ($prev_min||1);
print " min width $min_width oct ".sprintf('%o',$min_width)." at $min_pos factor $factor\n";
}
{
my $formula = (2*4**($level-1) + 1) / 3;
print " cf min formula $formula\n";
}
{
my $formula = (10*4**($level-1) - 1) / 3;
print " cf max formula $formula\n";
}
$prev_max = $max_width;
$prev_min = $min_width;
}
exit 0;
}
{
# min/max for level
require Math::PlanePath::QuadricCurve;
devel/quadric.pl view on Meta::CPAN
$min_width = $w;
$min_pos = "$x,$y n=$n (oct ".sprintf('%o',$n).")";
}
}
# print " max $max_width at $max_x,$max_y\n";
my $factor = $max_width / $prev_max;
print " min width $min_width oct ".sprintf('%o',$min_width)." at $min_pos factor $factor\n";
# print " max width $max_width oct ".sprintf('%o',$max_width)." at $max_pos factor $factor\n";
# print " cf formula ",(10*4**($level-1) - 1)/3,"\n";
# print " cf formula ",2* (4**($level-0) - 1)/3,"\n";
print " cf formula ",2*4**($level-1),"\n";
$prev_max = $max_width;
}
exit 0;
}
devel/r5-dragon.pl view on Meta::CPAN
# [a] [v0 v1 v2] -1 [v1]
# [b] = [v1 v2 v3] * [v2]
# [c] [v2 v3 v4] [v3]
$|=1;
my @array = (
54,90,150,250,422,714,1206,2042,3462
);
# @array = ();
# foreach my $k (5 .. 10) {
# push @array, R_formula(2*$k+1);
# }
# require MyOEIS;
# my $path = Math::PlanePath::R5DragonCurve->new;
# foreach my $k (0 .. 30) {
# my $value = MyOEIS::path_boundary_length($path, 5**$k,
# # side => 'left',
# );
# last if $value > 10_000;
# push @array, $value;
# print "$value,";
devel/r5-dragon.pl view on Meta::CPAN
foreach my $i (0 .. 10) { print wrong($i),","; }
print "\n";
print "recurrence_area815(): ";
foreach my $i (0 .. 10) { print recurrence_area815($i),","; }
print "\n";
print "recurrence_area43(): ";
foreach my $i (0 .. 10) { print recurrence_area43($i),","; }
print "\n";
print "formula_pow(): ";
foreach my $i (0 .. 10) { print formula_pow($i),","; }
print "\n";
print "recurrence_areaSU(): ";
foreach my $i (0 .. 10) { print recurrence_areaSU($i),","; }
print "\n";
print "recurrence_area2S(): ";
foreach my $i (0 .. 10) { print recurrence_area2S($i),","; }
print "\n";
exit 0;
devel/r5-dragon.pl view on Meta::CPAN
sub wrong {
my ($n) = @_;
if ($n <= 0) { return 0; }
if ($n == 1) { return 0; }
return 4*wrong($n-1) + 4*5**($n-2);
}
# A[n] = (5^k - 2*3^k + 1)/2
sub formula_pow {
my ($n) = @_;
return (5**$n - 2*3**$n + 1) / 2;
}
sub recurrence_area43 {
my ($n) = @_;
if ($n <= 0) { return 0; }
if ($n == 1) { return 0; }
return 4*recurrence_area43($n-1) - 3*recurrence_area43($n-2) + 4*5**($n-2);
}
devel/sierpinski-curve.pl view on Meta::CPAN
for (my $n = $n_hi; $n >= 0; $n-=2) {
my ($x,$y) = $path->n_to_xy($n);
push @points, [$x,$y];
}
### @points
my $polygon = Math::Geometry::Planar->new;
$polygon->points(\@points);
my $area = $polygon->area;
my $length = $polygon->perimeter;
my $formula = two_area_by_formula($level);
print "$level $length area $area $formula\n";
push @values, $area;
}
shift @values;
require Math::OEIS::Grep;
Math::OEIS::Grep->search(array => \@values);
exit 0;
sub two_area_by_formula {
my ($k) = @_;
return (13*4**$k - 7) / 3;
$k++;
return (27*4**($k-1) - 18*2**($k-1) -14 * 4**($k-1) + 18 * 2**($k-1) - 4) / 3 - 1;
return 9*4**($k-1) - 6*2**($k-1) + (-14 * 4**($k-1) + 18 * 2**($k-1) - 4) / 3 - 1;
return 9*4**($k-1) - 6*2**($k-1) - (14 * 4**($k-1) - 18 * 2**($k-1) + 4) / 3 - 1;
return 9*2**(2*$k-2) - 2*3*2**($k-1) + 1 - (14 * 4**($k-1) - 18 * 2**($k-1) + 4) / 3 - 2;
return (3*2**($k-1) - 1)**2 - (14 * 4**($k-1) - 18 * 2**($k-1) + 4) / 3 - 2;
return (3*2**($k-1) - 1)**2 - 4*(7 * 4**($k-1) - 9 * 2**($k-1) + 2) / 6 - 2;
return (3*2**($k-1) - 2 + 1)**2 - 4*((7 * 4**($k-1) - 9 * 2**($k-1) + 2) / 6) - 2;
return (3*2**($k-1) - 2 + 1)**2 - 4*area_by_formula($k-1) - 2;
}
}
{
# area under the curve
# (7*4^k - 9*2^k + 2)/6
#
require Math::Geometry::Planar;
my $path = Math::PlanePath::SierpinskiCurve->new;
my @values;
devel/sierpinski-curve.pl view on Meta::CPAN
my @points;
foreach my $n ($n_lo .. $n_hi) {
my ($x,$y) = $path->n_to_xy($n);
push @points, [$x,$y];
}
my $polygon = Math::Geometry::Planar->new;
$polygon->points(\@points);
my $area = $polygon->area;
my $length = $polygon->perimeter;
my $formula = area_by_formula($level);
print "$level $length area $area $formula\n";
push @values, $area;
}
shift @values;
require Math::OEIS::Grep;
Math::OEIS::Grep->search(array => \@values);
exit 0;
sub area_by_formula {
my ($k) = @_;
if ($k == 0) { return 0; }
return (7 * 4**$k - 9 * 2**$k + 2) / 6;
return (28 * 4**($k-1) - 18 * 2**($k-1) + 2) / 6;
{
return (
4**($k-1) * 6
- 4**($k-1) + 1
devel/sierpinski-curve.pl view on Meta::CPAN
}
{
# A(k) = 4*A(k-1) + Ytop(k) + 1 k >= 2
# A(1) = 2
# A(0) = 0
# Ytop(k) = 3*2^(k-1) - 2
# A(k) = 4*A(k-1) + 3*2^(k-1) - 1
#
if ($k == 1) { return 2; }
return 4*area_by_formula($k-1) + 3 * 2**($k-1) - 1;
}
}
}
{
# Y coordinate sequence
require Math::NumSeq::PlanePathCoord;
my $seq = Math::NumSeq::PlanePathCoord->new (planepath => 'SierpinskiCurve',
devel/theodorus.pl view on Meta::CPAN
exit 0;
}
{
{
package Math::Symbolic::Custom::MySimplification;
use base 'Math::Symbolic::Custom::Simplification';
# use Math::Symbolic::Custom::Pattern;
# my $formula = Math::Symbolic->parse_from_string("TREE_a * (TREE_b / TREE_c)");
# my $pattern = Math::Symbolic::Custom::Pattern->new($formula);
use Math::Symbolic::Custom::Transformation;
my $trafo = Math::Symbolic::Custom::Transformation::Group->new
(',',
'TREE_a * (TREE_b / TREE_c)' => '(TREE_a * TREE_b) / TREE_c',
'TREE_a * (TREE_b + TREE_c)' => 'TREE_a * TREE_b + TREE_a * TREE_c',
'(TREE_b + TREE_c) * TREE_a' => 'TREE_b * TREE_a + TREE_c * TREE_a',
# '(TREE_a / TREE_b) / TREE_c' => 'TREE_a / (TREE_b * TREE_c)',
examples/koch-svg.pl view on Meta::CPAN
# You should have received a copy of the GNU General Public License along
# with Math-PlanePath. If not, see <http://www.gnu.org/licenses/>.
# Usage: perl koch-svg.pl >output.svg
# perl koch-svg.pl LEVEL >output.svg
#
# Print SVG format graphics to standard output for a Koch snowflake curve of
# given LEVEL fineness. The default level is 4.
#
# The range of N values to plot follows the formulas in the
# Math::PlanePath::KochSnowflakes module POD.
#
# The svg output size is a fixed 300x300, but of course the point of svg is
# that it can be resized by a graphics viewer program.
use 5.006;
use strict;
use warnings;
use List::Util 'min';
use Math::PlanePath::KochSnowflakes;
lib/Math/NumSeq/PlanePathTurn.pm view on Meta::CPAN
use constant _NumSeq_Turn_Turn4_max => 3.5;
}
# { package Math::PlanePath::DragonMidpoint;
# }
{ package Math::PlanePath::AlternatePaper;
use constant _NumSeq_Turn_Turn4_min => 1; # left or right only
# A209615 is (-1)^e for each p^e prime=4k+3 or prime=2
# 3*3 mod 4 = 1 mod 4
# so picks out bit above lowest 1-bit, and factor -1 if an odd power-of-2
# which is the AlternatePaper turn formula
#
use constant _NumSeq_Turn_oeis_anum =>
{ 'arms=1' =>
{ LSR => 'A209615',
# OEIS-Catalogue: A209615 planepath=AlternatePaper turn_type=LSR
Right => 'A292077',
# OEIS-Catalogue: A292077 planepath=AlternatePaper turn_type=Right
# # Not quite, A106665 has OFFSET=0 cf first here i=1
# 'Left' => 'A106665', # turn, 1=left,0=right
lib/Math/PlanePath.pm view on Meta::CPAN
The sqrt(3) factor can be worked into a hypotenuse radial distance
calculation as follows if comparing distances from the origin.
hypot = sqrt(X*X + 3*Y*Y)
See for instance C<TriangularHypot> which is triangular points ordered by
this radial distance.
=head1 FORMULAS
The formulas section in the POD of each class describes some of the
calculations. This might be of interest even if the code is not.
=head2 Triangular Calculations
For a triangular lattice, the rotation formulas above allow calculations to
be done in the rectangular X,Y coordinates which are the inputs and outputs
of the PlanePath functions. Another way is to number vertically on a 60
degree angle with coordinates i,j,
...
* * * 2
* * * 1
* * * j=0
i=0 1 2
lib/Math/PlanePath.pm view on Meta::CPAN
this_dX,this_dY at int
next_dX,next_dY at int+1
at fractional N
dX = this_dX * (1-frac) + next_dX * frac
dY = this_dY * (1-frac) + next_dY * frac
This is combination of this_dX,this_dY and next_dX,next_dY in proportion to
the distances from positions N to int+1 and from int+1 to N+1.
The formulas can be rearranged to
dX = this_dX + frac*(next_dX - this_dX)
dY = this_dY + frac*(next_dY - this_dY)
which is like dX,dY at the integer position plus fractional part of a turn
or change to the next dX,dY.
=head2 N to dX,dY -- Self-Similar
For most of the self-similar paths such as C<HilbertCurve>, the change dX,dY
lib/Math/PlanePath/AlternatePaper.pm view on Meta::CPAN
Return C<(0, 2**$level)>, or for multiple arms return C<(0, $arms *
2**$level + ($arms-1))>.
This is the same as L<Math::PlanePath::DragonCurve/Level Methods>. Each
level is an unfold (on alternate sides left or right).
=back
=head1 FORMULAS
Various formulas for coordinates, lengths and area can be found in the
author's mathematical write-up
=over
L<http://user42.tuxfamily.org/alternate/index.html>
=back
=head2 Turn
lib/Math/PlanePath/AlternatePaper.pm view on Meta::CPAN
= / GRS(k) if k even
\ -GRS(k) if k odd
For dY reducing 2k+1 to k drops a 1 bit from the low end. If the second
lowest bit is also a 1 then they were a "11" bit pair which is lost from
GRS(k). The factor (-1)^k adjusts for that, being +1 if k even or -1 if k
odd.
=head2 dSum
From the dX and dY formulas above it can be seen that their sum is simply
GRS(N),
dSum = dX + dY = GRS(N)
The sum X+Y is a numbering of anti-diagonal lines,
| \ \ \
|\ \ \ \
| \ \ \ \
|\ \ \ \ \
lib/Math/PlanePath/AlternateTerdragon.pm view on Meta::CPAN
There are 3^level segments in a curve level, so 3^level+1 points numbered
from 0. For multiple arms there are arms*(3^level+1) points, numbered from
0 so n_hi = arms*(3^level+1)-1.
=back
=head1 FORMULAS
=cut
# Various formulas for coordinates, boundary, area and more can be found in
# the author's mathematical write-up
#
# =over
#
# L<http://user42.tuxfamily.org/terdragon/index.html>
#
# =back
#
# =head2 N to X,Y
#
lib/Math/PlanePath/ArchimedeanChords.pm view on Meta::CPAN
# is an upper bound, though a fairly slack one
#
#
# cf. arc length along the spiral r=a*theta with a=1/2pi
# arclength = (1/2) * a * (theta*sqrt(1+theta^2) + asinh(theta))
# = (1/4*pi) * (theta*sqrt(1+theta^2) + asinh(theta))
# and theta = 2*pi*r
# = (1/4*pi) * (4*pi^2*r^2 * sqrt(1+1/theta^2) + asinh(theta))
# = pi * (r^2 * sqrt(1+1/r^2) + asinh(theta)/(4*pi^2))
#
# and to compare to the circles formula
#
# = pi * (r*(r+1) * r/(r+1) * sqrt(1+1/r^2)
# + asinh(theta)/(4*pi^2))
#
# so it's smaller hence better upper bound. Only a little smaller than the
# squaring once get to 100 loops or so.
#
#
# not exact
sub rect_to_n_range {
lib/Math/PlanePath/ArchimedeanChords.pm view on Meta::CPAN
total = sum ------------
k=0 arcsin(1/2k)
and this is less than the chords along the spiral. Is there a good
polynomial over-estimate of arcsin, to become an under-estimate total,
without giving away so much?
=head2 Rectangle to N Range
For the C<rect_to_n_range()> upper bound, the current code takes the arc
length along with spiral with the usual formula
arc = 1/4pi * (theta*sqrt(1+theta^2) + asinh(theta))
Written in terms of the r radius (theta = 2pi*r) as calculated from the
biggest of the rectangle x,y corners,
arc = pi*r^2*sqrt(1+1/r^2) + asinh(2pi*r)/4pi
The arc length is longer than chords, so N=ceil(arc) is an upper bound for
the N range.
lib/Math/PlanePath/AztecDiamondRings.pm view on Meta::CPAN
Y<0 X>=0 d=X-Y N=(2d+2)*d+1 + Y
Y<0 X<0 d=X+Y N=(2d+4)*d+2 - Y
For example
Y=2 X=3 d=2+3=5 N=(2*5+2)*5+1 + 2 = 63
Y=2 X=-1 d=2-(-1)=3 N=2*3*3 - 2 = 16
Y=-1 X=4 d=4-(-1)=5 N=(2*5+2)*5+1 + -1 = 60
Y=-2 X=-3 d=-3+(-2)=-5 N=(2*-5+4)*-5+2 - (-2) = 34
The two XE<gt>=0 cases are the same N formula and can be combined with an
abs,
X>=0 d=X+abs(Y) N=(2d+2)*d+1 + Y
This works because at Y=0 the last line of one ring joins up to the start of
the next. For example N=11 to N=15,
15 2
\
14 1
lib/Math/PlanePath/BetaOmega.pm view on Meta::CPAN
----- -------------- -------------
0 0 0
1 0 0 1 = 1
2 -2 = -10 1 = 01
3 -2 = -010 5 = 101
4 -10 = -1010 5 = 0101
5 -10 = -01010 21 = 10101
6 -42 = -101010 21 = 010101
7 -42 = -0101010 85 = 1010101
The power of 4 divided by 3 formulas above for Ymin/Ymax have the effect of
producing alternating bit patterns like this.
For odd levels -Ymin/height approaches 1/3 and Ymax/height approaches 2/3,
ie. the start point is about 1/3 up the total extent. For even levels it's
the other way around, with -Ymin/height approaching 2/3 and Ymax/height
approaching 1/3.
=head2 Closed Curve
Wierum's idea for the curve is a closed square made from four betas,
lib/Math/PlanePath/CCurve.pm view on Meta::CPAN
if ($level == 2) { return (6,0); }
return (2, ($level == 0 ? 0 : 1));
}
my ($h, $rem) = _divrem($level, 2);
my $pow = 2**($h-1);
if ($rem) {
return (6*$pow, 7*$pow-4);
# return (2*S_formula($h) + 2*S_formula($h-1),
# Z_formula($h)/2 + Z_formula($h-1)
# + Z_formula($h)/2 + (S_formula($h)-S_formula($h-1)) );
} else {
return (7*$pow, 3*$pow-4);
# return (S_formula($h) + 2*S_formula($h-1) + S_formula($h)+(Z_formula($h-1)-Z_formula($h-2)),
# (Z_formula($h-1) + Z_formula($h-2)));
}
}
#------------------------------------------------------------------------------
{
my @_UNDOCUMENTED_level_to_hull_area = (0, 1/2, 2);
sub _UNDOCUMENTED_level_to_hull_area {
my ($self, $level) = @_;
if ($level < 3) {
if ($level < 0) { return undef; }
return $_UNDOCUMENTED_level_to_hull_area[$level];
}
my ($h, $rem) = _divrem($level, 2);
return 35*2**($level-4) - ($rem ? 13 : 10)*2**($h-1) + 2;
# if ($rem) {
# return 35*2**($level-4) - 13*$pow + 2;
#
# my $width = S_formula($h) + Z_formula($h)/2 + Z_formula($h-1)/2;
# my $ul = Z_formula($h-1)/2;
# my $ur = Z_formula($h)/2;
# my $bl = $width - Z_formula($h-1)/2 - S_formula($h-1);
# my $br = Z_formula($h-1)/2;
# return $width**2 - $ul**2/2 - $ur**2/2 - $bl**2/2 - $br**2/2;
#
# } else {
# return 35*2**($level-4) - 10*$pow + 2;
# return 0;
# return 35*2**($level-4) - 5*2**$h + 2;
#
# # my $width = S_formula($h) + Z_formula($h-1);
# # my $upper = Z_formula($h-1)/2;
# # my $lower = Z_formula($h-2)/2;
# # my $height = S_formula($h-1) + $upper + $lower;
# # return $width; # * $height - $upper*$upper - $lower*$lower;
# }
# }
}
}
#------------------------------------------------------------------------------
# levels
sub level_to_n_range {
lib/Math/PlanePath/CCurve.pm view on Meta::CPAN
| | | 1
*--* *--* | total
| | | height
* N=4^k N=0 * <-+ -> 1+1/4
| | | | | 1/4
*--*--* *--*--* <-+
^-----^ ^-----^
1/2 1 1/2 total width -> 2
The extent formulas can be found by considering the self-similar blocks.
The initial k=0 is a single line segment and all its extents are 0.
h[0] = 0
N=1 ----- N=0
l[0] = 0
w[0] = 0
Thereafter the replication overlap as
+-------+---+-------+
lib/Math/PlanePath/CCurve.pm view on Meta::CPAN
=over
=item C<($n_lo, $n_hi) = $path-E<gt>level_to_n_range($level)>
Return C<(0, 2**$level)>.
=back
=head1 FORMULAS
Some formulas and results can also be found in the author's mathematical
write-up
=over
L<http://user42.tuxfamily.org/c-curve/index.html>
=back
=head2 Direction
lib/Math/PlanePath/CCurve.pm view on Meta::CPAN
described in L<Math::PlanePath/N to dX,dY -- Fractional>.
# apply turn to make direction at Nint+1
turn = count_low_1_bits(N) - 1 # N integer part
dir = (dir - turn) mod 4 # direction at N+1
# adjust dx,dy by fractional amount in this direction
dx += Nfrac * (dir_to_dx[dir] - dx)
dy += Nfrac * (dir_to_dy[dir] - dy)
A small optimization can be made by working the "-1" of the turn formula
into a +90 degree rotation of the C<dir_to_dx[]> and C<dir_to_dy[]> parts by
swap and sign change,
turn_plus_1 = count_low_1_bits(N) # on N integer part
dir = (dir - turn_plus_1) mod 4 # direction-1 at N+1
# adjustment including extra +90 degrees on dir
dx -= $n*(dir_to_dy[dir] + dx)
dy += $n*(dir_to_dx[dir] - dy)
=head2 X,Y to N
The N values at a given X,Y can be found by taking terms low to high from
the complex number formula (the same as given above)
X+iY = b^k N = 2^k + 2^(k1) + 2^(k2) + ... in binary
+ b^k1 * i base b=1+i
+ b^k2 * i^2
+ ...
If the lowest term is b^0 then X+iY has X+Y odd. If the lowest term is not
b^0 but instead some power b^n then X+iY has X+Y even. This is because a
multiple of b=1+i,
lib/Math/PlanePath/CCurve.pm view on Meta::CPAN
|
-1,0 -- 0,0 -- 1,0
|
0,-1
It's not possible to wait for X=0,Y=0 to be reached because some dX,dY
directions will step infinitely among the four non-zeros. Only the case
X=dX,Y=dY is sure to reach 0,0.
The successive p decrements which rotate dX,dY by -90 degrees must end at p
== 0 mod 4 for highest term in the X+iY formula having i^0=1. This means
must end dX=1,dY=0 East. If this doesn't happen then there is no N for that
p direction.
The number of 1-bits in N is == p mod 4. So the order the N values are
obtained follows the order the p directions are attempted. In general the N
values will not be smallest to biggest N so a little sort is necessary if
that's desired.
It can be seen that sum X+Y is used for the bit calculation and then again
in the divide by 1+i. It's convenient to write the whole loop in terms of
lib/Math/PlanePath/CCurve.pm view on Meta::CPAN
The characteristic polynomial of these recurrences is
x^3 - 4x^2 + 6x - 4
= (x-2) * (x - (1-i)) * (x - (1+i))
=for GP-Test x^3 - 4*x^2 + 6*x - 4 == (x-2)*(x^2 - 2*x + 2)
=for GP-Test x^3 - 4*x^2 + 6*x - 4 == (x-2) * (x + (I-1)) * (x - (I+1))
So explicit formulas can be written in powers of the roots 2, 1-i and 1+i,
M0[k] = ( 2^k + (1-i)^k + (1+i)^k )/4 for k>=1
M1[k] = ( 2^k + i*(1-i)^k - i*(1+i)^k )/4
M2[k] = ( 2^k - (1-i)^k - (1+i)^k )/4
M3[k] = ( 2^k - i*(1-i)^k + i*(1+i)^k )/4
=cut
# GP-DEFINE M0pow(k) = if(k==0,1, (1/4)*(2^k + (1-I)^k + (1+I)^k))
# GP-DEFINE M1pow(k) = if(k==0,0, (1/4)*(2^k + I*(1-I)^k - I*(1+I)^k))
lib/Math/PlanePath/CellularRule190.pm view on Meta::CPAN
14 15 16 17 18 19 20 4
8 9 10 11 12 13 3
4 5 6 7 2
1 2 3 1
0 <- Y=0
-5 -4 -3 -2 -1 X=0 1 2 3 4 5
The effect is to push each N rightwards by 1, and wrapping around. So the
N=0,1,4,8,14,etc on the left were on the right of the default n_start=1.
This also has the effect of removing the +1 in the Nright formula given
above, so
Nleft = triangular(Y) + quartersquare(Y)
=head1 FUNCTIONS
See L<Math::PlanePath/FUNCTIONS> for behaviour common to all path classes.
=over 4
lib/Math/PlanePath/ChanTree.pm view on Meta::CPAN
6 | 4 7
5 | 3 8 N=6 is the 7/7 middle term
4 | 2 9
3 | 1 10
2 | 0 11
1 | 12
Y=0 |
+------------------------------
X=0 1 2 3 4 5 6 7
Each node has k children. The formulas for the children can be seen from
sample cases k=5 and k=6. A node X/Y descends to
k=5 k=6
1X+0Y / 2X+1Y 1X+0Y / 2X+1Y
2X+1Y / 3X+2Y 2X+1Y / 3X+2Y
3X+2Y / 2X+3Y 3X+2Y / 3X+3Y
2X+3Y / 1X+2Y 3X+3Y / 2X+3Y
1X+2Y / 0X+1Y 2X+3Y / 1X+2Y
1X+2Y / 0X+1Y
lib/Math/PlanePath/ChanTree.pm view on Meta::CPAN
parent = floor(N/k) when n_start=1
For other C<n_start> adjust before and after to an C<n_start=1> basis,
parent = floor((N-(Nstart-1)) / k) + Nstart-1
For example in the default k=0 Nstart=1 the parent of N=3 is
floor((3-(1-1))/3)=1.
The post-adjustment can be worked into the formula with (k-1)*(Nstart-1)
similar to the children above,
parent = floor((N + (k-1)*(Nstart-1)) / k)
But the first style is more convenient to compare to see that N is past the
top nodes and therefore has a parent.
N-(Nstart-1) >= k to check N is past top-nodes
=head2 N Root
lib/Math/PlanePath/ComplexMinus.pm view on Meta::CPAN
=item C<($n_lo, $n_hi) = $path-E<gt>level_to_n_range($level)>
Return C<(0, 2**$level - 1)>, or with C<realpart> option return C<(0,
$norm**$level - 1)> where norm=realpart^2+1.
=back
=head1 FORMULAS
Various formulas and pictures etc for the i-1 case can be found in the
author's long mathematical write-up (section "Complex Base i-1")
=over
L<http://user42.tuxfamily.org/dragon/index.html>
=back
=head2 X,Y to N
A given X,Y representing X+Yi can be turned into digits of N by successive
complex divisions by i-r. Each digit of N is a real remainder 0 to r*r
inclusive from that division.
The base formula above is
X+Yi = a[n]*b^n + ... + a[2]*b^2 + a[1]*b + a[0]
and want the a[0]=digit to be a real 0 to r*r inclusive. Subtracting a[0]
and dividing by b will give
(X+Yi - digit) / (i-r)
= - (X-digit + Y*i) * (i+r) / norm
= (Y - (X-digit)*r)/norm
+ i * - ((X-digit) + Y*r)/norm
lib/Math/PlanePath/ComplexMinus.pm view on Meta::CPAN
first norm^k many points, ie. N=0 to N=norm^k-1 inclusive, is calculated in
=over
William J. Gilbert, "The Fractal Dimension of Sets Derived From Complex
Bases", Canadian Mathematical Bulletin, volume 29, number 4, 1986.
L<http://www.math.uwaterloo.ca/~wgilbert/Research/GilbertFracDim.pdf>
=back
The boundary formula is a 3rd-order recurrence. For the twindragon case it
is
for realpart=1
boundary[k] = boundary[k-1] + 2*boundary[k-3]
= 4, 6, 10, 18, 30, 50, 86, 146, 246, 418, ... (2*A003476)
4 + 2*x + 4*x^2
generating function ---------------
1 - x - 2*x^3
lib/Math/PlanePath/ComplexMinus.pm view on Meta::CPAN
C -> B
starting from
A = 2*r
B = 2
C = 2 - 2*r
total boundary = A+B+C
For i-1 realpart=1 these A,B,C are already in the form of a recurrence
A-E<gt>A+2*C, B-E<gt>A, C-E<gt>B, per the formula above. For other real
parts a little matrix rearrangement turns the A,B,C parts into recurrence
boundary[k] = boundary[k-1] * (2*r - 1)
+ boundary[k-2] * (norm - 2*r)
+ boundary[k-3] * norm
starting from
boundary[0] = 4 # single square cell
boundary[1] = 2*norm + 2 # oblong of norm many cells
boundary[2] = 2*(norm-1)*(r+2) + 4
lib/Math/PlanePath/ComplexPlus.pm view on Meta::CPAN
=item C<($n_lo, $n_hi) = $path-E<gt>level_to_n_range($level)>
Return C<(0, 2**$level - 1)>, or for 2 arms return C<(0, 2 * 2**$level -
1)>. With the C<realpart> option return C<(0, $norm**$level - 1)> where
norm=realpart^2+1.
=back
=head1 FORMULAS
Various formulas and pictures etc for the i+1 case can be found in the
author's long mathematical write-up (section "Complex Base i+1")
=over
L<http://user42.tuxfamily.org/dragon/index.html>
=back
=head1 OEIS
lib/Math/PlanePath/Corner.pm view on Meta::CPAN
Then the X,Y coordinates are
if (Off < 0) then X=d+Off, Y=d
if (Off >= 0) then X=d, Y=d-Off
=head2 X,Y to N
For a given X,Y the bigger of X or Y determines the d gnomon.
If YE<gt>=X then X,Y is on the horizontal part. At X=0 have N=StartN(d) per
the Start(N) formula above, and any further X is an offset from there.
if Y >= X then
d=Y
N = StartN(d) + X
= Y^2 + 1 + X
Otherwise if YE<lt>X then X,Y is on the vertical part. At Y=0 N is the last
point on the gnomon, and one back from the start of the following gnomon,
if Y <= X then
lib/Math/PlanePath/Diagonals.pm view on Meta::CPAN
X^2 + 3X + 2XY + Y + Y^2
N = ------------------------ + Nstart
2
(X+Y)^2 + 3X + Y
= ---------------- + Nstart (using one square)
2
=head2 N to X,Y
The above formula N=d*(d+1)/2 can be solved for d as
d = floor( (sqrt(8*N+1) - 1)/2 )
# with n_start=0
For example N=12 is d=floor((sqrt(8*12+1)-1)/2)=4 as that N falls in the
fifth diagonal. Then the offset from the Y axis NY=d*(d-1)/2 is the X
position,
X = N - d*(d+1)/2
Y = X - d
In the code, fractional N is handled by imagining each diagonal beginning
0.5 back from the Y axis. That's handled by adding 0.5 into the sqrt, which
is +4 onto the 8*N.
d = floor( (sqrt(8*N+5) - 1)/2 )
# N>=-0.5
The X and Y formulas are unchanged, since N=d*(d-1)/2 is still the Y axis.
But each diagonal d begins up to 0.5 before that and therefore X extends
back to -0.5.
=cut
# GP-DEFINE \\ diagonal length
# GP-DEFINE d(n) = floor( (sqrt(8*n+1) - 1)/2 );
# GP-Test vector(16,n,n--; d(n)) == \
# GP-Test [0, 1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 5]
# GP-DEFINE X(n) = my(d=d(n)); n - d*(d+1)/2;
lib/Math/PlanePath/DiagonalsOctant.pm view on Meta::CPAN
}
$d = int( _sqrtint($d)/2 );
### $d
# remainder positive or negative relative to the start of the following
# diagonal
#
$n -= $d*($d+1) + 1;
### remainder: $n
# $n first in formulas to preserve n=BigFloat when d=integer is BigInt
#
if ($self->{'direction'} eq 'up') {
if (2*$n >= -1) {
return (-$n + $d,
$n + $d);
} else {
return (-$n - 1,
$n + 2*$d);
}
} else {
lib/Math/PlanePath/DiagonalsOctant.pm view on Meta::CPAN
d = floor sqrt(N-0.5)
= int( sqrt(int(4*$n)-2)/2 )
Taking /2 out of the sqrt helps with C<Math::BigInt> which circa Perl 5.14
doesn't inter-operate with flonums very well.
In any case N-Nstart is an offset into two diagonals, the first of length d
many points and the second d+1. For example d=3 starting Y=5 for points
N=10,11,12 followed by Y=6 N=13,14,15,16.
The formulas are simplified by calculating a remainder relative to the
second diagonal, so it's negative for the first and positive for the second,
Nrem = N - (d*(d+1)+1)
d*(d+1)+1 is 1 past the pronic numbers when end each first diagonal, as
described above. In any case for example d=3 is relative to N=13 making
Nrem=-3,-2,-1 or Nrem=0,1,2,3.
To include the preceding 0.5 in the second diagonal simply means reckoning
NremE<gt>=-0.5 as belonging to the second. In that base
lib/Math/PlanePath/DivisibleColumns.pm view on Meta::CPAN
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
=head1 FORMULAS
=head2 Rectangle to N Range
The cumulative divisor count up to and including a given X column can be
calculated from the fairly well-known sqrt formula, a sum from 1 to sqrt(X).
S = floor(sqrt(X))
/ i=S \
numdivs cumulative = 2 * | sum floor(X/i) | - S^2
\ i=1 /
This means the N range for 0 to X can be calculated without working out all
each column count up to X. In the current code if column counts have been
worked out then they're used, otherwise this formula.
=head1 OEIS
This pattern is in Sloane's Online Encyclopedia of Integer Sequences in the
following forms,
=over
L<http://oeis.org/A061017> (etc)
lib/Math/PlanePath/DragonCurve.pm view on Meta::CPAN
10--11,7---6 3---2
| |
16 13---12 0---1
| |
15---14 <- Wmin = -(2^2 - 1)/3 = -1
^ ^Lmin = Wmin = -1
|
Lmax = (7*2^2 - 4)/6 = 4
The formulas are all integer values, but the fractions 7/6, 1/3 and 2/3 show
the limits as the level increases. If scaled so that length Lend=2^k is
reckoned as 1 unit then Lmax extends 1/6 past the end, Lmin and Wmin extend
1/3, and Wmax extends across 2/3.
+--------+ --
| - | 1/6 total length
|| | | = 1/6+1+1/3 = 3/2
|| E | --
|| |
|| |
lib/Math/PlanePath/DragonCurve.pm view on Meta::CPAN
2**$level + ($arms-1))>.
There are 2^level segments in a curve level, so 2^level+1 points numbered
from 0. For multiple arms there are arms*(2^level+1) points, numbered from
0 so n_hi = arms*(2^level+1)-1.
=back
=head1 FORMULAS
Various formulas for coordinates, lengths and area can be found in the
author's long mathematical write-up
=over
L<http://user42.tuxfamily.org/dragon/index.html>
=back
=cut
lib/Math/PlanePath/GcdRationals.pm view on Meta::CPAN
"rows_reverse"
"diagonals_down"
"diagonals_up"
=back
=head1 FORMULAS
=head2 X,Y to N -- Rows
The defining formula above for X,Y can be inverted to give i,j and N. This
calculation doesn't notice if X,Y have a common factor, so a coprime(X,Y)
test must be made separately if necessary (for C<xy_to_n()> it is).
X/Y = g-1 + (i/g)/(j/g)
The g-1 integer part is recovered by a division X divide Y,
X = quot*Y + rem division by Y rounded towards 0
where 0 <= rem < Y
unless Y=1 in which case use quot=X-1, rem=1
g-1 = quot
g = quot+1
The Y=1 special case can instead be left as the usual kind of division
quot=X,rem=0, so 0E<lt>=remE<lt>Y. This will give i=0 which is outside the
intended 1E<lt>=iE<lt>=j range, but j is 1 bigger and the combination still
gives the correct N. It's as if the i=g,j=g point at the end of a row is
moved to i=0,j=g+1 just before the start of the next row. If only N is of
interest not the i,j then it can be left rem=0.
Equating the denominators in the X/Y formula above gives j by
Y = j/g the definition above
j = g*Y
= (quot+1)*Y
j = X+Y-rem per the division X=quot*Y+rem
And equating the numerators gives i by
X = (g-1)*Y + i/g the definition above
lib/Math/PlanePath/GcdRationals.pm view on Meta::CPAN
g = ceil(X/Y)
= cquot for division X=cquot*Y - crem
But division in most programming languages is towards 0 or towards
-infinity, not upwards towards +infinity.
=head2 X,Y to N -- Rows Reverse
For pairs_order="rows_reverse", the horizontal i is reversed to j-i+1. This
can be worked into the triangular part of the N formula as
Nrrev = (j-i+1) + j*(j-1)/2 for 1<=i<=j
= j*(j+1)/2 - i + 1
The Y=1 case described above cannot be left to go through with rem=0 giving
i=0 and j+1 since the reversal j-i+1 is then not correct. Either use rem=1
as described, or if not then compensate at the end,
if r=0 then j -= 2 adjust
Nrrev = j*(j+1)/2 - i + 1 same Nrrev as above
lib/Math/PlanePath/GcdRationals.pm view on Meta::CPAN
j-=2 so that j=6-2=4 gives the desired Nrrev=4*5/2-0+1=11 (per the table in
L</Rows Reverse> above).
=cut
# No, not quite
#
# =head2 Rectangle N Range -- Rows
#
# An over-estimate of the N range can be calculated just from the X,Y to N
# formula above.
#
# Within a row N increases with increasing X, so for a rectangle the minimum
# is in the left column and the maximum in the right column.
#
# Within a column N values increase until reaching the end of a "g" wedge,
# then drop down a bit. So the maximum is either the top-right corner of the
# rectangle, or the top of the next lower wedge, ie. smaller Y but bigger g.
# Conversely the minimum is either the bottom right of the rectangle, or the
# bottom of the next higher wedge, ie. smaller g but bigger Y. (Is that
# right?)
( run in 0.853 second using v1.01-cache-2.11-cpan-26ccb49234f )