File Coverage

blib/lib/Math/NumSeq/LucasNumbers.pm
Criterion Covered Total %
statement 118 125 94.4
branch 29 36 80.5
condition 8 12 66.6
subroutine 21 21 100.0
pod 9 11 81.8
total 185 205 90.2


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             package Math::NumSeq::LucasNumbers;
19 2     2   13535 use 5.004;
  2         9  
  2         91  
20 2     2   13 use strict;
  2         4  
  2         75  
21              
22 2     2   11 use vars '$VERSION','@ISA';
  2         3  
  2         139  
23             $VERSION = 71;
24 2     2   1699 use Math::NumSeq::Base::Sparse;
  2         6  
  2         65  
25             @ISA = ('Math::NumSeq::Base::Sparse');
26              
27 2     2   10 use Math::NumSeq;
  2         4  
  2         89  
28             *_is_infinite = \&Math::NumSeq::_is_infinite;
29             *_to_bigint = \&Math::NumSeq::_to_bigint;
30              
31 2     2   652 use Math::NumSeq::Fibonacci;
  2         5  
  2         483  
32             *_blog2_estimate = \&Math::NumSeq::Fibonacci::_blog2_estimate;
33             *_bit_split_hightolow = \&Math::NumSeq::Fibonacci::_bit_split_hightolow;
34              
35              
36             # uncomment this to run the ### lines
37             # use Smart::Comments;
38              
39              
40             # use constant name => Math::NumSeq::__('Lucas Numbers');
41             sub description {
42 4     4 1 19 my ($self) = @_;
43 4 100 100     20 if (ref $self && $self->i_start == 0) {
44 1         5 return Math::NumSeq::__('Lucas numbers 2, 1, 3, 4, 7, 11, 18, 29, etc, being L(i+1) = L(i) + L(i-1) starting from 2,1. This is the same recurrence as the Fibonacci numbers, but a different starting point.');
45             } else {
46 3         12 return Math::NumSeq::__('Lucas numbers 1, 3, 4, 7, 11, 18, 29, etc, being L(i+1) = L(i) + L(i-1) starting from 1,3. This is the same recurrence as the Fibonacci numbers, but a different starting point.');
47             }
48             }
49              
50             # negatives at i odd negative, otherwise minimum 1 at i=1
51             sub values_min {
52 2     2 1 17 my ($self) = @_;
53 2         11 my $i = $self->i_start;
54 2 100 66     18 if ($i <= 0 && $i % 2 == 0) {
55             # i even, increase to make i odd and i<=1
56 1         3 $i += 1;
57             }
58 2         6 return $self->ith($i);
59             }
60              
61             sub characteristic_increasing {
62 5     5 0 9 my ($self) = @_;
63             # not increasing at i=0 value=2 then i=1 value=1
64 5         40 return ($self->i_start >= 1);
65             }
66             sub characteristic_increasing_from_i {
67 1     1 0 3 my ($self) = @_;
68 1         5 return _max($self->i_start,1);
69             }
70 2     2   12 use constant characteristic_integer => 1;
  2         3  
  2         103  
71 2     2   10 use constant default_i_start => 1;
  2         4  
  2         1784  
72              
73             sub _max {
74 1     1   3 my $ret = shift;
75 1         5 while (@_) {
76 1         1 my $next = shift;
77 1 50       6 if ($next > $ret) {
78 0         0 $ret = $next;
79             }
80             }
81 1         6 return $ret;
82             }
83              
84             #------------------------------------------------------------------------------
85             # cf A000285 starting 1,4
86             # A022086 starting 0,3
87             # A022087 starting 0,4
88             # A022095 starting 1,5
89             # A022130 starting 4,9
90             # A080023 closest to log(phi), 2,3,4,7,11
91             # A169985 nearestint(phi^n) 1,2,3,4,7,11,18
92             {
93             my %oeis_anum = (
94             # OEIS-Catalogue array begin
95             0 => 'A000032', # i_start=0 # Lucas starting at 2,1,3,...
96             1 => 'A000204', # # Lucas starting at 1,3,...
97             # OEIS-Catalogue array end
98             );
99             sub oeis_anum {
100 2     2 1 11 my ($self) = @_;
101 2         7 return $oeis_anum{$self->i_start};
102             }
103             }
104              
105             #------------------------------------------------------------------------------
106              
107             # the biggest f0 for which both f0 and f1 fit into a UV, and which therefore
108             # for the next step will require BigInt
109             #
110             my $uv_limit;
111             my $uv_i_limit;
112             {
113             # Float integers too in 32 bits ?
114             # my $max = 1;
115             # for (1 .. 256) {
116             # my $try = $max*2 + 1;
117             # ### $try
118             # if ($try == 2*$max || $try == 2*$max+2) {
119             # last;
120             # }
121             # $max = $try;
122             # }
123             my $max = ~0;
124              
125             # f1+f0 > i
126             # f0 > i-f1
127             # check i-f1 as the stopping point, so that if i=UV_MAX then won't
128             # overflow a UV trying to get to f1>=i
129             #
130             my $i = 0;
131             my $f0 = 1;
132             my $f1 = 3;
133             my $prev_f0;
134             while ($f0 <= $max - $f1) {
135             $prev_f0 = $f0;
136             ($f1,$f0) = ($f1+$f0,$f1);
137             $i++;
138             }
139             ### $prev_f0
140             ### $f0
141             ### $f1
142             ### ~0 : ~0
143              
144             $uv_limit = $prev_f0;
145             $uv_i_limit = $i;
146             ### $uv_limit
147             ### $uv_i_limit
148             ### assert: __PACKAGE__->ith($uv_i_limit) == $uv_limit
149             };
150              
151             #------------------------------------------------------------------------------
152              
153             sub rewind {
154 11     11 1 562 my ($self) = @_;
155 11         50 $self->seek_to_i($self->i_start);
156             }
157             sub seek_to_i {
158 26     26 1 334 my ($self, $i) = @_;
159 26         62 $self->{'i'} = $i;
160 26 50       81 if (abs($i) >= $uv_i_limit) {
161 0         0 $i = Math::NumSeq::_to_bigint($i);
162             }
163 26         69 $self->{'f0'} = $self->ith($i);
164 26         73 $self->{'f1'} = $self->ith($i+1);
165              
166             # or perhaps
167             # ($self->{'f0'}, $self->{'f1'}) = $self->ith_pair($i);
168             }
169              
170             sub next {
171 109     109 1 1355 my ($self) = @_;
172             ### LucasNumbers next(): "f0=$self->{'f0'}, f1=$self->{'f1'}"
173              
174 109         298 (my $ret,
175             $self->{'f0'},
176             $self->{'f1'})
177             = ($self->{'f0'},
178             $self->{'f1'},
179             $self->{'f0'}+$self->{'f1'});
180             ### $ret
181              
182 109 50       257 if ($ret == $uv_limit) {
183 0         0 $self->{'f0'} = Math::NumSeq::_to_bigint($self->{'f0'});
184 0         0 $self->{'f1'} = Math::NumSeq::_to_bigint($self->{'f1'});
185             }
186              
187 109         318 return ($self->{'i'}++, $ret);
188             }
189              
190             # F[4] = (F[2]+L[2])^2/2 - 3*F[2]^2 - 2*(-1)^2
191             # = (1+3)^2/4 - 3*1^2 - 2
192             # = 16/4 - 3 - 2
193              
194             # F[3] = ((F[1]+L[1])^2 - 2*(-1)^1)/4 + F[1]^2
195             # = ((1+3)^2 - -2)/4 + 1^2
196             # = (16 + 2)/4 + 1
197             # = (16 + 2)/4 + 1
198              
199             # F[3] = (F[1]+L[1])^2/4 + F[1]^2
200             # = (1+3)^2/4 + 1^2
201             # = 16/4 + 1
202             # = 5
203              
204             # -8,5,-3,2,-1,1, 0, 1,1,2,3,5,8,13,21
205             # -11,7,-4,3,-1, 2, 1,3,4,7,11,18,29
206             #
207              
208             sub ith {
209 90     90 1 2425 my ($self, $i) = @_;
210             ### LucasNumbers ith(): $i
211              
212 90 50       235 if (_is_infinite($i)) {
213 0         0 return $i;
214             }
215 90         141 $i = int($i);
216 90 100       181 if ($i == 0) {
217 11         31 return 2;
218             }
219              
220             # automatic BigInt if not another bignum class
221 79   66     385 my $to_bigint = ($i >= $uv_i_limit || $i <= -$uv_i_limit && ! ref $i);
222             ### $to_bigint
223              
224 79         99 my $lowzeros = 0;
225 79         201 until ($i % 2) {
226 52         56 $lowzeros++;
227 52         137 $i /= 2;
228             }
229              
230 79         174 my ($L0) = $self->ith_pair($i);
231             ### $L0
232 79 100       171 if ($to_bigint) {
233             ### to bigint ...
234 4         16 $L0 = _to_bigint($L0);
235             ### $L0
236             }
237              
238             ### apply lowzeros: $lowzeros
239 79         241 my $add = -2;
240 79         193 while ($lowzeros--) {
241 52         80 $L0 *= $L0;
242 52         723 $L0 -= $add;
243 52         1250 $add = 2;
244             ### lowzeros L0: "$L0 " . ref $L0
245             }
246              
247             ### final ith() ...
248             ### L0: "$L0"
249              
250 79         237 return $L0;
251              
252              
253             # {
254             # # cf plain interative
255             # $i--;
256             # my $f0 = ($i * 0) + 1; # inherit bignum 1
257             # my $f1 = $f0 + 2; # inherit bignum 3
258             # while (--$i > 0) {
259             # $f0 += $f1;
260             #
261             # unless (--$i > 0) {
262             # return $f0;
263             # }
264             # $f1 += $f0;
265             # }
266             # return $f1;
267             # }
268             }
269              
270             sub ith_pair {
271 112     112 1 2547 my ($self, $i) = @_;
272             ### LucasNumbers ith_pair(): $i
273              
274 112 50       455 if (_is_infinite($i)) {
275 0         0 return ($i,$i);
276             }
277 112         181 $i = int($i);
278              
279 112         290 my $neg;
280 112 100       237 if ($i < 0) {
281 22         32 $i = -1 - $i; # L[-k] = L[-1-k] * (-1)^k
282 22         35 $neg = 1;
283             }
284              
285 112         120 my ($Lk, $Lkplus1);
286 112 100       218 if ($i == 0) {
287 9         11 $Lk = 2; # L[0] = 2
288 9         12 $Lkplus1 = 1; # L[1] = 1
289              
290             } else {
291             # initial k=1
292 103 50 33     284 my $zero = ($i >= $uv_i_limit && ! ref $i
293             ? _to_bigint(0) # automatic BigInt if not another bignum class
294             : $i * 0); # inherit bignum $i
295              
296 103         141 $Lk = 1 + $zero; # L[k] = L[1] = 1
297 103         120 $Lkplus1 = 3 + $zero; # L[k+1] = L[2] = 3
298 103         125 my $add = -2; # 2*(-1)^k
299              
300              
301 103         290 my @bits = _bit_split_hightolow($i);
302             ### @bits
303 103         152 shift @bits; # drop high 1-bit
304              
305 103         422 while (@bits) {
306             ### at: "Lk=$Lk Lkplus1=$Lkplus1 bit=$bits[0]"
307 111         137 $Lk *= $Lk;
308 111         142 $Lk -= $add; # L[2k] = L[k]^2 - 2*(-1)^k
309              
310 111         141 $Lkplus1 *= $Lkplus1;
311 111         293 $Lkplus1 += $add; # L[2k+2] = L[k+1]^2 + 2*(-1)^k
312              
313             # L[2k+1] = L[2k+2] - L[2k]
314 111         144 my $Lmid = $Lkplus1 - $Lk; # L[2k+1] = L[2k+2] - L[2k]
315              
316 111 100       373 if (shift @bits) { # high to low
317 57         76 $Lk = $Lmid; # L[2k+1], L[2k+2]
318 57         330 $add = -2;
319             } else {
320 54         69 $Lkplus1 = $Lmid; # L[2k], L[2k+1]
321 54         143 $add = 2;
322             }
323             }
324             }
325             ### loop final: "Lk=$Lk Lkplus1=$Lkplus1"
326              
327 112 100       253 if ($neg) {
328 22         52 ($Lk,$Lkplus1) = ($Lkplus1, $Lk);
329 22 100       50 if ($i % 2) {
330 4         7 $Lkplus1 = -$Lkplus1;
331             } else {
332 18         25 $Lk = -$Lk;
333             }
334             }
335 112         288 return ($Lk, $Lkplus1);
336             }
337              
338 2     2   12 use constant 1.02 _PHI => (1 + sqrt(5)) / 2;
  2         33  
  2         389  
339              
340             sub value_to_i_estimate {
341 33     33 1 699 my ($self, $value) = @_;
342 33 50       82 if (_is_infinite($value)) {
343 0         0 return $value;
344             }
345 33 100       6604 if ($value <= 0) {
346 14         36 return 0;
347             }
348 19 100       392 if (defined (my $blog2 = _blog2_estimate($value))) {
349             # i ~= log2(L(i)) / log2(phi)
350             # with log2(x) = log(x)/log(2)
351 2         700 return int($blog2 * (1 / (log(_PHI)/log(2))));
352             }
353              
354             # i ~= log(L(i)) / log(phi)
355 17         53 return int(log($value) * (1/log(_PHI)));
356             }
357              
358             1;
359             __END__