Any plans to incorporate polynomial rings modulo an irreversible polynomial?

[aklitzke]

Hi!

Great repo so far, I love the direction and believe this will be a very useful repo in cryptography! I've been looking at using feanor-math for secret sharing over finite polynomial rings (not fields), similar to what is done in https://eprint.iacr.org/2019/872 and a few other places. However, I'm running into a limitation. Hoping to check in with the dev(s) and see if I'm missing anything.

Consider for example the finite polynomial ring Z_4[x] / x^2+x+1. This ring is composed of elements (0+0x), (1+0x), ..., (3+3x). 4 is not prime, so elements generally cannot be inverted. However, this is not true for every element. Notice that while (2+2x) has no multiplicative inverse,(2+3x) DOES have a multiplicative inverse: (2+3x)^-1=(1+3x). To demonstrate:

(2+3x)(1+3x) = (2 + x + x^2)
(2 + x + x^2) mod (1 + x + x^2) = 1

What I'd like to be able to do with feanor-math is calculate the multiplicative inverse for elements in polynomial rings modulo some irreducible polynomial. So, given (2+3x) and (x^2+x+1) calculate (1+3x).

The way I've been doing this for non-polynomial Z_n fields is to use checked_left_div(1,...). However, doing checked_left_div(1,(2+3x)) on a PolyRing will simply return None. This makes sense, as the algorithm backing it does not consider the reduction modulo (x^2+x+1). Let me know if there's something I'm missing that would allow for this to be calculated, but it doesn't seem possible without changes to the library itself.

Adding polynomial rings modulo some irreducible polynomial doesn't seem terribly complicated to add. I toyed around with making a custom RingBase. However to calculate inverses you'd have to re-implement DivisibilityRing with a significantly more complex algorithm for checked_left_div, which started getting fairly in-the-weeds.

Let me know if there is something I'm missing that could enable this or if there are any plans to implement this behavior!

[FeanorTheElf]

Thank you for your post! The ring you are looking for basically exists as FreeAlgebra resp. FreeAlgebraImpl. However, it currently does not implement DivisibilityRing - mainly since the implementation is indeed quite involved. This is certainly planned for the future, and has quite high priority.

A few details on division in that ring:

• there are multiple ways to implement it. In the situation that the whole ring is finite and you only want to invert units, the easiest way seems to be to raise an element to the power of the size of the unit group minus one. This might be a valid workaround in your case (i.e. in Z_4[x] / (x^2+x+1), we always have for invertible a that a^-1 = a^3)
• An option that would work in any ring (but might be slower in some cases) would be based on linear algebra. This is the solution that I want to implement, but I have to decide on an approach to matrices and matrix inversion first - this is currently one of the most important missing features in the library
• There might be other ways, e.g. based on automorphisms or gcd-computations, that can be faster in certain cases

If this is an urgent requirement for you, I'd be happy to provide an ad-hoc implementation of the first approach, it should only be a few lines of code (I hope).

[aklitzke]

Thanks for the response! I took a look at FreeAlgebra, will play around with it! Good to hear you have divisibility planned for the future. I definitely appreciate your offer for an ad-hoc solution but your time is probably best spent on the right solution for the library. However a discussion of the math behind division in that ring would be helpful, if nothing else just for my own edification.

In my case the whole ring is finite and I do only want to invert known units, so a power solution could work. However I'm a bit confused about the underlying math and was hoping you could help out. Through some trial and error I see that for Z_4[x] / (x^2+x+1) any unit raised to the power of (6n-1 where n is any positive integer) will give its inverse. In other words, a^-1=a^5=a^11=a^17=...However, it isn't clear to me where the 6 is coming from or how to generalize this for other Z_p^k[x]/f rings.

You mentioned a^3, and at first that made intuitive sense as it is -1 in Z_4, but that didn't actually work when I tested it. You also mentioned 'size of the unit group', but there are 12 units in Z_4[x]/(x^2+x+1), not 6. I looked through some textbooks and papers, but struggled to answer why does a^(6n-1)=a^-1 for Z_4[x]/(x^2+x+1)

[FeanorTheElf]

The reason is that the set of units in R = Z_n[X]/(f) (or more generally in any ring) is a group (usually denoted by R*). Now the inverse of a unit is the inverse w.r.t. the group, and in any group, we have a^|R*| = 1 (basically this is Lagrange's theorem). This means that a * a^(|R*| - 1) = 1, so a^(|R*| - 1) is the inverse of a. Hence, it is left to find the order of the unit group R*, which is indeed 6 and not 3 (sorry, my mistake).

Edit The size of the unit group is 12. I should really do my calculations more carefully...

[aklitzke]

Got it. And I'm guessing ^(6n-1) works because, despite that particular ring having an order of 12, every element in the unit group happens to also be part of a smaller subgroup of order 6. 12 still works, 6 is just faster because of a fluke with this particular ring. Assuming this would not necessarily hold for other rings.

So we can raise elements to |R*|-1 to calculate inverses. So the question that remains is, how do we calculate |R*|?

I did some reading to try to answer this myself, saw on wiki and other sources that: Given R=GR[p^s, m], |R*| = G1 x G2 where for p=2 and s>=3 that G1=p^m-1 and G2=2*2^(s-2)*2^((s-1)*(m-1)). Doing some rough math for GR(2^8, 2), and that gives us 3*2^14. That's a big exponentiation!

My toolchain didn't work to exponentiate up to ^2^14, but I did some brute-forcing up to ~^2^10 and was able to determine a^(3*2^(7+n)-1) for any n>=0 works to invert GR[2^8, 2]. As 3*2^7 is a factor of |R*|=3*2^14, I'm guessing the R* can again be split into subgroups for this particular ring.

How might I do this for a degree of m=3? The same math for GR[2^8, 3] puts us at |R*| = 3*2^21. Seems very inefficient to perform that exponentiation. Is this where the other approaches you mentioned might come in? Or am I perhaps misunderstanding how to calculate |R*|?

[FeanorTheElf]

Yes, I think you are completely right regarding the 12 and 6. Your formulas for |R*| also look pretty solid.

However, I do not think that an exponentiation to 3*2^21 is difficult. Are you using FreeAlgebraImpl::pow()? That should do the job easily. The following code took less than a millisecond on my system:

let base_ring = Zn::new(256);
let modulo = base_ring.int_hom();
let ring = FreeAlgebraImpl::new(base_ring, [modulo.map(-1), modulo.map(-1)], default_memory_provider!());
let a = ring.from_canonical_basis([modulo.map(43), modulo.map(123)].into_iter());
let potential_inverse = ring.pow(ring.clone_el(&a), (3 << 21) - 1);
ring.println(&ring.mul_ref(&a, &potential_inverse));

Using RingStore::pow_gen(), you can even compute powers to exponents that go beyond the range of usize.

[aklitzke]

Oh, beautiful! Not realizing exponentiation could be that much faster, I was indeed doing rough prototyping in Mathematica. That solves that issue, thanks!

However, I'm now running into an issue using FreeAlgebraImpl. If a^n==1 then we'd expect a*a^(n-1)==1. However, that isn't the case using FreeAlgebraImpl over GR[2^8, 3] for some reason. Not sure if I set up the code wrong, or if I'm thinking about the math wrong? See the following minimal example:

Output:

is unit: (255θ + 255θ^2)^(640) => (1)
inverse: (255θ + 255θ^2)^(639) => (76 + 168θ + 33θ^2)
product: (255θ + 255θ^2)*(76 + 168θ + 33θ^2) => (201 + 158θ + 246θ^2)
thread 'minimal' panicked at src/main.rs:359:9:
Assertion failed: 1 != 201 + 158θ + 246θ^2

Code:

let base_ring = Zn::<256>::RING;
let modulo = base_ring.int_hom();
let ring = FreeAlgebraImpl::new(base_ring, [modulo.map(-1), modulo.map(-1), modulo.map(-1)], default_memory_provider!());
let a = ring.from_canonical_basis([modulo.map(0), modulo.map(255), modulo.map(255)].into_iter());
let n = 640;
let pow = ring.pow(ring.clone_el(&a), n);
if ring.eq_el(&pow, &ring.one()) {
    println!("is unit: ({})^({}) => ({})", ring.format(&a), n, ring.format(&pow));
    let inverse = ring.pow(ring.clone_el(&a), n-1);
    println!("inverse: ({})^({}) => ({})", ring.format(&a), n-1, ring.format(&inverse));
    let product = ring.mul_ref(&inverse, &a);
    println!("product: ({})*({}) => ({})", ring.format(&a), ring.format(&inverse), ring.format(&product));
    assert_el_eq!(&ring, &ring.one(), &product);
}

Worth noting when I check (255x+255x^2)^640 in Mathematica the product rule applies successfully, but I get a different result:

PolynomialMod[(255x+255x^2)^(640), x^3+x^2+1, Modulus->256]
121+191x+249x^2
PolynomialMod[(255 x+255 x^2)*(255 x+255 x^2)^(639), x^3+x^2+1, Modulus->256]
121+191x+249x^2

[FeanorTheElf]

That was indeed a bug in the multiplication (it only occurs for extensions of degree > 2). I have fixed it, and created a new version of feanor-math (1.4.0). It also includes division for FreeAlgebraImpl, so let me know if this is taken care of now!

Many thanks, I really enjoyed discussing this with you!

[aklitzke]

I've had the chance to throw more at the library and it's been working great for me for a few days now! I've been able to successfully do SSS over several polynomial rings, which runs the gamut of powers, multiplication, inversion, subtraction, and addition. Have used rings as large as GR(2^8,32) and it is still very fast. Thanks for fixing that bug and helping out with some of the math!