File Coverage

blib/lib/Math/PlanePath/CellularRule54.pm
Criterion Covered Total %
statement 102 119 85.7
branch 25 36 69.4
condition 6 9 66.6
subroutine 21 23 91.3
pod 5 5 100.0
total 159 192 82.8


line stmt bran cond sub pod time code
1             # Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018 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             # math-image --path=CellularRule54 --all --scale=10
19             # math-image --path=CellularRule54 --all --output=numbers --size=132x50
20              
21             package Math::PlanePath::CellularRule54;
22 2     2   34 use 5.004;
  2         7  
23 2     2   11 use strict;
  2         4  
  2         48  
24              
25 2     2   8 use vars '$VERSION', '@ISA';
  2         4  
  2         111  
26             $VERSION = 127;
27 2     2   11 use Math::PlanePath;
  2         3  
  2         131  
28             @ISA = ('Math::PlanePath');
29             *_divrem = \&Math::PlanePath::_divrem;
30             *_sqrtint = \&Math::PlanePath::_sqrtint;
31              
32             use Math::PlanePath::Base::Generic
33 2     2   13 'round_nearest';
  2         3  
  2         123  
34              
35             # uncomment this to run the ### lines
36             #use Smart::Comments;
37              
38              
39 2         107 use constant parameter_info_array =>
40             [ Math::PlanePath::Base::Generic::parameter_info_nstart1(),
41 2     2   12 ];
  2         41  
42              
43 2     2   11 use constant class_y_negative => 0;
  2         11  
  2         90  
44 2     2   10 use constant n_frac_discontinuity => .5;
  2         5  
  2         161  
45             sub x_negative_at_n {
46 0     0 1 0 my ($self) = @_;
47 0         0 return $self->n_start + 1;
48             }
49 2     2   12 use constant sumxy_minimum => 0; # triangular X>=-Y so X+Y>=0
  2         4  
  2         95  
50 2     2   12 use constant diffxy_maximum => 0; # triangular X<=Y so X-Y<=0
  2         3  
  2         94  
51 2     2   11 use constant dx_maximum => 4;
  2         4  
  2         111  
52 2     2   12 use constant dy_minimum => 0;
  2         4  
  2         116  
53 2     2   24 use constant dy_maximum => 1;
  2         5  
  2         92  
54 2     2   12 use constant absdx_minimum => 1;
  2         4  
  2         107  
55 2     2   12 use constant dsumxy_maximum => 4; # straight East dX=+4
  2         4  
  2         79  
56 2     2   11 use constant ddiffxy_maximum => 4; # straight East dX=+4
  2         3  
  2         104  
57 2     2   12 use constant dir_maximum_dxdy => (-1,0); # supremum, West and dY=+1 up
  2         3  
  2         1517  
58              
59              
60             #------------------------------------------------------------------------------
61              
62             sub new {
63 1     1 1 9 my $self = shift->SUPER::new (@_);
64 1 50       9 if (! defined $self->{'n_start'}) {
65 0         0 $self->{'n_start'} = $self->default_n_start;
66             }
67 1         6 return $self;
68             }
69             # left add
70             # even y=0 0 1
71             # 2 1 2
72             # 4 3 3
73             # 6 6 4
74             # left = y/2*(y/2+1)/2
75             # = y*(y+2)/8 of 4-cell figures
76             # inverse y = -1 + sqrt(2 * $n + -1)
77             #
78             # left add
79             # odd y=1 0 3
80             # 3 3 6
81             # 5 9 9
82             # 7 18 12
83             # left = 3*(y-1)/2*((y-1)/2+1)/2
84             # = 3*(y-1)*(y+1)/8 of 4-cell figures
85             #
86             # nbase y even = y*(y+2)/8 + 3*((y+1)-1)*((y+1)+1)/8
87             # = [ y*(y+2) + 3*y*(y+2) ] / 8
88             # = y*(y+2)/2
89             # y=0 nbase=0
90             # y=2 nbase=4
91             # y=4 nbase=12
92             # y=6 nbase=24
93             #
94             # nbase y odd = 3*(y-1)*(y+1)/8 + (y+1)*(y+3)/8
95             # = (y+1) * (3y-3 + y+3)/8
96             # = (y+1)*4y/8
97             # = y*(y+1)/2
98             # y=1 nbase=1
99             # y=3 nbase=6
100             # y=5 nbase=15
101             # y=7 nbase=28
102             # inverse y = -1/2 + sqrt(2 * $n + -7/4)
103             # = sqrt(2n-7/4) - 1/2
104             # = (2*sqrt(2n-7/4) - 1)/2
105             # = (sqrt(4n-7)-1)/2
106             #
107             # dual
108             # d = [ 0, 1, 2, 3 ]
109             # N = [ 1, 5, 13, 25 ]
110             # N = (2 d^2 + 2 d + 1)
111             # = ((2*$d + 2)*$d + 1)
112             # d = -1/2 + sqrt(1/2 * $n + -1/4)
113             # = sqrt(1/2 * $n + -1/4) - 1/2
114             # = [ 2*sqrt(1/2 * $n + -1/4) - 1 ] / 2
115             # = [ sqrt(4/2 * $n + -4/4) - 1 ] / 2
116             # = [ sqrt(2*$n - 1) - 1 ] / 2
117             #
118              
119             sub n_to_xy {
120 25     25 1 188 my ($self, $n) = @_;
121             ### CellularRule54 n_to_xy(): $n
122              
123 25         41 $n = $n - $self->{'n_start'}; # to N=0 basis
124 25         33 my $frac;
125             {
126 25         40 my $int = int($n);
  25         56  
127 25         35 $frac = $n - $int;
128 25         34 $n = $int; # BigFloat int() gives BigInt, use that
129 25 50       55 if (2*$frac >= 1) { # $frac>=0.5 and BigInt friendly
130 0         0 $frac -= 1;
131 0         0 $n += 1;
132             }
133             # -0.5 <= $frac < 0.5
134             ### assert: $frac >= -0.5
135             ### assert: $frac < 0.5
136             }
137              
138 25 100       46 if ($n < 0) {
139 3         9 return;
140             }
141              
142             # d is the two-row group number, d=2*y, where n belongs
143             # start of the two-row group is nbase = 2 d^2 + 2 d starting from N=0
144             #
145 22         50 my $d = int((_sqrtint(2*$n+1) - 1) / 2);
146 22         42 $n -= (2*$d + 2)*$d; # remainder within two-row
147             ### $d
148             ### remainder: $n
149 22 100       44 if ($n <= $d) {
150             # d+1 many points in the Y=0,2,4,6 etc even row, spaced 4*n apart
151 7         10 $d *= 2; # y=2*d
152 7         20 return ($frac + 4*$n - $d,
153             $d);
154             } else {
155             # 3*d many points in the Y=1,3,5,7 etc odd row, using 3 in 4 cells
156 15         23 $n -= $d+1; # remainder 0 upwards into odd row
157 15         20 $d = 2*$d+1; # y=2*d+1
158 15         33 my ($q) = _divrem($n,3);
159 15         37 return ($frac + $n + $q - $d,
160             $d);
161             }
162             }
163              
164             sub xy_to_n {
165 496     496 1 2430 my ($self, $x, $y) = @_;
166 496         884 $x = round_nearest ($x);
167 496         923 $y = round_nearest ($y);
168             ### CellularRule54 xy_to_n(): "$x,$y"
169              
170 496 100 66     1869 if ($y < 0
      100        
171             || $x < -$y
172             || $x > $y) {
173 240         439 return undef;
174             }
175 256         380 $x += $y;
176             ### x centred: $x
177 256 100       432 if ($y % 2) {
178             ### odd row, 3 in 4 ...
179 136 100       239 if (($x % 4) == 3) {
180 28         56 return undef;
181             }
182 108         332 return $x - int($x/4) + $y*($y+1)/2 + $self->{'n_start'};
183             } else {
184             ## even row, sparse ...
185 120 100       204 if ($x % 4) {
186 84         167 return undef;
187             }
188 36         101 return $x/4 + $y*($y+2)/2 + $self->{'n_start'};
189             }
190             }
191              
192             # not exact
193             sub rect_to_n_range {
194 0     0 1 0 my ($self, $x1,$y1, $x2,$y2) = @_;
195             ### CellularRule54 rect_to_n_range(): "$x1,$y1, $x2,$y2"
196              
197 0 0       0 ($x1,$y1, $x2,$y2) = _rect_for_V ($x1,$y1, $x2,$y2)
198             or return (1,0); # rect outside pyramid
199              
200 0         0 my $zero = ($x1 * 0 * $y1 * $x2 * $y2); # inherit bignum
201              
202             # nbase y even y*(y+2)/2
203             # nbase y odd y*(y+1)/2
204             # y even end (y+1)*(y+2)/2
205             # y odd end (y+1)*(y+3)/2
206              
207 0         0 $y2 += 1;
208             return (# even/odd left end
209             $zero + $y1*($y1 + 2-($y1%2))/2 + $self->{'n_start'},
210              
211             # even/odd right end
212 0         0 $zero + $y2*($y2 + 2-($y2%2))/2 + $self->{'n_start'} - 1);
213             }
214              
215             # Return ($x1,$y1, $x2,$y2) which is the rectangle part chopped to the top
216             # row entirely within the pyramid V and the bottom row partly within.
217             #
218             sub _rect_for_V {
219 222     222   427 my ($x1,$y1, $x2,$y2) = @_;
220             ### _rect_for_V(): "$x1,$y1, $x2,$y2"
221              
222 222         499 $y1 = round_nearest ($y1);
223 222         457 $y2 = round_nearest ($y2);
224 222 100       448 if ($y1 > $y2) { ($y1,$y2) = ($y2,$y1); } # swap to y1<=y2
  50         94  
225              
226 222 50       475 unless ($y2 >= 0) {
227             ### rect all negative, no N ...
228 0         0 return;
229             }
230 222 50       390 unless ($y1 >= 0) {
231             # increase y1 to zero, including negative infinity discarded
232 0         0 $y1 = 0;
233             }
234              
235 222         398 $x1 = round_nearest ($x1);
236 222         410 $x2 = round_nearest ($x2);
237 222 100       432 if ($x1 > $x2) { ($x1,$x2) = ($x2,$x1); } # swap to x1<=x2
  111         184  
238 222         323 my $neg_y2 = -$y2;
239              
240             # \ /
241             # y2 \ / +-----
242             # \ / |
243             # \ /
244             # \/ x1
245             #
246             # \ /
247             # ----+ \ / y2
248             # | \ /
249             # \ /
250             # x2 \/
251             #
252 222 50 33     687 if ($x1 > $y2 # off to the right
253             || $x2 < $neg_y2) { # off to the left
254             ### rect all off to the left or right, no N
255 0         0 return;
256             }
257              
258             # \ / x2
259             # \ +------+ y2
260             # \ | / |
261             # \ +------+
262             # \/
263             #
264 222 50       419 if ($x2 > $y2) {
265             ### top-right beyond pyramid, reduce ...
266 0         0 $x2 = $y2;
267             }
268              
269             #
270             # x1 \ /
271             # y2 +--------+ / y2
272             # | \ | /
273             # +--------+/
274             # \/
275             #
276 222 50       405 if ($x1 < $neg_y2) {
277             ### top-left beyond pyramid, increase ...
278 0         0 $x1 = $neg_y2;
279             }
280              
281             # \ | /
282             # \ |/
283             # \ /| |
284             # y1 \ / +-------+
285             # \/ x1
286             #
287             # \| /
288             # \ /
289             # |\ /
290             # -------+ \ / y1
291             # x2 \/
292             #
293             # in both of the following y1=x2 or y1=-x2 leaves y1<=y2 because have
294             # already established some part of the rectangle is in the V shape
295             #
296 222 50       480 if ($x1 > $y1) {
    50          
297             ### x1 off to the right, so y1 row is outside, increase y1 ...
298 0         0 $y1 = $x1;
299              
300             } elsif ((my $neg_x2 = -$x2) > $y1) {
301             ### x2 off to the left, so y1 row is outside, increase y1 ...
302 0         0 $y1 = $neg_x2;
303             }
304              
305             # values ordered
306             ### assert: $x1 <= $x2
307             ### assert: $y1 <= $y2
308              
309             # top row x1..x2 entirely within pyramid
310             ### assert: $x1 >= -$y2
311             ### assert: $x2 <= $y2
312              
313             # bottom row x1..x2 some part within pyramid
314             ### assert: $x1 <= $y1
315             ### assert: $x2 >= -$y1
316              
317 222         732 return ($x1,$y1, $x2,$y2);
318             }
319              
320             1;
321             __END__