Same Matrix Construction gives two different results

[Sihuang Hu]

I just used the following commands to construct an matrix:

q = 2
r = 2
s = 2
n = s*(r+1)
d = 3
P = matrix(ZZ, r+2, r+2, lambda u, k: Krawtchouk(r+1, q, k, u))

But sage gives me two different results at different times:

[ 1  3  3  1]
[ 1  1 -1 -1]
[ 1 -1 -1  1]
[ 1 -3  3 -1]

and

[ 1  3  3  1]
[ 1  0 -3 -1]
[ 1 -3  3  1]
[ 1 -3  3 -5]

The upper matrix is the right one. But sometimes I get the lower matrix.
Anyone knows why this happens? Thanks :)

[Dima Pasechnik]

No idea (unless you create another function named Krawtchouk, or something along these lines, e.g. a function used in the it). The whole implementation is

def Krawtchouk(n,q,l,i):
    """
    Compute ``K^{n,q}_l(i)``, the Krawtchouk polynomial:
    see :wikipedia:`Kravchuk_polynomials`.
    It is given by

    .. math::

        K^{n,q}_l(i)=\sum_{j=0}^l (-1)^j(q-1)^{(l-j)}{i \choose j}{n-i \choose l-j}
    """
    from sage.rings.arith import binomial
    # Use the expression in equation (55) of MacWilliams & Sloane, pg 151
    # We write jth term = some_factor * (j-1)th term
    kraw = jth_term = (q-1)**l * binomial(n, l) # j=0
    for j in range(1,l+1):
        jth_term *= -q*(l-j+1)*(i-j+1)/((q-1)*j*(n-j+1))
        kraw += jth_term
    return kraw

[Jori Mäntysalo]

On Sun, 8 Nov 2015, Sihuang Hu wrote:

> I just used the following commands to construct an matrix:
>
> q = 2
> r = 2
> s = 2
> n = s*(r+1)
> d = 3
> P = matrix(ZZ, r+2, r+2, lambda u, k: Krawtchouk(r+1, q, k, u))
>
> But sage gives me two different results at different times:

I was not able to reproduce the error. I tried in two system, about
500,000 tests together. If there is a real heisenbug, it must depend on
architechture etc. So Linux or Mac, what processor...?

--
Jori Mäntysalo

[William Stein]

It could also be flaky hardware. Does your computer ever crash?

William

>
> --
> Jori Mäntysalo



--
William (http://wstein.org)

[Sihuang Hu]

Under the Sagemath cloud and Windows, I always get

[ 1  3  3  1]
[ 1  1 -1 -1]
[ 1 -1 -1  1]
[ 1 -3  3 -1]

Under the linux version, I sometimes get

[ 1  3  3  1]
[ 1  0 -3 -1]
[ 1 -3  3  1]
[ 1 -3  3 -5]

I found that when I paste the following command

q = 2
r = 2
s = 2
n = s*(r+1)
d = 3
P = matrix(ZZ, r+2, r+2, lambda u, k: Krawtchouk(r+1, q, k, u))

into sage command line, I will get

[ 1  3  3  1]
[ 1  1 -1 -1]
[ 1 -1 -1  1]
[ 1 -3  3 -1]

When I load these command from a file "***.py", I will get

[ 1  3  3  1]
[ 1  0 -3 -1]
[ 1 -3  3  1]
[ 1 -3  3 -5]

在 2015年11月9日星期一 UTC+2上午7:34:59，jori.ma...@uta.fi写道：

[Sihuang Hu]

I installed Ubuntu 14.04 LTS 64-bit last weekend. I don't remember there are any crash.

在 2015年11月9日星期一 UTC+2上午7:42:40，William写道：

[Nils Bruin]

On Sunday, November 8, 2015 at 10:53:31 PM UTC-8, Sihuang Hu wrote:

When I load these command from a file "***.py", I will get
[ 1  3  3  1]
[ 1  0 -3 -1]
[ 1 -3  3  1]
[ 1 -3  3 -5]

Ah, that's a very important detail: that indicates that something that the preprocessor does makes a difference. Indeed:

sage: Krawtchouk(3,2,3,3)
-1
sage: Krawtchouk(int(3),int(2),int(3),int(3))
-5

That's a bug, and one that should be pretty straightforward to solve.

[Nils Bruin]

On Sunday, November 8, 2015 at 11:13:24 PM UTC-8, Nils Bruin wrote:

That's a bug, and one that should be pretty straightforward to solve.

Indeed, it's exceedingly silly. In Krawtchouk we execute:

        jth_term *= -q*(l-j+1)*(i-j+1)/((q-1)*j*(n-j+1))

If we'de doing that division with Sage Integer then we'll be getting rational numbers. With python int this amounts to "//", so will likely give other answers.

Workaround: make sure you input Sage Integers (i.e., if you absolutely have to use a ".py" file, write Integer(1) instead of 1.

Fix:
-    for j in range(1,l+1):
+    for j in srange(1,l+1):


 

[Sihuang Hu]

Great. Now I load the following commands:

q = 2
r = 2
s = 2
n = s*(r+1)
d = 3

P = matrix(ZZ, r+2, r+2, lambda u, k: Krawtchouk(r+1, q, ZZ(k), ZZ(u)))

I got:
[ 1  3  3  1]
[ 1  1 -1 -1]
[ 1 -1 -1  1]
[ 1 -3  3 -1]

So this is the right one? I'm a little bit confused.

在 2015年11月9日星期一 UTC+2上午9:19:45，Nils Bruin写道：

[Jori Mäntysalo]

On Sun, 8 Nov 2015, Nils Bruin wrote:

> Ah, that's a very important detail: that indicates that something that the
> preprocessor does makes a difference.

Todo to someone: add something like
http://doc.sagemath.org/html/en/faq/faq-usage.html#what-exactly-does-sage-do-when-i-type-0-6-2
to FAQ so that in future we could just give a link as an answer.

--
Jori Mäntysalo

[Dima Pasechnik]

On Monday, 9 November 2015 07:19:45 UTC, Nils Bruin wrote:

On Sunday, November 8, 2015 at 11:13:24 PM UTC-8, Nils Bruin wrote:

That's a bug, and one that should be pretty straightforward to solve.

Indeed, it's exceedingly silly. In Krawtchouk we execute:

        jth_term *= -q*(l-j+1)*(i-j+1)/((q-1)*j*(n-j+1))

If we'de doing that division with Sage Integer then we'll be getting rational numbers. With python int this amounts to "//", so will likely give other answers.

I wrote that code a while ago, as an internal function, and the reviewer (hi Basu!) insisted that I made it public. But me and he missed the argument checking....

Dima
 

[P Purkayastha]

Hi Dima!

Yes, we clearly overlooked that bug because we also tested on Sage
itself, and never from just python. An easy fix to this would be the
following in my opinion, so that everything gets coerced to Sage
integers, irrespective of the input:

Q = Integer(q)
kraw = jth_term = (Q-1)**l * binomial(n, l) # And replace all q with
Q. everywhere else

Regards,
basu.

[Daniel Krenn]

On 2015-11-09 08:06, P Purkayastha wrote:
> [...] so that everything gets coerced to Sage

> integers, irrespective of the input:
>
> Q = Integer(q)

FYI, this is a conversion; use
Q = ZZ.coerce(q)
instead.

Daniel

[Dima Pasechnik]

I've opened http://trac.sagemath.org/ticket/19561 to deal with it.
Dima

On Monday, 9 November 2015 07:19:45 UTC, Nils Bruin wrote:

[Dima Pasechnik]

ready for review