File Coverage

blib/lib/Math/PlanePath/WythoffArray.pm
Criterion Covered Total %
statement 43 122 35.2
branch 0 22 0.0
condition 4 16 25.0
subroutine 15 19 78.9
pod 7 7 100.0
total 69 186 37.1


line stmt bran cond sub pod time code
1             # Copyright 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021 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             # Classic Sequences
20             # http://oeis.org/classic.html
21             #
22             # Clark Kimberling
23             # http://faculty.evansville.edu/ck6/integer/intersp.html
24             #
25             # cf A175004 similar to wythoff but rows recurrence
26             # r(n-1)+r(n-2)+1 extra +1 in each step
27             # floor(n*phi+2/phi)
28             #
29             # cf Stolarsky round_nearest(n*phi)
30             # A035506 stolarsky by diagonals
31             # A035507 inverse
32             # A007067 stolarsky first column
33              
34             # Maybe:
35             # my ($x,$y) = $path->pair_to_xy($a,$b);
36             # Return the $x,$y which has ($a,$b).
37             # Advance $a,$b if before start of row.
38              
39             # Carlitz and Hoggatt "Fibonacci Representations", Fibonacci Quarterly,
40             # volume 10, number 1, January 1972
41             # http://www.fq.math.ca/10-1.html
42             # http://www.fq.math.ca/Scanned/10-1/carlitz1.pdf
43              
44              
45             package Math::PlanePath::WythoffArray;
46 1     1   1071 use 5.004;
  1         3  
47 1     1   5 use strict;
  1         2  
  1         40  
48             #use List::Util 'max';
49             *max = \&Math::PlanePath::_max;
50              
51 1     1   5 use vars '$VERSION', '@ISA';
  1         1  
  1         56  
52             $VERSION = 129;
53 1     1   717 use Math::PlanePath;
  1         2  
  1         51  
54             *_sqrtint = \&Math::PlanePath::_sqrtint;
55             @ISA = ('Math::PlanePath');
56              
57             use Math::PlanePath::Base::Generic
58 1         44 'is_infinite',
59 1     1   13 'round_nearest';
  1         2  
60             use Math::PlanePath::Base::Digits
61 1     1   542 'bit_split_lowtohigh';
  1         2  
  1         95  
62              
63             # uncomment this to run the ### lines
64             #use Smart::Comments;
65              
66              
67 1         57 use constant parameter_info_array =>
68             [ { name => 'x_start',
69             display => 'X start',
70             type => 'integer',
71             default => 0,
72             width => 3,
73             description => 'Starting X coordinate.',
74             },
75             { name => 'y_start',
76             display => 'Y start',
77             type => 'integer',
78             default => 0,
79             width => 3,
80             description => 'Starting Y coordinate.',
81             },
82 1     1   7 ];
  1         2  
83              
84 1     1   5 use constant default_n_start => 1;
  1         2  
  1         41  
85 1     1   5 use constant class_x_negative => 0;
  1         2  
  1         42  
86 1     1   6 use constant class_y_negative => 0;
  1         2  
  1         104  
87              
88             sub x_minimum {
89 1     1 1 10 my ($self) = @_;
90 1         5 return $self->{'x_start'};
91             }
92             sub y_minimum {
93 1     1 1 2 my ($self) = @_;
94 1         3 return $self->{'y_start'};
95             }
96 1     1   7 use constant absdx_minimum => 1;
  1         13  
  1         68  
97 1     1   6 use constant dir_maximum_dxdy => (3,-1); # N=4 to N=5 dX=3,dY=-1
  1         2  
  1         917  
98              
99             #------------------------------------------------------------------------------
100              
101             sub new {
102 3     3 1 689 my $self = shift->SUPER::new(@_);
103 3   100     24 $self->{'x_start'} ||= 0;
104 3   100     11 $self->{'y_start'} ||= 0;
105 3         7 return $self;
106             }
107              
108             sub xy_is_visited {
109 0     0 1   my ($self, $x, $y) = @_;
110             return ((round_nearest($x) >= $self->{'x_start'})
111 0   0       && (round_nearest($y) >= $self->{'y_start'}));
112             }
113              
114             #------------------------------------------------------------------------------
115             # 4 | 12 20 32 52 84 136 220 356 576 932 1508
116             # 3 | 9 15 24 39 63 102 165 267 432 699 1131
117             # 2 | 6 10 16 26 42 68 110 178 288 466 754
118             # 1 | 4 7 11 18 29 47 76 123 199 322 521
119             # Y=0 | 1 2 3 5 8 13 21 34 55 89 144
120             # +-------------------------------------------------------
121             # X=0 1 2 3 4 5 6 7 8 9 10
122             # 13,8,5,3,2,1
123             # 4 = 3+1 -> 1
124             # 6 = 5+1 -> 2
125             # 9 = 8+1 -> 3
126             # 12 = 8+3+1 -> 3+1=4
127             # 14 = 13+1 -> 5
128              
129             sub n_to_xy {
130 0     0 1   my ($self, $n) = @_;
131             ### WythoffArray n_to_xy(): $n
132              
133 0 0         if ($n < 1) { return; }
  0            
134 0 0 0       if (is_infinite($n) || $n == 0) { return ($n,$n); }
  0            
135              
136             {
137             # fractions on straight line between integer points
138 0           my $int = int($n);
  0            
139 0 0         if ($n != $int) {
140 0           my $frac = $n - $int; # inherit possible BigFloat/BigRat
141 0           my ($x1,$y1) = $self->n_to_xy($int);
142 0           my ($x2,$y2) = $self->n_to_xy($int+1);
143 0           my $dx = $x2-$x1;
144 0           my $dy = $y2-$y1;
145 0           return ($frac*$dx + $x1, $frac*$dy + $y1);
146             }
147 0           $n = $int;
148             }
149              
150             # f1+f0 > i
151             # f0 > i-f1
152             # check i-f1 as the stopping point, so that if i=UV_MAX then won't
153             # overflow a UV trying to get to f1>=i
154             #
155 0           my @fibs;
156             {
157 0           my $f0 = ($n * 0); # inherit bignum 0
  0            
158 0           my $f1 = $f0 + 1; # inherit bignum 1
159 0           while ($f0 <= $n-$f1) {
160 0           ($f1,$f0) = ($f1+$f0,$f1);
161 0           push @fibs, $f1; # starting $fibs[0]=1
162             }
163             }
164             ### @fibs
165              
166             # indices into fib[] which are the Fibonaccis adding up to $n
167 0           my @indices;
168 0           for (my $i = $#fibs; $i >= 0; $i--) {
169             ### at: "n=$n f=".$fibs[$i]
170 0 0         if ($n >= $fibs[$i]) {
171 0           push @indices, $i;
172 0           $n -= $fibs[$i];
173             ### sub: "$fibs[$i] to n=$n"
174 0           --$i;
175             }
176             }
177             ### @indices
178              
179             # X is low index, ie. how many low 0 bits in Zeckendorf form
180 0           my $x = pop @indices;
181             ### $x
182              
183             # Y is indices shifted down by $x and 2 more
184 0           my $y = 0;
185 0           my $shift = $x+2;
186 0           foreach my $i (@indices) {
187             ### y add: "ishift=".($i-$shift)." fib=".$fibs[$i-$shift]
188 0           $y += $fibs[$i-$shift];
189             }
190             ### $shift
191             ### $y
192              
193 0           return ($x+$self->{'x_start'},$y+$self->{'y_start'});
194             }
195              
196             # phi = (sqrt(5)+1)/2
197             # (y+1)*phi = (y+1)*(sqrt(5)+1)/2
198             # = ((y+1)*sqrt(5)+(y+1))/2
199             # = (sqrt(5*(y+1)^2)+(y+1))/2
200             #
201             # from x=0,y=0
202             # N = floor((y+1)*Phi) * Fib(x+2) + y*Fib(x+1)
203             #
204             sub xy_to_n {
205 0     0 1   my ($self, $x, $y) = @_;
206             ### WythoffArray xy_to_n(): "$x, $y"
207              
208 0           $x = round_nearest($x) - $self->{'x_start'};
209 0           $y = round_nearest($y) - $self->{'y_start'};
210 0 0 0       if ($x < 0 || $y < 0) {
211 0           return undef;
212             }
213              
214 0           my $zero = $x * 0 * $y;
215 0           $x += 2;
216 0 0         if (is_infinite($x)) { return $x; }
  0            
217 0 0         if (is_infinite($y)) { return $y; }
  0            
218              
219 0           my @bits = bit_split_lowtohigh($x);
220             ### @bits
221 0           pop @bits; # discard high 1-bit
222              
223 0           my $yplus1 = $zero + $y+1; # inherit bigint from $x perhaps
224              
225             # spectrum(Y+1) so Y,Ybefore are notional two values at X=-2 and X=-1
226 0           my $ybefore = int((_sqrtint(5*$yplus1*$yplus1) + $yplus1) / 2);
227             ### $ybefore
228              
229             # k=1, Fk1=F[k-1]=0, Fk=F[k]=1
230 0           my $Fk1 = $zero;
231 0           my $Fk = 1 + $zero;
232              
233 0           my $add = -2;
234 0           while (@bits) {
235             ### remaining bits: @bits
236             ### Fk1: "$Fk1"
237             ### Fk: "$Fk"
238              
239             # two squares and some adds
240             # F[2k+1] = 4*F[k]^2 - F[k-1]^2 + 2*(-1)^k
241             # F[2k-1] = F[k]^2 + F[k-1]^2
242             # F[2k] = F[2k+1] - F[2k-1]
243             #
244 0           $Fk *= $Fk;
245 0           $Fk1 *= $Fk1;
246 0           my $F2kplus1 = 4*$Fk - $Fk1 + $add;
247 0           $Fk1 += $Fk; # F[2k-1]
248 0           my $F2k = $F2kplus1 - $Fk1;
249              
250 0 0         if (pop @bits) { # high to low
251 0           $Fk1 = $F2k; # F[2k]
252 0           $Fk = $F2kplus1; # F[2k+1]
253 0           $add = -2;
254             } else {
255             # $Fk1 is F[2k-1] already
256 0           $Fk = $F2k; # F[2k]
257 0           $add = 2;
258             }
259             }
260              
261             ### final pair ...
262             ### Fk1: "$Fk1"
263             ### Fk: "$Fk"
264             ### @bits
265              
266 0           return ($Fk*$ybefore + $Fk1*$y);
267             }
268              
269             # exact
270             sub rect_to_n_range {
271 0     0 1   my ($self, $x1,$y1, $x2,$y2) = @_;
272             ### WythoffArray rect_to_n_range(): "$x1,$y1 $x2,$y2"
273              
274 0           $x1 = round_nearest($x1);
275 0           $y1 = round_nearest($y1);
276 0           $x2 = round_nearest($x2);
277 0           $y2 = round_nearest($y2);
278              
279 0 0         ($x1,$x2) = ($x2,$x1) if $x1 > $x2;
280 0 0         ($y1,$y2) = ($y2,$y1) if $y1 > $y2;
281              
282 0 0 0       if ($x2 < $self->{'x_start'} || $y2 < $self->{'y_start'}) {
283             ### all outside first quadrant ...
284 0           return (1, 0);
285             }
286              
287             # bottom left into first quadrant
288 0           $x1 = max($x1, $self->{'x_start'});
289 0           $y1 = max($y1, $self->{'y_start'});
290              
291 0           return ($self->xy_to_n($x1,$y1), # bottom left
292             $self->xy_to_n($x2,$y2)); # top right
293             }
294              
295             1;
296             __END__