File Coverage

blib/lib/Math/NumSeq/PolignacObstinate.pm

Criterion	Covered	Total	%
statement	69	70	98.5
branch	11	12	91.6
condition	7	9	77.7
subroutine	16	16	100.0
pod	3	3	100.0
total	106	110	96.3

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


# cf
# http://mathworld.wolfram.com/dePolignacsConjecture.html
#   consecutive primes difference is every even number infinitely many times

# 1,
# 127, 149, 251, 331,
# 337, 373, 509, 599,
# 701, 757, 809, 877,
# 905, 907, 959, 977,
# 997,

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

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

use Math::NumSeq::Primes; # for primes_list()
use Math::Prime::XS 'is_prime';

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


# use constant name => Math::NumSeq::__('Polignac Obstinate');
use constant description => Math::NumSeq::__('Odd integers N not representable as prime+2^k.');
use constant values_min => 1;
use constant characteristic_increasing => 1;
use constant characteristic_integer => 1;
use constant i_start => 1;

# cf A133122 - demanding k>0, so 1,3,127,... as 3 no longer representable
#    A065381 - primes not p+2^k k>=0, filtering from odds to primes
#
use constant oeis_anum => 'A006285'; # with k>=0 so 1,127,...


sub rewind {
  my ($self) = @_;
  $self->{'i'} = $self->i_start;
  $self->{'string'} = '';
  vec($self->{'string'},3/2,1) = 1;  # 2+2^0=3
  $self->{'done'} = -1;
  _resieve ($self, 20);
}

sub _resieve {
  my ($self, $hi) = @_;
  ### _resieve() ...

  $self->{'hi'} = $hi;
  my $sref = \$self->{'string'};
  vec($$sref,$hi,1) = 0;  # pre-extend
  my @primes = Math::NumSeq::Primes::_primes_list (3, $hi-1);
  for (my $power = 2; $power < $hi; $power *= 2) {
    foreach my $p (@primes) {
      if ((my $v = $p + $power) > $hi) {
        last;
      } else {
        vec($$sref,$v/2,1) = 1;
      }
    }
  }
}

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

  my $v = $self->{'done'};
  my $sref = \$self->{'string'};
  my $hi = $self->{'hi'};

  for (;;) {
    ### consider: "v=".($v+1)."  cf done=$self->{'done'}"
    if (($v+=2) > $hi) {
      _resieve ($self,
                $hi = ($self->{'hi'} *= 2));
    }
    unless (vec($$sref,$v/2,1)) {
      return ($self->{'i'}++,
              $self->{'done'} = $v);
    }
  }
}

sub pred {
  my ($self, $value) = @_;
  ### Obstinate pred(): $value

  unless ($value >= 0 && $value <= 0xFFFF_FFFF) {
    return undef;
  }
  if ($value != int($value)
      || $value < 127
      || ($value % 2) == 0) {
    return ($value == 1);
  }
  $value = "$value"; # numize Math::BigInt for speed

  # Maybe an is_any_prime(...)
  for (my $power = 2; $power < $value; $power *= 2) {
    if (is_prime($value - $power)) {
      return 0;
    }
  }
  return 1;
}

1;
__END__