File Coverage

blib/lib/Math/NumSeq/ProthNumbers.pm

Criterion	Covered	Total	%
statement	86	90	95.5
branch	20	22	90.9
condition	4	6	66.6
subroutine	17	17	100.0
pod	6	6	100.0
total	133	141	94.3

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

package Math::NumSeq::ProthNumbers;
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;

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


# use constant name => Math::NumSeq::__('Proth Numbers');
use constant description => Math::NumSeq::__('Proth numbers k*2^n+1 for odd k and k < 2^n.');
use constant i_start => 1;
use constant values_min => 3;
use constant characteristic_increasing => 1;
use constant characteristic_integer => 1;

# cf A157892 - value of k
#    A157893 - value of n
#    A080076 - Proth primes
#    A134876 - how many Proth primes for given n
#
#    A002253 3*2^n+1 is prime
#    A002254 5*2^n+1
#    A032353 7*2^n+1
#    A002256 9*2^n+1
#
#    A086341 - 2^x+/-1 is sqrt of Proth squares
#
use constant oeis_anum => 'A080075';


# $uv_limit is the biggest power 4**k which fits in a UV
#
my $uv_limit = do {
  my $max = ~0;

  # limit <= floor(max/4) means 4*limit <= max is still good
  my $limit = 4;
  while ($limit <= ($max >> 2)) {
    $limit <<= 2;
  }

  ### ProthNumbers UV limit ...
  ### $limit
  ### limit: sprintf('%#X',$limit)

  $limit
};

sub rewind {
  my ($self) = @_;
  $self->{'value'} = 1;
  $self->{'inc'} = 2;
  $self->{'limit'} = 4;
  $self->{'i'} = $self->i_start;
}
sub next {
  my ($self) = @_;
  ### ProthNumbers next(): sprintf('value=%b inc=%b limit=%b', $self->{'value'}, $self->{'inc'}, $self->{'limit'})

  my $value = ($self->{'value'} += $self->{'inc'});
  if ($value >= $self->{'limit'}) {
    ### bigger 2^n now ...
    if ($self->{'limit'} == $uv_limit) {
      ### go to BigInt
      $self->{'value'} = Math::NumSeq::_to_bigint($self->{'value'});
      $self->{'inc'}   = Math::NumSeq::_to_bigint($self->{'inc'});
      $self->{'limit'} = Math::NumSeq::_to_bigint($self->{'limit'});
    }
    $self->{'inc'} *= 2;
    $self->{'limit'} *= 4;
  }
  return ($self->{'i'}++, $value);
}

# seek value at i-1, inc and limit ready to make next value
#
# i     i         value       inc     limit           next value
#  1     1             1       10         100               11
#  2    10            11       10         100              101
#  3    11           101      100       10000             1001
#  4   100          1001      100       10000             1101
#  5   101          1101      100       10000            10001
#  6   110         10001     1000     1000000            11001
#  7   111         11001     1000     1000000           100001
#  8  1000        100001     1000     1000000           101001
#  9  1001        101001     1000     1000000           110001
# 10  1010        110001     1000     1000000           111001
# 11  1011        111001     1000     1000000          1000001
# 12  1100       1000001    10000   100000000          1010001
# 13  1101       1010001    10000   100000000          1100001
# 14  1110       1100001    10000   100000000          1110001
# 15  1111       1110001    10000   100000000         10000001
# 16 10000      10000001    10000   100000000         10010001
#
sub seek_to_i {
  my ($self, $i) = @_;
  $self->{'i'} = $i;

  my ($pow) = _round_down_pow($i, 2);  # i = 1zxxx
  $i -= $pow;                          # without high bit, so i=zxxx

  # $zpow is the z bit position, or 1 if $i==1
  my $zpow = ($pow >= 2 ? $pow/2 : 0);

  if ($i >= $zpow) {
    # z=1 so return value=1xxx 0 0001, extra low 0-bit
    $pow *= 2;
  } else {
    # z=0 so return value=1xxx 0001, turn on 1-bit in i
    $i += $zpow;    # i=0xxx becomes i=1xxx
  }

  if ($pow*$pow >= $uv_limit) {
    $pow = Math::NumSeq::_to_bigint($pow);
  }
  $self->{'value'} = $pow*$i + 1;
  $self->{'inc'}   = $pow;
  $self->{'limit'} = $pow*$pow;
}

sub ith {
  my ($self, $i) = @_;
  ### ProthNumbers ith(): $i

  $i += 1;
  my ($pow) = _round_down_pow($i, 2);  # i = 1zxxx
  $i -= $pow;                          # without high bit, so i=zxxx

  my $zpow = $pow/2;
  if ($i >= $zpow) { # test z bit
    # z=1 so return value=1xxx 0 0001, extra low 0-bit
    $pow *= 2;
  } else {
    # z=1 so return value=1xxx 0001, turn on 1-bit in i
    $i += $zpow;    # i=0xxx becomes i=1xxx
  }
  return $i*$pow + 1;
}

sub pred {
  my ($self, $value) = @_;
  ### ProthNumbers pred(): $value
  {
    my $int = int($value);
    if ($value != $int) { return 0; }
    $value = $int;
  }
  unless ($value >= 3
          && ($value % 2)
          && ! _is_infinite($value)) {
    return 0;
  }

  # ENHANCE-ME: maybe BigInt as_bin() and string match position of lowest 1
  # ENHANCE-ME: maybe (v^(v-1))+1 for lowest 1, ask if v < m*m
  #
  my $pow = ($value*0) + 2; # inherit bignum 2
  for (;;) {
    ### at: "$value   $pow"
    $value = int($value/2);
    if ($value < $pow) {
      return 1;
    }
    if ($value % 2) {
      return ($value < $pow);
    }
    $pow *= 2;
  }
}

# use 41/29 ~= sqrt(2) in case $value is a Math::BigInt, since can't
# multiply BigInt*float
#
sub value_to_i_estimate {
  my ($self, $value) = @_;
  if ($value < 0) { return 0; }
  my $est = int(sqrt(int($value)) * 41 / 29);
}

1;
__END__