File Coverage

blib/lib/Math/PlanePath/PythagoreanTree.pm

Criterion	Covered	Total	%
statement	239	391	61.1
branch	93	190	48.9
condition	49	82	59.7
subroutine	33	61	54.1
pod	25	25	100.0
total	439	749	58.6

# Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018 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 .


# math-image --path=PythagoreanTree --all --scale=3

# http://sunilchebolu.wordpress.com/pythagorean-triples-and-the-integer-points-on-a-hyperboloid/

# http://www.math.uconn.edu/~kconrad/blurbs/ugradnumthy/pythagtriple.pdf
#
# http://www.math.ou.edu/~dmccullough/teaching/pythagoras1.pdf
# http://www.math.ou.edu/~dmccullough/teaching/pythagoras2.pdf
#
# http://www.microscitech.com/pythag_eigenvectors_invariants.pdf
#


package Math::PlanePath::PythagoreanTree;
use 5.004;
use strict;
use Carp 'croak';

use vars '$VERSION', '@ISA';
$VERSION = 127;
use Math::PlanePath;
*_divrem = \&Math::PlanePath::_divrem;
*_sqrtint = \&Math::PlanePath::_sqrtint;
@ISA = ('Math::PlanePath');

#use List::Util 'min','max';
*min = \&Math::PlanePath::_min;
*max = \&Math::PlanePath::_max;

use Math::PlanePath::Base::Generic
  'is_infinite',
  'round_nearest';
use Math::PlanePath::Base::Digits
  'round_down_pow',
  'digit_split_lowtohigh',
  'digit_join_lowtohigh';
use Math::PlanePath::GrayCode;

# uncomment this to run the ### lines
# use Smart::Comments;

use constant class_x_negative => 0;
use constant class_y_negative => 0;
use constant tree_num_children_list => (3); # complete ternary tree
use constant tree_n_to_subheight => undef; # complete tree, all infinity

use constant parameter_info_array =>
  [ { name            => 'tree_type',
      share_key       => 'tree_type_uadfb',
      display         => 'Tree Type',
      type            => 'enum',
      default         => 'UAD',
      choices         => ['UAD','UArD','FB','UMT'],
    },
    { name            => 'coordinates',
      share_key       => 'coordinates_abcpqsm',
      display         => 'Coordinates',
      type            => 'enum',
      default         => 'AB',
      choices         => ['AB','AC','BC','PQ', 'SM','SC','MC',
                          # 'BA'
                          # 'UV',  # q from x=y diagonal down at 45-deg
                          # 'RS','ST',  # experimental
                         ],
    },
    { name            => 'digit_order',
      display         => 'Digit Order',
      type            => 'enum',
      default         => 'HtoL',
      choices         => ['HtoL','LtoH'],
    },
  ];

{
  my %UAD_coordinates_always_right = (PQ => 1,
                                      AB => 1,
                                      AC => 1);
  sub turn_any_left {
    my ($self) = @_;
    return ! ($self->{'tree_type'} eq 'UAD'
              && $UAD_coordinates_always_right{$self->{'coordinates'}});
  }
}
{
  my %UAD_coordinates_always_left = (BC => 1);
  sub turn_any_right {
    my ($self) = @_;
    return ! ($self->{'tree_type'} eq 'UAD'
              && $UAD_coordinates_always_left{$self->{'coordinates'}});
  }
}
{
  my %UMT_coordinates_any_straight = (BC => 1,  # UMT at N=5
                                      PQ => 1); # UMT at N=5
  sub turn_any_straight {
    my ($self) = @_;
    return ($self->{'tree_type'} eq 'UMT'
            && $UMT_coordinates_any_straight{$self->{'coordinates'}});
  }
}


#------------------------------------------------------------------------------
{
  my %coordinate_minimum = (A => 3,
                            B => 4,
                            C => 5,
                            P => 2,
                            Q => 1,
                            S => 3,
                            M => 4,
                           );
  sub x_minimum {
    my ($self) = @_;
    return $coordinate_minimum{substr($self->{'coordinates'},0,1)};
  }
  sub y_minimum {
    my ($self) = @_;
    return $coordinate_minimum{substr($self->{'coordinates'},1)};
  }
}
{
  my %diffxy_minimum = (PQ => 1, # octant X>=Y+1 so X-Y>=1
                       );
  sub diffxy_minimum {
    my ($self) = @_;
    return $diffxy_minimum{$self->{'coordinates'}};
  }
}
{
  my %diffxy_maximum = (AC => -2, # C>=A+2 so X-Y<=-2
                        BC => -1, # C>=B+1 so X-Y<=-1
                        SM => -1, # S
                        SC => -2, # S
                        MC => -1, # M
                       );
  sub diffxy_maximum {
    my ($self) = @_;
    return $diffxy_maximum{$self->{'coordinates'}};
  }
}
{
  my %absdiffxy_minimum = (PQ => 1,
                           AB => 1, # X=Y never occurs
                           BA => 1, # X=Y never occurs
                           AC => 2, # C>=A+2 so abs(X-Y)>=2
                           BC => 1,
                           SM => 1, # X=Y never occurs
                           SC => 2, # X<=Y-2
                           MC => 1, # X=Y never occurs
                          );
  sub absdiffxy_minimum {
    my ($self) = @_;
    return $absdiffxy_minimum{$self->{'coordinates'}};
  }
}
use constant gcdxy_maximum => 1;  # no common factor

{
  my %absdx_minimum = ('AB,UAD' => 2,
                       'AB,FB'  => 2,
                       'AB,UMT' => 2,

                       'AC,UAD' => 2,
                       'AC,FB'  => 2,
                       'AC,UMT' => 2,

                       'BC,UAD' => 4,  # at N=37
                       'BC,FB'  => 4,  # at N=2 X=12,Y=13
                       'BC,UMT' => 4,  # at N=2 X=12,Y=13

                       'PQ,UAD' => 0,
                       'PQ,FB'  => 0,
                       'PQ,UMT'  => 0,

                       'SM,UAD' => 1,
                       'SM,FB'  => 1,
                       'SM,UMT' => 2,

                       'SC,UAD' => 1,
                       'SC,FB'  => 1,
                       'SC,UMT' => 1,

                       'MC,UAD' => 3,
                       'MC,FB'  => 3,
                       'MC,UMT' => 1,
                      );
  sub absdx_minimum {
    my ($self) = @_;
    return $absdx_minimum{"$self->{'coordinates'},$self->{'tree_type'}"} || 0;
  }
}
{
  my %absdy_minimum = ('AB,UAD' => 4,
                       'AB,FB'  => 4,
                       'AB,UMT' => 4,

                       'AC,UAD' => 4,
                       'AC,FB'  => 4,
                       'BC,UAD' => 4,
                       'BC,FB'  => 4,
                       'PQ,UAD' => 0,
                       'PQ,FB'  => 1,

                       'SM,UAD' => 3,
                       'SM,FB'  => 3,
                       'SM,UMT' => 1,

                       'SC,UAD' => 4,
                       'SC,FB'  => 4,
                       'MC,UAD' => 4,
                       'MC,FB'  => 4,
                      );
  sub absdy_minimum {
    my ($self) = @_;
    return $absdy_minimum{"$self->{'coordinates'},$self->{'tree_type'}"} || 0;
  }
}

{
  my %dir_minimum_dxdy = (# AB apparent minimum dX=16,dY=8
                          'AB,UAD' => [16,8],
                          'AC,UAD' => [1,1],  # it seems
                          # 'BC,UAD' => [1,0], # infimum
                          # 'SM,UAD' => [1,0], # infimum
                          # 'SC,UAD' => [1,0], # N=255 dX=7,dY=0
                          # 'MC,UAD' => [1,0], # infimum

                          # 'SM,FB' => [1,0], # infimum
                          # 'SC,FB' => [1,0], # infimum
                          # 'SM,FB' => [1,0], # infimum

                          'AB,UMT' => [6,12], # it seems

                          # N=ternary 1111111122 dx=118,dy=40
                          # in general dx=3*4k-2 dy=4k
                          'AC,UMT' => [3,1], # infimum
                          #
                          # 'BC,UMT' => [1,0], # N=31 dX=72,dY=0
                          'PQ,UMT' => [1,1], # N=1
                          'SM,UMT' => [1,0],  # infiumum dX=big,dY=3
                          'SC,UMT' => [3,1],  # like AC
                          # 'MC,UMT' => [1,0],  # at N=31
                         );
  sub dir_minimum_dxdy {
    my ($self) = @_;
    return @{$dir_minimum_dxdy{"$self->{'coordinates'},$self->{'tree_type'}"}
               || [1,0] };
  }
}
{
  # AB apparent maximum dX=-6,dY=-12 at N=3
  # AC apparent maximum dX=-6,dY=-12 at N=3 same
  # PQ apparent maximum dX=-1,dY=-1
  my %dir_maximum_dxdy = ('AB,UAD'   => [-6,-12],
                          'AC,UAD'   => [-6,-12],
                          # 'BC,UAD' => [0,0],
                          'PQ,UAD'   => [-1,-1],
                          # 'SM,UAD' => [0,0],   # supremum
                          # 'SC,UAD' => [0,0],   # supremum
                          # 'MC,UAD' => [0,0],   # supremum

                          # 'AB,FB'  => [0,0],
                          # 'AC,FB'  => [0,0],
                          'BC,FB'    => [1,-1],
                          # 'PQ,FB'  => [0,0],
                          # 'SM,FB'  => [0,0],   # supremum
                          # 'SC,FB'  => [0,0],   # supremum
                          # 'MC,FB'  => [0,0],   # supremum

                          # N=ternary 1111111122 dx=118,dy=-40
                          # in general dx=3*4k-2 dy=-4k
                          'AB,UMT' => [3,-1], # supremum
                          #
                          'AC,UMT' => [-10,-20], # at N=9 apparent maximum
                          # 'BC,UMT' => [0,0],  # apparent approach
                          'PQ,UMT' => [1,-1], # N=2
                          # 'SM,UMT' => [0,0],  # supremum dX=big,dY=-1
                          'SC,UMT' => [-3,-5], # apparent approach
                          # 'MC,UMT' => [0,0], # supremum dX=big,dY=-small
                         );
  sub dir_maximum_dxdy {
    my ($self) = @_;
    return @{$dir_maximum_dxdy{"$self->{'coordinates'},$self->{'tree_type'}"}
               || [0,0]};
  }
}

#------------------------------------------------------------------------------

sub _noop {
  return @_;
}
my %xy_to_pq = (AB => \&_ab_to_pq,
                AC => \&_ac_to_pq,
                BC => \&_bc_to_pqa, # ignoring extra $a return
                PQ => \&_noop,
                SM => \&_sm_to_pq,
                SC => \&_sc_to_pq,
                MC => \&_mc_to_pq,
                UV => \&_uv_to_pq,
                RS => \&_rs_to_pq,
                ST => \&_st_to_pq,
               );
my %pq_to_xy = (AB => \&_pq_to_ab,
                AC => \&_pq_to_ac,
                BC => \&_pq_to_bc,
                PQ => \&_noop,
                SM => \&_pq_to_sm,
                SC => \&_pq_to_sc,
                MC => \&_pq_to_mc,
                UV => \&_pq_to_uv,
                RS => \&_pq_to_rs,
                ST => \&_pq_to_st,
               );

my %tree_types = (UAD  => 1,
                  UArD => 1,
                  FB   => 1,
                  UMT  => 1);
my %digit_orders = (HtoL  => 1,
                    LtoH => 1);
sub new {
  my $self = shift->SUPER::new (@_);
  {
    my $digit_order = ($self->{'digit_order'} ||= 'HtoL');
    $digit_orders{$digit_order}
      || croak "Unrecognised digit_order option: ",$digit_order;
  }
  {
    my $tree_type = ($self->{'tree_type'} ||= 'UAD');
    $tree_types{$tree_type}
      || croak "Unrecognised tree_type option: ",$tree_type;
  }
  {
    my $coordinates = ($self->{'coordinates'} ||= 'AB');
    $self->{'xy_to_pq'} = $xy_to_pq{$coordinates}
      || croak "Unrecognised coordinates option: ",$coordinates;
    $self->{'pq_to_xy'} = $pq_to_xy{$coordinates};
  }
  return $self;
}

sub n_to_xy {
  my ($self, $n) = @_;
  ### PythagoreanTree n_to_xy(): $n

  if ($n < 1) { return; }
  if (is_infinite($n)) { return ($n,$n); }

  {
    my $int = int($n);
    if ($n != $int) {
      my $frac = $n - $int;  # inherit possible BigFloat/BigRat
      my ($x1,$y1) = $self->n_to_xy($int);
      my ($x2,$y2) = $self->n_to_xy($int+1);
      my $dx = $x2-$x1;
      my $dy = $y2-$y1;
      return ($frac*$dx + $x1, $frac*$dy + $y1);
    }
  }

  return &{$self->{'pq_to_xy'}}(_n_to_pq($self,$n));
}

# maybe similar n_to_rsquared() as C^2=(P^2+Q^2)^2
sub n_to_radius {
  my ($self, $n) = @_;

  if (($self->{'coordinates'} eq 'AB'
       || $self->{'coordinates'} eq 'BA'
       || $self->{'coordinates'} eq 'SM')
      && $n == int($n)) {
    if ($n < 1) { return undef; }
    if (is_infinite($n)) { return $n; }
    my ($p,$q) = _n_to_pq($self,$n);
    return $p*$p + $q*$q;  # C=P^2+Q^2
  }

  return $self->SUPER::n_to_radius($n);
}

sub _n_to_pq {
  my ($self, $n) = @_;

  my $ndigits = _n_to_digits_lowtohigh($n);
  ### $ndigits

  if ($self->{'tree_type'} eq 'UArD') {
    Math::PlanePath::GrayCode::_digits_to_gray_reflected($ndigits,3);
    ### gray: $ndigits
  }
  if ($self->{'digit_order'} eq 'HtoL') {
    @$ndigits = reverse @$ndigits;
    ### reverse: $ndigits
  }

  my $zero = $n * 0;

  my $p = 2 + $zero;
  my $q = 1 + $zero;

  if ($self->{'tree_type'} eq 'FB') {
    ### FB ...

    foreach my $digit (@$ndigits) {  # high to low, possibly $digit=undef
      ### $p
      ### $q
      ### $digit

      if ($digit) {
        if ($digit == 1) {
          $q = $p-$q;                   # (2p, p-q)  M2
          $p *= 2;
        } else {
          # ($p,$q) = (2*$p, $p+$q);
          $q += $p;                     # (p+q, 2q)  M3
          $p *= 2;
        }
      } else { # $digit == 0
        # ($p,$q) = ($p+$q, 2*$q);
        $p += $q;                       # (p+q, 2q)  M1
        $q *= 2;
      }
    }
  } elsif ($self->{'tree_type'} eq 'UMT') {
    ### UMT ...

    foreach my $digit (@$ndigits) {  # high to low, possibly $digit=undef
      ### $p
      ### $q
      ### $digit

      if ($digit) {
        if ($digit == 1) {
          $q = $p-$q;                 # (2p, p-q)  M2
          $p *= 2;
        } else { # $digit == 2
          $p += 3*$q;                 # T
          $q *= 2;
        }
      } else { # $digit == 0
        # ($p,$q) = ($p+$q, 2*$q);
        ($p,$q) = (2*$p-$q, $p);      # "U" = (2p-q, p)
      }
    }
  } else {
    ### UAD or UArD ...
    ### assert: $self->{'tree_type'} eq 'UAD' || $self->{'tree_type'} eq 'UArD'

    # # Could optimize high zeros as repeated U
    # # high zeros as repeated U: $depth-scalar(@$ndigits)
    # # U^0 = p,    q
    # # U^1 = 2p-q, p          eg. P=2,Q=1 is 2*2-1,2 = 3,2
    # # U^2 = 3p-2q, 2p-q      eg. P=2,Q=1 is 3*2-2*1,2*2-1 = 4,3
    # # U^3 = 4p-3q, 3p-2q
    # # U^k = (k+1)p-kq, kp-(k-1)q   for k>=2
    # #     = p + k*(p-q), k*(p-q)+q
    # # and with initial p=2,q=1
    # # U^k = 2+k, 1+k
    # #
    # $q = $depth - $#ndigits + $zero;  # count high zeros + 1
    # $p = $q + 1 + $zero;

    foreach my $digit (@$ndigits) {  # high to low, possibly $digit=undef
      ### $p
      ### $q
      ### $digit

      if ($digit) {
        if ($digit == 1) {
          ($p,$q) = (2*$p+$q, $p);      # "A" = (2p+q, p)
        } else {
          $p += 2*$q;                   # "D" = (p+2q, q)
        }
      } else { # $digit==0
        ($p,$q) = (2*$p-$q, $p);        # "U" = (2p-q, p)
      }
    }

  }

  ### final pq: "$p, $q"

  return ($p, $q);
}

# _n_to_digits_lowtohigh() returns an arrayref $ndigits which is a list of
# ternary digits 0,1,2 from low to high which are the position of $n within
# its row of the tree.
# The length of the array is the depth.
#
# depth N  N%3      2*N-1   (N-2)/3*2+1
#   0   1   1         1         1/3
#   1   2   2         3          1
#   2   5   2         9          3
#   3   14  2        27          9
#   4   41  2        81         27       28 + (28/2-1) = 41
#
# (N-2)/3*2+1 rounded down to pow=3^k gives depth=k+1 and base=pow+(pow+1)/2
# is the start of the row base=1,2,5,14,41 etc.
#
# An easier calculation is 2*N-1 rounded down to pow=3^d gives depth=d and
# base=2*pow-1, but 2*N-1 and 2*pow-1 might overflow an integer.  Though
# just yet round_down_pow() goes into floats and so doesn't preserve 64-bit
# integer.  So the technique here helps 53-bit float integers, but not right
# up to 64-bits.
#
sub _n_to_digits_lowtohigh {
  my ($n) = @_;
  ### _n_to_digits_lowtohigh(): $n

  my @ndigits;
  if ($n >= 2) {
    my ($pow) = _divrem($n-2, 3);
    ($pow, my $depth) = round_down_pow (2*$pow+1, 3);
    ### $depth
    ### base: $pow + ($pow+1)/2
    ### offset: $n - $pow - ($pow+1)/2
    @ndigits = digit_split_lowtohigh ($n - $pow - ($pow+1)/2, 3);
    push @ndigits, (0) x ($depth - $#ndigits);   # pad to $depth with 0s
  }
  ### @ndigits
  return \@ndigits;


  # {
  #   my ($pow, $depth) = round_down_pow (2*$n-1, 3);
  #
  #   ### h: 2*$n-1
  #   ### $depth
  #   ### $pow
  #   ### base: ($pow + 1)/2
  #   ### rem n: $n - ($pow + 1)/2
  #
  #   my @ndigits = digit_split_lowtohigh ($n - ($pow+1)/2,  3);
  #   $#ndigits = $depth-1;   # pad to $depth with undefs
  #   ### @ndigits
  #
  #   return \@ndigits;
  # }
}

#------------------------------------------------------------------------------
# xy_to_n()

# Nrow(depth+1) - Nrow(depth)
#   = (3*pow+1)/2 - (pow+1)/2
#   = (3*pow + 1 - pow - 1)/2
#   = (2*pow)/2
#   = pow
#
sub xy_to_n {
  my ($self, $x, $y) = @_;
  $x = round_nearest ($x);
  $y = round_nearest ($y);
  ### PythagoreanTree xy_to_n(): "$x, $y"

  my ($p,$q) = &{$self->{'xy_to_pq'}}($x,$y)
    or return undef;    # not a primitive A,B,C

  unless ($p >= 2 && $q >= 1) {          # must be P > Q >= 1
    return undef;
  }
  if (is_infinite($p)) {
    return $p;  # infinity
  }
  if (is_infinite($q)) {
    return $q;  # infinity
  }
  if ($p%2 == $q%2) {  # must be opposite parity, not same parity
    return undef;
  }

  my @ndigits;  # low to high
  if ($self->{'tree_type'} eq 'FB') {
    for (;;) {
      unless ($p > $q && $q >= 1) {
        return undef;
      }
      last if $q <= 1 && $p <= 2;

      if ($q % 2) {
        ### q odd, p even, digit 1 or 2 ...
        $p /= 2;
        if ($q > $p) {
          ### digit 2, M3 ...
          push @ndigits, 2;
          $q -= $p;  # opp parity of p, and < new p
        } else {
          ### digit 1, M2 ...
          push @ndigits, 1;
          $q = $p - $q;  # opp parity of p, and < p
        }
      } else {
        ### q even, p odd, digit 0, M1 ...
        push @ndigits, 0;
        $q /= 2;
        $p -= $q;  # opp parity of q
      }
      ### descend: "$q / $p"
    }

  } elsif ($self->{'tree_type'} eq 'UMT') {
    for (;;) {
      ### at: "p=$p q=$q"
      my $qmod2 = $q % 2;
      unless ($p > $q && $q >= 1) {
        return undef;
      }
      last if $q <= 1 && $p <= 2;

      if ($p < 2*$q) {
        ($p,$q) = ($q, 2*$q-$p);  # U
        push @ndigits, 0;
      } elsif ($qmod2) {
        $p /= 2;   # M2
        $q = $p - $q;
        push @ndigits, 1;
      } else {
        $q /= 2;    # T
        $p -= 3*$q;
        push @ndigits, 2;
      }
    }

  } else {
    ### UAD or UArD ...
    ### assert: $self->{'tree_type'} eq 'UAD' || $self->{'tree_type'} eq 'UArD'
    for (;;) {
      ### $p
      ### $q
      if ($q <= 0 || $p <= 0 || $p <= $q) {
        return undef;
      }
      last if $q <= 1 && $p <= 2;

      if ($p > 2*$q) {
        if ($p > 3*$q) {
          ### digit 2 ...
          push @ndigits, 2;
          $p -= 2*$q;
        } else {
          ### digit 1
          push @ndigits, 1;
          ($p,$q) = ($q, $p - 2*$q);
        }

      } else {
        ### digit 0 ...
        push @ndigits, 0;
        ($p,$q) = ($q, 2*$q-$p);
      }
      ### descend: "$q / $p"
    }
  }
  ### @ndigits

  if ($self->{'digit_order'} eq 'LtoH') {
    @ndigits = reverse @ndigits;
    ### unreverse: @ndigits
  }
  if ($self->{'tree_type'} eq 'UArD') {
    Math::PlanePath::GrayCode::_digits_from_gray_reflected(\@ndigits,3);
    ### ungray: @ndigits
  }

  my $zero = $x*0*$y;
  ### offset: digit_join_lowtohigh(\@ndigits,3,$zero)
  ### depth: scalar(@ndigits)
  ### Nrow: $self->tree_depth_to_n($zero + scalar(@ndigits))

  return ($self->tree_depth_to_n($zero + scalar(@ndigits))
          + digit_join_lowtohigh(\@ndigits,3,$zero)); # offset into row
}

# numprims(H) = how many with hypot < H
# limit H->inf  numprims(H) / H -> 1/2pi
#
# not exact
sub rect_to_n_range {
  my ($self, $x1,$y1, $x2,$y2) = @_;
  ### PythagoreanTree rect_to_n_range(): "$x1,$y1  $x2,$y2"

  $x1 = round_nearest ($x1);
  $y1 = round_nearest ($y1);
  $x2 = round_nearest ($x2);
  $y2 = round_nearest ($y2);

  my $zero = ($x1 * 0 * $y1 * $x2 * $y2);  # inherit bignum

  ($x1,$x2) = ($x2,$x1) if $x1 > $x2;
  ($y1,$y2) = ($y2,$y1) if $y1 > $y2;
  ### x2: "$x2"
  ### y2: "$y2"

  if ($self->{'coordinates'} eq 'BA') {
    ($x2,$y2) = ($y2,$x2);
  }
  if ($self->{'coordinates'} eq 'SM') {
    if ($x2 > $y2) {   # both max
      $y2 = $x2;
    } else {
      $x2 = $y2;
    }
  }

  if ($self->{'coordinates'} eq 'PQ') {
    if ($x2 < 2 || $y2 < 1) {
      return (1,0);
    }
    # P > Q so reduce y2 to at most x2-1
    if ($y2 >= $x2) {
      $y2 = $x2-1;    # $y2 = min ($y2, $x2-1);
    }

    if ($y2 < $y1) {
      ### PQ y range all above X=Y diagonal ...
      return (1,0);
    }
  } else {
    # AB,AC,BC, SM,SC,MC
    if ($x2 < 3 || $y2 < 0) {
      return (1,0);
    }
  }

  my $depth;
  if ($self->{'tree_type'} eq 'FB') {
    ### FB ...
    if ($self->{'coordinates'} eq 'PQ') {
      $x2 *= 3;
    }
    my ($pow, $exp) = round_down_pow ($x2, 2);
    $depth = 2*$exp;
  } else {
    ### UAD or UArD, and UMT ...
    if ($self->{'coordinates'} eq 'PQ') {
      ### PQ ...
      # P=k+1,Q=k diagonal N=100..000 first of row is depth=P-2
      # anything else in that X=P column is smaller depth
      $depth = $x2 - 2;
    } else {
      my $xdepth = int (($x2+1) / 2);
      my $ydepth = int (($y2+31) / 4);
      $depth = min($xdepth,$ydepth);
    }
  }
  ### depth: "$depth"
  return (1, $self->tree_depth_to_n_end($zero+$depth));
}

#------------------------------------------------------------------------------
use constant tree_num_roots => 1;

sub tree_n_children {
  my ($self, $n) = @_;
  unless ($n >= 1) {
    return;
  }
  $n *= 3;
  return ($n-1, $n, $n+1);
}
sub tree_n_num_children {
  my ($self, $n) = @_;
  return ($n >= 1 ? 3 : undef);
}
sub tree_n_parent {
  my ($self, $n) = @_;
  unless ($n >= 2) {
    return undef;
  }
  return int(($n+1)/3);
}
sub tree_n_to_depth {
  my ($self, $n) = @_;
  ### PythagoreanTree tree_n_to_depth(): $n
  unless ($n >= 1) {
    return undef;
  }
  my ($pow, $depth) = round_down_pow (2*$n-1, 3);
  return $depth;
}

sub tree_depth_to_n {
  my ($self, $depth) = @_;
  return ($depth >= 0
          ? (3**$depth + 1)/2
          : undef);
}
# (3^(d+1)+1)/2-1 = (3^(d+1)-1)/2
sub tree_depth_to_n_end {
  my ($self, $depth) = @_;
  return ($depth >= 0
          ? (3**($depth+1) - 1)/2
          : undef);
}
sub tree_depth_to_n_range {
  my ($self, $depth) = @_;
  if ($depth >= 0) {
    my $n_lo = (3**$depth + 1) / 2;  # same as tree_depth_to_n()
    return ($n_lo, 3*$n_lo-2);
  } else {
    return;
  }
}
sub tree_depth_to_width {
  my ($self, $depth) = @_;
  return ($depth >= 0
          ? 3**$depth
          : undef);
}

#------------------------------------------------------------------------------

# Maybe, or abc_to_pq() perhaps with two of three values.
#
# @EXPORT_OK = ('ab_to_pq','pq_to_ab');
#
# =item C<($p,$q) = Math::PlanePath::PythagoreanTree::ab_to_pq($a,$b)>
#
# Return the P,Q coordinates for C<$a,$b>.  As described above this is
#
#     P = sqrt((C+A)/2)    where C=sqrt(A^2+B^2)
#     Q = sqrt((C-A)/2)
#
# The returned P,Q are integers PE=0,QE=0, but the further
# conditions for the path (namely PEQE=1 and no common factor) are
# not enforced.
#
# If P,Q are not integers or if BE0 then return an empty list.  This
# ensures A,B is a Pythagorean triple, ie. that C=sqrt(A^2+B^2) is an
# integer, but it might not be a primitive triple and might not have A odd B
# even.
#
# =item C<($a,$b) = Math::PlanePath::PythagoreanTree::pq_to_ab($p,$q)>
#
# Return the A,B coordinates for C<$p,$q>.  This is simply
#
#     $a = $p*$p - $q*$q
#     $b = 2*$p*$q
#
# This is intended for use with C<$p,$q> satisfying PEQE=1 and no
# common factor, but that's not enforced.


# a=p^2-q^2, b=2pq, c=p^2+q^2
# Done as a=(p-q)*(p+q) for one multiply instead of two squares, and to work
# close to a=UINT_MAX.
#
sub _pq_to_ab {
  my ($p, $q) = @_;
  return (($p-$q)*($p+$q), 2*$p*$q);
}

# C=(p-q)^2+B for one squaring instead of two.
# Also possible is C=(p+q)^2-B, but prefer "+B" so as not to round-off in
# floating point if (p+q)^2 overflows an integer.
sub _pq_to_bc {
  my ($p, $q) = @_;
  my $b = 2*$p*$q;
  $p -= $q;
  return ($b, $p*$p+$b);
}

# a=p^2-q^2, b=2pq, c=p^2+q^2
# Could a=(p-q)*(p+q) to avoid overflow if p^2 exceeds an integer as per
# _pq_to_ab(), but c overflows in that case anyway.
sub _pq_to_ac {
  my ($p, $q) = @_;
  $p *= $p;
  $q *= $q;
  return ($p-$q, $p+$q);
}

# a=p^2-q^2, b=2pq, c=p^2+q^2
# a
#  p^2-q^2 < 2pq
#  p^2 + 2pq - q^2 < 0
#  (p+q)^2 - 2*q^2 < 0
#  (p+q + sqrt(2)*q)*(p+q - sqrt(2)*q) < 0
#  (p+q - sqrt(2)*q) < 0
#  p + (1-sqrt(2))*q < 0
#  p < (sqrt(2)-1)*q
#
sub _pq_to_sc {
  my ($p, $q) = @_;
  my $b = 2*$p*$q;
  my $p_plus_q = $p + $q;
  $p -= $q;
  return (min($p_plus_q*$p, $b),  # A = P^2-Q^2 = (P+Q)*(P-Q)
          $p*$p+$b);              # C = P^2+Q^2 = (P-Q)^2 + 2*P*Q
}
sub _pq_to_mc {
  my ($p, $q) = @_;
  my $b = 2*$p*$q;
  my $p_plus_q = $p + $q;
  $p -= $q;
  return (max($p_plus_q*$p, $b),  # A = P^2-Q^2 = (P+Q)*(P-Q)
          $p*$p+$b);              # C = P^2+Q^2 = (P-Q)^2 + 2*P*Q
}
sub _pq_to_sm {
  my ($p, $q) = @_;
  my ($a, $b) = _pq_to_ab($p,$q);
  return ($a < $b ? ($a, $b) : ($b, $a));
}

# u = p+q, v=p-q
# at given p, vertical q
# u=p,v=p on diagonal then p+q,p-q is diagonal down
# so mirror p axis to x=y diagonal and measure down diagonal from there
sub _pq_to_uv {
  my ($p, $q) = @_;
  return ($p+$q, $p-$q);
}

# r = b+c = 2pq+p^2+q^2 = (p+q)^2
# s = c-a = p^2+q^2 - (p^2-q^2) = 2*q^2
sub _pq_to_rs {
  my ($p, $q) = @_;
  return (($p+$q)**2, 2*$q*$q);
}

# s = c-a = p^2+q^2 - (p^2-q^2) = 2*q^2
# t = a+b-c = p^2-q^2 + 2pq - (p^2+q^2) = 2pq-2q^2 = 2(p-q)q
sub _pq_to_st {
  my ($p, $q) = @_;
  my $q2 = 2*$q;
  return ($q2*$q, ($p-$q)*$q2);
}

#------------------------------------------------------------------------------

# a = p^2 - q^2
# b = 2pq
# c = p^2 + q^2
#
# q = b/2p
# a = p^2 - (b/2p)^2
#   = p^2 - b^2/4p^2
# 4ap^2 = 4p^4 - b^2
# 4(p^2)^2 - 4a*p^2 - b^2 = 0
# p^2 = [ 4a +/- sqrt(16a^2 + 16*b^2) ] / 2*4
#     = [ a +/- sqrt(a^2 + b^2) ] / 2
#     = (a +/- c) / 2   where c=sqrt(a^2+b^2)
# p = sqrt((a+c)/2)    since c>a
#
# a = (a+c)/2 - q^2
# q^2 = (a+c)/2 - a
#     = (c-a)/2
# q = sqrt((c-a)/2)
#
# if c^2 = a^2+b^2 is a perfect square then a,b,c is a pythagorean triple
# p^2 = (a+c)/2
#     = (a + sqrt(a^2+b^2))/2
# 2p^2 = a + sqrt(a^2+b^2)
#
# p>q so a>0
# a+c even is a odd, c odd or a even, c even
# if a odd then c=a^2+b^2 is opp of b parity, must have b even to make c+a even
# if a even then c=a^2+b^2 is same as b parity, must have b even to c+a even
#
# a=6,b=8 is c=sqrt(6^2+8^2)=10
# a=0,b=4 is c=sqrt(0+4^4)=4 p^2=(a+c)/2 = 2 not a square
# a+c even, then (a+c)/2 == 0,1 mod 4 so a+c==0,2 mod 4
#
sub _ab_to_pq {
  my ($a, $b) = @_;
  ### _ab_to_pq(): "A=$a, B=$b"

  unless ($b >= 4 && ($a%2) && !($b%2)) {   # A odd, B even
    return;
  }

  # This used to be $c=hypot($a,$b) and check $c==int($c), but libm hypot()
  # on Darwin 8.11.0 is somehow a couple of bits off being an integer, for
  # example hypot(57,176)==185 but a couple of bits out so $c!=int($c).
  # Would have thought hypot() ought to be exact on integer inputs and a
  # perfect square sum :-(.  Check for a perfect square by multiplying back
  # instead.
  #
  # The condition is "$csquared != $c*$c" with operands that way around
  # since the other way is bad for Math::BigInt::Lite 0.14.
  #
  my $c;
  {
    my $csquared = $a*$a + $b*$b;
    $c = _sqrtint($csquared);
    ### $csquared
    ### $c
    # since A odd and B even should have C odd, but floating point rounding
    # might prevent that
    unless ($csquared == $c*$c) {
      ### A^2+B^2 not a perfect square ...
      return;
    }
  }
  return _ac_to_pq($a,$c);
}

sub _bc_to_pqa {
  my ($b, $c) = @_;
  ### _bc_to_pqa(): "B=$b C=$c"

  unless ($c > $b && $b >= 4 && !($b%2) && ($c%2)) {  # B even, C odd
    return;
  }

  my $a;
  {
    my $asquared = $c*$c - $b*$b;
    unless ($asquared > 0) {
      return;
    }
    $a = _sqrtint($asquared);
    ### $asquared
    ### $a
    unless ($asquared == $a*$a) {
      return;
    }
  }

  # If $c is near DBL_MAX can have $a overflow to infinity, leaving A>C.
  # _ac_to_pq() will detect that.
  my ($p,$q) = _ac_to_pq($a,$c) or return;
  return ($p,$q,$a);
}

sub _ac_to_pq {
  my ($a, $c) = @_;
  ### _ac_to_pq(): "A=$a C=$c"

  unless ($c > $a && $a >= 3 && ($a%2) && ($c%2)) {  # A odd, C odd
    return;
  }
  $a = ($a-1)/2;
  $c = ($c-1)/2;
  ### halved to: "a=$a c=$c"

  my $p;
  {
    # If a,b,c is a triple but not primitive then can have psquared not an
    # integer.  Eg. a=9,b=12 has c=15 giving psquared=(9+15)/2=12 is not a
    # perfect square.  So notice that here.
    #
    my $psquared = $c+$a+1;
    $p = _sqrtint($psquared);
    ### $psquared
    ### $p
    unless ($psquared == $p*$p) {
      ### P^2=A+C not a perfect square ...
      return;
    }
  }

  my $q;
  {
    # If a,b,c is a triple but not primitive then can have qsquared not an
    # integer.  Eg. a=15,b=36 has c=39 giving qsquared=(39-15)/2=12 is not a
    # perfect square.  So notice that here.
    #
    my $qsquared = $c-$a;
    $q = _sqrtint($qsquared);
    ### $qsquared
    ### $q
    unless ($qsquared == $q*$q) {
      return;
    }
  }

  # Might have a common factor between P,Q here.  Eg.
  #     A=27 = 3*3*3, B=36 = 4*3*3
  #     A=45 = 3*3*5, B=108 = 4*3*3*3
  #     A=63, B=216
  #     A=75 =3*5*5  B=100 = 4*5*5
  #     A=81, B=360
  #
  return ($p, $q);
}

sub _sm_to_pq {
  my ($s, $m) = @_;
  unless ($s < $m) {
    return;
  }
  return _ab_to_pq($s % 2
                   ? ($s,$m)    # s odd is A
                   : ($m,$s));  # s even is B
}


# s^2+m^2=c^2
# if s odd then a=s
# ac_to_pq
# b = 2pq check isn't smaller than s
#
# p^2=(c+a)/2
# q^2=(c-a)/2

sub _sc_to_pq {
  my ($s, $c) = @_;
  my ($p,$q);
  if ($s % 2) {
    ($p,$q) = _ac_to_pq($s,$c)     # s odd is A
      or return;
    if ($s > 2*$p*$q) { return; }  # if s>B then s is not the smaller one
  } else {
    ($p,$q,$a) = _bc_to_pqa($s,$c)   # s even is B
      or return;
    if ($s > $a) { return; }         # if s>A then s is not the smaller one
  }
  return ($p,$q);
}

sub _mc_to_pq {
  my ($m, $c) = @_;
  ### _mc_to_pq() ...
  my ($p,$q);
  if ($m % 2) {
    ### m odd is A ...
    ($p,$q) = _ac_to_pq($m,$c)
      or return;
    if ($m < 2*$p*$q) { return; }   # if m
  } else {
    ### m even is B ...
    ($p,$q,$a) = _bc_to_pqa($m,$c)
      or return;
    ### $a
    if ($m < $a) { return; }         # if m
  }
  return ($p,$q);
}

# u = p+q, v=p-q
# u+v=2p   p = (u+v)/2
# u-v=2q   q = (u-v)/2
sub _uv_to_pq {
  my ($u, $v) = @_;
  return (($u+$v)/2, ($u-$v)/2);
}

# r = (p+q)^2
# s = 2*q^2 so   q = sqrt(r/2)
sub _rs_to_pq {
  my ($r, $s) = @_;

  return if $s % 2;
  $s /= 2;
  return unless $s >= 1;
  my $q = _sqrtint($s);
  return unless $q*$q == $s;

  return unless $r >= 1;
  my $p_plus_q = _sqrtint($r);
  return unless $p_plus_q*$p_plus_q == $r;

  return ($p_plus_q - $q, $q);
}

# s = 2*q^2
# t = a+b-c = p^2-q^2 + 2pq - (p^2+q^2) = 2pq-2q^2 = 2(p-q)q
#
# p=2,q=1  s=2  t=2.1.1=2
#
sub _st_to_pq {
  my ($s, $t) = @_;

  ### _st_to_pq(): "$s, $t"
  return if $s % 2;
  $s /= 2;
  return unless $s >= 1;
  my $q = _sqrtint($s);
  ### $q
  return unless $q*$q == $s;

  return if $t % 2;
  $t /= 2;
  ### rem: $t % $q
  return if $t % $q;
  $t /= $q;  # p-q

  ### pq: ($t+$q).", $q"

  return ($t+$q, $q);
}

1;
__END__