File Coverage

blib/lib/Math/PlanePath/CoprimeColumns.pm

Criterion	Covered	Total	%
statement	129	144	89.5
branch	38	48	79.1
condition	28	41	68.2
subroutine	24	25	96.0
pod	5	5	100.0
total	224	263	85.1

# Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019 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=CoprimeColumns --all --scale=10
# math-image --path=CoprimeColumns --output=numbers --all

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

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

use Math::PlanePath::Base::Generic
  'is_infinite',
  'round_nearest';

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


use constant default_n_start => 0;
use constant class_x_negative => 0;
use constant class_y_negative => 0;
use constant n_frac_discontinuity => .5;

use constant x_minimum => 1;
use constant y_minimum => 1;
               use constant diffxy_minimum => 0; # octant Y<=X so X-Y>=0
use constant gcdxy_maximum => 1;  # no common factor

use constant dx_minimum => 0;
use constant dx_maximum => 1;
use constant dir_maximum_dxdy => (1,-1); # South-East

use constant parameter_info_array =>
  [
   { name        => 'direction',
     share_key   => 'direction_updown',
     display     => 'Direction',
     type        => 'enum',
     default     => 'up',
     choices     => ['up','down'],
     choices_display => ['Down','Up'],
     description => 'Number points upwards or downwards in the columns.',
   },
   Math::PlanePath::Base::Generic::parameter_info_nstart0(),
  ];

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

sub new {
  my $self = shift->SUPER::new (@_);
  $self->{'direction'} ||= 'up';
  if (! defined $self->{'n_start'}) {
    $self->{'n_start'} = $self->default_n_start;
  }
  return $self;
}


# shared with DiagonalRationals
@_x_to_n = (0,0,1);
sub _extend {
  ### _extend(): $#_x_to_n
  my $x = $#_x_to_n;
  push @_x_to_n, $_x_to_n[$x] + _totient($x);

  # if ($x > 2) {
  #   if (($x & 3) == 2) {
  #     $x >>= 1;
  #     $next_n += $_x_to_n[$x] - $_x_to_n[$x-1];
  #   } else {
  #     $next_n +=
  #   }
  # }
  ### last x: $#_x_to_n
  ### second last: $_x_to_n[$#_x_to_n-2]
  ### last: $_x_to_n[$#_x_to_n-1]
  ### diff: $_x_to_n[$#_x_to_n-1] - $_x_to_n[$#_x_to_n-2]
  ### totient of: $#_x_to_n - 2
  ### totient: _totient($#_x_to_n-2)
  ### assert: $_x_to_n[$#_x_to_n-1] - $_x_to_n[$#_x_to_n-2] == _totient($#_x_to_n-2)
}

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

  $n = $n - $self->{'n_start'}; # to N=0 basis, and warn on undef

  # $n<-0.5 is ok for Math::BigInt circa Perl 5.12, it seems
  if (2*$n < -1) {
    return;
  }
  if (is_infinite($n)) {
    return ($n,$n);
  }

  my $frac;
  {
    my $int = int($n);
    $frac = $n - $int; # -.5 <= $frac < 1
    $n = $int;  # BigFloat int() gives BigInt, use that

    if (2*$frac >= 1) {
      $frac--;
      $n += 1;
      # now -.5 <= $frac < .5
    }
    ### $n
    ### $frac
    ### assert: 2*$frac >= -1
    ### assert: 2*$frac < 1
  }

  my $x = 1;
  for (;;) {
    while ($x > $#_x_to_n) {
      _extend();
    }
    if ($_x_to_n[$x] > $n) {
      $x--;
      last;
    }
    $x++;
  }
  $n -= $_x_to_n[$x];
  ### $x
  ### n base: $_x_to_n[$x]
  ### n next: $_x_to_n[$x+1]
  ### remainder: $n

  my $y = 1;
  for (;;) {
    if (_coprime($x,$y)) {
      if (--$n < 0) {
        last;
      }
    }
    if (++$y >= $x) {
      ### oops, not enough in this column ...
      return;
    }
  }

  $y += $frac;
  if ($x >= 2 && $self->{'direction'} eq 'down') {
    $y = $x - $y;
  }
  return ($x, $y);
}

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

sub xy_to_n {
  my ($self, $x, $y) = @_;
  ### CoprimeColumns xy_to_n(): "$x,$y"
  $x = round_nearest ($x);
  $y = round_nearest ($y);
  if (is_infinite($x)) { return $x; }
  if (is_infinite($y)) { return $y; }
  if ($x < 1
      || $y < 1
      || $y >= $x+($x==1)   # Y
      || ! _coprime($x,$y)) {
    return undef;
  }

  if ($x >= 2 && $self->{'direction'} eq 'down') {
    $y = $x - $y;
  }

  while ($#_x_to_n < $x) {
    _extend();
  }
  my $n = $_x_to_n[$x];
  ### base n: $n
  if ($y != 1) {
    foreach my $i (1 .. $y-1) {
      if (_coprime($x,$i)) {
        $n += 1;
      }
    }
  }
  return $n + $self->{'n_start'};
}

# Asymptotically
#     phisum(x) ~ 1/(2*zeta(2)) * x^2 + O(x ln x)
#               = 3/pi^2 * x^2 + O(x ln x)
# or by Walfisz
#     phisum(x) ~ 3/pi^2 * x^2 + O(x * (ln x)^(2/3) * (ln ln x)^4/3)
#
# but want an upper bound, so that for a given X at least enough N is
# covered ...
#
# Note: DiagonalRationals depends on this working only to column resolution.
# not exact
sub rect_to_n_range {
  my ($self, $x1,$y1, $x2,$y2) = @_;
  ### CoprimeColumns rect_to_n_range(): "$x1,$y1 $x2,$y2"

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

  if ($x2 < 1 || $y2 < 1
      # bottom right corner above X=Y diagonal, except X=1,Y=1 included
      || ($y1 >= $x2 + ($x2 == 1))) {
    ### outside ...
    return (1, 0);
  }
  if (is_infinite($x2)) {
    return ($self->{'n_start'}, $x2);
  }

  while ($#_x_to_n <= $x2) {
    _extend();
  }

  ### rect use xy_to_n at: "x=".($x2+1)." y=1"
  if ($x1 < 0) { $x1 = 0; }
  return ($_x_to_n[$x1]  + $self->{'n_start'},
          $_x_to_n[$x2+1] - 1 + $self->{'n_start'});

  # asympototically ?
  # return ($self->{'n_start'}, $self->{'n_start'} + .304*$x2*$x2 + 20);
}

# A000010
sub _totient {
  my ($x) = @_;
  my $count = (1                            # y=1 always
               + ($x > 2 && ($x&1))         # y=2 if $x odd
               + ($x > 3 && ($x % 3) != 0)  # y=3
               + ($x > 4 && ($x&1))         # y=4 if $x odd
              );
  for (my $y = 5; $y < $x; $y++) {
    $count += _coprime($x,$y);
  }
  return $count;
}

# code here only uses X>=Y but allow for any X,Y>=0 for elsewhere
sub _coprime {
  my ($x, $y) = @_;
  #### _coprime(): "$x,$y"

  if ($y > $x) {
    if ($x <= 1) {
      ### result yes ...
      return 1;
    }
    $y %= $x;
  }

  for (;;) {
    if ($y <= 1) {
      ### result: ($y == 1)
      return ($y == 1);
    }
    ($x,$y) = ($y, $x % $y);
  }
}

1;
__END__