File Coverage

blib/lib/Math/NumSeq/Pell.pm
Criterion Covered Total %
statement 67 75 89.3
branch 13 18 72.2
condition 2 3 66.6
subroutine 17 18 94.4
pod 5 5 100.0
total 104 119 87.3


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::Pell;
19 2     2   7013 use 5.004;
  2         5  
20 2     2   6 use strict;
  2         2  
  2         36  
21              
22 2     2   6 use vars '$VERSION', '@ISA';
  2         2  
  2         97  
23             $VERSION = 72;
24 2     2   953 use Math::NumSeq::Base::Sparse;
  2         3  
  2         67  
25             @ISA = ('Math::NumSeq::Base::Sparse');
26              
27 2     2   8 use Math::NumSeq 7; # v.7 for _is_infinite()
  2         19  
  2         82  
28             *_is_infinite = \&Math::NumSeq::_is_infinite;
29             *_to_bigint = \&Math::NumSeq::_to_bigint;
30              
31             # uncomment this to run the ### lines
32             # use Smart::Comments;
33              
34              
35             # a(r+s) = a(r)*a(s+1) + a(r-1)*a(s)
36             #
37             # P[2k+1] = P[k]^2 + P[k+1]^2
38              
39             # C[k]^2 - 8*P[k]^2 = 4(-1)^n
40              
41             # use constant name => Math::NumSeq::__('Pell Numbers');
42 2     2   7 use constant description => Math::NumSeq::__('The Pell numbers 0, 1, 2, 5, 12, 29, 70, etc, being P[k]=2*P[k-1]+P[k-2] starting from 0,1.');
  2         57  
  2         8  
43 2     2   7 use constant i_start => 0;
  2         2  
  2         66  
44 2     2   6 use constant values_min => 0;
  2         2  
  2         63  
45 2     2   5 use constant characteristic_increasing => 1;
  2         2  
  2         59  
46 2     2   6 use constant characteristic_integer => 1;
  2         1  
  2         69  
47              
48             # cf A001333 cont frac numerators, being P[n]+P[n-1]
49             # A002203 Pell companion
50             # A077985 P[-n] negatives
51             # A099011 Pell pseudoprimes
52             # Pell(N) == kronecker(2,N) mod N for all primes and some pseudos
53             # A048739 cumulative Pell, starting at value=1
54             #
55 2     2   19 use constant oeis_anum => 'A000129'; # pell
  2         3  
  2         900  
56              
57              
58             #------------------------------------------------------------------------------
59             # the biggest f0 for which both f0 and f1 fit into a UV, and which therefore
60             # for the next step will require BigInt
61             #
62             my $uv_limit;
63             my $uv_i_limit = -1; # index of $prev_f0
64             {
65             # Float integers too in 32 bits ?
66             # my $max = 1;
67             # for (1 .. 256) {
68             # my $try = $max*2 + 1;
69             # ### $try
70             # if ($try == 2*$max || $try == 2*$max+2) {
71             # last;
72             # }
73             # $max = $try;
74             # }
75             my $max = ~0;
76              
77             # 2*f1+f0 > max
78             # f0 > max-2*f1
79             # check max-2*f1 as the stopping point, so that if i=UV_MAX then won't
80             # overflow a UV trying to get to f1>=i
81             #
82             my $f0 = 0;
83             my $f1 = 1;
84             my $prev_f0;
85             while ($f1 <= ($max>>1) && $f0 <= $max - 2*$f1) {
86             $prev_f0 = $f0;
87             ($f1,$f0) = (2*$f1+$f0,$f1);
88             $uv_i_limit++;
89             }
90              
91             ### Pell UV limit ...
92             ### $prev_f0
93             ### $f0
94             ### $f1
95             ### ~0 : ~0
96              
97             $uv_limit = $prev_f0;
98              
99             ### $uv_limit
100             ### $uv_i_limit
101             ### ith: __PACKAGE__->ith($uv_i_limit)
102              
103             __PACKAGE__->ith($uv_i_limit) == $uv_limit
104             or warn "Oops, wrong uv_i_limit";
105             }
106              
107             sub seek_to_i {
108 18     18 1 215 my ($self, $i) = @_;
109             ### Pell rewind() ...
110 18         28 ($self->{'f0'}, $self->{'f1'}) = $self->ith_pair($i);
111 18         25 $self->{'i'} = $i;
112             }
113             sub rewind {
114 5     5 1 404 my ($self) = @_;
115             ### Pell rewind() ...
116 5         21 $self->{'i'} = $self->i_start;
117 5         5 $self->{'f0'} = 0;
118 5         11 $self->{'f1'} = 1;
119             }
120             sub next {
121 225     225 1 1682 my ($self) = @_;
122             (my $ret,
123             $self->{'f0'},
124             $self->{'f1'})
125             = ($self->{'f0'},
126             $self->{'f1'},
127 225         328 $self->{'f0'} + 2*$self->{'f1'});
128              
129 225 50       307 if ($ret == $uv_limit) {
130             ### go to bigint f1 ...
131 0         0 $self->{'f1'} = _to_bigint($self->{'f1'});
132             }
133              
134 225         287 return ($self->{'i'}++, $ret);
135             }
136              
137             # P[k-2] = P[k] - 2*P[k-1]
138             sub _UNTESTED_prev {
139 0     0   0 my ($self) = @_;
140             ($self->{'f0'},
141             $self->{'f1'},
142             my $ret)
143             = ($self->{'f0'} + 2*$self->{'f1'},
144             $self->{'f0'},
145 0         0 $self->{'f1'});
146              
147 0 0       0 if (abs($ret) == $uv_limit) {
148             ### go to bigint f1 ...
149 0         0 $self->{'f1'} = _to_bigint($self->{'f1'});
150             }
151              
152 0         0 return (--$self->{'i'}, $ret);
153             }
154              
155             # P[1] = 1
156             # P[0] = 0
157             # P[-1] = 1 so 1+2*0 == 1
158             # P[-2] = -2 so -2+2*1 == 0
159             # P[-3] = 5 so 5+2*-2 == 1
160             # so P[-i] = P[i] when i even, or -P[i] when i odd
161             #
162             sub ith {
163 116     116 1 1816 my ($self, $i) = @_;
164             ### ith(): $i
165 116 50       157 if (_is_infinite($i)) {
166 0         0 return $i;
167             }
168              
169 116         87 my $neg;
170 116 100       141 if ($i < 0) {
171 16         18 $i = -$i;
172 16         16 $neg = ($i % 2 == 0);
173             }
174             ### $neg
175              
176 116         83 my $f0 = ($i * 0); # inherit bignum 0
177 116         90 my $f1 = $f0 + 1; # inherit bignum 1
178              
179 116 100 66     161 if ($i > $uv_i_limit && ! ref $f0) {
180             ### automatic BigInt as not another bignum class ...
181 2         10 $f0 = _to_bigint($f0);
182 2         155 $f1 = _to_bigint($f1);
183             }
184              
185             # ENHANCE-ME: use one of the powering algorithms
186 116         182 while ($i-- > 0) {
187             ### at: "i=$i $f0, $f1"
188 1139         24461 ($f0,$f1) = ($f1, $f0 + 2*$f1);
189             }
190             ### final: "f0=$f0, f1=$f1"
191 116 100       492 return ($neg ? -$f0 : $f0);
192             }
193              
194             # P[i] = ( (1+sqrt(2))^i - (1-sqrt(2))^i ) / (2*sqrt(2))
195             # log(P[i]) ~= i*log(1+sqrt(2)) - log(2*sqrt(2))
196             # i = (log(P[i]) + log(2*sqrt(2))) / log(1+sqrt(2))
197              
198 2     2   1131 use Math::NumSeq::Fibonacci;
  2         2  
  2         238  
199             *_blog2_estimate = \&Math::NumSeq::Fibonacci::_blog2_estimate;
200              
201             sub value_to_i_estimate {
202 27     27 1 339 my ($self, $value) = @_;
203             ### Pell value_to_i_estimate(): "$value"
204              
205 27 50       45 if (_is_infinite($value)) {
206 0         0 return $value;
207             }
208 27 100       337 if ($value <= 0) {
209 8         9 return 0;
210             }
211              
212 19 100       145 if (defined (my $blog2 = _blog2_estimate($value))) {
213             ### $blog2
214 1         322 return int( ($blog2 + (log(2*sqrt(2))/log(2)))
215             / (log(1+sqrt(2))/log(2)) );
216             }
217 18         52 return int( (log($value) + log(2*sqrt(2)))
218             / log(1+sqrt(2)) );
219             }
220              
221             1;
222             __END__