A common point among some Discourse topics on Lagrange interpolation is the non-recommendation due to the Runge’s phenomenon (Edit: or numerical stability).

However, that is a problem on real numbers, not on finite fields. Also, weighting or rearranging samples makes nonsense in basic use cases.

Lagrange interpolation is supported by several packages, but some assume real numbers.
As far as I know, SpecialPolynomials.jl is the most suitable package. Although I have no major complaints on SpecialPolynomials.jl, the implementation uses the barycentric form. (The form in itself is fine, but the Lagrange object requires allocations each time nodes are updated.) Also, the Newton interpolation there is somewhat lacking in scalability, since the intermediate data is kept in a Matrix.

Of course, I can customize the code for my own use, but I prefer to reuse existing code.
If you know of any other good packages, please let me know.

Thanks for the introduction to LaurentPolynomials.jl and FiniteFields.jl.

I don’t understand about the internals of LaurentPolynomials.jl, is it really a Lagrange interpolation? based on the Newton form?

Of course, the interpolated result can be obtained if the underlying polynomial is known, but for my purposes, the emphasis is on the process of calculating from the nodes rather than the interpolated result.

I am not sure what is your question. Yes, the function does Lagrange interpolation, computing
a polynomial which takes the values vals at the points pts.

Aside: Polynomial interpolation is not the problem that causes Runge’s phenomenon, it’s equally spaced points. If you use e.g. Chebyshev nodes, then you can interpolate to very high orders without issue. (Note that Runge phenomena have nothing with roundoff errors, so it’s not completely obvious to me whether there is an analogue for finite fields … I’m not sure how you analyze convergence rates of polynomial interpolation for that case.)

A separate issue is that the Lagrange interpolation formula is numerically unstable. Even with Chebyshev points, this is why people use e.g. the Barycentric interpolation formula instead (see e.g. Trefethen (2004)). Of course, this doesn’t matter in exact arithmetic, e.g. for finite fields.

My challenge in the background is to implement security protocols based on existing algorithms in IoT devices (with GF(2^8)). It would be implemented in C or Rust, or possibly assembly language.

I would like to simulate and verify the algorithms at higher level with Julia.
The reason I (i.e., the algorithm I am trying to implement) use Lagrange interpolation is that there is no need to know the polynomial underlying.

The LaurentPolynomials.jl and FiniteFields.jl are definitely one of the solutions.

However, the family is intentionally separated from the Polynomials.jl and GaloisFields.jl family.

There is also the Nemo.jl family, although with a slightly different target.

I would like to implement the algorithms in Julia with T <: Number and minimize direct dependencies on specific package families. However, compatibility is a practical problem.

(Personally, I find it beneficial that elements of GF(2^8) can be represented by 1 byte for low-level implementation.)