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


line stmt bran cond sub pod time code
1             # Copyright 2010, 2011, 2012, 2013, 2014 Kevin Ryde
2              
3             # This file is part of Math-NumSeq.
4             #
5             # Math-NumSeq is free software; you can redistribute it and/or modify
6             # it under the terms of the GNU General Public License as published by the
7             # Free Software Foundation; either version 3, or (at your option) any later
8             # version.
9             #
10             # Math-NumSeq is distributed in the hope that it will be useful, but
11             # WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
12             # or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
13             # for more details.
14             #
15             # You should have received a copy of the GNU General Public License along
16             # with Math-NumSeq. If not, see .
17              
18              
19              
20             # cf A005894 centered tetrahedral numbers. (2*n+1)*(n^2+n+3)/3
21             # A005906 truncated tetrahedral numbers. (n+1)*(23*n^2+19*n+6)/6
22             # A015219 odd tetrahedrals (4n+1)(4n+2)(4n+3)/6
23             # A015220 even tetrahedrals
24              
25              
26             package Math::NumSeq::Tetrahedral;
27 2     2   7175 use 5.004;
  2         5  
28 2     2   7 use strict;
  2         2  
  2         46  
29              
30 2     2   6 use vars '$VERSION','@ISA';
  2         3  
  2         111  
31             $VERSION = 72;
32              
33 2     2   364 use Math::NumSeq;
  2         2  
  2         37  
34 2     2   324 use Math::NumSeq::Base::IterateIth;
  2         3  
  2         53  
35             @ISA = ('Math::NumSeq::Base::IterateIth',
36             'Math::NumSeq');
37              
38 2     2   318 use Math::NumSeq::Cubes;
  2         3  
  2         73  
39             *_cbrt_floor = \&Math::NumSeq::Cubes::_cbrt_floor;
40              
41             # uncomment this to run the ### lines
42             #use Smart::Comments;
43              
44              
45             # use constant name => Math::NumSeq::__('Tetrahedral Numbers');
46 2     2   8 use constant description => Math::NumSeq::__('The tetrahedral numbers 0, 1, 4, 10, 20, 35, 56, 84, 120, etc, i*(i+1)*(i+2)/6.');
  2         2  
  2         5  
47 2     2   6 use constant default_i_start => 0;
  2         2  
  2         67  
48 2     2   11 use constant values_min => 0;
  2         2  
  2         60  
49 2     2   6 use constant characteristic_increasing => 1;
  2         2  
  2         68  
50 2     2   6 use constant characteristic_integer => 1;
  2         2  
  2         62  
51 2     2   6 use constant oeis_anum => 'A000292'; # tetrahedrals
  2         2  
  2         532  
52              
53             # could next() by increment
54             # 0
55             # 1 +1
56             # 4 +3 +2
57             # 10 +6 +3
58             # 20 +10 +4
59             # 35 +15 +5
60             # 56 +21 +6
61             # 84 +28 +7
62             # 120 +36 +8
63             #
64             # T(i) = i*(i+1)*(i+2)/6
65             # = i*(i^2 + 3i + 2)/6
66             # = (i^3 + 3i^2 + 2i)/6
67              
68             sub rewind {
69 5     5 1 351 my ($self) = @_;
70 5         15 $self->{'i'} = $self->i_start;
71             }
72             sub _UNTESTED__seek_to_value {
73 0     0   0 my ($self, $value) = @_;
74 0         0 $self->seek_to_i($self->value_to_i_ceil($value));
75             }
76              
77             sub ith {
78 93     93 1 68 my ($self, $i) = @_;
79 93         160 return $i*($i+1)*($i+2)/6;
80             }
81              
82             sub pred {
83 250     250 1 557 my ($self, $value) = @_;
84             ### Tetrahedral pred(): $value
85              
86 250         145 $value *= 6;
87 250         262 my $i = _cbrt_floor($value);
88 250         286 return ($i*($i+1)*($i+2) == $value);
89             }
90              
91             # Cubic root formula
92             # 6*T(i) = i^3 + 3i^2 + 2i
93             # i^3 + 3i^2 + 2i - value = 0
94             # subst j=i+1, i=j-1
95             # (j-1)^3 + 3(j-1)^2 + 2(j-1) - value = 0
96             # j^3-3j^2+3j-1 + 3j^2-6j+3 + 2j-2 - value = 0
97             # j^3 + (-3+3)^2 + (3-6+2)j + (-1+3+-2) - value = 0
98             # j^3 - j - value = 0 p=-1 q=-v
99             #
100             # x^3+px+q=0
101             # x=a-b
102             # a^3 - 3*b*a^2 + 3*b^2*a - b^3 + p(a-b) + q = 0
103             # a^3 - b^3 + 3ab(-a+b) + p(a-b) + q = 0
104             # a^3 - b^3 + (-3ab+p)(a-b) + q = 0
105             # and 3ab=p so -3ab+p=0
106             # a^3 - b^3 + q = 0
107             # mul (3a)^3
108             # 27a^6 - (3ab)^3 + 27q*a^3 = 0
109             # 27a^6 - p^3 + 27q*a^3 = 0
110             # 27a^6 + 27q*a^3 - p^3 = 0
111             # 27(a^3)^2 + 27q*(a^3) - p^3 = 0 quadratic in a^3
112             # a^3 = (-27q + sqrt((27q)^2 + 4*27*p^3)) / 2*27
113             # = (-q + sqrt(q^2 + 4*p^3/27)) / 2
114             # 3ab=p b=p/3a
115             # b^3 = (-27q + sqrt((27q)^2 - 4*27*p^3)) / 2*27
116             #
117             # v=56=6*7*8/6
118             # p=-1 q=-v
119             # a^3 = (v + sqrt(v^2 - 4/27)) / 2
120             # = (56 + sqrt(56^2 - 4/27))/2
121             # a = 3.825847303806096100878703127
122             # b = -1/3a
123             #
124             # j^3 - j - 6*value = 0
125             # a^3 - 3*b*a^2 + 3*b^2*a - b^3 + -(a-b) - 6v = 0
126             # a^3 - b^3 + 3ab(-a + b) + -(a-b) - 6v = 0
127             # a^3 - b^3 - 3ab(a-b) + -(a-b) - 6v = 0
128             # a^3 - b^3 + (-1-3ab)*(a-b) - 6v = 0
129             # -1-3ab=0 3ab=-1
130             # a^3 - b^3 - 6v = 0
131             # 27a^6 - (3ab)^3 - 27*6v*a^3 = 0
132             # 27a^6 - (-1)^3 - 27*6v*a^3 = 0
133             # 27(a^3)^2 - 27*6v*(a^3) - (-1)^3 = 0
134             # 27(a^3)^2 - 27*6v*(a^3) + 1 = 0
135             # (a^3)^2 - 6v*(a^3) + 1/27 = 0
136             # a^3 = (6v + sqrt((6v)^2 - 4/27))/2
137             # = (6*56+sqrt((6*56)^2 - 4/27))/2
138             # 3ab=-1
139             # b=-1/3a
140             #
141             #
142             # 6*T(i) = i^3 + 3i^2 + 2i
143             # estimate i=cbrt(6*value)
144             # (i+1)^3 = i^3 + 3i^2 + 3i + 1 is bigger than T(i)
145             #
146             # v just below a cube so
147             # v=x^3-1
148             # then cbrt gives x
149             # T(x+1) = (x+1)^3 + 3*(x+1)^2 + 2*(x+1)
150             # = x^3 + 6*x^2 + 11*x + 6
151              
152              
153             # 6*value = i*(i+1)*(i+2)
154             # = i^3 + 3*i^2 + 2*i
155             # so i^3 < 6T(i) < (i+1)^3
156             #
157             sub value_to_i_estimate {
158 18     18 1 332 my ($self, $value) = @_;
159 18         26 return _cbrt_floor(6*$value);
160             }
161             sub value_to_i_floor {
162 80     80 1 161 my ($self, $value) = @_;
163             ### value_to_i_floor(): "$value"
164              
165 80         80 $value *= 6;
166 80 100       802 if ($value >= 0) {
167 59         859 my $i = _cbrt_floor($value);
168 59 100       1070 if ($i*($i+1)*($i+2) <= $value) {
169 32         2184 return $i;
170             } else {
171 27         44 return $i-1;
172             }
173             } else {
174             # secret undocumented negatives ...
175              
176 21         17 $value = abs($value);
177 21         38 my $i = _cbrt_floor($value);
178              
179             ### $i
180             ### prod: $i*($i+1)*($i+2)
181             ### value*6: "$value"
182              
183 21 100       37 if ($i*($i+1)*($i+2) >= $value) {
184 17         45 return -2-$i;
185             } else {
186 4         14 return -3-$i;
187             }
188             }
189             }
190              
191             1;
192             __END__