1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
|
/* 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)$
|
