I have some Maple code below that I'm trying to convert to Sage (which is some kind of Python) code.
This is the Maple code
restart:
m:=15:
a:=array(1..m):
eqn:=array(1..m):
expr:=sum(a['i']*u^('i'-1),'i'=1..m)-
product((1-u^'i')^(-a['i']),'i'=1..m):
for i from 1 to m do eqn[i]:=coeff(series(expr,u=0,m+2),u,i-1);
od:
sols:=solve({seq(eqn[i],i=1..m)},{seq(a[i],i=1..m)}):
assign(sols):
print(a);
This is the output for this code:
[1, 1, 2, 4, 9, 20, 48, 115, 286, 719, 1842, 4766, 12486, 32973, 87811]
This is the Sage code I have already:
u = var('u')
m = 15
a = [var(f"a_{i}") for i in range(1,m+1)]
eqn = [var(f"e_{i}") for i in range(1,m+1)]
expr = sum(a[i]*u^(i-1) for i in range(1,m)) - product((1-u^i)^(-a[i]) for i in range(1,m))
print(expr)
for i in range(0,m):
# eqn[i] = taylor(expr,u,0,m+2).coefficients(sparse=False)[i]
When I uncomment the code in the for-loop, I get IndexError: list index out of range.
I tried the CodeGeneration tool:
CodeGeneration:-Python(i -> coeff(series(expr, u = 0, m + 2), u, i - 1), output = embed);
CodeGeneration:-Python(sols -> assign(sols), output = embed);
This however gives me
Warning, the function names {coeff, series} are not recognized in the target language
Warning, the function names {assign} are not recognized in the target language
And it gives as output thus no useful code, since coeff and series don't exist:
cg0 = lambda i: coeff(series(expr, u == 0, m + 2), u, i - 1)
cg2 = lambda sols: assign(sols)
The question is now: what are the equivalent expressions for coeff, series and assign?
To give an answer, one has to understand first what happens behind the maple code, the solve the problem more or less in an optimal manner in sage. Here is what can be done in sage, to look as close as possible to the maple way.
m = 15
r = [1..m] # range to be used in the sequel
u = var('u');
a = dict(zip(r, var('a', n=m+1, latex_name='a')[1:]))
expr = sum ( a[j]*u^(j-1) for j in r ) \
- prod( (1 - u^j)^(-a[j]) for j in r )
expr = taylor(expr, u, 0, m)
sols = solve( [ diff(expr, u, j-1).subs(u == 0) for j in r ]
, list(a.values())
, solution_dict=True)
for sol in sols:
print(sol)
And there is only one solution, and the one print is:
{a1: 1, a2: 1, a3: 2, a4: 4, a5: 9, a6: 20, a7: 48, a8: 115, a9: 286
, a10: 719, a11: 1842, a12: 4766, a13: 12486, a14: 32973, a15: 87811}
Here are some words about the decision to use the "range" list (which is a list object, not a range object)... In the maple code, a was something that accepted a[1] and a[2] and so on. These "symbols" are variables, if not specified in a different way. But there is no a[0]. My choice was to use
var('a', n=m+1, latex_name='a')
which is a list of "variables", well it is the list starting with a0, so i am taking it from a1 to the end via (slicing)
var('a', n=m+1, latex_name='a')[1:]
and then, in order to have something like a[1] pointing to a1 i am using a dictionary in sage. So i am zipping / pairing the list r with the above sliced part of the var's. Then is want the dictionary that solidifies this pairing.
Then in order to have the u-coefficients, i am using the poor man's coefficient extraction, obtained by differentiation w.r.t. u, as many times as needed, then setting the u variable to zero.
Why is the above not working "somehow" by directly "taking coefficients"? Because in sage, expr is an expression.
sage: type(expr)
<class 'sage.symbolic.expression.Expression'>
And there is no way in general to get the coefficients of such an "expression". Instead, if we arrange to work in polynomial rings, then there is still a chance to proceed using the wanted method. Please compare the above with:
m = 15
r = [1..m] # range to be used in the sequel
R.<a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15> = \
PolynomialRing(QQ)
a = dict(zip(r, R.gens()))
S.<u> = PowerSeriesRing(R, default_prec=m)
expr = sum ( a[j]*u^(j-1) for j in r ) \
- exp( sum( -a[j]*log(1 - u^j) for j in r ) )
J = R.ideal(expr.coefficients())
J.variety() # solution(s) to the R-coefficient equations from above over QQ, base ring of R
And that J.variety() delivers...
sage: J.variety()
[{a15: 87811, a14: 32973, a13: 12486, a12: 4766, a11: 1842, a10: 719
, a9: 286, a8: 115, a7: 48, a6: 20, a5: 9, a4: 4, a3: 2, a2: 1, a1: 1}]
(Results were manually adjusted.)