roots(algorithm='numpy') does not work in arbitrary precision

[zimmermann6]

Consider the following example:

sage: R.<x> = PolynomialRing(ComplexField(3322))
sage: p=x^4+54*x^2+154
sage: z=p.roots(algorithm='pari')
sage: e=p-mul([x-z[i][0] for i in range(4)])
sage: n(max(abs(e.coeffs()[i]) for i in range(0,e.degree()+1)))
6.08883742371831e-999

This is ok. Compare now with:

sage: R.<x> = PolynomialRing(ComplexField(3322))
sage: p=x^4+54*x^2+154
sage: z=p.roots(algorithm='numpy')
sage: e=p-mul([x-z[i][0] for i in range(4)])
sage: n(max(abs(e.coeffs()[i]) for i in range(0,e.degree()+1)))
6.06533797844328e-14

Clearly the precision given by numpy is only 14 digits, not 1000
digits as expected.

CC: @rbeezer

Component: numerical

Author: Mike Hansen

Reviewer: Paul Zimmermann

Merged: sage-4.6.alpha2

Issue created by migration from https://trac.sagemath.org/ticket/8054

[jasongrout]

comment:1

Numpy does not do arbitrary precision things. So I suppose we should just document this.

[zimmermann6]

comment:2

Numpy does not do arbitrary precision things. So I suppose we should just document this.

still not done in 4.3.1.

[mwhansen]

comment:3

Attachment: trac_8054.patch.gz

[mwhansen]

Author: Mike Hansen

[zimmermann6]

comment:4

Mike, the warning message is just printed once, i.e., if one calls roots(algorithm='numpy') twice,
it is not printed the second time. Is it wanted?

Paul

[zimmermann6]

Reviewer: Paul Zimmermann

[zimmermann6]

comment:5

All test pass (with Sage 4.4.4). Thus I give a positive review.

[qed777]

comment:6

Could someone rebase the patch against 4.6.alpha1 when it's released (soon) and restore the positive review?

[qed777]

Work Issues: rebas

[jasongrout]

comment:7

If someone is rebasing it, they might also correct a typo:

Note that one should not use NumPy when wanting high precicion -> precision

[zimmermann6]

comment:8

by the way, is there a documentation somewhere explaining how to rebase a patch?

Paul

[qed777]

Changed work issues from rebas to rebase

[qed777]

comment:9

Replying to @zimmermann6:

by the way, is there a documentation somewhere explaining how to rebase a patch?

One way is to use Mercurial queues:

$ cd SAGE_ROOT/devel/sage
$ hg qimport https://github.com/sagemath/sage-prod/files/10647755/trac_8054.patch.gz
adding trac_8054.patch to series file
$ hg qpush
applying trac_8054.patch
patching file sage/rings/polynomial/polynomial_element.pyx
Hunk #2 FAILED at 4256
1 out of 4 hunks FAILED -- saving rejects to file sage/rings/polynomial/polynomial_element.pyx.rej
patch failed, unable to continue (try -v)
patch failed, rejects left in working dir
errors during apply, please fix and refresh trac_8054.patch

The file polynomial_element.pyx.rej is a diff of the changes that Mercurial was unable to apply. Editing polynomial_element.pyx by hand to incorporate the leftover changes and then doing

$ hg qrefresh
$ hg export qtip > /path/to/trac_8054.2.patch
$ hg qpop
$ hg qdelete trac_8054.patch

should write out an updated patch and restore the original state of the repository. If you choose not to delete the patch from the queue, e.g., for doctesting, then for convenience you can use

$ hg qrename trac_8054.2.patch

to rename it in the queue. After this, hg qseries, hg qapplied, and hg qunapplied will use the updated name (until you qdelete the patch).

There's more on using queues in the Developer's Guide, but as far as I can tell, there's no mention of rebasing patches. Let us know if you have questions.

[zimmermann6]

comment:10

Replying to @qed777:

thank you. I will try when 4.6.alpha1 is available. Do you know when?

Paul

[qed777]

comment:11

Replying to @zimmermann6:

Replying to @qed777:
thank you. I will try when 4.6.alpha1 is available. Do you know when?

No problem. I hope I haven't made any mistakes!

I very recently announced 4.6.alpha1 on sage-release.

I've cc'd Rob Beezer about rebasing patches.

[zimmermann6]

Attachment: trac_8054.2.patch.gz

original patch rebased on 4.6.alpha1

[zimmermann6]

comment:12

Replying to @qed777:

Could someone rebase the patch against 4.6.alpha1 when it's released (soon) and restore the positive review?

done.

Paul

[rbeezer]

comment:13

I'll try to send Paul some rebase instructions that he might look over and then and maybe grow those into something for the Developer's Guide.

Rob

[qed777]

comment:14

The rebased patch applies cleanly to 4.6.alpha1. But testing just the changed file, I get

sage -t -long "devel/sage-main/sage/rings/polynomial/polynomial_element.pyx"
**********************************************************************
Error: TAB character found.

         [13.2 s]
 
----------------------------------------------------------------------
The following tests failed:


        sage -t -long "devel/sage-main/sage/rings/polynomial/polynomial_element.pyx"
Total time for all tests: 13.2 seconds

There's a stray tab in this line

+       ``algorithm='numpy'`` with high-precision types.)                       

(in the patch).

[qed777]

Changed work issues from rebase to none

[qed777]

comment:15

Replying to @rbeezer:

I'll try to send Paul some rebase instructions that he might look over and then and maybe grow those into something for the Developer's Guide.

Thanks, Rob. I wasn't sure whether this is in the Guide. I should reread it in full, soon.

[zimmermann6]

Attachment: trac_8054.3.patch.gz

[zimmermann6]

comment:16

sorry, I've replaced the tab. Apply only the last patch.

Paul

[qed777]

Merged: sage-4.6.alpha2

[qed777]

comment:18

Replying to @qed777:

Replying to @zimmermann6:

by the way, is there a documentation somewhere explaining how to rebase a patch?

One way is to use Mercurial queues:

Another possible way:

$ hg import --no-commit filename.patch

which just updates the working directory. Then you can check the .rej file, make changes, re-commit, and export anew.