#!/usr/bin/perl -w
# Copyright 2012 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.004;
use strict;
use List::Util 'min', 'max';
# uncomment this to run the ### lines
use Smart::Comments;
{
# min/max for level
$|=1;
require Math::PlanePath::ComplexRevolving;
my $path = Math::PlanePath::ComplexRevolving->new;
my $prev_max = 1;
my @min = (1);
for (my $level = 1; $level < 25; $level++) {
my $n_start = 2**($level-1);
my $n_end = 2**$level;
my $min_hypot = 128*$n_end*$n_end;
my $min_x = 0;
my $min_y = 0;
my $min_pos = '';
my $max_hypot = 0;
my $max_x = 0;
my $max_y = 0;
my $max_pos = '';
# print "level $level n=$n_start .. $n_end\n";
foreach my $n ($n_start .. $n_end) {
my ($x,$y) = $path->n_to_xy($n);
my $h = $x*$x + $y*$y;
if ($h < $min_hypot) {
$min_hypot = $h;
$min_pos = "$x,$y";
}
if ($h > $max_hypot) {
$max_hypot = $h;
$max_pos = "$x,$y";
}
}
# print "$min_hypot,";
$min[$level] = $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 / $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;
}
{
require Math::PlanePath::ComplexRevolving;
require Image::Base::Text;
my $realpart = 2;
my $radix = $realpart*$realpart + 1;
my %seen;
my $isize = 20;
my $image = Image::Base::Text->new (-width => 2*$isize+1,
-height => 2*$isize+1);
foreach my $n (0 .. $radix**6) {
my $x = 0;
my $y = 0;
my $bx = 1;
my $by = 0;
foreach my $digit (digits($n,$radix)) {
if ($digit) {
$x += $digit * $bx;
$y += $digit * $by;
($bx,$by) = (-$by,$bx); # (bx+by*i)*i = bx*i - by, rotate +90
}
# (bx,by) = (bx + i*by)*(i+$realpart)
#
($bx,$by) = ($realpart*$bx - $by, $bx + $realpart*$by);
}
my $dup = ($seen{"$x,$y"}++ ? " dup" : "");
printf "%4d %2d,%2d%s\n", $n, $x,$y, $dup;
if ($x > -$isize && $x < $isize
&& $y > -$isize && $y < $isize) {
$image->xy($x+$isize,$y+$isize,'*');
}
}
$image->xy(0+$isize,0+$isize,'+');
$image->save_fh(\*STDOUT);
exit 0;
sub digits {
my ($n, $radix) = @_;
my @ret;
while ($n) {
push @ret, $n % $radix;
$n = int($n/$radix);
}
return @ret;
}
}