Add support for more general boundary constraints

### Feature description

As it currently stands, the most general type of boundary constraints are of the form $\frac{f(X) - b(X)}{Z(X)}$ where

1. $f(X)$ is the trace column polynomial,
2. $b(X)$ is a polynomial specifying the boundary values enforced on the aforementioned column,
3. $Z(X)$ is a zerofier polynomial on the enforcement domain associated to the boundary constraint. Currently, this is either $X - g^0$ (i.e., `first`) or $X - g^{-2}$ (i.e., `last`) where we assume one exemption point and $g$ is the trace domain generator.

However, there are situations where we want to have boundary constraints which involve more than just one trace column, say $f_0(X)$ and $f_1(X)$ and a constraint like $\frac{f_0(X) \cdot f_1(X) - 1}{X - g^{-2}}$.

### Why is this feature needed?

This is needed in order to be able to write some constraints where the interaction between components happen through different buses and the consistency of the interactions is enforced through the above kind of boundary constraints.

The issue describes extending boundary constraints in an AIR compiler. Start by reading the code that handles boundary constraints, likely in a module like `constraints` or `boundary`. Look for the current implementation of single-column constraints and the polynomial division logic. Understand how the zerofier polynomials are constructed. Then, design a way to represent multi-column constraints and integrate them into the constraint system. Testing will involve creating new AIR scripts that use these generalized constraints and verifying the generated proofs.

Written by the indexing model from the issue text.