Weyl algebra division by zero error in characteristic $p>0$.

[nielsthl]

In reducing polynomials to standard form in the Weyl algebra, the commutation rule $d x = x d + 1$ is applied. This ultimately leads to computation of binomial coefficients in the coefficient ring. In positive characteristic this fails in the kernel:

A = ZZ/5[x, d, WeylAlgebra => { x => d}]
d^20*x^10

resulting in

stdio:2:4:(3): error: division by zero

Suggested fix (courtesy of GPT 5.2 and @antonleykin): the culprit seems to be buried in the binomial function in

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;
	}

Suggested fix using https://en.wikipedia.org/wiki/Lucas%27s_theorem:

ring_elem WeylAlgebra::binomial(int top, int bottom) const
{
  // Assumption: 0 <= bottom <= top.
  if (bottom == 0) return K_->from_long(1);
  if (bottom == 1) return K_->from_long(top);

  const int p = K_->characteristic();  // (or however the coeff ring exposes it)

  // In char p, avoid division by 0 when bottom >= p
  if (p > 0 && bottom >= p)
    {
      ring_elem result = K_->from_long(1);

      int n = top;
      int k = bottom;

      while (k > 0)
        {
          int ni = n % p;
          int ki = k % p;

          if (ki > ni)
            {
              K_->remove(result);
              return K_->from_long(0);
            }

          // Safe: 0 <= ki <= ni < p, so denominators 1..(p-1) never hit 0
          ring_elem term = binomial(ni, ki); // this call will *not* re-enter Lucas
          ring_elem tmp  = K_->mult(result, term);

          K_->remove(result);
          K_->remove(term);
          result = tmp;

          n /= p;
          k /= p;
        }
      return result;
    }

  // existing code path (works in char 0, and in char p when bottom < p)
  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 in char 0; safe in char p if bottom < p
      K_->remove(c);
    }
  return result;
}

[d-torrance]

Thanks for the report! Would you like to open a pull request with the fix?

[nielsthl]

Thanks for the report! Would you like to open a pull request with the fix?

I will do my best (first PR to M2).