Implementing Number-theoretic functions #1114
Description
Fortran currently is well-supported in real-analysis problems, but the number-theory field is badly neglected (see #365 ). It would be nice if we had a more generally useful mathematical computing language that way users don't have to switch to another language just to do number-theory. (I personally do this myself, most of my Real arithmetic calculations is done in Fortran, but I always have to either switch to my Rust libraries for even basic number-theory or write Fortran bindings to them).
I propose that we add some of the following functionality. I can implement this during May-August 2026, but I want feedback and
interest first before I implement it given how much work it would be to efficiently work around the unfortunate restrictions of the Fortran standard (no unsigned integers, or 128-bit arithmetic).
This is just a preliminary proposal, some of the interfaces might be unintuitive for users, and I'm flexible with changing it.
Although Fortranners skew more mathematical than most programming communities, so the proposal might be close to what many users
are familiar with.
The formatting goes: (proposed_name) ({alternate_name}) (function_signature) (input_output_domains). Here Z and R refer to generic
signed integer and floating-point types respectively. N refers to a positive only Z. Since all additive inverses are mapped to their positive equivalent, all ring arithmetic modulo n is treated so that -n = n
Minimum Proposed functions
add_residue {mod_add}(in x, in y, in n) -> N
sub_residue {mod_sub}(in x, in y, in n) -> N
add_inverse (in a, in n) -> N
product_residue {mod_mul}(in x, in y, in n) -> N
quadratic_residue {mod_sqr}(in x, in n) -> N
exp_residue {mod_pow}(in a,in b,in n, out res) -> Error a^b mod n where a,b,n are in Z res \in N
optional legendre_symbol(in a, in n, out res ) -> Error
optional jacobi_symbol(in a, in n, out res) -> Error
Kronecker symbol is equivalent to all the above and is nearly as fast as the Jacobi symbol
Kronecker_symbol(in a, in n) -> {-1,0,1}
Gcd(a,b) -> N a,b \in Z
gcdv(A) -> N a is an array of Z elements
lcm(in a, in b, out prod) -> Error a,b \in Z
lcmv(in a, out prod) -> Error
factor(n) -> FactorizationType (prime_1,pow_1,...,etc)
This will likely use the Pollard-Brent algorithm
unit_cardinality {euler totient}(in n)-> N
least_exponent {carmichael_totient}(in n)-> N
extended_gcd(in a, in b,out x, out y) a,b,x,y \in Z
ord {multiplicative order}(a,n, out order) -> Error a,n \in Z, order \in N
prime_pi_bounds(in n, out lo,out hi) n \in Z, lo,hi \in N
n,lo,hi are in Real64 and or Real128
The reason why we implement this for reals is that estimates are mostly useful for very large intervals, no reason to restrict ourself to n < 2^63
And we can "secretly" use Axler's bounds which are far superior to the well-known n/ln(n) estimate that most analysis uses
prime_pi_bounds_real(in n,out lo, out hi) n,lo,hi \in R
mul_inverse(in a, in n, out inv) -> Error a,n, inv \in Z inv \in N
is_prime(n) n \in Z where a prime is only the product of a unit of Z {1,-1} and itself
This will use the BPSW test
I propose a non-standard representation of an integer factorisation
-1 stores the negative case,1 the positive case. If we factor an element of {-1,0,1} then we simply store it with no prime factors
FactorisationType
"unit" {-1,0,1}
primes
powers
Functions
These are mostly to exploit faster computation given an already factored integer
unit_cardinality(in x)
least_exponent
gcd(in x,in y) x,y FactorisationType
gcd_z(in x, in y) x is FactorisationType, y \in Z
lcm(in x, in y) x,y FactorisationType
lcm_z(in x, in y) x is a FactorisationType, y \in Z
Implement exponentiation since we can store much larger integers essentially free of charge
pow(in x,in p) x is FactorisationType, p is a Natural
product(in x, in y) x,y FactorisationType
Additional functions
This are additional functions, types that are less known/used by programmers but extensively used
by mathematicians. Residue Classes/Quotient rings in particular are the dominant structure used in actual research.
ResidueClass type
Allocated array of residue classes plus the ring
Functions
initialisers
from_units(in n, out RClass) -> Error (the number of units of a ring N can easily exceed memory, so we can restrict it to something like 1 trillion elements)
Calculate a ResidueClass whose elements are in both ResidueClasses
coprime_unify(in x, in y) -> ResidueClass
Coprime units are guaranteed to have solutions so we can avoid allocation and update X inplace
coprime_unify_inplace(in out x,in y)
Noncoprime CRT is just as useful, but does not have guaranteed solutions so we cannot update x inplace
noncoprime_unify(in x, in y) -> ResidueClass
Residue-Class Roots (modular sqrts)
Power Residue Symbol
Gaussian Integers
Generalised Quadratic Integers and possibly other integers