libmath-planepath-perl 113-1 / usr / share / perl5 / Math / PlanePath / FractionsTree.pm

# Copyright 2011, 2012, 2013 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/>.


# PowerPart has mostly square-free for X/Y > 1/2, then wedge of mostly
# multiple of 4, then mostly multiple of 16, then wedge of higher powers
# of 2.  Similar in AYT.



package Math::PlanePath::FractionsTree;
use 5.004;
use strict;

use vars '$VERSION', '@ISA';
$VERSION = 113;
use Math::PlanePath;
@ISA = ('Math::PlanePath');

use Math::PlanePath::Base::Generic
  'is_infinite',
  'round_nearest';
use Math::PlanePath::Base::Digits
  'bit_split_lowtohigh',
  'digit_join_lowtohigh';
use Math::PlanePath::RationalsTree;

use Math::PlanePath::CoprimeColumns;
*_coprime = \&Math::PlanePath::CoprimeColumns::_coprime;

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


use constant class_x_negative => 0;
use constant class_y_negative => 0;
use constant x_minimum => 1;
use constant y_minimum => 2;
use constant diffxy_maximum => -1; # upper octant X<=Y-1 so X-Y<=-1
use constant gcdxy_maximum => 1;  # no common factor
use constant tree_num_children_list => (2); # complete binary tree
use constant tree_n_to_subheight => undef; # complete tree, all infinity

use constant parameter_info_array =>
  [
   { name       => 'tree_type',
     share_key  => 'tree_type_fractionstree',
     display    => 'Tree Type',
     type       => 'enum',
     default    => 'Kepler',
     choices    => ['Kepler'],
   },
  ];

use constant dir_maximum_dxdy => (-2, -(sqrt(5)+1)); # phi


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

sub new {
  my $self = shift->SUPER::new(@_);
  $self->{'tree_type'} ||= 'Kepler';
  $self->{'n_start'} = 1; # for RationalsTree sharing
  return $self;
}

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

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

  # what to do for fractional $n?
  {
    my $int = int($n);
    if ($n != $int) {
      ### frac ...
      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;
      ### x1,y1: "$x1, $y1"
      ### x2,y2: "$x2, $y2"
      ### dx,dy: "$dx, $dy"
      ### result: ($frac*$dx + $x1).', '.($frac*$dy + $y1)
      return ($frac*$dx + $x1, $frac*$dy + $y1);
    }
    $n = $int;
  }

  my $zero = ($n * 0);  # inherit bignum 0
  my $one = $zero + 1;  # inherit bignum 1

  # my $tree_type = $self->{'tree_type'};
  # if ($tree_type eq 'Kepler')

  {
    ### Kepler tree ...

    #       X/Y
    #     /     \
    # X/(X+Y)  Y/(X+Y)
    #
    # (1 0) (x) = ( x )     (a b) (1 0) = (a+b b)   digit 0
    # (1 1) (y)   (x+y)     (c d) (1 1)   (c+d d)
    #
    # (0 1) (x) = ( y )     (a b) (0 1) = (b a+b)   digit 1
    # (1 1) (y)   (x+y)     (c d) (1 1)   (d c+d)

    my @bits = bit_split_lowtohigh($n);
    pop @bits;  # drop high 1 bit

    my $a = $one;     # initial  (1 0)
    my $b = $zero;    #          (0 1)
    my $c = $zero;
    my $d = $one;
    while (@bits) {
      ### at: "($a $b)"
      ### at: "($c $d)"
      ### $digit

      if (shift @bits) {      # low to high
        ($a,$b) = ($b, $a+$b);
        ($c,$d) = ($d, $c+$d);
      } else {
        $a += $b;
        $c += $d;
      }
    }
    ### final: "($a $b)"
    ### final: "($c $d)"

    # (a b) (1) = (a+b)
    # (c d) (2)   (c+d)
    return ($a+2*$b, $c+2*$d);
  }
}

sub xy_is_visited {
  my ($self, $x, $y) = @_;
  $x = round_nearest ($x);
  $y = round_nearest ($y);
  if ($x < 1 || $y < 2 || $x >= $y || ! _coprime($x,$y)) {
    return 0;
  }
  return 1;
}

sub xy_to_n {
  my ($self, $x, $y) = @_;
  $x = round_nearest ($x);
  $y = round_nearest ($y);
  ### FractionsTree xy_to_n(): "$x,$y   $self->{'tree_type'}"

  if ($x < 1 || $y < 2 || $x >= $y) {
    return undef;
  }
  if (is_infinite($x)) { return $x; }
  if (is_infinite($y)) { return $y; }

  my $zero = $x * 0 * $y;   # inherit bignum 0

  #       X/Y
  #     /     \
  # X/(X+Y)  Y/(X+Y)
  #
  # (x,y) <- (x, y-x)  nbit 0
  # (x,y) <- (y-x, x)  nbit 1
  #
  my @nbits;   # low to high
  for (;;) {
    ### at: "$x,$y n=$n"

    if ($y <= 2) {
      if ($x == 1 && $y == 2) {
        push @nbits, 1;  # high bit
        return digit_join_lowtohigh(\@nbits, 2, $zero);
      } else {
        return undef;
      }
    }

    ($y -= $x)          # (X,Y) <- (X, Y-X)
      || return undef;  # common factor if had X==Y
    if ($x > $y) {
      ($x,$y) = ($y,$x);
      push @nbits, 1;
    } else {
      push @nbits, 0;
    }
  }
}


# not exact
sub rect_to_n_range {
  my ($self, $x1,$y1, $x2,$y2) = @_;
  ### rect_to_n_range()

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

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


  #   |    /
  #   |   / x1
  #   |  /  +-----y2
  #   | /   |
  #   |/    +-----
  #
  if ($x2 < 1 || $y2 < 2 || $x1 >= $y2) {
    ### no values, rect outside upper octant ...
    return (1,0);
  }

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

  if ($x2 >= $y2) { $x2 = $y2-1; }
  if ($x1 < 1) { $x1 = 1; }
  if ($y1 < 2) { $y1 = 2; }

  # big x2, small y1
  # big y2, small x1
  # my $level = _bingcd_max ($y2,$x1);
  ### $level

  my $level = $y2;
  return (1, ($zero+2) ** $level);
}

sub _bingcd_max {
  my ($x,$y) = @_;
  ### _bingcd_max(): "$x,$y"

  if ($x < $y) { ($x,$y) = ($y,$x) }

  ### div: int($x/$y)
  ### bingcd: int($x/$y) + $y

  return int($x/$y) + $y + 1;
}

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

# Same structure as RationalsTree
*tree_n_children     = \&Math::PlanePath::RationalsTree::tree_n_children;
*tree_n_num_children = \&Math::PlanePath::RationalsTree::tree_n_num_children;
*tree_n_parent       = \&Math::PlanePath::RationalsTree::tree_n_parent;
*tree_n_to_depth     = \&Math::PlanePath::RationalsTree::tree_n_to_depth;
*tree_depth_to_n     = \&Math::PlanePath::RationalsTree::tree_depth_to_n;
*tree_depth_to_n_end = \&Math::PlanePath::RationalsTree::tree_depth_to_n_end;
*tree_depth_to_n_range=\&Math::PlanePath::RationalsTree::tree_depth_to_n_range;
*tree_depth_to_width = \&Math::PlanePath::RationalsTree::tree_depth_to_width;


1;
__END__

=for stopwords eg Ryde OEIS ie Math-PlanePath coprime Harmonices Mundi octant onwards Aiton

=head1 NAME

Math::PlanePath::FractionsTree -- fractions by tree

=head1 SYNOPSIS

 use Math::PlanePath::FractionsTree;
 my $path = Math::PlanePath::FractionsTree->new (tree_type => 'Kepler');
 my ($x, $y) = $path->n_to_xy (123);

=head1 DESCRIPTION

This path enumerates fractions X/Y in the range 0 E<lt> X/Y E<lt> 1 and in
reduced form, ie. X and Y having no common factor, using a method by
Johannes Kepler.

Fractions are traversed by rows of a binary tree which effectively
represents a coprime pair X,Y by subtraction steps of a subtraction-only
form of Euclid's greatest common divisor algorithm which would prove X,Y
coprime.  The steps left or right are encoded/decoded as an N value.

=head2 Kepler Tree

X<Kepler, Johannes>The default and only tree currently is by Johannes Kepler.

=over

"Harmonices Mundi", Book III.  Excerpt of translation by Aiton, Duncan and
Field at L<http://ndirty.cute.fi/~karttu/Kepler/a086592.htm>

=back

In principle similar bit reversal etc variations as in C<RationalsTree>
would be possible.

    N=1                             1/2
                              ------   ------
    N=2 to N=3             1/3               2/3
                          /    \            /   \
    N=4 to N=7         1/4      3/4      2/5      3/5
                       | |      | |      | |      | |
    N=8 to N=15     1/5  4/5  3/7 4/7  2/7 5/7  3/8 5/8

A node descends as

          X/Y
        /     \
    X/(X+Y)  Y/(X+Y)

Kepler described the tree as starting at 1, ie. 1/1, which descends to two
identical 1/2 and 1/2.  For the code here a single copy starting from 1/2 is
used.

Plotting the N values by X,Y is as follows.  Since it's only fractions
X/YE<lt>1, ie. XE<lt>Y, all points are above the X=Y diagonal.  The unused
X,Y positions are where X and Y have a common factor.  For example X=2,Y=6
have common factor 2 so is never reached.

    12  |    1024                  26        27                1025
    11  |     512   48   28   22   34   35   23   29   49  513     
    10  |     256        20                  21       257          
     9  |     128   24        18   19        25  129               
     8  |      64        14        15        65                    
     7  |      32   12   10   11   13   33                         
     6  |      16                  17                              
     5  |       8    6    7    9                                   
     4  |       4         5                                        
     3  |       2    3                                             
     2  |       1                                                  
     1  |
    Y=0 |   
         ----------------------------------------------------------
          X=0   1    2    3    4    5    6    7    8    9   10   11

The X=1 vertical is the fractions 1/Y at the left end of each tree row,
which is

    Nstart=2^level

The diagonal X=Y-1, fraction K/(K+1), is the second in each row, at
N=Nstart+1.  That's the maximum X/Y in each level.

The N values in the upper octant, ie. above the line Y=2*X, are even and
those below that line are odd.  This arises since XE<lt>Y so the left leg
X/(X+Y) E<lt> 1/2 and the right leg Y/(X+Y) E<gt> 1/2.  The left is an even
N and the right an odd N.

=head1 FUNCTIONS

See L<Math::PlanePath/FUNCTIONS> for behaviour common to all path classes.

=over

=item C<$path = Math::PlanePath::FractionsTree-E<gt>new ()>

Create and return a new path object.

=item C<$n = $path-E<gt>n_start()>

Return 1, the first N in the path.

=item C<($n_lo, $n_hi) = $path-E<gt>rect_to_n_range ($x1,$y1, $x2,$y2)>

Return a range of N values which occur in a rectangle with corners at
C<$x1>,C<$y1> and C<$x2>,C<$y2>.  The range is inclusive.

For reference, C<$n_hi> can be quite large because within each row there's
only one new 1/Y fraction.  So if X=1 is included then roughly C<$n_hi =
2**max(x,y)>.

=back

=head2 Tree Methods

X<Complete binary tree>Each point has 2 children, so the path is a complete
binary tree.

=over

=item C<@n_children = $path-E<gt>tree_n_children($n)>

Return the two children of C<$n>, or an empty list if C<$n E<lt> 1>
(before the start of the path).

This is simply C<2*$n, 2*$n+1>.  The children are C<$n> with an extra bit
appended, either a 0-bit or a 1-bit.

=item C<$num = $path-E<gt>tree_n_num_children($n)>

Return 2, since every node has two children, or return C<undef> if
C<$nE<lt>1> (before the start of the path).

=item C<$n_parent = $path-E<gt>tree_n_parent($n)>

Return the parent node of C<$n>, or C<undef> if C<$n E<lt>= 1> (the top of
the tree).

This is simply C<floor($n/2)>, stripping the least significant bit from
C<$n> (undoing what C<tree_n_children()> appends).

=item C<$depth = $path-E<gt>tree_n_to_depth($n)>

Return the depth of node C<$n>, or C<undef> if there's no point C<$n>.  The
top of the tree at N=1 is depth=0, then its children depth=1, etc.

The structure of the tree with 2 nodes per point means the depth is simply
floor(log2(N)), so for example N=4 through N=7 are all depth=2.

=back

=head2 Tree Descriptive Methods

=over

=item C<$num = $path-E<gt>tree_num_children_minimum()>

=item C<$num = $path-E<gt>tree_num_children_maximum()>

Return 2 since every node has 2 children, making that both the minimum and
maximum.

=item C<$bool = $path-E<gt>tree_any_leaf()>

Return false, since there are no leaf nodes in the tree.

=back

=head1 OEIS

The trees are in Sloane's Online Encyclopedia of Integer Sequences in the
following forms

=over

L<http://oeis.org/A020651> (etc)

=back

    tree_type=Kepler
      A020651  - X numerator (RationalsTree AYT denominators)
      A086592  - Y denominators
      A086593  - X+Y sum, and every second denominator
      A020650  - Y-X difference (RationalsTree AYT numerators)

The tree descends as X/(X+Y) and Y/(X+Y) so the denominators are two copies
of X+Y time after the initial 1/2.  A086593 is every second, starting at 2,
eliminating the duplication.  This is also the sum X+Y, from value 3
onwards, as can be seen by thinking of writing a node as the X+Y which would
be the denominators it descends to.

=head1 SEE ALSO

L<Math::PlanePath>,
L<Math::PlanePath::RationalsTree>,
L<Math::PlanePath::CoprimeColumns>,
L<Math::PlanePath::PythagoreanTree>

L<Math::NumSeq::SternDiatomic>,
L<Math::ContinuedFraction>

=head1 HOME PAGE

L<http://user42.tuxfamily.org/math-planepath/index.html>

=head1 LICENSE

Copyright 2011, 2012, 2013 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/>.

=cut