Math-PlanePath
view release on metacpan or search on metacpan
devel/c-curve.pl view on Meta::CPAN
# 3------2------1
# h=0 n=2*4^h=2 boundary=4
# h=1 n=2*4^h=8 boundary=16 area=12
# h=2 n=2*4^h=32 boundary=40
# h=3 n=2*4^h=128 boundary=88
# total 4*(2^h) + 4*(2*2^h - 2)
# = 3*2^h-8
# A182461 whole except 4 a(n) = a(n-1)*2+8 16,40,88,
# A131128 3*2^n - 4 half
# A033484 3*2^n - 2 quarter
# A153893 3*2^n - 1 eighth h>=1
require MyOEIS;
my @values;
foreach my $k (0 .. 16) {
my $n_end = 2**$k;
my $h = int($k/2);
my ($n1, $n2) = ($k % 2 ? diagonal_4k_axis_n_ends($h) : width_4k_axis_n_ends($h));
my ($x1,$y1) = $path->n_to_xy ($n1);
my ($x2,$y2) = $path->n_to_xy ($n2);
my $points = MyOEIS::path_boundary_points_ft($path, $n_end,
$x1,$y1, $x2,$y2,
side => 'right',
dir => $h,
);
if (@$points < 30) {
print "k=$k from N=$n1 $x1,$y1 to N=$n2 $x2,$y2\n";
print " ",points_str($points),"\n";
}
my $boundary = 2 * (scalar(@$points) - 1);
my $area;
if (@$points > 2) {
my $planar = Math::Geometry::Planar->new;
$planar->points($points);
$area = 2 * $planar->area;
} else {
$area = 0;
}
# push @values, $boundary;
push @values, $area;
print "$h B=$boundary A=$area n=$n1 xy=$x1,$y1 to n=$n2 xy=$x2,$y2 limit $n_end\n";
}
print join(',',@values),"\n";
shift @values;
shift @values;
Math::OEIS::Grep->search(array => \@values);
exit 0;
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;
}
{
# triangle area parts by individual recurrence
# e[k+8] = e[k+7] + 2*e[k+6] - e[k+4] + e[k+3] + 2*e[k+1] + 4*e[k]
# 1,0,0,0,0,0,0,2,6,10,22,40,80,156,308,622,1242,2494,4994,9988,19988,39952,79904,159786
my @e = (1,0,0,0,0,0,0,2);
foreach my $k (0 .. 15) {
push @e, $e[-1] + 2*$e[-2] - $e[-4] + $e[-5] + 2*$e[-7] + 4*$e[-8];
}
print join(",",@e),"\n";
exit 0;
}
{
# area parts by a-z recurrence
#
# a 0,0,0,2,4,8,16,30,60,116,232,466,932,1872,3744,7494
# c 0,1,1,1,1,2,4,8,18,39,79,159,315,628,1250,2494
# e 1,0,0,0,0,0,0,2,6,10,22,40,80,156,308,622
# g 0,0,0,0,0,0,0,2,2,6,10,20,40,76,156,310
# b 0,0,1,1,3,5,10,20,38,78,155,311,625,1247,2500,4994
# d 0,0,0,1,1,3,5,10,20,38,78,155,311,625,1247,2500
# f 0,0,0,0,1,1,3,5,10,20,38,78,155,311,625,1247
# h 0,0,0,0,0,1,1,3,5,10,20,38,78,155,311,625
# i 0,0,0,0,0,0,1,1,3,5,10,20,38,78,155,311
#
# [4]
# [2]
# [0]
# [1]
# [-1]
# [0]
# [2]
# [1]
# x^8 - (4*x^7 + 2*x^6 + 0*x^5 + 1*x^4 + -1*x^3 + 0*x^2 + 2*x + 1)
#
# 2,6,10,22,40,80,156,308,622,1242,2494,4994,9988,19988,39952,79904,159786,319550,639122,1278222,2556512,5113048,10226116,20452300,40904486
#
# 4*2 + 2*6 + 0*10 + 22 - 40 + 0*80 + 2*156 + 308
# a*x^2*g(x) + b*x*g(x) - g(x) = initial
# (-2 - 4*x)/(-1 + 1*x + 2*x^2 + 0*x^3 - x^4 + x^5 + 0*x^6 + 2*x^7 + 4*x^8)
#
#
my (@a,@b,@c,@d,@e,@f,@g,@h,@i);
my $a = 0;
my $b = 0;
my $c = 0;
devel/c-curve.pl view on Meta::CPAN
if ($k < 10) { return $want[$k]; }
return ($B->($k-4) * 2
+ $B->($k-3) * -4
+ $B->($k-2) * 1
+ $B->($k-1) * 2);
};
# $B = sub {
# my ($k) = @_;
# return MyOEIS::path_boundary_length($path, 2**$k);
# };
$|=1;
foreach my $k (0 .. $#want) {
my $want = $want[$k];
my $got = $B->($k);
my $diff = $want - $got;
print "$k $want $got $diff\n";
# print "$got,";
}
exit 0;
}
{
# left boundary vs recurrence
# L[k] = 4*L[k-1] - 5*L[k-2] + 2*L[k-3] k >= 6
# x^3 - 4*x^2 + 5*x - 2 = (x-1)^2 * (x-2) so a*2^k + b*k + c
#
# explicit L[2*h] = 55/4 * 2^h + 28*h - 130 # h>=3
# explicit L[2*h+1] = 78/4 * 2^h + 28*h - 116
# my @want = (1,4,16,64,202,450,918,1826,3614,7162,14230,28338,56526,112874,225542,450850,901438,1802586);
my @want = (2,8,32,124,308,648,1300,2576,5100,10120,20132,); # left 2k+1
my $L;
$L = sub {
my ($k) = @_;
if ($k < 6) { return $want[$k]; }
return ($L->($k-3) * 2
+ $L->($k-2) * -5
+ $L->($k-1) * 4);
};
$L = sub {
my ($k) = @_;
if ($k < 3) { return $want[$k]; }
return 78/4*2**$k + 28*$k - 116;
};
# $L = sub {
# my ($k) = @_;
# return MyOEIS::path_boundary_length($path, 2*4**$k, side => 'left');
# };
$|=1;
foreach my $k (0 .. $#want) {
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
return 10*2**$h - 4*$h - 8; # yes
return 2*2**$h + (2*$k-2)*(2**$h-1) - 4*($h-2)*2**$h - 8; # yes
{
my $r = 0;
foreach my $i (1 .. $h-1) { # yes
$r += $i * 2**$i;
}
return 2*2**$h + (2*$k-2)*(2**$h-1) - 4*$r;
}
{
my $r = 0;
foreach my $i (0 .. $h-1) {
$r += (2*$k-2 - 4*$i) * 2**$i;
}
return 2*2**$h + $r
}
{
my $r = 0;
my $pow = 1;
while ($k >= 3) {
### t: 2*$k-2
$r += (2*$k-2) * $pow;
$pow *= 2;
$k -= 2;
}
return $r + 2*$pow;
}
} else {
my $h = $k/2;
{
return 7*2**$h - 4*$h - 6; # yes
return (2*$k-1) * 2**$h - 2*$k + 2 - 4*(($h-1-1)*2**($h-1+1) + 2);
}
{
# right[k] = 2k-2 + 2*right[k-2] termwise, yes
my $r = 0;
foreach my $i (0 .. $h-1) {
$r += $i*2**$i;
}
return (2*$k-1) * 2**$h - 2*$k + 2 - 4*$r;
}
{
# right[k] = 2k-2 + 2*right[k-2] termwise, yes
my $r = 0;
my $pow = 1;
while ($k > 0) {
$r += (2*$k-2) * $pow;
$pow *= 2;
$k -= 2;
}
return $r + $pow;
}
return ($h-2) *2**$h;
}
};
my ($Scum_recurrence, $Zcum_recurrence, $R_recurrence);
$Scum_recurrence = sub {
my ($k) = @_;
if ($k == 0) { return 0; }
if ($k == 1) { return 0; }
return 2*$Scum_recurrence->($k-2) + $k-1; # yes
return $Zcum_recurrence->($k-1); # yes
};
$Zcum_recurrence = sub {
my ($k) = @_;
if ($k == 0) { return 0; }
if ($k == 1) { return 1; }
return 2*$Zcum_recurrence->($k-2) + $k; # yes
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;
my ($dx,$dy) = (1,0);
my ($target_x,$target_y) = $path->n_to_xy($n_end);
until ($x == $target_x && $y == $target_y) {
### at: "$x, $y $dx,$dy"
($dx,$dy) = ($dy,-$dx); # rotate -90
if (path_xy_is_visited_within ($path, $x+$dx,$y+$dy, $n_end)) {
$z++;
} else {
($dx,$dy) = (-$dy,$dx); # rotate +90
if (path_xy_is_visited_within ($path, $x+$dx,$y+$dy, $n_end)) {
$s++;
} else {
($dx,$dy) = (-$dy,$dx); # rotate +90
$z++;
path_xy_is_visited_within ($path, $x+$dx,$y+$dy, $n_end) or die;
}
}
$x += $dx;
$y += $dy;
}
return ($s, $z);
}
sub path_xy_is_visited_within {
my ($path, $x,$y, $n_end) = @_;
my @n_list = $path->xy_to_n_list($x,$y);
foreach my $n (@n_list) {
if ($n <= $n_end) {
return 1;
}
}
return 0;
}
}
{
# diagonal N endpoints search
my @values;
foreach my $k (0 .. 10) {
my ($n1, $n2) = diagonal_4k_axis_n_ends($k);
my ($x1,$y1) = $path->n_to_xy ($n1);
my ($x2,$y2) = $path->n_to_xy ($n2);
foreach (1 .. $k) {
($x1,$y1) = ($y1,-$x1); # rotate -90
($x2,$y2) = ($y2,-$x2); # rotate -90
}
push @values, $n2;
printf "$n1 xy=$x1,$y1 $n2 xy=$x2,$y2 %b %b\n", $n1, $n2;
( run in 0.975 second using v1.01-cache-2.11-cpan-df04353d9ac )