Re: [sympy] simplify bug with powers > 20
28 views
Skip to first unread message
Aaron Meurer
Sep 18, 2012, 2:09:35 PM9/18/12
to sy...@googlegroups.com
This is indeed a bug!
The problem appears to be in cancel(). It creates two polynomials,
Poly(1317378933230047068160*sqrt(5) + 2945748836994210856960, sqrt(5),
domain='ZZ') Poly(120352542776360960*sqrt(5) + 269116466014453760,
sqrt(5), domain='ZZ'), and applies Euclid's algorithm on them to find
the gcd, and then polynomial division. The problem is that the
generator sqrt(5) is an algebraic, so certain polynomials in it will
be identically zero. But Poly does not know this, so somewhere in gcd
or div it comes across such a polynomial and assumes that it is not
zero, leading to the incorrect result. If you want, you may try
applying the algorithm by hand as an exercise to see where it goes
awry. You might even be able to figure out why 21 is the magic number
where it starts to go bad.
If you use expand() and then radsimp(), you get the correct result.
The only way to really fix this would be to prevent cancel, and really
Poly in general, from creating generators from algebraic numbers.
This would unfortunately break a lot of things that otherwise work
just fine (because they never come across these identically zero
polynomials in the generators).
I should also point out that if you create the polynomials using an
algebraic extension, it works correctly:
In [20]: F, G = Poly(1317378933230047068160*sqrt(5) +
2945748836994210856960, sqrt(5), extension=sqrt(5)),
Poly(120352542776360960*sqrt(5) + 269116466014453760, sqrt(5),
extension=sqrt(5))
In [21]: F.cancel(G)Out[21]: (1, Poly(1317378933230047068160, sqrt(5),
domain='QQ<sqrt(5)>'), Poly(120352542776360960, sqrt(5),
domain='QQ<sqrt(5)>'))
In [22]: 1317378933230047068160/120352542776360960
Out[22]: 10946.0
But unfortunately, this can get quite slow if you have more than one extension.
Would you mind opening an issue in our issue tracker for this:
http://code.google.com/p/sympy/issues/list.
Aaron Meurer
On Tue, Sep 18, 2012 at 10:38 AM, Oscar <oscar.j....@gmail.com> wrote:
> Hi all,
>
> Having just discovered Binet's formula for the nth fibonacci number I
> thought
> I'd give it a try in sympy and stumbled across a possible bug.
>
> Binet's formula, although it has a number of irrational quantities, should
> always yield an integer. I found that sympy was able to compute this
> correctly
> up to the 20th Fibonacci number and that started returning irrational
> quantities when I asked it to simplify the expression.
>
> Demonstration using a fresh checkout of sympy master in a freshly created,
> otherwise empty Python 2.7.1 virtual env:
>
>>>> import sympy
>>>> n = sympy.Symbol('n')
>>>> Phi = (1 + sympy.sqrt(5)) / 2
>>>> fibn = (Phi**n - (-Phi)**(-n)) / sympy.sqrt(5)
>>>> for i in range(1, 25):
> ... print '%ith Fibonacci number:' % i, fibn.subs(n, i).simplify()
> ...
> 1th Fibonacci number: 1
> 2th Fibonacci number: 1
> 3th Fibonacci number: 2
> 4th Fibonacci number: 3
> 5th Fibonacci number: 5
> 6th Fibonacci number: 8
> 7th Fibonacci number: 13
> 8th Fibonacci number: 21
> 9th Fibonacci number: 34
> 10th Fibonacci number: 55
> 11th Fibonacci number: 89
> 12th Fibonacci number: 144
> 13th Fibonacci number: 233
> 14th Fibonacci number: 377
> 15th Fibonacci number: 610
> 16th Fibonacci number: 987
> 17th Fibonacci number: 1597
> 18th Fibonacci number: 2584
> 19th Fibonacci number: 4181
> 20th Fibonacci number: 6765
> 21th Fibonacci number: -4498 + sqrt(5) <--- should be 10946
> 22th Fibonacci number: -10108 + sqrt(5)
> 23th Fibonacci number: -6785 + sqrt(5)
> 24th Fibonacci number: sqrt(5) + 14490
>>>> fibn.subs(n, 21)
> sqrt(5)*(-1/(-sqrt(5)/2 - 1/2)**21 + (1/2 + sqrt(5)/2)**21)/5
>>>> fibn.subs(n, 21).simplify()
> -4498 + sqrt(5)
>>>> fibn.subs(n, 21).simplify().evalf()
> -4495.76393202250
>>>> fibn.subs(n, 21).evalf()
> 10946.0000000000
>
> That last line shows that the expression is correct before it is
> simplify()ed.
>
> Is this a known problem? I'd be happy to help fix this but I have little
> idea
> where to start. Is there any significance in the fact that it goes wrong at
> powers > 20?
>
> Oscar
>
> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To view this discussion on the web visit
> https://groups.google.com/d/msg/sympy/-/gH8z3L4TgbUJ.
> To post to this group, send email to sy...@googlegroups.com.
> To unsubscribe from this group, send email to
> sympy+un...@googlegroups.com.
> For more options, visit this group at
> http://groups.google.com/group/sympy?hl=en.
The problem appears to be in cancel(). It creates two polynomials,
Poly(1317378933230047068160*sqrt(5) + 2945748836994210856960, sqrt(5),
domain='ZZ') Poly(120352542776360960*sqrt(5) + 269116466014453760,
sqrt(5), domain='ZZ'), and applies Euclid's algorithm on them to find
the gcd, and then polynomial division. The problem is that the
generator sqrt(5) is an algebraic, so certain polynomials in it will
be identically zero. But Poly does not know this, so somewhere in gcd
or div it comes across such a polynomial and assumes that it is not
zero, leading to the incorrect result. If you want, you may try
applying the algorithm by hand as an exercise to see where it goes
awry. You might even be able to figure out why 21 is the magic number
where it starts to go bad.
If you use expand() and then radsimp(), you get the correct result.
The only way to really fix this would be to prevent cancel, and really
Poly in general, from creating generators from algebraic numbers.
This would unfortunately break a lot of things that otherwise work
just fine (because they never come across these identically zero
polynomials in the generators).
I should also point out that if you create the polynomials using an
algebraic extension, it works correctly:
In [20]: F, G = Poly(1317378933230047068160*sqrt(5) +
2945748836994210856960, sqrt(5), extension=sqrt(5)),
Poly(120352542776360960*sqrt(5) + 269116466014453760, sqrt(5),
extension=sqrt(5))
In [21]: F.cancel(G)Out[21]: (1, Poly(1317378933230047068160, sqrt(5),
domain='QQ<sqrt(5)>'), Poly(120352542776360960, sqrt(5),
domain='QQ<sqrt(5)>'))
In [22]: 1317378933230047068160/120352542776360960
Out[22]: 10946.0
But unfortunately, this can get quite slow if you have more than one extension.
Would you mind opening an issue in our issue tracker for this:
http://code.google.com/p/sympy/issues/list.
Aaron Meurer
On Tue, Sep 18, 2012 at 10:38 AM, Oscar <oscar.j....@gmail.com> wrote:
> Hi all,
>
> Having just discovered Binet's formula for the nth fibonacci number I
> thought
> I'd give it a try in sympy and stumbled across a possible bug.
>
> Binet's formula, although it has a number of irrational quantities, should
> always yield an integer. I found that sympy was able to compute this
> correctly
> up to the 20th Fibonacci number and that started returning irrational
> quantities when I asked it to simplify the expression.
>
> Demonstration using a fresh checkout of sympy master in a freshly created,
> otherwise empty Python 2.7.1 virtual env:
>
>>>> import sympy
>>>> n = sympy.Symbol('n')
>>>> Phi = (1 + sympy.sqrt(5)) / 2
>>>> fibn = (Phi**n - (-Phi)**(-n)) / sympy.sqrt(5)
>>>> for i in range(1, 25):
> ... print '%ith Fibonacci number:' % i, fibn.subs(n, i).simplify()
> ...
> 1th Fibonacci number: 1
> 2th Fibonacci number: 1
> 3th Fibonacci number: 2
> 4th Fibonacci number: 3
> 5th Fibonacci number: 5
> 6th Fibonacci number: 8
> 7th Fibonacci number: 13
> 8th Fibonacci number: 21
> 9th Fibonacci number: 34
> 10th Fibonacci number: 55
> 11th Fibonacci number: 89
> 12th Fibonacci number: 144
> 13th Fibonacci number: 233
> 14th Fibonacci number: 377
> 15th Fibonacci number: 610
> 16th Fibonacci number: 987
> 17th Fibonacci number: 1597
> 18th Fibonacci number: 2584
> 19th Fibonacci number: 4181
> 20th Fibonacci number: 6765
> 21th Fibonacci number: -4498 + sqrt(5) <--- should be 10946
> 22th Fibonacci number: -10108 + sqrt(5)
> 23th Fibonacci number: -6785 + sqrt(5)
> 24th Fibonacci number: sqrt(5) + 14490
>>>> fibn.subs(n, 21)
> sqrt(5)*(-1/(-sqrt(5)/2 - 1/2)**21 + (1/2 + sqrt(5)/2)**21)/5
>>>> fibn.subs(n, 21).simplify()
> -4498 + sqrt(5)
>>>> fibn.subs(n, 21).simplify().evalf()
> -4495.76393202250
>>>> fibn.subs(n, 21).evalf()
> 10946.0000000000
>
> That last line shows that the expression is correct before it is
> simplify()ed.
>
> Is this a known problem? I'd be happy to help fix this but I have little
> idea
> where to start. Is there any significance in the fact that it goes wrong at
> powers > 20?
>
> Oscar
>
> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To view this discussion on the web visit
> https://groups.google.com/d/msg/sympy/-/gH8z3L4TgbUJ.
> To post to this group, send email to sy...@googlegroups.com.
> To unsubscribe from this group, send email to
> sympy+un...@googlegroups.com.
> For more options, visit this group at
> http://groups.google.com/group/sympy?hl=en.
Oscar
Sep 18, 2012, 6:14:16 PM9/18/12
to sy...@googlegroups.com
On Tuesday, 18 September 2012 19:09:57 UTC+1, Aaron Meurer wrote:
This is indeed a bug!
The problem appears to be in cancel(). It creates two polynomials,
Poly(1317378933230047068160*sqrt(5) + 2945748836994210856960, sqrt(5),
domain='ZZ') Poly(120352542776360960*sqrt(5) + 269116466014453760,
sqrt(5), domain='ZZ'), and applies Euclid's algorithm on them to find
the gcd, and then polynomial division. The problem is that the
generator sqrt(5) is an algebraic, so certain polynomials in it will
be identically zero. But Poly does not know this, so somewhere in gcd
or div it comes across such a polynomial and assumes that it is not
zero, leading to the incorrect result. If you want, you may try
applying the algorithm by hand as an exercise to see where it goes
awry. You might even be able to figure out why 21 is the magic number
where it starts to go bad.
If you use expand() and then radsimp(), you get the correct result.
The only way to really fix this would be to prevent cancel, and really
Poly in general, from creating generators from algebraic numbers.
This would unfortunately break a lot of things that otherwise work
just fine (because they never come across these identically zero
polynomials in the generators).
It turns out to be a lot simpler than that. There was a typo in the return value from dup_zz_heu_gcd() in the euclidtools module. I've sent a pull request here:
I've no idea why it only occurs for all fibonacci numbers above the 20th as I don't fully understand the heuristic algorithm. I'm pretty confident about the typo though.
Also, it takes a really long time to run the test suite on my computer (I went to the shop and it was still going when I got back). Do people really run that every time before committing?
Oscar
Aaron Meurer
Sep 18, 2012, 8:50:08 PM9/18/12
to sy...@googlegroups.com
because Poly doesn't ever check for algebraic relations on the
generators, even trivial ones (like sqrt(5)**2 == 5).
>
> I've no idea why it only occurs for all fibonacci numbers above the 20th as
> I don't fully understand the heuristic algorithm. I'm pretty confident about
> the typo though.
>
> Also, it takes a really long time to run the test suite on my computer (I
> went to the shop and it was still going when I got back). Do people really
> run that every time before committing?
that is hanging?
Most people honestly run only those parts of the test suite that they
know might be affect (so e.g., if you change something in the stats
module, then you only need to run the stats tests). But that's only if
you have a good idea of what might or might not be affected by your
change. And anyway, for the polys, all of SymPy uses it, so you
really do have to run the whole tests to see.
Aaron Meurer
>
> Oscar
>
> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To view this discussion on the web visit
Reply all
Reply to author
Forward
0 new messages