libmath-planepath-perl 117-1 / usr / share / perl5 / Math / PlanePath / CfracDigits.pm

# Copyright 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/>.


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

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

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::RationalsTree;
*_xy_to_quotients = \&Math::PlanePath::RationalsTree::_xy_to_quotients;

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

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


use constant parameter_info_array =>
  [ { name      => 'radix',
      share_key => 'radix2_min1',
      display   => 'Radix',
      type      => 'integer',
      minimum   => 1,
      default   => 2,
      width     => 3,
    },
  ];

use constant n_start => 0;
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

# FIXME: believe this is right, but check N+1 always changes Y
sub absdy_minimum {
  my ($self) = @_;
  return ($self->{'radix'} < 3 ? 0 : 1);
}

# radix=1 N=1       has dir4=0
# radix=2 N=5628    has dir4=0 dx=9,dy=0
# radix=3 N=1189140 has dir4=0 dx=1,dy=0
# radix=4 N=169405  has dir4=0 dx=2,dy=0
# always eventually 0 ?
# use constant dir_minimum_dxdy => (1,0);  # the default

# radix=1 N=4    dX=1,dY=-1 for dir4=3.5
# radix=2 N=4413 dX=9,dY=-1
# radix=3 N=9492 dX=3,dY=-1
# ENHANCE-ME: suspect believe approaches 360 degrees, eventually, but proof?
# use constant dir_maximum_dxdy => (0,0);  # the default


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

sub new {
  my $self = shift->SUPER::new (@_);
  unless ($self->{'radix'} && $self->{'radix'} >= 1) {
    $self->{'radix'} = 2;
  }
  return $self;
}

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

  if ($n < 0) { return; }
  if (is_infinite($n)) { return ($n,$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 $radix = $self->{'radix'};
  my $zero = ($n * 0);  # inherit bignum 0
  my $x = $zero;
  my $y = 1 + $zero;    # inherit bignum 1

  foreach my $q (_n_to_quotients_bottomtotop($n,$radix,$zero)) {  # bottom to top
    ### at: "$x,$y   q=$q"

    # 1/(q + X/Y) = 1/((qY+X)/Y)
    #             = Y/(qY+X)
    ($x,$y) = ($y, $q*$y + $x);
  }

  ### return: "$x,$y"
  return ($x,$y);
}

# Return a list of quotients bottom to top.  The base3 digits of N are split
# by "3" delimiters and the parts adjusted so the first bottom-most q>=2 and
# the rest q>=1.  The values are ready to be used as continued fraction
# terms.
#
sub _n_to_quotients_bottomtotop {
  my ($n, $radix, $zero) = @_;
  ### _n_to_quotients_bottomtotop(): $n

  my $radix_plus_1 = $radix + 1;
  my @ret;
  my @group;
  foreach my $digit (_digit_split_1toR_lowtohigh($n,$radix_plus_1)) {
    if ($digit == $radix_plus_1) {
      ### @group
      push @ret, _digit_join_1toR_destructive(\@group, $radix, $zero) + 1;
      @group = ();
    } else {
      push @group, $digit;
    }
  }
  ### final group: @group
  push @ret, _digit_join_1toR_destructive(\@group, $radix, $zero) + 1;

  $ret[0] += 1;  # bottom-most is +2 rather than +1

  ### _n_to_quotients_bottomtotop result: @ret
  return @ret;
}

# Return a list of digits 1 <= d <= R which is $n written in $radix, low to
# high digits.
sub _digit_split_1toR_lowtohigh {
  my ($n, $radix) = @_;
  ### assert: $radix >= 1
  ### assert: $n >= 0

  if ($radix == 1) {
    return (1) x $n;
  }
  my @digits = digit_split_lowtohigh($n,$radix);

  # mangle 0 -> R
  my $borrow = 0;
  foreach my $digit (@digits) {   # low to high
    if ($borrow = (($digit -= $borrow) <= 0)) {  # modify array contents
      $digit += $radix;
    }
  }
  if ($borrow) {
    ### assert: $digits[-1] == $radix
    pop @digits;
  }

  return @digits;
}

sub _digit_join_1toR_destructive {
  my ($aref, $radix, $zero) = @_;
  ### assert: $radix >= 1

  if ($radix == 1) {
    return scalar(@$aref);
  }

  # mangle any digit==$radix down to digit=0
  my $carry = 0;
  foreach my $digit (@$aref) {   # low to high
    if ($carry = (($digit += $carry) >= $radix)) {  # modify array contents
      $digit -= $radix;
    }
  }
  if ($carry) {
    push @$aref, 1;
  }

  ### _digit_join_1toR_destructive() result: digit_join_lowtohigh($aref, $radix, $zero)
  return digit_join_lowtohigh($aref, $radix, $zero);
}

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

sub xy_to_n {
  my ($self, $x, $y) = @_;
  $x = round_nearest ($x);
  $y = round_nearest ($y);
  ### CfracDigits xy_to_n(): "$x,$y"

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

  my @quotients = _xy_to_quotients($x,$y)
    or return undef;  # $x,$y have a common factor
  ### @quotients

  # drop initial 0 integer part
  ### assert: $quotients[0] == 0
  shift @quotients;

  return _cfrac_join_toptobottom(\@quotients,
                                 $self->{'radix'},
                                 $x * 0 * $y);   # inherit bignum 0
}

# $aref is a list of continued fraction quotients from top-most to
# bottom-most.  There's no initial integer term in $aref.  Each quotient is
# q >= 1 except the bottom-most which q-1 and so also >=1.
#
sub _cfrac_join_toptobottom {
  my ($aref, $radix, $zero) = @_;
  ### _cfrac_join_toptobottom(): $aref

  my @digits;
  foreach my $q (reverse @$aref) {
    ### assert: $q >= 1
    push @digits, _digit_split_1toR_lowtohigh($q-1, $radix), $radix+1;
  }
  pop @digits;  # no high delimiter
  ### groups digits 1toR: @digits
  return _digit_join_1toR_destructive(\@digits, $radix+1, $zero);
}


# X/Y = F[k]/F[k+1] quotients all 1
# N = all delimiter digits R,R,...,R
#   = 1222...2221
#   = R^k + 2*(R^k+1)/(R-1) - 1
#   = (RR^k - R^k + 2R^k + 2 - R + 1) / (R-1)
#   = (RR^k + R^k - R + 3) / (R-1)
#   = ((R+1)R^k - R + 3) / (R-1)
# take high as "12" = R+2
# k = log(Y)/log(phi)
# N = (R+2) * R ** k
# N = Y ** (log(R)/log(phi))
#
# 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
  my $radix = $self->{'radix'};

  return (0,
          ($radix+3)
          * ($radix+1 + $zero) ** ($radix == 1
                                   ? $y2
                                   : _log_phi_estimate($y2,$radix)));
}

# Return an estimate of log base phi of $x, that being log($x)/log(phi),
# where phi=(1+sqrt(5))/2 the golden ratio.
#
sub _log_phi_estimate {
  my ($x) = @_;
  my ($pow,$exp) = round_down_pow ($x, 2);
  return int ($exp * (log(2) / log((1+sqrt(5))/2)));
}

1;
__END__

=for stopwords eg Ryde OEIS ie Math-PlanePath coprime octant onwards decrement Shallit radix-1 Radix radix HCS 10www w's

=head1 NAME

Math::PlanePath::CfracDigits -- continued fraction terms encoded by digits

=head1 SYNOPSIS

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

=head1 DESCRIPTION

X<Shallit, Jeffrey>This path enumerates reduced fractions
S<0 E<lt> X/Y E<lt> 1> with X,Y no common factor using a method by Jeffrey
Shallit encoding continued fraction terms in digit strings, as per

=over

Jeffrey Shallit, "Number Theory and Formal Languages", part 3,
L<https://cs.uwaterloo.ca/~shallit/Papers/ntfl.ps>

=back

Fractions up to a given denominator are covered by roughly N=den^2.28.  This
is a much smaller N range than the run-length encoding in C<RationalsTree>
and C<FractionsTree> (but is more than C<GcdRationals>).

=cut

# math-image --path=CfracDigits --output=numbers_xy --all --size=78x17

=pod

    15  |    25  27      91          61 115         307     105 104
    14  |    23      48      65             119     111     103
    13  |    22  24  46  29  66  59 113 120 101 109  99  98
    12  |    17              60     114              97
    11  |    16  18  30  64  58 112 118 102  96  95
    10  |    14      28             100      94
     9  |    13  15      20  38      36  35
     8  |     8      21      39      34
     7  |     7   9  19  37  33  32
     6  |     5              31
     5  |     4   6  12  11
     4  |     2      10
     3  |     1   3
     2  |     0
     1  |
    Y=0 |
         ----------------------------------------------------------
        X=0   1   2   3   4   5   6   7   8   9  10  11  12  13  14

A fraction S<0 E<lt> X/Y E<lt> 1> has a finite continued fraction of the
form

                      1
    X/Y = 0 + ---------------------
                            1
              q[1] + -----------------
                                  1
                     q[2] + ------------
                         ....
                                      1
                            q[k-1] + ----
                                     q[k]

    where each  q[i] >= 1
    except last q[k] >= 2


The terms are collected up as a sequence of integers E<gt>=0 by subtracting
1 from each and 2 from the last.

    # each >= 0
    q[1]-1,  q[2]-1, ..., q[k-2]-1, q[k-1]-1, q[k]-2

These integers are written in base-2 using digits 1,2.  A digit 3 is written
between each term as a separator.

    base2(q[1]-1), 3, base2(q[2]-1), 3, ..., 3, base2(q[k]-2)

If a term q[i]-1 is zero then its base-2 form is empty and there's adjacent
3s in that case.  If the high q[1]-1 is zero then a bare high 3, and if the
last q[k]-2 is zero then a bare final 3.  If there's just a single term q[1]
and q[1]-2=0 then the string is completely empty.  This occurs for X/Y=1/2.

The resulting string of 1s,2s,3s is reckoned as a base-3 value with digits
1,2,3 and the result is N.  All possible strings of 1s,2s,3s occur
(including the empty string) and so all integers NE<gt>=0 correspond
one-to-one with an X/Y fraction with no common factor.

Digits 1,2 in base-2 means writing an integer in the form

    d[k]*2^k + d[k-1]*2^(k-1) + ... + d[2]*2^2 + d[1]*2 + d[0]
    where each digit d[i]=1 or 2

Similarly digits 1,2,3 in base-3 which is used for N,

    d[k]*3^k + d[k-1]*3^(k-1) + ... + d[2]*3^2 + d[1]*3 + d[0]
    where each digit d[i]=1, 2 or 3

This is not the same as the conventional binary and ternary radix
representations by digits 0,1 or 0,1,2 (ie. 0 to radix-1).  The effect of
digits 1 to R is to change any 0 digit to instead R and decrement the value
above that position to compensate.

=head2 Axis Values

N=0,1,2,4,5,7,etc in the X=1 column is integers with no digit 0s in ternary.
N=0 is considered no digits at all and so no digit 0.  These points are
fractions 1/Y which are a single term q[1]=Y-1 and hence no "3" separators,
only a run of digits 1,2.  These N values are also those which are the same
when written in digits 0,1,2 as when written in digits 1,2,3, since there's
no 0s or 3s.

N=0,3,10,11,31,etc along the diagonal Y=X+1 are integers which are ternary
"10www..." where the w's are digits 1 or 2, so no digit 0s except the
initial "10".  These points Y=X+1 points are X/(X+1) with continued fraction

                     1
    X/(X+1) =  0 + -------
                        1
                   1 + ---
                        X

so q0=1 and q1=X, giving N="3,X-1" in digits 1,2,3, which is N="1,0,X-1" in
normal ternary.  For example N=34 is ternary "1021" which is leading "10"
and then X-1=7 ternary "21".

=head2 Radix

The optional C<radix> parameter can select another base for the continued
fraction terms, and corresponding radix+1 for the resulting N.  The default
is radix=2 as described above.  Any integer radixE<gt>=1 can be selected.
For example,

=cut

# math-image --path=CfracDigits,radix=5 --output=numbers_xy --all --size=78x17

=pod

    radix => 5

    11  |    10   30  114  469   75  255 1549 1374  240  225
    10  |     9       109                1369       224
     9  |     8   24        74  254       234  223
     8  |     7        78       258        41
     7  |     5   18   73  253  228   40
     6  |     4                  39
     5  |     3   12   42   38
     4  |     2        37
     3  |     1    6
     2  |     0
     1  |
    Y=0 |
         ----------------------------------------------------
        X=0   1    2    3    4    5    6    7    8    9   10

The X=1 column is integers with no digit 0 in base radix+1, so in radix=5
means no 0 digit in base-6.

=head2 Radix 1

The radix=1 case encodes continued fraction terms using only digit 1, which
means runs of q many "1"s to add up to q, and then digit "2" as separator.

    N =  11111 2 1111 2 ... 2 1111 2 11111     base2 digits 1,2
         \---/   \--/         \--/   \---/
         q[1]-1  q[2]-1     q[k-1]-1 q[k]-2

which becomes in plain binary

    N = 100000  10000   ...  10000  011111     base2 digits 0,1
        \----/  \---/        \---/  \----/
         q[1]    q[2]       q[k-1]  q[k]-1

Each "2" becomes "0" in plain binary and carry +1 into the run of 1s above
it.  That carry propagates through those 1s, turning them into 0s, and stops
at the "0" above them (which had been a "2").  The low run of 1s from q[k]-2
has no "2" below it and is therefore unchanged.

=cut

# math-image --path=CfracDigits,radix=1 --output=numbers_xy --all --size=60x12

=pod

    radix => 1

    11  |   511  32  18  21  39  55  29  26  48 767
    10  |   255      17              25     383
     9  |   127  16      19  27      24 191
     8  |    63      10      14      95
     7  |    31   8   9  13  12  47
     6  |    15              23
     5  |     7   4   6  11
     4  |     3       5
     3  |     1   2
     2  |     0
     1  |
    Y=0 |
         -------------------------------------------
        X=0   1   2   3   4   5   6   7   8   9  10

The result is similar to L<Math::PlanePath::RationalsTree/HCS Continued
Fraction>.  But the lowest run is "0111" here, instead of "1000" as it is in
the HCS.  So N-1 here, and a flip (Y-X)/X to map from X/YE<lt>1 here to
instead all rationals for the HCS tree.  For example

    CfracDigits radix=1       RationalsTree tree_type=HCS

    X/Y = 5/6                 (Y-X)/X = 1/5
    is at                     is at
    N = 23 = 0b10111          N = 24 = 0b11000
                ^^^^                      ^^^^

=head1 FUNCTIONS

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

=over

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

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

Create and return a new path object.

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

Return 0, the first N in the path.

=back

=head1 OEIS

Entries in Sloane's Online Encyclopedia of Integer Sequences related to this
path include

=over

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

=back

    radix=1
      A071766    X coordinate (numerator), except extra initial 1

    radix=2 (the default)
      A032924    N in X=1 column, ternary no digit 0 (but lacking N=0)

    radix=3
      A023705    N in X=1 column, base-4 no digit 0 (but lacking N=0)

    radix=4
      A023721    N in X=1 column, base-5 no digit 0 (but lacking N=0)

    radix=10
      A052382    N in X=1 column, decimal no digit 0 (but lacking N=0)

=head1 SEE ALSO

L<Math::PlanePath>,
L<Math::PlanePath::FractionsTree>,
L<Math::PlanePath::CoprimeColumns>

L<Math::PlanePath::RationalsTree>,
L<Math::PlanePath::GcdRationals>,
L<Math::PlanePath::DiagonalRationals>

L<Math::ContinuedFraction>

=head1 HOME PAGE

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

=head1 LICENSE

Copyright 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/>.

=cut