File Coverage

blib/lib/Math/PlanePath/CfracDigits.pm
Criterion Covered Total %
statement 124 150 82.6
branch 24 40 60.0
condition 5 18 27.7
subroutine 25 30 83.3
pod 7 7 100.0
total 185 245 75.5


line stmt bran cond sub pod time code
1             # Copyright 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019 Kevin Ryde
2              
3             # This file is part of Math-PlanePath.
4             #
5             # Math-PlanePath 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-PlanePath 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-PlanePath. If not, see .
17              
18              
19             package Math::PlanePath::CfracDigits;
20 1     1   1208 use 5.004;
  1         4  
21 1     1   6 use strict;
  1         2  
  1         25  
22              
23 1     1   5 use vars '$VERSION', '@ISA';
  1         2  
  1         69  
24             $VERSION = 128;
25 1     1   737 use Math::PlanePath;
  1         3  
  1         44  
26             @ISA = ('Math::PlanePath');
27              
28             use Math::PlanePath::Base::Generic
29 1         64 'is_infinite',
30 1     1   7 'round_nearest';
  1         2  
31             use Math::PlanePath::Base::Digits
32 1         73 'round_down_pow',
33             'digit_split_lowtohigh',
34 1     1   531 'digit_join_lowtohigh';
  1         4  
35              
36 1     1   683 use Math::PlanePath::RationalsTree;
  1         2  
  1         42  
37             *_xy_to_quotients = \&Math::PlanePath::RationalsTree::_xy_to_quotients;
38              
39 1     1   7 use Math::PlanePath::CoprimeColumns;
  1         2  
  1         46  
40             *_coprime = \&Math::PlanePath::CoprimeColumns::_coprime;
41              
42             # uncomment this to run the ### lines
43             #use Smart::Comments;
44              
45              
46 1         54 use constant parameter_info_array =>
47             [ { name => 'radix',
48             share_key => 'radix2_min1',
49             display => 'Radix',
50             type => 'integer',
51             minimum => 1,
52             default => 2,
53             width => 3,
54             },
55 1     1   6 ];
  1         3  
56              
57 1     1   6 use constant n_start => 0;
  1         2  
  1         40  
58 1     1   5 use constant class_x_negative => 0;
  1         2  
  1         56  
59 1     1   6 use constant class_y_negative => 0;
  1         2  
  1         55  
60 1     1   7 use constant x_minimum => 1;
  1         2  
  1         51  
61 1     1   6 use constant y_minimum => 2;
  1         2  
  1         57  
62 1     1   7 use constant diffxy_maximum => -1; # upper octant X<=Y-1 so X-Y<=-1
  1         2  
  1         61  
63 1     1   6 use constant gcdxy_maximum => 1; # no common factor
  1         3  
  1         1280  
64              
65             # FIXME: believe this is right, but check N+1 always changes Y
66             sub absdy_minimum {
67 0     0 1 0 my ($self) = @_;
68 0 0       0 return ($self->{'radix'} < 3 ? 0 : 1);
69             }
70              
71             # radix=1 N=1 has dir4=0
72             # radix=2 N=5628 has dir4=0 dx=9,dy=0
73             # radix=3 N=1189140 has dir4=0 dx=1,dy=0
74             # radix=4 N=169405 has dir4=0 dx=2,dy=0
75             # always eventually 0 ?
76             # use constant dir_minimum_dxdy => (1,0); # the default
77              
78             # radix=1 N=4 dX=1,dY=-1 for dir4=3.5
79             # radix=2 N=4413 dX=9,dY=-1
80             # radix=3 N=9492 dX=3,dY=-1
81             # ENHANCE-ME: suspect believe approaches 360 degrees, eventually, but proof?
82             # use constant dir_maximum_dxdy => (0,0); # the default
83              
84             sub turn_any_straight {
85 0     0 1 0 my ($self) = @_;
86 0         0 return ($self->{'radix'} != 1); # radix=1 never straight
87             }
88              
89             sub _UNDOCUMENTED__turn_any_left_at_n {
90 0     0   0 my ($self) = @_;
91 0         0 return $self->{'radix'} + 1;
92             }
93             sub _UNDOCUMENTED__turn_any_right_at_n {
94 0     0   0 my ($self) = @_;
95 0         0 return $self->{'radix'};
96             }
97              
98              
99             #------------------------------------------------------------------------------
100              
101             sub new {
102 3     3 1 604 my $self = shift->SUPER::new (@_);
103 3 50 33     19 unless ($self->{'radix'} && $self->{'radix'} >= 1) {
104 3         6 $self->{'radix'} = 2;
105             }
106 3         8 return $self;
107             }
108              
109             sub n_to_xy {
110 1     1 1 10 my ($self, $n) = @_;
111             ### CfracDigits n_to_xy(): "$n"
112              
113 1 50       4 if ($n < 0) { return; }
  0         0  
114 1 50       5 if (is_infinite($n)) { return ($n,$n); }
  0         0  
115              
116             {
117 1         3 my $int = int($n);
  1         2  
118 1 50       5 if ($n != $int) {
119             ### frac ...
120 0         0 my $frac = $n - $int; # inherit possible BigFloat/BigRat
121 0         0 my ($x1,$y1) = $self->n_to_xy($int);
122 0         0 my ($x2,$y2) = $self->n_to_xy($int+1);
123 0         0 my $dx = $x2-$x1;
124 0         0 my $dy = $y2-$y1;
125             ### x1,y1: "$x1, $y1"
126             ### x2,y2: "$x2, $y2"
127             ### dx,dy: "$dx, $dy"
128             ### result: ($frac*$dx + $x1).', '.($frac*$dy + $y1)
129 0         0 return ($frac*$dx + $x1, $frac*$dy + $y1);
130             }
131 1         2 $n = $int;
132             }
133              
134 1         3 my $radix = $self->{'radix'};
135 1         2 my $zero = ($n * 0); # inherit bignum 0
136 1         2 my $x = $zero;
137 1         2 my $y = 1 + $zero; # inherit bignum 1
138              
139 1         3 foreach my $q (_n_to_quotients_bottomtotop($n,$radix,$zero)) { # bottom to top
140             ### at: "$x,$y q=$q"
141              
142             # 1/(q + X/Y) = 1/((qY+X)/Y)
143             # = Y/(qY+X)
144 4         8 ($x,$y) = ($y, $q*$y + $x);
145             }
146              
147             ### return: "$x,$y"
148 1         3 return ($x,$y);
149             }
150              
151             # Return a list of quotients bottom to top. The base3 digits of N are split
152             # by "3" delimiters and the parts adjusted so the first bottom-most q>=2 and
153             # the rest q>=1. The values are ready to be used as continued fraction
154             # terms.
155             #
156             sub _n_to_quotients_bottomtotop {
157 6     6   336 my ($n, $radix, $zero) = @_;
158             ### _n_to_quotients_bottomtotop(): $n
159              
160 6         10 my $radix_plus_1 = $radix + 1;
161 6         10 my @ret;
162             my @group;
163 6         14 foreach my $digit (_digit_split_1toR_lowtohigh($n,$radix_plus_1)) {
164 21 100       31 if ($digit == $radix_plus_1) {
165             ### @group
166 7         15 push @ret, _digit_join_1toR_destructive(\@group, $radix, $zero) + 1;
167 7         15 @group = ();
168             } else {
169 14         25 push @group, $digit;
170             }
171             }
172             ### final group: @group
173 6         13 push @ret, _digit_join_1toR_destructive(\@group, $radix, $zero) + 1;
174              
175 6         13 $ret[0] += 1; # bottom-most is +2 rather than +1
176              
177             ### _n_to_quotients_bottomtotop result: @ret
178 6         13 return @ret;
179             }
180              
181             # Return a list of digits 1 <= d <= R which is $n written in $radix, low to
182             # high digits.
183             sub _digit_split_1toR_lowtohigh {
184 28     28   704 my ($n, $radix) = @_;
185             ### assert: $radix >= 1
186             ### assert: $n >= 0
187              
188 28 50       57 if ($radix == 1) {
189 0         0 return (1) x $n;
190             }
191 28         65 my @digits = digit_split_lowtohigh($n,$radix);
192              
193             # mangle 0 -> R
194 28         43 my $borrow = 0;
195 28         47 foreach my $digit (@digits) { # low to high
196 61 100       123 if ($borrow = (($digit -= $borrow) <= 0)) { # modify array contents
197 26         37 $digit += $radix;
198             }
199             }
200 28 100       49 if ($borrow) {
201             ### assert: $digits[-1] == $radix
202 6         8 pop @digits;
203             }
204              
205 28         74 return @digits;
206             }
207              
208             sub _digit_join_1toR_destructive {
209 18     18   31 my ($aref, $radix, $zero) = @_;
210             ### assert: $radix >= 1
211              
212 18 50       112 if ($radix == 1) {
213 0         0 return scalar(@$aref);
214             }
215              
216             # mangle any digit==$radix down to digit=0
217 18         26 my $carry = 0;
218 18         30 foreach my $digit (@$aref) { # low to high
219 34 100       68 if ($carry = (($digit += $carry) >= $radix)) { # modify array contents
220 16         25 $digit -= $radix;
221             }
222             }
223 18 100       30 if ($carry) {
224 6         8 push @$aref, 1;
225             }
226              
227             ### _digit_join_1toR_destructive() result: digit_join_lowtohigh($aref, $radix, $zero)
228 18         38 return digit_join_lowtohigh($aref, $radix, $zero);
229             }
230              
231             sub xy_is_visited {
232 0     0 1 0 my ($self, $x, $y) = @_;
233 0         0 $x = round_nearest ($x);
234 0         0 $y = round_nearest ($y);
235 0   0     0 return (! ($x < 1 || $y < 2 || $x >= $y)
236             && _coprime($x,$y));
237             }
238              
239             sub xy_to_n {
240 1     1 1 118 my ($self, $x, $y) = @_;
241 1         4 $x = round_nearest ($x);
242 1         2 $y = round_nearest ($y);
243             ### CfracDigits xy_to_n(): "$x,$y"
244              
245 1 50       3 if (is_infinite($x)) { return $x; }
  0         0  
246 1 50       3 if (is_infinite($y)) { return $y; }
  0         0  
247 1 50 33     11 if ($x < 1 || $y < 2 || $x >= $y) {
      33        
248 0         0 return undef;
249             }
250              
251 1 50       4 my @quotients = _xy_to_quotients($x,$y)
252             or return undef; # $x,$y have a common factor
253             ### @quotients
254              
255             # drop initial 0 integer part
256             ### assert: $quotients[0] == 0
257 1         3 shift @quotients;
258              
259             return _cfrac_join_toptobottom(\@quotients,
260 1         4 $self->{'radix'},
261             $x * 0 * $y); # inherit bignum 0
262             }
263              
264             # $aref is a list of continued fraction quotients from top-most to
265             # bottom-most. There's no initial integer term in $aref. Each quotient is
266             # q >= 1 except the bottom-most which q-1 and so also >=1.
267             #
268             sub _cfrac_join_toptobottom {
269 5     5   306 my ($aref, $radix, $zero) = @_;
270             ### _cfrac_join_toptobottom(): $aref
271              
272 5         8 my @digits;
273 5         9 foreach my $q (reverse @$aref) {
274             ### assert: $q >= 1
275 12         26 push @digits, _digit_split_1toR_lowtohigh($q-1, $radix), $radix+1;
276             }
277 5         14 pop @digits; # no high delimiter
278             ### groups digits 1toR: @digits
279 5         12 return _digit_join_1toR_destructive(\@digits, $radix+1, $zero);
280             }
281              
282              
283             # X/Y = F[k]/F[k+1] quotients all 1
284             # N = all delimiter digits R,R,...,R
285             # = 1222...2221
286             # = R^k + 2*(R^k+1)/(R-1) - 1
287             # = (RR^k - R^k + 2R^k + 2 - R + 1) / (R-1)
288             # = (RR^k + R^k - R + 3) / (R-1)
289             # = ((R+1)R^k - R + 3) / (R-1)
290             # take high as "12" = R+2
291             # k = log(Y)/log(phi)
292             # N = (R+2) * R ** k
293             # N = Y ** (log(R)/log(phi))
294             #
295             # not exact
296             sub rect_to_n_range {
297 2     2 1 222 my ($self, $x1,$y1, $x2,$y2) = @_;
298             ### rect_to_n_range()
299              
300 2         13 $x1 = round_nearest ($x1);
301 2         5 $y1 = round_nearest ($y1);
302 2         4 $x2 = round_nearest ($x2);
303 2         5 $y2 = round_nearest ($y2);
304              
305 2 50       6 ($x1,$x2) = ($x2,$x1) if $x1 > $x2;
306 2 50       5 ($y1,$y2) = ($y2,$y1) if $y1 > $y2;
307             ### $x2
308             ### $y2
309              
310             # | /
311             # | / x1
312             # | / +-----y2
313             # | / |
314             # |/ +-----
315             #
316 2 50 33     13 if ($x2 < 1 || $y2 < 2 || $x1 >= $y2) {
      33        
317             ### no values, rect outside upper octant ...
318 0         0 return (1,0);
319             }
320              
321 2         5 my $zero = ($x1 * 0 * $y1 * $x2 * $y2); # inherit bignum
322 2         3 my $radix = $self->{'radix'};
323              
324 2 50       7 return (0,
325             ($radix+3)
326             * ($radix+1 + $zero) ** ($radix == 1
327             ? $y2
328             : _log_phi_estimate($y2,$radix)));
329             }
330              
331             # Return an estimate of log base phi of $x, that being log($x)/log(phi),
332             # where phi=(1+sqrt(5))/2 the golden ratio.
333             #
334             sub _log_phi_estimate {
335 2     2   4 my ($x) = @_;
336 2         5 my ($pow,$exp) = round_down_pow ($x, 2);
337 2         10 return int ($exp * (log(2) / log((1+sqrt(5))/2)));
338             }
339              
340             1;
341             __END__