File Coverage

blib/lib/Math/NumSeq/AlgebraicContinued.pm
Criterion Covered Total %
statement 98 105 93.3
branch 14 22 63.6
condition n/a
subroutine 16 16 100.0
pod 3 3 100.0
total 131 146 89.7


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::AlgebraicContinued;
19 1     1   920 use 5.004;
  1         2  
20 1     1   3 use strict;
  1         2  
  1         18  
21 1     1   3 use Carp;
  1         1  
  1         52  
22              
23 1     1   3 use vars '$VERSION', '@ISA';
  1         1  
  1         46  
24             $VERSION = 72;
25 1     1   3 use Math::NumSeq 7; # v.7 for _is_infinite()
  1         9  
  1         43  
26             @ISA = ('Math::NumSeq');
27             *_is_infinite = \&Math::NumSeq::_is_infinite;
28             *_to_bigint = \&Math::NumSeq::_to_bigint;
29              
30 1     1   628 use Math::NumSeq::Cubes;
  1         2  
  1         27  
31              
32             # uncomment this to run the ### lines
33             #use Smart::Comments;
34              
35              
36             # use constant name => Math::NumSeq::__('Algebraic Continued Fraction');
37 1     1   4 use constant description => Math::NumSeq::__('Continued fraction expansion of an algebraic number, such as cube root or nth root.');
  1         1  
  1         3  
38 1     1   3 use constant default_i_start => 0;
  1         2  
  1         30  
39 1     1   3 use constant characteristic_smaller => 1;
  1         1  
  1         30  
40 1     1   3 use constant characteristic_increasing => 0;
  1         1  
  1         60  
41             # use constant characteristic_continued_fraction => 1;
42              
43 1         3 use constant parameter_info_array =>
44             [
45             {
46             name => 'expression',
47             display => Math::NumSeq::__('Expression'),
48             type => 'string',
49             width => 20,
50             default => 'cbrt 2',
51             choices => ['cbrt 2','7throot 123'],
52             description => Math::NumSeq::__('Expression to take continued fraction. Can be "cbrt 123", "7throot 123", etc'),
53             },
54 1     1   3 ];
  1         1  
55              
56             #------------------------------------------------------------------------------
57              
58             # A002945 to A002949 are OFFSET=1, unlike i_start=0 here
59             # (OFFSET=0 is rumoured to be the preferred style for continued fractions)
60             #
61             my @oeis_anum;
62              
63             $oeis_anum[3]
64             = {
65             # OEIS-Catalogue array begin
66              
67             '2,0,0,-1,i_start=1' => 'A002945', # expression=cbrt2 i_start=1
68             '3,0,0,-1,i_start=1' => 'A002946', # expression=cbrt3 i_start=1
69             '4,0,0,-1,i_start=1' => 'A002947', # expression=cbrt4 i_start=1
70             '5,0,0,-1,i_start=1' => 'A002948', # expression=cbrt5 i_start=1
71             '6,0,0,-1,i_start=1' => 'A002949', # expression=cbrt6 i_start=1
72             # undef, # expression=cbrt7
73             # undef, # expression=cbrt8
74             '9,0,0,-1' => 'A010239', # expression=cbrt9
75             '10,0,0,-1' => 'A010240', # expression=cbrt10
76             '11,0,0,-1' => 'A010241', # expression=cbrt11
77             '12,0,0,-1' => 'A010242', # expression=cbrt12
78             '13,0,0,-1' => 'A010243', # expression=cbrt13
79             '14,0,0,-1' => 'A010244', # expression=cbrt14
80             '15,0,0,-1' => 'A010245', # expression=cbrt15
81             '16,0,0,-1' => 'A010246', # expression=cbrt16
82             '17,0,0,-1' => 'A010247', # expression=cbrt17
83             '18,0,0,-1' => 'A010248', # expression=cbrt18
84             '19,0,0,-1' => 'A010249', # expression=cbrt19
85             '20,0,0,-1' => 'A010250', # expression=cbrt20
86             '21,0,0,-1' => 'A010251', # expression=cbrt21
87             '22,0,0,-1' => 'A010252', # expression=cbrt22
88             '23,0,0,-1' => 'A010253', # expression=cbrt23
89             '24,0,0,-1' => 'A010254', # expression=cbrt24
90             '25,0,0,-1' => 'A010255', # expression=cbrt25
91             '26,0,0,-1' => 'A010256', # expression=cbrt26
92             # '27,0,0,-1' => undef, # 27
93             '28,0,0,-1' => 'A010257', # expression=cbrt28
94             '29,0,0,-1' => 'A010258', # expression=cbrt29
95             '30,0,0,-1' => 'A010259', # expression=cbrt30
96             '31,0,0,-1' => 'A010260', # expression=cbrt31
97             '32,0,0,-1' => 'A010261', # expression=cbrt32
98             '33,0,0,-1' => 'A010262', # expression=cbrt33
99             '34,0,0,-1' => 'A010263', # expression=cbrt34
100             '35,0,0,-1' => 'A010264', # expression=cbrt35
101             '36,0,0,-1' => 'A010265', # expression=cbrt36
102             '37,0,0,-1' => 'A010266', # expression=cbrt37
103             '38,0,0,-1' => 'A010267', # expression=cbrt38
104             '39,0,0,-1' => 'A010268', # expression=cbrt39
105             '40,0,0,-1' => 'A010269', # expression=cbrt40
106             '41,0,0,-1' => 'A010270', # expression=cbrt41
107             '42,0,0,-1' => 'A010271', # expression=cbrt42
108             '43,0,0,-1' => 'A010272', # expression=cbrt43
109             '44,0,0,-1' => 'A010273', # expression=cbrt44
110             '45,0,0,-1' => 'A010274', # expression=cbrt45
111             '46,0,0,-1' => 'A010275', # expression=cbrt46
112             '47,0,0,-1' => 'A010276', # expression=cbrt47
113             '48,0,0,-1' => 'A010277', # expression=cbrt48
114             '49,0,0,-1' => 'A010278', # expression=cbrt49
115             '50,0,0,-1' => 'A010279', # expression=cbrt50
116             '51,0,0,-1' => 'A010280', # expression=cbrt51
117             '52,0,0,-1' => 'A010281', # expression=cbrt52
118             '53,0,0,-1' => 'A010282', # expression=cbrt53
119             '54,0,0,-1' => 'A010283', # expression=cbrt54
120             '55,0,0,-1' => 'A010284', # expression=cbrt55
121             '56,0,0,-1' => 'A010285', # expression=cbrt56
122             '57,0,0,-1' => 'A010286', # expression=cbrt57
123             '58,0,0,-1' => 'A010287', # expression=cbrt58
124             '59,0,0,-1' => 'A010288', # expression=cbrt59
125             '60,0,0,-1' => 'A010289', # expression=cbrt60
126             '61,0,0,-1' => 'A010290', # expression=cbrt61
127             '62,0,0,-1' => 'A010291', # expression=cbrt62
128             '63,0,0,-1' => 'A010292', # expression=cbrt63
129             # '64,0,0,-1' => undef, # 64
130             '65,0,0,-1' => 'A010293', # expression=cbrt65
131             '66,0,0,-1' => 'A010294', # expression=cbrt66
132             '67,0,0,-1' => 'A010295', # expression=cbrt67
133             '68,0,0,-1' => 'A010296', # expression=cbrt68
134             '69,0,0,-1' => 'A010297', # expression=cbrt69
135             '70,0,0,-1' => 'A010298', # expression=cbrt70
136             '71,0,0,-1' => 'A010299', # expression=cbrt71
137             '72,0,0,-1' => 'A010300', # expression=cbrt72
138             '73,0,0,-1' => 'A010301', # expression=cbrt73
139             '74,0,0,-1' => 'A010302', # expression=cbrt74
140             '75,0,0,-1' => 'A010303', # expression=cbrt75
141             '76,0,0,-1' => 'A010304', # expression=cbrt76
142             '77,0,0,-1' => 'A010305', # expression=cbrt77
143             '78,0,0,-1' => 'A010306', # expression=cbrt78
144             '79,0,0,-1' => 'A010307', # expression=cbrt79
145             '80,0,0,-1' => 'A010308', # expression=cbrt80
146             '81,0,0,-1' => 'A010309', # expression=cbrt81
147             '82,0,0,-1' => 'A010310', # expression=cbrt82
148             '83,0,0,-1' => 'A010311', # expression=cbrt83
149             '84,0,0,-1' => 'A010312', # expression=cbrt84
150             '85,0,0,-1' => 'A010313', # expression=cbrt85
151             '86,0,0,-1' => 'A010314', # expression=cbrt86
152             '87,0,0,-1' => 'A010315', # expression=cbrt87
153             '88,0,0,-1' => 'A010316', # expression=cbrt88
154             '89,0,0,-1' => 'A010317', # expression=cbrt89
155             '90,0,0,-1' => 'A010318', # expression=cbrt90
156             '91,0,0,-1' => 'A010319', # expression=cbrt91
157             '92,0,0,-1' => 'A010320', # expression=cbrt92
158             '93,0,0,-1' => 'A010321', # expression=cbrt93
159             '94,0,0,-1' => 'A010322', # expression=cbrt94
160             '95,0,0,-1' => 'A010323', # expression=cbrt95
161             '96,0,0,-1' => 'A010324', # expression=cbrt96
162             '97,0,0,-1' => 'A010325', # expression=cbrt97
163             '98,0,0,-1' => 'A010326', # expression=cbrt98
164             '99,0,0,-1' => 'A010327', # expression=cbrt99
165             '100,0,0,-1' => 'A010328', # expression=cbrt100
166              
167             # OEIS-Catalogue array end
168             };
169              
170             $oeis_anum[4]
171             = {
172             # OEIS-Catalogue array begin
173             '2,0,0,0,-1,i_start=1' => 'A179613', # expression=4throot2 i_start=1
174             '3,0,0,0,-1,i_start=1' => 'A179615', # expression=4throot3 i_start=1
175             '5,0,0,0,-1,i_start=1' => 'A179616', # expression=4throot5 i_start=1
176             '91,0,0,0,-10,i_start=1' => 'A093876', # expression=4throot9.1 i_start=1
177             # OEIS-Catalogue array end
178             };
179              
180             $oeis_anum[5]
181             = {
182             # OEIS-Catalogue array begin
183             '2,0,0,0,0,-1,i_start=1' => 'A002950', # expression=5throot2 i_start=1
184             '3,0,0,0,0,-1,i_start=1' => 'A003117', # expression=5throot3 i_start=1
185             '4,0,0,0,0,-1,i_start=1' => 'A003118', # expression=5throot4 i_start=1
186             '5,0,0,0,0,-1,i_start=1' => 'A002951', # expression=5throot5 i_start=1
187             # OEIS-Catalogue array end
188             };
189              
190             sub oeis_anum {
191 2     2 1 7 my ($self) = @_;
192 2         2 my $key = join(',',@{$self->{'orig_poly'}});
  2         7  
193 2 50       112 if (my $i_start = $self->i_start) {
194             ### $i_start
195 0         0 $key .= ",i_start=$i_start"; # if non-zero
196             }
197             ### $key
198 2         7 return $oeis_anum[$self->{'root'}]->{$key};
199             }
200              
201             #------------------------------------------------------------------------------
202             #
203             # (aC+b)/(cC+d) - j >= 0
204             # (aC+b) - j*(cC+d) >= 0
205             # aC + b - jcC - jd >= 0
206             # (a-jc)C >= (jd-b)
207             # CC*(a-jc)^3 >= (jd-b)^3
208             # CC*(a-jc)^3 - (jd-b)^3 >= 0
209             # (-d^3 - c^3*CC)*j^3
210             # + (3*b*d^2 + 3*c^2*a*CC)*j^2
211             # + (-3*b^2*d - 3*c*a^2*CC)*j
212             # + (b^3 + a^3*CC)
213             # poly
214             # p = -d^3 - c^3*CC
215             # q = 3*b*d^2 + 3*c^2*a*CC
216             # r = -3*b^2*d - 3*c*a^2*CC
217             # s = b^3 + a^3*CC
218             # initial a=1,b=0,c=0,d=1
219             # p = -1
220             # q = 0
221             # r = 0
222             # s = 1*CC
223             #
224             # new = 1 / ((aC+b)/(cC+d) - j)
225             # = 1 / (((aC+b) - j*(cC+d))/(cC+d))
226             # = (cC+d) / ((aC+b) - j*(cC+d))
227             # = (cC+d) / ((a-jc)*C + b-jd)
228             # new p = -(b-jd)^3 - (a-jc)^3*CC
229             # = (d^3 + c^3*CC)*j^3 + (-3*b*d^2 - 3*c^2*a*CC)*j^2 + (3*b^2*d + 3*c*a^2*CC)*j + (-b^3 - a^3*CC)
230             # = -p*j^3 + (-3*b*d^2 - 3*c^2*a*CC)*j^2 + (3*b^2*d + 3*c*a^2*CC)*j + (-b^3 - a^3*CC)
231             # = d^3*j^3 + c^3*CC*j^3
232             # + (-3*b*d^2*j^2 - 3*c^2*a*CC*j^2)
233             # + (3*b^2*d*j + 3*c*a^2*j*CC)
234             # + (-b^3 - a^3*CC)
235             # new q = 3*d*(b-jd)^2 + 3*(a-jc)^2*c*CC
236             # new r = -3*d^2*(b-jd) - 3*(a-jc)*c^2*CC
237             # new s = d^3 + c^3*CC
238             # = -p
239              
240             # x = root of p,q,r,s
241             # shift to (x+j)
242             # p*(x+j)^3 + q*(x+j)^2 + r*(x+j) + s
243             # reverse and negate for -1/x
244             # new s = -p
245             # new r = -(3*p*j + q)
246             # new q = -(3*p*j^2 + 2*q*j + r)
247             # new p = -(p*j^3 + q*j^2 + r*j + s)
248             #
249             # o*(x+j)^4 + p*(x+j)^3 + q*(x+j)^2 + r*(x+j) + s
250             # low = s + r*j + q*j^2 + p*j^3 + o*j^4
251             # next = r + q*2*j + p*3j^2 + o*4*j^3
252             # next = q + p*3j + o* 6 j^2 1,3,6,10,15,21,28
253             # next = p + o*4j 1,4,10,20,35,56
254             # high = o 1,5,15,35,70,126
255             #
256             # bin(n,m) = n!/m!(n-m)!
257             # bin(n+1,m+1) = (n+1)!/(m+1)!(n+1-m-1)!
258             # = (n+1)/(m+1) * n!/m!/(n-m)!
259             # bin(10,3)=120 120*11/4 = 330
260             # bin(11,4)=330 = 120*11/4
261             #
262             #-------------
263             # perfect cube C=2
264             # new = (C+0)/(0+1) - 2
265             # = (0+1) / ((1-2*0)C + 0-2*1)
266             # = (0+1) / (1C + -2)
267             # p = -(-2)^3 - 1^3 * CC
268             # = 0
269              
270             # nthroot(num/den)
271             # f(x) = num/den - x^n
272             # = num - den*x^n
273             sub _nthroot_to_poly {
274 2     2   2 my ($power, $num, $den) = @_;
275 2         6 my $zero = Math::NumSeq::_to_bigint(0);
276 2         138 return (_to_bigint("$num"),
277             ($zero) x ($power-1),
278             - _to_bigint("$den"));
279             }
280              
281             my %name_to_root = (sqrt => 2,
282             cbrt => 3);
283             sub rewind {
284 6     6 1 1697 my ($self) = @_;
285 6         25 $self->{'i'} = $self->i_start;
286              
287 6 100       18 if (! $self->{'orig_poly'}) {
288 2         3 my $str = $self->{'expression'};
289 2         6 $str =~ s/^\s+//;
290 2         6 $str =~ s/\s+$//;
291 2         4 $str =~ s/\s+/ /g;
292 2         2 my ($root, $operand);
293 2 50       12 if ($str =~ /^(sqrt|cbrt|(\d+)throot) ?\(? ?(\d+(\.\d*)?|\d*\.\d+) ?\)?$/) {
294 2 50       6 $root = (defined $2 ? $2 : $name_to_root{$1});
295 2         4 $operand = $3;
296             } else {
297 0         0 croak "Unrecognised expression: ",$str;
298             }
299             ### $root
300             ### $operand
301              
302 2         4 my ($num, $den);
303 2 50       5 if ($operand =~ /^(\d*)\.(\d*)$/) {
304 0         0 $num = $1.$2;
305 0         0 $den = '1' . ('0' x length($2));
306             } else {
307 2         3 $num = $operand;
308 2         2 $den = 1;
309             }
310              
311 2         6 $self->{'orig_poly'} = [ _nthroot_to_poly($root,$num,$den) ];
312             ### orig_poly join(',',@{$self->{'orig_poly'}})
313              
314             # if root<1 then initial continued fraction term is 0
315 2 50       94 $self->{'values_min'} = (_eval($self->{'orig_poly'},1) < 0 ? 0 : 1);
316              
317 2         166 $self->{'root'} = $root;
318 2         4 $self->{'operand'} = $operand;
319             }
320              
321 6         8 $self->{'poly'} = [ @{$self->{'orig_poly'}} ]; # copy array
  6         34  
322             }
323              
324             sub _eval {
325 304     304   294 my ($poly, $x) = @_;
326 304         228 my $ret = 0;
327 304         340 foreach my $coeff (reverse @$poly) { # high to low
328 1216         57797 $ret *= $x;
329 1216         61811 $ret += $coeff;
330             }
331 304         11242 return $ret;
332             }
333              
334             sub next {
335 94     94 1 1015 my ($self) = @_;
336             ### AlgebraicContinued next(): "poly=".join(',',@{$self->{'poly'}})
337              
338 94         82 my $poly = $self->{'poly'};
339 94 50       166 if ($poly->[-1] == 0) {
340             ### poly high zero, perfect power ...
341 0         0 return;
342             }
343              
344             # doubling probe for p(lo) >= 0 by p(hi) < 0
345 94         6884 my $lo = $self->{'values_min'}; # 0 or 1
346 94         71 my $hi = 2;
347 94         143 while (_eval($poly,$hi) >= 0) {
348             ### assert: _eval($poly,$lo) >= 0
349             ### lohi: "$lo,$hi eval "._eval($poly,$lo)." "._eval($poly,$hi)
350 104 50       8287 if ($hi == 0x4000_0000) {
351             # ENHANCE-ME: are terms ever bignums ?
352 0         0 $hi = _to_bigint($hi);
353             }
354 104         207 ($lo,$hi) = ($hi,2*$hi);
355             }
356              
357             # binary search for smallest c with poly(c) < 0
358 94         6922 my $c;
359 94         77 for (;;) {
360 198         304 $c = int(($lo+$hi)/2);
361              
362             ### lohi: "$lo,$hi poly "._eval($poly,$lo)." "._eval($poly,$hi)
363             ### c: "$c"
364             ### assert: _eval($poly,$lo) >= 0
365             ### assert: _eval($poly,$hi) < 0
366              
367 198 100       313 if ($c == $lo) {
368 94         91 last;
369             }
370 104 100       130 if (_eval($poly,$c) >= 0) {
371 34         2625 $lo = $c;
372             } else {
373 70         5170 $hi = $c;
374             }
375             }
376             ### c: "$c"
377              
378             # column = j row=j-i
379             # binomial(j+1,j-i+1)
380             # factor (n+1)/(m+1) = (j+2)/(j-i+2)
381             #
382 94         134 my @new = @$poly;
383             {
384 94         70 my $cpow = _to_bigint($c); # c^i
  94         198  
385 94         1789 foreach my $i (1 .. $#$poly) {
386 282         12674 my $t = $cpow;
387 282         435 foreach my $j ($i .. $#$poly-1) {
388             ### term: "t=$t coeff=$poly->[$j] next mul ".($j+1)." div ".($j+1-$i)
389 282         7505 $new[$j-$i] += $t * $poly->[$j];
390 282         21741 $t *= $j+1;
391             ### assert: $t % ($j+1-$i) == 0
392 282         20004 $t /= $j+1-$i;
393             }
394 282         14137 $new[$#$poly-$i] += $t * $poly->[-1];
395 282         22683 $cpow *= $c;
396             }
397             }
398 94         7092 @$poly = map{-$_} reverse @new;
  376         3416  
399              
400 94 50       1474 if ($poly->[-1] > 0) {
401 0         0 die "Oops, AlgebraicContinued poly not negative";
402             }
403              
404 94         7308 return ($self->{'i'}++, $c);
405              
406              
407             # $self->{'p'} = -($p*$c**3 + $q*$c**2 + $r*$c + $s);
408             # $self->{'q'} = -(3*$p*$c**2 + 2*$q*$c + $r);
409             # $self->{'r'} = -(3*$p*$c + $q);
410             # $self->{'s'} = -$p;
411             }
412              
413             1;
414             __END__