File Coverage

blib/lib/Math/NumSeq/LuckyNumbers.pm

Criterion	Covered	Total	%
statement	65	65	100.0
branch	8	8	100.0
condition			n/a
subroutine	14	14	100.0
pod	2	2	100.0
total	89	89	100.0

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


# http://oeis.org/wiki/Lucky_numbers
#
# Gardiner, R. Lazarus, N. Metropolis and S. Ulam,"On certain sequences of
# integers defined by sieves", Mathematics Magazine 29:3 (1955),
# pp. 117-122.


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

use vars '$VERSION','@ISA';
$VERSION = 71;

use Math::NumSeq 7; # v.7 for _is_infinite()
@ISA = ('Math::NumSeq');
*_is_infinite = \&Math::NumSeq::_is_infinite;

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


# use constant name => Math::NumSeq::__('Lucky Numbers');
use constant description => Math::NumSeq::__('Sieved out multiples according to the sequence itself.');
use constant values_min => 1;
use constant i_start => 1;
use constant characteristic_increasing => 1;
use constant characteristic_integer => 1;

#------------------------------------------------------------------------------
# cf A145649 - 0,1 characteristic of Lucky numbers
#    A050505 - complement, the non-Lucky numbers
#
#    A007951 - ternary sieve, dropping 3rd, 6th, 9th, etc
#    1,2,_,4,5,_,7,8,_,10,11,_,12,13,_,14,15,_
#                              ^9th
#    1,2,4,5,7,8,10,11,14,16,17,19,20,22,23,25,28,29,31,32,34,35,37,38,41,
#
use constant oeis_anum => 'A000959';

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

# i=22
# mod: '8,12,14,20,24'
# int(21/24)=0 to i=21
# int(21/20)=1 to i=22
# int(22/14)=1 to i=23     two_pos=2 three_pos=1
# int(23/12)=1 to i=24
# int(24/8)=3 to i=27
# int(27/6)=4 to i=31
# int(31/2)=15 to i=46
# 2*46+1=93
#
# i+1 then twos advance for +1 more = +2
# threes advance
# i>3*mod[pos]  i-=2 pos++
# i-2>3*mod[pos+1]
#     3*2

sub rewind {
  my ($self) = @_;
  $self->{'i'}         = $self->i_start;
  $self->{'mod'}       = [ 8 ];
  $self->{'mod_pos'}   = 0;
  $self->{'one_pos'}   = 0;
  $self->{'two_pos'}   = 0;
  $self->{'three_pos'} = 0;
}

my @twelve = (1,3,
              7,9,
              13,15,
              21,
              25,27,
              31,33,
              37
             );

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

  my $ret_i = $self->{'i'}++;
  my $i = $ret_i - 1;
  my $mod = $self->{'mod'};
  if ($i >= $mod->[-1]) {
    ### extend mod ...
    my $mod_i = $#$mod + 4;
    my $mod_pos = $self->{'mod_pos'};
    if ($mod_i >= $mod->[$mod_pos]) {
      $self->{'mod_pos'} = ++$mod_pos;
    }
    push @$mod, _ith_by_mod($mod, $mod_i, $mod_pos) - 1;
  }

  # my $one_pos = $#$mod;
  my $two_pos = $self->{'two_pos'};

  ### mod: join(',',@{$self->{'mod'}})
  ### at: "i=$i one_pos=$#$mod two_pos=$self->{'two_pos'} diff=".($#$mod - $two_pos)

  $i += $#$mod - $two_pos;
  ### after ones: "i=$i cmp ".(2*$mod->[$two_pos])

  if ($i > 2*$mod->[$two_pos]) {
    ### advance two ...
    $self->{'two_pos'} = ++$two_pos;
    $i--;
  }
  ### after twos: "i=$i  two_pos=$two_pos  three_pos=$self->{'three_pos'}"

  my $three_pos = $self->{'three_pos'};
  $i += 2*($two_pos - $three_pos);
  ### increased i with twos: $i

  ### cmp: "i=$i three_pos=$three_pos  cmp ".(3*$mod->[$three_pos])
  if ($i > 3*$mod->[$three_pos]) {
    ### advance three ...
    $self->{'three_pos'} = ++$three_pos;
    $i -= 2;
  }

  return ($ret_i,
          _ith_by_mod($mod, $i, $three_pos));
}

sub _ith_by_mod {
  my ($mod, $i, $pos) = @_;
  ### _ith_by_mod(): "i=$i pos=$pos is mods=".join(',',@{$mod}[0..$pos-1])

  while (--$pos >= 0) {
    ### at: "pos=$pos i=$i"
    $i += int($i / $mod->[$pos]);
  }
  ### final: "i=$i  result=".(int($i/12)*21 + $twelve[$i%12])
  return int($i/12)*42 + $twelve[$i%12];
}

# i~=value/log(value)
#
use Math::NumSeq::Primes;
*value_to_i_estimate = \&Math::NumSeq::Primes::value_to_i_estimate;

1;
__END__