Fix binomial coefficient computation in Weyl algebras over ZZ/p by nielsthl · Pull Request #4087 · Macaulay2/M2

nielsthl · 2026-01-11T19:26:10Z

Fixes division-by-zero errors when computing binomial coefficients
in characteristic p for indices >= p by implementing Lucas's theorem.

This PR addresses issue #4082.

Changes to weylalg.cpp:

Add Lucas's theorem implementation for computing C(n,k) mod p
Recursively decompose n and k into base-p digits
Use the identity: C(n,k) ≡ ∏ C(n_i, k_i) (mod p)
Preserves original implementation for char 0 and small indices

Add comprehensive test suite in weyl-binomial.m2:

Tests over QQ with large exponents (up to 50)
Tests over ZZ/p for 25+ different primes
Verifies Lucas's theorem patterns (exponents up to 342)
Tests multivariate Weyl algebras
Includes Singular nc_algebra compatibility verification
All Part E tests use exponents > 20 to test beyond precomputed tables

Fixes division-by-zero errors when computing binomial coefficients in characteristic p for indices >= p by implementing Lucas's theorem. This PR addresses issue Macaulay2#4082. Changes to weylalg.cpp: - Add Lucas's theorem implementation for computing C(n,k) mod p - Recursively decompose n and k into base-p digits - Use the identity: C(n,k) ≡ ∏ C(n_i, k_i) (mod p) - Preserves original implementation for char 0 and small indices Add comprehensive test suite in weyl-binomial.m2: - Tests over QQ with large exponents (up to 50) - Tests over ZZ/p for 25+ different primes - Verifies Lucas's theorem patterns (exponents up to 342) - Tests multivariate Weyl algebras - Includes Singular nc_algebra compatibility verification - All Part E tests use exponents > 20 to test beyond precomputed tables Co-Authored-By: Claude Sonnet 4.5 <noreply@anthropic.com>

d-torrance · 2026-01-11T23:50:41Z

Thank you!

However, this is causing a test error from the Complexes package:

i1 : -- test source: ../../../../Macaulay2/packages/Complexes/FreeResolutionTests.m2:574:5-589:3
     
     -*
       restart
       needsPackage "Complexes"
     *-  
       -- Weyl algebras
       R = ZZ/101[s,t, Degrees => {0,0}]/(s^2-t^2-1)

o1 = R

o1 : QuotientRing

i2 :   S = R[x,y,Dx,Dy, WeylAlgebra => {{x,Dx}, {y,Dy}}, Degrees => {1, 1, -1, -1}, Join => false]

o2 = S

o2 : PolynomialRing, 2 differential variable(s)

i3 :   I = ideal(x*Dx + s*t, y*Dx)

o3 = ideal (x*Dx + s*t, y*Dx)

o3 : Ideal of S

i4 :   assert isHomogeneous I

i5 :   gbTrace=3

o5 = 3

i6 :   gens gb I

   -- registering gb 1 at 0x7efe7e8e3a80

   -- [gb]{2}(2)mm{3}(1)m{4}(4)mooo{5}(2)oonumber of (nonminimal) gb elements = 5
   -- number of monomials                = 9
   -- #reduction steps = 7
   -- #spairs done = 9
   -- ncalls = 3
   -- nloop = 6
   -- nsaved = 0
   -- 
o6 = | sty (t3+t)y yDx xDx+st |

             1      4
o6 : Matrix S  <-- S

i7 :   m = gens I

o7 = | xDx+st yDx |

             1      2
o7 : Matrix S  <-- S

i8 :   -- TODO: once git issue as #2789 is fixed this will give an error, so change it to just call the free resolution.
       assert try (F = freeResolution I; false) else true

   -- registering SchreyerOrder 0 at 0x7efe7ea4f450

   -- registering gb 2 at 0x7efe7e8e38c0

   -- [gb]{2}(2)mm{3}(1)m{4}(4)mKilled

nielsthl · 2026-01-12T09:51:41Z

The code seems to stall because of a dangling pointer somehow (coming from K_->mult and my lack of understanding of the M2 code). An interesting observation is that the stable version of the M2 engine fails on the same example when we remove try:

needsPackage "Complexes"
R = ZZ/101[s,t, Degrees => {0,0}]/(s^2-t^2-1)
S = R[x,y,Dx,Dy, WeylAlgebra => {{x,Dx}, {y,Dy}}, Degrees => {1, 1, -1, -1}, Join => false]
I = ideal(x*Dx + s*t, y*Dx)
gbTrace=3
gens gb I
m = gens I
F = freeResolution I

gives

-- [gb]z{5}(2)stdio:10:18:(3): error: division by zero

Seems to be related to #2789. Supposing that there is no dangling pointer or similar in the PR, perhaps one should take another look at #2789 - @mikestillman

nielsthl · 2026-01-13T09:26:52Z

Just to be on the same page: I am waiting for feedback on #2789 - the stable version gives a division by zero error there. My PR gives non-termination. The test associated with this is embedded in

assert try (F = freeResolution I; false) else true

and therefore the error in the stable version passes this test. Is non-termination expected in #2789? @mikestillman

mahrud · 2026-01-13T16:56:31Z

M2/Macaulay2/tests/normal/weyl-binomial.m2

+-- Test: dx^3 * x^3 = x^3*dx^3 + 9*x^2*dx^2 + 18*x*dx + 6
+assert(dx^3 * x^3 == x^3*dx^3 + 9*x^2*dx^2 + 18*x*dx + 6)


Writing the test as a comment and as an assertion is redundant.

Will look into cleaning up test file. Any comments on the issue in #2789 related to the PR?

mahrud · 2026-01-13T17:00:06Z

M2/Macaulay2/tests/normal/weyl-binomial.m2

+Rp3e1 = (ZZ/3)[x, dx, WeylAlgebra => {x => dx}];
+f3e1 = dx^25 * x^25;
+assert(f3e1 == x^25*dx^25 + x^24*dx^24);
+<< "  ✓ p=3: dx^25*x^25" << endl


These printed lines are also redundant.

mahrud · 2026-01-13T17:02:36Z

M2/Macaulay2/tests/normal/weyl-binomial.m2

+--------------------------------------------------------------------------------
+-- Summary
+--------------------------------------------------------------------------------
+
+<< endl;
+<< "========================================" << endl;
+<< "All Weyl algebra binomial tests passed!" << endl;
+<< "========================================" << endl;
+<< endl
+<< "Tested:" << endl
+<< "  - QQ coefficients with exponents up to 50" << endl
+<< "  - ZZ/p for p = 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101" << endl
+<< "  - ALL Part E tests use exponents > 20 (testing outside precomputed table bounds)" << endl
+<< "  - Largest tested exponents: 342 (p³-1 pattern), 124 (p³-1 pattern), 120 (p²-1 pattern)" << endl
+<< "  - Lucas's theorem verification at very large exponents (up to 342)" << endl
+<< "  - Coefficient extraction and verification across fields" << endl
+<< "  - Multivariate Weyl algebras with exponents > 20" << endl
+<< "  - Singular nc_algebra compatibility tests (all exponents > 20):" << endl
+<< "    * Coefficient patterns for n = p+1, p+2, p+3" << endl
+<< "    * Multiples of p (n = kp) verification" << endl
+<< "    * Asymmetric exponent cases (up to 60)" << endl
+<< "    * Composite exponents with precise coefficient verification" << endl
+<< "    * Multivariate patterns matching Singular's nc_algebra" << endl
+<< "    * Special cases m = p^2 - 1 and m = p^3 - 1" << endl
+<< "    * Pattern periodicity verification across power levels" << endl
+<< "    * Very large symmetric exponents (up to 75)" << endl
+
+exit 0


This section is also redundant.

mahrud · 2026-01-13T17:07:50Z

M2/Macaulay2/tests/normal/weyl-binomial.m2

Having more unit tests is always nice but is a 1091 line unit test file really necessary?

This is a courtesy of Claude Code. Will look into cleaning up.

mahrud · 2026-01-13T17:08:44Z

M2/Macaulay2/e/weylalg.cpp

    }
  return result;
 }
+*/


Is there a reason you leave the old code as comment and duplicate it just a few lines above rather than deleting it?

mahrud · 2026-01-13T17:17:00Z

M2/Macaulay2/e/weylalg.cpp

+          ring_elem tmp  = K_->mult(result, term);
+
+          K_->remove(result);
+          K_->remove(term);
+          result = tmp;


What's the point of the extra variable tmp here? Would result = K_->mult(result, term); not suffice?

This should align with memory allocation/deallocation in the old code:

M2/M2/Macaulay2/e/weylalg.cpp

Lines 227 to 246 in 54c4a83

ring_elem WeylAlgebra::binomial(int top, int bottom) const

{

// This should be located elsewhere

// Assumption: bottom <= top, top >= 0, bottom >= 0.

if (bottom == 0) return K_->from_long(1);

if (bottom == 1) return K_->from_long(top);

if (top <= binomtop) return K_->from_long(binomtable[top][bottom]);

ring_elem result = K_->from_long(1);

for (int a = 0; a < bottom; a++)

{

ring_elem b = K_->from_long(top - a);

ring_elem result1 = K_->mult(result, b);

K_->remove(result);

K_->remove(b);

ring_elem c = K_->from_long(a + 1);

result = K_->divide(result1, c); // exact

K_->remove(c);

}

return result;

}

It seems that K_->mult returns a "pointer" to a new allocated object.

nielsthl · 2026-02-09T09:24:41Z

I am still waiting for an answer to #4087 (comment) @d-torrance @mahrud @mikestillman

Without this input the PR is stalled.

mikestillman · 2026-02-17T16:57:31Z

I'm looking into the noticed bug(s)...

nielsthl mentioned this pull request Jan 12, 2026

bug in computing syzygy module for tower rings involving the Weyl algebra in the presence of Schreyer orders #2789

Open

mahrud reviewed Jan 13, 2026

View reviewed changes

mahrud mentioned this pull request Jan 13, 2026

AI/LLM contribution policy #4066

Open

mahrud reviewed Jan 13, 2026

View reviewed changes

d-torrance requested a review from mikestillman February 10, 2026 17:08

		-- Test: dx^3 * x^3 = x^3dx^3 + 9x^2dx^2 + 18x*dx + 6
		assert(dx^3 * x^3 == x^3dx^3 + 9x^2dx^2 + 18x*dx + 6)

	ring_elem WeylAlgebra::binomial(int top, int bottom) const
	{
	// This should be located elsewhere
	// Assumption: bottom <= top, top >= 0, bottom >= 0.
	if (bottom == 0) return K_->from_long(1);
	if (bottom == 1) return K_->from_long(top);
	if (top <= binomtop) return K_->from_long(binomtable[top][bottom]);
	ring_elem result = K_->from_long(1);
	for (int a = 0; a < bottom; a++)
	{
	ring_elem b = K_->from_long(top - a);
	ring_elem result1 = K_->mult(result, b);
	K_->remove(result);
	K_->remove(b);
	ring_elem c = K_->from_long(a + 1);
	result = K_->divide(result1, c); // exact
	K_->remove(c);
	}
	return result;
	}

Conversation

nielsthl commented Jan 11, 2026

Uh oh!

d-torrance commented Jan 11, 2026

Uh oh!

nielsthl commented Jan 12, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

nielsthl commented Jan 13, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

mahrud Jan 13, 2026

Choose a reason for hiding this comment

Uh oh!

nielsthl Jan 13, 2026

Choose a reason for hiding this comment

Uh oh!

mahrud Jan 13, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

mahrud Jan 13, 2026

Choose a reason for hiding this comment

Uh oh!

mahrud Jan 13, 2026

Choose a reason for hiding this comment

Uh oh!

nielsthl Jan 13, 2026

Choose a reason for hiding this comment

Uh oh!

mahrud Jan 13, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

nielsthl Jan 13, 2026

Choose a reason for hiding this comment

Uh oh!

mahrud Jan 13, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

nielsthl Jan 13, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

nielsthl commented Feb 9, 2026

Uh oh!

mikestillman commented Feb 17, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

nielsthl commented Jan 12, 2026 •

edited

Loading

nielsthl commented Jan 13, 2026 •

edited

Loading

mahrud Jan 13, 2026 •

edited

Loading

mahrud Jan 13, 2026 •

edited

Loading

mahrud Jan 13, 2026 •

edited

Loading

nielsthl Jan 13, 2026 •

edited

Loading