Opened 4 years ago
Closed 19 months ago
#20145 closed defect (fixed)
Hilbert series bug
Reported by:  stumpc5  Owned by:  

Priority:  major  Milestone:  sage8.4 
Component:  commutative algebra  Keywords:  Hilbert series, polynomial ring 
Cc:  Merged in:  
Authors:  Simon King  Reviewers:  Travis Scrimshaw 
Report Upstream:  Reported upstream. Developers acknowledge bug.  Work issues:  
Branch:  a1545cc (Commits)  Commit:  a1545cc6c843fcdb85cbca3c85ca8ce49de28616 
Dependencies:  #26243  Stopgaps: 
Description (last modified by )
UPSTREAM: https://www.singular.unikl.de:8005/trac/ticket/750
This works correctly:
sage: n=4;m=10;P = PolynomialRing(QQ,n*m,"x"); x = P.gens(); M = Matrix(n,x) sage: M [ x0 x1 x2 x3 x4 x5 x6 x7 x8 x9] [x10 x11 x12 x13 x14 x15 x16 x17 x18 x19] [x20 x21 x22 x23 x24 x25 x26 x27 x28 x29] [x30 x31 x32 x33 x34 x35 x36 x37 x38 x39] sage: I = P.ideal(M.minors(2)) sage: I.hilbert_series().numerator() 84*t^3  108*t^2  27*t  1 sage: factor(I.hilbert_series().denominator()) (t  1)^13
but if I increase the number of columns by one, I get
sage: n=4;m=11;P = PolynomialRing(QQ,n*m,"x"); x = P.gens(); M = Matrix(n,x) sage: M [ x0 x1 x2 x3 x4 x5 x6 x7 x8 x9 x10] [x11 x12 x13 x14 x15 x16 x17 x18 x19 x20 x21] [x22 x23 x24 x25 x26 x27 x28 x29 x30 x31 x32] [x33 x34 x35 x36 x37 x38 x39 x40 x41 x42 x43] sage: I = P.ideal(M.minors(2)) sage: I.hilbert_series().numerator() 120*t^33  3465*t^32 + 48180*t^31  429374*t^30 + 2753520*t^29  13522410*t^28 + 52832780*t^27  168384150*t^26 + 445188744*t^25  987193350*t^24 + 1847488500*t^23 + 1372406746*t^22  403422496*t^21  8403314*t^20  471656596*t^19 + 1806623746*t^18 + 752776200*t^17 + 752776200*t^16  1580830020*t^15 + 1673936550*t^14  1294246800*t^13 + 786893250*t^12  382391100*t^11 + 146679390*t^10  42299400*t^9 + 7837830*t^8  172260*t^7  468930*t^6 + 183744*t^5  39270*t^4 + 5060*t^3  330*t^2 + 1 sage: factor(I.hilbert_series().denominator()) (t  1)^44
which is clearly not correct, as it can be checked with
sage: from sage.libs.singular.function import singular_function sage: hilb = singular_function("hilb") sage: std = singular_function("std") sage: hilb(std(I)) // ** _ is no standard basis ... // 1 t^0 // 30 t^1 // 135 t^2 // 120 t^3 // dimension (proj.) = 13 // degree (proj.) = 286
Change History (66)
comment:1 Changed 4 years ago by
comment:2 Changed 4 years ago by
The correct answer, as given in the second computation hilb(std(I))
, is
( 120*t^3 + 135*t^2 + 30*t + 1 ) / (t1)^14
which is not what Sage outputs. (But the total degree 11
is still correct...)
comment:3 followups: ↓ 4 ↓ 5 Changed 4 years ago by
well, for m=10 you only get the output above due to the fact that numerator has a factor (t1)^{27} that gets cancelled.
And for m=11 the numerator is irreducible, so no cancellation happens. Does this mean that it's only due to this cancellation accident in the m=10 case you get the same output in Sage as you expect?
Numerator is computed by Singular and so it should not give any bug, I guess...
comment:4 in reply to: ↑ 3 ; followup: ↓ 6 Changed 4 years ago by
Replying to dimpase:
well, for m=10 you only get the output above due to the fact that numerator has a factor (t1)^{27} that gets cancelled.
I am not sure I follow what you mean  in the m=10 case, the (afaik) correct Hilbert series is
( 84*t^3 + 108*t^2 + 27*t + 1 ) / (1  t)^13
which is what I get above. But in the m=11 case, the (afaik) correct Hilbert series is
( 120*t^3 + 135*t^2 + 30*t + 1 ) / (1  t)^14
which is not what I get above. In particular, the numerator output of singular should be divisible by (1t)^30
, but  as you write  it is irreducible instead.
comment:5 in reply to: ↑ 3 Changed 4 years ago by
Replying to dimpase:
Numerator is computed by Singular and so it should not give any bug, I guess...
At least luisfe on sagedevel (https://groups.google.com/d/msg/sagedevel/tyTDd5dEDGw/8FxH9hUzBwAJ) seems to confirm this to be a bug in Singular
...
comment:6 in reply to: ↑ 4 ; followup: ↓ 7 Changed 4 years ago by
Replying to stumpc5:
Replying to dimpase:
well, for m=10 you only get the output above due to the fact that numerator has a factor (t1)^{27} that gets cancelled.
I am not sure I follow what you mean
You can read the code of I.hilbert_series
to see that it takes I.hilbert_numerator()
, which is computed by Singular (again, read the code), and divides it by (t1)^{P.ngens()}.
 in the m=10 case, the (afaik) correct Hilbert series is
( 84*t^3 + 108*t^2 + 27*t + 1 ) / (1  t)^13which is what I get above. But in the m=11 case, the (afaik) correct Hilbert series is
( 120*t^3 + 135*t^2 + 30*t + 1 ) / (1  t)^14which is not what I get above. In particular, the numerator output of singular should be divisible by
(1t)^30
, but  as you write  it is irreducible instead.
So it is a bug in Singular, indeed.
comment:7 in reply to: ↑ 6 ; followup: ↓ 8 Changed 4 years ago by
Replying to dimpase:
Replying to stumpc5:
Replying to dimpase:
You can read the code of
I.hilbert_series
to see that it takesI.hilbert_numerator()
, which is computed by Singular (again, read the code), and divides it by (t1)^{P.ngens()}.
I actually did read the code, but
 I only see "Singular is doing the computation" in
sage.libs.singular.function_factory.ff.hilb
,  the output is then going to Sage in a way I don't know, and finally
 there is a polynomial division.
So a priori, the bug can be in either of these parts.
Anyway, it seems to be a bug in Singular, so I suggest to leave this ticket open until that is fixed, and then add a doctest for checking...
I would also appreciate if someone knowing how this Hilbert numerator is computed in Singular would send a bug report there (or shows me how to reproduce the bug in Singular).
comment:8 in reply to: ↑ 7 ; followup: ↓ 9 Changed 4 years ago by
 Report Upstream changed from N/A to Not yet reported upstream; Will do shortly.
Replying to stumpc5:
Anyway, it seems to be a bug in Singular, so I suggest to leave this ticket open until that is fixed, and then add a doctest for checking...
I would also appreciate if someone knowing how this Hilbert numerator is computed in Singular would send a bug report there (or shows me how to reproduce the bug in Singular).
I think you know that in the output of hilb(std(I)) (the part you omitted and indicated by ...) is the numerator that should be divisible by the appropriate power of t1, but it is not.
As to writing up the corresponding Singular commands to show this... I can do this for you, although I didn't touch it for 5+ years...
comment:9 in reply to: ↑ 8 ; followup: ↓ 10 Changed 4 years ago by
Replying to dimpase:
Replying to stumpc5: I think you know that in the output of hilb(std(I)) (the part you omitted and indicated by ...) is the numerator that should be divisible by the appropriate power of t1, but it is not.
Oh, I didn't know, thanks. If that is the not correct Singular output, how is the stuff below
// 1 t^0 // 30 t^1 // 135 t^2 // 120 t^3 // dimension (proj.) = 13
then again correct (saying that it gives the numerator divided by the maximal power of (t1)
?
comment:10 in reply to: ↑ 9 ; followup: ↓ 11 Changed 4 years ago by
Replying to stumpc5:
Replying to dimpase:
Replying to stumpc5: I think you know that in the output of hilb(std(I)) (the part you omitted and indicated by ...) is the numerator that should be divisible by the appropriate power of t1, but it is not.
Oh, I didn't know, thanks. If that is the not correct Singular output, how is the stuff below
// 1 t^0 // 30 t^1 // 135 t^2 // 120 t^3 // dimension (proj.) = 13then again correct (saying that it gives the numerator divided by the maximal power of
(t1)
?
It's all a bit Greek to me now, and https://www.singular.unikl.de/Manual/latest/sing_310.htm#SEC349 does not help much.
Perhaps there is more than one way to compute this info, and Singular does something leading to the correct answer, even though the data it supplies to Sage is not correct.
You seem to know a reference for the correct answers to this question for all m, no? This would be handy in the bug report upstream.
comment:11 in reply to: ↑ 10 Changed 4 years ago by
Replying to dimpase:
You seem to know a reference for the correct answers to this question for all m, no? This would be handy in the bug report upstream.
This Hilbert series is known (though maybe implicit as some h vector of a simplicial complex), but I don't have a reference at hand. I will try to get one.
The reason it seemed suspicious is that I conjecture a term order for which its initial ideal is the StanleyReisner ideal of a simplicial sphere I was studying for completely different reasons. And this example was spit out to contradict that conjecture. But, as often, it was the code that's buggy and and not the conjectural description. That's anyway not yet in a form to be made public.
comment:12 Changed 4 years ago by
One would try to compute the same with Macaulay2. Did you try this?
comment:13 Changed 4 years ago by
Yes, I did  here it is the output again (with a few long outputs ...
'ed):
i31 : R = QQ[x11,x12,x13,x14,x15,x16,x17,x18,x19,x110,x111,x21,x22,x23,x24,x25,x26,x27,x28,x29,x210,x211,x31,x32,x33,x34,x35,x36,x37,x38,x39,x310,x311,x41,x42,x43,x44,x45,x46,x47,x48,x49,x410,x411] o31 = R o31 : PolynomialRing i32 : M = matrix{{x11,x12,x13,x14,x15,x16,x17,x18,x19,x110,x111},{x21,x22,x23,x24,x25,x26,x27,x28,x29,x210,x211},{x31,x32,x33,x34,x35,x36,x37,x38,x39,x310,x311},{x41,x42,x43,x44,x45,x46,x47,x48,x49,x410,x411}} o32 =  x11 x12 x13 x14 x15 x16 x17 x18 x19 x110 x111   x21 x22 x23 x24 x25 x26 x27 x28 x29 x210 x211   x31 x32 x33 x34 x35 x36 x37 x38 x39 x310 x311   x41 x42 x43 x44 x45 x46 x47 x48 x49 x410 x411  4 11 o32 : Matrix R < R i33 : I = minors(2,M) o33 = ... o33 : Ideal of R i34 : s = hilbertSeries (R/I) o34 = ... o34 : Expression of class Divide i35 : denominator s 44 o35 = (1  T) o35 : Expression of class Product i36 : n = numerator s o36 = ... o36 : ZZ[T] i37 : factor n 30 2 3 o37 = ( 1 + T) (1 + 30T + 135T + 120T ) o37 : Expression of class Product
comment:14 Changed 4 years ago by
by the way, here is what they mean by the 1st and by the 2nd Hilbert series. https://www.singular.unikl.de/DEMO_HTML/Examples/Hilbert/index.html
Indeed, this alone seems enough for a bug report. (save the Singular code :))
comment:15 Changed 4 years ago by
Okay, I am sending a bug report to Singular now...
comment:16 Changed 4 years ago by
any progress on this?
comment:17 Changed 4 years ago by
You find the bug report from Feb at https://www.singular.unikl.de:8005/trac/ticket/750, but no one seems to have looked at it so far...
comment:18 Changed 4 years ago by
 Report Upstream changed from Not yet reported upstream; Will do shortly. to Reported upstream. No feedback yet.
comment:19 Changed 2 years ago by
As rechecked at http://www.singular.unikl.de:8002/trac/ticket/750#comment:1, this appears to be an interface bug, and no bug directly in singular.
comment:20 Changed 2 years ago by
Looks like an overflow error
sage: QQ['z'](hilb(I, 2))*(1z)**30QQ['z'](hilb(I, 1)) // ** _ is no standard basis // ** _ is no standard basis 4294967296*z^22 + 4294967296*z^21  4294967296*z^20 + 4294967296*z^19  4294967296*z^18 sage: factor(4294967296) 2^32
comment:21 Changed 2 years ago by
On the singular ticket they observe that what is shown here as the 2nd Hilbert series is actually their 1st Hilbert series.
comment:22 Changed 2 years ago by
But the long singular answer should factorise with a factor (1t)^30
, no ? And it does not!
comment:23 Changed 2 years ago by
So (I double checked) this is a bug in singular indeed. Someone should reopen http://www.singular.unikl.de:8002/trac/ticket/750
comment:24 Changed 2 years ago by
More precisely, the difference between the good result and the singular result is
(2^32) * t^18 * (t^4  t^3 + t^2  t + 1)
in singular 4.1.0 (which is the version used by sage and available online)
comment:25 Changed 2 years ago by
Okay, I can follow your computation. Is this problem now just a coincidence, or how should it be related to the bug I reported?
The bug I had was that the Sage Hilbert series is the singular second Hilbert series. This is given correctly by the computation 4x10, while the (as you say wrongly computed) first Hilbert series was given for the 4x11 example.
comment:26 Changed 2 years ago by
Sage only ever calls the "first Hilbert series". If the result of singular for hilb(I, 1) was correct, the rational function in sage would simplify and its numerator would be the same result as the "second Hilbert series" of singular. *The problem is a bug in Singular*.
Now we could "fix" this bug by calling instead hilb(I, 2) and changing the sage code, but I would not call that a good solution.
comment:27 Changed 2 years ago by
oh, I see, sorry for my ignorance.
comment:28 Changed 2 years ago by
 Report Upstream changed from Reported upstream. No feedback yet. to Reported upstream. Developers acknowledge bug.
https://www.singular.unikl.de:8005/trac/ticket/750 has been reopened
comment:29 Changed 2 years ago by
 Description modified (diff)
comment:30 followup: ↓ 34 Changed 19 months ago by
Note that in the current version 4.1.1p2.p0 of Singular, there is not a wrong result due to integer overflow, but an honest "int overflow error" is raised.
In my computations of group cohomology, I suffered from these overflow errors as well, which is why I opened #26243 and provided an implementation of Hilbert series computations that goes beyond Singular's limitations.
My implementation is not much slower than Singular in the examples in which Singular works. However, in #26243, I do not propose to use my implementation as default way to compute Hilbert series (in particular, I didn't touch the .hilbert_numerator()
method of multipolynomial ideals). However, if you consider to use it in the .hilbert_*()
methods by default in order to fix *this* ticket: Go ahead!
EDIT: The failing example from here became doctest for the code from #26243.
comment:31 followup: ↓ 32 Changed 19 months ago by
comment:32 in reply to: ↑ 31 Changed 19 months ago by
comment:33 Changed 19 months ago by
 Dependencies set to #26243
comment:34 in reply to: ↑ 30 Changed 19 months ago by
Replying to SimonKing:
However, in #26243, I do not propose to use my implementation as default way to compute Hilbert series (in particular, I didn't touch the
.hilbert_numerator()
method of multipolynomial ideals). However, if you consider to use it in the.hilbert_*()
methods by default in order to fix *this* ticket: Go ahead!
+1 to this. I may work on it this coming week if no one else volunteers, but be advised that various obligations mean I'm slow.
comment:35 followup: ↓ 36 Changed 19 months ago by
If I see that correctly, the following functions need to be taken care of (all in sage.rings.polynomial.multi_polynomial_ideal):
hilbert_series
: Here we can use what I calledhilbert_poincare_series
in sage.rings.hilbert (name according to both the book of Greuel and Pfister, and to the article of Bigatti).hilbert_numerator
: Here we can use what I calledfirst_hilbert_series
in sage.rings.hilbert (name according to the book of Greuel and Pfister).hilbert_polynomial
: This needs to be reimplemented.
If I understand correctly, one wouldn't provide degree weights for hilbert_polynomial
, and I think it is easy enough to follow Definition 5.1.4 in the GreuelPfister book.
comment:36 in reply to: ↑ 35 ; followup: ↓ 37 Changed 19 months ago by
Replying to SimonKing:
If I see that correctly, the following functions need to be taken care of (all in sage.rings.polynomial.multi_polynomial_ideal):
hilbert_series
: Here we can use what I calledhilbert_poincare_series
in sage.rings.hilbert (name according to both the book of Greuel and Pfister, and to the article of Bigatti).hilbert_numerator
: Here we can use what I calledfirst_hilbert_series
in sage.rings.hilbert (name according to the book of Greuel and Pfister).
Thanks. I'll aim for that.
hilbert_polynomial
: This needs to be reimplemented.
I have done that before, so I should be able to do it for this. :)
If I understand correctly, one wouldn't provide degree weights for
hilbert_polynomial
...
Yes, the Hilbert series is nonpolynomial for nonstandard weights, and I understand that finding formulas for this is an area of research that no one is actually involved in, unless someone extended the work of Massimo Caboara's Master's student about 45 years ago.
comment:37 in reply to: ↑ 36 Changed 19 months ago by
Replying to john_perry:
Replying to SimonKing:
If I see that correctly, the following functions need to be taken care of (all in sage.rings.polynomial.multi_polynomial_ideal):
hilbert_series
: Here we can use what I calledhilbert_poincare_series
in sage.rings.hilbert (name according to both the book of Greuel and Pfister, and to the article of Bigatti).hilbert_numerator
: Here we can use what I calledfirst_hilbert_series
in sage.rings.hilbert (name according to the book of Greuel and Pfister).Thanks. I'll aim for that.
hilbert_polynomial
: This needs to be reimplemented.I have done that before, so I should be able to do it for this. :)
I'm almost done and will probably submit it later today. But of course I'd be glad about a review.
comment:38 Changed 19 months ago by
 Branch set to u/SimonKing/Hilbert_series_bug
comment:39 Changed 19 months ago by
 Commit set to dcc3f1709927a1ea622e42b38543e42ed3fb0f16
 Status changed from new to needs_review
With the new commit, the stuff from #26243 is used for all Hilbert series computations. The tests in sage/rings/polynomial pass, so, I guess it is "needs review" (I trust that the patchbot will run full tests).
For the record, here is a timing that comes out of some of my cohomology computations, in a size that Singular can manage:
sage: I Ideal (a_1_0^2, a_1_0*b_1_1, b_2_2*a_1_0, b_2_3*a_1_0, b_2_1*b_1_1, a_1_0*b_3_4, a_1_0*b_3_5, a_1_0*b_3_6, b_2_1*b_2_2, b_1_1*b_3_5, b_4_10*a_1_0, b_2_1*b_3_4, b_2_1*b_3_6, b_2_2*b_3_5, b_4_10*b_1_1 + b_2_3*b_3_4 + b_2_3^2*b_1_1 + b_2_2*b_3_6, a_1_0*b_5_17, b_2_1*b_4_10 + b_2_1*b_2_3^2, b_3_4^2 + b_1_1^3*b_3_6 + b_2_2*b_1_1*b_3_4 + b_2_2*b_2_3*b_1_1^2 + b_2_2^3, b_3_4*b_3_5, b_3_4*b_3_6 + b_1_1*b_5_17 + b_2_3*b_1_1*b_3_6 + b_2_3*b_1_1*b_3_4 + b_2_2^2*b_2_3, b_3_5^2 + b_2_1*b_2_3^2, b_3_5*b_3_6, b_3_6^2 + b_2_3*b_1_1*b_3_6 + b_2_2*b_2_3^2 + c_4_11*b_1_1^2, b_2_1*b_5_17, b_4_10*b_3_4 + b_2_3*b_1_1^2*b_3_6 + b_2_3^2*b_3_4 + b_2_2*b_5_17 + b_2_2*b_2_3*b_3_6 + b_2_2*b_2_3^2*b_1_1, b_4_10*b_3_5 + b_2_3^2*b_3_5, b_4_10*b_3_6 + b_2_3*b_5_17 + b_2_3^2*b_3_4 + b_2_2*b_2_3*b_3_6 + b_2_2*c_4_11*b_1_1, b_4_10^2 + b_2_3^2*b_1_1*b_3_6 + b_2_3^4 + b_2_2*b_2_3*b_4_10 + b_2_2^2*c_4_11, b_3_4*b_5_17 + b_2_3*b_1_1*b_5_17 + b_2_3^2*b_1_1*b_3_6 + b_2_3^2*b_1_1*b_3_4 + b_2_2*b_1_1*b_5_17 + b_2_2^2*b_4_10 + c_4_11*b_1_1^4, b_3_5*b_5_17, b_3_6*b_5_17 + b_2_3^2*b_1_1*b_3_6 + b_2_2*b_2_3*b_4_10 + c_4_11*b_1_1*b_3_4 + b_2_3*c_4_11*b_1_1^2, b_4_10*b_5_17 + b_2_3^3*b_3_6 + b_2_3^3*b_3_4 + b_2_2*b_2_3*b_5_17 + b_2_2*b_2_3^2*b_3_6 + b_2_3*c_4_11*b_1_1^3 + b_2_2*c_4_11*b_3_4 + b_2_2*b_2_3*c_4_11*b_1_1, b_5_17^2 + b_2_3^3*b_1_1*b_3_6 + b_2_2*b_2_3*b_1_1*b_5_17 + b_2_2*b_2_3^2*b_1_1*b_3_6 + b_2_2*b_2_3^4 + b_2_2^2*b_2_3*b_4_10 + c_4_11*b_1_1^3*b_3_6 + b_2_3*c_4_11*b_1_1^4 + b_2_3^2*c_4_11*b_1_1^2 + b_2_2*c_4_11*b_1_1*b_3_4 + b_2_2*b_2_3*c_4_11*b_1_1^2 + b_2_2^3*c_4_11) of Multivariate Polynomial Ring in b_2_1, b_2_2, b_2_3, b_4_10, c_4_11, a_1_0, b_1_1, b_3_4, b_3_5, b_3_6, b_5_17 over Finite Field of size 2 sage: G = I.groebner_basis() # caching the Gröbner bases
Now, with the new commit:
sage: %timeit I.hilbert_polynomial() The slowest run took 9.05 times longer than the fastest. This could mean that an intermediate result is being cached. 1000 loops, best of 3: 686 µs per loop sage: I.hilbert_polynomial() 2/3*t^3 + 4*t^2 + 16/3*t + 1
Without the new commit (i.e., based on using Singular):
sage: %timeit I.hilbert_polynomial() The slowest run took 8.71 times longer than the fastest. This could mean that an intermediate result is being cached. 100 loops, best of 3: 2.33 ms per loop sage: I.hilbert_polynomial() 2/3*t^3 + 4*t^2 + 16/3*t + 1
Last 10 new commits:
7698bb6  Some code optimizations in Hilbert series computation

8ae070b  Use flint polynomial boilerplate for Hilbert series computation

d997abe  Fix a corner case in Hilbert series computation

5af3439  Fix the erroneous fix of some corner cases

138f85b  Turn some functions of polynomial.hilbert to methods of ETuple

e900222  Hilbert series: Remove unnecessary bound checks, simplify code branches

f0e5c85  Add documentation to new ETuple methods

2d29b9b  Using slices, custom median, and formatting/P#P8 tweaks to Hilbert code.

07ff2e1  Reverting use of get_exp for __getitem__.

dcc3f17  Do not use Singular for Hilbert series computations

comment:40 followup: ↓ 41 Changed 19 months ago by
Impressive speedup! Where does it come from? A better algorithm?
comment:41 in reply to: ↑ 40 Changed 19 months ago by
Replying to dimpase:
Impressive speedup! Where does it come from? A better algorithm?
I am not totally sure. We have to distinguish: hilbert_numerator
and hilbert_series
on the one hand, and hilbert_polynomial
on the other hand.
The computation of "first Hilbert series" (also known as "Hilbert numerator") and "HilbertPoincaré series" (also known as "Hilbert series") relies on a new implementation introduced at #26243. The original motivation was to allow for really big computations (>70 variables, >2000 ideal generators) that would fail in Singular because of its 32bit integer limitation. But in the end, it was also about speed: The new implementation beats Singular in medium size examples.
The algorithm for the HilbertPoincaré series isn't exactly new, but differs from what is done in Singular as follows. Let I
be a monomial ideal in a multivariate polynomial Ring R
and let m
be a monomial of degree d
. Then, I think all implementations for the computation of the HilbertPoincaré series HP(R/I)
rely on the equation HP(R/I) = HP(R/(I+m))+t^d*HP(R/(I:(m)))
. The differences lie in the choice of m. Some systems choose m as one of the generators of I and read the above equation backwards: HP(R/(I+m)) = HP(R/I)  t^d*HP(R/(I:(m)))
. Singular chooses m to be a variable that appears in at least one generator of degree >1 of I. I think CoCoA uses the gcd of two randomly chosen generators. My implementation chooses m to be the variable that appears most frequently in the generators of I, to the power of the median of all exponents of that variable in the generators of I.
I was told by John Perry that this choice is known in the literature, but I don't have a reference. Anyway, I think part of the speedup is due to good data structures and careful usage of Cython.
However, all of the above doesn't apply here. With the example of the above timing for hilbert_polynomial
, the timings for hilbert_numerator
and hilbert_series
have hardly changed. With Singular:
sage: %timeit I.hilbert_numerator() The slowest run took 4.14 times longer than the fastest. This could mean that an intermediate result is being cached. 1000 loops, best of 3: 270 µs per loop sage: %timeit I.hilbert_series() The slowest run took 5.27 times longer than the fastest. This could mean that an intermediate result is being cached. 1000 loops, best of 3: 375 µs per loop
With the new implementation:
sage: %timeit I.hilbert_numerator() 1000 loops, best of 3: 280 µs per loop sage: %timeit I.hilbert_series() The slowest run took 5.73 times longer than the fastest. This could mean that an intermediate result is being cached. 1000 loops, best of 3: 367 µs per loop
So, in this example, the change really is in the computation of hilbert_polynomial
. Since my code is just straight forward, I really don't know why the old code was so much slower.
comment:42 Changed 19 months ago by
 Commit changed from dcc3f1709927a1ea622e42b38543e42ed3fb0f16 to 631824276e86028e86931364c97b46f88e325eff
Branch pushed to git repo; I updated commit sha1. New commits:
6318242  Fix a sign in Hilbert polynomial computation

comment:43 Changed 19 months ago by
By playing with a bigger example (where I
is the ideal defined in the attachment to #26243), I found that I got the sign wrong, meaning that depending on the example I got the correct polynomial or its negative. The new commit fixes it.
With the new commit and the example from #26243, I get consistent results (i.e., the Hilbert polynomial seems to return the coefficients of the HilbertPoincaré series for all but finitely many degrees):
sage: P = PowerSeriesRing(QQ,'t',default_prec=50) sage: hp = I.hilbert_series() sage: hs = P(hp.numerator())/P(hp.denominator()) sage: hs 1 + 77*t + 436*t^2 + 1396*t^3 + 3446*t^4 + 7194*t^5 + 13376*t^6 + 22856*t^7 + 36626*t^8 + 55806*t^9 + 81644*t^10 + 115516*t^11 + 158926*t^12 + 213506*t^13 + 281016*t^14 + 363344*t^15 + 462506*t^16 + 580646*t^17 + 720036*t^18 + 883076*t^19 + 1072294*t^20 + 1290346*t^21 + 1540016*t^22 + 1824216*t^23 + 2145986*t^24 + 2508494*t^25 + 2915036*t^26 + 3369036*t^27 + 3874046*t^28 + 4433746*t^29 + 5051944*t^30 + 5732576*t^31 + 6479706*t^32 + 7297526*t^33 + 8190356*t^34 + 9162644*t^35 + 10218966*t^36 + 11364026*t^37 + 12602656*t^38 + 13939816*t^39 + 15380594*t^40 + 16930206*t^41 + 18593996*t^42 + 20377436*t^43 + 22286126*t^44 + 24325794*t^45 + 26502296*t^46 + 28821616*t^47 + 31289866*t^48 + 33913286*t^49 + O(t^50) sage: p = I.hilbert_polynomial() sage: p(t=1) # the only exception 86 sage: p(t=2) 436 sage: p(t=3) 1396 sage: p(t=4) 3446 sage: p(t=5) 7194 sage: p(t=20) 1072294 sage: p(t=40) 15380594
comment:44 Changed 19 months ago by
 Commit changed from 631824276e86028e86931364c97b46f88e325eff to ae4a2fc0526e9e9a88bf010fd9aceacdefed8b2b
Branch pushed to git repo; I updated commit sha1. New commits:
ae4a2fc  Add consistency check for Hilbert polynomial and fix a corner case

comment:45 Changed 19 months ago by
It was needed to specialcase zerodimensional ideals. Problem was: In that case, the pole order of the HilbertPoincaré series at t=1 is zero, which means one was dividing by the factorial of 1, which gave an error. In fact, if the HilbertPoincaré series is a polynomial, then all but finitely many slices of R/I are of dimension zero, which means that the Hilbert polynomial is zero.
While I was at it: In the first commit on this ticket, I thought I needed to amend the sign, but it turned out that I don't need. Just to be on the safe side, I am now asserting that the leading coefficient of the Hilbert polynomial is positive (it has to be, because if t is sufficiently large then the Hilbert polynomial is a natural number).
comment:46 followup: ↓ 47 Changed 19 months ago by
Some microoptimizations: if len(denom):
> if denom:
, you could use for n, coeff in enumerate(second_hilbert)
, and you do not need to pass a list to sum
and prod
(which is slower because it has to create and populate the list).
Some stylistic things: PEP8 spaces are my friend and can be your friend too :P, and error messages in Python are causes (no initial capital letter and no period/fullstop).
I think it would be good to still have a way to access the Singular implementation through an algorithm
keyword (of course, defaulting to 'sage'
). (Is this there, and I missed it?) That way you can implement more robust tests by comparing the results and perhaps there are some cases where Singular does win that someone may need frequently.
comment:47 in reply to: ↑ 46 ; followup: ↓ 48 Changed 19 months ago by
Replying to tscrim:
Some microoptimizations:
if len(denom):
>if denom:
,
OK. I was not sure if a factorisation object behaves the same as a list.
I think it would be good to still have a way to access the Singular implementation through an
algorithm
keyword (of course, defaulting to'sage'
). (Is this there, and I missed it?) That way you can implement more robust tests by comparing the results and perhaps there are some cases where Singular does win that someone may need frequently.
That's a very good idea! So, I should merge in the old implementation.
comment:48 in reply to: ↑ 47 Changed 19 months ago by
Replying to SimonKing:
Replying to tscrim:
Some microoptimizations:
if len(denom):
>if denom:
,OK. I was not sure if a factorisation object behaves the same as a list.
This is true in basically every case I've come across:
sage: f = factor(123521321) sage: f 7 * 11 * 997 * 1609 sage: %timeit bool(len(f)) The slowest run took 35.73 times longer than the fastest. This could mean that an intermediate result is being cached. 1000000 loops, best of 3: 334 ns per loop sage: %timeit bool(f) The slowest run took 34.95 times longer than the fastest. This could mean that an intermediate result is being cached. 1000000 loops, best of 3: 314 ns per loop
comment:49 Changed 19 months ago by
 Commit changed from ae4a2fc0526e9e9a88bf010fd9aceacdefed8b2b to 1505cceba6f52b43b00d92b09fa087a8eaff07fb
Branch pushed to git repo; I updated commit sha1. New commits:
1505cce  Keep Singular as optional backend for Hilbert series computations; some PEP8

comment:50 Changed 19 months ago by
I changed the (old!) error message to lower case and removed the fullstop. And I hope the style is now enough "spacey"...
Moreover, I added an optional parameter algorithm
.
However, one thing is very very strange now. A doctest that previously worked is now still giving the correct result, but also showing an odd comment; this does NOT happen when running the test interactively:
File "src/sage/rings/polynomial/multi_polynomial_ideal.py", line 2506, in sage.rings.polynomial.multi_polynomial_ideal.?.hilbert_series Failed example: I.hilbert_series() Expected: (t^4 + t^3 + t^2 + t + 1)/(t^2  2*t + 1) Got: skipping text from `parameter` (t^4 + t^3 + t^2 + t + 1)/(t^2  2*t + 1)
WTH??
New commits:
1505cce  Keep Singular as optional backend for Hilbert series computations; some PEP8

comment:51 Changed 19 months ago by
Apparently the comment is an artefact from a previous libsingular error:
sage: n = 4; m = 11; P = PolynomialRing(QQ, n * m, "x"); x = P.gens(); M = Matrix(n, x) sage: Minors = P.ideal(M.minors(2)) sage: hp = Minors.hilbert_polynomial(); hp 1/21772800*t^13 + 61/21772800*t^12 + 1661/21772800*t^11 + 26681/21772800*t^10 + 93841/7257600*t^9 + 685421/7257600*t^8 + 1524809/3110400*t^7 + 39780323/21772800*t^6 + 6638071/1360800*t^5 + 12509761/1360800*t^4 + 2689031/226800*t^3 + 1494509/151200*t^2 + 12001/2520*t + 1 sage: Minors.hilbert_polynomial(algorithm = 'singular')  RuntimeError Traceback (most recent call last) RuntimeError: Error raised calling singular function Exception RuntimeError: RuntimeError('Error raised calling singular function',) in 'sage.libs.singular.function.LibraryCallHandler.handle_call' ignored  RuntimeError Traceback (most recent call last) <ipythoninput2bfe31fe9a36b> in <module>() > 1 Minors.hilbert_polynomial(algorithm = 'singular') /home/king/Sage/git/sage/local/lib/python2.7/sitepackages/sage/rings/polynomial/multi_polynomial_ideal.pyc in __call__(self, *args, **kwds) 295 if not R.base_ring().is_field(): 296 raise ValueError("Coefficient ring must be a field for function '%s'."%(self.f.__name__)) > 297 return self.f(self._instance, *args, **kwds) 298 299 require_field = RequireField /home/king/Sage/git/sage/local/lib/python2.7/sitepackages/sage/rings/polynomial/multi_polynomial_ideal.pyc in hilbert_polynomial(self, algorithm) 2466 hilbPoly = sage.libs.singular.function_factory.ff.poly__lib.hilbPoly 2467 > 2468 hp = hilbPoly(self) 2469 t = ZZ['t'].gen() 2470 fp = ZZ(len(hp)1).factorial() /home/king/Sage/git/sage/local/lib/python2.7/sitepackages/sage/libs/singular/function.pyx in sage.libs.singular.function.SingularFunction.__call__ (build/cythonized/sage/libs/singular/function.cpp:15394)() 1328 if not (isinstance(ring, MPolynomialRing_libsingular) or isinstance(ring, NCPolynomialRing_plural)): 1329 raise TypeError("Cannot call Singular function '%s' with ring parameter of type '%s'"%(self._name,type(ring))) > 1330 return call_function(self, args, ring, interruptible, attributes) 1331 1332 def _instancedoc_(self): /home/king/Sage/git/sage/local/lib/python2.7/sitepackages/sage/libs/singular/function.pyx in sage.libs.singular.function.call_function (build/cythonized/sage/libs/singular/function.cpp:17499)() 1526 if errorreported: 1527 errorreported = 0 > 1528 raise RuntimeError("error in Singular function call %r:\n%s" % 1529 (self._name, "\n".join(error_messages))) 1530 RuntimeError: error in Singular function call 'hilbPoly': int overflow in hilb 1 error occurred in or before poly.lib::hilbPoly line 58: ` intvec v=hilb(I,2);` expected intvecexpression. type 'help intvec;' sage: P.<x,y,z> = PolynomialRing(QQ) sage: I = Ideal([x^3*y^2 + 3*x^2*y^2*z + y^3*z^2 + z^5]) sage: I.hilbert_series() skipping text from `parameter` (t^4 + t^3 + t^2 + t + 1)/(t^2  2*t + 1) sage: P.<x,y,z> = PolynomialRing(QQ) sage: I = Ideal([x^3*y^2 + 3*x^2*y^2*z + y^3*z^2 + z^5]) sage: I.hilbert_series() (t^4 + t^3 + t^2 + t + 1)/(t^2  2*t + 1)
What shall I do to fix it?
comment:52 Changed 19 months ago by
The easiest solution would be to remove the test that demonstrates the shortcomings of Singular (maybe that would also be better for my karma). Do you agree it can be removed?
comment:53 followup: ↓ 54 Changed 19 months ago by
Rather than remove it, I would mark it as either # known bug (20145)
or # not tested
(with a quick explanation why it is not tested).
comment:54 in reply to: ↑ 53 ; followup: ↓ 55 Changed 19 months ago by
Replying to tscrim:
Rather than remove it, I would mark it as either
# known bug (20145)
or# not tested
(with a quick explanation why it is not tested).
"not tested" might be an option. If I understand correctly, "known bug" would still mean that it is tested, right? That would be bad, because of the sideeffect of the error.
In any case, I will open a new ticket so that the sideeffect can be fixed. And I would prefer to remove the test for now and introduce it at the new ticket (showing that the sideeffect is fixed).
comment:55 in reply to: ↑ 54 Changed 19 months ago by
Replying to SimonKing:
Replying to tscrim:
Rather than remove it, I would mark it as either
# known bug (20145)
or# not tested
(with a quick explanation why it is not tested)."not tested" might be an option. If I understand correctly, "known bug" would still mean that it is tested, right? That would be bad, because of the sideeffect of the error.
I don't remember offhand. However, it is simple enough to test.
In any case, I will open a new ticket so that the sideeffect can be fixed. And I would prefer to remove the test for now and introduce it at the new ticket (showing that the sideeffect is fixed).
I would prefer # not tested
(unless that is what you meant by remove). At least that gives an example to the user where it doesn't work and the side effect is the only reason it would not be tested.
comment:56 Changed 19 months ago by
comment:57 Changed 19 months ago by
 Commit changed from 1505cceba6f52b43b00d92b09fa087a8eaff07fb to abe9c18f3441e29ed8d45258642dda7784092a94
Branch pushed to git repo; I updated commit sha1. New commits:
abe9c18  Remove a test because of a sideeffect

comment:58 Changed 19 months ago by
I hope that the last commit complies with your suggestion: The test is still present, but marked # not tested
, and it is commented that this is because of a sideeffect that is investigated at #26300.
comment:59 Changed 19 months ago by
 Branch changed from u/SimonKing/Hilbert_series_bug to u/tscrim/hilbert_series_bug20145
 Commit changed from abe9c18f3441e29ed8d45258642dda7784092a94 to ea5a851616b8c84fda8032d9f9c1219e50d03a99
 Milestone changed from sage7.1 to sage8.4
 Reviewers set to Travis Scrimshaw
Yes, that is good. I made a few other tweaks to match PEP8 and our coding/doc conventions in Sage. If my changes are good, then positive review.
New commits:
aac2eed  Merge branch 'u/SimonKing/Hilbert_series_bug' of git://trac.sagemath.org/sage into u/tscrim/hilbert_series_bug20145

ea5a851  A few reviewer tweaks.

comment:60 Changed 19 months ago by
 Status changed from needs_review to positive_review
I trust that the tests still pass (didn't check myself). Your changes look good to me. However, although I did read PEP8, I don't understand why you added blank spaces around some operators and removed blank spaces around others. And where can I read about doc formatting conventions of Sage (in particular when it is about optional parameters)?
comment:61 followup: ↓ 62 Changed 19 months ago by
PEP8 allows you to not write spaces for higher priority operators. So like in the exponentiation, that is higher than the addition, and thus, I don't need the spaces around the **
. Of course, that is my own stylistic preference showing. :) For this:
 def hilbert_series(self, grading = None, algorithm = 'sage'): + def hilbert_series(self, grading=None, algorithm='sage'):
 return self.hilbert_numerator(algorithm = 'singular') / (1  t)**n + return self.hilbert_numerator(algorithm='singular') / (1  t)**n
PEP8 says arguments should not have space around the =
.
Some of these things are not directly specified in the Sage doc conventions, but done as a general acceptedish convention. At least, an optional
to me means it does not have a default value (in some appropriate sense of the term for a user).
comment:62 in reply to: ↑ 61 ; followup: ↓ 64 Changed 19 months ago by
Replying to tscrim:
PEP8 allows you to not write spaces for higher priority operators. So like in the exponentiation, that is higher than the addition, and thus, I don't need the spaces around the
**
. Of course, that is my own stylistic preference showing. :)
My intention (although I quite often do typos) is to put spaces around...
 assignments of variables (but not of optional arguments):
bla = foo(bar)
 comparisons:
if bla <= foo:
+
and
 after (but not before) a comma or semicolon
and to not put spaces around
*
,/
and**
, because they have higher priority than+
and
:(2*x  1)**4)
 assignments of optional arguments:
foo(bla, bar=4)
Actually I dislike (4 * x  1) ** 3
, I find it difficult to read.
Would the above comply to PEP8?
comment:63 Changed 19 months ago by
 Commit changed from ea5a851616b8c84fda8032d9f9c1219e50d03a99 to a1545cc6c843fcdb85cbca3c85ca8ce49de28616
 Status changed from positive_review to needs_review
Branch pushed to git repo; I updated commit sha1 and set ticket back to needs_review. New commits:
a8d02dd  Fixing doctest output in palp_database.py due to signs.

7b2ed50  Merge branch 'u/tscrim/hilbert_functions26243' of git://trac.sagemath.org/sage into u/tscrim/hilbert_functions26243

5637e07  Fixing INPUT:: to INPUT: in hilbert.pyx.

a1545cc  Merge branch 'u/tscrim/hilbert_functions26243' of git://trac.sagemath.org/sage into u/tscrim/hilbert_series_bug20145

comment:64 in reply to: ↑ 62 Changed 19 months ago by
Replying to SimonKing:
Replying to tscrim:
PEP8 allows you to not write spaces for higher priority operators. So like in the exponentiation, that is higher than the addition, and thus, I don't need the spaces around the
**
. Of course, that is my own stylistic preference showing. :)My intention (although I quite often do typos) is to put spaces around...
 assignments of variables (but not of optional arguments):
bla = foo(bar)
 comparisons:
if bla <= foo:
+
and
 after (but not before) a comma or semicolon
and to not put spaces around
*
,/
and**
, because they have higher priority than+
and
:(2*x  1)**4)
 assignments of optional arguments:
foo(bla, bar=4)
Actually I dislike
(4 * x  1) ** 3
, I find it difficult to read.Would the above comply to PEP8?
By my understanding, yes. Of course, I don't always follow it as PEP8 are more guidelines, but I do try to stay as close to it as I can.
comment:65 Changed 19 months ago by
 Status changed from needs_review to positive_review
Last commit was a trivial fix of a doctest and merging in a fix of a bad block from #26243. I'm allowing myself to reset to positive.
comment:66 Changed 19 months ago by
 Branch changed from u/tscrim/hilbert_series_bug20145 to a1545cc6c843fcdb85cbca3c85ca8ce49de28616
 Resolution set to fixed
 Status changed from positive_review to closed
Where is the error?
44 is the number of generators of P in
n=4;m=11;P = PolynomialRing(QQ,n*m,"x");
, and Sage simply returns(1t)^P.ngens()
as the denominator. Is this wrong mathematically?