IMPL: generic and complex integration

[GComitini]

This draft PR implements integration of complex functions of one real variable. To do so, it first makes integration generic enough and then implements it for both real and complex functions. One of the guiding principles of this PR was not to change the pre-existing APIs. Of course, I also didn't change the current implementation of real integrals.

Theory

Since the integral of a complex function can be computed by splitting it into its real and imaginary part and separately computing two real integrals, all the integration methods currently defined in peroxide can be easily lifted to complex integration (with some qualifications needed for Gauss-Kronrod quadrature, which I'll discuss later). It would make little sense to re-define quadratures from scratch for complex number. The most logic way to proceed is making quadratures generic enough that they apply to both real and complex functions, and then implement the quadratures for both.

As to what type should be the subject of the generics, instead of defining a RealOrComplex trait or using Num from the num crate, I decided to go with an implementation which would make it possible to easily extend integration to types like vectors or matrices. Just like complex numbers, by linearity, the integral of vector/matrix functions of one real variable can be computed by evaluating a certain number of real integrals. By not committing to "numbers", but using more general traits instead, integration is easily implemented for these types too.

Implementation

How general? I'd say as general as the single quadratures needs them to be. For this reason I introduced three new traits: NCIntegrable, GLKIntegrable and GKIntegrable. These represent, respectively, types (which are the codomain of functions) which can be integrated using the Newton-Cotes quadrature, the Gauss-Legendre or the Kronrod quadrature, and the Gauss-Kronrod quadrature. Gauss-Legendre and Kronrod in theory should be kept separate, but computing these quadratures requires performing the very same operations (modulo tables), so in practice they coincide. Each of these traits is fairly easy to implement and only defines functions/constants which are actually needed to compute the quadratures. Keeping them separate leaves room, in the future, to implement integration for types which may not be integrable with some of the quadratures mentioned above, but that can still be integrated using the others. All traits are public so that users of the library can implement them for their own types.

Having defined the three traits, the building blocks of peroxide's integration routines, that is, fn newton_cotes_quadrature, fn gauss_kronrod_quadrature, fn unit_gauss_legendre_quadrature and fn unit_kronrod_quadrature are generalized to make use of those traits. All other functions are just wrappers around these, so they don't need more changes (asides from having to move some commutative operands from left to right to deal with the stupidity of rust being unable to overload operations on the left). In particular, switching to generic functions the way I did does not break the API and one will be able to use, for complex functions, the same library functions used for real ones.

Details (not all of them)

This PR introduces two associated types, NCIntegrable::{NodeY,PolynomialY}. These generalize node_y: Vec<f64> and p: Polynomial for Newton-Cotes integration. Implementing NCIntegrable requires defining how to compute a NodeY and a PolynomialY, and how to integrate the latter and evaluate it at a point. For types like complex numbers, vectors, matrices, etc., these operations should simply be duplicated over the codomain's degrees of freedom (e.g. real/imaginary part for complex numbers, single entries in vectors/matrices). Evaluation then rearranges the result in the correct format (e.g. the real and imaginary part are independently evaluated and then combined to obtain a single complex number; for vector/matrices the entries would be evaluated independently and then combined into a vector/matrix of the same format etc.)

Gauss-Kronrod quadrature works with a single tolerance, whereas the codomain of complex/vector/matrix integrals has 2 or more degrees of freedom. What to do with the tolerance depends on the implementation, and this is reflected in the existence of a is_within_tol trait function whose purpose is to say whether a certain value (real/complex number, vector, matrix, etc.) is within the specified tolerance. is_within_tol has a default implementation, which declares the value to be within tolerance if its gk_norm is. The gk_norm is the value used to compute relative tolerances and generalizes the absolute value of real numbers. For complex numbers it is taken to be equal to their modulus, and is_within_tol is left to its default. What this means is this PR computes the Gauss-Kronrod quadrature of complex functions by applying its tolerance to the modulus of G-K. One alternative could be using a different gk_norm, like max(|z.re|, |z.im|).

Notes

Some of the changes introduced by this MR are whitespace changes. This is because the crate is not formatted using cargo fmt. I'd ask you to please do so, as it makes contributing to the project much easier. I'd then rebase my changes on the formatted versions.

One example of using the functionality provided by this PR is here.

[Axect]

I sincerely apologize for the delay in reviewing your PR. I've been swamped with work recently, and I regret not getting to this sooner.

After carefully examining your detailed description and code changes, I must say I'm thoroughly impressed. Your approach to implementing complex function integration is both elegant and well-thought-out.

The introduction of NCIntegrable, GKIntegrable, and GLKIntegrable traits to generically implement operations for both f64 and C64 is a brilliant design choice. I particularly appreciate how you've managed to extend the functionality to include complex numbers without affecting the existing f64 operations. This demonstrates a deep understanding of the library's architecture and Rust's type system.

Your code design is clean and efficient, leaving little room for improvement. The way you've structured the changes maintains API compatibility while significantly expanding the library's capabilities.

However, I do have a couple of questions:

1. Regarding the NCIntegrable trait, while 'NodeY' seems like an appropriate name, 'PolynomialY' feels a bit ambiguous.
   Could you elaborate on your reasoning behind this name choice?
2. I noticed that you've included Clone but not Copy in the GKIntegrable trait. Was this a deliberate decision to allow for future extensions? I'd be interested in hearing your thoughts on this.

Once again, thank you for your outstanding contribution. Your work will significantly enhance Peroxide's functionality, and we're excited to merge it once we've addressed these minor points.

Oh, one last thing: Could you please move your work to the '35-support-complex-number' branch? I've set it up for the Complex feature. Thanks!

[GComitini]

Thanks for the appreciation @Axect!

1. Regarding the NCIntegrable trait, while 'NodeY' seems like an appropriate name, 'PolynomialY' feels a bit ambiguous. Could you elaborate on your reasoning behind this name choice?

I thought I bit about this but I couldn't come up with a better naming choice that was also concise enough. In the end I settled with PolynomialY in order not to reuse Polynomial and since its choice depends on the choice of the codomain of integration, Y. However, I am not committed to this name, and I'm open to other possibilities. For example, I just realized that we could even just call it Polynomial, since the only way to access that type is by typing NCIntegrable::Polynomial for a specific NCIntegrable all the more (e.g. f64::Polynomial or C64::Polynomial), so that the choice of the codomain is already made explicit by the namespace and there should be no overlap with the pre-existing Polynomial type. I just checked and the code builds with the substitution PolynomialY -> Polynomial.

2. I noticed that you've included Clone but not Copy in the GKIntegrable trait. Was this a deliberate decision to allow for future extensions? I'd be interested in hearing your thoughts on this.

Yes, it was in fact intentional. While for complex and real numbers (i.e. C64 and f64) Copy would have been appropriate, integration the way I implemented it in this PR is easily extended to vectors and matrices (even if I haven't tried to actually do so). These types won't be able to implement the Copy trait, whereas they will be able to implement Clone.

[GComitini]

Actually, on a second thought, we could rename NCIntegrable::Polynomial as NCIntegrable::NCPolynomial, since that type defines a polynomial structure specifically used for Newton-Cotes integration.

[Axect]

Thank you for your detailed explanations. They've addressed my questions thoroughly.

1. I agree with renaming to NCIntegrable::NCPolynomial. It's clear and specific to Newton-Cotes integration.

2. Your reasoning for using Clone instead of Copy in GKIntegrable makes perfect sense. It's a good choice for future extensibility.

Your suggestions resolve all my concerns. Please implement these naming changes, then we can proceed with merging your PR into the '35-support-complex-number' branch.

Great work on this contribution. Your design choices are well-considered and will significantly enhance the library.

[GComitini]

Done! How shall we proceed with the docs for the new traits and the small changes that may be needed to the old docs?

[Axect]

Thank you for implementing the changes. I've reviewed the modifications, and they look good.

I've also confirmed that all tests are passing successfully.

Regarding the documentation work you mentioned, let's handle that as a separate PR after we merge these changes. This approach will help us keep the current PR focused and easier to review.

I'll be changing the target of this PR from the master branch to the '35-support-complex-number' branch to align with our development workflow.

Once again, I appreciate your valuable contribution. I apologize for any delays in PR reviews lately - my workload has increased significantly. I'll do my best to review PRs as soon as I can find time, and I appreciate your understanding.

[GComitini]

Thanks for the merge @Axect!

Regarding the documentation work you mentioned, let's handle that as a separate PR after we merge these changes. This approach will help us keep the current PR focused and easier to review.

ACK.