Math-PlanePath
view release on metacpan or search on metacpan
lib/Math/PlanePath/CCurve.pm view on Meta::CPAN
# + sqrt(2)*(2*Z[h-1]/2 + 2*Z[h-2]/2)
# = 2*S[h] + 2*S[h-1] + Z[h-1]-Z[h-2] + sqrt(2) * (Z[h-1] + Z[h-2])
# = 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)
# = 3*2^h + 2*2^(h-1)-2 - 2*2^(h-2) + 2 + sqrt(2) * (3*2^(h-1) - 4)
# = 3*2^h + 2^(h-1) + sqrt(2) * (3*2^(h-1) - 4)
# = 7*2^(h-1) + sqrt(2) * (3*2^(h-1) - 4)
# = 7*sqrt(2)^(2h-2) + sqrt(2) * (3*sqrt(2)^(2h-2) - 4)
# = 7*sqrt(2)^(k-2) + sqrt(2) * (3*sqrt(2)^(k-2) - 4)
# = 7*sqrt(2)^(k-2) + sqrt(2)*3*sqrt(2)^(k-2) - 4*sqrt(2)
# = 7*sqrt(2)^(k-2) + 3*sqrt(2)*sqrt(2)^(k-2) - 4*sqrt(2)
# = (7 + 3*sqrt(2))*sqrt(2)^(k-2) - 4*sqrt(2)
#
# S[2]=4
# 11--10--7,9--6---5 Z[1]=2 k=4 h=2
# | | |
# 13--12 8 4---3 4 + 2*2 + 4+(2-0) = 14
# | | S[1]=2 (2+0) = 2
# 14 2
# | |
# 15---16 0---1 Z[0] = 0
#
# k odd
# S[h]
# ----
# Z[h-1] / \ middle Z[h]
# S[h-1] | \
# \ \
# | S[h]
# |
# \ / Z[h-1]
# --
# S[h-1]
#
# 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] )
# = 2*S[h] + 2*S[h-1] + sqrt(2)*( Z[h] + Z[h-1] + S[h]-S[h-1] )
# = 2*2^h + 2*2^(h-1) + sqrt(2)*( 2*2^h-2 + 2*2^(h-1)-2 + 2^h - 2^(h-1) )
# = 3*2^h + sqrt(2)*( 3*2^h + 2^(h-1) - 4 )
# = 3*2^h + sqrt(2)*( 7*2^(h-1) - 4 )
sub _UNDOCUMENTED_level_to_hull_boundary {
my ($self, $level) = @_;
my ($a, $b) = $self->_UNDOCUMENTED_level_to_hull_boundary_sqrt2($level)
or return undef;
return $a + $b*sqrt(2);
}
sub _UNDOCUMENTED_level_to_hull_boundary_sqrt2 {
my ($self, $level) = @_;
if ($level <= 2) {
if ($level < 0) { return; }
if ($level == 2) { return (6,0); }
return (2, ($level == 0 ? 0 : 1));
}
my ($h, $rem) = _divrem($level, 2);
my $pow = 2**($h-1);
if ($rem) {
return (6*$pow, 7*$pow-4);
# return (2*S_formula($h) + 2*S_formula($h-1),
# Z_formula($h)/2 + Z_formula($h-1)
# + Z_formula($h)/2 + (S_formula($h)-S_formula($h-1)) );
} else {
return (7*$pow, 3*$pow-4);
# return (S_formula($h) + 2*S_formula($h-1) + S_formula($h)+(Z_formula($h-1)-Z_formula($h-2)),
# (Z_formula($h-1) + Z_formula($h-2)));
}
}
#------------------------------------------------------------------------------
{
my @_UNDOCUMENTED_level_to_hull_area = (0, 1/2, 2);
sub _UNDOCUMENTED_level_to_hull_area {
my ($self, $level) = @_;
if ($level < 3) {
if ($level < 0) { return undef; }
return $_UNDOCUMENTED_level_to_hull_area[$level];
}
my ($h, $rem) = _divrem($level, 2);
return 35*2**($level-4) - ($rem ? 13 : 10)*2**($h-1) + 2;
# if ($rem) {
# return 35*2**($level-4) - 13*$pow + 2;
#
# my $width = S_formula($h) + Z_formula($h)/2 + Z_formula($h-1)/2;
# my $ul = Z_formula($h-1)/2;
# my $ur = Z_formula($h)/2;
# my $bl = $width - Z_formula($h-1)/2 - S_formula($h-1);
# my $br = Z_formula($h-1)/2;
# return $width**2 - $ul**2/2 - $ur**2/2 - $bl**2/2 - $br**2/2;
#
# } else {
# return 35*2**($level-4) - 10*$pow + 2;
# return 0;
# return 35*2**($level-4) - 5*2**$h + 2;
#
# # my $width = S_formula($h) + Z_formula($h-1);
# # my $upper = Z_formula($h-1)/2;
# # my $lower = Z_formula($h-2)/2;
# # my $height = S_formula($h-1) + $upper + $lower;
# # return $width; # * $height - $upper*$upper - $lower*$lower;
# }
# }
}
}
#------------------------------------------------------------------------------
# levels
sub level_to_n_range {
my ($self, $level) = @_;
return (0, 2**$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_up_pow ($n, 2);
return $exp;
}
#------------------------------------------------------------------------------
1;
__END__
=for stopwords eg Ryde Math-PlanePath ie OEIS dX,dY dX combinatorial Ramus th zig zags stairstep Duvall Keesling vy Preprint tilings optimization
=head1 NAME
Math::PlanePath::CCurve -- Levy C curve
=head1 SYNOPSIS
use Math::PlanePath::CCurve;
my $path = Math::PlanePath::CCurve->new;
my ($x, $y) = $path->n_to_xy (123);
=head1 DESCRIPTION
This is an integer version of the "C" curve by LE<233>vy.
=over
"Les Courbes Planes ou Gauches et les Surfaces ComposE<233>e de Parties
Semblables au Tout", Journal de l'E<201>cole Polytechnique, July 1938 pages
227-247 and October 1938 pages 249-292
L<http://gallica.bnf.fr/ark:/12148/bpt6k57344323/f53.image>,
L<http://gallica.bnf.fr/ark:/12148/bpt6k57344820/f3.image>
=back
It spirals anti-clockwise, variously crossing and overlapping itself. The
construction is straightforward but various measurements like how many
distinct points are quite complicated.
11-----10-----9,7-----6------5 3
| | |
13-----12 8 4------3 2
lib/Math/PlanePath/CCurve.pm view on Meta::CPAN
+ ...
Notice the i power is not the bit position k, but rather how many 1-bits are
above the position.
=head2 Level Ranges 4^k
The X,Y extents of the path through to Nlevel=2^k can be expressed as a
width and height measured relative to the endpoints.
*------------------* <-+
| | |
*--* *--* | height h[k]
| | |
* N=4^k N=0 * <-+
| | | | | below l[k]
*--*--* *--*--* <-+
^-----^ ^-----^
width 2^k width
w[k] w[k] Extents to N=4^k
<------------------------>
total width = 2^k + 2*w[k]
N=4^k is on either the X or Y axis and for the extents here it's taken
rotated as necessary to be horizontal. k=2 N=4^2=16 shown above is already
horizontal. The next level k=3 N=64=4^3 would be rotated -90 degrees to be
horizontal.
The width w[k] is measured from the N=0 and N=4^k endpoints. It doesn't
include the 2^k length between those endpoints. The two ends are symmetric
so the extent is the same at each end.
h[k] = 2^k - 1 0,1,3,7,15,31,etc
w[k] = / 0 for k=0
\ 2^(k-1) - 1 for k>=1 0,0,1,3,7,15,etc
l[k] = / 0 for k<=1
\ 2^(k-2) - 1 for k>=2 0,0,0,1,3,7,etc
The initial N=0 to N=64 shown above is k=3. h[3]=7 is the X=-7 horizontal.
l[3]=1 is the X=1 horizontal. w[3]=3 is the vertical Y=3, and also Y=-11
which is 3 below the endpoint N=64 at Y=8.
Expressed as a fraction of the 2^k distance between the endpoints the
extents approach total 2 wide by 1.25 high,
*------------------* <-+
| | | 1
*--* *--* | total
| | | height
* N=4^k N=0 * <-+ -> 1+1/4
| | | | | 1/4
*--*--* *--*--* <-+
^-----^ ^-----^
1/2 1 1/2 total width -> 2
The extent formulas can be found by considering the self-similar blocks.
The initial k=0 is a single line segment and all its extents are 0.
h[0] = 0
N=1 ----- N=0
l[0] = 0
w[0] = 0
Thereafter the replication overlap as
+-------+---+-------+
| | | |
+------+ | | +------+
| | D | | C | | B | | <-+
| +-------+---+-------+ | | 2^(k-1)
| | | | | previous
| | | | | level ends
| E | | A | <-+
+------+ +------+
^---------------^
2^k this level ends
w[k] = max (h[k-1], w[k-1]) # right of A,B
h[k] = 2^(k-1) + max (h[k-1], w[k-1]) # above B,C,D
l[k] = max w[k-1], l[k-1]-2^(k-1) # below A,E
Since h[k]=2^(k-1)+w[k] have S<h[k] E<gt> w[k]> for kE<gt>=1 and with the
initial h[0]=w[k]=0 have h[k]E<gt>=w[k] always. So the max of those two
is h.
h[k] = 2^(k-1) + h[k-1] giving h[k] = 2^k-1 for k>=1
w[k] = h[k-1] giving w[k] = 2^(k-1)-1 for k>=1
The max for l[k] is always w[k-1] as l[k] is never big enough that the parts
B-C and C-D can extend down past their 2^(k-1) vertical position.
(l[0]=w[0]=0 and thereafter by induction l[k]E<lt>=w[k].)
l[k] = w[k-1] giving l[k] = 2^(k-2)-1 for k>=2
=head2 Repeated Points
The curve crosses itself and can repeat X,Y positions up to 4 times. The
first double, triple and quadruple points are at
visits X,Y N
------ ------- ----------------------
2 -2, 3 7, 9
3 18, -7 189, 279, 281
4 -32, 55 1727, 1813, 2283, 2369
=cut
# binary
# 2 -10, 11 111, 1001
# 3 2
# 3 10010, -111 10111101, 100010111, 100011001
# 6 5 4
# 4 -100000, 110111 11010111111, 11100010101,
# 100011101011, 100101000001
# 9, 6, 7, 4
lib/Math/PlanePath/CCurve.pm view on Meta::CPAN
# The S<"_ _ _"> line shown which is part of the 24-pattern above but omitted
# here. This line is at Y=2^k. The extents described above mean that it
# extends down to Y=2^k - h[k] = 2^k-(2^k-1)=1, so it visits some points in
# row Y=1 and higher. Omitting the curve means there are YE<gt>=1 not visited
# 4 times. Similarly YE<lt>=-1 and XE<lt>-1 and XE<gt>=+1.
=pod
=head1 FUNCTIONS
See L<Math::PlanePath/FUNCTIONS> for the behaviour common to all path
classes.
=over 4
=item C<$path = Math::PlanePath::CCurve-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 = $path-E<gt>xy_to_n ($x,$y)>
Return the point number for coordinates C<$x,$y>. If there's nothing at
C<$x,$y> then return C<undef>. If C<$x,$y> is visited more than once then
return the smallest C<$n> which visits it.
=item C<@n_list = $path-E<gt>xy_to_n_list ($x,$y)>
Return a list of N point numbers at coordinates C<$x,$y>. If there's
nothing at C<$x,$y> then return an empty list.
A given C<$x,$y> is visited at most 4 times so the returned list is at most
4 values.
=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, 2**$level)>.
=back
=head1 FORMULAS
Some formulas and results can also be found in the author's mathematical
write-up
=over
L<http://user42.tuxfamily.org/c-curve/index.html>
=back
=head2 Direction
The direction or net turn of the curve is the count of 1 bits in N,
direction = count_1_bits(N) * 90degrees
For example N=11 is binary 1011 has three 1 bits, so direction 3*90=270
degrees, ie. to the south.
This bit count is because at each power-of-2 position the curve is a copy of
the lower bits but turned +90 degrees, so +90 for each 1-bit.
For powers-of-2 N=2,4,8,16, etc, there's only a single 1-bit so the
direction is always +90 degrees there, ie. always upwards.
=head2 Turn
At each point N the curve can turn in any direction: left, right, straight,
or 180 degrees back. The turn is given by the number of low 0-bits of N,
turn right = (count_low_0_bits(N) - 1) * 90degrees
For example N=8 is binary 0b100 which is 2 low 0-bits for turn=(2-1)*90=90
degrees to the right.
When N is odd there's no low zero bits and the turn is always (0-1)*90=-90
to the right, so every second turn is 90 degrees to the left.
=head2 Next Turn
The turn at the point following N, ie. at N+1, can be calculated by counting
the low 1-bits of N,
next turn right = (count_low_1_bits(N) - 1) * 90degrees
For example N=11 is binary 0b1011 which is 2 low one bits for
nextturn=(2-1)*90=90 degrees to the right at the following point, ie. at
N=12.
This works simply because low 1-bits like ..0111 increment to low 0-bits
..1000 to become N+1. The low 1-bits at N are thus the low 0-bits at N+1.
=head2 N to dX,dY
C<n_to_dxdy()> is implemented using the direction described above. For
integer N the count mod 4 gives the direction for dX,dY.
dir = count_1_bits(N) mod 4
dx = dir_to_dx[dir] # table 0 to 3
dy = dir_to_dy[dir]
For fractional N the direction at int(N)+1 can be obtained from the
direction at int(N) and the turn at int(N)+1, which is the low 1-bits of N
per L</Next Turn> above. Those two directions can then be combined as
described in L<Math::PlanePath/N to dX,dY -- Fractional>.
# apply turn to make direction at Nint+1
turn = count_low_1_bits(N) - 1 # N integer part
dir = (dir - turn) mod 4 # direction at N+1
# adjust dx,dy by fractional amount in this direction
dx += Nfrac * (dir_to_dx[dir] - dx)
dy += Nfrac * (dir_to_dy[dir] - dy)
A small optimization can be made by working the "-1" of the turn formula
into a +90 degree rotation of the C<dir_to_dx[]> and C<dir_to_dy[]> parts by
swap and sign change,
turn_plus_1 = count_low_1_bits(N) # on N integer part
dir = (dir - turn_plus_1) mod 4 # direction-1 at N+1
# adjustment including extra +90 degrees on dir
dx -= $n*(dir_to_dy[dir] + dx)
dy += $n*(dir_to_dx[dir] - dy)
=head2 X,Y to N
The N values at a given X,Y can be found by taking terms low to high from
the complex number formula (the same as given above)
X+iY = b^k N = 2^k + 2^(k1) + 2^(k2) + ... in binary
+ b^k1 * i base b=1+i
+ b^k2 * i^2
+ ...
If the lowest term is b^0 then X+iY has X+Y odd. If the lowest term is not
b^0 but instead some power b^n then X+iY has X+Y even. This is because a
multiple of b=1+i,
X+iY = (x+iy)*(1+i)
= (x-y) + (x+y)i
so X=x-y Y=x+y
sum X+Y = 2x is even if X+iY a multiple of 1+i
So the lowest bit of N is found by
bit = (X+Y) mod 2
If bit=1 then a power i^p is to be subtracted from X+iY. p is how many
1-bits are above that point, and this is not yet known. It represents a
direction to move X,Y to put it on an even position. It's also the
direction of the step N-2^l to N, where 2^l is the lowest 1-bit of N.
The reduction should be attempted with p commencing as each of the four
possible directions N,S,E,W. Some or all will lead to an N. For quadrupled
points (such as X=-32, Y=55 described above) all four will lead to an N.
for p 0 to 3
dX,dY = i^p # directions [1,0] [0,1] [-1,0] [0,-1]
loop until X,Y = [0,0] or [1,0] or [-1,0] or [0,1] or [0,-1]
{
bit = X+Y mod 2 # bits of N from low to high
if bit == 1 {
X -= dX # move to "even" X+Y == 0 mod 2
Y -= dY
(dX,dY) = (dY,-dX) # rotate -90 as for p-1
}
X,Y = (X+Y)/2, (Y-X)/2 # divide (X+iY)/(1+i)
}
if not (dX=1 and dY=0)
wrong final direction, try next p
if X=dX and Y=dY
further high 1-bit for N
found an N
if X=0 and Y=0
found an N
The "loop until" ends at one of the five points
0,1
|
-1,0 -- 0,0 -- 1,0
|
0,-1
It's not possible to wait for X=0,Y=0 to be reached because some dX,dY
directions will step infinitely among the four non-zeros. Only the case
X=dX,Y=dY is sure to reach 0,0.
The successive p decrements which rotate dX,dY by -90 degrees must end at p
== 0 mod 4 for highest term in the X+iY formula having i^0=1. This means
must end dX=1,dY=0 East. If this doesn't happen then there is no N for that
p direction.
The number of 1-bits in N is == p mod 4. So the order the N values are
obtained follows the order the p directions are attempted. In general the N
values will not be smallest to biggest N so a little sort is necessary if
that's desired.
It can be seen that sum X+Y is used for the bit calculation and then again
in the divide by 1+i. It's convenient to write the whole loop in terms of
sum S=X+Y and difference D=Y-X.
for dS = +1 or -1 # four directions
for dD = +1 or -1 #
S = X+Y
D = Y-X
loop until -1 <= S <= 1 and -1 <= D <= 1 {
bit = S mod 2 # bits of N from low to high
if bit == 1 {
S -= dS # move to "even" S+D == 0 mod 2
D -= dD
(dS,dD) = (dD,-dS) # rotate -90
}
(S,D) = (S+D)/2, (D-S)/2 # divide (S+iD)/(1+i)
}
if not (dS=1 and dD=-1)
wrong final direction, try next dS,dD direction
if S=dS and D=dD
further high 1-bit for N
found an N
if S=0 and D=0
found an N
The effect of S=X+Y, D=Y-D is to rotate by -45 degrees and use every second
point of the plane.
D= 2 X=0,Y=2 . rotate -45
D= 1 X=0,Y=1 . X=1,Y=2 .
D= 0 X=0,Y=0 . X=1,Y=1 . X=2,Y=2
D=-1 X=1,Y=0 . X=2,Y=1 .
D=-2 X=2,Y=0 .
S=0 S=1 S=2 S=3 S=4
The final five points described above are then in a 3x3 block at the origin.
The four in-between points S=0,D=1 etc don't occur so range tests
-1E<lt>=SE<lt>=1 and -1E<lt>=DE<lt>=1 can be used.
S=-1,D=1 . S=1,D=1
. S=0,D=0 .
S=-1,D=-1 . S=1,D=-1
lib/Math/PlanePath/CCurve.pm view on Meta::CPAN
# GP-Test vector(length(M0samples),k,M0half(k-1)) == M0samples
# GP-Test vector(length(M1samples),k,M1half(k-1)) == M1samples
# GP-Test vector(length(M2samples),k,M2half(k-1)) == M2samples
# GP-Test vector(length(M3samples),k,M3half(k-1)) == M3samples
=pod
The counts can be calculated in two ways. Firstly they satisfy mutual
recurrences. Each adds the preceding rotated M.
M0[k+1] = M0[k] + M3[k] initially M0[0] = 1 (N=0 to N=1)
M1[k+1] = M1[k] + M0[k] M1[0] = 0
M2[k+1] = M2[k] + M1[k] M2[0] = 0
M3[k+1] = M3[k] + M2[k] M3[0] = 0
Geometrically this can be seen from the way each level extends by a copy of
the previous level rotated +90,
7---6---5 Easts in N=0 to 8
| | = Easts in N=0 to 4
8 4---3 + Wests in N=0 to 4
| since N=4 to N=8 is
2 the N=0 to N=4 rotated +90
|
0---1
For the bits in N, level k+1 introduces a new bit either 0 or 1. In M0[k+1]
the a 0-bit is count M0[k] the same direction, and when a 1-bit is M3[k]
since one less bit mod 4. Similarly the other counts.
Some substitutions give 3rd order recurrences
for k >= 4
M0[k] = 4*M0[k-1] - 6*M0[k-2] + 4*M0[k-3] initial 1,1,1,1
M1[k] = 4*M1[k-1] - 6*M1[k-2] + 4*M1[k-3] initial 0,1,2,3
M2[k] = 4*M2[k-1] - 6*M2[k-2] + 4*M2[k-3] initial 0,0,1,3
M3[k] = 4*M3[k-1] - 6*M3[k-2] + 4*M3[k-3] initial 0,0,0,1
=cut
# GP-DEFINE M0rec(k) = if(k<4,1, 4*M0rec(k-1) - 6*M0rec(k-2) + 4*M0rec(k-3))
# GP-DEFINE M1rec(k) = if(k<4,k, 4*M1rec(k-1) - 6*M1rec(k-2) + 4*M1rec(k-3))
# GP-DEFINE M2rec(k) = if(k<2,0, if(k==2,1, if(k==3,3, 4*M2rec(k-1) - 6*M2rec(k-2) + 4*M2rec(k-3))))
# GP-DEFINE M3rec(k) = if(k<3,0, if(k==3,1, 4*M3rec(k-1) - 6*M3rec(k-2) + 4*M3rec(k-3)))
# GP-Test vector(20,k,M0rec(k-1)) == vector(20,k,M0half(k-1))
# GP-Test vector(20,k,M1rec(k-1)) == vector(20,k,M1half(k-1))
# GP-Test vector(20,k,M2rec(k-1)) == vector(20,k,M2half(k-1))
# GP-Test vector(20,k,M3rec(k-1)) == vector(20,k,M3half(k-1))
=pod
The characteristic polynomial of these recurrences is
x^3 - 4x^2 + 6x - 4
= (x-2) * (x - (1-i)) * (x - (1+i))
=for GP-Test x^3 - 4*x^2 + 6*x - 4 == (x-2)*(x^2 - 2*x + 2)
=for GP-Test x^3 - 4*x^2 + 6*x - 4 == (x-2) * (x + (I-1)) * (x - (I+1))
So explicit formulas can be written in powers of the roots 2, 1-i and 1+i,
M0[k] = ( 2^k + (1-i)^k + (1+i)^k )/4 for k>=1
M1[k] = ( 2^k + i*(1-i)^k - i*(1+i)^k )/4
M2[k] = ( 2^k - (1-i)^k - (1+i)^k )/4
M3[k] = ( 2^k - i*(1-i)^k + i*(1+i)^k )/4
=cut
# GP-DEFINE M0pow(k) = if(k==0,1, (1/4)*(2^k + (1-I)^k + (1+I)^k))
# GP-DEFINE M1pow(k) = if(k==0,0, (1/4)*(2^k + I*(1-I)^k - I*(1+I)^k))
# GP-DEFINE M2pow(k) = if(k==0,0, (1/4)*(2^k - (1-I)^k - (1+I)^k))
# GP-DEFINE M3pow(k) = if(k==0,0, (1/4)*(2^k - I*(1-I)^k + I*(1+I)^k))
# GP-Test vector(50,k,M0pow(k-1)) == vector(50,k,M0half(k-1))
# GP-Test vector(50,k,M1pow(k-1)) == vector(50,k,M1half(k-1))
# GP-Test vector(50,k,M2pow(k-1)) == vector(50,k,M2half(k-1))
# GP-Test vector(50,k,M3pow(k-1)) == vector(50,k,M3half(k-1))
=pod
The complex numbers 1-i and 1+i are 45 degree lines clockwise and
anti-clockwise respectively. The powers turn them in opposite directions so
the imaginary parts always cancel out. The remaining real parts can be had
by a half power h=floor(k/2) which is the magnitude abs(1-i)=sqrt(2)
projected onto the real axis. The sign selector d(n) above is whether the
positive or negative part of the real axis, or zero when at the origin.
The second way to calculate is the combinatorial interpretation that per
L</Direction> above the direction is count_1_bits(N) mod 4 so East segments
are all N values with count_1_bits(N) == 0 mod 4, ie. N with 0, 4, 8, etc
many 1-bits. The number of ways to have those bit counts within total k
bits is k choose 0, 4, 8 etc.
M0[k] = /k\ + /k\ + ... + / k\ m = floor(k/4)
\0/ \4/ \4m/
M1[k] = /k\ + /k\ + ... + / k \ m = floor((k-1)/4)
\1/ \5/ \4m+1/
M2[k] = /k\ + /k\ + ... + / k \ m = floor((k-2)/4)
\2/ \6/ \4m+2/
M3[k] = /k\ + /k\ + ... + / k \ m = floor((k-3)/4)
\3/ \7/ \4m+3/
=cut
# GP-DEFINE M0sum(k) = sum(i=0,floor(k/4), binomial(k, 4*i))
# GP-DEFINE M1sum(k) = sum(i=0,floor(k/4), binomial(k, 4*i+1))
# GP-DEFINE M2sum(k) = sum(i=0,floor(k/4), binomial(k, 4*i+2))
# GP-DEFINE M3sum(k) = sum(i=0,floor(k/4), binomial(k, 4*i+3))
# GP-Test vector(length(M0samples),k,M0sum(k-1)) == M0samples
# GP-Test vector(length(M1samples),k,M1sum(k-1)) == M1samples
# GP-Test vector(length(M2samples),k,M2sum(k-1)) == M2samples
# GP-Test vector(length(M3samples),k,M3sum(k-1)) == M3samples
=pod
The power forms above are cases of the identity by Ramus for sums of
binomial coefficients in arithmetic progression like this. (See Knuth
volume 1 section 1.2.6 exercise 30 for a form with cosines resulting from
( run in 0.958 second using v1.01-cache-2.11-cpan-df04353d9ac )