See the last line, x, c1 || c2. We're joining two carry bits there because it seems that they are never 1 at the same time, but that remains to be proven. Exhaustive tests seem to pass, though.

In Julia, MyType{T, V} <: MyType{T}, so in some sense the first type argument is more "basic" than the second one. MLUInt/MLInt have arguments {N, T}, so the number of limbs is first, then the type. This is made by analogy with NTuple, for which this may be more logical, but for multi-limb integers the order may be better reversed.

Although, to be fair, it is exclusively a cosmetic change, since it is no problem to generalize the first type argument explicitly when necessary: MyType{Int, Bool} <: MyType{T, Bool} where T.

Useful for HEAAN. For reference see:

• HEAAN project (rns.jl)
• Richard Crandall & Carl Pomerance, "Prime Numbers", A2.1.7 or A9.5.26
• Ananda Mohan, "Residue Number Systems"

The question is, can we cram the moduli into the type like ModUInt does? In particular, how fast it will be?

Also check out RNS conversion based on core functions (see especially Abtahi et al, 2005, 10.1016/j.camwa.2005.03.008). Also see other ways in "Residue Number Systems"; and RNS to integer conversion using reconstruction coefficients in Montgomery form (allows us to keep inside the modulus range, but apparently quite slow as compared to straightforward reconstruction coefficient multiplication).

Specifically, for MLUInt. If we don't have to calculate the remainder/the quotient, will it be faster?

Currently it is calculated just as 0 - x, which can probably be improved.

Currently mulmod() uses divrem() which can be quite slow.

Two possible variants:

• modulus is a power of 2 (see HEAAN project)
• moduli with low numbers of significant bits
• moduli close to the typemax() (see e.g. efficient multiplication https://hal-lirmm.ccsd.cnrs.fr/file/index/docid/106470/filename/D547.PDF)

Currently it is a very crude implementation using Julia's bitstring(). Although it is called only during NTT plan creation, and the plans are cached, so its performance is not too important.

A number of paper works with "centralized" modulo numbers, that is for a modulus q they lie in [-q÷2, q÷2]. For example, see Shuhong Gao's FHE scheme (and SGFHE project).

Those would be useful to have, but one needs to check first if Julia provides any constant-time guarantees for its functions on builtin types. If not, this may be achieved by linking a C library, but that would defeat the purpose of easily-modifiable all-Julia implementation.

from_montgomery()'s performance depends on whether we use mulmod_montgomery(x, 1) or mulmod_montgomery(1, x) - the former works faster with MLUInt{3, UInt32}, and the latter with MLUInt{2, UInt64}. Need to run tests on different combinations, determine the pattern, and make several implementations.

The branch adk_mul contains an implementation of arbitrary degree Karatsuba (ADK), in polynomials.jl::adk_mul().

Sources:

1. "Missing a trick: Karatsuba variations" by M. Scott (https://miracl.com/assets/pdf-downloads/sk-1.pdf)
2. "Generalizations of the Karatsuba Algorithm for Efficient Implementations" by A. Weimerskirch and C. Paar (https://eprint.iacr.org/2006/224.pdf)

It is faster than mul_naive() and even mul_karatsuba for small polynomial sizes (still N^2 complexity though, so Karatsuba wins for larger sizes). It also does not require the polynomial size to be a power of 2.

It may be possible to get rid of the additional buffer requirement, in which case it will become even faster.

The title says it all.

One problem I can see is what if the highmost coefficient isn't a multiple of the highmost coefficient of the divisor? For example, (3x^4 + 1) ÷ (2x^2 + 1). Convert to floating point coefficients (possibly the same as Julia does it, that is have / and ÷ operations)? Disallow divisors with the highmost coefficient other than one (this will work for the purposes of polynomial moduli, at least, since that's the case for all of them)?

At the moment of writing it the primary focus were negacyclic convolutions, so nussbaumer_mul_cyclic() never got unrolled, but doing so will speed it up.

Currently limbs in MLUInt/MLInt are stored as big-endian, because that's how multi-limb algorithms are generally described in the literature (e.g. in TAOCP or Handbook of Applied Cryptography). Perhaps, little-endian would be a more convenient format - limbs are already little-endian (on most systems where this library actually works), and it will help visually compare MLUInt/MLInt numbers with builtin types, or initialize them with predefined constants.

Currently it is super-slow, even streamlined from mod(MLUInt{N, T}, MLUInt{N, T}) (10x the times for BigInt), but it would be very useful for RNS decomposition (issue #12). Need to see if we can lift NTL or https://github.com/niekbouman/ctbignum implementation.

According to Knuth it loops at most twice, so unrolling it may speed the function up.

See the LogProof project for an example.

Need to support constructions like

rand(MLUInt)
rand(MLInt)
rand(a:b) # where `a`, `b` are of type `MLUInt`/`MLInt`

Originally they were removed because of combinatorial explosion when new types are added, and currently one has to unpack ModUInt and MgModUInt using value() before conversion.

It may be possible to build a sensible conversion system that won't double in size when, e.g. ModInt and MgModInt are added (see issue #10).

Currently it transforms the argument from Montgomery representation first. Can we determine if the value is odd or not without the transformation?

See e.g. R. Crandall and C. Pomerance, "Prime Numbers", Algorithm 9.5.2

_benchmark_distribution() can be replaced by setup= keyword of @benchmark. Instead of standard deviation, standard error of the mean should be printed.

Benchmark results should be looked at after this replacement - there is some overhead in @benchmark, plus an incorrect choice between using symbols vs splicing arguments in this macro can change the results drastically.

To reduce the error, tests can be designed better - e.g. for divrem of multi-limb integers check several cases of n limbs by m limbs, instead of taking a random over the whole range (which, naturally ends up as maximum limbs by maximum limbs in most cases).

Now that Polynomial has a proper modulus field, support for arbitrary moduli (in addition to the current x^N - 1/x^N + 1) can be added. The Polynomial class may be split into two, non-modulo and modulo ones (same as we have now with integers).

Will need issue #6 implemented.

Currently implemented: MLUInt % built_in_type. Remaining:

• MLInt % built_in_type
• MLUInt/MLInt % MLUInt/MLInt
• built_in_type % MLUInt/MLInt

See the line a = submod(a, p2, m). submod() requires both of its first two arguments to be lower than m, but I am not sure how to prove that a satisfies that.

On the other hand, exhaustive tests with MLUInt{2} pass, so perhaps it's okay?

Useful, for example, in RNS decomposition (issue #12).