\\ Kevin Ryde, November 2021
\\
\\ Usage: gp <a349009.gp
\\
\\ Calculate terms of A349009, being decimal digits of the area "HAf" of the
\\ convex hull around the R5 dragon fractal.
\\
\\ The method is the formulas presented in the sequence, with a scale factor
\\ 10^len so that the desired decimals are the integer digits.  The formulas
\\ are lightly adapted by moving "scale * 2/25" inside HAgrowf, since it is
\\ linear HAgrowf(c*z) = c*HAgrowf(z) for real c>=0.
\\
\\ The loop end condition uses the sum upper bound given in the sequence.
\\ When "omitted" terms are too small to affect the integer part of the
\\ scaled HAf, then the integer part and hence desired digits are certain.

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

b = 1+2*I;

\\ return z rotated if necessary to put it in the 1st or 3rd quadrant
RotQ(z) = if(sign(real(z))==sign(imag(z)), z, z*I);

\\ return complex z with the imaginary part sheared by doubling
ShearIm(z) = z + I*imag(z);

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

HAgrowf(z) = MinReIm(ShearIm(RotQ(z)));

\\ Return a vector of the first len many decimal digits of A349009 = HAf
a_vector(len) =
{
  my(scale = 10^len,       \\ to get len many decimal digits
     j = 0,
     HAf = scale * 17/25,  \\ start with constant part
     B = scale * 2/25,     \\ scale * 2/25 / b^j
     omitted = scale);     \\ scale * 1/sqrt(5)^j
  while(1,
        j++;
        if(j%2==0,
           omitted /= 5;
           \\ omitted == scale/5^(j/2) || error();
           if(frac(HAf) + omitted < 1,
              return(digits(floor(HAf)))));
        B /= b;
        HAf += 3*HAgrowf(      B)
             +   HAgrowf((4+I)*B)) ;
}

{
  my(v=a_vector(24));
  print("A349009 digits");
  print1("  ");
  for(i=1,#v, print1(v[i],", "));
  print("...");
  print("  = 0.",fromdigits(v),"...");
}

print("end");