#!/usr/bin/perl -w
# Copyright 2011, 2012, 2013, 2014 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/>.
use 5.010;
use strict;
use warnings;
use Math::PlanePath::KochSquareflakes;
# uncomment this to run the ### lines
# use Smart::Comments;
{
# area
#
# start diag[0] = 0
# start straight[0] = 4
# diag[n+1] = 2*straight[n] + 2*diag[n]
# straight[n+1] = 2*straight[n] + 2*diag[n]
#
#
require Math::Geometry::Planar;
my $path = Math::PlanePath::KochSquareflakes->new;
my @values;
my $prev_a_log = 0;
my $prev_len_log = 0;
foreach my $level (1 .. 7) {
my $n_level = (4**($level+1) - 1) / 3;
my $n_end = $n_level + 4**$level;
my @points;
foreach my $n ($n_level .. $n_end) {
my ($x,$y) = $path->n_to_xy($n);
push @points, [$x,$y];
}
### @points
my $polygon = Math::Geometry::Planar->new;
$polygon->points(\@points);
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;
use lib 'xt'; require MyOEIS;
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
#
# diag: 1, 3, 10, 34, 116, 396, 1352, 4616, 15760, 53808, 183712, 627232
# A007052 a(n+1) = 4*a(n) - 2*a(n-1), starting 1,3
# A007070 max horiz dist from ring start pos 4,14,48,164 side width
# A206374 N of the max position 2,9,37,149 corner
# A003480 X of the max position 2,7,24,82 last
# A007052 max vert dist from ring start pos 3,10,34,116 height
# A072261 N of the max Y position 7,29,117,469 Y axis
# A007052 Y of the max position 3,10,34,116
my $path = Math::PlanePath::KochSquareflakes->new;
my @values;
my $coord = 1;
foreach my $level (1 .. 8) {
my $nstart = (4**($level+1) - 1) / 3;
my $nend = $nstart + 4**$level;
my @start = $path->n_to_xy($nstart);
my $max_offset = 0;
my $max_offset_n = $nstart;
my $max_offset_c = $start[$coord];
foreach my $n ($nstart .. $nend) {
my @this = $path->n_to_xy($n);
my $offset = abs ($this[$coord] - $start[$coord]);
if ($offset > $max_offset) {
$max_offset = $offset;
$max_offset_n = $n;
$max_offset_c = $this[$coord];
}
}
push @values, $max_offset;
print "level $level start=$start[$coord] max offset $max_offset at N=$max_offset_n (of $nstart to $nend) Y=$max_offset_c\n";
}
use lib 'xt'; require MyOEIS;
Math::OEIS::Grep->search(array => \@values);
exit 0;
}
{
# X or Y coordinate of first point of ring
# X or Y coord: 1, 2,7,24,82,280,
# A003480 1,2,7 OFFSET=0
# A020727 2,7
#
# cf A006012 same recurrence, start 1,2
my $path = Math::PlanePath::KochSquareflakes->new;
my @values;
foreach my $level (1 .. 12) {
my $nstart = (4**($level+1) - 1) / 3;
my ($x,$y) = $path->n_to_xy($nstart);
push @values, -$y;
}
use lib 'xt'; require MyOEIS;
Math::OEIS::Grep->search(array => \@values);
exit 0;
}
{
# Xstart power
# Xstart = b^level
# b = Xstart^(1/level)
#
# D = P^2-4Q = 4^2-4*-2 = 24
# sqrt(24) = 4.898979485566356196394568149
#
my $path = Math::PlanePath::KochSquareflakes->new;
my $prev = 1;
foreach my $level (1 .. 12) {
my $nstart = (4**($level+1) - 1) / 3;
my ($xstart,$ystart) = $path->n_to_xy($nstart);
$xstart = -$xstart;
my $f = $xstart / $prev;
# my $b = $xstart ** (1/($level+1));
print "level=$level xstart=$xstart f=$f\n";
$prev = $xstart;
}
print "\n";
exit 0;
}
{
my @horiz = (1);
my @diag = (1);
foreach my $i (0 .. 10) {
$horiz[$i+1] = 2*$horiz[$i] + 2*$diag[$i];
$diag[$i+1] = $horiz[$i] + 2*$diag[$i];
$i++;
}
print "horiz: ",join(', ',@horiz),"\n";
print "diag: ",join(', ',@diag),"\n";
exit 0;
}