File Coverage

blib/lib/Math/NumSeq/RepdigitRadix.pm

Criterion	Covered	Total	%
statement	105	108	97.2
branch	35	40	87.5
condition	5	8	62.5
subroutine	18	18	100.0
pod	6	6	100.0
total	169	180	93.8

# Copyright 2010, 2011, 2012, 2013, 2014 Kevin Ryde

# This file is part of Math-NumSeq.
#
# Math-NumSeq 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-NumSeq 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-NumSeq.  If not, see .


# ENHANCE-ME: ith() might be a touch faster than next() now, perhaps
# something sieve/flag in next()


package Math::NumSeq::RepdigitRadix;
use 5.004;
use strict;

use vars '$VERSION', '@ISA';
$VERSION = 71;
use Math::NumSeq;
@ISA = ('Math::NumSeq');
*_is_infinite = \&Math::NumSeq::_is_infinite;

use Math::NumSeq::NumAronson;
*_round_down_pow = \&Math::NumSeq::NumAronson::_round_down_pow;

use Math::Factor::XS 0.40 'factors'; # version 0.40 for factors() on BigInt

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


# use constant name => Math::NumSeq::__('Repdigit Radix');
use constant description => Math::NumSeq::__('First radix in which i is a repdigit (at most base=i-1 since otherwise "11" would always give i).');
use constant characteristic_smaller => 1;
use constant characteristic_increasing => 0;
use constant characteristic_integer => 1;
use constant characteristic_value_is_radix => 1;
sub values_min {
  my ($self) = @_;
  return ($self->i_start >= 3 ? 2 : 0);
}
sub i_start {
  my ($self) = @_;
  return $self->{'i_start'} || 0;
}

# smallest base in which n is a repdigit, starting n=3
sub oeis_anum { 'A059711' }
# OEIS-Catalogue: A059711 i_start=3


# d * (b^3 + b^2 + b + 1) = i
# b^3 + b^2 + b + 1 = i/d
# (b+1)^3 = b^3 + 3b^2 + 3b + 1
#
# (b-1) * (b^3 + b^2 + b + 1) = b^4 - 1
#
# 8888 base 9 = 6560
# 1111 base 10

# b^2 + b + 1 = k
# (b+0.5)^2 + .75 = k
# (b+0.5)^2 = (k-0.75)
# b = sqrt(k-0.75)-0.5;

sub rewind {
  my ($self) = @_;
  $self->{'i'} = $self->i_start;
  $self->{'ones'}   = [ undef, undef, 7 ];
  $self->{'digits'} = [ undef, undef, 1 ];
  ### rewind to: $self
}

# (r+1)^2 + (r+1) + 1
#   = r^2 + 2r + 1 + r +1 + 1
#   = r^2 + 3r + 3
#   = (r + 3)*r + 3

#            0  1  2
my @small = (2, 0, 0);

sub next {
  my ($self) = @_;
  ### RepdigitRadix next(): $self->{'i'}

  my $i = $self->{'i'}++;
  my $ones = $self->{'ones'};
  my $digits = $self->{'digits'};

  if ($i < 3) {
    return ($i, $small[$i]);
  }

  for (my $radix = 2; ; $radix++) {
    ### $radix
    ### ones: $ones->[$radix]
    ### digit: $digits->[$radix]

    my $one;
    if ($radix > $#$ones) {
      ### maybe extend array: $radix
      $one = $radix + 1;  # or three digits ... ($radix + 1) *
      unless ($one <= $i) {
        ### not repdigit of 3 digits in any radix, take as 2 digits ...
        return ($i, $i-1);
      }
      $ones->[$radix] = $one;
      $digits->[$radix] = 1;

    } else {
      $one = $ones->[$radix];
    }

    my $repdigit = $one * $digits->[$radix];
    while ($repdigit < $i) {
      my $digit = ++$digits->[$radix];
      if ($digit >= $radix) {
        $digit = $digits->[$radix] = 1;
        $one = $ones->[$radix] = ($one * $radix + 1);
      }
      $repdigit = $one * $digit;
    }
    ### consider repdigit: $repdigit
    if ($repdigit == $i) {
      ### found radix: $radix
      return ($i, $radix);
    }
  }
}

# d=r-1
# v = d*(r^(len-1)+...+1)
#   = r^len-1
# dlimit = nthroot(v+1, len)
#
# q=v/d
# q=r^(len-1)+...+1
# q > r^(len-1)    # when len>2
# r < nthroot(q, len-1)

sub ith {
  my ($self, $i) = @_;
  ### RepdigitRadix ith(): $i
  if ($i < 0) {
    $i = abs($i);
  }
  if ($i < 3) {
    return $small[$i];
  }
  if ($i > 0xFFFF_FFFF) {
    return undef;
  }
  if (_is_infinite($i)) {
    return $i; # nan
  }

  my @factors = reverse (1, factors($i));
  ### @factors

  my ($pow, $len) = _round_down_pow ($i, 2);
  $len++;
  ### initial len: $len

  my $r_found;
  for ( ; $len >= 3; $len--) {
    my $d_limit = (defined $r_found ? $r_found-1 : _nth_root_floor($i+1,$len));
    ### $len
    ### $d_limit

    foreach my $d (grep {$_<=$d_limit} @factors) {  # descending order
      ### try d: $d

      # if ($d > $d_limit) {
      #   ### stop for d > d_limit ...
      #   last;
      # }

      my $q = $i / $d;
      my $r = _nth_root_floor($q,$len-1);
      ### $q
      ### $r
      ### ones: ($r**$len - 1) / ($r-1)

      if (defined $r_found && $r >= $r_found) {
        ### stop at r >= r_found ...
        last;
      }

      if ($r <= $d) {
        ### r smaller than d ...
        # since d>=1 this also excludes r<2
        next;
      }

      if ($q == ($r**$len - 1) / ($r-1)) {
        $r_found = $r;
      }
    }
  }

  my $d_limit = (defined $r_found ? $r_found-1 : int(sqrt($i+1)));
  foreach my $d (grep {$_<=$d_limit} @factors) {  # descending order
    ### try d: $d

    # v = d*(r+1)
    # v/d = r+1
    # r = v/d - 1
    #
    my $r = $i/$d - 1;
    ### $r

    if (defined $r_found && $r >= $r_found) {
      ### stop at r >= r_found ...
      last;
    }

    if ($r <= $d) {
      ### r smaller than d ...
      # since d>=1 this also excludes r<2
      next;
    }

    $r_found = $r;
  }

  return (defined $r_found ? $r_found : $i-1);


  # for (my $radix = 2; ; $radix++) {
  #   ### $radix
  #
  #   my $one = $radix + 1;  # ... or 3 digits 111 ($radix + 1) *
  #   unless ($one <= $i) {
  #     ### stop at ones too big not a 3-digit repdigit: $one
  #     return $i-1;
  #   }
  #   ### $one
  #
  #   do {
  #     if ($one == $i) {
  #       return $radix;
  #     }
  #     foreach my $digit (2 .. $radix-1) {
  #       ### $digit
  #       if ((my $repdigit = $digit * $one) <= $i) {
  #         if ($repdigit == $i) {
  #           return $radix;
  #         }
  #       }
  #     }
  #   } while (($one = $one * $radix + 1) <= $i);
  # }
}

# value = root^power
# log(value) = power*log(root)
# log(root) = log(value)/power
# root = exp(log(value)/power)
#
sub _nth_root_floor {
  my ($value, $power) = @_;
  my $root = int (exp (log($value)/$power));
  if ($root**$power > $value) {
    return $root-1;
  }
  if (($root+1)**$power < $value) {
    return $root+1;
  }
  return $root;
}

# R^2+R+1
# R=65 "111"=4291
#
# Does every radix occur?  Is it certain that at least one repdigit in base
# R is not a repdigit in anything smaller?
#
# sub pred {
#   my ($self, $value) = @_;
#   return ($value == int($value)
#           && ($value == 0 || $value >= 2));
# }

1;
__END__