Source file evaluations.ml
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
module Fr = Bls12_381.Fr
module Stubs = struct
type fr = Fr.t
type fr_array = Fr_carray.t
(** [add res a b size_a size_b] writes the result of polynomial addition
of [a] and [b] using the evaluation representation in [res], where
- [a] is evaluated on [domain_a] of size [size_a]
- [b] is evaluated on [domain_b] of size [size_b]
requires:
- [size a = size_a]
- [size b = size_b]
- [size res = min (size_a, size_b)]
- [res], [a] and [b] are either pairwise disjoint or equal
- [size_b mod size_a = 0] *)
external add : fr_array -> fr_array -> fr_array -> int -> int -> unit
= "caml_bls12_381_polynomial_polynomial_evaluations_add_stubs"
[@@noalloc]
(** [mul_arrays res eval_evallen_comp_power_powlen size_res nb_evals] writes
the result of computing [p₁(gᶜ₁·x)ᵐ₁·p₂(gᶜ₂·x)ᵐ₂·…·pₖ(gᶜₖ·x)ᵐₖ] using
the evaluation representation in [res], where
- [eval_evallen_comp_power_powlen.[i] = (pᵢ, size_p_i, cᵢ, mᵢ, size_bits_m_i)]
- a polynomial [pᵢ] is evaluated on [domain_p_i] of size [size_p_i]
- [cᵢ] is a composition_gx; it computes [pᵢ(gᶜᵢ·x)] instead of [pᵢ(x)],
where [g] is a primitive [size_p_i]-th root of unity
- [mᵢ] is a power in [pᵢ(x)ᵐᵢ]
- [size_bits_m_i] is the *exact* number of bits in [mᵢ]
- [nb_evals] is a number of evaluations, i.e., [i = 1..nb_evals]
requires:
- [size res = size_res]
- [size eval_evallen_comp_power_powlen = nb_evals]
- [size p_i = size_p_i]
- [size_bits m_i = size_bits_m_i]
- [size_p_i mod size_res = 0]
- [res] and [p_i] are disjoint *)
external mul_arrays :
fr_array ->
(fr_array * int * int * Bytes.t * int) array ->
int ->
int ->
unit = "caml_bls12_381_polynomial_polynomial_evaluations_mul_arrays_stubs"
[@@noalloc]
(** [linear_arrays res eval_evallen_coeff_comp add_constant size_res nb_evals]
writes the result of computing
[λ₁·p₁(gᶜ₁·x) + λ₂·p₂(gᶜ₂·x) + … + λₖ·pₖ(gᶜₖ·x) + add_constant] using
the evaluation representation in [res], where
- [eval_evallen_coeff_comp.[i] = (pᵢ, size_p_i, λᵢ, cᵢ)]
- a polynomial [pᵢ] is evaluated on [domain_p_i] of size [size_p_i]
- [cᵢ] is a composition_gx; it computes [pᵢ(gᶜᵢ·x)] instead of [pᵢ(x)],
where [g] is a primitive [size_p_i]-th root of unity
- [λᵢ] is a coefficient in [λᵢ·pᵢ(x)]
- [nb_evals] is a number of evaluations, i.e., [i = 1..nb_evals]
requires:
- [size res = size_res]
- [size eval_evallen_coeff_comp = nb_evals]
- [size p_i = size_p_i]
- [size_p_i mod size_res = 0]
- [res] and [p_i] are disjoint *)
external linear_arrays :
fr_array -> (fr_array * int * fr * int) array -> fr -> int -> int -> unit
= "caml_bls12_381_polynomial_polynomial_evaluations_linear_arrays_stubs"
[@@noalloc]
(** [fft_inplace p domain log log_degree] computes Fast Fourier Transform.
It converts the coefficient representation of a polynomial [p] to
the evaluation representation
requires:
- [size p = size domain]
- [size domain = 2^log]
- [domain = [one; g; ..; g^{n-1}]] where [g] is a primitive
[n]-th root of unity and [n = 2^log] (as done by {!Domain.Stubs.compute_domain}) *)
external fft_inplace :
fr_array -> domain:fr_array -> log:int -> log_degree:int -> unit
= "caml_bls12_381_polynomial_fft_inplace_on_stubs"
(** [ifft_inplace p domain log] computes Inverse Fast Fourier Transform.
It converts the evaluation representation of a polynomial [p] to
the coefficient representation
requires:
- [size p = size domain]
- [size domain = 2^log]
- [domain = [one; g; ..; g^{n-1}]] where [g] is a primitive
[n]-th root of unity and [n = 2^log] (as done by {!Domain.Stubs.compute_domain}) *)
external ifft_inplace : fr_array -> domain:fr_array -> log:int -> unit
= "caml_bls12_381_polynomial_ifft_inplace_on_stubs"
[@@noalloc]
(** [dft_inplace coefficients domain inverse length] computes the Fourier Transform.
requires:
- [size domain = size coefficients = length]
- [length <= 2^10]
- [length != 2^k]
- if [inverse = true] then the inverse dft is performed
- [domain = [one; g; ..; g^{n-1}]] where [g] is a primitive
[n]-th root of unity (as done by {!Domain.Stubs.compute_domain}) *)
external dft_inplace : fr_array -> fr_array -> bool -> int -> unit
= "caml_bls12_381_polynomial_dft_stubs"
[@@noalloc]
(** [fft_prime_factor_algorithm_inplace coefficient (domain1, domain1_length) (domain2, domain2_length) inverse]
computes the Fast Fourier Transform following
{{: https://en.wikipedia.org/wiki/Prime-factor_FFT_algorithm }the Prime-factor FFT algorithm}.
requires:
- [size domain1 = domain1_length]
- [size domain2 = domain2_length]
- [size domain1] and [size domain2] to be coprime
- if for some k [size domain1 != 2^k] then [size domain1 <= 2^10]
- if for some k [size domain2 != 2^k] then [size domain2 <= 2^10]
- [size coefficients = domain1_length * domain2_length]
- if [inverse = true] then the inverse fft is performed
- [domain = [one; g; ..; g^{n-1}]] where [g] is a primitive
[n]-th root of unity (as done by {!Domain.Stubs.compute_domain}) *)
external fft_prime_factor_algorithm_inplace :
fr_array -> fr_array * int -> fr_array * int -> bool -> unit
= "caml_bls12_381_polynomial_prime_factor_algorithm_fft_stubs"
[@@noalloc]
end
module type Evaluations_sig = sig
type scalar
type polynomial
type t [@@deriving repr]
(** [init size f degree] initializes an evaluation of a polynomial of the given
[degree]. *)
val init : int -> (int -> scalar) -> degree:int -> t
(** [of_array (d, e)] creates a value of type [t] from the evaluation
representation of a polynomial [e] of degree [d], i.e, it converts an OCaml
array to a C array *)
val of_array : int * scalar array -> t
(** [to_array] converts a C array to an OCaml array *)
val to_array : t -> scalar array
(** [string_of_eval e] returns the string representation of evaluation [e] *)
val string_of_eval : t -> string
type domain
(** [of_domain d] converts [d] from type {!type:domain} to type {!type:t}.
Note: [of_domain d] doesn't create a copy of [d] *)
val of_domain : domain -> t
(** [to_domain d] converts [d] from type {!type:t} to type {!type:domain}.
Note: [to_domain d] doesn't create a copy of [d] *)
val to_domain : t -> domain
(** [zero] returns the evaluation representation of the zero polynomial *)
val zero : t
(** [is_zero p] checks whether a polynomial [p] is the zero polynomial *)
val is_zero : t -> bool
(** [degree] returns the degree of a polynomial. Returns [-1] for the zero polynomial *)
val degree : t -> int
(** [length e] returns the size of domain where a polynomial is evaluated, or equally,
the size of a C array where evaluation [e] is stored *)
val length : t -> int
(** [create len] returns the evaluation representation of a zero polynomial of size [len] *)
val create : int -> t
(** [copy ?res a] returns a copy of evaluation [a]. The function writes the result in [res]
if [res] has the correct size and allocates a new array for the result otherwise *)
val copy : ?res:t -> t -> t
(** [get p i] returns the [i]-th element of a given array [p] *)
val get : t -> int -> scalar
(** [get_inplace p i res] copies the [i]-th element of a given array [p] in res *)
val get_inplace : t -> int -> scalar -> unit
(** [mul_by_scalar] computes muliplication of a polynomial by a blst_fr element *)
val mul_by_scalar : scalar -> t -> t
(** [mul_c] computes [p₁(gᶜ₁·x)ᵐ₁·p₂(gᶜ₂·x)ᵐ₂·…·pₖ(gᶜₖ·x)ᵐₖ], where
- [pᵢ = List.nth evaluations i]
- [mᵢ = List.nth powers i]
- [cᵢ = List.nth (fst composition_gx) i]
- [n = snd composition_gx] is the order of generator, i.e., [gⁿ = 1]
The function writes the result in [res] if [res] has the correct size (= min (size pᵢ))
and allocates a new array for the result otherwise
Note: [res] and [pᵢ] are disjoint *)
val mul_c :
?res:t ->
evaluations:t list ->
?composition_gx:int list * int ->
?powers:int list ->
unit ->
t
(** [linear_c] computes [λ₁·p₁(gᶜ₁·x) + λ₂·p₂(gᶜ₂·x) + … + λₖ·pₖ(gᶜₖ·x) + add_constant], where
- [pᵢ = List.nth evaluations i]
- [λᵢ = List.nth linear_coeffs i]
- [cᵢ = List.nth (fst composition_gx) i]
- [n = snd composition_gx] is the order of generator, i.e., [gⁿ = 1]
The function writes the result in [res] if [res] has the correct size (= min (size pᵢ))
and allocates a new array for the result otherwise
Note: [res] and [pᵢ] are disjoint *)
val linear_c :
?res:t ->
evaluations:t list ->
?linear_coeffs:scalar list ->
?composition_gx:int list * int ->
?add_constant:scalar ->
unit ->
t
(** [linear_with_powers p s] computes [∑ᵢ sⁱ·p.(i)]. This function is more efficient
than [linear] + [powers] for evaluations of the same size *)
val linear_with_powers : t list -> scalar -> t
(** [add ?res a b] computes polynomial addition of [a] and [b]. The function writes
the result in [res] if [res] has the correct size (= min (size (a, b))) and allocates
a new array for the result otherwise
Note: [res] can be equal to either [a] or [b] *)
val add : ?res:t -> t -> t -> t
(** [equal a b] checks whether a polynomial [a] is equal to a polynomial [b]
Note: [equal] is defined as restrictive equality, i.e., the same polynomial
evaluated on different domains are said to be different *)
val equal : t -> t -> bool
(** [evaluation_fft domain p] converts the coefficient representation of
a polynomial [p] to the evaluation representation. [domain] can be obtained
using {!Domain.build}
Note:
- size of domain must be a power of two
- degree of polynomial must be strictly less than the size of domain *)
val evaluation_fft : domain -> polynomial -> t
(** [interpolation_fft domain p] converts the evaluation representation of
a polynomial [p] to the coefficient representation. [domain] can be obtained
using {!Domain.build}
Note:
- size of domain must be a power of two
- size of a polynomial must be equal to size of domain *)
val interpolation_fft : domain -> t -> polynomial
val interpolation_fft2 : domain -> scalar array -> polynomial
(** [dft domain polynomial] converts the coefficient representation of
a polynomial [p] to the evaluation representation. [domain] can be obtained
using {!Domain.build}.
Computes the discrete Fourier transform in time O(n^2)
where [n = size domain].
requires:
- [size domain] to divide Bls12_381.Fr.order - 1
- [size domain != 2^k]
- [degree polynomial < size domain] *)
val dft : domain -> polynomial -> t
(** [idft domain t] converts the evaluation representation of
a polynomial [p] to the coefficient representation. [domain] can be obtained
using {!Domain.build}.
Computes the inverse discrete Fourier transform in time O(n^2)
where [n = size domain].
requires:
- [size domain] to divide Bls12_381.Fr.order - 1
- [size domain != 2^k]
- [size domain = size t] *)
val idft : domain -> t -> polynomial
(** [evaluation_fft_prime_factor_algorithm domain1 domain2 p] converts the
coefficient representation of a polynomial [p] to the evaluation representation.
[domain] can be obtained using {!Domain.build}.
See {{: https://en.wikipedia.org/wiki/Prime-factor_FFT_algorithm }the Prime-factor FFT algorithm}.
requires:
- [size domain1 * size domain2] to divide Bls12_381.Fr.order - 1
- [size domain1] and [size domain2] to be coprime
- if for some k [size domain1 != 2^k] then [size domain1 <= 2^10]
- if for some k [size domain2 != 2^k] then [size domain2 <= 2^10]
- [degree polynomial < size domain1 * size domain2] *)
val evaluation_fft_prime_factor_algorithm :
domain1:domain -> domain2:domain -> polynomial -> t
(** [interpolation_fft_prime_factor_algorithm domain1 domain2 t] converts the
evaluation representation of a polynomial [p] to the coefficient representation.
[domain] can be obtained using {!Domain.build}.
See {{: https://en.wikipedia.org/wiki/Prime-factor_FFT_algorithm }the Prime-factor FFT algorithm}.
requires:
- [size domain1 * size domain2] to divide Bls12_381.Fr.order - 1
- [size domain1] and [size domain2] to be coprime
- if for some k [size domain1 != 2^k] then [size domain1 <= 2^10]
- if for some k [size domain2 != 2^k] then [size domain2 <= 2^10]
- [size t = size domain1 * size domain2] *)
val interpolation_fft_prime_factor_algorithm :
domain1:domain -> domain2:domain -> t -> polynomial
end
module Evaluations_impl = struct
module Domain_c = Domain.Stubs
module Domain = Domain.Domain_unsafe
module Polynomial_c = Polynomial.Stubs
module Polynomial = Polynomial.Polynomial_unsafe
type scalar = Fr.t
type polynomial = Polynomial.t
type t = int * Fr_carray.t [@@deriving repr]
type domain = Domain.t
let init size f ~degree = (degree, Fr_carray.init size f)
let of_array (d, p) =
if d < -1 then
raise @@ Invalid_argument "make_evaluation: degree must be >= -1" ;
if Array.length p <= d then
raise
@@ Invalid_argument "make_evaluation: array must be longer than degree" ;
let res = Fr_carray.of_array p in
(d, res)
let to_array (_d, e) = Fr_carray.to_array e
let to_carray (_, e) = e
let string_of_eval (d, e) =
Printf.sprintf "%d : [%s]" d Polynomial.(to_string (of_carray e))
let of_domain domain =
let d = Domain.to_carray domain in
(1, d)
let allocate = Fr_carray.allocate
let to_domain (_, eval) = Domain.of_carray eval
let zero = (-1, allocate 1)
let degree (d, _) = d
let length (_, e) = Fr_carray.length e
let create n = (-1, allocate n)
let is_zero (d, _e) =
d = -1
let allocate_for_res res length_result =
match res with
| Some (_, res) when Fr_carray.length res = length_result -> res
| _ -> allocate length_result
let copy ?res (d, e) =
let len = Fr_carray.length e in
let res = allocate_for_res res len in
Fr_carray.blit e ~src_off:0 res ~dst_off:0 ~len ;
(d, res)
let get (_, eval) i = Fr_carray.get eval i
let get_inplace (_, eval) i res = Fr_carray.get_inplace eval i res
let mul_by_scalar lambda (d, e) =
let len = Fr_carray.length e in
let res = allocate len in
Polynomial_c.mul_by_scalar res lambda e len ;
(d, res)
let mul_c ?res ~evaluations ?composition_gx ?powers () =
let len_evaluations = List.length evaluations in
let composition_gx =
match composition_gx with
| Some composition_gx -> composition_gx
| None -> (List.init len_evaluations (fun _ -> 0), 1)
in
let powers =
match powers with
| Some powers -> powers
| None -> List.init len_evaluations (fun _ -> 1)
in
let domain_len = snd composition_gx in
assert (domain_len > 0) ;
assert (len_evaluations > 0) ;
assert (List.compare_length_with (fst composition_gx) len_evaluations = 0) ;
assert (List.compare_length_with powers len_evaluations = 0) ;
assert (List.for_all (fun power -> power > 0) powers) ;
let length_result =
List.fold_left min Int.max_int @@ List.map length evaluations
in
let res = allocate_for_res res length_result in
if List.exists is_zero evaluations then (
Fr_carray.erase res length_result ;
(-1, res))
else
let degree_result =
List.fold_left2
(fun acc d pow -> acc + (d * pow))
0
(List.map degree evaluations)
powers
in
if degree_result >= length_result then
raise
(Invalid_argument
(Printf.sprintf
"mul_evaluations : evaluation is too short (length=%d) for \
expected result size %d"
length_result
(degree_result + 1)))
else
let list_array =
List.map2
(fun (evaluation, pow) composition ->
let pow = Z.of_int pow in
let pow_bytes = Z.to_bits pow |> Bytes.unsafe_of_string in
let pow_len = Z.numbits pow in
let l = length evaluation in
let rescale_composition = composition * l / domain_len in
(snd evaluation, l, rescale_composition, pow_bytes, pow_len))
(List.combine evaluations powers)
(fst composition_gx)
in
let nb_evals = List.length evaluations in
Stubs.mul_arrays res (Array.of_list list_array) length_result nb_evals ;
(degree_result, res)
let constant p c =
for i = 0 to Fr_carray.length p - 1 do
Fr_carray.set p c i
done
let linear_c ?res ~evaluations ?linear_coeffs ?composition_gx
?(add_constant = Fr.zero) () =
let len_evaluations = List.length evaluations in
let linear_coeffs =
match linear_coeffs with
| Some linear_coeffs -> linear_coeffs
| None -> List.init len_evaluations (fun _ -> Fr.(copy one))
in
let composition_gx =
match composition_gx with
| Some composition_gx -> composition_gx
| None -> (List.init len_evaluations (fun _ -> 0), 1)
in
let domain_len = snd composition_gx in
assert (domain_len > 0) ;
assert (len_evaluations > 0) ;
assert (List.compare_length_with linear_coeffs len_evaluations = 0) ;
assert (List.compare_length_with (fst composition_gx) len_evaluations = 0) ;
let list_eval_coeff_composition =
List.map2
(fun (eval, coeff) composition ->
let rescale_composition = composition * length eval / domain_len in
(eval, coeff, rescale_composition))
(List.combine evaluations linear_coeffs)
(fst composition_gx)
|> List.filter (fun (eval, _, _) -> not (is_zero eval))
in
match list_eval_coeff_composition with
| [] ->
let length_result =
List.fold_left min Int.max_int @@ List.map length evaluations
in
let res = allocate_for_res res length_result in
constant res add_constant ;
let degree_result = if Fr.is_zero add_constant then -1 else 0 in
(degree_result, res)
| _ :: _ ->
let length_result =
List.fold_left min Int.max_int
@@ List.map
(fun (eval, _, _) -> length eval)
list_eval_coeff_composition
in
let degree_result =
List.fold_left max 0
@@ List.map
(fun (eval, _, _) -> degree eval)
list_eval_coeff_composition
in
let nb_evals = List.length list_eval_coeff_composition in
let array_eval_coeff_composition =
List.map
(fun (eval, linear_coeff, composition) ->
(snd eval, length eval, linear_coeff, composition))
list_eval_coeff_composition
|> Array.of_list
in
let res = allocate_for_res res length_result in
Stubs.linear_arrays
res
array_eval_coeff_composition
add_constant
length_result
nb_evals ;
(degree_result, res)
let add ?res e1 e2 =
let d1 = fst e1 in
let d2 = fst e2 in
let l1 = length e1 in
let l2 = length e2 in
if d1 = -1 then
let res = allocate_for_res res l2 in
copy ~res:(d2, res) e2
else if d2 = -1 then
let res = allocate_for_res res l1 in
copy ~res:(d1, res) e1
else
let deg_result = max d1 d2 in
let length_result = min l1 l2 in
let res = allocate_for_res res length_result in
Stubs.add res (snd e1) (snd e2) l1 l2 ;
(deg_result, res)
let linear_with_powers evals coeff =
let nb_evals = List.length evals in
assert (nb_evals > 0) ;
let eval_lenghts = List.map length evals in
let eval0_length = List.hd eval_lenghts in
let is_equal_size = List.for_all (Int.equal eval0_length) eval_lenghts in
if is_equal_size then (
let length_result = eval0_length in
let degree_result = List.fold_left max (-1) @@ List.map degree evals in
if degree_result = -1 then create length_result
else
let res = allocate length_result in
let evals =
List.map (fun (_, e) -> (e, Fr_carray.length e)) evals
|> Array.of_list
in
Polynomial_c.linear_with_powers res coeff evals nb_evals ;
(degree_result, res))
else
let coeffs = Fr_carray.powers nb_evals coeff |> Array.to_list in
linear_c ~evaluations:evals ~linear_coeffs:coeffs ()
let equal (deg1, eval1) (deg2, eval2) =
if deg1 <> deg2 || Fr_carray.(length eval1 <> length eval2) then false
else Polynomial.(equal (of_carray eval1) (of_carray eval2))
let evaluation_fft_internal : Domain.t -> polynomial -> Fr_carray.t =
fun domain poly ->
let degree = Polynomial.degree poly in
let log_degree = Z.log2up (Z.of_int (degree + 1)) in
let domain = Domain.to_carray domain in
let n_domain = Fr_carray.length domain in
let poly = Polynomial.to_carray poly in
let log = Z.(log2up @@ of_int n_domain) in
if not (Helpers.is_power_of_two n_domain) then
raise @@ Invalid_argument "Size of domain should be a power of 2." ;
if not (degree < n_domain) then
raise
@@ Invalid_argument
"Degree of poly should be strictly less than domain size." ;
let res = allocate n_domain in
Fr_carray.blit poly ~src_off:0 res ~dst_off:0 ~len:(degree + 1) ;
Stubs.fft_inplace res ~domain ~log ~log_degree ;
res
let evaluation_fft : domain -> polynomial -> t =
fun domain poly ->
let d = Polynomial.degree poly in
let domain_length = Domain.length domain in
if d = -1 then (d, allocate domain_length)
else
let res = evaluation_fft_internal domain poly in
(d, res)
let interpolation_fft_internal : Domain.t -> Fr_carray.t -> polynomial =
fun domain coefficients ->
let domain = Domain.to_carray domain in
let n_domain = Fr_carray.length domain in
let log = Z.(log2up @@ of_int n_domain) in
if not (Helpers.is_power_of_two n_domain) then
raise @@ Invalid_argument "Size of domain should be a power of 2." ;
let n_coefficients = Fr_carray.length coefficients in
if not (n_coefficients = n_domain) then
raise @@ Invalid_argument "Size of coefficients should be same as domain." ;
Stubs.ifft_inplace coefficients ~domain ~log ;
Polynomial.of_carray coefficients
let interpolation_fft : domain -> t -> polynomial =
fun domain (d, evaluation) ->
if d = -1 then Polynomial.zero
else
let length_res = Domain.length domain in
let rescaled_eval = allocate length_res in
Domain_c.rescale
rescaled_eval
evaluation
length_res
(Fr_carray.length evaluation) ;
interpolation_fft_internal domain rescaled_eval
let interpolation_fft2 : Domain.t -> scalar array -> polynomial =
fun domain coefficients ->
interpolation_fft_internal domain (Fr_carray.of_array coefficients)
let dft domain polynomial =
let length = Domain.length domain in
if length > 1 lsl 10 then
raise @@ Invalid_argument "Domain size must be <= 2^10." ;
if Helpers.is_power_of_two length then
raise @@ Invalid_argument "Domain size must not be a power of two" ;
let d = Polynomial.degree polynomial in
let polynomial = Polynomial.to_carray polynomial in
if not (d < length) then
raise
@@ Invalid_argument
"Degree of poly should be strictly less than domain size." ;
let evaluations = allocate length in
Fr_carray.blit polynomial ~src_off:0 evaluations ~dst_off:0 ~len:(d + 1) ;
Stubs.dft_inplace evaluations (Domain.to_carray domain) false length ;
(d, evaluations)
let idft domain (_, evaluations) =
let length = Domain.length domain in
if length > 1 lsl 10 then
raise @@ Invalid_argument "Domain size must be <= 2^10." ;
if Helpers.is_power_of_two length then
raise @@ Invalid_argument "Domain size must not be a power of two" ;
if not (length = Fr_carray.length evaluations) then
raise @@ Invalid_argument "Size of coefficients should be same as domain." ;
let coefficients = Fr_carray.copy evaluations in
Stubs.dft_inplace coefficients (Domain.to_carray domain) true length ;
Polynomial.of_carray coefficients
let evaluation_fft_prime_factor_algorithm ~domain1 ~domain2 polynomial =
let domain1_length = Domain.length domain1 in
let domain2_length = Domain.length domain2 in
if
(not (Helpers.is_power_of_two domain1_length))
&& domain1_length > 1 lsl 10
then
raise
@@ Invalid_argument "Domain of non power of 2 length must be <= 2^10." ;
if
(not (Helpers.is_power_of_two domain2_length))
&& domain2_length > 1 lsl 10
then
raise
@@ Invalid_argument "Domain of non power of 2 length must be <= 2^10." ;
if not Z.(gcd (of_int domain1_length) (of_int domain2_length) = one) then
raise @@ Invalid_argument "Size of domains must be coprime." ;
let n_domain = domain1_length * domain2_length in
let d = Polynomial.degree polynomial in
let coefficients = Polynomial.to_carray polynomial in
if not (d < n_domain) then
raise
@@ Invalid_argument
"Degree of poly should be strictly less than domain size." ;
let res = allocate n_domain in
if d = -1 then (d, res)
else
let domain1 = Domain.to_carray domain1 in
let domain2 = Domain.to_carray domain2 in
Fr_carray.blit coefficients ~src_off:0 res ~dst_off:0 ~len:(d + 1) ;
Stubs.fft_prime_factor_algorithm_inplace
res
(domain1, domain1_length)
(domain2, domain2_length)
false ;
(d, res)
let interpolation_fft_prime_factor_algorithm ~domain1 ~domain2 (d, evaluations)
=
let domain1_length = Domain.length domain1 in
let domain2_length = Domain.length domain2 in
if
(not (Helpers.is_power_of_two domain1_length))
&& domain1_length > 1 lsl 10
then
raise
@@ Invalid_argument "Domain of non power of 2 length must be <= 2^10." ;
if
(not (Helpers.is_power_of_two domain2_length))
&& domain2_length > 1 lsl 10
then
raise
@@ Invalid_argument "Domain of non power of 2 length must be <= 2^10." ;
if not Z.(gcd (of_int domain1_length) (of_int domain2_length) = one) then
raise @@ Invalid_argument "Size of domains must be coprime." ;
let n_domain = domain1_length * domain2_length in
let n_evaluations = Fr_carray.length evaluations in
if not (n_evaluations = n_domain) then
raise @@ Invalid_argument "Size of coefficients should be same as domain." ;
if d = -1 then Polynomial.zero
else
let domain1 = Domain.to_carray domain1 in
let domain2 = Domain.to_carray domain2 in
let coefficients = Fr_carray.copy evaluations in
Stubs.fft_prime_factor_algorithm_inplace
coefficients
(domain1, domain1_length)
(domain2, domain2_length)
true ;
Polynomial.of_carray coefficients
end
module type Evaluations_unsafe_sig = sig
include Evaluations_sig
(** [to_carray t] converts [t] from type {!type:t} to type {!type:Fr_carray.t}
Note: [to_carray t] doesn't create a copy of [t] *)
val to_carray : t -> Fr_carray.t
end
module Evaluations_unsafe :
Evaluations_unsafe_sig
with type scalar = Bls12_381.Fr.t
and type domain = Domain.t
and type polynomial = Polynomial.t =
Evaluations_impl
include (
Evaluations_unsafe :
Evaluations_sig
with type t = Evaluations_unsafe.t
and type scalar = Evaluations_unsafe.scalar
and type domain = Evaluations_unsafe.domain
and type polynomial = Evaluations_unsafe.polynomial)