Macaulay2 / M2

The primary source code repository for Macaulay2, a system for computing in commutative algebra, algebraic geometry and related fields.
https://macaulay2.com
345 stars 230 forks source link

(trim, Ideal) should work correctly in single variable polynomial rings #1273

Open kschwede opened 4 years ago

kschwede commented 4 years ago

When working over a PID, trim should work correctly, even if the ideal is not homogeneous.

i1 : R = ZZ/7[t];

i2 : I = ideal(t^3*(t+1), t*(t-1))

             4    3   2
o2 = ideal (t  + t , t  - t)

o2 : Ideal of R

i3 : trim I

                2
o3 = ideal (t, t  - t)

o3 : Ideal of R

Obviousy the gcd command on the generators works ok, but presumably this should work via trim as well.

This probably isn't a priority for version #1.16, but it should be easy to implement.

pzinn commented 4 years ago

I personally use this on my inhomogeneous ideals. it will in particular work for a single variable.

-- simplifying inhomogeneous ideals
simplify = x -> (
    R:=ring x;
    local u;
    flag:=true;
    y:=new MutableList from first entries gens x;
    while flag do (
    flag=false; 
    scan(#y, i-> (
        u=y#i;
        y#i=0_R;
        y#i=u%(ideal toList y);
        if y#i =!= u then flag=true;
        ));
    );
    trim ideal toList y -- just to remove zeroes
    );
kschwede commented 4 years ago

That's very nice and I am definitely taking it :-). Although maybe then there should be a simplify function in Macaulay2 if we don't want this behavior implemented into trim (I imagine for Euclidean domains though, the gcd function might be faster if it's implemented at a lower level, athough I didn't look at the code).

pzinn commented 2 years ago

I realize this still hasn't been fixed. it's a bit ridiculous that M2 still can't simplify this:

i1 : R=QQ[x]

o1 = R

o1 : PolynomialRing

i2 : ideal(3*x^4-x^3+x^2-x-2,3*x^4+2*x^3+6*x+4)

              4    3    2            4     3
o2 = ideal (3x  - x  + x  - x - 2, 3x  + 2x  + 6x + 4)

o2 : Ideal of R

i3 : trim oo

              3    2             4    3    2
o3 = ideal (3x  - x  + 7x + 6, 3x  - x  + x  - x - 2)

o3 : Ideal of R
mikestillman commented 2 years ago

I agree!

In general, how should we trim ideals which are not homogeneous?

pzinn commented 11 months ago

This is something that regularly bothers me. I decided to take a look at modules2.m2 and see how trim is implemented. @mahrud can you explain the logic behind trimHelper? looks like trim uses hooks, but actually it doesn't, since everything goes through trimHelper anyway...

mahrud commented 11 months ago

Not everything goes through trimHelper: https://github.com/Macaulay2/M2/blob/86d6c70f8746e6b13d95c53ad20d880e0d60a8dc/M2/Macaulay2/packages/LocalRings.m2#L354-L358 Hooks were added to accommodate a strategy for trimming modules over local rings, so the existing implementation in Core should have been added as the Default strategy. Except that the existing implementation already had a Strategy option with two possible values, Complement and null, so they needed to be separated into two hooks, one accessible with Strategy => Complement and one with a new name, Strategy => Inhomogeneous.

I could have split the code into two, but wanted to change as few lines as possible, so renamed the previous implementation trimHelper and added two hooks that call trimHelper with specific strategy values: https://github.com/Macaulay2/M2/blob/4be4f70a827facd35d4b2eb09f3069f7a0e2702e/M2/Macaulay2/m2/matrix2.m2#L247-L248

A similar thing happened with most other functions that needed new strategies for local rings.

I think you could easily add a new hook (after Complement and Inhomogeneous, so that it is attempted first) which does what you suggest when the ring is a PID and otherwise returns null, so the existing strategies are attempted, so something like this:

diff --git a/M2/Macaulay2/m2/matrix2.m2 b/M2/Macaulay2/m2/matrix2.m2
index 989c38855e..830bdd6a85 100644
--- a/M2/Macaulay2/m2/matrix2.m2
+++ b/M2/Macaulay2/m2/matrix2.m2
@@ -246,6 +246,18 @@ trimHelper = ((opts, M) -> (

 addHook((trim, Module), Strategy => Inhomogeneous, (opts, M) -> trimHelper(opts ++ {Strategy => Inhomogeneous}, M))
 addHook((trim, Module), Strategy => Complement,    (opts, M) -> trimHelper(opts ++ {Strategy => Complement},    M))
+addHook((trim, Module), Strategy => "PID",
+    (opts, M) -> (
+       R := ring M;
+       -- check that R is a principal ideal domain:
+       if not instance(R, PolynomialRing) or dim R > 1
+       -- TODO: add other cases
+       then return null;
+       -- modules are direct sums of cyclic modules.
+       -- ideals are generated by the gcd of their generators.
+       -- TODO: implementation here
+       ))
pzinn commented 11 months ago

sure I can do that. I still think the trimHelper makes the code hard to parse. I'm still not sure I understand these two lines:

addHook((trim, Module), Strategy => Inhomogeneous, (opts, M) -> trimHelper(opts ++ {Strategy => Inhomogeneous}, M))
addHook((trim, Module), Strategy => Complement,    (opts, M) -> trimHelper(opts ++ {Strategy => Complement},    M))
pzinn commented 11 months ago

actually, it seems strategy 3 of minimalPresentation (modules2.m2 L207) already implements pretty much what you suggest.

pzinn commented 11 months ago

one thing I can never figure up working with modules: When one defines a module M as an image, how can one easily access the embedding map from M to ambient M? the best I could come up with is gens M * map(cover M,M,1) which is kind of ugly since the second map is not well-defined...

pzinn commented 11 months ago

maybe this could just replace the current "Inhomogeneous" strategy? which doesn't do much anyway. edit: oh, is that strategy ever called anyway? help Inhomogeneous suggests that for gb computations this option is used as default, but for trim, AFAICT, it's never used unless explicitly introduced with Strategy=>Inhomogeneous...

mahrud commented 11 months ago

I agree that trimHelper is very convoluted and hard to read -- I didn't write it, only gave it a name, and would like to see it simplified.

The syntax addHook(key, Strategy => label, func) simply means hang a new strategy named label that runs func and is accessible via key. In this case:

key   = (trim, Module) -- meaning this hook is used by trim(Module)
label = Complement     -- just a way to distinguish this strategy
func  = (opts, M) -> trimHelper(opts ++ {Strategy => Complement}, M)

In particular, you can use the key and label to find func using either:

code hooks((trim, Module), Strategy => Complement)
code hooks( trim, Module,  Strategy => Complement)

(side note: I just remembered that locate and lookup don't work for hooks, would be nice to add) and if you wish you can run func on a module M like this:

runHooks((trim, Module), Strategy => Complement, (options trim, M))

actually, it seems strategy 3 of minimalPresentation (modules2.m2 L207) already implements pretty much what you suggest.

It must be wrong:

R = QQ[x]
I = ideal(3*x^4-x^3+x^2-x-2,3*x^4+2*x^3+6*x+4)
M = image gens I
N = prune M
rank M -- gives 1
rank N -- gives 2!!
N.cache.pruningMap -- gives 0!!!

When one defines a module M as an image, how can one easily access the embedding map from M to ambient M?

It's hard to remember, but this should do the trick: inducedMap(ambient M, M)

pzinn commented 11 months ago

I'll look into the prune bug you found

pzinn commented 11 months ago

I suspect smithNormalForm. in particular, the fact that it gives different results depending on the option ChangeMatrix is highly suspicious, but the documentation is not very clear on that.

pzinn commented 11 months ago

There seem to be at least two distinct issues:

i1 : R=QQ[x]

o1 = R

o1 : PolynomialRing

i2 : f=random(R^{1,2,3,4},R^{0,1})

o2 = {-1} | 7x  1   |
     {-2} | 8x2 7x  |
     {-3} | 6x3 4x2 |
     {-4} | 9x4 3x3 |

             4       2
o2 : Matrix R  <--- R

i3 : smithNormalForm f

o3 = (| 41x2 0 |, | -7x   1     0  0  |, {0}  | -1 0 |)
      | 0    1 |  | 1     0     0  0  |  {-1} | 7x 1 |
      | 0    0 |  | -10x2 -22x  41 0  |
      | 0    0 |  | -39x3 -12x2 0  41 |

o3 : Sequence

i4 : smithNormalForm (f,ChangeMatrix=>{false,true})

o4 = (| x2 0 |, {0}  | -1 0 |)
      | 0  1 |  {-1} | 7x 1 |

o4 : Sequence

the first one is that the result of gb depends (mildly) on the value of ChangeMatrix. annoying. that explains the constants being there or not on the diagonal. the other one (more serious) is the wrong size of the matrix in the second calculation.

pzinn commented 11 months ago

Fatal bug:

i1 : R=QQ[x]

o1 = R

o1 : PolynomialRing

i2 : smithNormalForm matrix{{x^3+2},{x^3+x^2}}

o2 = (| 2 |, | x2+1  -x2+x-2     |, {3} | 1 |)
      | 0 |  | x3+x2 -x3-2       |
      | 0 |  | x4-x2 -x4+x3-2x+2 |

the result doesn't even have the right size

pzinn commented 11 months ago

I don't know how to fix the smithNormalForm bug, I made a separate issue. The best I could come up with is a one-line workaround so that prune works normally: (modules2.m2 L220)

rows := first \ piv | toList(rank target f..<rank target g);
pzinn commented 11 months ago

that being said, this whole prune hook should probably be reviewed once smithNormalForm is fixed, since it's clearly out of sync with the code of the latter -- for example it doesn't trust the first matrix to be diagonal, even though it definitely should be.

pzinn commented 11 months ago

can finally be closed!