File Coverage

blib/lib/Math/NumSeq/SqrtContinued.pm
Criterion Covered Total %
statement 104 108 96.3
branch 13 16 81.2
condition 6 9 66.6
subroutine 20 20 100.0
pod 6 6 100.0
total 149 159 93.7


line stmt bran cond sub pod time code
1             # Copyright 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             # cf Knuth volume 2 Seminumerical Algorithms section 4.5.3 exercise 12.
20             #
21              
22             package Math::NumSeq::SqrtContinued;
23 2     2   14185 use 5.004;
  2         10  
  2         144  
24 2     2   13 use strict;
  2         5  
  2         84  
25              
26 2     2   10 use vars '$VERSION', '@ISA';
  2         4  
  2         150  
27             $VERSION = 71;
28 2     2   558 use Math::NumSeq 7; # v.7 for _is_infinite()
  2         35  
  2         107  
29             @ISA = ('Math::NumSeq');
30             *_is_infinite = \&Math::NumSeq::_is_infinite;
31              
32 2     2   14 use List::Util 'min','max';
  2         5  
  2         163  
33              
34 2     2   2650 use Math::NumSeq::Squares;
  2         7  
  2         73  
35 2     2   3201 use Math::NumSeq::SqrtContinuedPeriod;
  2         7  
  2         71  
36              
37             # uncomment this to run the ### lines
38             #use Smart::Comments;
39              
40              
41             # use constant name => Math::NumSeq::__('Sqrt Continued Fraction');
42 2     2   10 use constant description => Math::NumSeq::__('Continued fraction expansion of a square root.');
  2         5  
  2         7  
43 2     2   10 use constant default_i_start => 0;
  2         4  
  2         65  
44 2     2   11 use constant characteristic_smaller => 1;
  2         2  
  2         70  
45 2     2   10 use constant characteristic_increasing => 0;
  2         3  
  2         73  
46             # use constant characteristic_continued_fraction => 1;
47              
48 2     2   510 use Math::NumSeq::SqrtDigits;
  2         7  
  2         231  
49 2         26 use constant parameter_info_array =>
50             [
51             Math::NumSeq::SqrtDigits->parameter_info_hash->{'sqrt'},
52 2     2   17 ];
  2         5  
53              
54             #------------------------------------------------------------------------------
55              
56             # http://oeis.org/index/Con#confC
57             #
58             my @oeis_anum = (
59             # A010171 to A010175 have OFFSET=1, unlike the rest
60             # OFFSET=0, but still include them in the catalogue for now
61              
62             # OEIS-Catalogue array begin
63             undef, # sqrt=0
64             undef, # sqrt=1
65             'A040000', # sqrt=2
66             'A040001', # sqrt=3
67             undef, # sqrt=4
68             'A040002', # sqrt=5
69             'A040003', # sqrt=6
70             'A010121', # sqrt=7
71             'A040005', # sqrt=8
72             undef, # sqrt=9
73              
74             'A040006', # sqrt=10
75             'A040007', # sqrt=11
76             'A040008', # sqrt=12
77             'A010122', # sqrt=13
78             'A010123', # sqrt=14
79             'A040011', # sqrt=15
80             undef, # sqrt=16
81             'A040012', # sqrt=17
82             'A040013', # sqrt=18
83             'A010124', # sqrt=19
84              
85             'A040015', # sqrt=20
86             'A010125', # sqrt=21
87             'A010126', # sqrt=22
88             'A010127', # sqrt=23
89             'A040019', # sqrt=24
90             undef, # sqrt=25
91             'A040020', # sqrt=26
92             'A040021', # sqrt=27
93             'A040022', # sqrt=28
94             'A010128', # sqrt=29
95              
96             'A040024', # sqrt=30
97             'A010129', # sqrt=31
98             'A010130', # sqrt=32
99             'A010131', # sqrt=33
100             'A010132', # sqrt=34
101             'A040029', # sqrt=35
102             undef, # sqrt=36
103             'A040030', # sqrt=37
104             'A040031', # sqrt=38
105             'A040032', # sqrt=39
106              
107             'A040033', # sqrt=40
108             'A010133', # sqrt=41
109             'A040035', # sqrt=42
110             'A010134', # sqrt=43
111             'A040037', # sqrt=44
112             'A010135', # sqrt=45
113             'A010136', # sqrt=46
114             'A010137', # sqrt=47
115             'A040041', # sqrt=48
116             undef, # sqrt=49
117              
118             'A040042', # sqrt=50
119             'A040043', # sqrt=51
120             'A010138', # sqrt=52
121             'A010139', # sqrt=53
122             'A010140', # sqrt=54
123             'A010141', # sqrt=55
124             'A040048', # sqrt=56
125             'A010142', # sqrt=57
126             'A010143', # sqrt=58
127             'A010144', # sqrt=59
128              
129             'A040052', # sqrt=60
130             'A010145', # sqrt=61
131             'A010146', # sqrt=62
132             'A040055', # sqrt=63
133             undef, # sqrt=64
134             'A040056', # sqrt=65
135             'A040057', # sqrt=66
136             'A010147', # sqrt=67
137             'A040059', # sqrt=68
138             'A010148', # sqrt=69
139              
140             'A010149', # sqrt=70
141             'A010150', # sqrt=71
142             'A040063', # sqrt=72
143             'A010151', # sqrt=73
144             'A010152', # sqrt=74
145             'A010153', # sqrt=75
146             'A010154', # sqrt=76
147             'A010155', # sqrt=77
148             'A010156', # sqrt=78
149             'A010157', # sqrt=79
150              
151             'A040071', # sqrt=80
152             undef, # sqrt=81
153             'A040072', # sqrt=82
154             'A040073', # sqrt=83
155             'A040074', # sqrt=84
156             'A010158', # sqrt=85
157             'A010159', # sqrt=86
158             'A040077', # sqrt=87
159             'A010160', # sqrt=88
160             'A010161', # sqrt=89
161              
162             'A040080', # sqrt=90
163             'A010162', # sqrt=91
164             'A010163', # sqrt=92
165             'A010164', # sqrt=93
166             'A010165', # sqrt=94
167             'A010166', # sqrt=95
168             'A010167', # sqrt=96
169             'A010168', # sqrt=97
170             'A010169', # sqrt=98
171             'A010170', # sqrt=99
172              
173             undef, # sqrt=100
174             undef, # sqrt=101, is 10, 20,20,rep
175             undef, # sqrt=102, is 10, 10,20,10,20,rep
176             'A010171', # sqrt=103
177             undef, # sqrt=104, is 10, 5,20,5,20,rep
178             undef, # sqrt=105
179             'A010172', # sqrt=106
180             'A010173', # sqrt=107
181             'A010174', # sqrt=108
182             'A010175', # sqrt=109
183              
184             undef, # sqrt=110
185             'A010176', # sqrt=111
186             'A010177', # sqrt=112
187             'A010178', # sqrt=113
188             'A010179', # sqrt=114
189             'A010180', # sqrt=115
190             'A010181', # sqrt=116
191             'A010182', # sqrt=117
192             'A010183', # sqrt=118
193             'A010184', # sqrt=119
194              
195             undef, # sqrt=120
196             undef, # sqrt=121
197             undef, # sqrt=122
198             undef, # sqrt=123
199             'A010185', # sqrt=124
200             'A010186', # sqrt=125
201             'A010187', # sqrt=126
202             'A010188', # sqrt=127
203             'A010189', # sqrt=128
204             'A010190', # sqrt=129
205              
206             undef, # sqrt=130
207             'A010191', # sqrt=131
208             undef, # sqrt=132
209             'A010192', # sqrt=133
210             'A010193', # sqrt=134
211             'A010194', # sqrt=135
212             'A010195', # sqrt=136
213             'A010196', # sqrt=137
214             'A010197', # sqrt=138
215             'A010198', # sqrt=139
216              
217             'A010199', # sqrt=140
218             'A010200', # sqrt=141
219             'A010201', # sqrt=142
220             undef, # sqrt=143
221             undef, # sqrt=144
222             undef, # sqrt=145
223             undef, # sqrt=146
224             undef, # sqrt=147
225             undef, # sqrt=148
226             'A010202', # sqrt=149
227              
228             undef, # sqrt=150
229             'A010203', # sqrt=151
230             undef, # sqrt=152
231             'A010204', # sqrt=153
232             'A010205', # sqrt=154
233             undef, # sqrt=155
234             undef, # sqrt=156
235             'A010206', # sqrt=157
236             'A010207', # sqrt=158
237             'A010208', # sqrt=159
238              
239             'A010209', # sqrt=160
240             'A010210', # sqrt=161
241             'A010211', # sqrt=162
242             'A010212', # sqrt=163
243             undef, # sqrt=164
244             'A010213', # sqrt=165
245             'A010214', # sqrt=166
246             'A010215', # sqrt=167
247             undef, # sqrt=168
248             undef, # sqrt=169
249              
250             undef, # sqrt=170
251             undef, # sqrt=171
252             'A010216', # sqrt=172
253             'A010217', # sqrt=173
254             'A010218', # sqrt=174
255             'A010219', # sqrt=175
256             'A010220', # sqrt=176
257             'A010221', # sqrt=177
258             'A010222', # sqrt=178
259             'A010223', # sqrt=179
260              
261             undef, # sqrt=180
262             'A010224', # sqrt=181
263             undef, # sqrt=182
264             'A010225', # sqrt=183
265             'A010226', # sqrt=184
266             'A010227', # sqrt=185
267             'A010228', # sqrt=186
268             'A010229', # sqrt=187
269             'A010230', # sqrt=188
270             'A010231', # sqrt=189
271              
272             'A010232', # sqrt=190
273             'A010233', # sqrt=191
274             'A010234', # sqrt=192
275             'A010235', # sqrt=193
276             'A010236', # sqrt=194
277             undef, # sqrt=195
278             undef, # sqrt=196
279             undef, # sqrt=197
280             undef, # sqrt=198
281             'A010237', # sqrt=199
282             # OEIS-Catalogue array end
283             );
284              
285             sub oeis_anum {
286 1     1 1 6 my ($self) = @_;
287 1         9 return $oeis_anum[$self->{'sqrt'}];
288             }
289              
290             #------------------------------------------------------------------------------
291              
292             sub values_min {
293 200     200 1 829 my ($self) = @_;
294 200         383 _values_min_max($self);
295 200         750 return $self->{'values_min'};
296             }
297             sub values_max {
298 199     199 1 749 my ($self) = @_;
299 199         420 _values_min_max($self);
300 199         445 return $self->{'values_max'};
301             }
302             sub _values_min_max {
303 399     399   516 my ($self) = @_;
304 399 100       1152 return if defined $self->{'values_min'};
305              
306 200   100     1223 my $period = ($self->{'period'}
307             ||= Math::NumSeq::SqrtContinuedPeriod->ith($self->{'sqrt'}));
308 200         365 my $values_min = $self->{'root'};
309 200         485 my $values_max = $self->{'root'};
310              
311 200         272 my $sqrt = $self->{'sqrt'};
312 200         254 my $root = $self->{'root'};
313 200         246 my $p = $root;
314 200         372 my $q = $sqrt - $root*$root;
315 200         424 while ($period-- > 0) {
316 1086         1899 my $value = int (($root + $p) / $q);
317 1086         1276 $p = $value*$q - $p;
318 1086         1275 $q = ($sqrt - $p*$p) / $q;
319 1086         1612 $values_min = min($values_min, $value);
320 1086         3041 $values_max = max($values_max, $value);
321             }
322 200         336 $self->{'values_min'} = $values_min;
323 200         1198 $self->{'values_max'} = $values_max;
324             }
325              
326             #------------------------------------------------------------------------------
327              
328             # V = floor[ (P+sqrt(S))/Q ]
329             #
330             # (P+sqrt(S))/Q = V + 1/x
331             # 1/x = (P+sqrt(S) - VQ)/Q
332             # x = Q/(P+sqrt(S) - VQ)
333             # = Q/( sqrt(S) + (P-VQ))
334             # = Q*( sqrt(S) - (P-VQ)) / ( S - (P-VQ)^2)
335             # newP = VQ-P
336             # newQ = (S - (P-VQ)^2)/Q
337             # = (S- (P^2 - 2PVQ + VVQQ))/Q
338             # = (S - P^2 + 2PVQ - VVQQ)/Q
339             # = (S - P^2)/Q + (2PVQ - VVQQ)/Q
340             # = (S - P^2)/Q + 2PV - VVQ
341             #
342             # T = (S-P^2)/Q
343             # newQ = T + 2PV - VVQ
344             # newT = (S-newP^2)/newQ
345             # = (S-VQ+P)/(T + 2PV - VVQ)
346             #
347             sub rewind {
348 202     202 1 827 my ($self) = @_;
349 202         728 $self->{'i'} = $self->i_start;
350              
351 202         381 my $sqrt = $self->{'sqrt'};
352 202 50       377 if ($sqrt <= 0) {
353 0         0 $self->{'a'} = 0;
354             } else {
355             # ENHANCE-ME: 'root' and 'perfect_square' one-off in new()
356 202         454 my $root = $self->{'root'} = sqrt($sqrt);
357 202         291 my $int = int($root);
358 202 100       738 if ($root == $int) {
359 13         24 $self->{'perfect_square'} = 1;
360 13         42 $self->{'P'} = $root;
361             } else {
362 189         363 $self->{'P'} = 0;
363 189         365 $self->{'Q'} = 1;
364 189         558 $self->{'root'} = $int;
365             }
366             }
367             }
368             sub next {
369 18638     18638 1 105292 my ($self) = @_;
370             ### SqrtContinued next() ...
371              
372 18638         28257 my $p = $self->{'P'};
373 18638         19004 my $value;
374 18638 100       38769 if ($self->{'perfect_square'}) {
375 26 100       42 if (defined $p) {
376 13         37 delete $self->{'P'};
377 13         35 return (1, $p);
378             } else {
379             # perfect square no more terms
380 13         36 return;
381             }
382             }
383              
384             # always "+ 1" to round up because sqrt() is not an integer so the
385             # numerator is not divisible by the denominator
386             #
387 18612         26918 my $q = $self->{'Q'};
388 18612         29933 $value = int (($self->{'root'} + $p) / $q);
389              
390             ### $p
391             ### $q
392             ### $value
393              
394 18612         22094 $p -= $value*$q;
395 18612         25372 $self->{'P'} = -$p;
396 18612         30343 $self->{'Q'} = ($self->{'sqrt'} - $p*$p) / $q;
397              
398             ### assert: $self->{'P'} >= 0
399             ### assert: $self->{'Q'} >= 0
400             ### assert: $self->{'P'} <= $self->{'root'}
401             ### assert: $self->{'Q'} <= 2*$self->{'root'}+1
402             ### assert: (($self->{'P'} * $self->{'P'} - $self->{'sqrt'}) % $self->{'Q'}) == 0
403              
404 18612         62877 return ($self->{'i'}++, $value);
405             }
406              
407             # initial
408             # P=0 Q=1
409             # value = (root+P)/Q = root
410             # P=value*Q = root
411             # Q = (S - P*P)/Q = S-P*P
412             sub ith {
413 21     21 1 89 my ($self, $i) = @_;
414              
415 21         38 my $root = $self->{'root'};
416 21 100       52 if ($i == 0) {
417 2         9 return $root;
418             }
419              
420 19 50 33     81 if ($self->{'perfect_square'} || _is_infinite($i)) {
421 0         0 return undef;
422             }
423              
424 19   66     59 my $period = ($self->{'period'}
425             ||= Math::NumSeq::SqrtContinuedPeriod->ith($self->{'sqrt'}));
426 19         32 $i = ($i - 1) % $period;
427              
428 19         33 my $sqrt = $self->{'sqrt'};
429 19         26 my $p = $root;
430 19         29 my $q = $sqrt - $root*$root;
431 19         23 for (;;) {
432 19         34 my $value = int (($root + $p) / $q);
433 19 50       46 if (--$i < 0) {
434 19         69 return $value;
435             }
436 0           $p = $value*$q - $p;
437 0           $q = ($sqrt - $p*$p) / $q;
438             }
439             }
440              
441             1;
442             __END__