File Coverage

blib/lib/Math/PlanePath/KochSquareflakes.pm
Criterion Covered Total %
statement 212 245 86.5
branch 68 90 75.5
condition 13 14 92.8
subroutine 23 25 92.0
pod 5 5 100.0
total 321 379 84.7


line stmt bran cond sub pod time code
1             # Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020 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 it
6             # under the terms of the GNU General Public License as published by the Free
7             # 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::KochSquareflakes;
20 1     1   7642 use 5.004;
  1         8  
21 1     1   4 use strict;
  1         2  
  1         18  
22              
23 1     1   4 use vars '$VERSION', '@ISA';
  1         2  
  1         53  
24             $VERSION = 128;
25 1     1   551 use Math::PlanePath;
  1         2  
  1         41  
26             @ISA = ('Math::PlanePath');
27             *_divrem = \&Math::PlanePath::_divrem;
28              
29             use Math::PlanePath::Base::Generic
30 1         37 'is_infinite',
31 1     1   5 'round_nearest';
  1         1  
32             use Math::PlanePath::Base::Digits
33 1         69 'round_down_pow',
34 1     1   368 'digit_split_lowtohigh';
  1         1  
35              
36             # uncomment this to run the ### lines
37             #use Devel::Comments;
38              
39              
40 1     1   6 use constant n_frac_discontinuity => 0;
  1         1  
  1         62  
41              
42 1         39 use constant parameter_info_array =>
43             [ { name => 'inward',
44             display => 'Inward',
45             type => 'boolean',
46             default => 0,
47             description => 'Whether to direct the sides of the square inward, rather than outward.',
48 1     1   6 } ];
  1         1  
49              
50 1     1   5 use constant x_negative_at_n => 1;
  1         1  
  1         40  
51 1     1   6 use constant y_negative_at_n => 1;
  1         1  
  1         34  
52 1     1   5 use constant sumabsxy_minimum => 1;
  1         1  
  1         39  
53 1     1   6 use constant rsquared_minimum => 0.5; # minimum X=0.5,Y=0.5
  1         1  
  1         52  
54              
55             # jump across rings is South-West, so
56 1     1   5 use constant dx_maximum => 1;
  1         2  
  1         49  
57 1     1   5 use constant dy_maximum => 1;
  1         2  
  1         33  
58 1     1   5 use constant dsumxy_maximum => 2; # diagonal NE
  1         1  
  1         42  
59 1     1   5 use constant ddiffxy_maximum => 2;
  1         1  
  1         42  
60 1     1   5 use constant ddiffxy_minimum => -2;
  1         2  
  1         36  
61 1     1   4 use constant dir_maximum_dxdy => (1,-1); # South-East
  1         2  
  1         46  
62              
63             # N=1,2,3,4 gcd(1/2,1/2) = 1/2
64 1     1   6 use constant gcdxy_minimum => 1/2;
  1         1  
  1         42  
65              
66 1     1   6 use constant turn_any_straight => 0; # never straight
  1         1  
  1         1155  
67              
68              
69             #------------------------------------------------------------------------------
70              
71             # level 0 inner square
72             # sidelen = 4^level
73             # ring points 4*4^level
74             # Nend = 4 * [ 1 + ... + 4^level ]
75             # = 4 * (4^(level+1) - 1) / 3
76             # = (4^(level+2) - 4) / 3
77             # Nstart = Nend(level-1) + 1
78             # = (4^(level+1) - 4) / 3 + 1
79             # = (4^(level+1) - 4 + 3) / 3
80             # = (4^(level+1) - 1) / 3
81             #
82             # level Nstart Nend
83             # 0 (4-1)/3=1 (16-4)/3=12/3=4
84             # 1 (16-1)/3=15/3=5 (64-4)/3=60/3=20
85             # 2 (64-1)/3=63/3=21 (256-4)/3=252/3=84
86             # 3 (256-1)/3=255/3=85
87             #
88              
89             sub n_to_xy {
90 190     190 1 5300 my ($self, $n) = @_;
91             ### KochSquareflakes n_to_xy(): $n
92 190 50       331 if ($n < 1) { return; }
  0         0  
93              
94 190         203 my $frac;
95             {
96 190         208 my $int = int($n);
  190         237  
97 190         249 $frac = $n - $int;
98 190         215 $n = $int; # BigFloat int() gives BigInt, use that
99             }
100              
101             # (4^(level+1) - 1) / 3 = N
102             # 4^(level+1) - 1 = 3*N
103             # 4^(level+1) = 3*N+1
104             #
105 190         428 my ($pow,$level) = round_down_pow (3*$n + 1, 4);
106             ### $level
107             ### $pow
108 190 50       366 if (is_infinite($level)) { return ($level,$level); }
  0         0  
109              
110             # Nstart = (4^(level+1)-1)/3 with $power=4^(level+1) here
111             #
112 190         353 $n -= ($pow-1)/3;
113              
114             ### base: ($pow-1)/3
115             ### next base would be: (4*$pow-1)/3
116             ### n remainder from base: $n
117              
118 190         242 my $sidelen = $pow/4;
119 190         351 (my $rot, $n) = _divrem ($n, $sidelen); # high part is rot
120             ### $sidelen
121             ### n remainder: $n
122             ### $rot
123              
124             ### assert: $n>=0
125             ### assert: $n < 4 ** $level
126              
127 190         306 my @horiz = (1);
128 190         228 my @diag = (1);
129 190         221 my $i = 0;
130 190         341 while (--$level > 0) {
131 386         593 $horiz[$i+1] = 2*$horiz[$i] + 2*$diag[$i];
132 386         554 $diag[$i+1] = $horiz[$i] + 2*$diag[$i];
133             ### horiz: $horiz[$i+1]
134             ### diag: $diag[$i+1]
135 386         551 $i++;
136             }
137              
138             ### horiz: join(', ',@horiz)
139             ### $i
140 190         346 my $x
141             = my $y
142             = ($n * 0) + $horiz[$i]/-2; # inherit bignum
143 190 100       311 if ($rot & 1) {
144 87         155 ($x,$y) = (-$y,$x);
145             }
146 190 100       282 if ($rot & 2) {
147 87         114 $x = -$x;
148 87         109 $y = -$y;
149             }
150 190         205 $rot *= 2;
151              
152 190         239 my $inward = $self->{'inward'};
153 190         396 my @digits = digit_split_lowtohigh($n,4);
154              
155 190         303 while ($i > 0) {
156 386         486 $i--;
157 386   100     693 my $digit = $digits[$i] || 0;
158              
159 386         436 my ($dx, $dy, $drot);
160 386 100       681 if ($digit == 0) {
    100          
    100          
    50          
161 149         155 $dx = 0;
162 149         150 $dy = 0;
163 149         155 $drot = 0;
164             } elsif ($digit == 1) {
165 72 100       101 if ($rot & 1) {
166 16         19 $dx = $diag[$i];
167 16         21 $dy = $diag[$i];
168             } else {
169 56         73 $dx = $horiz[$i];
170 56         62 $dy = 0;
171             }
172 72 100       96 $drot = ($inward ? 1 : -1);
173             } elsif ($digit == 2) {
174 72 100       100 if ($rot & 1) {
175 16 100       25 if ($inward) {
176 8         11 $dx = $diag[$i];
177 8         11 $dy = $diag[$i] + $horiz[$i];
178             } else {
179 8         9 $dx = $diag[$i] + $horiz[$i];
180 8         11 $dy = $diag[$i];
181             }
182             } else {
183 56         73 $dx = $horiz[$i] + $diag[$i];
184 56         65 $dy = $diag[$i];
185 56 100       86 unless ($inward) { $dy = -$dy; }
  28         33  
186             }
187 72 100       90 $drot = ($inward ? -1 : 1);
188             } elsif ($digit == 3) {
189 93 100       132 if ($rot & 1) {
190 16         24 $dx = $dy = $diag[$i] + $horiz[$i];
191             } else {
192 77         93 $dx = $horiz[$i] + 2*$diag[$i];
193 77         96 $dy = 0;
194             }
195 93         104 $drot = 0;
196             }
197             ### delta: "$dx,$dy rot=$rot drot=$drot"
198              
199 386 100       559 if ($rot & 2) {
200 165         240 ($dx,$dy) = (-$dy,$dx);
201             }
202 386 100       563 if ($rot & 4) {
203 165         188 $dx = -$dx;
204 165         176 $dy = -$dy;
205             }
206             ### delta with rot: "$dx,$dy"
207              
208 386         410 $x += $dx;
209 386         423 $y += $dy;
210 386         822 $rot += $drot;
211             }
212              
213             {
214 190         212 my $dx = $frac;
  190         215  
215 190 100       268 my $dy = ($rot & 1 ? $frac : 0);
216 190 100       259 if ($rot & 2) {
217 87         136 ($dx,$dy) = (-$dy,$dx);
218             }
219 190 100       279 if ($rot & 4) {
220 87         97 $dx = -$dx;
221 87         98 $dy = -$dy;
222             }
223 190         206 $x = $dx + $x;
224 190         222 $y = $dy + $y;
225             }
226              
227 190         453 return ($x,$y);
228             }
229              
230             my @inner_to_n = (1,2,4,3);
231              
232             sub xy_to_n {
233 880     880 1 45318 my ($self, $x, $y) = @_;
234             ### KochSquareflakes xy_to_n(): "$x, $y"
235              
236             # +/- 0.75
237 880 50 100     2245 if (4*$x < 3 && 4*$y < 3 && 4*$x >= -3 && 4*$y >= -3) {
      100        
      66        
238 0         0 return $inner_to_n[($x >= 0) + 2*($y >= 0)];
239             }
240              
241 880         1468 $x = round_nearest($x);
242 880         1289 $y = round_nearest($y);
243              
244             # quarter curve segment and high digit
245 880         1059 my $n;
246             {
247 880         949 my $negx = -$x;
  880         974  
248 880 100       1409 if (($y > 0 ? $x > $y : $x >= $y)) {
    100          
249             ### below leading diagonal ...
250 440 100       587 if ($negx > $y) {
251             ### bottom quarter ...
252 220         243 $n = 1;
253             } else {
254             ### right quarter ...
255 220         237 $n = 2;
256 220         320 ($x,$y) = ($y, $negx); # rotate -90
257             }
258             } else {
259             ### above leading diagonal
260 440 100       562 if ($y > $negx) {
261             ### top quarter ...
262 220         233 $n = 3;
263 220         250 $x = $negx; # rotate 180
264 220         282 $y = -$y;
265             } else {
266             ### right quarter ...
267 220         236 $n = 4;
268 220         384 ($x,$y) = (-$y, $x); # rotate +90
269             }
270             }
271             }
272 880         995 $y = -$y;
273             ### rotate to: "$x,$y n=$n"
274              
275 880 50       1320 if (is_infinite($x)) {
276 0         0 return $x;
277             }
278 880 50       1485 if (is_infinite($y)) {
279 0         0 return $y;
280             }
281              
282 880         1212 my @horiz;
283             my @diag;
284 880         967 my $horiz = 1;
285 880         977 my $diag = 1;
286 880         921 for (;;) {
287 1664         1986 push @horiz, $horiz;
288 1664         1770 push @diag, $diag;
289 1664         1823 my $offset = $horiz+$diag;
290 1664         1827 my $nextdiag = $offset + $diag; # horiz + 2*diag
291             ### $horiz
292             ### $diag
293             ### $offset
294             ### $nextdiag
295              
296 1664 100       2356 if ($y <= $nextdiag) {
297             ### found level at: "top=$nextdiag vs y=$y"
298 880         1003 $y -= $offset;
299 880         1008 $x += $offset;
300 880         1052 last;
301             }
302 784         837 $horiz = 2*$offset; # 2*horiz+2*diag
303 784         938 $diag = $nextdiag;
304             }
305             ### base subtract to: "$x,$y"
306              
307              
308 880 100       1339 if ($self->{'inward'}) {
309 440         502 $y = -$y;
310             ### inward invert to: "$x,$y"
311             }
312              
313             ### origin based side: "$x,$y horiz=$horiz diag=$diag with levels ".scalar(@horiz)
314              
315             # loop 4*1, 4*4, 4*4^2 etc, extra +1 on the digits to include that in the sum
316             #
317 880         1017 my $slope;
318 880         1359 while (@horiz) {
319             ### at: "$x,$y slope=".($slope||0)." n=$n"
320 1664         2286 $horiz = pop @horiz;
321 1664         1809 $diag = pop @diag;
322 1664         1783 $n *= 4;
323              
324 1664 100       2123 if ($slope) {
325 336 100       440 if ($y < $diag) {
326             ### slope digit 0 ...
327 152         238 $n += 1;
328             } else {
329 184         200 $x -= $diag;
330 184         195 $y -= $diag;
331             ### slope not digit 0, move to: "$x,$y"
332              
333 184 100       235 if ($y < $horiz) {
334             ### digit 1 ...
335 80         88 $n += 2;
336 80         115 ($x,$y) = ($y, -$x); # rotate -90
337 80         185 $slope = 0;
338             } else {
339 104         112 $y -= $horiz;
340             ### slope not digit 1, move to: "$x,$y"
341              
342 104 100       143 if ($x < $horiz) {
343             ### digit 2 ...
344 48         51 $n += 3;
345 48         85 $slope = 0;
346              
347             } else {
348             ### digit 3 ...
349 56         63 $n += 4;
350 56         90 $x -= $horiz;
351             }
352             }
353             }
354              
355             } else {
356 1328 100       1658 if ($x < $horiz) {
357             ### digit 0 ...
358 384         569 $n += 1;
359             } else {
360 944         1046 $x -= $horiz;
361             ### not digit 0, move to: "$x,$y"
362              
363 944 100       1177 if ($x < $diag) {
364             ### digit 1 ...
365 280         315 $n += 2;
366 280         443 $slope = 1;
367             } else {
368 664         725 $x -= $diag;
369             ### not digit 1, move to: "$x,$y"
370              
371 664 100       824 if ($x < $diag) {
372             ### digit 2 ...
373 280         293 $n += 3;
374 280         297 $slope = 1;
375 280         589 ($x,$y) = ($diag-$y, $x); # offset and rotate +90
376              
377             } else {
378             ### digit 3 ...
379 384         401 $n += 4;
380 384         588 $x -= $diag;
381             }
382             }
383             }
384             }
385             }
386             ### final: "$x,$y n=$n"
387              
388 880 100 100     1827 if ($x == 0 && $y == 0) {
389 160         293 return $n;
390             } else {
391 720         1222 return undef;
392             }
393             }
394              
395             # not exact
396             sub rect_to_n_range {
397 0     0 1 0 my ($self, $x1,$y1, $x2,$y2) = @_;
398             ### KochSquareflakes rect_to_n_range(): "$x1,$y1 $x2,$y2"
399              
400 0         0 foreach ($x1,$y1, $x2,$y2) {
401 0 0       0 if (is_infinite($_)) {
402 0         0 return (0, $_);
403             }
404 0         0 $_ = abs(round_nearest($_));
405             }
406 0 0       0 if ($x1 > $x2) { ($x1,$x2) = ($x2,$x1); }
  0         0  
407 0 0       0 if ($y1 > $y2) { ($y1,$y2) = ($y2,$y1); }
  0         0  
408 0 0       0 my $max = ($x2 > $y2 ? $x2 : $y2);
409              
410             # Nend = 4 * [ 1 + ... + 4^level ]
411             # = 4 + 16 + ... + 4^(level+1)
412             #
413 0         0 my $horiz = 4;
414 0         0 my $diag = 3;
415 0         0 my $nhi = 4;
416 0         0 for (;;) {
417 0         0 $nhi += 1;
418 0         0 $nhi *= 4;
419 0         0 my $nextdiag = $horiz + 2*$diag;
420 0 0       0 if (($self->{'inward'} ? $horiz : $nextdiag) >= 2*$max) {
    0          
421 0         0 return (1, $nhi);
422             }
423 0         0 $horiz = $nextdiag + $horiz; # 2*$horiz + 2*$diag;
424 0         0 $diag = $nextdiag;
425             }
426             }
427              
428              
429             #------------------------------------------------------------------------------
430             # Nstart = (4^(k+1) - 1)/3
431             # Nend = Nstart(k+1) - 1
432             # = (4*4^(k+1) - 1)/3 - 1
433             # = (4*4^(k+1) - 1 - 3)/3
434             # = (4*4^(k+1) - 4)/3
435             # = 4*(4^(k+1) - 1)/3
436             # = 4*Nstart(k)
437              
438             sub level_to_n_range {
439 10     10 1 784 my ($self, $level) = @_;
440 10         25 my $n_lo = (4**($level+1) - 1)/3;
441 10         25 return ($n_lo, 4*$n_lo);
442             }
443             sub n_to_level {
444 0     0 1   my ($self, $n) = @_;
445 0 0         if ($n < 1) { return undef; }
  0            
446 0 0         if (is_infinite($n)) { return $n; }
  0            
447 0           my ($pow,$exp) = round_down_pow (3*$n + 1, 4);
448 0           return $exp-1;
449             }
450              
451             #------------------------------------------------------------------------------
452             1;
453             __END__