\\ Kevin Ryde, November 2021
\\
\\ Usage: gp <a349008.gp
\\
\\ Calculate terms of A349008, being area "HA" of the convex hull around
\\ R5 dragon curve expansion level n.
\\
\\ The implementation follows the formulas in the sequence except powers of
\\ 5 and powers of b are maintained instead of say 5^(n-2-j) each time.
\\ (This is a speedup, but by less than it might be thought actually.)

default(strictargs,1);
default(recover,0);

b = 1+2*I;

\\ return z rotated if necessary so it's in quadrant 1 or 3
RotQ(z) = if(sign(real(z))==sign(imag(z)), z, z*I);

\\ return complex z with the real part sheared by doubling
ShearRe(z) = z + real(z);

MinReIm(z) = min(abs(real(z)),abs(imag(z)));

HAgrow(z) = MinReIm(ShearRe(RotQ(z)));

a(n) =
{
  if(n==0, 0,     \\ initial exceptions, then general case formula
     n==1, 2,
     my(B=1,        \\ b^j
        F=5^(n-2)); \\ 5^(n-2-j)
     17*F - 1
     + sum(j=1,n-2,
           F/=5; B*=b;
           2*((3*F-1)*HAgrow(      B)
              +  F   *HAgrow((4-I)*B))));
}

{
  print("A349008 hull areas, for n>=0");
  print1("  ");
  for(n=0,9, print1(a(n),", "));
  print("...\n");
}

\\-----------------------------------------------------------------------------
\\ Fractal Limit
\\
\\ As noted in the sequence, a(n) grows as
\\
\\         a(n)
\\     lim ----  ->  A349009 = "HAf" = 0.9761640029...
\\         5^n
\\
\\ Limit HAf is the area of the hull around the R5 dragon fractal.
\\
\\ The formulas in the limit A349009 differ from here in that they can omit
\\ the "-1" from "17*F-1" and "3*F-1" since /5^n means they -> 0.  (Other
\\ apparent differences in A349009 are just some conjugating so that b^j/5^j
\\ = conj(1/b^j) can be written without the conj.)

{
  print("Print some a(n)/5^n, but with denominators 10^n to show");
  print("how decimals -> A349009 = 0.9761640029...");

  my(n_list=[1,2,3,4,5,6, 10,15,20,30,40]);
  for(i=1,#n_list,
     my(n = n_list[i],
        rat = a(n) / 5^n,
        den = 10^n,
        num = rat*den);
     num/den == a(n)/5^n || error();
     printf("  n=%2d:  %s / 10^%d\n", n, num, n));
}

\\-----------------------------------------------------------------------------
print("end");