File Coverage

blib/lib/Math/NumSeq/ReRound.pm

Criterion	Covered	Total	%
statement	71	75	94.6
branch	12	16	75.0
condition	3	3	100.0
subroutine	18	18	100.0
pod	4	4	100.0
total	108	116	93.1

# 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::ReRound;
use 5.004;
use strict;
use POSIX 'ceil';
use List::Util 'max';

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

use Math::NumSeq;
use Math::NumSeq::Base::IterateIth;
@ISA = ('Math::NumSeq::Base::IterateIth',
        'Math::NumSeq');
*_is_infinite = \&Math::NumSeq::_is_infinite;

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

# use constant name => Math::NumSeq::__('Repeated Rounding');
use constant description => Math::NumSeq::__('Repeated rounding up to multiple of N, N-1, ..., 1.');
use constant values_min => 1; # at i=1
use constant i_start => 1;
use constant characteristic_increasing => 1;
use constant characteristic_integer => 1;

use constant parameter_info_array =>
  [
   { name    => 'extra_multiples',
     display => Math::NumSeq::__('Extra Multiples'),
     type    => 'integer',
     default => '0',
     minimum => 0,
     width   => 2,
      description => Math::NumSeq::__('Extra multiples to round up at each stage.'),
   },
  ];

#------------------------------------------------------------------------------
# cf A007952 sieve+1
#    A099361 sieve by primes
#    A099204 sieve by primes
#    A099207 sieve by primes
#    A099243 sieve by primes
#    A002960 square sieve
#
#    A056533 sieve out even every 2nd,4th,6th,etc
#    A039672 sieve fibonacci style i+j prev terms
#    A056530 Flavius Josephus after 2nd round
#    A056531 Flavius Josephus after 4th round
#    A119446 cf A100461
#
#    A113749 k multiples
#
my @oeis_anum
  = (
     # OEIS-Catalogue array begin
     'A002491',   #                     # Mancala stones
     'A000960',   # extra_multiples=1   # Flavius Josephus
     'A112557',   # extra_multiples=2   # more Mancala stones ...
     'A112558',   # extra_multiples=3   # more Mancala stones ...
     'A113742',   # extra_multiples=4   # more Mancala stones ...
     'A113743',   # extra_multiples=5   # more Mancala stones ...
     'A113744',   # extra_multiples=6   # more Mancala stones ...
     'A113745',   # extra_multiples=7   # more Mancala stones ...
     'A113746',   # extra_multiples=8   # more Mancala stones ...
     'A113747',   # extra_multiples=9   # more Mancala stones ...
     'A113748',   # extra_multiples=10  # more Mancala stones ...
     # OEIS-Catalogue array end
    );
sub oeis_anum {
  my ($self) = @_;
  return $oeis_anum[$self->{'extra_multiples'}];
}

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


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

  if (_is_infinite($i)) {
    return $i;  # don't loop forever if $i is +infinity
  }

  my $extra_multiples = $self->{'extra_multiples'};
  ### $extra_multiples

  for (my $m = $i-1; $m >= 1; $m--) {
    ### add: (-$i % $m) + $extra_multiples*$m

    $i += (-$i % $m) + $extra_multiples*$m;
  }
  return $i;
}

# 1,3,7,13,19
# 2->2+1=3
# 3->4+2=6->6+1=7
# 4->6+m3=9->10+m2=12->12+m1=13
#
# next = prev + (-prev mod m) + k*m
# next-k*m = prev + (-prev mod m)

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

  my $extra_multiples = $self->{'extra_multiples'};

  if (_is_infinite($value)) {
    return undef;
  }
  if ($value <= 1 || $value != int($value)) {
    return ($value == 1);
  }

  # special case m=1 stepping down to an even number
  if (($value -= $extra_multiples) % 2) {
    return 0;
  }

  my $m = 2;
  while ($value > $m) {
    ### at: "value=$value  m=$m"

    if (($value -= $extra_multiples*$m) <= 0) {
      ### no, negative: $value
      return 0;
    }
    ### subtract to: "value=$value"

    ### rem: "modulus=".($m+1)." rem ".($value%($m+1))
    my $rem;
    if (($rem = ($value % ($m+1))) == $m) {
      ### no, remainder: "rem=$rem  modulus=".($m+1)
      return 0;
    }

    $value -= $rem;
    $m++;
  }

  ### final ...
  ### $value
  ### $m

  return ($value == $m);
}

# # v1.02 for leading underscore
use constant 1.02 _PI => 4*atan2(1,1); # similar to Math::Complex pi()

# value = m*1+m*2+m*3+...+m*i
# value = m*i*(i+1)/2
# i*i + i - 2v/m = 0
# i = (-1 + sqrt(1 + 4*2v/m)) / 2
#   = -1/2 + sqrt(1 + 4*2v/m)/2
#   = -1/2 + sqrt(4 + 2v/m)
# i -> sqrt(2v/m)
#
# as extra_multiples dominates the estimate changes from
#     extra_multiples==0         est = sqrt(pi*value)
# up to
#     extra_multiples==m large   est = sqrt(2*value/m)
#
# What formula that progressively morphs pi to 2/m ?
#
sub value_to_i_estimate {
  my ($self, $value) = @_;
  if ($value < 0) { return 0; }

  if ($self->{'extra_multiples'} == 0) {
    # use sqrt(pi)~=296/167 to preserve Math::BigInt
    return int( (sqrt(int($value)) * 296) / 167);
  } else {
    return int(sqrt(int (2 * $value / $self->{'extra_multiples'})));

    # return int(sqrt(int ($value * _PI/($self->{'extra_multiples'}+1)**2))
    #            + ($self->{'extra_multiples'} > 0
    #               ? -0.5 + sqrt(4 + 2*$value/$self->{'extra_multiples'})
    #               : 0));
  }



  # # large extra_multiples
  # return int(sqrt(int (2 * $value
  #                      / ($self->{'extra_multiples'}+1))));
  #
  # return int(sqrt(int (($value * 355)
  #                      / (113 * ($self->{'extra_multiples'}+1)))));
  #
  #
  # # extra_multiples==14
  # return int(sqrt(int (_PI * $value)) * 429/2048); #        429=13*11*3
  #
  # # extra_multiples==12
  # return int(sqrt(int (_PI * $value)) * 462/2048); # 231/1024) # 11*7*3=231
  #
  # # extra_multiples==10
  # return int(sqrt(int (_PI * $value))) * 126/512; # 63/256;  # 7*3*2=126
  #
  # # extra_multiples==8
  # return int(sqrt(int (_PI * $value))) * 35/128;             # 7*5=35
  #
  # # extra_multiples==6
  # return int(sqrt(int (_PI * $value))) * 10/32; # 5/16       # 5*2=10
  #
  # # extra_multiples==4
  # return int(sqrt(int (_PI * $value))) * 3/8;
  #
  # # extra_multiples==2
  # return int(sqrt(int (_PI * $value))) * 1/2;
  #
  # # extra_multiples==0
  # return int(sqrt(int (_PI * $value))) * 1/1;
  #
  #  2 1
  #  4 11
  #  6 1010
  #  8 100011
  # 10 1111110
  # 12 111001110
  # 14 110101101

  #
  #
  # # extra_multiples==11
  # return int(sqrt(int (1/_PI * $value)) * 1024/1386); # 512/693);
  #
  # # extra_multiples==9
  # return int(sqrt(int (1/_PI * $value)) * 256/315);
  #
  # # extra_multiples==7
  # return int(sqrt(int (1/_PI * $value)) * 64/70); # 32/35;
  #
  # # extra_multiples==5
  # return int(sqrt(int (1/_PI * $value))) * 16/15;
  #
  # # extra_multiples==3
  # return int(sqrt(int (1/_PI * $value))) * 4/3;
  #
  # # extra_multiples==1
  # return int(sqrt(int (1/_PI * $value ))) * 2/1;

}

1;
__END__