Add MultiPolynomialIdeal.groebner_basis("libgiac")

[malb]

Giac is an optional package and supports fast computation of Gröbner bases. Let's make that easily accessible.

This trac depends on #18841 in the following points: 1) a sig_count error detection in the doctest. 2) the output of the Katsura3 doctest in multi_polynomial_ideal.py was generated with giac-1.2.0-19 as provided in #18841. (Earlier versions provides a reversed output)

[frederichan]

Thank you. NB: Currently giacpy is just an optional spkg not a sage library so there is no sage.lib.giac

I don't know if there is a suitable place on sage git to expand the so called upstream spkg, I would be happy to put it in public so that sage developper could also modify it.

[malb]

giacpy being optional is no problem. I'd suggest to add sage.libs.libgiac which imports the right stuff from giacpy and raises an error if that fails, e.g.

def groebner_basis_libgiac(gens, epsilon=0, prot=False, <other args>):
   from giacpy import libgiac
   F = libgiac(gens)
   giacsettings.proba_epsilon=epsilon
   # get ring from gens, check term order etc.
   B = F.gbasis([P.gens()],'revlex') 
   # convert result to PolynomialSequence and return

[frederichan]

Sorry I am struggling with git (I am failing to create a branch), could you please start a branch with your function so that I can pull?

NB: you need to import giacsettings also.

[frederichan]

ouf! this time I have succeed to push a trivial thing. I can start to work

---

New commits:

​6c54aac

start groebner basis branch with libgiac

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​9e02904	Merge branch 'public/giacpyGB' of trac.sagemath.org:sage into public/giacpy-groebner
​0988f8d	rename sage.libs.libgiac to sage.libs.giac
​0b695a0	poc of groebner_basis_libgiac

[malb]

I added a proof of concept implementation

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​f57602f	Merge branch 'develop' of trac.sagemath.org:sage into t/public/giacpyGB
​2fa05ba	Merge branch 'public/giacpyGB' of trac.sagemath.org:sage into t/public/giacpyGB
​fab096e	test import, gens, rings. change conversion to sage

[frederichan]

I have noticed that with giac1.2.0-13 (spkg version) it may happen that the giac answer contains a 0 term, but it is fixed in giac1.2.0-15. May be updating the spkg is the best solution?

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​55fff07	drop unused code and imports
​5f95d87	don't check for list type explicitly
​b8725bc	obey PEP8 style conventions
​86bd0a4	use decorators for preserving global settings

[malb]

I cleaned it up a bit more. It still needs tests, documentation and a way to shut Giac up. Also, does Giac have an option to print a protocol of the GB computation (most systems do).

Finally, we get this doctest failure:

File "src/sage/libs/giac.py", line 19, in sage.libs.giac
Failed example:
    sig_on_count()
Expected:
    0
Got:
    1

which suggests there's a bug somewhere.

[malb]

Updating Giac should be a separate ticket.

[frederichan]

Thank you. I am still working on this. I found a way for the info level.

[frederichan]

The sig_count pb is from giacpy 0.4.8 Pygen __iter__, it could/should? be fixed in a separate trac? we can avoid it with: return PolynomialSequence(P, tuple(gb_giac), immutable=True) instead of just gb_giac. I don't remark time loss with tuple, should we adopt this?

is there a way to substitute efficiently variables in PolynomialRing? (ex: 'i' will give use troubles in giac) if not I will just add some protection with an error message

[frederichan]

sorry I wanted to say: I don't remark time loss with tuple. (so I have edited my comment:14)

[malb]

Replying to frederichan:

The sig_count pb is from giacpy 0.4.8 Pygen __iter__, it could/should? be fixed in a separate trac?

yes, this should be a separate ticket but this ticket IMHO should depend on that ticket. That is, this ticket should wait until the fix is in.

we can avoid it with: return PolynomialSequence(P, tuple(gb_giac), immutable=True) instead of just gb_giac. I remark time loss with tuple, should we adopt this?

You mean instead of PolynomialSequence(gb_giac, P, immutable=True)? I don't think we should work around a bug which should be fixed instead.

is there a way to substitute efficiently variables in PolynomialRing? (ex: 'i' will give use troubles in giac) if not I will just add some protection with an error message

You can try but it's not very fast:

sage: P.<a,b,c,i> = QQ[]
sage: R.<x,y,z,w> = QQ[]
sage: f = P.random_element()
sage: R(f)
-2*y^2 - 1/4*y*z + 6*z*w - 2*z - 59/5
sage: %timeit R(f)
10000 loops, best of 3: 77.3 µs per loop

Maybe check P.gens_dict().keys() for problematic names and only switch over if there is a problem?

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​6d810f6	add prot and restore global proba
​f725910	Merge branch 'public/giacpyGB' of trac.sagemath.org:sage into t/public/giacpyGB
​cc81303	add prot feature + fixes in doctests
​2f65302	add and change doctests. Add a threads option

[frederichan]

giac gbasis have also a rur option to obtain a Rational Univariate Representation of a 0 dimensional ideal. Are you also interested by a more convient interface to this or should we leave it for later?

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​425b8dd

add doctests

[malb]

I'm a one-feature-per-ticket kind of guy, so I'd put that rur enhancement in a separate ticket.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​756d223	wrap long lines as recommended by style guide
​7076dd0	follow PEP8 style guide
​28e2db9	make error message more concise

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​061f277

Merge branch 'develop' of trac.sagemath.org:sage into giacpyGB

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​81b68fa	add groebner_basis("giac") option
​8bf1320	we have namespaces, use them

[frederichan]

Dear martin, thank you for the modifs. About the doctests I don't trust even the first 2 characters: // of the message // Giac ... I have not yet tested this giac spkg 1.2.0-19 on os X to see but last time I got pbs with the begining of this message on OSX versus linux. so I think that in a doctest it is safer to do it in 2 steps: sage: gb = J.groebner_basis('giac') # optional - giacpy, random sage: gb # optional - giacpy

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​5d23482

fix-up doctest

[malb]

I guess this ticket is done then (?)

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​b590a50

blacklist 'e' for variable names

[frederichan]

may be I should also add some doc in multi_polynomial_ideal.py or pointer to the sage.libs.giac.groebner_basis doc

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​5603124

add doc in sage.rings.polynomial.multi_polynomial_ideal.groebner_basis

[malb]

Perhaps

sage: I9.groebner_basis("giac",proba_epsilon=1e-7) # optional - giacpy

should be tagged long?

[frederichan]

On a 2012 (2core i5-2435M CPU @ 2.40GHz) notebook with 4Go of RAM:

sage: A9=PolynomialRing(QQ,9,'x') # optional - giacpy
sage: I9=sage.rings.ideal.Katsura(A9) # optional - giacpy
sage: time I9.groebner_basis("giac",proba_epsilon=1e-7) # optional - giacpy
// Giac share root-directory:/home/fred/dev/sage/develop/sage.run/local/share/giac/
// Unable to find keyword file /home/fred/dev/sage/develop/sage.run/local/share/giac/doc/fr/keywords
// Giac share root-directory:/home/fred/dev/sage/develop/sage.run/local/share/giac/
Help file /home/fred/dev/sage/develop/sage.run/local/share/giac/doc/fr/aide_cas not found
Added 0 synonyms
Running a probabilistic check for the reconstructed Groebner basis. If successfull, error probability is less than 1e-07 and is estimated to be less than 10^-133. Use proba_epsilon:=0 to certify (this takes more time).
CPU times: user 3.29 s, sys: 19.9 ms, total: 3.31 s
Wall time: 3.23 s
Polynomial Sequence with 143 Polynomials in 9 Variables

[malb]

Three seconds is long by sage doctest standards. We're running thousands of doctests.

[frederichan]

OK so some others in giac.py need it also. I see that sometime there is an indication of the time. would # optional giacpy, long time (3s on a 2012 notebook)

be OK or is it too long.

[malb]

It's fine but not mandatory to add this (I'd say)

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​36620a9

add long time tags to some doctests

[frederichan]

to change I have tried it on OSX(10.7) with #18841 and got no problems.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​24a7254	Merge branch 'develop' of trac.sagemath.org:sage into t/public/giacpyGB
​f3a235b	update docstring in cachefunc.pyx (groebner_basis call)

[frederichan]

The zero ideal gives a segfault

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​02cf665

add an is_zero ideal test before calling giac gbasis

[malb]

I think this ticket is done, i.e. I'm positively reviewing your contributions. If you review mine and agree with them, we can set this to positive review.

[frederichan]

OK, I thought we had to wait for a third person. I will look if I can reproduce this doctest error: ​http://patchbot.sagemath.org/log/18824/Ubuntu/14.04/x86_64/3.13.0-61-generic/librae/2015-08-17%2014:47:28%20+0000?short do you think it could be related to this patch?

[frederichan]

I am positively reviewing your changes.

I was not able to reproduce the previous patchbot doctest error with sage 6.9.beta2.