Help with sympy, please

[Dick Moores]

from __future__ import division
Here's what I'm trying to do, but using sympy:
=============================
from math import e
n = 1
prod = 1
k = 0
while k < 1000:
k += 1
term = (e**(1.0/n))/(e**(1.0/(n+1)))
prod *= term
n += 2
print prod, term
================================
Output:
1.99950018746 1.00000025013 (the limit is 2)


============================
from sympy import *
from sympy import Rational as R
from sympy.numerics import *

prec = 50
Float.setdps(prec)
e = evalf(E)
n = 1
k = 0
prod = evalf(R(1,1))
while k < 1000:
k += 1
n = evalf(R(n,1))
term = (e**(1/n))/(e**(1/(n+1)))
prod *= term
n += 2
print prod, term
===========================

This gets:
Traceback (most recent call last):
File "E:\PythonWork\Untitled 5.py", line 20, in <module>
term = (e**(1/n))/(e**(1/(n+1)))
TypeError: unsupported operand type(s) for /: 'int' and 'Float'

Thanks,

Dick Moores

 

[Dennis Lee Bieber]

This gets:
Traceback (most recent call last):
File "E:\PythonWork\Untitled 5.py", line 20, in <module>
term = (e**(1/n))/(e**(1/(n+1)))
TypeError: unsupported operand type(s) for /: 'int' and 'Float'

Seems self-explanatory... try using 1.0 rather than 1
--
Wulfraed Dennis Lee Bieber KD6MOG
(e-mail address removed) (e-mail address removed)
HTTP://wlfraed.home.netcom.com/
(Bestiaria Support Staff: (e-mail address removed))
HTTP://www.bestiaria.com/

 

[Dick Moores]

Seems self-explanatory... try using 1.0 rather than 1

term = (e**(1.0/n))/(e**(1.0/(n+1)))
TypeError: unsupported operand type(s) for /: 'float' and 'Float'

Dick

 

[Fredrik Johansson]

term = (e**(1.0/n))/(e**(1.0/(n+1)))
TypeError: unsupported operand type(s) for /: 'float' and 'Float'

Hi Dick, I recognize you from python-list, where you had a question
about mpmath.

Your code still won't work if you convert the numbers to Floats
because the Float type in sympy.numerics does not implement ** for
fractional numbers. You could use the exp function in
sympy.numerics.functions instead to compute e**x.

Basically, sympy.numerics is an old version of mpmath. The
sympy.numerics module is not very well integrated in SymPy, slower
than mpmath, and has a couple bugs that have subsequently been fixed
in mpmath. In sympycore (http://code.google.com/p/sympycore/), we're
using the latest version of mpmath and integrating it directly into
the symbolic engine; it will be much more robust and user-friendly. It
will hopefully not be long until we merge the improvements we've done
in the sympycore project with the main SymPy branch.

Fredrik

 

[Dick Moores]

Hi Dick, I recognize you from python-list, where you had a question
about mpmath.

Your code still won't work if you convert the numbers to Floats
because the Float type in sympy.numerics does not implement ** for
fractional numbers. You could use the exp function in
sympy.numerics.functions instead to compute e**x.

Thanks, Fredrik, but I get the same error using either exp or power:
============================
from __future__ import division
from sympy import *
from sympy import Rational as R
from sympy.numerics import *
from sympy.numerics.functions import power, exp

prec = 50
Float.setdps(prec)
e = evalf(E)
n = 1
m = n + 1
k = 0
n = evalf(R(n,1))
m = evalf(R(m,1))
prod = evalf(R(1,1))
prec = 50
Float.setdps(prec)
e = evalf(E)
n = 1
k = 0
prod = evalf(R(1,1))
while k < 1000:
k += 1
n = evalf(R(n,1))
term = (exp(1/n))/(exp(1/(n+1)))
prod *= term
n += 2
print prod, term
=============================
term = (exp(1/n))/(exp(1/(n+1)))
TypeError: unsupported operand type(s) for /: 'int' and 'Float'

And also if I use sympy.numerics.functions.power in that line, as
term = power(e,(1/n))/power(e,(1/(n+1)))

I get:

term = power(e,(1/n))/power(e,(1/(n+1)))
TypeError: unsupported operand type(s) for /: 'int' and 'Float'

Dick

 

[Dennis Lee Bieber]

TypeError: unsupported operand type(s) for /: 'float' and 'Float'

<blink><blink>

What restrictive system can't even coerce a Python float to their
own internal data type?

Does that library have any functions for type conversion? Float(1)
(or Float(1.0) ) perhaps?
--
Wulfraed Dennis Lee Bieber KD6MOG
(e-mail address removed) (e-mail address removed)
HTTP://wlfraed.home.netcom.com/
(Bestiaria Support Staff: (e-mail address removed))
HTTP://www.bestiaria.com/

 

[Dick Moores]

Basically, sympy.numerics is an old version of mpmath. The
sympy.numerics module is not very well integrated in SymPy, slower
than mpmath, and has a couple bugs that have subsequently been fixed
in mpmath. In sympycore (http://code.google.com/p/sympycore/), we're
using the latest version of mpmath and integrating it directly into
the symbolic engine; it will be much more robust and user-friendly. It
will hopefully not be long until we merge the improvements we've done
in the sympycore project with the main SymPy branch.

Fredrik

OK, I tried mpmath again, and to my surprise, it went well!

===================================
#!/usr/bin/env python
#coding=utf-8
from mpmath import *
mpf.dps = 50
n = 1
k = 0
prod = mpf(1)
while k < 100000:
k += 1
term = exp(1.0/n)/exp(1.0/(n+1))
prod *= term
n += 2
print prod, term
======================================
Output:
1.9999950000187499635016028080844735182389158683797
1.0000000000250001250004074790133889386806610626172

Thanks,

Dick

 

[Fredrik Johansson]

Thanks, Fredrik, but I get the same error using either exp or power:
============================
from __future__ import division
from sympy import *
from sympy import Rational as R
from sympy.numerics import *
from sympy.numerics.functions import power, exp

prec = 50
Float.setdps(prec)
e = evalf(E)
n = 1
m = n + 1
k = 0
n = evalf(R(n,1))
m = evalf(R(m,1))
prod = evalf(R(1,1))
prec = 50
Float.setdps(prec)
e = evalf(E)
n = 1
k = 0
prod = evalf(R(1,1))
while k < 1000:
k += 1
n = evalf(R(n,1))
term = (exp(1/n))/(exp(1/(n+1)))
prod *= term
n += 2
print prod, term
=============================
term = (exp(1/n))/(exp(1/(n+1)))
TypeError: unsupported operand type(s) for /: 'int' and 'Float'

And also if I use sympy.numerics.functions.power in that line, as
term = power(e,(1/n))/power(e,(1/(n+1)))

I get:

term = power(e,(1/n))/power(e,(1/(n+1)))
TypeError: unsupported operand type(s) for /: 'int' and 'Float'

Try not using "from __future__ import division".

Also,

n = evalf(R(n,1))

won't work, because the arguments to Rational must be ints.

Fredrik

 

[Fredrik Johansson]

<blink><blink>

What restrictive system can't even coerce a Python float to their
own internal data type?

Does that library have any functions for type conversion? Float(1)
(or Float(1.0) ) perhaps?

Coercion works fine; the problem is that Float implements __div__ but
not __truediv__. It works if you don't import division from
__future__. If you do that, even Float / Float division ceases to
work. That is definitely a bug, but not as bad a bug as not being able
to coerce floats at all ;-)

Fredrik

 

[Fredrik Johansson]

OK, I tried mpmath again, and to my surprise, it went well!

===================================
#!/usr/bin/env python
#coding=utf-8
from mpmath import *
mpf.dps = 50
n = 1
k = 0
prod = mpf(1)
while k < 100000:
k += 1
term = exp(1.0/n)/exp(1.0/(n+1))
prod *= term
n += 2
print prod, term
======================================
Output:
1.9999950000187499635016028080844735182389158683797
1.0000000000250001250004074790133889386806610626172

You're getting slightly wrong results, though, because 1.0/n and
1.0/(n+1) just performs regular float division with ~16-digit
precision when n is a Python int. You should convert either 1.0 or n
to an mpf before dividing.

Fredrik

 

[Dick Moores]

Try not using "from __future__ import division".

Didn't help.

Also,

n = evalf(R(n,1))

won't work, because the arguments to Rational must be ints.

But both n and 1 are ints, aren't they?

Anyway, you've by now seen my new success with mpmath. I'll abandon
syspy for the time being.

Dick

 

[Dick Moores]

You're getting slightly wrong results, though, because 1.0/n and
1.0/(n+1) just performs regular float division with ~16-digit
precision when n is a Python int. You should convert either 1.0 or n
to an mpf before dividing.

Ah, yes. That was my original mistake with using mpmath. How's this?:

=============================
#!/usr/bin/env python
#coding=utf-8
import psyco
psyco.full()
from mpmath import *
mpf.dps = 50
n = mpf(1)
k = 0
prod = mpf(1)
while k < 100000:
k += 1
term = exp(1.0/n)/exp(1.0/(n+1))
prod *= term
n += 2
print prod, term
===============================
Output:
1.9999950000187499635415917971337956346129920295869
1.0000000000250001250009375062500416669401059407666

Dick