File Coverage

blib/lib/Math/PlanePath/DivisibleColumns.pm

Criterion	Covered	Total	%
statement	112	172	65.1
branch	26	66	39.3
condition	7	39	17.9
subroutine	20	24	83.3
pod	7	7	100.0
total	172	308	55.8

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


# A006218 - cumulative count of divisors
#
#   Dirichlet:
#   n * (log(n) + 2*gamma - 1) + O(sqrt(n))  gamma=0.57721... Euler-Mascheroni
#
#   n * (log(n) + 2*gamma - 1) + O(log(n)*n^(1/3))
#
#   Chandrasekharan: bounds
#   n log(n) + (2 gamma - 1) n - 4 sqrt(n) - 1
#   <= a(n) <=
#   n log(n) + (2 gamma - 1) n + 4 sqrt(n)
#
# a(n)=2 * sum[ i=1 to floor(sqrt(n)) of floor(n/i) ] - floor(sqrt(n))^2
#
# cf A003988,A010766 - triangle with values floor(i/j)
#
# http://mathworld.wolfram.com/DirichletDivisorProblem.html
#
# compile-command: "math-image --path=DivisibleColumns --all"
#
# math-image --path=DivisibleColumns --output=numbers --all
#

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

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

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

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

use constant parameter_info_array =>
  [ { name      => 'divisor_type',
      share_key => 'divisor_type_allproper',
      display   => 'Divisor Type',
      type      => 'enum',
      choices   => ['all','proper'],
      default   => 'all',
      description => 'Divisor type, with "proper" meaning divisors d
    },
    # { 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(),
  ];

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

# X=2,Y=1 when proper
# X=1,Y=1 when not
sub x_minimum {
  my ($self) = @_;
  return ($self->{'proper'} ? 2 : 1);
}
use constant y_minimum => 1;

sub diffxy_minimum {
  my ($self) = @_;
  # octant Y<=X so X-Y>=0
  return ($self->{'proper'} ? 1 : 0);
}

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


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

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

  my $divisor_type = ($self->{'divisor_type'} ||= 'all');
  $self->{'proper'} = ($divisor_type eq 'proper');  # bool

  $self->{'direction'} ||= 'up';

  if (! defined $self->{'n_start'}) {
    $self->{'n_start'} = $self->default_n_start;
  }
  return $self;
}

my @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] + _count_divisors($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]
  ### divisors of: $#x_to_n - 2
  ### divisors: _count_divisors($#x_to_n-2)
  ### assert: $x_to_n[$#x_to_n-1] - $x_to_n[$#x_to_n-2] == _count_divisors($#x_to_n-2)
}

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

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

  # $n<-0.5 works with Math::BigInt circa Perl 5.12, it seems
  if ($n < -0.5) {
    return;
  }
  if (is_infinite($n)) {
    return ($n,$n);
  }

  my $frac;
  {
    my $int = int($n);
    if ($n == $int) {
      $frac = 0;
    } else {
      $frac = $n - $int; # -.5 <= $frac < 1
      $n = $int;         # BigFloat int() gives BigInt, use that
      if ($frac > .5) {
        $frac--;
        $n += 1;
        # now -.5 <= $frac < .5
      }
    }
    ### $n
    ### n: "$n"
    ### $frac
    ### assert: $frac >= -.5
    ### assert: $frac < .5
  }
  my $proper = $self->{'proper'} || 0;  # cannot add false '' to BigInt
  ### $proper

  my $x;
  if ($proper) {
    $x = 2;
    ### proper adjusted n: $n
  } else {
    $x = 1;
  }

  for (;;) {
    while ($x > $#x_to_n) {
      _extend();
    }
    $n += $proper;
    ### consider: "n=$n x=$x  x_to_n=".$x_to_n[$x]
    if ($x_to_n[$x] > $n) {
      $x--;
      last;
    }
    $x++;
  }
  $n -= $x_to_n[$x];
  $n -= $proper;
  ### $x
  ### x_to_n: $x_to_n[$x]
  ### x_to_n next: $x_to_n[$x+1]
  ### remainder: $n

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

# Feturn a count of the number of integers dividing $x, including 1 and $x
# itself.   Cf Math::Factor::XS maybe.
sub _count_divisors {
  my ($x) = @_;
  my $ret = 1;
  unless ($x % 2) {
    my $count = 1;
    do {
      $x /= 2;
      $count++;
    } until ($x % 2);
    $ret *= $count;
  }
  my $limit = _sqrtint($x);
  for (my $d = 3; $d <= $limit; $d+=2) {
    unless ($x % $d) {
      my $count = 1;
      do {
        $x /= $d;
        $count++;
      } until ($x % $d);
      $limit = sqrt($x);
      $ret *= $count;
    }
  }
  if ($x > 1) {
    $ret *= 2;
  }
  return $ret;
}

sub xy_is_visited {
  my ($self, $x, $y) = @_;
  $x = round_nearest ($x);
  $y = round_nearest ($y);
  return ($y >= 1
          && ($self->{'proper'}
              ? $x >= 2 && $y <= int($x/2)
              : $x >= 1 && $y <= $x)
          && ($x%$y) == 0);
}

sub xy_to_n {
  my ($self, $x, $y) = @_;
  ### DivisibleColumns 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; }

  my $proper = $self->{'proper'};
  if ($proper) {
    if ($x < 2
        || $y < 1
        || $y > int($x/2)
        || ($x%$y)) {
      return undef;
    }
  } else {
    if ($x < 1
        || $y < 1
        || $y > $x
        || ($x%$y)) {
      return undef;
    }
  }

  while ($#x_to_n < $x) {
    _extend();
  }
  ### x_to_n: $x_to_n[$x]

  my $n = $x_to_n[$x] - ($proper ? $x-1 : 1);
  ### base n: $n

  for (my $i = 1+$proper; $i <= $y; $i++) {
    unless ($x % $i) {
      $n += 1;
    }
  }
  return $n + $self->{'n_start'};
}

# not exact
sub rect_to_n_range {
  my ($self, $x1,$y1, $x2,$y2) = @_;
  ### DivisibleColumns 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;

  ### rounded ...
  ### $x2
  ### $y2

  if ($self->{'proper'}) {
    if ($x2 < 2            # rect all negative
        || $y2 < 1         # rect all negative
        || 2*$y1 > $x2) {  # rect all above X=2Y octant
      ### outside proper divisors ...
      return (1, 0);
    }
    if ($x1 < 2) { $x1 = 2; }
  } else {
    if ($x2 < 1           # rect all negative
        || $y2 < 1        # rect all negative
        || $y1 > $x2) {   # rect all above X=Y diagonal
      ### outside all divisors ...
      return (1, 0);
    }
    if ($x1 < 1) { $x1 = 1; }
  }
  if (is_infinite($x2)) {
    return ($self->{'n_start'}, $x2);
  }

  my ($n_lo, $n_hi);
  if ($x1 <= $#x_to_n) {
    $n_lo = $x_to_n[$x1];
  } else {
    $n_lo = _count_divisors_cumulative($x1-1);
  }
  if ($x2 < $#x_to_n) {
    $n_hi = $x_to_n[$x2+1];
  } else {
    $n_hi = _count_divisors_cumulative($x2);
  }
  $n_hi -= 1;

  ### rect at: "x=".($x2+1)." x_to_n=".($x_to_n[$x2+1]||'none')

  if ($self->{'proper'}) {
    $n_lo -= $x1-1;
    $n_hi -= $x2;
  }
  return ($n_lo + $self->{'n_start'},
          $n_hi + $self->{'n_start'});
}

# Return a total count of all the divisors of all the integers 1 to $x
# inclusive.
sub _count_divisors_cumulative {
  my ($x) = @_;
  my $total = 0;
  my $limit = _sqrtint($x);
  foreach my $i (1 .. $limit) {
    $total += int($x/$i);
  }
  return 2*$total - $limit*$limit;
}

1;
__END__