view release on metacpan or  search on metacpan

lib/Math/NumSeq/Tetrahedral.pm  view on Meta::CPAN

$VERSION = 75;

use Math::NumSeq;
use Math::NumSeq::Base::IterateIth;
@ISA = ('Math::NumSeq::Base::IterateIth',
        'Math::NumSeq');

use Math::NumSeq::Cubes;
*_cbrt_floor = \&Math::NumSeq::Cubes::_cbrt_floor;

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


# use constant name => Math::NumSeq::__('Tetrahedral Numbers');
use constant description => Math::NumSeq::__('The tetrahedral numbers 0, 1, 4, 10, 20, 35, 56, 84, 120, etc, i*(i+1)*(i+2)/6.');
use constant default_i_start => 0;
use constant values_min => 0;
use constant characteristic_increasing => 1;
use constant characteristic_integer => 1;
use constant oeis_anum => 'A000292'; # tetrahedrals

# could next() by increment
# 0
# 1    +1
# 4    +3  +2
# 10   +6  +3
# 20  +10  +4
# 35  +15  +5
# 56  +21  +6
# 84  +28  +7
# 120 +36  +8
#
# T(i) = i*(i+1)*(i+2)/6
#      = i*(i^2 + 3i + 2)/6
#      = (i^3 + 3i^2 + 2i)/6

sub rewind {
  my ($self) = @_;
  $self->{'i'} = $self->i_start;
}
sub _UNTESTED__seek_to_value {
  my ($self, $value) = @_;
  $self->seek_to_i($self->value_to_i_ceil($value));
}

sub ith {
  my ($self, $i) = @_;
  return $i*($i+1)*($i+2)/6;
}

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

  $value *= 6;
  my $i = _cbrt_floor($value);
  return ($i*($i+1)*($i+2) == $value);
}

# Cubic root formula
# 6*T(i) = i^3 + 3i^2 + 2i
# i^3 + 3i^2 + 2i - value = 0
# subst j=i+1, i=j-1
# (j-1)^3 + 3(j-1)^2 + 2(j-1) - value = 0
# j^3-3j^2+3j-1 + 3j^2-6j+3 + 2j-2 - value = 0
# j^3 + (-3+3)^2 + (3-6+2)j + (-1+3+-2) - value = 0
# j^3 - j - value = 0   p=-1 q=-v
#
# x^3+px+q=0
#  x=a-b
# a^3 - 3*b*a^2 + 3*b^2*a - b^3 + p(a-b) + q = 0
# a^3 - b^3 + 3ab(-a+b) + p(a-b) + q = 0
# a^3 - b^3 + (-3ab+p)(a-b) + q = 0
#  and 3ab=p so -3ab+p=0
# a^3 - b^3 + q = 0
#  mul (3a)^3
# 27a^6 - (3ab)^3 + 27q*a^3 = 0
# 27a^6 - p^3 + 27q*a^3 = 0
# 27a^6 + 27q*a^3 - p^3 = 0
# 27(a^3)^2 + 27q*(a^3) - p^3 = 0    quadratic in a^3
# a^3 = (-27q + sqrt((27q)^2 + 4*27*p^3)) / 2*27
#     = (-q + sqrt(q^2 + 4*p^3/27)) / 2
# 3ab=p b=p/3a
# b^3 = (-27q + sqrt((27q)^2 - 4*27*p^3)) / 2*27
#
# v=56=6*7*8/6
# p=-1 q=-v
# a^3 = (v + sqrt(v^2 - 4/27)) / 2
#     = (56 + sqrt(56^2 - 4/27))/2
# a = 3.825847303806096100878703127
# b = -1/3a
# 
# j^3 - j - 6*value = 0
# a^3 - 3*b*a^2 + 3*b^2*a - b^3 + -(a-b) - 6v = 0
# a^3 - b^3 + 3ab(-a + b) + -(a-b) - 6v = 0
# a^3 - b^3 - 3ab(a-b) + -(a-b) - 6v = 0
# a^3 - b^3 + (-1-3ab)*(a-b) - 6v = 0
# -1-3ab=0 3ab=-1
# a^3 - b^3 - 6v = 0
# 27a^6 - (3ab)^3 - 27*6v*a^3 = 0
# 27a^6 - (-1)^3 - 27*6v*a^3 = 0
# 27(a^3)^2 - 27*6v*(a^3) - (-1)^3 = 0
# 27(a^3)^2 - 27*6v*(a^3) + 1 = 0
# (a^3)^2 - 6v*(a^3) + 1/27 = 0
# a^3 = (6v + sqrt((6v)^2 - 4/27))/2
#     = (6*56+sqrt((6*56)^2 - 4/27))/2
# 3ab=-1
# b=-1/3a
#
# 
# 6*T(i)   = i^3 + 3i^2 + 2i
# estimate i=cbrt(6*value)
# (i+1)^3 = i^3 + 3i^2 + 3i + 1 is bigger than T(i)
#
# v just below a cube so
# v=x^3-1
# then cbrt gives x
# T(x+1) = (x+1)^3 + 3*(x+1)^2 + 2*(x+1)
#        = x^3 + 6*x^2 + 11*x + 6