\\ Copyright 2021 Kevin Ryde
\\
\\ This file 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.
\\
\\ This file 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 this file. If not, see .
\\ Usage: gp =0.
\\
\\ The most efficient (so far) single-term function is b_by_LtoH(m).
\\ But its LtoH_state_table[] is obscure until working up to how "remainder"
\\ bits can be accumulated in a state while reducing m.
default(strictargs,1);
default(recover,0);
\\-----------------------------------------------------------------------------
\\ Morphism
\\
\\ The morphism used by A073057 expands by replacing each term
\\
\\ 1 -> 1,2
\\ 2 -> 1,3
\\ 3 -> 4,2
\\ 4 -> 4,3
\\
morphism_table = {[ [1,2], \\ 1
[1,3], \\ 2
[4,2], \\ 3
[4,3] ]}; \\ 4
\\ v is a vector of terms, return it expanded "reps" many times
morphism_expand_vector(v,reps=1) =
{
for(i=1,reps,
v=concat(vecextract(morphism_table,v)));
v;
}
\\ s is an integer 1..4 and bit is 0 or 1.
\\ Expand s by the morphism and return the first or second term
\\ according as bit=0 or bit=1 respectively.
morphism_bit(s,bit) = morphism_table[s][bit+1];
{
print1("Morphism: ");
for(s=1,4,
printf(" %d->%d,%d ",
s, morphism_bit(s,0), morphism_bit(s,1)));
}
\\-----------------------------------------------------------------------------
\\ A073057 by Appending
\\
\\ The following a_vector() is a direct implementation of the definition:
\\ append the morphism expansion of the sequence so-far.
\\ a_vector(k) returns a vector of sequence terms after k expansions.
\\ k=0 is the initial [1,2,3,4].
\\ The vector returned has v[n] = A073057(n).
a_vector(k) =
{
my(v=[1,2,3,4]);
for(i=1,k, v=concat(v, morphism_expand_vector(v)));
v;
}
{
print("Expand and append: ");
for(k=0,2,
my(v=a_vector(k));
print1(" k=",k,": ");
for(n=1,#v, print1(v[n], if(n==#v,"",",")));
print());
}
\\ The length of each expansion is 2* the original and is appended, so new
\\ length 3* each, starting from len(0) = 4 initial terms.
\\
\\ len(k) = 4*3^k (A003946 after its first "1")
\\
len(k) = 4*3^k;
for(k=0,6, #a_vector(k)==len(k) || error());
{
print1("len(k) = ");
for(k=0,5, print1(len(k),", "));
print("...");
print();
}
\\-----------------------------------------------------------------------------
\\ Morphism Blocks
\\
\\ Another way to think of the appending is in terms of "Mblocks" (morphism
\\ blocks). Take
\\
\\ Mblock(0) = 1,2,3,4 initial terms
\\
\\ Then Mblock(p) is those initial terms expanded by the morphism p times.
\\ Each expansion doubles the length so Mlen(p) = 4*2^p.
\\ return a vector of terms 1,2,3,4 expanded p many times
Mblock(p) = morphism_expand_vector([1,2,3,4],p);
\\ length of Mblock(p) vector
Mlen(p) = 4*2^p;
for(p=0,6, #Mblock(p)==Mlen(p) || error());
{
print("Mblocks: (morphism expansions of 1,2,3,4)");
for(p=0,3,
print1(" Mblock("p") = ");
my(v=Mblock(p));
for(i=1,#v,print1(v[i],if(i==#v,"",",")));
print());
print();
}
\\ Each Mblock(p) in the sequence so-far expands to Mblock(p+1).
\\
\\ k=1 so far ---> to k=2 by append
\\ Mblock(0), Mblock(1) Mblock(0+1), Mblock(1+1)
\\ /------------------------\ /-----------------------------\
\\ Mblock(0) Mblock(1) Mblock(1) Mblock(2)
\\ 1,2,3,4 1,2,1,3,4,2,4,3 1,2,1,3,4,2,4,3 1,2,1,3,1,...
\\ j=0 j=1 j=2 j=3
\\
\\ The number of Mblocks doubles each time. Number them starting from j=0.
\\ Then block j is Mblock(p) expansion level
\\
\\ p = hammingweight(j)
\\
\\ j = 0 1 2 3 4 5 6 7 8 15
\\ p = 0,1, 1,2, 1,2,2,3, 1,2,2,3,2,3,3,4, ...
\\
\\ k=0 ^
\\ k=1 |-|
\\ k=2 |------|
\\ k=3 |---------------|
\\ k=4 |--------------------------------| 2^4-1 = 16 blocks 0..15
\\
\\ Each expansion level k comprises 2^k Mblocks j = 0 .. 2^k-1 inclusive.
\\ Each append puts a new highest 1-bit on each j, and each Mblock expands
\\ to p+1, hence p = hammingweight(j).
\\
\\ a_vector() is enough for calculating a vector of terms, but this Mblock
\\ structure shows sub-ranges which are useful in calculations below.
\\-----------------------------------------------------------------------------
\\ Recurrence
\\
\\ In any sort of block replication sequence, one way to calculate an
\\ individual a(n) is to consider where it replicated from, and any
\\ transformation applied.
\\
\\ Here it's convenient to work with an index m starting m=0. Take this as
\\ sequence b(m) with b(m) = a(m+1).
\\
\\ Sequence level k=0 expanded and appended gives k=1 as follows
\\
\\ k=0 +--- m=len(1) +--- m=len(2)
\\ /-----\ v v
\\ m = 0 1 2 3 4 5 6 7 8 9 10 11 12
\\ b(m) = 1,2,3,4, 1,2, 1,3, 4,2, 4,3 ...
\\ from m 0 0 1 1 2 2 3 3
\\ k=0 \-----/
\\ k=1 \----------------------------/
\\ k=2 \-------------------------------- ... --/
\\
\\ Terms m=4,5 expanded from m=0; terms m=6,7 expanded from m=1; and so on.
\\ In general, for len(k) <= m < len(k+1), how far m is after len(k)
\\ determines the m' where m expanded from, and a "bit" 0 or 1 for whether
\\ to take the first or second term of expansion of term b(m'),
\\
\\ m = len(k) + 2*m' + bit, where m in range len(k) <= m < len(k+1)
\\
\\ / m-len(k) \
\\ m' = floor| -------- | and bit = (m-len(k)) mod 2
\\ \ 2 /
\\
\\ Since len(k) is even, "bit" is unchanged by subtracting len(k), so also
\\
\\ bit = m mod 2
\\
\\ return k for which len(k) <= m < len(k+1), for m >= 4
m_came_from_k(m) = logint(m>>2,3);
\\ return a 2-element vector [mdash,bit] which is where m expanded from
m_came_from_and_bit(m) =
{
m>=4 || error();
my(k=m_came_from_k(m));
m -= len(k);
[m>>1, bitand(m,1)];
}
b_by_recurrence(m) =
{
m >= 0 || error();
if(m<4, return(m+1)); \\ initial m=0..3 values b(m) = 1..4
my(mdash,bit); [mdash,bit] = m_came_from_and_bit(m);
morphism_bit(b_by_recurrence(mdash), bit);
}
{
\\ check b_by_recurrence(m) against a_vector()
my(v=a_vector(7));
vector(#v,m,m--; b_by_recurrence(m)) == v || error();
}
{
print("Recurrence");
my(limit=21);
print1(" m = 0 1 2 3 ");
for(m=4,limit, printf("%2d ", m));
print("...");
print1(" from ");
for(m=4,limit,
my(mdash,bit); [mdash,bit] = m_came_from_and_bit(m);
printf("%2d ", mdash));
print("...");
print1(" bit ");
for(m=4,limit,
my(mdash,bit); [mdash,bit] = m_came_from_and_bit(m);
printf("%2d ", bit));
print("...");
print();
}
\\-----------------------------------------------------------------------------
\\ Offset into Mblock
\\
\\ Another approach for an individual term is to determine the Mblock(p)
\\ which m falls within, and the distance r into that block.
\\
\\ This is implicit in b_by_recurrence(). Successive "m_came_from" track
\\ back through the expansions of the block containing m until initial
\\ 1,2,3,4.
\\
\\ A slow way to find the Mblock start is just search up by block number j,
\\ knowing block levels p=hammingweight(j) and lengths Mlen(p).
m_to_Mpr_by_slow_search(m) =
{
my(j=0);
while(1,
m>=0 || error();
my(p=hammingweight(j),
L=Mlen(p));
if(m < L,
return([p,m]));
m -= L; \\ successively reduced, so will end as remainder r
j++);
}
\\ Having found the block, the block has length Mlen(p) and the distance of
\\ m into it is
\\
\\ 0 <= r < Mlen(p) = 2^(p+2)
\\
\\ For example Mblock p=0 has,
\\
\\ r = 0 1 2 3 offset into block
\\ Mblock(0) = 1,2,3,4
\\
\\ Or Mblock p=1 has,
\\
\\ r = 0 1 2 3 4 5 6 7 offset into block
\\ Mblock(1) = 1,2,1,3,4,2,4,3
\\
\\ Taking r as exactly p+2 bits (with high 0s as necessary), the highest
\\ 2 bits decide the initial 1,2,3,4 and then the rest from high to low are
\\ morphism expansions.
\\
\\ high low
\\ p+1 p p-1 ... 0 bit position
\\ r = h h x x
\\ \-----/ \----------/
\\ initial morphism
\\ 1,2,3,4 steps
\\ return the r'th term of Mblock(p), where first term is at r=0
Mblock_term(p,r) =
{
(0 <= r && r < Mlen(p)) || error();
my(s = r>>p + 1); \\ highest 2 bits for initial 1,2,3,4
forstep(i=p-1,0,-1, \\ then other bits high to low
s = morphism_bit(s,bittest(r,i)));
s;
}
b_by_Mpr_slow_search(m) =
{
my(p,r); [p,r] = m_to_Mpr_by_slow_search(m);
Mblock_term(p,r);
}
{
\\ check b_by_Mpr_slow_search(m) against a_vector()
my(v=a_vector(5));
vector(#v,m,m--; b_by_Mpr_slow_search(m)) == v || error();
}
\\ High to low by bits of r is the nature of any morphism which has
\\ fixed-width expansions (a uniform morphism).
\\
\\ For comparison, b_by_recurrence() holds bits of r on the call stack, in
\\ each "bit" local variable, waiting to be applied with morphism_bit().
\\ b_by_Mpr_slow_search() instead calculates all r in an integer.
\\-----------------------------------------------------------------------------
\\ Powers
\\
\\ In b_by_Mpr_slow_search(), the successive Mblock lengths are 4* Gould's
\\ sequence A001316 = 1,2,2,4,2,4,4,8,... = 2^hammingweight(j). Summing
\\ terms of that sequence is A006046 and 4* gives m = Mstartpos(j) which is
\\ the start location m of the j'th Mblock in the sequence.
\\
\\ M. F. Hasler's recurrence in A006046 shows the nature of the summing.
\\ Expanding this recurrence repeatedly is
\\
\\ A006046(n) = 3^k[1] + 2*3^k[2] + ... + 2^(p-1)*3^k[p]
\\
\\ where n = 2^k[1] + 2^k[2] + ... + 2^k[p]
\\ is the binary representation of n,
\\ with descending exponents k[1] > k[2] > ... > k[p]
\\
\\ Or equivalently, replace the 1-bits of n with successive powers of 2 and
\\ interpret the result as ternary.
\\
\\ n = 83 = 1 0 1 0 0 1 1 binary
\\ A006046(n) = 1 0 2 0 0 4 8 ternary
\\
\\ These are not ordinary ternary digits 0,1,2 of course, but are to be
\\ treated like digits with usual factors 3^position. PARI fromdigits()
\\ allows such larger-than-normal digits.
A006046(n) =
{
my(v=binary(n),t=1);
for(i=1,#v, if(v[i], v[i]=t; t<<=1));
fromdigits(v,3);
}
{
vector(20,n,n--; A006046(n)) \\ sample data from A006046
== [0,1,3,5,9,11,15,19,27,29,33,37,45,49,57,65,81,83,87,91] || error();
my(n=83);
binary(n) == [1,0,1, 0,0,1,1] || error();
A006046(n) == fromdigits([1,0,2, 0,0,4,8],3) || error();
}
\\ m = Mstartpos(j) = 4*A006046(j) is the m which is the first term of
\\ Mblock(j),
\\ high low
\\ for j = 2^k[1] + 2^k[2] + ... + 2^k[p]
\\ descending k[1] > k[2] > ... > k[p]
\\
\\ Mstartpos(j) = 4*3^k[1] + 8*3^k[2] + ... + 2^(p+1)*3^k[p]
\\
Mstartpos(j) = 4*A006046(j);
{
print1("Mlen(p) = ");
for(p=0,5, print1(Mlen(p),", "));
print("... (4*2^p)");
print1("Gould's = ");
for(j=0,7, print1(2^hammingweight(j),", "));
print("... (2^hammingweight(j), A001316)");
print1("A006046 = ");
for(j=0,8, print1(A006046(j),", "));
print("... (cumulative Gould's)");
print1("Mstartpos = ");
for(j=0,8, print1(Mstartpos(j),", "));
print("... (4*A006046)");
}
\\-----------------------------------------------------------------------------
\\ Powers and Offset
\\
\\ An efficient version of m_to_Mpr_by_slow_search() can be made by finding
\\ the unique representation of m as Mstartpos of a block and offset r into
\\ that block.
\\
\\ m = 4*3^k[1] + 8*3^k[2] + ... + 2^(p+1)*3^k[p] + r
\\
\\ where descending exponents k[1] > k[2] > ... > k[d] >= 0
\\ and "remainder" 0 <= r < 2^(p+2)
m_to_Mpr(m) =
{
my(p=0);
while(m >= 1<<(p+2),
my(k=logint(m>>(p+2),3));
m -= 2^(p+2) * 3^k; \\ successively reduced, ends as remainder r
p++);
[p,m];
}
b_by_Mpr(m) =
{
\\ same as b_by_Mpr_slow_search() but the fast m_to_Mpr()
my(p,r); [p,r] = m_to_Mpr(m);
Mblock_term(p,r);
}
{
\\ check b_by_Mpr(m) against a_vector()
my(v=a_vector(6));
vector(#v,m,m--; b_by_Mpr(m)) == v || error();
}
\\ m_to_Mpr() discards the k exponents it finds. If desired they could be
\\ made into a list (of p many) and so as to give the powers representation
\\ of m. But to calculate b(m), just p and r are enough.
\\-----------------------------------------------------------------------------
\\ State Machine High to Low
\\
\\ Mblock_term() treats the most significant 2 bits of r separately to get
\\ initial s for the morphism. An alternative is to have them as additional
\\ "states" (5 and 6) in a state machine which implements the morphism.
\\
\\ s=6 is the start state, for no bits of r yet seen. Its 0-bit case goes
\\ to state 1 which suits for next bit to go to s = 1 or 2. Or s=6 1-bit
\\ case goes to state s=5 which goes to s = 3 or 4 to complete the 1,2,3,4
\\ initial values.
HtoL_state_table = {[1, 2; \\ \
1, 3; \\ | same as morphism_table[]
4, 2; \\ |
4, 3; \\ /
3, 4; \\ state 5
1, 5]}; \\ state 6 start
{
\\ check first 4 rows of matrix same as morphism_table[]
for(s=1,4,
HtoL_state_table[s,] == morphism_table[s] || error());
}
\\ state = 1 to 6, bit = 0 or 1
HtoL_transition(state,bit) = HtoL_state_table[state,bit+1];
\\ same as Mblock_term() but using HtoL state machine
Mblock_term_by_HtoL(p,r) =
{
my(s=6);
forstep(i=p+1,0,-1, \\ from high to low
s = HtoL_transition(s, bittest(r,i)));
s;
}
\\ Bits of r put through the state machine starting at s=6.
b_by_HtoL(m) = \
my(p,r); [p,r] = m_to_Mpr(m); Mblock_term_by_HtoL(p,r);
{
\\ check b_by_HtoL(m) against a_vector()
my(v=a_vector(7));
vector(#v,m,m--; b_by_HtoL(m)) == v || error();
for(p=0,8,
for(r=0,Mlen(p)-1,
Mblock_term(p,r) == Mblock_term_by_HtoL(p,r) || error()));
}
\\-----------------------------------------------------------------------------
\\ State Machine Low to High
\\
\\ Within m_to_Mpr(), bits of r are determined from least to most significant.
\\ Some state machine manipulations can reverse HtoL to an equivalent machine
\\ taking bits in low to high order.
\\
\\ LtoH is bigger, but allows tighter code (and code where the reductions
\\ part is more like successive m_came_from_and_bit() actually).
\\
\\ See a073057-diagrams.pdf for a picture of LtoH, which ASCII art would not
\\ do justice to. The picture there or the numbers here are both frankly
\\ obscure without knowing their purpose. But it's the sort of thing where
\\ pre-calculating some table data can reduce code.
LtoH_state_table = {[
\\ generated by tools/gen-states-gp.pl short
1,2; \\ 1 2a -> 2a,3a
5,14; \\ 2 3a -> 2b,3c
8,11; \\ 3 4a -> 1b,4c
4,3; \\ 4 1a -> 1a,4a
1,9; \\ 5 2b -> 2a,2c
8,11; \\ 6 3g -> 1b,4c
18,10; \\ 7 4g -> 3d,3b
4,12; \\ 8 1b -> 1a,1c
5,14; \\ 9 2c -> 2b,3c
5,13; \\ 10 3b -> 2b,2d
19,15; \\ 11 4c -> 4d,4b
8,11; \\ 12 1c -> 1b,4c
18,10; \\ 13 2d -> 3d,3b
18,10; \\ 14 3c -> 3d,3b
8,16; \\ 15 4b -> 1b,1d
19,15; \\ 16 1d -> 4d,4b
5,13; \\ 17 2e -> 2b,2d
26,22; \\ 18 3d -> 3f,3e
23,27; \\ 19 4d -> 4f,4e
8,16; \\ 20 1e -> 1b,1d
18,17; \\ 21 2f -> 3d,2e
18,17; \\ 22 3e -> 3d,2e
23,24; \\ 23 4f -> 4f,1f
19,20; \\ 24 1f -> 4d,1e
0,0; \\ 25 unused
26,21; \\ 26 3f -> 3f,2f
19,20 \\ 27 4e -> 4d,1e
\\ end generated
]};
LtoH_transition(state,bit) = LtoH_state_table[state,bit+1];
\\ On ending in a state like "3f", the final result is b(m) = 3.
\\ Letters like "f" distinguish among various states of same result.
\\ States are numbered so that
\\
\\ b(m) = state mod 4 = result which is why unused state 25
\\ start state = 4 + (m mod 4)
\\ and otherwise arbitrary
\\
\\ The starting state is by explicitly taking the two lowest bits of r,
\\ which are also the two lowest bits of m. Further bits of r (bit
\\ positions 2 and above) are then transitions.
b_by_LtoH(m) =
{
my(state=(m%4)+4); \\ start state 4,5,6,7 by two low bits
m>>=2; \\ rest of m will be successively reduced
while(m,
m-=3^logint(m,3);
state=LtoH_transition(state,m%2);
m>>=1);
state%4 + 1;
}
{
my(v=a_vector(6));
vector(#v,m,m--; b_by_LtoH(m)) == v || error();
}
\\ In all these various reductions, it might we wondered how efficient
\\ k=logint(x,3) followed by 3^k is for finding powers of 3.
\\ The alternative would be some comparisons: Run down through powers
\\ T = 3^k and check m>=T to see if m is big enough to use this power.
\\
\\ Logs would be expected to be good if there's big gaps between T powers
\\ needed. Is that likely? Might imagine small gaps most of the time,
\\ depending on expected inputs.
\\
\\ An initial T can be any convenient over-estimate 3^k, as long as sure
\\ m < 3^(k+1). k_over_estimate() uses logint(m,2), which is just the bit
\\ length, and converts to power of 3 by factor
\\
\\ 12/19 = 0.6315... > log(2)/log(3) = 0.6309...
\\
\\ Any convenient similar rational would suit (a power-of-2 denominator
\\ would be multiply and shift). The next closer rational is 53/84.
k_over_estimate(m) = if(m==0,0, logint(m,2)*12\19+1);
b_by_LtoH_pow_steps(m) =
{
my(state=(m%4)+4);
m>>=2;
my(T=3^k_over_estimate(m));
m < 3*T || error("oops, over-estimate T too small, at m>>2 = ",m);
while(T>=1,
if(m>=T, \\ if big enough to use this T
m-=T;
state=LtoH_transition(state,m%2);
m>>=1);
T \= 3); \\ an exact division actually
state%4 + 1;
}
{
my(v=a_vector(6));
vector(#v,m,m--; b_by_LtoH_pow_steps(m)) == v || error();
}
\\ Another alternative would be to keep m in ternary digits all the time, so
\\ that m>=T just means the next highest non-0 ternary digit. But a
\\ digit-by-digit division by 2 in ternary is likely to be both slower an
\\ inconvenient (unless special tricks!).
\\-----------------------------------------------------------------------------
\\ Low to High State Machine Notes
\\
\\ In the LtoH diagram in the PDF, it can be noted the states are in two
\\ halves. The upper (as drawn) is values b(m) = 1,4 and lower b(m) = 2,3.
\\ This separation is since the morphism always expands
\\
\\ x -> [1 or 4], [2 or 3]
\\
\\ So that in Mblock(1) and above, the even r terms are 1 or 4 and the odd r
\\ terms are 2 or 3. Three bits of r reach Mblock(1) and above, so that in
\\ LtoH three bits are enough to decide between 1,4 and 2,3 (two bits for
\\ the starting state then one more).
\\
\\ States 3g and 4g exist to allow for Mblock(0) not following the even,odd
\\ pattern,
\\
\\ v-----v----- even r, would be each 1 or 4 in general
\\ Mblock(0) = 1, 2, 3, 4
\\ ^-----^-- odd r, would be each 2 or 3 in general
\\
\\ States 3g and 4g could be eliminated by noticing in the code when m<=3.
\\ But the reason for having a state machine is the converse: let the states
\\ handle cases so uniformity in the code.
\\ If the morphism is treated as a 4-state machine, starting at 1, then its
\\ reversal is two copies of the same 4 states and same transitions, but
\\ different accepting states, and which is simpler than LtoH.
\\
\\ Roughly speaking, LtoH is complicated since the HtoL start is not
\\ determined only at the highest 2 bits of r. LtoH remembers enough about
\\ the lower bits to allow for any stopping point to be the highest 2 bits.
\\-----------------------------------------------------------------------------
\\ Equal Proportions of 1,2,3,4
\\
\\ In the morphism expansion, four terms 1,2,3,4 expand to 1,1,2,2,3,3,4,4,
\\ though not in that order. So if started from equal numbers of 1s=2s=3s=4s
\\ then each expansion maintains that equality.
\\
\\ Mblock(0) = 1,2,3,4 has equal proportions, and hence so does every
\\ subsequent Mblock(p), and so the whole sequence has:
\\
\\ among the first Mstartpos(j) many terms of the sequence, for any j:
\\ number of 1s = number of 2s = number of 3s = number of 4s
\\
\\ These are also the only places with equality, since nowhere within an
\\ Mblock has equality. One way to see this is to interpret terms as
\\ segment directions in the alternate paperfolding curve,
\\
\\ 2
\\ ^ terms as alternate paperfolding
\\ | segment directions E,N,S,W
\\ 4 <---*---> 1
\\ |
\\ v
\\ 3
\\
\\ The morphism expansion is the alternate paperfolding curve directions
\\ expansion. Starting at the origin "O", the initial terms Mblock(0) =
\\ 1,2,3,4 are the following L shape return to the origin,
\\
\\ *
\\ ^ | initial
\\ 2| |3 Mblock(0) = 1,2,3,4
\\ 1 | v (which occurs just once)
\\ O -----> *
\\ <-----
\\ 4
\\
\\ One morphism expansion for Mblock(1) is a figure-of-8,
\\
\\ 4 1
\\ * <----- * -----> *
\\ | ^ ^ | Mblock(1) = 1,2, 1,3, 4,2, 4,3
\\ |3 2| |2 |3 figure of 8
\\ v 1 | | 4 v
\\ O -----> * <----- *
\\
\\ A square of sub-curves is what I call a "twin alternate" (my alternate
\\ paperfolding curve write-up). "Twin" here refers to the way two curves
\\ 1 and 4 back-to-back mesh perfectly, and on one expansion become a
\\ square. Similarly 2 and 3. This is more familiar in the dragon curve
\\ where two curves back-to-back are the twindragon.
\\
\\ Curve extents in the alternate paperfolding curve show that the
\\ figure-of-8 does not revisit the origin O again until the end. Nowhere
\\ within an Mblock has 1s=2s=3s=4s, because such equality would mean
\\
\\ 2s = 3s so point r on the "X" axis of the figure-of-8
\\ 1s = 4s so point r on the "Y" axis of the figure-of-8
\\ and both occurring at some r means the origin O
\\
\\ See the accompanying diagrams PDF for some bigger sample pictures.
vector_is_all_equal(v) = #Set(v)==1;
{
\\ check counts 1s=2s=3s=4s are if and only if m = Mstartpos(j)
my(k_limit=7,
v=a_vector(k_limit),
counts=[0,0,0,0],
j=0,
expect=Mstartpos(j), \\ next num terms to expect counts equal
n=0); \\ number of terms now in counts[]
while(n<#v,
vector_is_all_equal(counts) == (n==expect) || error();
if(n==expect,
expect=Mstartpos(j++)); \\ new expected num terms
n++;
counts[v[n]]++);
\\ the end of a_vector() is the end of its Mblock too
n == Mstartpos(j) || error();
vector_is_all_equal(counts) == (n==expect) || error();
}
\\-----------------------------------------------------------------------------
print("end");
\\ LocalWords: OEIS paperfolding halvings GRS fxtbook www org pdf
\\ LocalWords: Heighway Harter Joerg Arndt's Hasler's Ryde
\\ LocalWords: bittest logint hammingweight fromdigits
\\ LocalWords: len Mblock Mblocks Mlen Mpr Mstartpos HtoL LtoH
\\ LocalWords: dir twindragon PARI gp