File Coverage

blib/lib/Math/PlanePath/GcdRationals.pm

Criterion	Covered	Total	%
statement	137	159	86.1
branch	35	48	72.9
condition	12	23	52.1
subroutine	27	30	90.0
pod	7	7	100.0
total	218	267	81.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 .


# A003989 diagonals from (1,1)
# A109004        0,1,1,2,1,2,3,1,1,3,4,1,2,1,4,5,1,1,1,1
#  gcd by diagonals (0,0)=0
#                   (1,0)=1 (0,1)=1
#                   (2,0)=2 (1,1)=1 (0,2)=2
# A050873 gcd rows n>=1, k=1..n
#            1,1,2,1,1,3,1,2,1,4,1,1,1,1,5,1,2,3,2,1,6,1,1,1,
# add        0,1,0,1,1,0,1,1,1,0,1,1,1,1,0,1,1,1,1,1,0  A023532 0 at m(m+3)/2
# IntXY      1,0,2,0,0,3,0,1,0,4,0,0,0,0,5,0,1,2,1,0,6,
# IntXY+1    2,1,3,1,1,4,1,2,1,5,1,1,1,1,6,1,2,3,2,1,7
# diff       1,0,1,0,0,1,0,0,0,1,0,0,0,0,1,0,0,0,0,0,1   A023531
# A178340  1,2,1,3,1,1,4,1,2,1,5,1,1,1,1,6,1,2,3,2,1,7,1,1 Bernoulli
#   T(n,m) = A003989(n-m+1,m) m>=1, except when factor cancels

# diagonals_down even/odd in wedges, and other modulo

# math-image --path=GcdRationals --expression='i<30*31/2?i:0' --text --size=40
# math-image --path=GcdRationals --output=numbers --expression='i<100?i:0'
# math-image --path=GcdRationals --all --output=numbers

# Y = v = j/g
# X = (g-1)*v + u
#   = (g-1)*j/g + i/g
#   = ((g-1)*j + i)/g

# j=5  11 ...
# j=4  7 8 9 10
# j=3  4 5 6
# j=2  2 3
# j=1  1
#
# N = (1/2 d^2 - 1/2 d + 1)
#   = (1/2*$d**2 - 1/2*$d + 1)
#   = ((1/2*$d - 1/2)*$d + 1)
# j = 1/2 + sqrt(2 * $n + -7/4)
#   = [ 1 + 2*sqrt(2 * $n + -7/4) ] /2
#   = [ 1 + sqrt(8*$n -7) ] /2
#

# Primes
# i=3*a,j=3*b
# N=3*a*(3*b-1)/2


package Math::PlanePath::GcdRationals;
use 5.004;
use strict;
use Carp 'croak';
#use List::Util 'min','max';
*min = \&Math::PlanePath::_min;
*max = \&Math::PlanePath::_max;

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

use Math::PlanePath::Base::Generic
  'is_infinite',
  'round_nearest';
*_divrem = \&Math::PlanePath::_divrem;

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 => 1;
use constant gcdxy_maximum => 1;  # no common factor

use constant parameter_info_array =>
  [ { name        => 'pairs_order',
      display     => 'Pairs Order',
      type        => 'enum',
      default     => 'rows',
      choices     => ['rows','rows_reverse','diagonals_down','diagonals_up'],
      choices_display => ['Rows',
                          'Rows Reverse',
                          'Diagonals Down',
                          'Diagonals Up'],
      description => 'Order in the i,j pairs.',
    } ];

sub absdy_minimum {
  my ($self) = @_;
  return ($self->{'pairs_order'} eq 'diagonals_down'
          ? 1
          : 0);
}

{
  my %dir_minimum_dxdy
    = (rows           => [1,0],  # N=4 to N=5 horiz
       rows_reverse   => [1,0],  # N=1 to N=2 horiz
       diagonals_down => [0,1],  # N=1 to N=2 vertical, nothing less
       diagonals_up   => [1,0],  # N=4 to N=5 horiz
      );
  sub dir_minimum_dxdy {
    my ($self) = @_;
    return @{$dir_minimum_dxdy{$self->{'pairs_order'}}};
  }
}
{
  my %dir_maximum_dxdy
    = (rows           => [1,-1], # N=2 to N=3 SE diagonal
       rows_reverse   => [2,-1], # N=3 to N=4 dX=2,dY=-1
       diagonals_down => [1,-1], # N=5 to N=6 SE diagonal
       diagonals_up   => [2,-1], # N=9 to N=10 dX=2,dY=-1
      );
  sub dir_maximum_dxdy {
    my ($self) = @_;
    return @{$dir_maximum_dxdy{$self->{'pairs_order'}}};
  }
}

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

# all rationals X,Y >= 1 no common factor
use Math::PlanePath::DiagonalRationals;
*xy_is_visited = Math::PlanePath::DiagonalRationals->can('xy_is_visited');

sub new {
  my $self = shift->SUPER::new(@_);

  my $pairs_order = ($self->{'pairs_order'} ||= 'rows');
  (($self->{'pairs_order_n_to_xy'}
    = $self->can("_pairs_order__${pairs_order}__n_to_xy"))
   && ($self->{'pairs_order_xygr_to_n'}
       = $self->can("_pairs_order__${pairs_order}__xygr_to_n")))
    or croak "Unrecognised pairs_order: ",$pairs_order;

  return $self;
}

sub n_to_xy {
  my ($self, $n) = @_;
  ### GcdRationals 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;
      return ($frac*$dx + $x1, $frac*$dy + $y1);
    }
    $n = $int;
  }

  my ($x,$y) = $self->{'pairs_order_n_to_xy'}->($n);

  # if ($self->{'pairs_order'} eq 'rows'
  #     || $self->{'pairs_order'} eq 'rows_reverse') {
  #   $y = int((sqrt(8*$n-7) + 1) / 2);
  #   $x = $n - ($y - 1)*$y/2;
  #
  #   if ($self->{'pairs_order'} eq 'rows_reverse') {
  #     $x = $y - ($x-1);
  #   }
  #
  #   # require Math::PlanePath::PyramidRows;
  #   # my ($x,$y) = Math::PlanePath::PyramidRows->new(step=>1)->n_to_xy($n);
  #   # $x+=1;
  #   # $y+=1;
  #
  # } else {
  #   require Math::PlanePath::DiagonalsOctant;
  #   ($x,$y) = Math::PlanePath::DiagonalsOctant->new->n_to_xy($n);
  #   if ($self->{'pairs_order'} eq 'diagonals_up') {
  #     my $d = $x+$y;      # top 0,d measure diag down by x
  #     my $e = int($d/2);  # end e,d-e
  #     ($x,$y) = ($e-$x, $d - ($e-$x));
  #   }
  #   $x+=1;
  #   $y+=1;
  # }
  ### triangle: "$x,$y"

  my $gcd = _gcd($x,$y);
  $x /= $gcd;
  $y /= $gcd;

  ### $gcd
  ### reduced: "$x,$y"
  ### push out to x: $x + ($gcd-1)*$y

  return ($x + ($gcd-1)*$y, $y);
}

sub _pairs_order__rows__n_to_xy {
  my ($n) = @_;
  my $y = int( (_sqrtint(8*$n-7) + 1) / 2 );
  return ($n - ($y-1)*$y/2,
          $y);
}
sub _pairs_order__rows_reverse__n_to_xy {
  my ($n) = @_;
  my $y = int( (_sqrtint(8*$n-7) + 1) / 2 );
  return ($y*($y+1)/2 + 1 - $n,
          $y);
}
sub _pairs_order__diagonals_down__n_to_xy {
  my ($n) = @_;
  my $d = _sqrtint($n-1);   # eg. N=10 d=3
  $n -= $d*($d+1);          # eg. d=3 subtract 12
  if ($n > 0) {
    return ($n,
            2 - $n + 2*$d);
  } else {
    return ($n + $d,
            1 - $n + $d);
  }
}
sub _pairs_order__diagonals_up__n_to_xy {
  my ($n) = @_;
  my $d = _sqrtint($n-1);
  $n -= $d*($d+1);
  if ($n > 0) {
    return (-$n + $d + 2,
            $n + $d);
  } else {
    return (1 - $n,
            $n + 2*$d);
  }
}


# X=(g-1)*v+u
# Y=v
# u = x % y
# i = u*g
#   = (x % y)*g
#   = (x % y)*(floor(x/y)+1)
#
# Better:
#   g-1 = floor(x/y)
#   Y = j/g
#   X = ((g-1)*j + i)/g
#   j = Y*g
#   (g-1)*j + i = X*g
#   i = X*g - (g-1)*j
#     = X*g - (g-1)*Y*g
#   N = i + j*(j-1)/2
#     = X*g - (g-1)*Y*g + Y*g*(Y*g-1)/2
#     = X*g + Y*g * (-(g-1) + (Y*g-1)/2)    # but Y*g-1 may be odd
#     = X*g + Y*g * (Y*g-1 - (2g-2))/2
#     = X*g + Y*g * (Y*g-1 - 2g + 2))/2
#     = X*g + Y*g * (Y*g - 2g + 1))/2
#     = X*g + Y*g * ((Y-2)*g + 1) / 2
#     = g * [ X + Y*((Y-2)*g + 1) / 2 ]
#
#   N = X*g - (g-1)*Y*g + Y*g*(Y*g-1)/2
#     = [ 2*X*g - 2*(g-1)*Y*g + Y*g*(Y*g-1) ] / 2
#     = [ 2*X - 2*(g-1)*Y + Y*(Y*g-1) ] * g / 2
#     = [ 2*X + Y*(- 2*(g-1) + (Y*g-1)) ] * g / 2
#     = [ 2*X + Y*(-2g + 2 + Y*g - 1) ] * g / 2
#     = [ 2*X + Y*((Y-2)*g + 1) ] * g / 2
#     = X*g + [(Y-2)*g + 1]*Y*g/2
#
#  if Y and g both odd then (Y-2)*g+1 is odd+1 so even

# q=int(x/y)
# x = qy+r   qy=x-r
# r = x % y
# g-1 = q
# g = q+1
# g*y = (q+1)*y
#     = q*y + y
#     = x-r + y
#
#   N = X*g + Y*g * ((Y-2)*g + 1) / 2
#     = X*g + (X+Y-r) * ((Y-2)*g + 1) / 2
#     = X*g + (X+Y-r) * ((g*Y-2*g + 1) / 2
#     = X*g + (X+Y-r) * (((X+Y-r) - 2*g + 1) / 2
#     ... not much better

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

  if (is_infinite($x)) { return $x; }
  if (is_infinite($y)) { return $y; }
  if ($x < 1 || $y < 1 || ! _coprime($x,$y)) {
    return undef;
  }

  my ($p,$r) = _divrem ($x,$y);
  ### $x
  ### $y
  ### $p
  ### $r
  return $self->{'pairs_order_xygr_to_n'}->($x,$y,$p+1,$r);


  # my $g = int($x/$y) + 1;
  # ### g: "$g"
  # ### halve: ''.$y*(($y-2)*$g + 1)
  # return $self->{'pairs_order_xygr_to_n'}->($x,$y,$g);
}

sub _pairs_order__rows__xygr_to_n {
  my ($x,$y,$g,$r) = @_;
  ### j: $x+$y-$r
  ### i: $g*$r
  $x += $y;
  $x -= $r;  # j=X+Y-r
  return $x*($x-1)/2 + $g*$r;   # i=g*r
}

# i = X*g - (g-1)*g*Y
#   = X*g - (g-1)*(X+Y-r)
#   = X*g - g*(X+Y-r) + *(X+Y-r)
#   = X*g - g*X - g*Y + g*r + (X+Y-r)
#   = X*g - g*X - (X+Y-r) + g*r + (X+Y-r)
#   = g*r
#
# N = j-i+1 + j*(j-1)/2
#   = [2j-2i + 2 + $j*($j-1)] / 2
#   = [-2i + 2 + 2j+ j*(j-1)] / 2
#   = [-2i + 2 + j*(j-1+2)] / 2
#   = [-2i + 2 + j*(j+1)] / 2
#   = 1-i + j*(j+1)/2
#
sub _pairs_order__rows_reverse__xygr_to_n {
  my ($x,$y,$g,$r) = @_;
  $y += $x;
  $y -= $r;    # j = X+Y-r
  if ($r == 0) {
    # Case r=0 which is Y=1 becomes i=0 and that doesn't reverse to the
    # correct place by j-i+1.  Can either set $r=1,$g+=1 or leave $r==0
    # alone and adjust $y.
    $y -= 2;
  }
  return $y*($y+1)/2 - $r*$g + 1;
}

# d = (i-1)+(j-1)+1
#   = i+j-1
#   = rg + X+Y-r - 1
#   = X+Y + r*(g-1) - 1
# if r==0 Y==1 then r=1 g=X-1
# i = r*g = X-1
# j = X+Y-r = X+1-1 = X-1
# d = i+j-1
#   = 2X-2
# N = (d*d - (d%2))/4 + X-1
#   = ((2X-2)*(2X-2) - 0)/4 + X-1
#   = (X-1)^2 + X-1
#
sub _pairs_order__diagonals_down__xygr_to_n {
  my ($x,$y,$g,$r) = @_;

  $y += $x + $r*($g-1) - 1;   # d=X+Y + r*(g-1) - 1
  if ($r == 0) {
    $y *= 2;   # d=2*g-2
  }
  return ($y*$y - ($y % 2))/4 + $r*$g;
}
sub _pairs_order__diagonals_up__xygr_to_n {
  my ($x,$y,$g,$r) = @_;

  $y += $x + $r*($g-1);   # d=X+Y + r*(g-1)
  if ($r == 0) {
    $y = 2*$x - 1;
  }
  return ($y*$y - ($y % 2))/4 - $r*$g + 1;
}


# increase in rows, so right column
# in column increase within g wedge, then drop
#
# int(x2/y2) is slope of top of the wedge containing x2,y2
# g = int(x2/y2)+1 is the slope of the bottom of that wedge
# yw = floor(x2 / g) is the Y of that bottom
# N at x2,yw,g+1 is the top of the wedge underneath, bigger g smaller y
# or x2,y2,g is the top-right corner
#
# Eg.
# x=19 y=2 to 4
# g=int(19/4)+1=5
# yw=int(19/5)=3
# N(19,3,6)=
#
# at X=Y+1 g=2
# nhi = (y*((y-2)*g + 1) / 2 + x)*g
#     = (y*((y-2)*2 + 1) / 2 + y+1)*2
#     = (y*(2y-4 + 1) / 2 + y+1)*2
#     = (y*(2y-3) / 2 + y+1)*2
#     = y*(2y-3)  + 2y+2
#     = 2y^2 - 3y + 2y + 2
#     = 2y^2 - y + 2
#     = y*(2y-1) + 2

# 11  12  13  14      47  49  51  53     108 111 114 117     194 198 202 206
#  7       9      30      34      69      75     124     132     195     205
#  4   5      17  19      39  42      70  74     110 115     159 165     217
#  2       8      18      32      50      72      98     128     162     200
#  1   3   6  10  15  21  28  36  45  55  66  78  91 105 120 136 153 171 190

# 206=20*19/2+16  i=16,j=20 gcd=4
# 19,5 is slope=floor(19/5)=3 so g=4
#
# 205=20*19/2+15  i=15,j=20 gcd=5
# 19,4 is slope=floor(19/4)=4 so g=5
#
# 217=21*20/2 + 7, i=21,j=7  gcd=7
# 19,3 is slope=floor(19/3)=6 so g=7

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

  $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

  if ($x2 < 1 || $y2 < 1) {
    return (1, 0);  # outside quadrant
  }

  if ($x1 < 1) { $x1 = 1; }
  if ($y1 < 1) { $y1 = 1; }

  if ($self->{'pairs_order'} =~ /^diagonals/) {
    my $d = $x2 + max($x2,$y2);
    return (1, int($d*($d+($d%2)) / 4));  # N end of diagonal d
  }

  my $nhi;
  {
    my $c = max($x2,$y2);
    $nhi = _pairs_order__rows__xygr_to_n($c,$c,2,0);

    # my $rev = ($self->{'pairs_order'} eq 'rows_reverse');
    # my $slope = int($x2/$y2);
    # my $g = $slope + 1;
    #
    # # within top row
    # {
    #   my $x;
    #   if ($rev) {
    #     if ($slope > 0) {
    #       $x = max ($x1, $y2*$slope);  # left-most within this wedge
    #     } else {
    #       $x = $x1;  # top-left corner
    #     }
    #   } else {
    #     # pairs_order=rows
    #     $x = $x2;  # top-right corner
    #   }
    #   $nhi = $self->{'pairs_order_xygr_to_n'}->($x, $y2, $g, 0);
    #
    #   ### $slope
    #   ### $g
    #   ### x for hi: $x
    #   ### nhi for x,y2: $nhi
    # }
    #
    # # within x2 column, top of wedge below
    # #
    # my $yw = int(($x2+$g-1) / $g); # rounded up
    # if ($yw >= $y1) {
    #   $nhi = max ($nhi, $self->{'pairs_order_xygr_to_n'}->($x2,$yw,$g+1,0));
    #
    #   ### $yw
    #   ### nhi_wedge: $self->{'pairs_order_xygr_to_n'}->($x2,$yw,$g+1,0)
    # }
    #   my $yw = int($x2 / $g) - ($g==1);  # below X=Y diagonal when g==1
    #   if ($yw >= $y1) {
    #     $g = int($x2/$yw) + 1;  # perhaps went across more than one wedge
    #     $nhi = max ($nhi,
    #                 ($yw*(($yw-2)*($g+1) + 1) / 2 + $x2)*($g+1));
    #     ### $yw
    #     ### nhi_wedge: ($yw*(($yw-2)*($g+1) + 1) / 2 + $x2)*($g+1)
    #   }
  }

  my $nlo;
  {
    $nlo = _pairs_order__rows__xygr_to_n(1,$x1, 1, $x1-1);

    # my $g = int($x1/$y1) + 1;
    # $nlo = $self->{'pairs_order_xygr_to_n'}->($x1,$y1,$g,0);
    #
    # ### glo: $g
    # ### $nlo
    #
    # if ($g > 1) {
    #   my $yw = max (int($x1 / $g),
    #                 1);
    #   ### $yw
    #   if ($yw <= $y2) {
    #     $g = int($x1/$yw); # no +1, and perhaps up across more than one wedge
    #     $nlo = min ($nlo, $self->{'pairs_order_xygr_to_n'}->($x1,$yw,$g,0));
    #     ### glo_wedge: $g
    #     ### nlo_wedge: $self->{'pairs_order_xygr_to_n'}->($x1,$yw,$g,0)
    #   }
    # }
    # if ($nlo < 1) {
    #   $nlo = 1;
    # }
  }

  ### $nhi
  ### $nlo
  return ($nlo, $nhi);
}

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

  # bgcd() available in even the earliest Math::BigInt
  if ((ref $x && $x->isa('Math::BigInt'))
      || (ref $y && $y->isa('Math::BigInt'))) {
    return Math::BigInt::bgcd($x,$y);
  }

  $x = abs(int($x));
  $y = abs(int($y));
  unless ($x > 0) {
    return $y;  # gcd(0,y)=y for y>=0, giving gcd(0,0)=0
  }
  if ($y > $x) {
    $y %= $x;
  }
  for (;;) {
    ### assert: $x >= 1

    if ($y <= 1) {
      return ($y == 0
              ? $x   # gcd(x,0)=x
              : 1);  # gcd(x,1)=1
    }
    ($x,$y) = ($y, $x % $y);
  }
}



# # old code, rows only ...
# sub rect_to_n_range {
#   my ($self, $x1,$y1, $x2,$y2) = @_;
#   ### rect_to_n_range(): "$x1,$y1  $x2,$y2"
#
#   $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
#
#   if ($x2 < 1 || $y2 < 1) {
#     return (1, 0);  # outside quadrant
#   }
#
#   if ($x1 < 1) { $x1 = 1; }
#   if ($y1 < 1) { $y1 = 1; }
#
#   my $g = int($x2/$y2) + 1;
#   my $nhi = ($y2*(($y2-2)*$g + 1) / 2 + $x2)*$g;
#   ### ghi: $g
#   ### $nhi
#
#   my $yw = int($x2 / $g) - ($g==1);  # below X=Y diagonal when g==1
#   if ($yw >= $y1) {
#     $g = int($x2/$yw) + 1;  # perhaps went across more than one wedge
#     $nhi = max ($nhi,
#                 ($yw*(($yw-2)*($g+1) + 1) / 2 + $x2)*($g+1));
#     ### $yw
#     ### nhi_wedge: ($yw*(($yw-2)*($g+1) + 1) / 2 + $x2)*($g+1)
#   }
#
#   $g = int($x1/$y1) + 1;
#   my $nlo = ($y1*(($y1-2)*$g + 1) / 2 + $x1)*$g;
#
#   ### glo: $g
#   ### $nlo
#
#   if ($g > 1) {
#     $yw = max (int($x1 / $g),
#                1);
#     ### $yw
#     if ($yw <= $y2) {
#       $g = int($x1/$yw); # no +1, and perhaps up across more than one wedge
#       $nlo = min ($nlo,
#                   ($yw*(($yw-2)*$g + 1) / 2 + $x1)*$g);
#       ### glo_wedge: $g
#       ### nlo_wedge: ($yw*(($yw-2)*$g + 1) / 2 + $x1)*$g
#     }
#   }
#
#   return ($nlo, $nhi);
# }


1;
__END__