Optimize non native field mul for pseudo mersenne primes

[nicholas-mainardi]

This PR addresses issue #160 , introducing an optimization that allows to save from 30% to 40% constraints for the multiplication of non native field elements when the prime field is a pseudo-mersenne, i.e, the prime modulus p = 2^n - c for a small integer c.

This optimization allows to represent the result of a multiplication as a NonNativeFieldGagdet rather than a NonNativeFieldMulResultGadget, enabling a more efficient reduction. The private functions that deals with multiplication of field elements in a pseudo-mersenne field are added, and they are accordingly called by the publicly exposed multiplication functions in case the non native field is pseudo-mersenne. The definition of the prime field trait in algebra is modified to allow specifying the value of the small integer c for a pseudo-mersenne prime field.

This PR also fixes the estimation of the number of constraint for a multiplication performed in the find_parameters function, adding also the estimation in the case of a pseudo-mersenne field.

[DanieleDiBenedetto]

Will return on that to check the arithmetic, from a SW perspective maybe it's worth to explicitly separate the functions (and maybe put the ones for pseudo mersenne in another module ?)

[DanieleDiBenedetto]

@UlrichHaboeck75 what do you think ?

[UlrichHaboeck75]

I confirm. The factor 2 is missing. The correct formula should be

   surfeit' = log(2 * num_limbs * (num_add(L) + 1) * (num_add(R) + 1) + num_limbs).

This factor is introduced by 2^len(p) / p <= 2 in the limb bound of the quotient k.

[DanieleDiBenedetto]

Let's remove this TODO now no ?

[DanieleDiBenedetto]

One of the advantages of this optimization based on Mersenne is that it allows us to return a NonNativeFieldGadget as a result of mul instead of a NonNativeFieldMulResultGadget (on which we can only perform additions); here, however, to keep the interfaces consistent, we are still returning a NonNativeFieldMulResultGadget, obliging us to perform a reduce if we want to do another mul while it shouldn't be needed.
In future we will refactor everything so NonNativeFieldMulResultGadget won't exist anymore but, for this PR, would be nice to test if we are wasting constraints and how much.
I suggest to perform a test that does

let non_native_field_mul_result_gadget= a.mul(b)
let non_native_field_gadget = non_native_field_mul_result_gadget.reduce()
non_native_field_gadget .mul(c)

Using the function we expose to the user and another one:

let non_native_field_gadget = a.mul_mersenne(b)
non_native_field_gadget.mul_mersenne(c)

Using instead our internal function for Mersenne mul.
Let's count number of constraints in both cases and compare

[DanieleDiBenedetto]

As suspected, the strategy we are using now seems to be suboptimal. However, optimizing in a cleaner way requires most likely a huge refactoring (already envisioned by #152 anyway). I would suggest to postpone it when addressing #152.

[DanieleDiBenedetto]

Was it wrong before ? Why ?

[UlrichHaboeck75]

I am confused now, too. We implement the improvement using mul_by_constant(), which saves num_limbs(p) constraints. But we find an error in the previous computation, which neglects the number of constraints in of the conditional select, again num_limbs(p) many. Shouldn't we now end up with the old formula

     C = 3 * S + surfeit + num_limbs(p) + 1?

[nicholas-mainardi]

Yeah, you are right, I forgot to subtract num_limbs from the comment when I change the function to use mul_by_constant. Actually I realized that I have not yet updated the comment with the correct formulas I derived while fixing the find_parameters function, so I will update this comment with the correct estimated cost.

[nicholas-mainardi]

Done

[DanieleDiBenedetto]

you can simplify the ceiling in this way:

ceil(a/b) = (a + b - 1)/b

[DDT92]

This can be made a bit more precise:
If n is the number of additions over normal form of the initial gadget "x", then we have that
x < (n+1) * (2^bpl - 1)
which is equivalent to
x < { (n+1) / (2^bpl - 1) } * (2^bpl - 1)^2
It follows, called n' the new number of additions over normal form, n' = ceil ((n+1) / (2^bpl - 1) - 1)
In practice there will be no difference, since this will be always 0.

[nicholas-mainardi]

You are right about the n+1, I forgot it, thanks. But it is not clear why we should consider 2^bpl-1 as an upper bound for the limb rather than 2^bpl.

[UlrichHaboeck75]

I am not sure about the validity of the implemented formula. In terms of the limb bound factors num_adds + 1, shouldn't we have

   (num_adds + 1) = (num_adds(product) + 1) *  2^bits_per_limb,

and hence num_adds(product) = Floor[(num_adds + 1)/2^bits_per_limb] - 1?

[DDT92]

You are right about the n+1, I forgot it, thanks. But it is not clear why we should consider 2^bpl-1 as an upper bound for the limb rather than 2^bpl.

just because 2^bpl - 1 = (1...1)_2, but it does not make a big difference (in this case no difference in practice).

[DDT92]

I am not sure about the validity of the implemented formula. In terms of the limb bound factors num_adds + 1, shouldn't we have

   (num_adds + 1) = (num_adds(product) + 1) *  2^bits_per_limb,

and hence num_adds(product) = Floor[(num_adds + 1)/2^bits_per_limb] - 1?

num_adds(product) = Floor[(num_adds + 1)/2^bits_per_limb] - 1

This could be a negative number (in practice it will be -1), we need to put ceil instead of floor.

[nicholas-mainardi]

I would still use 2^bpl since, at least for a non-native field gagdet that is the result of a multiplication, n+1 should be divisible by 2^bpl, hence avoiding rounding.

[DDT92]

I would still use 2^bpl since, at least for a non-native field gagdet that is the result of a multiplication, n+1 should be divisible by 2^bpl, hence avoiding rounding.

It depends if you change the conversion from NonNativeFieldMulResultGadget -> NonNativeFieldGadget
according to my comment below. In that case (n+1) will be a multiple o (2^bpl - 1).

[nicholas-mainardi]

Uhm not sure, consider that n+1 is computed in mul_without_prereduce_for_pseudomersenne to be multiple of 2^bpl, it is not initialized from a NonNativeFieldMulResultGadget with the function you referred to.

[nicholas-mainardi]

Formula fixed in the code employing 2^bpl

[DanieleDiBenedetto]

Let's add a TODO about possibility to find some library supporting types bigger than 2^128 with const operation avoiding this additional complexity, or perhaps another solution at all @DDT92 @UlrichHaboeck75

[DanieleDiBenedetto]

Especially because this works only given our hardcoded choose of SURFEIT

[DDT92]

Let's add a TODO about possibility to find some library supporting types bigger than 2^128 with const operation avoiding this additional complexity, or perhaps another solution at all @DDT92 @UlrichHaboeck75

See my comment above

[UlrichHaboeck75]

The solution I propose here is the refactoring as described in #152, the different treatment of parameters (and hence the optimum finder).

[DanieleDiBenedetto]

As suspected, the strategy we are using now seems to be suboptimal. However, optimizing in a cleaner way requires most likely a huge refactoring (already envisioned by #152 anyway). I would suggest to postpone it when addressing #152.

[DanieleDiBenedetto]

Some minor comments, but overall fine for me

[DanieleDiBenedetto]

Let's remove this TODO now no ?

[UlrichHaboeck75]

Please update the comment.

[nicholas-mainardi]

I think the todo is still pending for non pseudo-mersenne fields, no?

[UlrichHaboeck75]

Right. Then let us clarify here?

[UlrichHaboeck75]

For the records: Let us add a reminder to overthink the treatment of the generic/special prime discrimination to issue #152. I wouldn't be surprised that we can do better once we drop mulresult gadgets.