is_square and sqrt for polynomials and fraction fields

[robertwb]

Implement is_square and sqrt for polynomials and fraction fields.

Only apply: attachment: trac_9094-sqrt-mderickx.patch

CC: @koffie @pjbruin @mstreng @mminzlaff

Component: algebra

Author: Robert Bradshaw, Maarten Derickx

Reviewer: John Cremona, Marco Streng, Robert Bradshaw

Merged: sage-4.7.alpha5

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

[robertwb]

Attachment: 9094-sqrt.patch.gz

[JohnCremona]

comment:2

Looks good (and useful). Patch applies fine to 4.4.3. I tested all sage/rings and had one failure (in a new doctest, on 64-bit ubuntu):

File "/storage/jec/sage-4.4.3/devel/sage-tests/sage/rings/polynomial/polynomial_element.pyx", line 1262:
    sage: (9 * f^10 * g^4).sqrt() == 3 * f^5 * g^2
Expected:
    True
Got:
    False

I have never been happy about doctests using random elements, however deterministic they are supposed to be.

I would have tested for self==0 at the start of the sqrt function, but that case will be caught quickly anyway, so no problem.

I have been working on a patch to add both all= and extend= parameters to the sqrt function for AA and QQbar, having found how awkward it is when the parameters are not uniform across types. Would it be worth adding an extend= parameter here, even if a NotImplementedError is raised when it is needed (i.e. for sqrt of a non-square when extend=True)?

So: needs work because of the doctest failure; up to you whether to do the extend= thing.

[JohnCremona]

Work Issues: doctest failure

[JohnCremona]

Reviewer: John Cremona

[JohnCremona]

Author: Robert Bradshaw

[pjbruin]

comment:5

is_square and sqrt can be done generally for fraction fields, given the corresponding functions for the base ring:

is_square (a / b) = is_square (a * b)
sqrt (a / b) = sqrt (a * b) / b

[koffie]

comment:7

On sage days 23 we discussed how to make this faster by avoiding factoring in polynomial rings. We are currently implementing this together with the remark of pbruin. When finnished the patch will be added to 9054 because it also fixes some other bugs on that ticket.

[koffie]

comment:9

This patch also leaks memory, probably because of bug #9129

t=get_memory_usage()
Pr.<x>=ZZ[]
for i in range(500):
    C=((x^2+1)*x+1)
    B=C^2
    print "memusage", get_memory_usage(t)
    time x=B.sqrt()

executing the above gives:

memusage 0.0
Time: CPU 0.00 s, Wall: 0.00 s
memusage 0.0
Time: CPU 0.00 s, Wall: 0.00 s
memusage 0.0
Time: CPU 0.00 s, Wall: 0.00 s
memusage 0.0
Time: CPU 0.00 s, Wall: 0.00 s
memusage 0.0
Time: CPU 0.02 s, Wall: 0.02 s
memusage 1.20703125
Time: CPU 0.18 s, Wall: 0.19 s
memusage 23.79296875
Time: CPU 2.33 s, Wall: 2.36 s
memusage 148.12890625
Time: CPU 31.70 s, Wall: 32.24 s
memusage 1534.01171875
^C
^C

Ps. if you want to check this, please don't let the loop run as long as I did ;). It might make your computer on the edge of crashing

[koffie]

comment:10

I just added a patch which should fix this also. It does it in a slightly different way. The sqrt() functions just are using the is_square(root=True) functionality (one should not have two square root finding algorithms in parralel). Also the sqrt() function has now the optional parameter extend. The is_square() function now also uses the trick mentioned by Peter the Bruin which is very general as it should be in this general base class.

The patch should be installed aplied to 4.4.4 instead of the other patch.

[mstreng]

Changed work issues from doctest failure to improve documentation, add more tests

[mstreng]

Changed reviewer from John Cremona to mstreng

[mstreng]

Changed author from Robert Bradshaw to mderickx

[mstreng]

comment:11

Looks good. No time for a complete review right now, but here are some remarks:

• Check spelling of your documentation, you probably don't mean "wether" (=castrated sheep), but "whether" (spell checker doesn't help that much here). "Requiered" has one "e" too much.

• Illustrate is_square(root = True) for the user by providing a square example and a non-square example.

• The internal variable name is_sqrt is confusing, remove the "t".

• You write "This code is quite general, it could even be implemented for all integral domains as soon as they have the is_square(root=True) option", don't you mean "This code is quite general, it works for all integral domains that have the is_square(root = True) option"?

• I think sqrt(self, extend = True, all = False, name=None ) should have tests of all combinations of options and squareness, i.e: sq/non-sq, extend True/False, all True/False = 222 = at least 8 examples. The extend=True for a non-square should also have examples with name=None and name=something.

• Perhaps include examples for more base fields?

• consider adding a few more spaces for readability thougout, e.g. around "="

• "if root==True:", why not write "if root:"?

[koffie]

comment:23

Only apply: trac_9094-sqrt-mderickx.2.patch

[mminzlaff]

comment:25

Since you told me you're going to rewrite the patch, I'm holding off on doctesting. Three short comments:

in structure/element.pyx

1. Line 1920 (): true --> True

2. Line 2002 (): all=False = all=True

3. The memory leak (from comment #9) is no leak after all. Note that you replace x in each iteration by C. If one doesn't do that, memory usage in each iteration is 0.0. :)

[robertwb]

comment:27

• is_square for fraction field elements now always computes a root, even if it's not needed. This can be quite expensive.

  • Why are you writing boilerplate wrappers for all of LazyNamedUnop? Just do is_square and sqrt if need be.

  • Using _parent is just fine, especially for subclasses. No need to incur the extra expense of a method call.

[mstreng]

comment:28

Replying to @robertwb:

• is_square for fraction field elements now always computes a root, even if it's not needed. This can be quite expensive.

Indeed, it should use root=root instead of root=True

• Why are you writing boilerplate wrappers for all of LazyNamedUnop? Just do is_square and sqrt if need be.

Hi Robert,

I can answer this question. This was a discussion on sage-devel. 4 options were discussed for dealing with these unitary operators for lazy numbers.

A) leaving them as they were

B) hard-coding them (as Maarten just did)

C) "fixing LazyNamedUnop to preserve documentation"

D) "populating these methods at class creation time
(rather than attribute lookup time, perhaps dynamically creating a
docstring for them)"

Option A was discarded: it lacks tab completion, and (as we saw in this ticket) leads to a big mess and a lot of confusion when a base class such as RingElement gets a default implementation for one of these unitary operators.

Options C and D were suggested by you, but you seemed to be the only one in the discussion who knew how to do them. Maarten said he was interested in them, and then nothing happened for months.

So I guess that in the end Maarten implemented the only remaining option, which leads to stable, transparent, understandable code that is unfortunately a bit long.

I think the extra time that Sage developers will need when extending or changing Maarten's long list of similar methods is not too bad if you compare it to the time that Sage developers will need when trying to understand more complicated code and having to debug it if something unforeseen happens. But that's just my opinion.

+1 from me for the way Maarten implemented it.

• Using _parent is just fine, especially for subclasses. No need to incur the extra expense of a method call.

Are there general guidelines for this? I can imagine _parent to be faster, but parent() to be more stable in case creative implementations are in base classes.

[koffie]

comment:29

New file added Only apply: trac_9094-sqrt-mderickx.patch

It's RLF is now boiler plate code free, works with tab completion and works with inheritance thanks to Thomas Kluiver here at sage days since he told me how to do it :).

Also added the other suggestions. Hopefully finally a final version??

[mstreng]

comment:30

Hi Maarten, could you add a doctest to __dir__ so that doctest coverage goes up instead of down? Just something simple like

sage: "log" in RLF(sqrt(8)).__dir__()
True

[koffie]

A patch instead of the other one, which I ofcourse think is better ;)

[koffie]

comment:31

Attachment: trac_9094-sqrt-mderickx.patch.gz

Done :)

[robertwb]

comment:32

Looks good.

[mstreng]

Changed reviewer from John Cremona, Marco Streng to John Cremona, Marco Streng, Robert Bradshaw

[jdemeyer]

Changed work issues from make general sqrt compatible with lazy reals, or go back to double code to none

[jdemeyer]

Merged: sage-4.7.alpha5