File Coverage

blib/lib/Math/PlanePath/CellularRule190.pm
Criterion Covered Total %
statement 101 110 91.8
branch 34 42 80.9
condition 2 3 66.6
subroutine 22 23 95.6
pod 5 5 100.0
total 164 183 89.6


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              
19             # http://mathworld.wolfram.com/ElementaryCellularAutomaton.html
20             #
21             # Loeschian numbers strips on the right ...
22              
23              
24             package Math::PlanePath::CellularRule190;
25 2     2   9262 use 5.004;
  2         13  
26 2     2   11 use strict;
  2         6  
  2         74  
27              
28 2     2   11 use vars '$VERSION', '@ISA';
  2         5  
  2         143  
29             $VERSION = 127;
30 2     2   690 use Math::PlanePath;
  2         4  
  2         106  
31             *_sqrtint = \&Math::PlanePath::_sqrtint;
32             @ISA = ('Math::PlanePath');
33              
34             use Math::PlanePath::Base::Generic
35 2     2   15 'round_nearest';
  2         3  
  2         98  
36              
37 2     2   454 use Math::PlanePath::CellularRule54;
  2         4  
  2         88  
38             *_rect_for_V = \&Math::PlanePath::CellularRule54::_rect_for_V;
39              
40             # uncomment this to run the ### lines
41             # use Smart::Comments;
42              
43              
44 2     2   13 use constant class_y_negative => 0;
  2         4  
  2         134  
45 2     2   12 use constant n_frac_discontinuity => .5;
  2         3  
  2         198  
46              
47 2         202 use constant parameter_info_array =>
48             [ { name => 'mirror',
49             display => 'Mirror',
50             type => 'boolean',
51             default => 0,
52             description => 'Mirror to "rule 246" instead.',
53             },
54             Math::PlanePath::Base::Generic::parameter_info_nstart1(),
55 2     2   13 ];
  2         3  
56              
57             sub x_negative_at_n {
58 0     0 1 0 my ($self) = @_;
59 0         0 return $self->n_start + 1;
60             }
61 2     2   16 use constant sumxy_minimum => 0; # triangular X>=-Y so X+Y>=0
  2         4  
  2         100  
62 2     2   13 use constant diffxy_maximum => 0; # triangular X<=Y so X-Y<=0
  2         4  
  2         102  
63 2     2   13 use constant dx_maximum => 2; # across gap
  2         4  
  2         111  
64 2     2   13 use constant dy_minimum => 0;
  2         4  
  2         97  
65 2     2   13 use constant dy_maximum => 1;
  2         3  
  2         97  
66 2     2   13 use constant absdx_minimum => 1;
  2         3  
  2         110  
67 2     2   12 use constant dsumxy_maximum => 2; # straight East dX=+2
  2         4  
  2         95  
68 2     2   11 use constant ddiffxy_maximum => 2; # straight East dX=+2
  2         17  
  2         97  
69 2     2   13 use constant dir_maximum_dxdy => (-1,0); # supremum, West except dY=+1
  2         5  
  2         1343  
70              
71              
72             #------------------------------------------------------------------------------
73              
74             sub new {
75 10     10 1 1832 my $self = shift->SUPER::new (@_);
76 10 100       42 if (! defined $self->{'n_start'}) {
77 7         26 $self->{'n_start'} = $self->default_n_start;
78             }
79 10         31 return $self;
80             }
81              
82             # 31 32 33 34 35 36 37 38 39 40
83             # 22 23 24 25 26 27 28 29 30
84             # 15 6 17 18 19 20 21
85             # 9 10 11 12 13 14
86             # 5 6 7 8
87             # 2 3 4
88             # 1
89             #
90             # even y = [ 0, 2, 4, 6 ]
91             # N = [ 1, 5, 15, 31 ]
92             # Neven = (3/4 y^2 + 1/2 y + 1)
93             # = (3y + 2)*y/4 + 1
94             # = ((3y + 2)*y + 4) /4
95             # = (3 (y/2)^2 + (y/2) + 1)
96             # = (3*(y/2) + 1)*(y/2) + 1
97             #
98             # odd y = [ 1, 3, 5,7 ]
99             # N = [ 2,9,22,41 ]
100             # Nodd = (3/4 y^2 + 1/2 y + 3/4)
101             # = ((3y+2)*y+ 3) / 4
102             #
103             # pair even d = [0,1,2,3]
104             # N = [ 1, 5, 15, 31 ]
105             # Npair = (3 d^2 + d + 1)
106             # d = -1/6 + sqrt(1/3 * $n + -11/36)
107             # = [ -1 + sqrt(1/3 * $n + -11/36)*6 ] / 6
108             # = [ -1 + sqrt(1/3 * $n*36 + -11/36*36) ] / 6
109             # = [ -1 + sqrt(12n-11) ] / 6
110             #
111             sub n_to_xy {
112 140     140 1 11343 my ($self, $n) = @_;
113             ### CellularRule190 n_to_xy(): $n
114              
115 140         254 $n = $n - $self->{'n_start'} + 1; # to N=1 basis, and warn if $n undef
116 140         198 my $frac;
117             {
118 140         198 my $int = int($n);
  140         205  
119 140         200 $frac = $n - $int;
120 140         193 $n = $int; # BigFloat int() gives BigInt, use that
121 140 50       367 if (2*$frac >= 1) {
122 0         0 $frac -= 1;
123 0         0 $n += 1;
124             }
125             # now -0.5 <= $frac < 0.5
126             ### assert: 2*$frac >= -1 || $n!=$n || $n+1==$n
127             ### assert: 2*$frac < 1 || $n!=$n || $n+1==$n
128             }
129              
130 140 100       296 if ($n < 1) {
131 6         20 return;
132             }
133              
134             # d is the two-row number, ie. d=2*y, where n belongs
135             # start of the two-row group is nbase = 3 d^2 + d + 1
136             #
137 134         350 my $d = int ((_sqrtint(12*$n-11) - 1) / 6);
138 134         244 $n -= ((3*$d + 1)*$d + 1); # remainder within two-row
139             ### $d
140             ### remainder: $n
141 134 100       268 if ($n <= 3*$d) {
142             # 3d+1 many points in the Y=0,2,4,6 etc even row
143 83         126 $d *= 2; # y=2*d
144 83 100       281 return ($frac + $n + int(($n + ($self->{'mirror'} ? 2 : 0))/3) - $d,
145             $d);
146             } else {
147             # 3*d many points in the Y=1,3,5,7 etc odd row, using 3 in 4 cells
148 51         75 $n -= 3*$d+1; # remainder 0 upwards into odd row
149 51         77 $d = 2*$d+1; # y=2*d+1
150 51         133 return ($frac + $n + int($n/3) - $d,
151             $d);
152             }
153             }
154              
155             sub xy_to_n {
156 1628     1628 1 17143 my ($self, $x, $y) = @_;
157 1628         3164 $x = round_nearest ($x);
158 1628         3106 $y = round_nearest ($y);
159             ### CellularRule190 xy_to_n(): "$x,$y"
160              
161 1628 100 66     4952 if ($y < 0 || $x > $y) {
162 240         459 return undef;
163             }
164              
165 1388         1935 $x += $y; # move to have x=0 the start of the row
166 1388 100       2417 if ($x < 0) {
167 240         447 return undef;
168             }
169              
170             ### x centred: $x
171 1148 100       2023 if ($y % 2) {
172             ### odd row, 3s from the start ...
173 496 100       872 if (($x % 4) == 3) {
174 56         112 return undef;
175             }
176             # 3y^2+2y-1 = (3y-1)*(y+1)
177             return ($x
178             - int($x/4)
179             + ((3*$y+2)*$y-1)/4
180 440         1381 + $self->{'n_start'});
181             } else {
182             ### even row, 3s then 1 sep, or mirror 1 sep start ...
183 652         953 my $mirror = $self->{'mirror'};
184 652 100       1347 if (($x % 4) == ($mirror ? 1 : 3)) {
    100          
185 56         118 return undef;
186             }
187             return ($x
188             - int(($x+($mirror ? 2 : 1))/4)
189             + (3*$y+2)*$y/4
190 596 100       2093 + $self->{'n_start'});
191             }
192             }
193              
194             # exact
195             sub rect_to_n_range {
196 222     222 1 8547 my ($self, $x1,$y1, $x2,$y2) = @_;
197             ### CellularRule190 rect_to_n_range(): "$x1,$y1, $x2,$y2"
198              
199 222 50       582 ($x1,$y1, $x2,$y2) = _rect_for_V ($x1,$y1, $x2,$y2)
200             or return (1,0); # rect outside pyramid
201              
202             # # inherit bignum (before collapsing some y1 to x1 etc)
203             # my $zero = ($x1 * 0 * $y1 * $x2 * $y2);
204              
205 222         401 my $mirror = $self->{'mirror'};
206 222         345 my $unincremented_x1 = $x1;
207              
208             # \+------+
209             # | | /
210             # |\ | /
211             # | \ | /
212             # | \ | /
213             # y1 +------+ /
214             # x1 \ /
215             # \/
216             #
217 222 100       605 if ($x1 < (my $neg_y1 = -$y1)) {
    100          
    50          
218             ### bottom-left outside, move across to: "$neg_y1,$y1"
219 11         15 $x1 = $neg_y1;
220              
221             # For the following blank checks a blank doesn't occur at the ends of a
222             # row, so when on a blank it's always possible to increment or decrement
223             # X to go to a non-blank -- as long as that adjacent space is within the
224             # rectangle.
225             #
226             } elsif ((($mirror ? $y1-$x1 : $x1+$y1) % 4) == 3) {
227             ### x1,y1 bottom left is on a blank: "x1+y1=".($x1+$y1)
228 0 0       0 if ($x1 < $x2) {
229             ### rect width >= 2, so increase x1 ...
230 0         0 $x1 += 1;
231             } else {
232             ### rect is a single column width==1, increase y1 ...
233 0 0       0 if (($y1 += 1) > $y2) {
234             ### rect was a single blank square, contains no N ...
235 0         0 return (1,0);
236             }
237             }
238             }
239              
240 222 100       588 if ((($mirror ? $y2-$x2 : $x2+$y2) % 4) == 3) {
    100          
241             ### x2,y2 top right is on a blank, decrement ...
242 12 50       26 if ($x2 > $unincremented_x1) {
243             ### rect width >= 2, so decrease x2 ...
244 12         19 $x2 -= 1;
245             } else {
246             ### rect is a single column width==1, decrease y2 ...
247 0         0 $y2 -= 1;
248              
249             # Can decrement without checking whether the rect is a single square.
250             # If the rect was a single blank square then the x1+y1 bottom left
251             # above detects and returns. And the bottom left blank check
252             # incremented y1 to leave a single square then that's a non-blank
253             # because there's no vertical blank pairs (they go on the diagonal).
254             ### assert $y2 >= $y1
255             }
256             }
257              
258             # At this point $x1,$y1 is a non-blank bottom left corner, and $x2,$y2
259             # is a non-blank top right corner, being the N lo to hi range.
260              
261             ### range: "bottom-right $x1,$y1 top-left $x2,$y2"
262 222         459 return ($self->xy_to_n ($x1,$y1),
263             $self->xy_to_n ($x2,$y2));
264             }
265              
266              
267             # old rect_to_n_range() of row endpoints
268             #
269             # # even right y = [ 0, 2, 4, 6 ]
270             # # N = [ 1,8,21,40 ]
271             # # Nright = (3/4 y^2 + 2 y + 1)
272             # # = (3 y^2 + 8 y + 4) / 4
273             # # = ((3y + 8)y + 4) / 4
274             # #
275             # # odd right y = [ 1, 3, 5, 7 ]
276             # # N = [ 4,14,30, 52 ]
277             # # Nright = (3/4 y^2 + 2 y + 5/4)
278             # # = (3 y^2 + 8 y + 5) / 4
279             # # = ((3y + 8)y + 5) / 4
280             # #
281             # # Nleft y even ((3y+2)*y + 4)/4
282             # # Nleft y odd ((3y+2)*y + 3)/4
283             # # Nright even ((3(y+1)+2)*(y+1) + 3)/4 - 1
284             # # = ((3y+3+2)*(y+1) + 3 - 4)/4
285             # # = ((3y+5)*(y+1) - 1)/4
286             # # = ((3y^2 + 8y + 5 - 1)/4
287             # # = ((3y^2 + 8y + 4)/4
288             # # = ((3y+8)y + 4)/4
289             # # = ((3y+2)(y+2)/4
290             # #
291             # ### $y1
292             # ### $y2
293             # $y2 += $zero;
294             # $y1 += $zero;
295             # return (((3*$y1 + 2)*$y1 + 4 - ($y1%2)) / 4, # even/odd Nleft
296             # ((3*$y2 + 8)*$y2 + 4 + ($y2%2)) / 4); # even/odd Nright
297              
298             1;
299             __END__