File Coverage

blib/lib/Math/NumSeq/Tetrahedral.pm

Criterion	Covered	Total	%
statement	57	59	96.6
branch	6	6	100.0
condition			n/a
subroutine	17	18	94.4
pod	5	5	100.0
total	85	88	96.5

# 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 A005894 centered tetrahedral numbers. (2*n+1)*(n^2+n+3)/3
#    A005906 truncated tetrahedral numbers. (n+1)*(23*n^2+19*n+6)/6
#    A015219 odd tetrahedrals (4n+1)(4n+2)(4n+3)/6
#    A015220 even tetrahedrals


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

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

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


# 6*value = i*(i+1)*(i+2)
#         = i^3 + 3*i^2 + 2*i
# so i^3 < 6T(i) < (i+1)^3
#
sub value_to_i_estimate {
  my ($self, $value) = @_;
  return _cbrt_floor(6*$value);
}
sub value_to_i_floor {
  my ($self, $value) = @_;
  ### value_to_i_floor(): "$value"

  $value *= 6;
  if ($value >= 0) {
    my $i = _cbrt_floor($value);
    if ($i*($i+1)*($i+2) <= $value) {
      return $i;
    } else {
      return $i-1;
    }
  } else {
    # secret undocumented negatives ...

    $value = abs($value);
    my $i = _cbrt_floor($value);

    ### $i
    ### prod: $i*($i+1)*($i+2)
    ### value*6: "$value"

    if ($i*($i+1)*($i+2) >= $value) {
      return -2-$i;
    } else {
      return -3-$i;
    }
  }
}

1;
__END__