File Coverage

blib/lib/Math/NumSeq/Factorials.pm
Criterion Covered Total %
statement 115 132 87.1
branch 28 38 73.6
condition 8 11 72.7
subroutine 22 24 91.6
pod 8 8 100.0
total 181 213 84.9


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             # http://www.luschny.de/math/factorial/approx/SimpleCases.html
20             #
21              
22              
23             package Math::NumSeq::Factorials;
24 2     2   481470 use 5.004;
  2         10  
  2         204  
25 2     2   16 use strict;
  2         3  
  2         78  
26              
27 2     2   10 use vars '$VERSION','@ISA';
  2         3  
  2         135  
28             $VERSION = 71;
29              
30 2     2   682 use Math::NumSeq;
  2         5  
  2         83  
31             @ISA = ('Math::NumSeq');
32             *_is_infinite = \&Math::NumSeq::_is_infinite;
33              
34 2     2   599 use Math::NumSeq::Fibonacci;
  2         5  
  2         114  
35             *_blog2_estimate = \&Math::NumSeq::Fibonacci::_blog2_estimate;
36              
37             # uncomment this to run the ### lines
38             #use Smart::Comments;
39              
40              
41             # use constant name => Math::NumSeq::__('Factorials');
42 2     2   14 use constant description => Math::NumSeq::__('The factorials 1, 2, 6, 24, 120, etc, 1*2*...*N.');
  2         5  
  2         8  
43 2     2   10 use constant values_min => 1;
  2         4  
  2         74  
44 2     2   12 use constant i_start => 0;
  2         3  
  2         79  
45 2     2   8 use constant characteristic_increasing => 1;
  2         5  
  2         243  
46 2     2   21 use constant characteristic_integer => 1;
  2         5  
  2         149  
47              
48             #------------------------------------------------------------------------------
49             # cf A006882 double a(n)=n*a(n-2), n*(n-2)*(n-4)*...*3*1 or *4*2
50             # A001147 double factorial 1*3*5*...*(2n-1) odd numbers, bisection
51             # A000165 double factorial 2*4*6*...*2n even numbers, bisection
52             #
53             # A007661 triple a(n)=n*a(n-3)
54             # A007662 quadruple a(n)=n*a(n-4)
55             # A047053 quad 4^n*n! quad on multiples of 4
56             # A007696 quad n=4k+1 products
57             # A001813 quad (2*n)!/n!
58             # A008545 quad n=4k+1 products
59             # A080500 squares n*(n-1)*(n-4)*(n-9)*(n-16)*(n-25)*...*1
60             #
61             # A001013 Jordan-Polya products of factorials
62             #
63             # A008906 n! num digits excl trailing zeros
64             # A027868 n! num trailing zeros, is power of 5
65             # A000966 n! never ends these 0s
66             #
67             # A008904 n! low non-zero
68             # A136690 base 3
69             # A136691 base 4
70             # A136692 base 5
71             # A136693 base 6
72             # A136694 base 7
73             # A136695 base 8
74             # A136696 base 9
75             # A136697 base 11
76             # A136698 base 12
77             # A136699 base 13
78             # A136700 base 14
79             # A136701 base 15
80             # A136702 base 16
81             #
82             # A008905 n! leading digit
83             # A136754 base 3
84             # A136755 base 4
85             # A136756 base 5
86             # A136757 base 6
87             # A136758 base 7
88             # A136759 base 8
89             # A136760 base 9
90             # A136761 base 11
91             # A136762 base 12
92             # A136763 base 13
93             # A136764 base 14
94             # A136765 base 15
95             # A136766 base 16
96              
97 2     2   9 use constant oeis_anum => 'A000142'; # factorials 1,1,2,6,24, including 0!==1
  2         5  
  2         107  
98              
99             #------------------------------------------------------------------------------
100              
101              
102 2     2   10 use constant 1.02; # for leading underscore
  2         43  
  2         195  
103 2         3 use constant _UV_I_LIMIT => do {
104 2         4 my $u = ~0 >> 1;
105 2         3 my $limit = 1;
106 2         6 my $i = 2;
107 2         112 for (; $i++; ) {
108 38 100       125 if ($u < $i) {
109             ### _UV_LIMIT stop before: "i=$i"
110 2         6 last;
111             }
112 36         38 $u -= ($u % $i);
113 36         33 $u /= $i;
114 36         54 $limit *= $i;
115             }
116             ### $limit
117             ### $i
118 2         230 $i
119 2     2   9 };
  2         5  
120             ### _UV_I_LIMIT: _UV_I_LIMIT()
121              
122 2         4 use constant _NV_LIMIT => do {
123 2         5 my $f = 1.0;
124 2         2 my $max;
125 2         4 for (;;) {
126 104         98 $max = $f;
127 104         100 my $l = 2.0*$f;
128 104         100 my $h = 2.0*$f+2.0;
129 104         92 $f = 2.0*$f + 1.0;
130 104         134 $f = sprintf '%.0f', $f;
131 104 100 66     395 last unless ($f < $h && $f > $l);
132             }
133             ### uv : ~0
134             ### 53 : 1<<53
135             ### $max
136 2         1783 $max
137 2     2   12 };
  2         9  
138              
139              
140              
141             sub rewind {
142 7     7 1 230 my ($self) = @_;
143             ### Factorials rewind()
144 7         38 $self->{'i'} = $self->i_start;
145 7         22 $self->{'f'} = 1;
146             }
147             sub seek_to_i {
148 20     20 1 642 my ($self, $i) = @_;
149 20         29 $self->{'i'} = $i;
150 20         48 $self->{'f'} = $self->ith($i-1);
151             }
152             sub _UNTESTED__seek_to_value {
153 0     0   0 my ($self, $value) = @_;
154 0         0 my $i = $self->{'i'} = $self->value_to_i_ceil($value);
155 0         0 $self->{'f'} = $self->ith($i);
156             }
157             sub next {
158 55     55 1 893 my ($self) = @_;
159             ### Factorials next() ...
160              
161 55         85 my $i = $self->{'i'}++;
162 55 50       94 if ($i == _UV_I_LIMIT) {
163 0         0 $self->{'f'} = Math::NumSeq::_to_bigint($self->{'f'});
164             }
165 55   100     211 return ($i, $self->{'f'} *= ($i||1));
166             }
167              
168             sub ith {
169 57     57 1 1045 my ($self, $i) = @_;
170             ### Factorials ith() ...
171              
172 57 50       176 if (_is_infinite($i)) {
173 0         0 return $i;
174             }
175              
176 57 100 66     244 if (! ref $i && $i >= _UV_I_LIMIT) {
177             # Plain index $i automatically use Math::BigInt when UV limit reached.
178             # Maybe should check if BigInt new enough to have bfac(), circa vers 1.60
179 3         1932 return Math::NumSeq::_bigint()->bfac($i);
180             }
181              
182 54         70 my $value = ($i*0) + 1; # inherit bignum 1
183 54         111 while ($i >= 2) {
184 168         167 $value *= $i;
185 168         296 $i -= 1;
186             }
187 54         144 return $value;
188             }
189              
190             sub pred {
191 1463     1463 1 7312 my ($self, $value) = @_;
192 1463         2521 return defined($self->value_to_i($value));
193             }
194             sub value_to_i {
195 1495     1495 1 5489 my ($self, $value) = @_;
196              
197             # NV inf or nan gets $value%$i=nan and nan==0 is false,
198             # but Math::BigInt binf()%$i=0 so would go into infinite loop
199             # hence explicit check against _is_infinite()
200             #
201 1495 50       3434 if (_is_infinite($value)) {
202 0         0 return undef;
203             }
204              
205 1495 100       5246 if ($value == 1) {
206 8         225 return 0;
207             }
208              
209 1487         2118 my $i = 1;
210 1487         1641 for (;;) {
211 2589 100       6461 if ($value <= 1) {
212 31 100       454 return ($value == 1 ? $i : undef);
213             }
214 2558         3913 $i++;
215 2558 100       4375 if (($value % $i) == 0) {
216 1102         4763 $value /= $i;
217             } else {
218 1456         4142 return undef;
219             }
220             }
221 0         0 return $i;
222             }
223              
224             sub value_to_i_floor {
225 52     52 1 169 my ($self, $value) = @_;
226 52 50       105 if (_is_infinite($value)) {
227 0         0 return $value;
228             }
229 52 100       1815 if ($value < 2) {
230 14         180 return $self->i_start;
231             }
232              
233             # "/" operator converts 64-bit UV to an NV and so loses bits, making the
234             # result come out 1 too small sometimes. Experimental switch to BigInt to
235             # keep precision. ENHANCE-ME: Maybe better _divrem().
236             #
237 38 50 66     517 if (! ref $value && $value > _NV_LIMIT) {
238 0         0 $value = Math::NumSeq::_to_bigint($value);
239             }
240              
241 38         46 my $i = 2;
242 38         42 for (;; $i++) {
243             ### $value
244             ### $i
245              
246 137 100       1982 if ($value < $i) {
247 38         418 return $i-1;
248             }
249 99         1086 $value = int($value/$i);
250             }
251             }
252              
253             # ENHANCE-ME: should be able to notice rounding in $value/$i divisions of
254             # value_to_i_floor(), rather than multiplying back.
255             #
256             sub _UNTESTED__value_to_i_ceil {
257 0     0   0 my ($self, $value) = @_;
258 0 0       0 if ($value < 0) { return 0; }
  0         0  
259 0         0 my $i = $self->value_to_i_floor($value);
260 0 0       0 if ($self->ith($i) < $value) {
261 0         0 $i += 1;
262             }
263 0         0 return $i;
264             }
265              
266              
267             #--------
268             # Stirling
269             # n! ~= sqrt(2pi*n) * binomial(n,e)^n
270             # n! ~= sqrt(2*Pi) * n^(n+1/2) / e^n
271             # log(i!) ~= i*log(i) - i
272             #
273             # f(x) = x*log(x) - x - t
274             # f'(x) = log(x)
275             # sub = f(x) / f'(x)
276             # = (x*log(x) - x - t) / log(x)
277             # = x - (x+t)/log(x)
278             # new = x - sub
279             # = x - (x - (x+t)/log(x))
280             # = (x+t)/log(x)
281             #
282             # start x=t
283             # new1 = 2t/log(t)
284             # new2 = (2t/log(t) + t) / log(2t/log(t))
285             # = (2t/log(t) + t) / (log(2t) - log(log(t)))
286             #
287             # log2(i!) = log(i!)/log(2)
288             # ~= (i*log(i) - i)/log(2)
289             # log2(i!)*log(2) ~= i*log(i) - i
290             #
291             #--------
292             # Gosper, approximating terms of Stirling series
293             #
294             # n! ~= sqrt((2n+1/3)pi) * n^n * e^-n
295             #
296             sub value_to_i_estimate {
297 16     16 1 269 my ($self, $value) = @_;
298             ### value_to_i_estimate: $value
299              
300 16 100       30 if ($value <= 1) {
301 10         22 return 0;
302             }
303 6 50       82 if ($value <= 3) {
304 0         0 return 1;
305             }
306              
307 6         83 my $t;
308 6 100       19 if (defined (my $blog2 = _blog2_estimate($value))) {
309 1         313 $t = $blog2 * log(2);
310             } else {
311 5         14 $t = log($value);
312             }
313              
314             # two steps of Newton's method
315 6         12 my $x = 2*$t/log($t);
316 6         18 return int(($x+$t)/log($x));
317              
318              
319             # # single step of Newton's method starting x=t
320             # # x-1 is a touch under the true i, so just int() down
321             # return int(2*$t/log($t));
322             #
323             # multiple steps
324             # my $x = $t;
325             # for (1 .. 10) {
326             # ### $x
327             # ### log: log($x)
328             # ### f: ($x*log($x)-$x - $t)
329             # ### fd: log($x)
330             # ### sub: ($x*log($x) - $x - $t)/log($x)
331             # ### new: ($x+$t)/log($x)
332             #
333             # $x = ($x+$t)/log($x);
334             # }
335             # return int($x)-1;
336             }
337              
338             1;
339             __END__