File Coverage

blib/lib/Math/NumSeq/Catalan.pm
Criterion Covered Total %
statement 94 97 96.9
branch 19 22 86.3
condition 2 3 66.6
subroutine 19 19 100.0
pod 7 8 87.5
total 141 149 94.6


line stmt bran cond sub pod time code
1             # Copyright 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             package Math::NumSeq::Catalan;
19 2     2   6780 use 5.004;
  2         5  
20 2     2   6 use strict;
  2         2  
  2         42  
21              
22 2     2   5 use vars '$VERSION','@ISA';
  2         2  
  2         80  
23             $VERSION = 72;
24              
25 2     2   341 use Math::NumSeq;
  2         7  
  2         65  
26             @ISA = ('Math::NumSeq');
27             *_is_infinite = \&Math::NumSeq::_is_infinite;
28              
29 2     2   369 use Math::NumSeq::Fibonacci;
  2         2  
  2         57  
30             *_blog2_estimate = \&Math::NumSeq::Fibonacci::_blog2_estimate;
31              
32             # uncomment this to run the ### lines
33             #use Smart::Comments;
34              
35              
36             # use constant name => Math::NumSeq::__('Catalan Numbers');
37 2     2   7 use constant values_min => 1;
  2         2  
  2         131  
38 2     2   7 use constant default_i_start => 0;
  2         3  
  2         76  
39 2     2   6 use constant characteristic_integer => 1;
  2         3  
  2         70  
40 2     2   6 use constant characteristic_non_decreasing => 1;
  2         3  
  2         301  
41             {
42             my %characteristic_increasing_from_i = (C => 1,
43             odd => 2);
44             sub characteristic_increasing_from_i {
45 2     2 0 3 my ($self) = @_;
46 2         8 return $characteristic_increasing_from_i{$self->{'values_type'}};
47             }
48             }
49             {
50             my %description = (C => Math::NumSeq::__('The Catalan numbers 1, 1, 2, 5, 14, 42, ... (2n)!/(n!*(n+1)!).'),
51             odd => Math::NumSeq::__('The odd part of the Catalan numbers 1, 1, 2, 5, 14, 42, ... (2n)!/(n!*(n+1)!).'),);
52             sub description {
53 4     4 1 8 my ($self) = @_;
54 4 100       10 return $description{ref $self ? $self->{'values_type'} : 'C'};
55             }
56             }
57              
58 2         7 use constant parameter_info_array =>
59             [ {
60             name => 'values_type',
61             share_key => 'values_type_Codd',
62             type => 'enum',
63             default => 'C',
64             choices => ['C',
65             'odd',
66             ],
67             choices_display => [Math::NumSeq::__('C'),
68             Math::NumSeq::__('Odd'),
69             ],
70             description => Math::NumSeq::__('The Catalan numbers, or just the odd part.'),
71             },
72 2     2   10 ];
  2         2  
73              
74             #------------------------------------------------------------------------------
75             # A048990 Catalans at even i
76             # A024492 Catalans at odd i
77             # A014137 Catalans cumulative
78             # A094639 Catalans squared cumulative
79             # A000984 central binomial coeff (2n)! / n!^2
80             # A048881 trailing zeros a(n) = A000120(n+1) - 1 = onebits(n+1) - 1
81             #
82             my %oeis_anum = (C => 'A000108',
83             odd => 'A098597', # Catalan odd part, divide out powers-of-2
84             # OEIS-Catalogue: A000108
85             # OEIS-Catalogue: A098597 values_type=odd
86             );
87             sub oeis_anum {
88 2     2 1 6 my ($self) = @_;
89 2         4 return $oeis_anum{$self->{'values_type'}};
90             }
91              
92              
93             #------------------------------------------------------------------------------
94              
95 2         6 use constant 1.02 _UV_I_LIMIT => do {
96 2         1 my $uv_max = ~0 >> 1;
97             ### $uv_max
98 2         2 my $value = 1;
99 2         2 my $i = 1;
100 2         6 for (; $i++; ) {
101             ### at: "i=$i value=$value"
102 66         41 my $mul = 2*(2*$i-1);
103 66         39 my $div = $i+1;
104 66 100       75 if ($value > ($uv_max - ($uv_max%$mul)) / $mul) {
105 2         3 last;
106             }
107 64         44 $value *= $mul;
108 64         71 $value /= $div;
109             }
110             ### _UV_I_LIMIT: $i
111             ### $value
112 2         857 $i
113 2     2   7 };
  2         23  
114              
115             # use constant _NV_LIMIT => do {
116             # my $f = 1.0;
117             # my $max;
118             # for (;;) {
119             # $max = $f;
120             # my $l = 2.0*$f;
121             # my $h = 2.0*$f+2.0;
122             # $f = 2.0*$f + 1.0;
123             # $f = sprintf '%.0f', $f;
124             # last unless ($f < $h && $f > $l);
125             # }
126             # ### uv : ~0
127             # ### 53 : 1<<53
128             # ### $max
129             # $max
130             # };
131              
132              
133             # C(0) = 0!/(0!*1!) = 1
134             # C(1) = 2!/(1!*2!) = 1
135             # C(2) = 4!/(2!*3!) = 4/2 = 2
136             sub rewind {
137 10     10 1 727 my ($self) = @_;
138             ### Catalan rewind()
139 10         26 $self->{'i'} = $self->i_start;
140 10         23 $self->{'f'} = 1;
141             }
142             sub seek_to_i {
143 51     51 1 2942 my ($self, $i) = @_;
144 51         45 $self->{'i'} = $i;
145 51         78 $self->{'f'} = $self->ith($i-1);
146             }
147             # sub _UNTESTED__seek_to_value {
148             # my ($self, $value) = @_;
149             # my $i = $self->{'i'} = $self->value_to_i_ceil($value);
150             # $self->{'f'} = $self->ith($i);
151             # }
152              
153             # C(i) = C(i-1) * 2i(2i-1) / i*(i+1)
154             # = C(i-1) * 2(2i-1) / (i+1)
155             # at i=0 mul 2*(2i+1)=2
156             # div 1
157             sub next {
158 247     247 1 1799 my ($self) = @_;
159             ### Catalan next() ...
160              
161 247         171 my $i = $self->{'i'}++;
162 247 50       286 if ($i == _UV_I_LIMIT) {
163 0         0 $self->{'f'} = Math::NumSeq::_to_bigint($self->{'f'});
164             }
165 247 100       253 if ($i) {
166 239 100       251 if ($self->{'values_type'} eq 'odd') {
167 97         79 $self->{'f'} *= (2*$i-1);
168 97         60 my $div = $i+1;
169 97         117 until ($div & 1) {
170 74         84 $div >>= 1;
171             }
172             ### next f: $self->{'f'} / $div
173             ### assert: ($self->{'f'} % $div) == 0
174 97         79 $self->{'f'} /= $div;
175              
176             } else {
177 142         125 $self->{'f'} *= 2*(2*$i-1);
178             ### next f: $self->{'f'} / ($i+1)
179             ### assert: ($self->{'f'} % ($i+1)) == 0
180 142         124 $self->{'f'} /= ($i+1);
181             }
182             }
183              
184 247         277 return ($i, $self->{'f'});
185             }
186              
187             sub ith {
188 158     158 1 3150 my ($self, $i) = @_;
189             ### Catalan ith(): $i
190              
191 158 50       223 if (_is_infinite($i)) {
192 0         0 return $i;
193             }
194              
195 158         121 my $value;
196 158 100 66     380 if (! ref $i && $i >= _UV_I_LIMIT) {
197             ### use bigint ...
198 2         7 $value = Math::NumSeq::_to_bigint(1);
199             } else {
200 156         123 $value = ($i*0) + 1; # inherit bignum 1
201             }
202              
203 158 100       231 if ($self->{'values_type'} eq 'odd') {
204 73         94 foreach my $k (1 .. $i) {
205 498         7873 $value *= (2*$k-1);
206 498         7104 my $div = $k+1;
207 498         612 until ($div & 1) { $div >>= 1 }
  437         503  
208             ### assert: ($value % $div) == 0
209 498         438 $value /= $div;
210             }
211             } else {
212 85         102 foreach my $k (1 .. $i) {
213 636         8009 $value *= 2*(2*$k-1);
214             ### assert: ($value % ($k+1)) == 0
215 636         7136 $value /= ($k+1);
216             }
217             }
218              
219             ### $value
220 158         393 return $value;
221             }
222              
223             # i=0 i=1 i=2 i=3 i=4 i=5
224             # 2*1 2*3 2*5 2*7 2*9
225             # C = * --- * --- * --- * --- * ---
226             # 2 3 4 5 6
227             # C=1 C=1 C=2 C=5 C=14 C=42
228             # =2*7 =14*3
229             #
230             # C(5) = 42 = 14 * 2*(2*5-1)/6
231             #
232             # sub pred {
233             # my ($self, $value) = @_;
234             # ### Catalan pred(): $value
235             #
236             # # NV inf or nan gets $value%$i=nan and nan==0 is false.
237             # # Math::BigInt binf()%$i=0 so would go into infinite loop
238             # # hence explicit check against _is_infinite()
239             # #
240             # if (_is_infinite($value)) {
241             # return undef;
242             # }
243             #
244             # for (my $i = 2; ; $i++) {
245             # ### at: "i=$i value=$value mul ".($i+1)." div ".(2*(2*$i-1))
246             # if ($value <= 1) {
247             # return ($value == 1);
248             # }
249             # $value *= ($i+1);
250             # my $div = 2*(2*$i-1);
251             # if ($value % $div) {
252             # ### not divisible, false: "value=$value div=$div"
253             # return 0;
254             # }
255             # $value /= $div;
256             # }
257             # }
258             # =item C<$bool = $seq-Epred($value)>
259             #
260             # Return true if C<$value> is a factorial, ie. equal to C<1*2*...*i> for
261             # some i.
262              
263              
264             # sub _UNTESTED__value_to_i {
265             # my ($self, $value) = @_;
266             #
267             # if (_is_infinite($value)) {
268             # return undef;
269             # }
270             # my $i = 1;
271             # for (;;) {
272             # if ($value <= 1) {
273             # return $i;
274             # }
275             # $i++;
276             # if (($value % $i) == 0) {
277             # $value /= $i;
278             # } else {
279             # return 0;
280             # }
281             # }
282             # }
283              
284             # sub _UNTESTED__value_to_i_floor {
285             # my ($self, $value) = @_;
286             # if (_is_infinite($value)) {
287             # return $value;
288             # }
289             # if ($value < 2) {
290             # return $self->i_start;
291             # }
292             #
293             # # "/" operator converts 64-bit UV to an NV and so loses bits, making the
294             # # result come out 1 too small sometimes. Experimental switch to BigInt to
295             # # keep precision.
296             # #
297             # if (! ref $value && $value > _NV_LIMIT) {
298             # $value = Math::NumSeq::_to_bigint($value);
299             # }
300             #
301             # my $i = 2;
302             # for (;; $i++) {
303             # ### $value
304             # ### $i
305             #
306             # $value *= ($i+1);
307             # my $mul = 2*(2*$i-1);
308             # if ($value < $mul) {
309             # return $i-1;
310             # }
311             # $value = int($value/$mul);
312             # }
313             # }
314              
315             # # ENHANCE-ME: should be able to notice rounding in $value/$i divisions of
316             # # value_to_i_floor(), rather than multiplying back.
317             # #
318             # sub _UNTESTED__value_to_i_ceil {
319             # my ($self, $value) = @_;
320             # if ($value < 0) { return 0; }
321             # my $i = $self->value_to_i_floor($value);
322             # if ($self->ith($i) < $value) {
323             # $i += 1;
324             # }
325             # return $i;
326             # }
327              
328              
329             #--------
330             # Stirling approximation to n!
331             # n! ~= sqrt(2pi*n) * binomial(n,e)^n
332             # log(i!) ~= i*log(i) - i
333             #
334             # noted by Dan Fux in A000108 gives
335             # C(n) ~= 4^n / (sqrt(pi*n)*(n+1))
336             #
337             # log((2i)!/(i!(i+1)!))
338             # ~= (2i*log(2i) - 2i) - (i*log(i) - i) - ((i+1)*log(i+1) - i+1)
339             # = 2i*log(2i) - 2i - i*log(i) + i - (i+1)*log(i+1) + i+1
340             # = 2i*log(2i) - i*log(i) - (i+1)*log(i+1) + 1
341             # = 2i*(log(2)+log(i)) - i*log(i) - (i+1)*log(i+1) + 1
342             # = 2i*log(2) + 2i*log(i) - i*log(i) - (i+1)*log(i+1) + 1
343             # = 2i*log(2) + (2i-i)*log(i) - (i+1)*log(i+1) + 1
344             # ~= 2i*log(2) + (2i-i-i-1)*log(i) + 1
345             # = 2i*log(2) - log(i) + 1
346             #
347             # f(x) = 2x*log(2) - log(x) + 1 - t
348             # f'(x) = 2log(2) - log(x)
349             # sub = f(x) / f'(x)
350             # = (2x*log(2) - log(x) + 1 - t) / (2log(2) - log(x))
351             # new = x - sub
352             # = x - (2x*log(2) - log(x) + 1 - t) / (2log(2) - log(x))
353             # = ( - x*log(x) + log(x) - 1 + t) / (2log(2) - log(x))
354             # = ((1-x)*log(x) - 1 + t) / (2log(2) - log(x))
355             #
356             # start x=t
357             # new1 =
358             # new2 =
359             #------
360             #
361             # f(x) = 4^x / (sqrt(Pi * x) * (x + 1)) - targ
362             # f'(x) = (((((((4 ^ x) * 1.38629436111989) * ((3.14 * x) ^ 0.5)) - ((4 ^ x) * ((0.5 * ((3.14 * x) ^ 0.5)) * (3.14 / (3.14 * x))))) / (3.14 * x)) * (1 + x)) - ((4 ^ x) / ((3.14 * x) ^ 0.5))) / ((1 + x) ^ 2)
363             # = ((((((4^x * 1.38629436111989) * sqrt(pi*x)) - (4^x * (0.5 * sqrt(pi*x) * 1/x))) / (pi*x)) * (1 + x)) - ((4^x) / (sqrt(pi*x)))) / ((1 + x) ^ 2)
364              
365             # ENHANCE-ME: slightly off for small values, but for big the 4^n dominates
366             sub value_to_i_estimate {
367 43     43 1 529 my ($self, $value) = @_;
368             ### value_to_i_estimate: $value
369              
370 43 100       64 if ($value <= 1) {
371 20         20 return 0;
372             }
373 23 50       144 if ($value <= 3) {
374 0         0 return 1;
375             }
376              
377 23         134 my $i = _blog2_estimate($value);
378 23 100       455 unless (defined $i) {
379 21         28 $i = log($value) * (1/log(2));
380             }
381 23         17 $i /= 2;
382              
383 23         26 return int($i);
384             }
385              
386             1;
387             __END__