This file is indexed.

/usr/share/perl5/Math/PlanePath/KochCurve.pm is in libmath-planepath-perl 117-1.

This file is owned by root:root, with mode 0o644.

The actual contents of the file can be viewed below.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
# Copyright 2011, 2012, 2013, 2014 Kevin Ryde

# This file is part of Math-PlanePath.
#
# Math-PlanePath is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by the
# Free Software Foundation; either version 3, or (at your option) any later
# version.
#
# Math-PlanePath is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
# or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
# for more details.
#
# You should have received a copy of the GNU General Public License along
# with Math-PlanePath.  If not, see <http://www.gnu.org/licenses/>.


# math-image --path=KochCurve --lines --scale=10
# math-image --path=KochCurve --all --scale=10

# continuous but nowhere differentiable
#
# Sur une courbe continue sans tangente, obtenue par une construction
# géométrique élémentaire
#
# http://www.nku.edu/~curtin/grenouille.html
# http://www.nku.edu/~curtin/koch_171.jpg
#
# Cesàro, "Remarques sur la courbe de von Koch." Atti della
# R. Accad. della Scienze fisiche e matem. Napoli 12, No. 15, 1-12,
# 1905. Reprinted as §228 in Opere scelte, a cura dell'Unione matematica
# italiana e col contributo del Consiglio nazionale delle ricerche, Vol. 2:
# Geometria, analisi, fisica matematica. Rome: Edizioni Cremonese,
# pp. 464-479, 1964.
#
# Thue-Morse count 1s mod 2 is net direction
# Toeplitz first diffs is turn sequence +1 or -1
#
# J. Ma and J.A. Holdener. When Thue-Morse Meets Koch. In Fractals:
# Complex Geometry, Patterns, and Scaling in Nature and Society, volume 13,
# pages 191-206, 2005.
# http://personal.kenyon.edu/holdenerj/StudentResearch/WhenThueMorsemeetsKochJan222005.pdf
#
# F.M. Dekking. On the distribution of digits in arithmetic sequences. In
# Seminaire de Theorie des Nombres de Bordeaux, volume 12, pages 3201-3212,
# 1983.
#



package Math::PlanePath::KochCurve;
use 5.004;
use strict;
use List::Util 'sum','first';

use vars '$VERSION', '@ISA';
$VERSION = 117;
use Math::PlanePath;
@ISA = ('Math::PlanePath');
*_divrem_mutate = \&Math::PlanePath::_divrem_mutate;

use Math::PlanePath::Base::Generic
  'is_infinite',
  'round_nearest';
use Math::PlanePath::Base::Digits
  'round_down_pow',
  'digit_split_lowtohigh',
  'digit_join_lowtohigh';

# uncomment this to run the ### lines
# use Smart::Comments;

# undocumented previous name for this ...
*_round_down_pow = \&round_down_pow;


use constant n_start => 0;
use constant class_x_negative => 0;
use constant class_y_negative => 0;
use constant diffxy_minimum => 0;  # X>=Y octant so X-Y>=0
use constant dx_minimum => -2;
use constant dx_maximum => 2;
use constant dy_minimum => -1;
use constant dy_maximum => 1;
*_UNDOCUMENTED__dxdy_list = \&Math::PlanePath::_UNDOCUMENTED__dxdy_list_six;
use constant absdx_minimum => 1; # never vertical
use constant dsumxy_minimum => -2; # diagonals
use constant dsumxy_maximum => 2;
use constant ddiffxy_minimum => -2;
use constant ddiffxy_maximum => 2;
use constant dir_maximum_dxdy => (1,-1); # South-East


#------------------------------------------------------------------------------

sub n_to_xy {
  my ($self, $n) = @_;
  ### KochCurve n_to_xy(): $n

  # secret negatives to -.5
  if (2*$n < -1) { return; }
  if (is_infinite($n)) { return ($n,$n); }

  my $x;
  my $y;
  {
    my $int = int($n);
    $x = 2 * ($n - $int);  # usually positive, but n=-0.5 gives x=-0.5
    $y = $x * 0;           # inherit possible bigrat 0
    $n = $int;             # BigFloat int() gives BigInt, use that
  }

  my $len = $y+1;  # inherit bignum 1
  foreach my $digit (digit_split_lowtohigh($n,4)) {
    ### at: "$x,$y  digit=$digit"

    if ($digit == 0) {

    } elsif ($digit == 1) {
      ($x,$y) = (($x-3*$y)/2 + 2*$len,     # rotate +60
                 ($x+$y)/2);

    } elsif ($digit == 2) {
      ($x,$y) = (($x+3*$y)/2 + 3*$len,    # rotate -60
                 ($y-$x)/2   + $len);

    } else {
      ### assert: $digit==3
      $x += 4*$len;
    }
    $len *= 3;
  }

  ### final: "$x,$y"
  return ($x,$y);
}

sub xy_to_n {
  my ($self, $x, $y) = @_;
  ### KochPeaks xy_to_n(): "$x, $y"

  $x = round_nearest ($x);
  $y = round_nearest ($y);
  if ($y < 0 || $x < 0 || (($x ^ $y) & 1)) {
    ### neg y or parity different ...
    return undef;
  }
  my ($len,$level) = round_down_pow(($x/2)||1, 3);
  ### $level
  ### $len
  if (is_infinite($level)) {
    return $level;
  }

  my $n = 0;
  foreach (0 .. $level) {
    $n *= 4;
    ### at: "level=$level len=$len   x=$x,y=$y  n=$n"
    if ($x < 3*$len) {
      if ($x < 2*$len) {
        ### digit 0 ...
      } else {
        ### digit 1 ...
        $x -= 2*$len;
        ($x,$y) = (($x+3*$y)/2,   # rotate -60
                   ($y-$x)/2);
        $n += 1;
      }
    } else {
      $x -= 4*$len;
      ### digit 2 or 3 to: "x=$x"
      if ($x < $y) {   # before diagonal
        ### digit 2...
        $x += $len;
        $y -= $len;
        ($x,$y) = (($x-3*$y)/2,     # rotate +60
                   ($x+$y)/2);
        $n += 2;
      } else {
        #### digit 3...
        $n += 3;
      }
    }
    $len /= 3;
  }
  ### end at: "x=$x,y=$y   n=$n"
  if ($x != 0 || $y != 0) {
    return undef;
  }
  return $n;
}

# level extends to x= 2*3^level
#                  level = log3(x/2)
#
# exact
sub rect_to_n_range {
  my ($self, $x1,$y1, $x2,$y2) = @_;
  ### KochCurve rect_to_n_range(): "$x1,$y1  $x2,$y2"

  $x1 = round_nearest ($x1);
  $x2 = round_nearest ($x2);
  $y1 = round_nearest ($y1);
  $y2 = round_nearest ($y2);
  if ($x1 > $x2) { ($x1,$x2) = ($x2,$x1); }
  if ($y1 > $y2) { ($y1,$y2) = ($y2,$y1); }

  if ($x2 < 0 || $y2 < 0
      || 3*$y1 > $x2 ) {   # above line Y=X/3
    return (1,0);
  }

  #        \
  #          \
  #       *    \
  #      / \     \
  # o-+-*   *-+-e  \
  # 0     3     6
  #
  # 3*Y+X/2 - (Y!=0)
  #
  #                  /
  #             *-+-*
  #              \
  #       *       *
  #      / \     /
  # o-+-*   *-+-*
  # 0     3     6   X/2
  #
  my ($len, $level) = round_down_pow ($x2/2, 3);
  return _rect_to_n_range_rot ($len, $level, 0, $x1,$y1, $x2,$y2);



  # (undef, my $level) = round_down_pow ($x2/2, 3);
  # ### $level
  # return (0, 4**($level+1)-1);
}


my @dir6_to_dx = (2, 1,-1,-2, -1, 1);
my @dir6_to_dy = (0, 1, 1, 0, -1,-1);
my @max_digit_to_rot = (1, -2, 1, 0);
my @min_digit_to_rot = (0, 1, -2, 1);
my @max_digit_to_offset = (-1, -1, -1, 2);

sub _rect_to_n_range_rot {
  my ($initial_len, $level_max, $initial_rot, $x1,$y1, $x2,$y2) = @_;
  ### KochCurve _rect_to_n_range_rot(): "$x1,$y1  $x2,$y2  len=$initial_len level=$level_max rot=$initial_rot"

  my ($rot, $len, $x, $y);
  my $overlap = sub {
    ### overlap: "$x,$y len=$len rot=$rot"

    if ($len == 1) {
      return ($x >= $x1 && $x <= $x2
              && $y >= $y1 && $y <= $y2);
    }
    my $len = $len / 3;

    if ($rot < 3) {
      if ($rot == 0) {
        #       *
        #      / \
        # o-+-*   *-+-.
        return ($y <= $y2               # bottom before end
                && $y+$len >= $y1
                && $x <= $x2
                && $x+6*$len > $x1);    # right before end, exclusive
      } elsif ($rot == 1) {
        #       .
        #      /
        # *-+-*
        #  \
        #   *  +-----
        #  /   |x1,y2
        # o
        return ($x <= $x2              # left before end
                && $y+3*$len > $y1     # top after start, exclusive
                && $y-$x <= $y2-$x1);  # diag before corner
      } else {
        # .    |x1,y1
        #  \   +-----
        #   *
        #  /
        # *-+-*
        #      \
        #       o
        return ($y <= $y2              # bottom before end
                && $x-3*$len <=$x2     # left before end
                && $y+$x >= $y1+$x1);  # diag after corner
      }
    } else {
      if ($rot == 3) {
        # .-+-*   *-+-o
        #      \ /
        #       *
        return ($y >= $y1              # top after start
                && $y-$len <= $y2      # bottom before end
                && $x >= $x1           # right after start
                && $x-6*$len < $x2);   # left before end, exclusive
      } elsif ($rot == 4) {
        # x2,y1|    o
        # -----+   /
        #         *
        #          \
        #       *-+-*
        #      /
        #     .
        return ($x >= $x1              # right after start
                && $y-3*$len < $y2     # bottom before end, exclusive
                && $y-$x >= $y1-$x2);  # diag after corner
      } else {
        #    o
        #     \
        #      *-+-*
        #         /
        #        *
        # -----+  \
        # x2,y2|   .
        return ($y >= $y1              # top after start
                && $x+3*$len >= $x1    # right after start
                && $y+$x <= $y2+$x2);  # diag before corner
      }
    }
  };

  my $zero = 0*$x1*$x2*$y1*$y2;
  my @lens = ($initial_len);
  my $n_hi;
  $rot = $initial_rot;
  $len = $initial_len;
  $x = $zero;
  $y = $zero;
  my @digits = (4);

  for (;;) {
    my $digit = --$digits[-1];
    ### max at: "digits=".join(',',@digits)."  xy=$x,$y   len=$len"

    if ($digit < 0) {
      pop @digits;
      if (! @digits) {
        ### nothing found to level_max ...
        return (1, 0);
      }
      ### end of digits, backtrack ...
      $len = $lens[$#digits];
      next;
    }

    my $offset = $max_digit_to_offset[$digit];
    $rot = ($rot - $max_digit_to_rot[$digit]) % 6;
    $x += $dir6_to_dx[$rot] * $offset * $len;
    $y += $dir6_to_dy[$rot] * $offset * $len;

    ### $offset
    ### $rot

    if (&$overlap()) {
      if ($#digits >= $level_max) {
        ### yes overlap, found n_hi ...
        ### digits: join(',',@digits)
        ### n_hi: _digit_join_hightolow (\@digits, 4, $zero)
        $n_hi = _digit_join_hightolow (\@digits, 4, $zero);
        last;
      }
      ### yes overlap, descend ...
      push @digits, 4;
      $len = ($lens[$#digits] ||= $len/3);
    } else {
      ### no overlap, next digit ...
    }
  }

  $rot = $initial_rot;
  $x = $zero;
  $y = $zero;
  $len = $initial_len;
  @digits = (-1);

  for (;;) {
    my $digit = ++$digits[-1];
    ### min at: "digits=".join(',',@digits)."  xy=$x,$y   len=$len"

    if ($digit > 3) {
      pop @digits;
      if (! @digits) {
        ### oops, n_lo not found to level_max ...
        return (1, 0);
      }
      ### end of digits, backtrack ...
      $len = $lens[$#digits];
      next;
    }

    ### $digit
    ### rot increment: $min_digit_to_rot[$digit]
    $rot = ($rot + $min_digit_to_rot[$digit]) % 6;

    if (&$overlap()) {
      if ($#digits >= $level_max) {
        ### yes overlap, found n_lo ...
        ### digits: join(',',@digits)
        ### n_lo: _digit_join_hightolow (\@digits, 4, $zero)
        return (_digit_join_hightolow (\@digits, 4, $zero),
                $n_hi);
      }
      ### yes overlap, descend ...
      push @digits, -1;
      $len = ($lens[$#digits] ||= $len/3);

    } else {
      ### no overlap, next digit ...
      $x += $dir6_to_dx[$rot] * $len;
      $y += $dir6_to_dy[$rot] * $len;
    }
  }
}

# $aref->[0] high digit
sub _digit_join_hightolow {
  my ($aref, $radix, $zero) = @_;
  my @lowtohigh = reverse @$aref;
  return digit_join_lowtohigh(\@lowtohigh, $radix, $zero);
}


my @digit_to_dir = (0, 1, -1, 0);
my @digit_to_nextturn = (1,  # digit=1 (with +1 for "next" N)
                         -2, # digit=2
                         1); # digit=3
sub n_to_dxdy {
  my ($self, $n) = @_;
  ### n_to_dxdy(): $n

  if ($n < 0) {
    return;  # first direction at N=0
  }
  if (is_infinite($n)) {
    return ($n,$n);
  }

  my $int = int($n);
  $n -= $int;
  my @ndigits = digit_split_lowtohigh($int,4);

  my $dir6 = sum(0, map {$digit_to_dir[$_]} @ndigits) % 6;
  my $dx = $dir6_to_dx[$dir6];
  my $dy = $dir6_to_dy[$dir6];

  if ($n) {
    # fraction part

    # lowest non-3 digit, or zero if all 3s (0 above high digit)
    $dir6 += $digit_to_nextturn[ first {$_!=3} @ndigits, 0 ];
    $dir6 %= 6;
    $dx += $n*($dir6_to_dx[$dir6] - $dx);
    $dy += $n*($dir6_to_dy[$dir6] - $dy);
  }
  return ($dx, $dy);
}

sub _UNTESTED__n_to_dir6 {
  my ($self, $n) = @_;
  if ($n < 0) {
    return undef;  # first direction at N=0
  }
  if (is_infinite($n)) {
    return ($n,$n);
  }
  return (sum (map {$digit_to_dir[$_]} digit_split_lowtohigh($n,4))
          || 0) # if empty
    % 6;
}

my @n_to_turn6 = (undef,
                      1,  # +60 degrees
                      -2, # -120 degrees
                      1); # +60 degrees
sub _UNTESTED__n_to_turn6 {
  my ($self, $n) = @_;
  if (is_infinite($n)) {
    return undef;
  }
  while ($n) {
    my $digit = _divrem_mutate($n,4);
    if ($digit) {
      # lowest non-zero digit
      return $n_to_turn6[$digit];
    }
  }
  return 0;
}
sub _UNTESTED__n_to_turn_LSR {
  my ($self, $n) = @_;
  my $turn6 = $self->_UNTESTED__n_to_turn6($n) || return undef;
  return ($turn6 > 0 ? 1 : -1);
}
sub _UNTESTED__n_to_turn_left {
  my ($self, $n) = @_;
  my $turn6 = $self->_UNTESTED__n_to_turn6($n) || return undef;
  return ($turn6 > 0 ? 1 : 0);
}
sub _UNTESTED__n_to_turn_right {
  my ($self, $n) = @_;
  my $turn6 = $self->_UNTESTED__n_to_turn6($n) || return undef;
  return ($turn6 < 0 ? 1 : 0);
}

#------------------------------------------------------------------------------
# levels

sub level_to_n_range {
  my ($self, $level) = @_;
  return (0, 4**$level);
}
sub n_to_level {
  my ($self, $n) = @_;
  if ($n < 0) { return undef; }
  if (is_infinite($n)) { return $n; }
  $n = round_nearest($n);
  my ($pow, $exp) = round_down_pow ($n-1, 4);
  return $exp + 1;
}


#------------------------------------------------------------------------------
1;
__END__

=for stopwords eg Ryde Helge von Koch Math-PlanePath Nlevel differentiable ie OEIS Xlevel floorlevel Nhi Nlo Ndigit Une thode trique mentaire tude de Certaines orie des Courbes Acta Arithmetica

=head1 NAME

Math::PlanePath::KochCurve -- horizontal Koch curve

=head1 SYNOPSIS

 use Math::PlanePath::KochCurve;
 my $path = Math::PlanePath::KochCurve->new;
 my ($x, $y) = $path->n_to_xy (123);

=head1 DESCRIPTION

X<Koch, Helge von>This is an integer version of the self-similar Koch curve,

=over 4

Helge von Koch, "Une ME<233>thode GE<233>omE<233>trique
E<201>lE<233>mentaire pour l'E<201>tude de Certaines Questions de la
ThE<233>orie des Courbes Planes", Acta Arithmetica, volume 30, 1906, pages
145-174.  L<http://archive.org/details/actamathematica11lefgoog>

=back

It goes along the X axis and makes triangular excursions upwards.

                               8                                   3
                             /  \
                      6---- 7     9----10                18-...    2
                       \              /                    \
             2           5          11          14          17     1
           /  \        /              \        /  \        /
     0----1     3---- 4                12----13    15----16    <- Y=0

     ^
    X=0   2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19

The replicating shape is the initial N=0 to N=4,

            *
           / \
      *---*   *---*

which is rotated and repeated 3 times in the same pattern to give sections
N=4 to N=8, N=8 to N=12, and N=12 to N=16.  Then that N=0 to N=16 is itself
replicated three times at the angles of the base pattern, and so on
infinitely.

The X,Y coordinates are arranged on a square grid using every second point,
per L<Math::PlanePath/Triangular Lattice>.  The result is flattened
triangular segments with diagonals at a 45 degree angle.

=head2 Level Ranges

Each replication adds 3 copies of the existing points and is thus 4 times
bigger, so if N=0 to N=4 is reckoned as level 1 then a given replication
level goes from

    Nstart = 0
    Nlevel = 4^level   (inclusive)

Each replication is 3 times the width.  The initial N=0 to N=4 figure is 6
wide and in general a level runs from

    Xstart = 0
    Xlevel = 2*3^level   at N=Nlevel

The highest Y is 3 times greater at each level similarly.  The peak is at
the midpoint of each level,

    Npeak = (4^level)/2
    Ypeak = 3^level
    Xpeak = 3^level

It can be seen that the N=6 point backtracks horizontally to the same X as
the start of its section N=4 to N=8.  This happens in the further
replications too and is the maximum extent of the backtracking.

The Nlevel is multiplied by 4 to get the end of the next higher level.  The
same 4*N can be applied to all points N=0 to N=Nlevel to get the same shape
but a factor of 3 bigger X,Y coordinates.  The in-between points 4*N+1,
4*N+2 and 4*N+3 are then new finer structure in the higher level.

=head2 Fractal

Koch conceived the curve as having a fixed length and infinitely fine
structure, making it continuous everywhere but differentiable nowhere.  The
code here can be pressed into use for that sort of construction for a given
level of granularity by scaling

    X/3^level
    Y/3^level

which makes it a fixed 2 wide by 1 high.  Or for unit-side equilateral
triangles then apply further factors 1/2 and sqrt(3)/2, as noted in
L<Math::PlanePath/Triangular Lattice>.

    (X/2) / 3^level
    (Y*sqrt(3)/2) / 3^level


=head2 Area

The area under the curve to a given level can be calculated from its
self-similar nature.  The curve at level+1 is 3 times wider and higher and
adds a triangle of unit area onto each line segment.  So reckoning the line
segment N=0 to N=1 as level=0 (which is area[0]=0),

    area[level] = 9*area[level-1] + 4^(level-1)
                = 4^(level-1) + 9*4^(level-2) + ... + 9^(level-1)*4^0

                  9^level - 4^level
                = -----------------
                          5

                = 0, 1, 13, 133, 1261, 11605, 105469, ...  (A016153)

The sides are 6 different angles.  The triangles added on the sides are
always the same shape either pointing up or down.  Base width=2 and height=1
gives area=1.

       *            *-----*   ^
      / \            \   /    | height=1
     /   \            \ /     |
    *-----*            *      v     triangle area = 2*1/2 = 1

    <-----> width=2

If the Y coordinates are stretched to make equilateral triangles then the
number of triangles is not changed and so the area increases by a factor of
the area of the equilateral triangle, sqrt(3)/4.

=head1 FUNCTIONS

See L<Math::PlanePath/FUNCTIONS> for behaviour common to all path classes.

=over 4

=item C<$path = Math::PlanePath::KochCurve-E<gt>new ()>

Create and return a new path object.

=item C<($x,$y) = $path-E<gt>n_to_xy ($n)>

Return the X,Y coordinates of point number C<$n> on the path.  Points begin
at 0 and if C<$n E<lt> 0> then the return is an empty list.

Fractional positions give an X,Y position along a straight line between the
integer positions.

=item C<($n_lo, $n_hi) = $path-E<gt>rect_to_n_range ($x1,$y1, $x2,$y2)>

The returned range is exact, meaning C<$n_lo> and C<$n_hi> are the smallest
and biggest in the rectangle.

=item C<$n = $path-E<gt>n_start()>

Return 0, the first N in the path.

=back

=head2 Level Methods

=over

=item C<($n_lo, $n_hi) = $path-E<gt>level_to_n_range($level)>

Return C<(0, 4**$level)>.

=back

=head1 FORMULAS

=head2 N to Turn

The curve always turns either +60 degrees or -120 degrees, it never goes
straight ahead.  In the base 4 representation of N the lowest non-zero digit
gives the turn.  The first turn is at N=1 so there's always a non-zero digit
in N.

   low digit
    base 4         turn
   ---------   ------------
      1         +60 degrees (left)
      2        -120 degrees (right)
      3         +60 degrees (left)

For example N=8 is 20 base 4, so lowest nonzero "2" means turn -120 degrees
for the next segment.

If the least significant digit is non-zero then it determines the turn,
making the base N=0 to N=4 shape.  If the least significant is zero then the
next level up is in control, eg. N=0,4,8,12,16, making a turn according to
the base shape again at that higher level.  The first and last segments of
the base shape are "straight" so there's no extra adjustment to apply in
those higher digits.

This base 4 digit rule is equivalent to counting low 0-bits.  A low base-4
digit 1 or 3 is an even number of low 0-bits and a low digit 2 is an odd
number of low 0-bits.

    count low 0-bits         turn
    ----------------     ------------
         even             +60 degrees (left)
         odd             -120 degrees (right)

For example N=8 in binary "1000" has 3 low 0-bits and 3 is odd so turn -120
degrees (right).

See L<Math::PlanePath::GrayCode/Turn> for a similar turn sequence arising
from binary Gray code.

=head2 N to Next Turn

The turn at N+1, ie the next turn, can be found from the base-4 digits by
considering how the digits of N change when 1 is added, and the low-digit
turn calculation is applied on those changed digits.

Adding 1 means low digit 0, 1 or 2 will become non-zero.  Any low 3s wrap
around to become low 0s.  So the turn at N+1 can be found from the digits of
N by seeking the lowest non-3

   lowest non-3       turn
    digit of N       at N+1
   ------------   ------------
        0          +60 degrees (left)
        1         -120 degrees (right)
        2          +60 degrees (left)

=head2 N to Direction

The total turn at a given N can be found by counting digits 1 and 2 in
base 4.

    direction = ((count of 1-digits in base 4)
                 - (count of 2-digits in base 4)) * 60 degrees

For example N=11 is "23" in base 4, so 60*(0-1) = -60 degrees.

In this formula the count of 1s and 2s can go past 360 degrees, representing
a spiralling around which occurs at progressively higher replication levels.
The direction can be taken mod 360 degrees, or the count mod 6, for a
direction 0 to 5 if desired.

=head2 N to abs(dX),abs(dY)

The direction expressed as abs(dX) and abs(dY) can be calculated simply from
N modulo 3.  abs(dX) is a repeating pattern 2,1,1 and abs(dY) repeating
0,1,1.

    N mod 3     abs(dX),abs(dY)
    -------     ---------------
       0             2,0            horizontal, East or West
       1             1,1            slope North-East or South-West
       2             1,1            slope North-West or South-East

This works because the direction calculation above corresponds to N mod 3.
Each N digit in base 4 becomes

    N digit
    base 4    direction add
    -------   -------------
       0            0
       1            1
       2           -1
       3            0

Notice that direction == Ndigit mod 3.  Then because 4==1 mod 3 the
power-of-4 for each digit reduces down to 1,

    N = 4^k * digit_k + ... 4^0 * digit_0
    N mod 3 = 1 * digit_k + ... 1 * digit_0
            = digit_k + ... digit_0
    same as
    direction = digit_k + ... + digit_0    taken mod 3

=head2 Rectangle to N Range -- Level

An easy over-estimate of the N values in a rectangle can be had from the
Xlevel formula above.  If XlevelE<gt>rectangleX then Nlevel is past the
rectangle extent.

    X = 2*3^level

so

    floorlevel = floor log_base_3(X/2)
    Nhi = 4^(floorlevel+1) - 1

For example a rectangle extending to X=13 has floorlevel =
floor(log3(13/2))=1 and so Nhi=4^(1+1)-1=15.

The rounding-down of the log3 ensures a point such as X=18 which is the
first in the next Nlevel will give that next level.  So
floorlevel=log3(18/2)=2 (exactly) and Nhi=4^(2+1)-1=63.

The worst case for this over-estimate is when rectangleX==Xlevel, ie. the
rectangle is just into the next level.  In that case Nhi is nearly a factor
4 bigger than it needs to be.

=head2 Rectangle to N Range -- Exact

The exact Nlo and Nhi in a rectangle can be found by searching along the
curve.  For Nlo search forward from the origin N=0.  For Nhi search backward
from the Nlevel over-estimate described above.

At a given digit position in the prospective N the sub-part of the curve
comprising the lower digits has a certain triangular extent.  If it's
outside the target rectangle then step to the next digit value, and to the
next of the digit above when past digit=3 (or below digit=0 when searching
backwards).

There's six possible orientations for the curve sub-part.  In the following
pictures "o" is the start and the surrounding lines show the triangular
extent.  There's just four curve parts shown in each, but these triangles
bound a sub-curve of any level.

   rot=0   -+-               +-----------------+
         --   --              - .-+-*   *-+-o -
       --   *   --             --    \ /    --
     --    / \    --             --   *   --
    - o-+-*   *-+-. -              --   --
   +-----------------+       rot=3   -+-

   rot=1
   +---------+               rot=4    /+
   |      . /                        / |
   |     / /                        / o|
   |*-+-* /                        / / |
   | \   /                        / *  |
   |  * /                        /   \ |
   | / /                        / *-+-*|
   |o /                        / /     |
   | /                        / .      |
   +/                        +---------+

   +\  rot=2                 +---------+
   | \                        \ o      |
   |. \                        \ \     |
   | \ \                        \ *-+-*|
   |  * \                        \   / |
   | /   \                        \ *  |
   |*-+-* \                        \ \ |
   |     \ \                        \ .|
   |      o \                rot=5   \ |
   +---------+                        \+

The "." is the start of the next sub-curve.  It belongs to the next digit
value and so can be excluded.  For rot=0 and rot=3 this means simply
shortening the X range permitted.  For rot=1 and rot=4 similarly the Y
range.  For rot=2 and rot=5 it would require a separate test.

Tight sub-part extent checking reduces the sub-parts which are examined, but
it works perfectly well with a looser check, such as a square box for the
sub-curve extents.  Doing that might be easier if the target region is not a
rectangle but instead some trickier shape.

=head1 OEIS

The Koch curve is in Sloane's Online Encyclopedia of Integer Sequences in
various forms,

=over

L<http://oeis.org/A035263> (etc)

=back

    A177702   abs(dX) from N=1 onwards, being 1,1,2 repeating
    A011655   abs(dY), being 0,1,1 repeating
    A035263   turn 1=left,0=right, by morphism
    A096268   turn 0=left,1=right, by morphism
    A056832   turn 1=left,2=right, by replicate and flip last
    A029883   turn +/-1=left,0=right, Thue-Morse first differences
    A089045   turn +/-1=left,0=right, by +/- something

    A003159   N positions of left turns, ending even number 0 bits
    A036554   N positions of right turns, ending odd number 0 bits

    A020988   number of left turns N=0 to N < 4^k, being 2*(4^k-1)/3
    A002450   number of right turns N=0 to N < 4^k, being (4^k-1)/3
    A016153   area under the curve, (9^n-4^n)/5

For reference, A217586 is not quite the same as A096268 right turn.  A217586
differs by a 0E<lt>-E<gt>1 flip at N=2^k due to different initial a(1)=1.

=head1 SEE ALSO

L<Math::PlanePath>,
L<Math::PlanePath::PeanoCurve>,
L<Math::PlanePath::HilbertCurve>,
L<Math::PlanePath::KochPeaks>,
L<Math::PlanePath::KochSnowflakes>,
L<Math::PlanePath::CCurve>

L<Math::Fractal::Curve>

=head1 HOME PAGE

L<http://user42.tuxfamily.org/math-planepath/index.html>

=head1 LICENSE

Copyright 2011, 2012, 2013, 2014 Kevin Ryde

Math-PlanePath is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the Free
Software Foundation; either version 3, or (at your option) any later
version.

Math-PlanePath is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for
more details.

You should have received a copy of the GNU General Public License along with
Math-PlanePath.  If not, see <http://www.gnu.org/licenses/>.

=cut