File Coverage

blib/lib/Math/PlanePath/CCurve.pm
Criterion Covered Total %
statement 138 192 71.8
branch 33 70 47.1
condition 20 22 90.9
subroutine 20 25 80.0
pod 8 8 100.0
total 219 317 69.0


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
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             package Math::PlanePath::CCurve;
20 1     1   2022 use 5.004;
  1         4  
21 1     1   6 use strict;
  1         1  
  1         36  
22 1     1   7 use List::Util 'min','max','sum';
  1         2  
  1         84  
23              
24 1     1   6 use vars '$VERSION', '@ISA';
  1         2  
  1         71  
25             $VERSION = 128;
26 1     1   847 use Math::PlanePath;
  1         2  
  1         30  
27 1     1   604 use Math::PlanePath::Base::NSEW;
  1         2  
  1         41  
28             @ISA = ('Math::PlanePath::Base::NSEW',
29             'Math::PlanePath');
30              
31             use Math::PlanePath::Base::Generic
32 1         49 'is_infinite',
33 1     1   6 'round_nearest';
  1         2  
34             use Math::PlanePath::Base::Digits
35 1         85 'round_up_pow',
36             'round_down_pow',
37             'bit_split_lowtohigh',
38             'digit_split_lowtohigh',
39 1     1   641 'digit_join_lowtohigh';
  1         3  
40             *_divrem = \&Math::PlanePath::_divrem;
41             *_divrem_mutate = \&Math::PlanePath::_divrem_mutate;
42              
43 1     1   714 use Math::PlanePath::KochCurve;
  1         3  
  1         50  
44             *_digit_join_hightolow = \&Math::PlanePath::KochCurve::_digit_join_hightolow;
45              
46             # uncomment this to run the ### lines
47             # use Smart::Comments;
48              
49              
50             # Not sure about this yet ... 2 or 4? With mirror images too 8 arms would
51             # fill the plane everywhere 4-visited points double-traversed segments.
52             # use constant parameter_info_array => [ { name => 'arms',
53             # share_key => 'arms_2',
54             # display => 'Arms',
55             # type => 'integer',
56             # minimum => 1,
57             # maximum => 2,
58             # default => 1,
59             # width => 1,
60             # description => 'Arms',
61             # } ];
62              
63 1     1   7 use constant n_start => 0;
  1         2  
  1         52  
64 1     1   5 use constant x_negative_at_n => 6;
  1         2  
  1         42  
65 1     1   7 use constant y_negative_at_n => 22;
  1         1  
  1         54  
66 1     1   5 use constant 1.02 _UNDOCUMENTED__dxdy_list_at_n => 7;
  1         13  
  1         1585  
67              
68              
69             #------------------------------------------------------------------------------
70              
71             sub new {
72 3     3 1 72 my $self = shift->SUPER::new(@_);
73 3   50     28 $self->{'arms'} = max(1, min(2, $self->{'arms'} || 1));
74 3         9 return $self;
75             }
76              
77              
78             sub n_to_xy {
79 1425     1425 1 27815 my ($self, $n) = @_;
80             ### CCurve n_to_xy(): $n
81              
82 1425 100       2806 if ($n < 0) { return; }
  1         3  
83 1424 50       2686 if (is_infinite($n)) { return ($n, $n); }
  0         0  
84              
85 1424         2502 my $zero = ($n * 0); # inherit bignum 0
86 1424         1923 my $x = $zero;
87 1424         2314 my $y = $zero;
88             {
89 1424         1757 my $int = int($n);
  1424         1785  
90 1424         1894 $x = $n - $int; # inherit possible BigFloat
91 1424         1990 $n = $int; # BigFloat int() gives BigInt, use that
92             }
93              
94             # initial rotation from arm number $n mod $arms
95 1424         2711 my $rot = _divrem_mutate ($n, $self->{'arms'});
96              
97 1424         2064 my $len = $zero+1;
98 1424         2619 foreach my $digit (digit_split_lowtohigh($n,4)) {
99             ### $digit
100              
101 6169 100       11806 if ($digit == 0) {
    100          
    100          
102 1194         2039 ($x,$y) = ($y,-$x); # rotate -90
103             } elsif ($digit == 1) {
104 1672         2149 $y -= $len; # at Y=-len
105             } elsif ($digit == 2) {
106 1654         2184 $x += $len; # at X=len,Y=-len
107 1654         2162 $y -= $len;
108             } else {
109             ### assert: $digit == 3
110 1649         2876 ($x,$y) = (2*$len - $y, # at X=2len,Y=-len and rotate +90
111             $x-$len);
112             }
113 6169         7670 $rot++; # to keep initial direction
114 6169         8495 $len *= 2;
115             }
116              
117 1424 100       2798 if ($rot & 2) {
118 224         321 $x = -$x;
119 224         295 $y = -$y;
120             }
121 1424 100       2398 if ($rot & 1) {
122 957         1700 ($x,$y) = (-$y,$x);
123             }
124              
125             ### final: "$x,$y"
126 1424         3160 return ($x,$y);
127             }
128              
129             # point N=2^(2k) at XorY=+/-2^k radius 2^k
130             # N=2^(2k-1) at X=Y=+/-2^(k-1) radius sqrt(2)*2^(k-1)
131             # radius = sqrt(2^level)
132             # R(l)-R(l-1) = sqrt(2^level) - sqrt(2^(level-1))
133             # = sqrt(2^level) * (1 - 1/sqrt(2))
134             # about 0.29289
135              
136             # len=1 extent of lower level 0
137             # len=4 extent of lower level 2
138             # len=8 extent of lower level 4+1 = 5
139             # len=16 extent of lower level 8+3
140             # len/2 + len/4-1
141              
142             my @digit_to_rot = (-1, 1, 0, 1);
143             my @dir4_to_dsdd = ([1,-1],[1,1],[-1,1],[-1,-1]);
144              
145             sub xy_to_n {
146 61     61 1 12168 return scalar((shift->xy_to_n_list(@_))[0]);
147             }
148             sub xy_to_n_list {
149 113     113 1 8392 my ($self, $x, $y) = @_;
150             ### CCurve xy_to_n(): "$x, $y"
151              
152 113         282 $x = round_nearest($x);
153 113         239 $y = round_nearest($y);
154 113         192 my $zero = $x*0*$y;
155              
156 113         224 ($x,$y) = ($x + $y, $y - $x); # sum and diff
157 113 50       225 if (is_infinite($x)) { return $x; }
  0         0  
158 113 50       250 if (is_infinite($y)) { return $y; }
  0         0  
159              
160 113         190 my @n_list;
161 113         203 foreach my $dsdd (@dir4_to_dsdd) {
162 452         779 my ($ds,$dd) = @$dsdd;
163             ### attempt: "ds=$ds dd=$dd"
164 452         608 my $s = $x; # sum X+Y
165 452         577 my $d = $y; # diff Y-X
166 452         588 my @nbits;
167              
168 452   100     1734 until ($s >= -1 && $s <= 1 && $d >= -1 && $d <= 1) {
      100        
      100        
169             ### at: "s=$s, d=$d nbits=".join('',reverse @nbits)
170 2852         1963 my $bit = $s % 2;
171 2852         1932 push @nbits, $bit;
172 2852 100       2292 if ($bit) {
173 510         647 $s -= $ds;
174 510         659 $d -= $dd;
175 510         797 ($ds,$dd) = ($dd,-$ds); # rotate -90
176             }
177              
178             # divide 1/(1+i) = (1-i)/(1^2 - i^2)
179             # = (1-i)/2
180             # so multiply (s + i*d) * (1-i)/2
181             # s = (s + d)/2
182             # d = (d - s)/2
183             #
184             ### assert: (($s+$d)%2)==0
185              
186             # this form avoids overflow near DBL_MAX
187 2852         1953 my $odd = $s % 2;
188 2852         1716 $s -= $odd;
189 2852         1693 $d -= $odd;
190 2852         1888 $s /= 2;
191 2852         1775 $d /= 2;
192 2852         12594 ($s,$d) = ($s+$d+$odd, $d-$s);
193             }
194              
195             # five final positions
196             # . 0,1 . ds,dd
197             # |
198             # -1,0--0,0--1,0
199             # |
200             # . 0,-1 .
201             #
202             ### end: "s=$s d=$d ds=$ds dd=$dd"
203              
204             # last step must be East dx=1,dy=0
205 452 100 100     1055 unless ($ds == 1 && $dd == -1) { next; }
  288         544  
206              
207 164 100 100     533 if ($s == $ds && $d == $dd) {
    100 66        
208 86         163 push @nbits, 1;
209             } elsif ($s != 0 || $d != 0) {
210 75         139 next;
211             }
212             # ended s=0,d=0 or s=ds,d=dd, found an N
213 89         252 push @n_list, digit_join_lowtohigh(\@nbits, 2, $zero);
214             ### found N: "$n_list[-1]"
215             }
216             ### @n_list
217 113         376 return sort {$a<=>$b} @n_list;
  39         131  
218             }
219              
220             # f = (1 - 1/sqrt(2) = .292
221             # 1/f = 3.41
222             # N = 2^level
223             # Rend = sqrt(2)^level
224             # Rmin = Rend / 2 maybe
225             # Rmin^2 = (2^level)/4
226             # N = 4 * Rmin^2
227             #
228             sub rect_to_n_range {
229 5     5 1 970 my ($self, $x1,$y1, $x2,$y2) = @_;
230             ### CCurve rect_to_n_range(): "$x1,$y1 $x2,$y2"
231              
232 5         18 $x1 = round_nearest ($x1);
233 5         13 $x2 = round_nearest ($x2);
234 5         10 $y1 = round_nearest ($y1);
235 5         11 $y2 = round_nearest ($y2);
236              
237 5 50       14 ($x1,$x2) = ($x2,$x1) if $x1 > $x2;
238 5 50       12 ($y1,$y2) = ($y2,$y1) if $y1 > $y2;
239              
240 5         13 my ($len,$level) = _rect_to_k ($x1,$y1, $x2,$y2);
241 5 50       14 if (is_infinite($level)) {
242 0         0 return (0, $level);
243             }
244 5         19 return (0, 4*$len*$len*$self->{'arms'} - 1);
245             }
246              
247             # N=16 is Y=4 away k=2
248             # N=64 is Y=-8+1=-7 away k=3
249             # N=256=4^4 is X=2^4=16-3=-7 away k=4
250             # dist = 2^k - (2^(k-2)-1)
251             # = 2^k - 2^(k-2) + 1
252             # = 4*2^(k-2) - 2^(k-2) + 1
253             # = 3*2^(k-2) + 1
254             # k=2 3*2^(2-2)+1=4 len=4^2=16
255             # k=3 3*2^(3-2)+1=7 len=4^3=64
256             # k=4 3*2^(4-2)+1=13
257             # 2^(k-2) = (dist-1)/3
258             # 2^k = (dist-1)*4/3
259             #
260             # up = 3*2^(k-2+1) + 1
261             # 2^(k+1) = (dist-1)*4/3
262             # 2^k = (dist-1)*2/3
263             #
264             # left = 3*2^(k-2+1) + 1
265             # 2^(k+1) = (dist-1)*4/3
266             # 2^k = (dist-1)*2/3
267             #
268             # down = 3*2^(k-2+1) + 1
269             # 2^(k+1) = (dist-1)*4/3
270             # 2^k = (dist-1)*2/3
271             #
272             # m=2 4*(2-1)/3=4/3=1
273             # m=4 4*(4-1)/3=4
274             sub _rect_to_k {
275 5     5   12 my ($x1,$y1, $x2,$y2) = @_;
276             ### _rect_to_k(): $x1,$y1
277              
278             {
279 5         26 my $m = max(abs($x1),abs($y1),abs($x2),abs($y2));
  5         21  
280 5 50       12 if ($m < 2) {
281 0         0 return (2, 1);
282             }
283 5 50       11 if ($m < 4) {
284 0         0 return (4, 2);
285             }
286             ### round_down: 4*($m-1)/3
287 5         21 my ($len, $k) = round_down_pow (4*($m-1)/3, 2);
288 5         16 return ($len, $k);
289             }
290              
291 0         0 my $len;
292 0         0 my $k = 0;
293              
294 0         0 my $offset = -1;
295 0         0 foreach my $m ($x2, $y2, -$x1, -$y1) {
296 0         0 $offset++;
297             ### $offset
298             ### $m
299 0 0       0 next if $m < 0;
300              
301 0         0 my ($len1, $k1);
302             # if ($m < 2) {
303             # $len1 = 1;
304             # $k1 = 0;
305             # } else {
306             # }
307              
308 0         0 ($len1, $k1) = round_down_pow (($m-1)/3, 2);
309 0 0       0 next if $k1 < $offset;
310 0         0 my $sub = ($offset-$k1) % 4;
311 0         0 $k1 -= $sub; # round down to k1 == offset mod 4
312              
313 0 0       0 if ($k1 > $k) {
314 0         0 $k = $k1;
315 0         0 $len = $len1 / 2**$sub;
316             }
317             }
318              
319             ### result: "k=$k len=$len"
320 0         0 return ($len, 2*$k);
321             }
322              
323              
324              
325             my @dir4_to_dx = (1,0,-1,0);
326             my @dir4_to_dy = (0,1,0,-1);
327              
328             # Not yet supporting "arms" parameter ...
329             sub n_to_dxdy {
330 2530     2530 1 49550 my ($self, $n) = @_;
331             ### n_to_dxdy(): $n
332              
333 2530 100       4445 if ($n < 0) { return; }
  3         18  
334 2527 50       4415 if (is_infinite($n)) { return ($n,$n); }
  0         0  
335 2527         4096 my $int = int($n);
336 2527         3406 $n -= $int; # $n fraction part
337              
338 2527         4451 my @digits = bit_split_lowtohigh($int);
339 2527   100     10986 my $dir = (sum(@digits)||0) & 3; # count of 1-bits
340 2527         3862 my $dx = $dir4_to_dx[$dir];
341 2527         3319 my $dy = $dir4_to_dy[$dir];
342              
343 2527 100       4059 if ($n) {
344             # apply fraction part $n
345              
346             # count low 1-bits is right turn of N+1, apply as dir-(turn-1) so decr $dir
347 8         17 while (shift @digits) {
348 8         15 $dir--;
349             }
350              
351             # this with turn=count-1 turn which is dir++ worked into swap and negate
352             # of dir4_to_dy parts
353 8         16 $dir &= 3;
354 8         21 $dx -= $n*($dir4_to_dy[$dir] + $dx); # with rot-90 instead of $dir+1
355 8         16 $dy += $n*($dir4_to_dx[$dir] - $dy);
356              
357             # this the equivalent with explicit dir++ for turn=count-1
358             # $dir++;
359             # $dir &= 3;
360             # $dx += $n*($dir4_to_dx[$dir] - $dx);
361             # $dy += $n*($dir4_to_dy[$dir] - $dy);
362             }
363              
364             ### result: "$dx, $dy"
365 2527         7178 return ($dx,$dy);
366             }
367              
368             #------------------------------------------------------------------------------
369             # k even
370             # S[h]
371             # ---------
372             # / \ Z[h-1]
373             # / \
374             # | | S[h-1]
375             # \ / Z[h-2]
376             # -- --
377             # Hb[k] = S[h] + 2*S[h-1] + S[h] + 2*(Z[h-1]/2 - Z[h-2]/2)
378             # + sqrt(2)*(2*Z[h-1]/2 + 2*Z[h-2]/2)
379             # = 2*S[h] + 2*S[h-1] + Z[h-1]-Z[h-2] + sqrt(2) * (Z[h-1] + Z[h-2])
380             # = 2*2^h + 2*2^(h-1) + 2*2^(h-1)-2 - (2*2^(h-2)-2) + sqrt(2) * (2*2^(h-1)-2 + 2*2^(h-2)-2)
381             # = 3*2^h + 2*2^(h-1)-2 - 2*2^(h-2) + 2 + sqrt(2) * (3*2^(h-1) - 4)
382             # = 3*2^h + 2^(h-1) + sqrt(2) * (3*2^(h-1) - 4)
383             # = 7*2^(h-1) + sqrt(2) * (3*2^(h-1) - 4)
384             # = 7*sqrt(2)^(2h-2) + sqrt(2) * (3*sqrt(2)^(2h-2) - 4)
385             # = 7*sqrt(2)^(k-2) + sqrt(2) * (3*sqrt(2)^(k-2) - 4)
386             # = 7*sqrt(2)^(k-2) + sqrt(2)*3*sqrt(2)^(k-2) - 4*sqrt(2)
387             # = 7*sqrt(2)^(k-2) + 3*sqrt(2)*sqrt(2)^(k-2) - 4*sqrt(2)
388             # = (7 + 3*sqrt(2))*sqrt(2)^(k-2) - 4*sqrt(2)
389             #
390             # S[2]=4
391             # 11--10--7,9--6---5 Z[1]=2 k=4 h=2
392             # | | |
393             # 13--12 8 4---3 4 + 2*2 + 4+(2-0) = 14
394             # | | S[1]=2 (2+0) = 2
395             # 14 2
396             # | |
397             # 15---16 0---1 Z[0] = 0
398             #
399              
400             # k odd
401             # S[h]
402             # ----
403             # Z[h-1] / \ middle Z[h]
404             # S[h-1] | \
405             # \ \
406             # | S[h]
407             # |
408             # \ / Z[h-1]
409             # --
410             # S[h-1]
411             #
412             # Hb[k] = 2*S[h] + 2*S[h-1] + sqrt(2)*( Z[h]/2 + Z[h-1] + Z[h]/2 + S[h]-S[h-1] )
413             # = 2*S[h] + 2*S[h-1] + sqrt(2)*( Z[h] + Z[h-1] + S[h]-S[h-1] )
414             # = 2*2^h + 2*2^(h-1) + sqrt(2)*( 2*2^h-2 + 2*2^(h-1)-2 + 2^h - 2^(h-1) )
415             # = 3*2^h + sqrt(2)*( 3*2^h + 2^(h-1) - 4 )
416             # = 3*2^h + sqrt(2)*( 7*2^(h-1) - 4 )
417              
418             sub _UNDOCUMENTED_level_to_hull_boundary {
419 0     0     my ($self, $level) = @_;
420 0 0         my ($a, $b) = $self->_UNDOCUMENTED_level_to_hull_boundary_sqrt2($level)
421             or return undef;
422 0           return $a + $b*sqrt(2);
423             }
424             sub _UNDOCUMENTED_level_to_hull_boundary_sqrt2 {
425 0     0     my ($self, $level) = @_;
426 0 0         if ($level <= 2) {
427 0 0         if ($level < 0) { return; }
  0            
428 0 0         if ($level == 2) { return (6,0); }
  0            
429 0 0         return (2, ($level == 0 ? 0 : 1));
430             }
431              
432 0           my ($h, $rem) = _divrem($level, 2);
433 0           my $pow = 2**($h-1);
434              
435 0 0         if ($rem) {
436 0           return (6*$pow, 7*$pow-4);
437              
438             # return (2*S_formula($h) + 2*S_formula($h-1),
439             # Z_formula($h)/2 + Z_formula($h-1)
440             # + Z_formula($h)/2 + (S_formula($h)-S_formula($h-1)) );
441              
442             } else {
443 0           return (7*$pow, 3*$pow-4);
444              
445             # return (S_formula($h) + 2*S_formula($h-1) + S_formula($h)+(Z_formula($h-1)-Z_formula($h-2)),
446             # (Z_formula($h-1) + Z_formula($h-2)));
447             }
448             }
449              
450             #------------------------------------------------------------------------------
451             {
452             my @_UNDOCUMENTED_level_to_hull_area = (0, 1/2, 2);
453              
454             sub _UNDOCUMENTED_level_to_hull_area {
455 0     0     my ($self, $level) = @_;
456              
457 0 0         if ($level < 3) {
458 0 0         if ($level < 0) { return undef; }
  0            
459 0           return $_UNDOCUMENTED_level_to_hull_area[$level];
460             }
461 0           my ($h, $rem) = _divrem($level, 2);
462 0 0         return 35*2**($level-4) - ($rem ? 13 : 10)*2**($h-1) + 2;
463              
464             # if ($rem) {
465             # return 35*2**($level-4) - 13*$pow + 2;
466             #
467             # my $width = S_formula($h) + Z_formula($h)/2 + Z_formula($h-1)/2;
468             # my $ul = Z_formula($h-1)/2;
469             # my $ur = Z_formula($h)/2;
470             # my $bl = $width - Z_formula($h-1)/2 - S_formula($h-1);
471             # my $br = Z_formula($h-1)/2;
472             # return $width**2 - $ul**2/2 - $ur**2/2 - $bl**2/2 - $br**2/2;
473             #
474             # } else {
475             # return 35*2**($level-4) - 10*$pow + 2;
476             # return 0;
477             # return 35*2**($level-4) - 5*2**$h + 2;
478             #
479             # # my $width = S_formula($h) + Z_formula($h-1);
480             # # my $upper = Z_formula($h-1)/2;
481             # # my $lower = Z_formula($h-2)/2;
482             # # my $height = S_formula($h-1) + $upper + $lower;
483             # # return $width; # * $height - $upper*$upper - $lower*$lower;
484             # }
485             # }
486             }
487             }
488              
489             #------------------------------------------------------------------------------
490             # levels
491              
492             sub level_to_n_range {
493 0     0 1   my ($self, $level) = @_;
494 0           return (0, 2**$level);
495             }
496             sub n_to_level {
497 0     0 1   my ($self, $n) = @_;
498 0 0         if ($n < 0) { return undef; }
  0            
499 0 0         if (is_infinite($n)) { return $n; }
  0            
500 0           $n = round_nearest($n);
501 0           my ($pow, $exp) = round_up_pow ($n, 2);
502 0           return $exp;
503             }
504              
505             #------------------------------------------------------------------------------
506             1;
507             __END__