search

sagemath / sage

Main repository of SageMath
https://www.sagemath.org
Other
1.44k stars 480 forks source link

Hilbert series bug #20145

Closed stumpc5 closed 6 years ago

stumpc5 commented 8 years ago

UPSTREAM: https://www.singular.uni-kl.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

Depends on #26243

Upstream: Reported upstream. Developers acknowledge bug.

Component: commutative algebra

Keywords: Hilbert series, polynomial ring

Author: Simon King

Branch/Commit: a1545cc

Reviewer: Travis Scrimshaw

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

dimpase commented 8 years ago
comment:1

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 (1-t)^P.ngens() as the denominator. Is this wrong mathematically?

stumpc5 commented 8 years ago
comment:2

The correct answer, as given in the second computation hilb(std(I)), is

( 120*t^3 + 135*t^2 + 30*t + 1 ) / (t-1)^14

which is not what Sage outputs. (But the total degree -11 is still correct...)

dimpase commented 8 years ago
comment:3

well, for m=10 you only get the output above due to the fact that numerator has a factor (t-1)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...

stumpc5 commented 8 years ago
comment:4

Replying to @dimpase:

well, for m=10 you only get the output above due to the fact that numerator has a factor (t-1)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 (1-t)^30, but - as you write - it is irreducible instead.

stumpc5 commented 8 years ago
comment:5

Replying to @dimpase:

Numerator is computed by Singular and so it should not give any bug, I guess...

At least luisfe on sage-devel (https://groups.google.com/d/msg/sage-devel/tyTDd5dEDGw/8FxH9hUzBwAJ) seems to confirm this to be a bug in Singular...

dimpase commented 8 years ago
comment:6

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 (t-1)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 (t-1)P.ngens().

-- 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 (1-t)^30, but - as you write - it is irreducible instead.

So it is a bug in Singular, indeed.

stumpc5 commented 8 years ago
comment:7

Replying to @dimpase:

Replying to @stumpc5:

Replying to @dimpase:

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 (t-1)P.ngens().

I actually did read the code, but

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).

dimpase commented 8 years ago
comment:8

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 t-1, 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...

dimpase commented 8 years ago

Upstream: Not yet reported upstream; Will do shortly.

stumpc5 commented 8 years ago
comment:9

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 t-1, 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 (t-1)?

dimpase commented 8 years ago
comment:10

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 t-1, 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 (t-1)?

It's all a bit Greek to me now, and https://www.singular.uni-kl.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.

stumpc5 commented 8 years ago
comment:11

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 Stanley-Reisner 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.

dimpase commented 8 years ago
comment:12

One would try to compute the same with Macaulay2. Did you try this?

stumpc5 commented 8 years ago
comment:13

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
dimpase commented 8 years ago
comment:14

by the way, here is what they mean by the 1st and by the 2nd Hilbert series. https://www.singular.uni-kl.de/DEMO_HTML/Examples/Hilbert/index.html

Indeed, this alone seems enough for a bug report. (save the Singular code :-))

stumpc5 commented 8 years ago
comment:15

Okay, I am sending a bug report to Singular now...

dimpase commented 8 years ago
comment:16

any progress on this?

stumpc5 commented 8 years ago
comment:17

You find the bug report from Feb at https://www.singular.uni-kl.de:8005/trac/ticket/750, but no one seems to have looked at it so far...

dimpase commented 8 years ago

Changed upstream from Not yet reported upstream; Will do shortly. to Reported upstream. No feedback yet.

stumpc5 commented 6 years ago
comment:19

As rechecked at http://www.singular.uni-kl.de:8002/trac/ticket/750#comment:1, this appears to be an interface bug, and no bug directly in singular.

fchapoton commented 6 years ago
comment:20

Looks like an overflow error

sage: QQ['z'](hilb(I, 2))*(1-z)**30-QQ['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
stumpc5 commented 6 years ago
comment:21

On the singular ticket they observe that what is shown here as the 2nd Hilbert series is actually their 1st Hilbert series.

fchapoton commented 6 years ago
comment:22

But the long singular answer should factorise with a factor (1-t)^30, no ? And it does not!

fchapoton commented 6 years ago
comment:23

So (I double checked) this is a bug in singular indeed. Someone should reopen ​http://www.singular.uni-kl.de:8002/trac/ticket/750

fchapoton commented 6 years ago
comment:24

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)

stumpc5 commented 6 years ago
comment:25

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.

fchapoton commented 6 years ago
comment:26

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.

stumpc5 commented 6 years ago
comment:27

oh, I see, sorry for my ignorance.

fchapoton commented 6 years ago

Changed upstream from Reported upstream. No feedback yet. to Reported upstream. Developers acknowledge bug.

fchapoton commented 6 years ago
comment:28

https://www.singular.uni-kl.de:8005/trac/ticket/750 has been re-opened

fchapoton commented 6 years ago

Description changed:

--- 
+++ 
@@ -1,3 +1,5 @@
+UPSTREAM: https://www.singular.uni-kl.de:8005/trac/ticket/750
+
 This works correctly:
simon-king-jena commented 6 years ago
comment:30

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.

simon-king-jena commented 6 years ago
comment:31

The code in #25243 has improved, it now beats Singular also speedwise.

So, I suggest to use #25243 here.

tscrim commented 6 years ago
comment:32

Replying to @simon-king-jena:

The code in #25243 has improved, it now beats Singular also speedwise.

So, I suggest to use #25243 here.

Typo: #26243 :P

simon-king-jena commented 6 years ago

Dependencies: #26243

johnperry-math commented 6 years ago
comment:34

Replying to @simon-king-jena:

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.

simon-king-jena commented 6 years ago
comment:35

If I see that correctly, the following functions need to be taken care of (all in sage.rings.polynomial.multi_polynomial_ideal):

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 Greuel-Pfister book.

johnperry-math commented 6 years ago
comment:36

Replying to @simon-king-jena:

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 called hilbert_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 called first_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 re-implemented.

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 non-polynomial for non-standard 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 4-5 years ago.

simon-king-jena commented 6 years ago
comment:37

Replying to @johnperry-math:

Replying to @simon-king-jena:

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 called hilbert_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 called first_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 re-implemented.

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.

simon-king-jena commented 6 years ago

Branch: u/SimonKing/Hilbert_series_bug

simon-king-jena commented 6 years ago

Commit: dcc3f17

simon-king-jena commented 6 years ago

Author: Simon King

simon-king-jena commented 6 years ago
comment:39

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:

7698bb6Some code optimizations in Hilbert series computation
8ae070bUse flint polynomial boilerplate for Hilbert series computation
d997abeFix a corner case in Hilbert series computation
5af3439Fix the erroneous fix of some corner cases
138f85bTurn some functions of polynomial.hilbert to methods of ETuple
e900222Hilbert series: Remove unnecessary bound checks, simplify code branches
f0e5c85Add documentation to new ETuple methods
2d29b9bUsing slices, custom median, and formatting/P#P8 tweaks to Hilbert code.
07ff2e1Reverting use of get_exp for __getitem__.
dcc3f17Do not use Singular for Hilbert series computations
dimpase commented 6 years ago
comment:40

Impressive speedup! Where does it come from? A better algorithm?

simon-king-jena commented 6 years ago
comment:41

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 "Hilbert-Poincaré 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 32-bit integer limitation. But in the end, it was also about speed: The new implementation beats Singular in medium size examples.

The algorithm for the Hilbert-Poincaré 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 Hilbert-Poincaré 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.

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

Changed commit from dcc3f17 to 6318242

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

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

6318242Fix a sign in Hilbert polynomial computation
simon-king-jena commented 6 years ago
comment:43

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 Hilbert-Poincaré 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
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

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

ae4a2fcAdd consistency check for Hilbert polynomial and fix a corner case
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

Changed commit from 6318242 to ae4a2fc