Opened 13 years ago

Closed 13 years ago

Last modified 13 years ago

#6237 closed defect (fixed)

repeated roots with roots(CDF, multiplicities=False)

Reported by: ncalexan Owned by: somebody
Priority: major Milestone: sage-4.3.1
Component: basic arithmetic Keywords: roots CDF multiplicities
Cc: Merged in: sage-4.3.1.alpha2
Authors: Alex Ghitza Reviewers: Karl-Dieter Crisman
Report Upstream: N/A Work issues:
Branch: Commit:
Dependencies: Stopgaps:

Status badges

Description

sage: pari('v')
v
sage: pari('u')
u
sage: u = QQ['u'].0
sage: v = QQ['u']['v'].0
sage: f = v^3 - u^7 + 2*u^3*v
sage: f.discriminant()
-27*u^14 - 32*u^9
sage: f.discriminant().roots(CDF, multiplicities=False)

[-1.03456371594,
 0,
 0,
 0,
 0,
 0,
 0,
 0,
 0,
 0,
 -0.31969776999 - 0.983928563571*I,
 -0.31969776999 + 0.983928563571*I,
 0.836979627962 - 0.608101294789*I,
 0.836979627962 + 0.608101294789*I]

Note the repetition of 0.

Attachments (1)

trac_6237.patch (2.1 KB) - added by AlexGhitza 13 years ago.

Download all attachments as: .zip

Change History (10)

comment:1 Changed 13 years ago by AlexGhitza

  • Authors set to Alex Ghitza
  • Report Upstream set to N/A
  • Status changed from new to needs_review

And if you run it with just roots(CDF) you get nine different occurrences of (0, 1), i.e. root 0 with multiplicity 1.

The attached patch fixes this.

comment:2 Changed 13 years ago by kcrisman

  • Status changed from needs_review to needs_info

This is fine overall, assuming the answer to the following question is yes.

            if output_complex:
                rts = sort_complex_numbers_for_display([L(root) for root in ext_rts])
            else:
                rts = [L(root.real()) for root in ext_rts if root.imag() == 0]

The first list gives a canonical ordering, so using

                rts_mult.append((rt, mult))
                j += mult

is okay, it won't append the same thing twice. Is that also true for the second list? I couldn't come up with an example that breaks it, but that just means I don't know much about how polynomials are represented internally. At any rate, it should probably be ordered, just to be on the safe side, if you use that way of finding multiplicities. Does Pari not compute multiplicities for this? Maxima's implementation does, though perhaps it doesn't do arbitrary precision. I assume numpy definitely doesn't do multiplicities.

comment:3 Changed 13 years ago by AlexGhitza

  • Status changed from needs_info to needs_work

Thanks for catching this.

From Pari's documentation for the function we're using:

polroots(pol,{flag = 0})

complex roots of the polynomial pol, given as a column vector where each root is repeated according to its multiplicity. [...]

The algorithm used is a modification of A.Sch¨nhage's root-finding algorithm, due to and implemented by X.Gourdon. Barring bugs, it is guaranteed to converge and to give the roots to the required accuracy.

There is no mention of the roots being sorted. I guess I could read the source code and find out, but I like you suggestion of sorting the Pari output anyway -- just in case they change the behaviour in the future.

Also from the above snippet, Pari indeed does not give the multiplicities, it just repeats each root the correct number of times.

I will replace the patch soon.

Changed 13 years ago by AlexGhitza

comment:4 Changed 13 years ago by AlexGhitza

  • Status changed from needs_work to needs_review

OK, the revised patch is up.

comment:5 Changed 13 years ago by kcrisman

  • Reviewers set to Karl-Dieter Crisman
  • Status changed from needs_review to positive_review

Looks good. Passes relevant tests.

I did find something else weird, perhaps just my install is somehow corrupt... Did you try the original thing from the description as well? It seems to be broken at the calculation of the discriminant both before and after the patch. I get an error message about PariError?(8). However, this does not appear on sage.math, so I assume something weird happened in my local install. But I thought I'd mention it in case you find it.

Anyway, positive review!

comment:6 follow-up: Changed 13 years ago by AlexGhitza

Did you try the whole thing, including pari('v') and pari('u') ?

It seems to work for me. Of course, if I had written the example I would have tried something like

sage: R.<u> = QQ[]
sage: S.<v> = R[]
sage: f = v^3 - u^7 + 2*u^3*v
sage: f.discriminant()

This does indeed fail with a PariError(8). And this is apparently documented in the docstring for discriminant. It's a shame that it doesn't work out of the box. In fact, I would even like something like the following to work:

sage: R.<u, v> = QQ[]
sage: f = v^3 - u^7 + 2*u^3*v
sage: f.discriminant(v)

Anyway, it's beyond the scope of this ticket. I might open an enhancement ticket about it.

comment:7 in reply to: ↑ 6 Changed 13 years ago by kcrisman

Replying to AlexGhitza:

Did you try the whole thing, including pari('v') and pari('u') ?

Yes. Again, oddly enough, the same thing did work when I tried it on sage.math, so it must be highly sensitive to something, though not sure what - maybe I had typed in something earlier that made it work/not work.

It seems to work for me. Of course, if I had written the example I would have tried something like

sage: R.<u> = QQ[]
sage: S.<v> = R[]
sage: f = v^3 - u^7 + 2*u^3*v
sage: f.discriminant()

This does indeed fail with a PariError(8). And this is apparently documented in the docstring for discriminant.

Huh. Well, glad to know it.

Anyway, it's beyond the scope of this ticket. I might open an enhancement ticket about it.

Go for it, though I won't be able to help on it.

comment:8 Changed 13 years ago by rlm

  • Merged in set to 4.3.1.alpha2
  • Resolution set to fixed
  • Status changed from positive_review to closed

comment:9 Changed 13 years ago by mvngu

  • Merged in changed from 4.3.1.alpha2 to sage-4.3.1.alpha2
Note: See TracTickets for help on using tickets.