File Coverage

blib/lib/Math/PlanePath/KochSnowflakes.pm
Criterion Covered Total %
statement 125 150 83.3
branch 35 56 62.5
condition 11 14 78.5
subroutine 26 27 96.3
pod 7 7 100.0
total 204 254 80.3


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             # math-image --path=KochSnowflakes --lines --scale=10
20             #
21             # area approaches sqrt(48)/10
22             # * height=sqrt(1-1/4)=sqrt(3)/2
23             # /|\ halfbase=1/2
24             # / | \ trianglearea = sqrt(3)/4
25             # *-----*
26             # segments = 3*4^level = 3,12,48,192,...
27             #
28             # with initial triangle area=1
29             # add a new triangle onto each side
30             # x,y scale by 3* so 9*area
31             #
32             # area[level+1] = 9*area[level] + segments
33             # = 9*area[level] + 3*4^level
34              
35             # area[0] = 1
36             # area[1] = 9*area[0] + 3 = 9 + 3 = 12
37             # area[2] = 9*area[1] + 3*4
38             # = 9*(9*1 + 3) + 3*4
39             # = 9*9 + 3*9 + 3*4 = 120
40             # area[3] = 9*area[2] + 3*4
41             # = 9*(9*9 + 3*9 + 3*4) + 3*4^2
42             # = 9^3 + 3*9^2 + 3*0*4 + 3*4^2
43              
44             # area[level+1]
45             # = 9^(level+1) + (9^(level+1) - 4^(level+1)) * 3/5
46             # = (5*9^(level+1) + 3*9^(level+1) - 3*4^(level+1)) / 5
47             # = (8*9^(level+1) - 3*4^(level+1)) / 5
48             #
49             # area[level] = (8*9^level - 3*4^level) / 5
50             # = 1,12,120,1128,10344,93864,847848
51             #
52             # .
53             # / \ area[0] = 1
54             # .---.
55             #
56             # .
57             # / \ area=[1] = 12 = 9*area[0] + 3*4^0
58             # .---.---.---.
59             # \ / \ / \ /
60             # .---.---.
61             # / \ / \ / \
62             # .---.---.---.
63             # \ /
64             # .
65             #
66             # area[level] / 9^level
67             # = (8*9^level / 9^level - 3*4^level / 9^level) / 5
68             # = (8 - 3*(4/0)^level)/5
69             # -> 8/5 as level->infinity
70              
71             # in integer coords
72             # initial triangle area
73             # * 2/3 1*2 / 2 = 1 unit
74             # / \
75             # *---* -1/3
76             # -1 +1
77             #
78             # so area[level] / (sqrt(3)/2)
79              
80              
81              
82             package Math::PlanePath::KochSnowflakes;
83 1     1   1507 use 5.004;
  1         4  
84 1     1   5 use strict;
  1         2  
  1         56  
85             #use List::Util 'max';
86             *max = \&Math::PlanePath::_max;
87              
88 1     1   7 use vars '$VERSION', '@ISA';
  1         1  
  1         83  
89             $VERSION = 127;
90 1     1   790 use Math::PlanePath;
  1         3  
  1         41  
91             @ISA = ('Math::PlanePath');
92              
93             use Math::PlanePath::Base::Generic
94 1         43 'is_infinite',
95 1     1   7 'round_nearest';
  1         2  
96             use Math::PlanePath::Base::Digits
97 1     1   580 'round_down_pow';
  1         2  
  1         54  
98 1     1   671 use Math::PlanePath::KochCurve;
  1         2  
  1         33  
99              
100             # uncomment this to run the ### lines
101             # use Smart::Comments;
102              
103              
104 1     1   7 use constant n_frac_discontinuity => 0;
  1         2  
  1         54  
105 1     1   6 use constant x_negative_at_n => 1;
  1         2  
  1         43  
106 1     1   6 use constant y_negative_at_n => 1;
  1         2  
  1         44  
107 1     1   6 use constant sumabsxy_minimum => 2/3; # minimum X=0,Y=2/3
  1         2  
  1         42  
108 1     1   5 use constant rsquared_minimum => 4/9; # minimum X=0,Y=2/3
  1         2  
  1         108  
109             # maybe: use constant radius_minimum => 2/3; # minimum X=0,Y=2/3
110              
111             # jump across rings is WSW slope 2, so following maximums
112 1     1   8 use constant dx_maximum => 2;
  1         1  
  1         44  
113 1     1   6 use constant dy_maximum => 1;
  1         2  
  1         37  
114 1     1   5 use constant dsumxy_maximum => 2;
  1         2  
  1         52  
115 1     1   6 use constant ddiffxy_maximum => 2;
  1         2  
  1         42  
116              
117 1     1   6 use constant absdx_minimum => 1; # never vertical
  1         1  
  1         56  
118 1     1   6 use constant dir_maximum_dxdy => (1,-1); # South-East
  1         2  
  1         57  
119              
120             # N=1 gcd(-1, -1/3) = 1/3
121             # N=2 gcd( 1, -1/3) = 1/3
122             # N=3 gcd( 0, 2/3) = 2/3
123 1     1   6 use constant gcdxy_minimum => 1/3;
  1         2  
  1         41  
124              
125 1     1   6 use constant turn_any_straight => 0; # never straight
  1         2  
  1         1009  
126              
127              
128             #------------------------------------------------------------------------------
129             sub new {
130 4     4 1 312 my $self = shift->SUPER::new (@_);
131 4   50     29 $self->{'sides'} ||= 3; # default
132 4         10 return $self;
133             }
134              
135             # N=1 to 3 3 of, level=1
136             # N=4 to 15 12 of, level=2
137             # N=16 to .. 48 of, level=3
138             #
139             # each loop = 3*4^level
140             #
141             # n_base = 1 + 3*4^0 + 3*4^1 + ... + 3*4^(level-1)
142             # = 1 + 3*[ 4^0 + 4^1 + ... + 4^(level-1) ]
143             # = 1 + 3*[ (4^level - 1)/3 ]
144             # = 1 + (4^level - 1)
145             # = 4^level
146             #
147             # each side = loop/3
148             # = 3*4^level / 3
149             # = 4^level
150             #
151             # 6 sides
152             # n_base = 1 + 2*3*4^0 + ...
153             # = 2*4^level - 1
154             # level = log4 (n+1)/2
155              
156             ### loop 1: 3* 4**1
157             ### loop 2: 3* 4**2
158             ### loop 3: 3* 4**3
159              
160             # sub _level_to_base {
161             # my ($level) = @_;
162             # return -3*$level + 4**($level+1) - 2;
163             # }
164             # ### level_to_base(1): _level_to_base(1)
165             # ### level_to_base(2): _level_to_base(2)
166             # ### level_to_base(3): _level_to_base(3)
167              
168             sub n_to_xy {
169 18     18 1 697 my ($self, $n) = @_;
170             ### KochSnowflakes n_to_xy(): $n
171 18 50       47 if ($n < 1) { return; }
  0         0  
172 18 50       39 if (is_infinite($n)) { return ($n,$n); }
  0         0  
173              
174 18         50 my $sides = $self->{'sides'};
175 18 50       63 my ($sidelen, $level) = round_down_pow (($sides == 6 ? ($n+1)/2 : $n),
176             4);
177 18 50       34 my $base = ($sides == 6 ? 2*$sidelen - 1 : $sidelen);
178 18         30 my $rem = $n - $base;
179              
180             ### $level
181             ### $base
182             ### $sidelen
183             ### $rem
184             ### assert: $n >= $base
185             ### assert: $rem >= 0
186             ### assert: $rem < $sidelen * $sides
187              
188 18         30 my $side = int($rem / $sidelen);
189             ### $side
190             ### $rem
191             ### subtract: $side*$sidelen
192 18         26 $rem -= $side*$sidelen;
193              
194             ### assert: $side >= 0 && $side < $sides
195              
196 18         47 my ($x, $y) = Math::PlanePath::KochCurve->n_to_xy ($rem);
197             ### $x
198             ### $y
199              
200 18 50       36 if ($sides == 3) {
201 18         28 my $len = 3**($level-1);
202 18 100       50 if ($side < 1) {
    50          
203             ### horizontal rightwards
204 12         40 return ($x - 3*$len,
205             -$y - $len);
206             } elsif ($side < 2) {
207             ### right slope upwards
208 0         0 return (($x-3*$y)/-2 + 3*$len, # flip vert and rotate +120
209             ($x+$y)/2 - $len);
210             } else {
211             ### left slope downwards
212 6         23 return ((-3*$y-$x)/2, # flip vert and rotate -120
213             ($y-$x)/2 + 2*$len);
214             }
215             } else {
216              
217             # 3
218             # 5-----4
219             # 4 / \
220             # / \ 2
221             # 6 o 3
222             # 5 \ . . /
223             # \ . . / 1
224             # 1-----2
225             # 0
226             # 7
227             #
228 0         0 my $len = 3**$level;
229 0         0 $x -= $len; # -y flip vert and offset
230 0         0 $y = -$y - $len;
231 0 0       0 if ($side >= 3) {
232             ### rotate 180 ...
233 0         0 $x = -$x; # rotate 180
234 0         0 $y = -$y;
235 0         0 $side -= 3;
236             }
237 0 0       0 if ($side >= 2) {
238             ### upper right slope upwards ...
239 0         0 return (($x+3*$y)/-2, # rotate +120
240             ($x-$y)/2);
241             }
242 0 0       0 if ($side >= 1) {
243             ### lower right slope upwards ...
244 0         0 return (($x-3*$y)/2, # rotate +60
245             ($x+$y)/2);
246             }
247             ### horizontal ...
248 0         0 return ($x,$y);
249             }
250             }
251              
252              
253             # N=1 overlaps N=5
254             # N=2 overlaps N=7
255             # +---------+ +---------+ Y=1.5
256             # | | | |
257             # | +---------+ | Y=7/6 = 1.166
258             # | | | |
259             # | * 13 | | * 11 | Y=1
260             # | | | |
261             # | | * 3 | | Y=2/3 = 0.666
262             # | | | |
263             # +---------+ +---------+ Y=0.5
264             # | |
265             # +---------+---------+---------+ Y=1/6 = 0.166
266             # | | O | | --Y=0
267             # | | | |
268             # | | | |
269             # | * 1 | | * 2 | Y=-1/3 = -0.333
270             # | | | |
271             # +---------+ +---------+ Y=-3/6 = -0.5
272             # | | | |
273             # +---------+ +---------+ Y=-5/6 = -0.833
274             # | | | |
275             # | * 5 | | * 7 | Y=-1
276             # | | | |
277             # | | | |
278             # +---------+ +---------+ Y=-1.5
279             #
280             sub xy_to_n {
281 40     40 1 3443 return scalar((shift->xy_to_n_list(@_))[0]);
282             }
283             sub xy_to_n_list {
284 58     58 1 1993 my ($self, $x, $y) = @_;
285             ### KochSnowflakes xy_to_n(): "$x, $y"
286              
287 58         134 $x = round_nearest ($x);
288 58 100       126 if (abs($x) <= 1) {
289 44 100       78 if ($x == 0) {
290 13         22 my $y6 = 6*$y;
291 13 100 66     44 if ($y6 >= 1 && $y6 < 7) {
292             # Y = 2/3-1/2=1/6 to 2/3+1/2=7/6
293 9         20 return 3;
294             }
295             } else {
296 31         47 my $y6 = 6*$y;
297 31 100 66     105 if ($y6 >= -5 && $y6 < 1) {
298             # Y = -1/3-1/2=-5/6 to -1/3+1/2=+1/6
299 24 100       76 return (1 + ($x > 0),
300             ($y6 < -3 ? (5+2*($x>0)) : ())); # 5 or 7 up to Y<-1/2
301             }
302             }
303             }
304              
305 25         48 $y = round_nearest ($y);
306 25 100       63 if (($x % 2) != ($y % 2)) {
307             ### diff parity...
308 9         18 return;
309             }
310              
311 16         20 my $high;
312 16 100 100     77 if ($x > 0 && $x >= -3*$y) {
    100 100        
313             ### right upper third n=2 ...
314 3         7 ($x,$y) = ((3*$y-$x)/2, # rotate -120 and flip vert
315             ($x+$y)/2);
316 3         7 $high = 2;
317             } elsif ($x <= 0 && 3*$y > $x) {
318             ### left upper third n=3 ...
319 3         8 ($x,$y) = (($x+3*$y)/-2, # rotate +120 and flip vert
320             ($y-$x)/2);
321 3         5 $high = 3;
322             } else {
323             ### lower third n=1 ...
324 10         16 $y = -$y; # flip vert
325 10         15 $high = 1;
326             }
327             ### rotate/flip is: "$x,$y"
328              
329 16 50       34 if ($y <= 0) {
330 0         0 return;
331             }
332              
333 16         40 my ($len,$level) = round_down_pow($y, 3);
334 16         31 $level += 1;
335             ### $level
336             ### $len
337 16 50       35 if (is_infinite($level)) {
338 0         0 return $level;
339             }
340              
341              
342 16         28 $y -= $len; # shift to Y=0 basis
343 16         22 $len *= 3;
344              
345             ### compare for end: ($x+$y)." >= 3*len=".$len
346 16 100       34 if ($x + $y >= $len) {
347             ### past end of this level, no points ...
348 3         9 return;
349             }
350 13         19 $x += $len; # shift to X=0 basis
351              
352 13         39 my $n = Math::PlanePath::KochCurve->xy_to_n($x, $y);
353              
354             ### plain curve on: ($x+3*$len).",".($y-$len)." n=".(defined $n && $n)
355             ### $high
356             ### high: (4**$level)*$high
357              
358 13 100       28 if (defined $n) {
359 10         34 return (4**$level)*$high + $n;
360             } else {
361 3         7 return;
362             }
363             }
364              
365             # level extends to x= +/- 3^level
366             # y= +/- 2*3^(level-1)
367             # = 2/3 * 3^level
368             # 1.5*y = 3^level
369             #
370             # ENHANCE-ME: share KochCurve segment checker to find actual min/max
371             #
372             # not exact
373             sub rect_to_n_range {
374 11     11 1 1266 my ($self, $x1,$y1, $x2,$y2) = @_;
375             ### KochSnowflakes rect_to_n_range(): "$x1,$y1 $x2,$y2"
376              
377 11         30 $x1 = round_nearest ($x1);
378 11         22 $y1 = round_nearest ($y1);
379 11         23 $x2 = round_nearest ($x2);
380 11         22 $y2 = round_nearest ($y2);
381              
382 11 100       26 ($x1,$x2) = ($x2,$x1) if $x1 > $x2;
383 11 100       34 ($y1,$y2) = ($y2,$y1) if $y1 > $y2;
384              
385             #
386             # |
387             # +------ . -----+
388             # |x1,y2 /|\ x2,y2|
389             # / | \
390             # / | \
391             # -----/---m---\-----
392             # / | \
393             # .-----------.
394             # |
395             # y1
396             # -------
397             #
398             # -y1 bottom horizontal
399             # (x2+y2)/2 right side
400             # (-x1+y2)/2 left side
401             # each giving a power of 3 of the level
402             #
403             ### right: ($x2+$y2)/2
404             ### left: (-$x1+$y2)/2
405             ### bottom: -$y1
406              
407 11         21 my $sides = $self->{'sides'};
408 11 50       49 my ($pow, $level) = round_down_pow (max ($sides == 6
409             ? ($x1/-2,
410             $x2/2,
411             -$y1,
412             $y2)
413             : (int(($x2+$y2)/2),
414             int((-$x1+$y2)/2),
415             -$y1)),
416             3);
417             ### $level
418             # end of $level is 1 before base of $level+1
419 11         36 return (1, 4**($level+2) - 1);
420             }
421              
422             #------------------------------------------------------------------------------
423             # Nstart = 4^k
424             # length = 3*4^k many points
425             # Nend = Nstart + length-1
426             # = 4^k + 3*4^k - 1
427             # = 4*4^k - 1
428             # = Nstart(k+1) - 1
429             sub level_to_n_range {
430 9     9 1 953 my ($self, $level) = @_;
431 9         17 my $pow = 4**$level;
432 9         30 return ($pow, 4*$pow-1);
433             }
434             sub n_to_level {
435 0     0 1   my ($self, $n) = @_;
436 0 0         if ($n < 1) { return undef; }
  0            
437 0 0         if (is_infinite($n)) { return $n; }
  0            
438 0           $n = round_nearest($n);
439 0 0         my ($sidelen, $level) = round_down_pow (($self->{'sides'} == 6 ? ($n+1)/2 : $n),
440             4);
441 0           return $level;
442             }
443              
444             #------------------------------------------------------------------------------
445             1;
446             __END__