Can't take 'value' of a Hilbert series

[kroner]

Hilbert series are returned as Expressions, but I would like to work with them as rational functions. The method 'value' is supposed to do this conversion. However, the numerator and denominator are returned in a Laurent polynomial ring, and such rings give an error when you try to form their fraction field.

R = QQ[x,y]
H = hilbertSeries ideal{x}
value H

To get around this requires jumping through some hoops, subbing the numerator and denominator into a polynomial ring before applying value.

S = ZZ[gens degreesRing R]
H = sub(numerator H,S)/sub(denominator H,S)
value H

[DanGrayson]

For some reason, the engine doesn't support such fraction rings. It should be possible to modify the routine for making fraction fields to use the workaround you present.

[mahrud]

Is there a fix for this yet? I've heard from others that this is a real pain to deal with. Somehow there is no way to have the numerator and denominator of the Hilbert series in a ring without inverses:

i12 : newRing(degreesRing 1, Inverses => false, Global => true, Degrees => {1})
stdio:17:1:(3): error: not all variables are > 1, and Global => true

[pzinn]

I use the workaround below. Alternatively, one can use my package Factor (still not PRed due to annoying incompatibilities with some code that uses factor) which completely redefines everything in a more sensible manner. I'll cook up a little demo on the web server soon.

debug Core

frac EngineRing := R -> if isField R then R else if R.?frac then R.frac else (
     o := options R;
     if o.WeylAlgebra =!= {} or R.?SkewCommutative
     then error "fraction field of non-commutative ring requested";
     if not factoryAlmostGood R then error "not implemented yet: fraction fields of polynomial rings over rings other than ZZ, QQ, or a finite field";
     local F;
     if o.Inverses then (
	 R1:=newRing(R,Inverses=>false,MonomialOrder=>GRevLex);
	 f:=map(R1,R); g:=map(R,R1);
	 R.frac = F = frac R1; -- !!!
	 F.baseRings=append(F.baseRings,R);
	 promote(R,F) := (x,F) -> (f numerator x)/(f denominator x);
	 oldnum := F#numerator; oldden := F#denominator;
	 numerator F := (x) -> g oldnum x;
	 denominator F := (x) -> g oldden x;
	 lift(F,R) := opts -> (f,R) -> if isUnit denominator f then numerator f*(denominator f)^(-1) else error "cannot lift given ring element";
	 fraction(R,R) := (x,y) -> (f (numerator x*denominator y))/(f (numerator y*denominator x));
	 return F;
	 );
     R.frac = F = new FractionField from rawFractionRing R.RawRing;
     F.frac = F;
     F.baseRings = append(R.baseRings,R);
     commonEngineRingInitializations F;
     factor F := options -> f -> factor numerator f / factor denominator f; -- options?
     toString F := x -> toString expression x;
     net F := x -> net expression x;
     baseName F := (f) -> (
	  if denominator f != 1 
	  then error "expected a generator"
	  else baseName numerator f);
     expression F := (f) -> expression numerator f / expression denominator f;
     numerator F := (f) -> new R from rawNumerator raw f;
     denominator F := (f) -> new R from rawDenominator raw f;
     fraction(F,F) := F / F := (x,y) -> if y != 0 then x//y else error "division by 0";
     fraction(R,R) := (r,s) -> new F from rawFraction(F.RawRing,raw r,raw s);
     F % F := (x,y) -> if y == 0 then x else 0_F;	    -- not implemented in the engine, for some reason
     F.generators = apply(generators R, m -> promote(m,F));
     if R.?generatorSymbols then F.generatorSymbols = R.generatorSymbols;
     if R.?generators then F.generators = apply(R.generators, r -> promote(r,F));
     if R.?generatorExpressions then F.generatorExpressions = (
	  R.generatorExpressions
	  -- apply(R.generatorExpressions,F.generators,(e,x)->new Holder2 from {e#0,x})
	  );
     if R.?indexSymbols then F.indexSymbols = applyValues(R.indexSymbols, r -> promote(r,F));
     if R.?indexStrings then F.indexStrings = applyValues(R.indexStrings, r -> promote(r,F));
     if R.?numallvars then F.numallvars=R.numallvars;
     F)

[mahrud]

I'm glad you responded, because I kind of wish there was a simpler solution purely involving expressions. For instance, is there a simple way to add two Hilbert series?

R = kk[x,y,z]
K = koszul vars R
a = hilbertSeries ker K.dd_1
b = hilbertSeries ker K.dd_2

a + b

Current output:

i5 : a+b

       2    3       3
     3T  - T       T
o5 = -------- + --------
            3          3
     (1 - T)    (1 - T)

o5 : Expression of class Sum

Output that I'd like:

           2
         3T
o12 = --------
             3
      (1 - T)

o12 : Expression of class Divide

I realize general simplifications are most likely best handled as rational functions, but operations involving Hilbert series of a fixed ring always involve the same denominator, so maybe this case should be handled separately.

[pzinn]

I don't think expressions is the solution since by definition you can't do mathematical operations on them. here's the demo, BTW: https://www.unimelb-macaulay2.cloud.edu.au/?user=xxx#chat the first version ("old way") is just the workaround above which I've randomly inserted inside Factor for no good reason. the "new way" involves using Factor's own redefinitions.

[pzinn]

BTW Factor hasn't been updated in a while, in particular not since factor.m2 (sigh, lowercase vs uppercase) was introduced. there's also an issue with monomial ordering that needs fixing.

[mahrud]

I don't think expressions is the solution since by definition you can't do mathematical operations on them.

What do you think about this?

Divide + Divide := (x,y) -> if denominator x === denominator y then new Divide from {numerator x + numerator y, denominator x} else new Sum from {x,y}

[d-torrance]

Would it make sense to have a HilbertSeries class that maybe uses expressions for its net method function but otherwise allows for more Hilbert series-specific functionality?

[pzinn]

What do you think about this?

Divide + Divide := (x,y) -> if denominator x === denominator y then new Divide from {numerator x + numerator y, denominator x} else new Sum from {x,y}

the reason I'm reluctant to add such rules for expressions is that expressions are used in the output process, and they're really meant to be output as they are, not transformed any further. yes, there are a couple of exceptions already in place (like adding 0 doesn't do anything to an expression), and in fact one day we may regret these exceptions (maybe one day there'll be need to display "0+0"). But imagine the next stage: how about defining the product of two expressions being their actual product if they happen to be ring elements... and now factor does nothing at all...

I like better my "Factor" solution since Hilbert series really are normal rational functions that should be usable as such.

[mahrud]

I don't understand the issue. You can define 0+0:

i1 : new Sum from {0,0}

o1 = 0 + 0

You can also define new Sum from {x, y} for two hilbert series if you really want to keep them separate. The question is what's a more useful output. Another example is Power: I really think this should be simplified:

i3 : (new Power from {2, 3}) * (new Power from {2, 5})

      3 5
o3 = 2 2

[pzinn]

what I meant was, expression 1 + expression 2 will give you 1+2, whereas expression 0 + expression 2 will give you 2. BTW I think we're digressing from the original issue...

[mahrud]

This feels incomplete without a method for simplifying expressions. Has this always been the idea behind expressions? If so, I think I agree with Doug that maybe Hilbert series shouldn't be expressions, but another type that allows manipulations and taking values.

[pzinn]

yes, another type, like an actual ring element

[mikestillman]

maybe there should be a type: FactoredRingElement, or FactoredFraction. Multiplication is easy, inverse too, addition would take out the common parts and add the rest, maybe there should be a simplify function which would factor the rest too.

This could be a useful type, including with fractional ideals, where we often have denominators which only have a few different irreducible factors, but potentially (reasonably) high powers of these.

[pzinn]

Actually this is more or less exactly what my package Factor does, cf #1645 (at the time I gave up on the PR because Macs are unable to distinguish factor.m2 from Factor.m2 -- was that ever resolved?) In particular it takes care of Hilbert series in a much more satisfactory way. You can test it on Macaulay2Web, see screenshot

[pzinn]

That being said, there's a subtle difference between what Factor does and what you @mikestillman suggest: Factor introduces a type FactorPolynomialRing which is a descendent of PolynomialRing. You suggest a type FactoredRingElement which is a descendant of RingElement. I'd have to think hard what is better. Relatedly, it used to be that factor.m2 would define factor at the level of each polynomial ring, whereas now it acts on RingElement. I could imagine rewriting Factor in a similar manner.

[mahrud]

at the time I gave up on the PR because Macs are unable to distinguish factor.m2 from Factor.m2 -- was that ever resolved?

Certainly Core scripts don't get confused with files elsewhere on the path anymore.