#!/usr/bin/perl -w
# Copyright 2012, 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::BaseCnv;
use List::MoreUtils;
use POSIX 'floor';
use Math::Libm 'M_PI', 'hypot';
use List::Util 'min', 'max';
use Math::PlanePath::R5DragonCurve;
use Math::BigInt try => 'GMP';
use Math::BigFloat;
use lib 'devel/lib';
use lib 'xt';
use MyOEIS;
# uncomment this to run the ### lines
# use Smart::Comments;
{
  # partial fractions
  require Math::Polynomial;
  Math::Polynomial->string_config({ascending=>1});
  # x^3/((x-1)*(2*x-1)*(x^2-x+1))
  require Math::Complex;
  my $b = Math::Complex->make(1,1);
  my @numerators = MyOEIS::polynomial_partial_fractions
    (Math::Polynomial->new(Math::Complex->make(  3,  -2),
                           Math::Complex->make(-10,   1),
                           Math::Complex->make( 11,   5),
                           Math::Complex->make(  2, -24),
                           Math::Complex->make(-18,  18),
                           Math::Complex->make(  8, -16),
                           Math::Complex->make( 16, -32),
                          ), # numerator
     # ((   16 - 32*I)*x^6
     #  + (  8 - 16*I)*x^5
     #  + (-18 + 18*I)*x^4
     #  + (  2 - 24*I)*x^3
     #  + ( 11 +  5*I)*x^2
     #  + (-10 +    I)*x
     #  + (  3 -  2*I))
     Math::Polynomial->new(1,-$b),            # 1-b*x
     Math::Polynomial->new(1,1)**2,              # 1+x
     Math::Polynomial->new(1,-2,2)**2,           # 1 - 2x + 2x^2
     Math::Polynomial->new(1, -$b, -2*$b**3), # 1 - b*x - 2*b^3*x^3
    );
  print "@numerators\n";
  # my @numerators = MyOEIS::polynomial_partial_fractions
  #   (Math::Polynomial->new(0,0,0,3), # numerator x^3
  #    Math::Polynomial->new(1,-1),     # 1-x
  #    Math::Polynomial->new(1,-2),     # 1-2x
  #    Math::Polynomial->new(1, -1, 1), # 1 - x + x^2
  #   );
  # print "@numerators\n";
  # my @numerators = MyOEIS::polynomial_partial_fractions
  #   (1640 * Math::Polynomial->new(1), # numerator 1
  #    Math::Polynomial->new(1,-1),     # 1-x
  #    Math::Polynomial->new(1, 1),     # 1+x
  #    Math::Polynomial->new(1, -2, 2), # 1 - 2*x + 2*x^2
  #   );
  # print "@numerators\n";
  # use Math::BigRat;
  # my $o = Math::BigRat->new(1);
  # $o = 1;
  # my @numerators = MyOEIS::polynomial_partial_fractions
  #   (1640 * Math::Polynomial->new(0*$o/10, 65*$o/10, 18*$o/10, 13*$o/10), # numerator 13*x^2 + 18*x + 65
  #    Math::Polynomial->new(25*$o, 6*$o, 1*$o),  # (25 + 6*x + x^2)
  #    Math::Polynomial->new(1*$o, -5*$o),        #  * (1-5*x)
  #    Math::Polynomial->new(1*$o, -1*$o));       # (1-x)
  # print "@numerators\n";
  # dragon dir N touching next level
  # p = (1-2*x^3)/(1-2*x-x^3+2*x^4)
  # (1-2*x^3)/((1-2*x)*(1-x)*(1+x+x^2)) * 21 == (18+12*x)/(1+x+x^2) + 3/(1-2*x)
  # p*21 == ((3 + 6*x + 12*x^2)/(1-x^3) + 3/(1-2*x)
  # p*21 == (-4 -5*x)/(1+x+x^2) + 7/(1-x) + 18/(1-2*x)
  # my @numerators = MyOEIS::polynomial_partial_fractions
  #   (21 * Math::Polynomial->new(1,0,0,-2), # numerator 1-2*x^3
  #     Math::Polynomial->new(1,0,0,-1),     # 1-x^3
  #    # Math::Polynomial->new(1,1,1),    # 1+x+x^2
  #    # Math::Polynomial->new(1,-1),     # 1-x
  #    Math::Polynomial->new(1,-2));    # 1-2x
  # print "@numerators\n";
  # # dragon JA[k] area
  # # x^4/ ((1 - x - 2*x^3)*(1-x)*(1-2*x))
  # my @numerators = MyOEIS::polynomial_partial_fractions
  #   (Math::Polynomial->new(1),         # numerator
  #    Math::Polynomial->new(1,-1,0,-2), # 1-x-2*x^3
  #    Math::Polynomial->new(1,-1));     # 1-x
  # print "@numerators\n";
  # # dragon A[k] area
  # # x^4/ ((1 - x - 2*x^3)*(1-x)*(1-2*x))
  # my @numerators = MyOEIS::polynomial_partial_fractions
  #   (Math::Polynomial->new(2),         # numerator
  #    Math::Polynomial->new(1,-1,0,-2), # 1-x-2*x^3
  #    Math::Polynomial->new(1,-2),      # 1-2*x
  #    Math::Polynomial->new(1,-1));     # 1-x
  # print "@numerators\n";
  # # dragon B[k]=R[k+1] total boundary
  # # (4 + 2 x + 4 x^2)/(1-x-2*x^3) +  (-2)/(1-x)
  # my @numerators = MyOEIS::polynomial_partial_fractions
  #   (Math::Polynomial->new(2,0,2),     # numerator reduced 2*x + 2*x^3
  #    Math::Polynomial->new(1,-1,0,-2), # 1-x-2*x^3
  #    Math::Polynomial->new(1,-1));     # 1-x
  # print "@numerators\n";
  # # dragon R right boundary
  # my @numerators = MyOEIS::polynomial_partial_fractions
  #   (Math::Polynomial->new(1,0,1,0,2),
  #    Math::Polynomial->new(1,-1,0,-2),
  #    Math::Polynomial->new(1,-1));
  # print "@numerators\n";
  exit 0;
}
{
  # convex hull
  # hull 8 new vertices
  require Math::Geometry::Planar;
  my $points = [ [0,0], [1,0], [0,0] ];
  $points = [ [Math::BigInt->new(0), Math::BigInt->new(0)],
              [Math::BigInt->new(1), Math::BigInt->new(0)],
              [Math::BigInt->new(0), Math::BigInt->new(0)] ];
  my $end_x = Math::BigInt->new(1);
  my $end_y = Math::BigInt->new(0);
  my $path = Math::PlanePath::R5DragonCurve->new;
  my $num_points_prev = 0;
  for (my $k = Math::BigInt->new(0);
       $k < 40;
       $k++) {
    my $angle = 0; # Math::BigFloat->new($end_y)->batan2(Math::BigFloat->new($end_x), 10);
    my $num_points = scalar(@$points);
    my $num_points_diff = $num_points - $num_points_prev;
    print "k=$k end=$end_x,$end_y a=$angle  $num_points diff=$num_points_diff\n";
    my @new_points = @$points;
    {
      my $p = Math::Geometry::Planar->new;
      $p->points(points_copy($points));
      $p->move (-$end_y, $end_x);
      push @new_points, @{$p->points};
      ### move 1: $p->points
    }
    {
      my $p = Math::Geometry::Planar->new;
      $p->points(points_copy($points));
      $p->move (2*-$end_y, 2*$end_x);
      push @new_points, @{$p->points};
      ### move 2: $p->points
    }
    {
      my $p = Math::Geometry::Planar->new;
      $p->points(points_copy($points));
      planar_rotate_plus90($p);
      push @new_points, @{$p->points};
      ### rot: $p->points
    }
    {
      my $p = Math::Geometry::Planar->new;
      $p->points(points_copy($points));
      planar_rotate_plus90($p);
      $p->move ($end_x + -$end_y, $end_y + $end_x);
      push @new_points, @{$p->points};
      ### rot move: $p->points
    }
    my $p = Math::Geometry::Planar->new;
    $p->points(\@new_points);
    $p = $p->convexhull2;
    $points = $p->points;
    ($end_x,$end_y) = ($end_x - 2*$end_y,
                       $end_y + 2*$end_x);
    $num_points_prev = $num_points;
    my ($x,$y) = $path->n_to_xy(5**($k+1));
    ### $end_y
    ### $y
    $x == $end_x or die;
    $y == $end_y or die;
  }
  exit 0;
  sub planar_rotate_plus90 {
    my ($planar) = @_;
    my $points = $planar->points;
    foreach my $p (@$points) {
      ($p->[0],$p->[1]) = (- $p->[1], $p->[0]);
    }
    return $planar;
  }
  sub points_copy {
    my ($points) = @_;
    return [ map {[$_->[0],$_->[1]]} @$points ];
  }
  # {
  #   my $pl = Math::Geometry::Planar->new;
  #   $pl->points($points);
  #   $pl->rotate(- atan2(2,1));
  #   $pl->scale(1/sqrt(5));
  #   $points = $pl->points;
  # }
}
{
  # extents  h->4/5 w->2/5
  #            1/sqrt(5)
  #           *--*              1/5 + 4/5 = 1
  # 2/sqrt(5) | /  1
  #           |/
  #           *
  #
  my $h = 0;
  my $w = 0;
  my $sum = 0;
  foreach my $k (0 .. 20) {
    print "$h $w $sum\n";
    $sum += (3/5)**$k;
    $h /= sqrt(5);
    $w /= sqrt(5);
    my $s = 1/sqrt(5);
    my $add = $s * 2/sqrt(5);
    ($h, $w) = ($h*2/sqrt(5) + $w*1/sqrt(5) + $add,
                $h*2/sqrt(5) + $w*1/sqrt(5));
  }
  exit 0;
}
{
  # min/max for level
  # radial extent
  #
  # dist0to5 = sqrt(1*1+2*2) = sqrt(5)
  #
  #   4-->5
  #   ^
  #   |
  #   3<--2
  #       ^
  #       |
  #   0-->1
  #
  # Rlevel = sqrt(5)^level + Rprev
  #        = sqrt(5) + sqrt(5)^2 + ... + sqrt(5)^(level-1) + sqrt(5)^level
  # if level 
  #        = sqrt(5) + sqrt(5)^2 + sqrt(5)*sqrt(5)^2 + ... 
  #        = sqrt(5) + (1+sqrt(5))*5^1 + (1+sqrt(5))*5^2 + ... 
  #        = sqrt(5) + (1+sqrt(5))* [ 5^1 + 5^2 + ... ]
  #        = sqrt(5) + (1+sqrt(5))* (5^k - 1)/4
  #        <= 5^k
  # Rlevel^2 <= 5^level
  require Math::BaseCnv;
  require Math::PlanePath::R5DragonCurve;
  my $path = Math::PlanePath::R5DragonCurve->new;
  my $prev_min = 1;
  my $prev_max = 1;
  for (my $level = 1; $level < 10; $level++) {
    my $n_start = 5**($level-1);
    my $n_end = 5**$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 $min_hypot   at $min_x,$min_y\n";
    # print "  max $max_hypot   at $max_x,$max_y\n";
    {
      my $factor = $min_hypot / $prev_min;
      my $min_hypot_5 = Math::BaseCnv::cnv($min_hypot,10,5);
      print "  min r^2 $min_hypot ${min_hypot_5}[5]  at $min_pos  factor $factor\n";
    }
    {
      my $factor = $max_hypot / $prev_max;
      my $max_hypot_5 = Math::BaseCnv::cnv($max_hypot,10,5);
      print "  max r^2 $max_hypot ${max_hypot_5}[5])  at $max_pos  factor $factor\n";
    }
    $prev_min = $min_hypot;
    $prev_max = $max_hypot;
  }
  exit 0;
}
{
  # boundary length between arms = 2*3^k
  #
  #         *---1    length=6
  #         |
  # 2   *---*---*
  # |   |   |   |
  # *---*   0---*
  #
  # T[0] = 2
  # T[k+1] = R[k] + T[k] + U[k]
  # T[k+1] = 4*3^k + T[k]
  #              i=k-1
  # T[k] = 2 +  sum    4*3^i
  #              i=0
  #      = 2 + 4*(3^k - 1)/(3-1)
  #      = 2 + 2*(3^k - 1)
  #      = 2*3^k
  my $arms = 2;
  my $path = Math::PlanePath::R5DragonCurve->new (arms => $arms);
  my @values;
  foreach my $k (0 .. 8) {
    my $n_limit = $arms * 5**$k + $arms-1;
    my $n_from = $n_limit-1;
    my $n_to = $n_limit;
    print "k=$k  n_limit=$n_limit\n";
    my $points = MyOEIS::path_boundary_points_ft ($path, $n_limit,
                                                  $path->n_to_xy($n_from),
                                                  $path->n_to_xy($n_to),
                                                  side => 'right',
                                                 );
    if (@$points < 10) {
      foreach my $p (@$points) {
        print " $p->[0],$p->[1]";
      }
      print "\n";
    }
    my $length = scalar(@$points) - 1;
    print "  length $length\n";
    push @values, $length;
  }
  shift @values;
  print join(',',@values),"\n";
  Math::OEIS::Grep->search(array => \@values);
  exit 0;
}
{
  # right boundary N
  my $path = Math::PlanePath::R5DragonCurve->new;
  my %non_values;
  my %n_values;
  my @n_values;
  my @values;
  foreach my $k (3){
    my $n_limit = 5**$k;
    print "k=$k  n_limit=$n_limit\n";
    foreach my $n (0 .. $n_limit-1) {
      $non_values{$n} = 1;
    }
    my $points = MyOEIS::path_boundary_points ($path, $n_limit,
                                               side => 'right',
                                              );
    ### $points
    for (my $i = 0; $i+1 <= $#$points; $i++) {
      my ($x,$y) = @{$points->[$i]};
      my ($x2,$y2) = @{$points->[$i+1]};
      # my @n_list = $path->xy_to_n_list($x,$y);
      my @n_list = path_xyxy_to_n($path, $x,$y, $x2,$y2);
      foreach my $n (@n_list) {
        delete $non_values{$n};
        if ($n <= $n_limit) { $n_values{$n} = 1; }
        my $n5 = Math::BaseCnv::cnv($n,10,5);
        my $pred = $path->_UNDOCUMENTED__n_segment_is_right_boundary($n);
        my $diff = $pred ? '' : '  ***';
        if ($k <= 4) { print "$n  $n5$diff\n"; }
      }
    }
    @n_values = keys %n_values;
    @n_values = sort {$a<=>$b} @n_values;
    my @non_values = keys %non_values;
    @non_values = sort {$a<=>$b} @non_values;
    my $count = scalar(@n_values);
    print "count $count\n";
    # push @values, $count;
    @values = @n_values;
    if ($k <= 4) {
      foreach my $n (@non_values) {
        my $pred = $path->_UNDOCUMENTED__n_segment_is_right_boundary($n);
        my $diff = $pred ? '  ***' : '';
        my $n5 = Math::BaseCnv::cnv($n,10,5);
        print "non $n  $n5$diff\n";
      }
    }
    # @values = @non_values;
    # print "func ";
    # foreach my $i (0 .. $count-1) {
    #   my $n = $path->_UNDOCUMENTED__right_boundary_i_to_n($i);
    #   my $n5 = Math::BaseCnv::cnv($n,10,5);
    #   print "$n,";
    # }
    # print "\n";
    print "vals ";
    foreach my $i (0 .. $count-1) {
      my $n = $values[$i];
      my $n5 = Math::BaseCnv::cnv($n,10,5);
      print "$n,";
    }
    print "\n";
  }
  @values = MyOEIS::first_differences(@values);
  shift @values;
  shift @values;
  shift @values;
  print join(',',@values),"\n";
  Math::OEIS::Grep->search(array => \@values);
  exit 0;
  sub path_xyxy_to_n {
    my ($path, $x1,$y1, $x2,$y2) = @_;
    ### path_xyxy_to_n(): "$x1,$y1, $x2,$y2"
    my @n_list = $path->xy_to_n_list($x1,$y1);
    ### @n_list
    my $arms = $path->arms_count;
    foreach my $n (@n_list) {
      my ($x,$y) = $path->n_to_xy($n + $arms);
      if ($x == $x2 && $y == $y2) {
        return $n;
      }
    }
    return;
  }
}
{
  my $C = sub {
    my ($k) = @_;
    return 3**$k - $k;    # A024024
  };
  my $E = sub {
    my ($k) = @_;
    return 3**$k + $k;    # A104743
  };
  my @values = map { $E->($_) } 0 .. 10;
  print join(',',@values),"\n";
  Math::OEIS::Grep->search(array => \@values);
  exit;
}
{
  # left boundary N
  my $path = Math::PlanePath::R5DragonCurve->new;
  my %non_values;
  my %n_values;
  my @n_values;
  my @values;
  foreach my $k (2) {
    my $n_limit = 3*5**$k;
    print "k=$k  n_limit=$n_limit\n";
    foreach my $n (0 .. $n_limit-1) {
      $non_values{$n} = 1;
    }
    my $points = MyOEIS::path_boundary_points ($path, $n_limit,
                                               side => 'left',
                                              );
    @$points = reverse @$points; # for left
    ### $points
    for (my $i = 0; $i+1 <= $#$points; $i++) {
      my ($x,$y) = @{$points->[$i]};
      my ($x2,$y2) = @{$points->[$i+1]};
      # my @n_list = $path->xy_to_n_list($x,$y);
      my @n_list = path_xyxy_to_n($path, $x,$y, $x2,$y2);
      foreach my $n (@n_list) {
        delete $non_values{$n};
        if ($n <= $n_limit) { $n_values{$n} = 1; }
        my $n5 = Math::BaseCnv::cnv($n,10,5);
        my $pred = $path->_UNDOCUMENTED__n_segment_is_left_boundary($n);
        my $diff = $pred ? '' : '  ***';
        if ($k <= 4) { print "$n  $n5$diff\n"; }
      }
    }
    @n_values = keys %n_values;
    @n_values = sort {$a<=>$b} @n_values;
    my @non_values = keys %non_values;
    @non_values = sort {$a<=>$b} @non_values;
    my $count = scalar(@n_values);
    print "count $count\n";
    # push @values, $count;
    @values = @n_values;
    if ($k <= 4) {
      foreach my $n (@non_values) {
        my $pred = $path->_UNDOCUMENTED__n_segment_is_left_boundary($n);
        my $diff = $pred ? '  ***' : '';
        my $n5 = Math::BaseCnv::cnv($n,10,5);
        print "non $n  $n5$diff\n";
      }
    }
    # @values = @non_values;
    # print "func ";
    # foreach my $i (0 .. $count-1) {
    #   my $n = $path->_UNDOCUMENTED__left_boundary_i_to_n($i);
    #   my $n5 = Math::BaseCnv::cnv($n,10,5);
    #   print "$n,";
    # }
    # print "\n";
    print "vals ";
    foreach my $i (0 .. $count-1) {
      my $n = $values[$i];
      my $n5 = Math::BaseCnv::cnv($n,10,5);
      print "$n5,";
    }
    print "\n";
  }
  # @values = MyOEIS::first_differences(@values);
  shift @values;
  shift @values;
  shift @values;
  print join(',',@values),"\n";
  Math::OEIS::Grep->search(array => \@values);
  exit 0;
}
{
  # recurrence
  # v3 = a*v0 + b*v1 + c*v2
  #  [v0 v1 v2] [a]   [v3]
  #  [v1 v2 v3] [b] = [v4]
  #  [v2 v3 v4] [c]   [v5]
  #  [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,";
  # }
  print "\n";
  array_to_recurrence_pari(\@array);
  print "\n";
  my @recurr = array_to_recurrence(\@array);
  print join(', ',@recurr),"\n";
  exit 0;
  sub array_to_recurrence_pari {
    my ($aref) = @_;
    my $order = int(scalar(@array)/2); # 2*order-1 = @array-1
    my $str = "m=[";
    foreach my $i (0 .. $order-1) {
      if ($i) { $str .= "; " }
      foreach my $j (0 .. $order-1) {
        if ($j) { $str .= "," }
        $str .= $aref->[$i+$j];
      }
    }
    $str .= "]\n";
    $str .= "v=[";
    foreach my $i ($order .. 2*$order-1) {
      if ($i > $order) { $str .= ";" }
      $str .= $aref->[$i];
    }
    $str .= "];";
    $str .= "(m^-1)*v\n";
    print $str;
    require IPC::Run;
    IPC::Run::run(['gp'],'<',\$str);
  }
  sub array_to_recurrence {
    my ($aref) = @_;
    # 2*order-1 = @array-1
    my $order = int(scalar(@array)/2);
    require Math::Matrix;
    my $m = Math::Matrix->new(map {[
                                    map { $array[$_]
                                        } $_ .. $_+$order-1
                                   ]}
                              0 .. $order-1);
    print $m;
    print $m->determinant,"\n";
    my $v = Math::Matrix->new(map {[ $array[$_] ]} $order .. 2*$order-1);
    print $v;
    $m = $m->invert;
    print $m;
    $v = $m*$v;
    print $v;
    return (map {$v->[$_][0]} reverse 0 .. $order-1);
  }
}
{
  # at N=29
  require Math::NumSeq::PlanePathDelta;
  require Math::PlanePath::R5DragonMidpoint;
  my $path = Math::PlanePath::R5DragonMidpoint->new;
  my $n = 29;
  my ($x,$y) = $path->n_to_xy($n);
  my ($dx,$dy) = $path->n_to_dxdy($n);
  my $tradius = Math::NumSeq::PlanePathCoord::_path_n_to_tradius($path,$n);
  my $next_tradius = Math::NumSeq::PlanePathCoord::_path_n_to_tradius($path,$n + $path->arms_count);
  my $dtradius = Math::NumSeq::PlanePathDelta::_path_n_to_dtradius($path,$n);
  print "$n  x=$x,y=$y   $dx,$dy  dtradius=$dtradius\n";
  print "   tradius $tradius to $next_tradius\n";
  exit 0;
}
{
  # first South step dY=-1 on Y axis
  require Math::PlanePath::R5DragonMidpoint;
  my $path = Math::PlanePath::R5DragonMidpoint->new;
  require Math::NumSeq::PlanePathDelta;
  my $seq = Math::NumSeq::PlanePathDelta->new (path => $path);
  my @values;
  my $n = 0;
 OUTER: for ( ; ; $n++) {
    my ($x,$y) = $path->n_to_xy($n);
    my ($dx,$dy) = $path->n_to_dxdy($n);
    if ($x == 0 && $dx == 0 && $dy == -($y < 0 ? -1 : 1)) {
      my $tradius = Math::NumSeq::PlanePathCoord::_path_n_to_tradius($path,$n);
      my $next_tradius = Math::NumSeq::PlanePathCoord::_path_n_to_tradius($path,$n + $path->arms_count);
      my $dtradius = Math::NumSeq::PlanePathDelta::_path_n_to_dtradius($path,$n);
      print "$n  $x,$y   $dx,$dy  dtradius=$dtradius\n";
      print "   tradius $tradius to $next_tradius\n";
      push @values, $n;
      last OUTER if @values > 20;
    }
  }
  print join(',',@values),"\n";
  Math::OEIS::Grep->search(array => \@values);
  exit 0;
}
{
  # any South step dY=-1 on Y axis
  # use Math::BigInt try => 'GMP';
  # use Math::BigFloat;
  require Math::PlanePath::R5DragonMidpoint;
  my $path = Math::PlanePath::R5DragonMidpoint->new;
  require Math::NumSeq::PlanePathDelta;
  my $seq = Math::NumSeq::PlanePathDelta->new (path => $path);
  my @values;
  my $x = 0;
  my $y = 0;
  # $x = Math::BigFloat->new($x);
  # $y = Math::BigFloat->new($y);
  OUTER: for ( ; ; $y++) {
    ### y: "$y"
    foreach my $sign (1,-1) {
      ### at: "$x, $y  sign=$sign"
      if (defined (my $n = $path->xy_to_n($x,$y))) {
        my ($dx,$dy) = $path->n_to_dxdy($n);
        ### dxdy: "$dx, $dy"
        if ($dx == 0 && $dy == $sign) {
          my $tradius = Math::NumSeq::PlanePathCoord::_path_n_to_tradius($path,$n);
          my $next_tradius = Math::NumSeq::PlanePathCoord::_path_n_to_tradius($path,$n + $path->arms_count);
          my $dtradius = Math::NumSeq::PlanePathDelta::_path_n_to_dtradius($path,$n);
          print "$n  $x,$y   $dx,$dy  dtradius=$dtradius\n";
          print "   tradius $tradius to $next_tradius\n";
          push @values, $y;
          last OUTER if @values > 20;
        }
      }
      $y = -$y;
    }
  }
  print join(',',@values),"\n";
  require Math::OEIS::Grep;
  Math::OEIS::Grep->search(array => \@values);
  exit 0;
}
{
  # boundary join 4,13,40,121,364
  # A003462 (3^n - 1)/2.
  require Math::PlanePath::R5DragonCurve;
  my $path = Math::PlanePath::R5DragonCurve->new;
  my @values;
  $| = 1;
  foreach my $exp (2 .. 6) {
    my $t_lo = 5**$exp;
    my $t_hi = 2*5**$exp - 1;
    my $count = 0;
    foreach my $n (0 .. $t_lo-1) {
      my ($x,$y) = $path->n_to_xy($n);
      my @n_list = $path->xy_to_n_list($x,$y);
      if (@n_list >= 2
          && $n_list[0] < $t_lo
          && $n_list[1] >= $t_lo
          && $n_list[1] < $t_hi) {
        $count++;
      }
    }
    push @values, $count;
    print "$count,";
  }
  print "\n";
  require Math::OEIS::Grep;
  Math::OEIS::Grep->search(array => \@values);
  exit 0;
}
{
  # overlaps
  require Math::PlanePath::R5DragonCurve;
  require Math::BaseCnv;
  my $path = Math::PlanePath::R5DragonCurve->new;
  my $width = 5;
  foreach my $n (0 .. 5**($width-1)) {
    my ($x,$y) = $path->n_to_xy($n);
    my @n_list = $path->xy_to_n_list($x,$y);
    next unless @n_list >= 2;
    if ($n_list[1] == $n) { ($n_list[0],$n_list[1]) = ($n_list[1],$n_list[0]); }
    my $n_list = join(',',@n_list);
    my @n5_list = map { sprintf '%*s', $width, Math::BaseCnv::cnv($_,10,5) } @n_list;
    print "$n5_list[0]  $n5_list[1]  ($n_list)\n";
  }
  exit 0;
}
{
  # tiling
  require Image::Base::Text;
  require Math::PlanePath::R5DragonCurve;
  my $path = Math::PlanePath::R5DragonCurve->new;
  my $width = 37;
  my $height = 21;
  my $image = Image::Base::Text->new (-width => $width,
                                      -height => $height);
  my $xscale = 3;
  my $yscale = 2;
  my $w2 = int(($width+1)/2);
  my $h2 = int($height/2);
  $w2 -= $w2 % $xscale;
  $h2 -= $h2 % $yscale;
  my $affine = sub {
    my ($x,$y) = @_;
    return ($x*$xscale + $w2,
            -$y*$yscale + $h2);
  };
  my ($n_lo, $n_hi) = $path->rect_to_n_range(-$w2/$xscale, -$h2/$yscale,
                                             $w2/$xscale, $h2/$yscale);
  print "n to $n_hi\n";
  foreach my $n ($n_lo .. $n_hi) {
    next if ($n % 5) == 2;
    my ($x,$y) = $path->n_to_xy($n);
    my ($next_x,$next_y) = $path->n_to_xy($n+1);
    foreach (1 .. 4) {
      $image->line ($affine->($x,$y),
                    $affine->($next_x,$next_y),
                    ($x==$next_x ? '|' : '-'));
      $image->xy ($affine->($x,$y),
                  '+');
      $image->xy ($affine->($next_x,$next_y),
                  '+');
      ($x,$y) = (-$y,$x); # rotate +90
      ($next_x,$next_y) = (-$next_y,$next_x); # rotate +90
    }
  }
  $image->xy ($affine->(0,0),
              'o');
  foreach my $x (0 .. $width-1) {
    foreach my $y (0 .. $height-1) {
      next unless $image->xy($x,$y) eq '+';
      if ($x > 0 && $image->xy($x-1,$y) eq ' ') {
        $image->xy($x,$y, '|');
      } elsif ($x < $width-1 && $image->xy($x+1,$y) eq ' ') {
        $image->xy($x,$y, '|');
      } elsif ($y > 0 && $image->xy($x,$y-1) eq ' ') {
        $image->xy($x,$y, '-');
      } elsif ($y < $height-1 && $image->xy($x,$y+1) eq ' ') {
        $image->xy($x,$y, '-');
      }
    }
  }
  $image->save('/dev/stdout');
  exit 0;
}
{
  # area recurrence
  foreach my $i (0 .. 10) {
    print recurrence($i),",";
  }
  print "\n";
  print "wrong():              ";
  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;
  # A[n+1] = 4*A[n] - 3*A[n-1] + 4*5^(n-1)
  # - A[n+1] + 4*A[n] + 4*5^(n-1) = 3*A[n-1]
  # 3*A[n-1] = - A[n+1] + 4*A[n] + 4*5^(n-1)
  # 3*A[n-2] = - A[n] + 4*A[n-1] + 4*5^(n-2)
  # D[n+1] = 4*A[n] - 3*A[n-1]            + 4*5^(n-1)
  #          -       (4*A[n-1] - 3*A[n-2] + 4*5^(n-2))
  #        = 4*A[n] - 3*A[n-1]            + 4*5^(n-1)
  #                 - 4*A[n-1] + 3*A[n-2] - 4*5^(n-2))
  #        = 4*A[n] - 3*A[n-1]            + 4*5^(n-1)
  #                 - 4*A[n-1] - A[n] + 4*A[n-1] + 4*5^(n-2) - 4*5^(n-2))
  #        = 4*A[n] - 3*A[n-1]            + 4*5^(n-1)
  #          - A[n]
  # D[n+1] = 4*A[n] - 3*A[n-1]            + 4*5^(n-1)
  #          - A[n]
  # D[n+1] = 3*A[n] - 3*A[n-1]            + 4*5^(n-1)
  # D[n+1] = 3*D[n] + 4*5^(n-1)
  #      = 4*A[n] - 7*A[n-1] + 3*A[n-2] + (4*5-4)*5^(n-2)
  #      = 4*A[n] - 7*A[n-1] + 3*A[n-2] + 16*5^(n-2)
  #      = 4*A[n] - 7*A[n-1] + A[n] + 4*A[n-1] + 4*5^(n-2) + 16*5^(n-2)
  #      = 3*A[n] - 3*A[n-1] + 20*5^(n-2)
  # 4*A[n] - 12*A[n-1] + 4 - 4*5^(n-1) = 0 ??
  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);
  }
  # A[n+1] = 8*A[n] - 15*A[n-1] + 4
  sub recurrence_area815 {
    my ($n) = @_;
    if ($n <= 0) { return 0; }
    if ($n == 1) { return 0; }
    return 8*recurrence_area815($n-1) - 15*recurrence_area815($n-2) + 4;
  }
  sub recurrence {
    my ($n) = @_;
    if ($n <= 0) { return 0; }
    if ($n == 1) { return 2; }
    return 8*recurrence($n-1) - 15*recurrence($n-2) + 2;
  }
  sub recurrence_area2S {
    my ($n) = @_;
    return 2*recurrence_S($n+1);
  }
  sub recurrence_areaSU {
    my ($n) = @_;
    return 4*recurrence_S($n) + 2*recurrence_U($n);
  }
  sub recurrence_S {
    my ($n) = @_;
    if ($n <= 0) { return 0; }
    if ($n == 1) { return 0; }
    return 2*recurrence_S($n-1) + recurrence_U($n-1);
  }
  sub recurrence_U {
    my ($n) = @_;
    if ($n <= 0) { return 0; }
    if ($n == 1) { return 0; }
    return recurrence_S($n-1) + 2*recurrence_U($n-1) + 2*5**($n-2);
  }
  # A(n)=a(n)*2
  # A(n)/2 = 8*A(n-1)/2 - 15*A(n-2)/2 + 2
  # A(n) = 8*A(n-1) - 15*A(n-2) + 4
}
{
  # arm xy modulus
  require Math::PlanePath::R5DragonMidpoint;
  my $path = Math::PlanePath::R5DragonMidpoint->new (arms => 4);
  my %dxdy_to_digit;
  my %seen;
  for (my $n = 0; $n < 6125; $n++) {
    my $digit = $n % 5;
    foreach my $arm (0 .. 3) {
      my ($x,$y) = $path->n_to_xy(4*$n+$arm);
      my $nb = int($n/5);
      my ($xb,$yb) = $path->n_to_xy(4*$nb+$arm);
      # (x+iy)*(1+2i) = x-2y + 2x+y
      ($xb,$yb) = ($xb-2*$yb, 2*$xb+$yb);
      my $dx = $xb - $x;
      my $dy = $yb - $y;
      my $dxdy = "$dx,$dy";
      my $show = "${dxdy}[$digit]";
      $seen{$x}{$y} = $show;
      if ($dxdy eq '0,0') {
      }
      # if (defined $dxdy_to_digit{$dxdy} && $dxdy_to_digit{$dxdy} != $digit) {
      #   die;
      # }
      $dxdy_to_digit{$dxdy} = $digit;
    }
  }
  foreach my $y (reverse -45 .. 45) {
    foreach my $x (-5 .. 5) {
      printf " %9s", $seen{$x}{$y}//'e'
    }
    print "\n";
  }
  ### %dxdy_to_digit
  exit 0;
}
{
  # Midpoint xy to n
  require Math::PlanePath::DragonMidpoint;
  require Math::BaseCnv;
  my @yx_adj_x = ([0,1,1,0],
                  [1,0,0,1],
                  [1,0,0,1],
                  [0,1,1,0]);
  my @yx_adj_y = ([0,0,1,1],
                  [0,0,1,1],
                  [1,1,0,0],
                  [1,1,0,0]);
  sub xy_to_n {
    my ($self, $x,$y) = @_;
    my $n = ($x * 0 * $y) + 0; # inherit bignum 0
    my $npow = $n + 1;         # inherit bignum 1
    while (($x != 0 && $x != -1) || ($y != 0 && $y != 1)) {
      # my $ax = ((($x+1) ^ ($y+1)) >> 1) & 1;
      # my $ay = (($x^$y) >> 1) & 1;
      # ### assert: $ax == - $yx_adj_x[$y%4]->[$x%4]
      # ### assert: $ay == - $yx_adj_y[$y%4]->[$x%4]
      my $y4 = $y % 4;
      my $x4 = $x % 4;
      my $ax = $yx_adj_x[$y4]->[$x4];
      my $ay = $yx_adj_y[$y4]->[$x4];
      ### at: "$x,$y  n=$n  axy=$ax,$ay  bit=".($ax^$ay)
      if ($ax^$ay) {
        $n += $npow;
      }
      $npow *= 2;
      $x -= $ax;
      $y -= $ay;
      ### assert: ($x+$y)%2 == 0
      ($x,$y) = (($x+$y)/2,   # rotate -45 and divide sqrt(2)
                 ($y-$x)/2);
    }
    ### final: "xy=$x,$y"
    my $arm;
    if ($x == 0) {
      if ($y) {
        $arm = 1;
        ### flip ...
        $n = $npow-1-$n;
      } else { #  $y == 1
        $arm = 0;
      }
    } else { # $x == -1
      if ($y) {
        $arm = 2;
      } else {
        $arm = 3;
        ### flip ...
        $n = $npow-1-$n;
      }
    }
    ### $arm
    my $arms_count = $self->arms_count;
    if ($arm > $arms_count) {
      return undef;
    }
    return $n * $arms_count + $arm;
  }
  foreach my $arms (4,3,1,2) {
    ### $arms
    my $path = Math::PlanePath::DragonMidpoint->new (arms => $arms);
    for (my $n = 0; $n < 50; $n++) {
      my ($x,$y) = $path->n_to_xy($n)
        or next;
      my $rn = xy_to_n($path,$x,$y);
      my $good = '';
      if (defined $rn && $rn == $n) {
        $good .= "good N";
      }
      my $n2 = Math::BaseCnv::cnv($n,10,2);
      my $rn2 = Math::BaseCnv::cnv($rn,10,2);
      printf "n=%d xy=%d,%d got rn=%d    %s\n",
        $n,$x,$y,
          $rn,
            $good;
    }
  }
  exit 0;
}
{
  # 2i+1 powers
  my $x = 1;
  my $y = 0;
  foreach (1 .. 10) {
    ($x,$y) = ($x - 2*$y,
                 $y + 2*$x);
    print "$x  $y\n";
  }
  exit 0;
}
{
  # turn sequence
  require Math::NumSeq::PlanePathTurn;
  my @want = (0);
  foreach (1 .. 5) {
    @want = map { $_ ? (0,0,1,1,1) : (0,0,1,1,0) } @want;
  }
  my @got;
  foreach my $i (1 .. @want) {
    push @got, calc_n_turn($i);
  }
  # my $seq = Math::NumSeq::PlanePathTurn->new (planepath => 'R5DragonCurve',
  #                                             turn_type => 'Right');
  # while (@got < @want) {
  #   my ($i,$value) = $seq->next;
  #   push @got, $value;
  # }
  my $got = join(',',@got);
  my $want = join(',',@want);
  print "$got\n";
  print "$want\n";
  if ($got ne $want) {
    die;
  }
  exit 0;
  # return 0 for left, 1 for right
  sub calc_n_turn {
    my ($n) = @_;
    $n or die;
    for (;;) {
      if (my $digit = $n % 5) {
        return ($digit >= 3 ? 1 : 0);
      }
      $n = int($n/5);
    }
  }
}