fcode(contract=True) fails interpreting indexed with arithmetic for indexes

[ThePauliPrinciple]

Taken from http://docs.sympy.org/latest/_modules/sympy/printing/fcode.html

from sympy import Eq, IndexedBase, Idx, fcode
len_y = 5
y = IndexedBase('y', shape=(len_y,))
t = IndexedBase('t', shape=(len_y,))
Dy = IndexedBase('Dy', shape=(len_y-1,))
i = Idx('i', len_y-1)
e=Eq(Dy[i], (y[i+1]-y[i])/(t[i+1]-t[i]))
print fcode(e.rhs, assign_to=e.lhs, contract=False)
print fcode(e.rhs, assign_to=e.lhs, contract=True)

I added the last line to explicitely write the loop structure. Unfortunately this does not work.

Here a minimal example:

e=Eq(Dy[i], y[i])
print(fcode(e.rhs, assign_to=e.lhs, contract=True))

works as expected

e=Eq(Dy[i], y[i+1])
print(fcode(e.rhs, assign_to=e.lhs, contract=True))

However does not. The error message seems to suggest determining the loop bounds is unsuccesfull, supposedly because of the arithmetic performed on the index.

See also
https://groups.google.com/forum/#!topic/sympy/sP6vhwtCELw

[manauref]

Has this been resolved? It's been many months since this was brought up.

This is by the way also broken with other code printers. I've tried with ccode and cxxcode and it is also broken.

A fix to this feature of Sympy, or alternative solutions would be greatly appreciated.

[oscarbenjamin]

@manauref maybe you could open a PR?

I have no experience of using any of the code printers and don't know the printing code at all but I am happy to give general advice.