File Coverage

blib/lib/Math/NumSeq/SqrtEngel.pm
Criterion Covered Total %
statement 64 67 95.5
branch 8 12 66.6
condition n/a
subroutine 16 17 94.1
pod 4 4 100.0
total 92 100 92.0


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             package Math::NumSeq::SqrtEngel;
19 2     2   7225 use 5.004;
  2         6  
20 2     2   5 use strict;
  2         2  
  2         40  
21              
22 2     2   5 use vars '$VERSION', '@ISA';
  2         3  
  2         89  
23             $VERSION = 72;
24 2     2   344 use Math::NumSeq;
  2         4  
  2         47  
25             @ISA = ('Math::NumSeq');
26              
27 2     2   316 use Math::NumSeq::Squares;
  2         8  
  2         49  
28              
29             # uncomment this to run the ### lines
30             #use Smart::Comments;
31              
32              
33             # use constant name => Math::NumSeq::__('Sqrt Engel Expansion');
34 2     2   7 use constant description => Math::NumSeq::__('Engel expansion for a square root.');
  2         1  
  2         6  
35 2     2   6 use constant characteristic_increasing => 0; # in general
  2         2  
  2         68  
36 2     2   5 use constant characteristic_non_decreasing => 1;
  2         3  
  2         61  
37 2     2   6 use constant characteristic_integer => 1;
  2         3  
  2         61  
38 2     2   6 use constant i_start => 1;
  2         2  
  2         60  
39              
40 2     2   359 use Math::NumSeq::SqrtDigits;
  2         3  
  2         113  
41             use constant parameter_info_array =>
42             [
43 2         12 Math::NumSeq::SqrtDigits->parameter_info_hash->{'sqrt'},
44 2     2   11 ];
  2         3  
45              
46 2     2   8 use constant values_min => 1;
  2         2  
  2         608  
47             sub values_max {
48 0     0 1 0 my ($self) = @_;
49 0 0       0 return ($self->Math::NumSeq::Squares::pred($self->{'sqrt'})
50             ? 1 # perfect square, only some 1s
51             : undef);
52             }
53              
54             # cf A028259 of phi=(sqrt(5)+1)/2 golden ratio
55             # A068388 of sqrt(3/2)
56             # A059178 cube root 2
57             # A059179 cube root 3
58             #
59             my @oeis_anum = (
60             # OEIS-Catalogue array begin
61             undef, # 0
62             undef, # 1
63             'A028254', # # 2
64             'A028257', # sqrt=3
65             undef, # 4
66             'A059176', # sqrt=5
67             undef, # 6
68             'A161368', # sqrt=7
69             undef, # 8
70             undef, # 9
71             'A059177', # sqrt=10
72             # OEIS-Catalogue array end
73             );
74             sub oeis_anum {
75 3     3 1 7 my ($self) = @_;
76 3         4 return $oeis_anum[$self->{'sqrt'}];
77             }
78              
79             # a/b + 1/(b*v) = sqrt(s)
80             # a + 1/v = b*sqrt(s)
81             # 1/v = b*sqrt(s) - a
82             # v = 1 / (b*sqrt(s) - a)
83             # = (b*sqrt(s) + a) / (b*sqrt(s) - a)*(b*sqrt(s) + a)
84             # = (b*sqrt(s) + a) / (b^2*s - a^2)
85             # = (sqrt(s*b^2) + a) / (s*b^2 - a^2)
86             # round up v
87             # sqrt(sb2) never an integer
88             # bigint sqrt rounds down so +1 to round up
89             # division add den-1 to round up
90             # so add 1+den-1=den means
91             # v = floor( (floor(sqrt(s*b^2)) + a) / (s*b^2 - a^2) ) + 1
92             #
93             # a/b + 1/bv < sqrt(s)
94             # a + 1/v < b*sqrt(s)
95             # a^2 + 2a/v + 1/v^2 < s*b^2
96             # a^2*v^2 + 2a*v + 1 < s*b^2*v^2
97             #
98             # a/b < sqrt(s)
99             # a < b*sqrt(s)
100             # a^2 < s*b^2
101             #
102             sub rewind {
103 29     29 1 1032 my ($self) = @_;
104 29         61 $self->{'i'} = $self->i_start;
105              
106 29         31 my $sqrt = $self->{'sqrt'};
107 29 50       51 if ($sqrt <= 0) {
108 0         0 $self->{'a'} = 0;
109             } else {
110 29         33 my $root = sqrt($sqrt);
111 29 100       54 if ($root == int($root)) {
112             # perfect square
113 9         9 $self->{'perfect_square'} = 1;
114 9         21 $self->{'a'} = $root;
115             } else {
116             # start 0/1 a=0,b=1, so sb2=s*b^2=s
117 20         45 $self->{'a'} = Math::NumSeq::_to_bigint(0);
118 20         1165 $self->{'sb2'} = Math::NumSeq::_to_bigint($sqrt);
119             }
120             }
121             }
122             sub next {
123 192     192 1 137874 my ($self) = @_;
124             ### SqrtEngel next() ...
125              
126 192         218 my $a = $self->{'a'};
127 192         137 my $value;
128 192 100       296 if ($self->{'perfect_square'}) {
129 22 100       28 if ($a) {
130 19         17 $value = 1;
131 19         19 $self->{'a'} -= 1;
132             } else {
133             # perfect square no more terms
134 3         6 return;
135             }
136             } else {
137             ### a: "$self->{'a'}"
138             ### sb2: "$self->{'sb2'}"
139             ### b: "".sqrt($self->{'sb2'} / $self->{'sqrt'})
140             ### assert: $self->{'a'} * $self->{'a'} <= $self->{'sb2'}
141             ### num: (sqrt($self->{'sb2'}) + $self->{'a'}).''
142             ### den: ($self->{'sb2'} - $self->{'a'}**2).''
143              
144             # always "+ 1" to round up because sqrt() is not an integer so the
145             # numerator is not divisible by the denominator
146             #
147 170         362 $value = (sqrt($self->{'sb2'}) + $a) / ($self->{'sb2'} - $a*$a) + 1;
148 170         53018 $self->{'a'} = $a*$value + 1;
149 170         17768 $self->{'sb2'} *= $value * $value;
150              
151             ### new value: "$value"
152             ### assert: $self->{'a'} * $self->{'a'} <= $self->{'sb2'}
153              
154 170 50       10658 if ($value <= ~0) {
155 170         9700 $value = $value->numify;
156             }
157             }
158              
159 189         1793 return ($self->{'i'}++, $value);
160             }
161              
162             # ### assert: $self->{'a'} ** 2 * $value ** 2 + 2 * $self->{'a'} * $value + 1 < $self->{'sb2'} * $value ** 2
163              
164              
165             # a/b + 1/bv + 1/bvw < sqrt(s)
166             # 1/bvw < sqrt(s) - a/b - 1/bv
167             # 1/w < bv*sqrt(s) - bv*a/b - bv/bv
168             # 1/w < bv*sqrt(s) - v*a - 1
169             # 1/w < v*(b*sqrt(s) - a) - 1
170             #
171             # 1/bv + 1/bvw
172             # = 1/bv + 1/bvw + 1/bw - 1/bw
173             # = 1/bw + 1/bvw + (1/bv - 1/bw)
174             # = 1/bw + 1/bvw + w/bvw - v/bvw
175             # = 1/bw + (1+w-v)/bvw
176             #
177             # 1/bv+1/bv1, should have had smaller v at previous stage
178             # 1/bv < sqrt 1/b(v-1) > sqrt
179             # sqrt-1/bv > 0
180             #
181             # 1/bv+1/bv(v-1)
182             # = 1/bv * (1 + 1/(v-1))
183             #
184             # 1/bv < t < 1/b(v-1)
185             # 1/bv + 1/bvw < t < 1/bv + 1/bv(w-1)
186             # 1/bvw < t-1/bv < 1/bv(w-1)
187             #
188             # sqrt(2) - (1 + 1/3 + 1/3*5) > 0
189             # sqrt(2) - (1 + 1/3 + 1/3*6)
190             # = sqrt(2) - 25/18
191             #
192             # R = t - (a/b + 1/bv) > 0
193             # bvR = bvt - (bva/b + bv/bv)
194             # = bvt - (va + 1) > 0
195             # bvwR = tbvw - (avw + w)
196             # bvt - (va + 1) > 0
197             # bvt > (va + 1)
198              
199             # t - (a/b + 1/b(v-1)) < 0
200             # bvt - (bva/b + bv/b(v-1)) < 0
201             # bvt - (va + v/(v-1)) < 0
202             # bvt < (va + v/(v-1))
203             #
204             # S = t - (a/b + 1/bv + 1/bvw) > 0
205             # bvwS = tbvw - (bvwa/b + bvw/bv + bvw/bvw)
206             # = tbvw - (vwa + w + 1) > 0
207              
208              
209              
210             1;
211             __END__