#!/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::KochSnowflakes;
# uncomment this to run the ### lines
# use Smart::Comments;
{
# num angles taking points pairwise
#
# one direction, so modulo 180
# can go outside: 3,12,102,906,7980
# not going outside: 3,12,90,666
#
# both directions
# can go outside: 3,18,155,1370
# not going outside: 3,18,124,917,6445
# 12, one direction always X increasing
# 0,-1 12->6 south
# 1,-3 12->7
# 1,-1 12->11
# 2,-1 14->7
# 3,-1 12->10
# 5,-1 14->9
# 1,0 14->13 east
# 5,1 3,1 2,1 1,1 1,3
#
# 12
# / \
# 14----13 11----10
# \ /
# 15 9
# \
# 4---- 5 7---- 8
# \ /
# 6
require Math::Geometry::Planar;
my $path = Math::PlanePath::KochSnowflakes->new;
my @values;
require Math::PlanePath::GcdRationals;
foreach my $level (0 .. 7) {
my $n_level = 4**$level;
my $n_end = 4**($level+1) - 1;
my @points;
foreach my $n ($n_level .. $n_end) {
my ($x,$y) = $path->n_to_xy($n);
push @points, [$x,$y];
}
sub points_equal {
my ($p1, $p2) = @_;
return (abs($p1->[0] - $p2->[0]) < 0.00001
&& abs($p1->[1] - $p2->[1]) < 0.00001);
}
# Return true if line $p1--$p2 is entirely within the snowflake.
# If $p1--$p2 intersects any boundary segment then it's not entirely
# within.
#
# -3,-1 ---- -1,-1
my $segment_inside = sub {
my ($p1,$p2) = @_;
### segment_inside: join(',',@$p1).' '.join(',',@$p2)
foreach my $i (0 .. $#points) {
my $pint = Math::Geometry::Planar::SegmentIntersection([$p1,$p2, $points[$i-1],$points[$i]]) || next;
### intersects: join(',',@{$points[$i-1]}).' '.join(',',@{$points[$i]}).' at '.join(',',@$pint)
if (points_equal($pint,$p1) || points_equal($pint,$p1)
|| points_equal($pint,$points[$i-1]) || points_equal($pint,$points[$i])) {
### is an endpoint ...
next;
}
if (Math::Geometry::Planar::Colinear([$p1,$pint,$points[$i]])) {
### but colinear ...
next;
}
### no, crosses boundary ...
return 0;
}
### yes, all inside ...
return 1;
};
my %seen;
foreach my $i (0 .. $#points) {
foreach my $j ($i+1 .. $#points) {
### pair: "$i -- $j"
my $p1 = $points[$i];
my $p2 = $points[$j];
my $dx = $p1->[0] - $p2->[0];
my $dy = $p1->[1] - $p2->[1];
my $rdx = $dx;
my $rdy = $dy;
if ($rdx < 0 || ($rdx == 0 && $rdy < 0)) {
$rdx = -$rdx;
$rdy = -$rdy;
}
my $g = Math::PlanePath::GcdRationals::_gcd($rdx,$rdy);
$rdx /= $g;
$rdy /= $g;
next if $seen{"$rdx,$rdy"};
next unless $segment_inside->($p1, $p2);
$seen{"$rdx,$rdy"} = 1;
}
}
my $count = scalar(keys %seen);
if ($count <= 12) {
my @strs = keys %seen;
@strs = sort strpoint_cmp_by_atan2 @strs;
print join(' ',@strs),"\n";
}
push @values, $count;
print "$level $count\n";
}
require Math::OEIS::Grep;
Math::OEIS::Grep->search(array => \@values);
exit 0;
sub strpoint_cmp_by_atan2 {
my ($ax,$ay) = split /,/, $a;
my ($bx,$by) = split /,/, $b;
return atan2($ay,$ax) <=> atan2($by,$bx);
}
}
{
# num acute angle turns, being num right hand turns
# A178789
my $path = Math::PlanePath::KochSnowflakes->new;
my @values;
require Math::NumSeq::PlanePathTurn;
foreach my $level (0 .. 8) {
my $n_level = 4**$level;
my $n_end = 4**($level+1) - 1;
my @x;
my @y;
foreach my $n ($n_level .. $n_end) {
my ($x,$y) = $path->n_to_xy($n);
push @x, $x;
push @y, $y;
}
my $count = 0;
foreach my $i (0 .. $#x) {
my $dx = $x[$i-1] - $x[$i-2];
my $dy = $y[$i-1] - $y[$i-2];
my $next_dx = $x[$i] - $x[$i-1];
my $next_dy = $y[$i] - $y[$i-1];
my $tturn6 = Math::NumSeq::PlanePathTurn::_turn_func_TTurn6($dx,$dy, $next_dx,$next_dy);
### $tturn6
if ($tturn6 == 2 || $tturn6 == 4) {
$count++;
}
}
$count = ($n_end - $n_level + 1) - $count; # non-acute ones
push @values, $count;
my $calc = 4**$level + 2;
print "$level $count $calc\n";
}
require Math::OEIS::Grep;
Math::OEIS::Grep->search(array => \@values);
exit 0;
}
{
# num lines at 60deg angles, A110593 2*3^n
# 3,6,18,54,162
# 18*3=54
#
my $path = Math::PlanePath::KochSnowflakes->new;
my @values;
require Math::NumSeq::PlanePathDelta;
foreach my $level (0 .. 7) {
my $n_level = 4**$level;
my $n_end = 4**($level+1) - 1;
my %seen;
foreach my $n ($n_level .. $n_end) {
my ($x,$y) = $path->n_to_xy($n);
my ($dx,$dy) = $path->n_to_dxdy($n);
if ($n == $n_end) {
my ($x2,$y2) = $path->n_to_xy($n_level);
$dx = $x2 - $x;
$dy = $y2 - $y;
}
my $tdir6 = Math::NumSeq::PlanePathDelta::_delta_func_TDir6($dx,$dy);
my $value;
if ($tdir6 % 3 == 0) {
$value = "H-$y";
} elsif ($tdir6 % 3 == 1) {
$value = "I-".($x-$y);
} else {
$value = "J-".($x+$y);
}
# print " $n $value\n";
$seen{$value} = 1;
}
my $count = scalar(keys %seen);
push @values, $count;
print "$level $count\n";
}
require Math::OEIS::Grep;
Math::OEIS::Grep->search(array => \@values);
exit 0;
}
{
# area
require Math::Geometry::Planar;
my $path = Math::PlanePath::KochSnowflakes->new;
my @values;
my $prev_a_log = 0;
my $prev_len_log = 0;
# area of an equilateral triangle of side length 1
use constant AREA_EQUILATERAL_1 => sqrt(3)/4;
foreach my $level (0 .. 7) {
my $n_level = 4**$level;
my $n_end = 4**($level+1) - 1;
my @points;
foreach my $n ($n_level .. $n_end) {
my ($x,$y) = $path->n_to_xy($n);
# ($x,$y) = triangular_to_unit_side($x,$y);
push @points, [$x,$y];
}
my $polygon = Math::Geometry::Planar->new;
$polygon->points(\@points);
my $a = $polygon->area;
my $len = $polygon->perimeter;
# $a /= 3**$len;
# $len /= sqrt(3)**$len;
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;
# $a /= AREA_EQUILATERAL_1();
print "$level $a\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 triangular_to_unit_side {
my ($x,$y) = @_;
return ($x/2, $y*(sqrt(3)/2));
}
}
{
# area
print sqrt(48)/10,"\n";
my $sum = 0;
foreach my $level (1 .. 10) {
my $area = (1 + 3/9*$sum) * sqrt(3)/4;
print "$level $area\n";
$sum += (4/9)**$level;
}
exit 0;
}
{
# X axis N increasing
my $path = Math::PlanePath::KochSnowflakes->new;
my $prev_n = 0;
foreach my $x (0 .. 10000000) {
my $n = $path->xy_to_n($x,0) // next;
if ($n < $prev_n) {
print "decrease N at X=$x N=$n prev_N=$prev_n\n";
}
$prev_n = $n;
}
}