147 lines
5.3 KiB
Plaintext
147 lines
5.3 KiB
Plaintext
/* 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)$
|