/* Print out the coefficients for the geodesic series in a format that
allows Math::polyval to be used. */
taylordepth:5$
simplenum:true$
count:0$
jtaylor(expr,var1,var2,ord):=expand(subst([zz=1],
    ratdisrep(taylor(subst([var1=zz*var1,var2=zz*var2],expr),zz,0,ord))))$
h(x):=if x=0 then 0 else block([n:0],while not(oddp(x)) do (x:x/2,n:n+1),n);
decorhex(x):=if x=0 then "0" else if abs(x)<10^12 then string(x)
else concat(if x<0 then "-" else "",
  "0x",block([obase:16,s],
    s:sdowncase(string(abs(x))),
    if substring(s,1,2) = "0" then s:substring(s,2),
    s))$
formatnum(x):=block([n,neg:is(x<0),s1,s2,ll],x:abs(x),n:h(x),
  ll:if x<2^31 then "" else "LL",
  s1:concat(if neg then "-" else "", decorhex(x), ll),
  s2:concat(if neg then "-(" else "", decorhex(x/2^n), ll, "<<", n,
    if neg then ")" else ""),
  if slength(s2) < slength(s1) then s2 else s1)$
formatnumx(x):=if simplenum then
concat(string(x),if abs(x) < 2^31 then "" else "LL") else
if abs(x)<2^63 then (if abs(x) < 2^24
  /* test used to be: abs(x)/2^h(abs(x)) < 2^24; but Visual Studio
  complains about truncation of __int64 to real even if trailing bits
  are zero */
  then formatnum(x)
  else concat(s:if x<0 then "-" else "","real(",formatnum(abs(x)),")"))
else
block([sb:if x<0 then "-reale(" else "reale(",se:")"], x:abs(x),
  se:concat(formatnum(x-floor(x/2^52)*2^52),se),x:floor(x/2^52),
  while x > 0 do (
  se:concat(formatnum(x-floor(x/2^52)*2^52),",",se),x:floor(x/2^52)),
  concat(sb,se))$
printterm(x,line):=block([lx:slength(x)+1,lline:slength(line)],
  count:count+1,
  x:concat(x,if simplenum then ", " else ","),
  if lline=0 then line:concat("      ",x)
  else (if lx+lline<80 then line:concat(line,x)
    else (print(line),line:concat("      ",x))),
  line)$
flushline(line):=(if slength(line)>0 then (print(line),line:""),line)$
polyprint(p, var, title):=block([linel:90,d,line:"",h],
  p:ratsimp(p),
  d:abs(denom(p)),
  p:expand(d*p),
  h:hipow(p,var),
  print(concat("      // ", title, ", polynomial in ", var, " of order ",h)),
  for k:h step -1 thru 0 do
  line:printterm(formatnumx(coeff(p,var,k)),line),
  line:printterm(formatnumx(d),line),
  flushline(line),
  done)$

value1(val,var,pow,dord):=block([tab2:"    ",linel:90,div],
  print(concat(tab2,"// Generated by Maxima on ",timedate())),
  div:if pow = 1 then "" else concat("/",pow),
  for nn:0 step pow thru maxpow do block([c],
    if nn = 0 then
    print(concat("#if ", macro, div, " == ",nn/pow))
    else
    print(concat("#elif ", macro, div, " == ",nn/pow)),
    count:0,
    print(concat(tab2,"static const real coeff[] = {")),
    block(
      [q:ratdisrep(taylor(ev(val),var,0,nn-dord)),line:""],
      if pow = 2 then (
        q:subst(var=sqrt(concat(var,2)),expand(q)),
        polyprint(q,concat(var,2),string(val)))
      else (polyprint(q,var,string(val))),
      line:flushline(line)),
    print(concat(tab2,"};  // count = ",count))),
  print("#else
#error", concat("\"Bad value for ", macro, "\""), "
#endif
"),
'done)$

array1(array,var,pow,dord):=block([tab2:"    ",linel:90],
  print(concat(tab2,"// Generated by Maxima on ",timedate())),
  for nn:0 thru maxpow do block([c],
    if nn = 0 then
    print(concat("#if ", macro, " == ",nn))
    else
    print(concat("#elif ", macro, " == ",nn)),
    count:0,
    print(concat(tab2,"static const real coeff[] = {")),
    for m:0 thru nn do if part(arrayapply(array,[m]),0) # array then block(
      [q:ratdisrep(taylor(arrayapply(array,[m]),var,0,nn-dord)),line:""],
      if pow = 2 then (
        q:subst(var=sqrt(concat(var,2)),expand(q/var^(m-dord))),
        polyprint(q,concat(var,2),concat(array, "[", m, "]/",var,"^",m-dord)))
      else (
        q:expand(q/var^(m-dord)),
        polyprint(q,var,concat(array, "[", m, "]/",var,"^",m-dord))),
      line:flushline(line)),
    print(concat(tab2,"};  // count = ",count))),
  print("#else
#error", concat("\"Bad value for ", macro, "\""), "
#endif
"),
'done)$

value2(val,var1,var2,dord):=block([tab2:"    ",linel:90],
  print(concat(tab2,"// Generated by Maxima on ",timedate())),
  for nn:0 thru maxpow do block([c],
    if nn = 0 then
    print(concat("#if ", macro, " == ",nn))
    else
    print(concat("#elif ", macro, " == ",nn)),
    count:0,
    print(concat(tab2,"static const real coeff[] = {")),
    block(
      [q:jtaylor(ev(val),n,eps,nn-dord),line:""],
      for j:nn-dord step -1 thru 0 do block(
        [p:coeff(q,var2,j)],
        polyprint(p,var1,concat(val, ", coeff of eps^",j))),
      line:flushline(line)),
    print(concat(tab2,"};  // count = ",count))),
  print("#else
#error", concat("\"Bad value for ", macro, "\""), "
#endif
"),
'done)$

array2(array,var1,var2,dord):=block([tab2:"    ",linel:90],
  print(concat(tab2,"// Generated by Maxima on ",timedate())),
  for nn:0 thru maxpow do block([c],
    if nn = 0 then
    print(concat("#if ", macro, " == ",nn))
    else
    print(concat("#elif ", macro, " == ",nn)),
    count:0,
    print(concat(tab2,"static const real coeff[] = {")),
    for m:0 thru nn-1 do if part(arrayapply(array,[m]),0) # array then block(
      [q:jtaylor(arrayapply(array,[m]),n,eps,nn-dord),line:""],
      for j:nn-dord step -1 thru m do block(
        [p:coeff(q,var2,j)],
        polyprint(p,var1,concat(array, "[", m ,"], coeff of eps^",j))),
      line:flushline(line)),
    print(concat(tab2,"};  // count = ",count))),
  print("#else
#error", concat("\"Bad value for ", macro, "\""), "
#endif
"),
'done)$