1//! Allocation-free, (best-effort) constant-time, finite field arithmetic for large integers.
2//!
3//! Unlike `std.math.big`, these integers have a fixed maximum length and are only designed to be used for modular arithmetic.
4//! Arithmetic operations are meant to run in constant-time for a given modulus, making them suitable for cryptography.
5//!
6//! Parts of that code was ported from the BSD-licensed crypto/internal/bigmod/nat.go file in the Go language, itself inspired from BearSSL.
7
8const std = @import("std");
9const builtin = @import("builtin");
10const crypto = std.crypto;
11const math = std.math;
12const mem = std.mem;
13const testing = std.testing;
14const assert = std.debug.assert;
15const Endian = std.builtin.Endian;
16
17// A Limb is a single digit in a big integer.
18const Limb = usize;
19
20// The number of reserved bits in a Limb.
21const carry_bits = 1;
22
23// The number of active bits in a Limb.
24const t_bits: usize = @bitSizeOf(Limb) - carry_bits;
25
26// A TLimb is a Limb that is truncated to t_bits.
27const TLimb = @Int(.unsigned, t_bits);
28
29const native_endian = builtin.target.cpu.arch.endian();
30
31// A WideLimb is a Limb that is twice as wide as a normal Limb.
32const WideLimb = struct {
33 hi: Limb,
34 lo: Limb,
35};
36
37/// Value is too large for the destination.
38pub const OverflowError = error{Overflow};
39
40/// Invalid modulus. Modulus must be odd.
41pub const InvalidModulusError = error{ EvenModulus, ModulusTooSmall };
42
43/// Exponentiation with a null exponent.
44/// Exponentiation in cryptographic protocols is almost always a sign of a bug which can lead to trivial attacks.
45/// Therefore, this module returns an error when a null exponent is encountered, encouraging applications to handle this case explicitly.
46pub const NullExponentError = error{NullExponent};
47
48/// Invalid field element for the given modulus.
49pub const FieldElementError = error{NonCanonical};
50
51/// Invalid representation (Montgomery vs non-Montgomery domain.)
52pub const RepresentationError = error{UnexpectedRepresentation};
53
54/// The set of all possible errors `std.crypto.ff` functions can return.
55pub const Error = OverflowError || InvalidModulusError || NullExponentError || FieldElementError || RepresentationError;
56
57/// An unsigned big integer with a fixed maximum size (`max_bits`), suitable for cryptographic operations.
58/// Unless side-channels mitigations are explicitly disabled, operations are designed to be constant-time.
59pub fn Uint(comptime max_bits: comptime_int) type {
60 comptime assert(@bitSizeOf(Limb) % 8 == 0); // Limb size must be a multiple of 8
61
62 return struct {
63 const Self = @This();
64 const max_limbs_count = @divCeil(max_bits, t_bits);
65
66 limbs_buffer: [max_limbs_count]Limb,
67 /// The number of active limbs.
68 limbs_len: usize,
69
70 /// Number of bytes required to serialize an integer.
71 pub const encoded_bytes = @divCeil(max_bits, 8);
72
73 /// Constant slice of active limbs.
74 fn limbsConst(self: *const Self) []const Limb {
75 return self.limbs_buffer[0..self.limbs_len];
76 }
77
78 /// Mutable slice of active limbs.
79 fn limbs(self: *Self) []Limb {
80 return self.limbs_buffer[0..self.limbs_len];
81 }
82
83 // Removes limbs whose value is zero from the active limbs.
84 fn normalize(self: Self) Self {
85 var res = self;
86 if (self.limbs_len < 2) {
87 return res;
88 }
89 var i = self.limbs_len - 1;
90 while (i > 0 and res.limbsConst()[i] == 0) : (i -= 1) {}
91 res.limbs_len = i + 1;
92 assert(res.limbs_len <= res.limbs_buffer.len);
93 return res;
94 }
95
96 /// The zero integer.
97 pub const zero: Self = .{
98 .limbs_buffer = @splat(0),
99 .limbs_len = max_limbs_count,
100 };
101
102 /// Creates a new big integer from a primitive type.
103 /// This function may not run in constant time.
104 pub fn fromPrimitive(comptime T: type, init_value: T) OverflowError!Self {
105 var x = init_value;
106 var out: Self = .{
107 .limbs_buffer = undefined,
108 .limbs_len = max_limbs_count,
109 };
110 for (&out.limbs_buffer) |*limb| {
111 limb.* = if (@bitSizeOf(T) > t_bits) @as(TLimb, @truncate(x)) else x;
112 x = math.shr(T, x, t_bits);
113 }
114 if (x != 0) {
115 return error.Overflow;
116 }
117 return out;
118 }
119
120 /// Converts a big integer to a primitive type.
121 /// This function may not run in constant time.
122 pub fn toPrimitive(self: Self, comptime T: type) OverflowError!T {
123 var x: T = 0;
124 var i = self.limbs_len - 1;
125 while (true) : (i -= 1) {
126 if (@bitSizeOf(T) >= t_bits and math.shr(T, x, @bitSizeOf(T) - t_bits) != 0) {
127 return error.Overflow;
128 }
129 x = math.shl(T, x, t_bits);
130 const v = math.cast(T, self.limbsConst()[i]) orelse return error.Overflow;
131 x |= v;
132 if (i == 0) break;
133 }
134 return x;
135 }
136
137 /// Encodes a big integer into a byte array.
138 pub fn toBytes(self: Self, bytes: []u8, comptime endian: Endian) OverflowError!void {
139 if (bytes.len == 0) {
140 if (self.isZero()) return;
141 return error.Overflow;
142 }
143 @memset(bytes, 0);
144 var shift: usize = 0;
145 var out_i: usize = switch (endian) {
146 .big => bytes.len - 1,
147 .little => 0,
148 };
149 for (0..self.limbs_len) |i| {
150 var remaining_bits = t_bits;
151 var limb = self.limbsConst()[i];
152 while (remaining_bits >= 8) {
153 bytes[out_i] |= math.shl(u8, @as(u8, @truncate(limb)), shift);
154 const consumed = 8 - shift;
155 limb >>= @as(u4, @truncate(consumed));
156 remaining_bits -= consumed;
157 shift = 0;
158 switch (endian) {
159 .big => {
160 if (out_i == 0) {
161 if (i != self.limbs_len - 1 or limb != 0) {
162 return error.Overflow;
163 }
164 return;
165 }
166 out_i -= 1;
167 },
168 .little => {
169 out_i += 1;
170 if (out_i == bytes.len) {
171 if (i != self.limbs_len - 1 or limb != 0) {
172 return error.Overflow;
173 }
174 return;
175 }
176 },
177 }
178 }
179 bytes[out_i] |= @as(u8, @truncate(limb));
180 shift = remaining_bits;
181 }
182 }
183
184 /// Creates a new big integer from a byte array.
185 pub fn fromBytes(bytes: []const u8, comptime endian: Endian) OverflowError!Self {
186 if (bytes.len == 0) return Self.zero;
187 var shift: usize = 0;
188 var out = Self.zero;
189 var out_i: usize = 0;
190 var i: usize = switch (endian) {
191 .big => bytes.len - 1,
192 .little => 0,
193 };
194 while (true) {
195 const bi = bytes[i];
196 out.limbs()[out_i] |= math.shl(Limb, bi, shift);
197 shift += 8;
198 if (shift >= t_bits) {
199 shift -= t_bits;
200 out.limbs()[out_i] = @as(TLimb, @truncate(out.limbs()[out_i]));
201 const overflow = math.shr(Limb, bi, 8 - shift);
202 out_i += 1;
203 if (out_i >= out.limbs_len) {
204 if (overflow != 0 or i != 0) {
205 return error.Overflow;
206 }
207 break;
208 }
209 out.limbs()[out_i] = overflow;
210 }
211 switch (endian) {
212 .big => {
213 if (i == 0) break;
214 i -= 1;
215 },
216 .little => {
217 i += 1;
218 if (i == bytes.len) break;
219 },
220 }
221 }
222 return out;
223 }
224
225 /// Returns `true` if both integers are equal.
226 pub fn eql(x: Self, y: Self) bool {
227 return crypto.timing_safe.eql([max_limbs_count]Limb, x.limbs_buffer, y.limbs_buffer);
228 }
229
230 /// Compares two integers.
231 pub fn compare(x: Self, y: Self) math.Order {
232 return crypto.timing_safe.compare(
233 Limb,
234 x.limbsConst(),
235 y.limbsConst(),
236 .little,
237 );
238 }
239
240 /// Returns `true` if the integer is zero.
241 pub fn isZero(x: Self) bool {
242 var t: Limb = 0;
243 for (x.limbsConst()) |elem| {
244 t |= elem;
245 }
246 return ct.eql(t, 0);
247 }
248
249 /// Returns `true` if the integer is odd.
250 pub fn isOdd(x: Self) bool {
251 return @as(u1, @truncate(x.limbsConst()[0])) != 0;
252 }
253
254 /// Adds `y` to `x`, and returns `true` if the operation overflowed.
255 pub fn addWithOverflow(x: *Self, y: Self) u1 {
256 return x.conditionalAddWithOverflow(true, y);
257 }
258
259 /// Subtracts `y` from `x`, and returns `true` if the operation overflowed.
260 pub fn subWithOverflow(x: *Self, y: Self) u1 {
261 return x.conditionalSubWithOverflow(true, y);
262 }
263
264 // Replaces the limbs of `x` with the limbs of `y` if `on` is `true`.
265 fn cmov(x: *Self, on: bool, y: Self) void {
266 for (x.limbs(), y.limbsConst()) |*x_limb, y_limb| {
267 x_limb.* = ct.select(on, y_limb, x_limb.*);
268 }
269 }
270
271 // Adds `y` to `x` if `on` is `true`, and returns `true` if the
272 // operation overflowed.
273 fn conditionalAddWithOverflow(x: *Self, on: bool, y: Self) u1 {
274 var carry: u1 = 0;
275 for (x.limbs(), y.limbsConst()) |*x_limb, y_limb| {
276 const res = x_limb.* + y_limb + carry;
277 x_limb.* = ct.select(on, @as(TLimb, @truncate(res)), x_limb.*);
278 carry = @truncate(res >> t_bits);
279 }
280 return carry;
281 }
282
283 // Subtracts `y` from `x` if `on` is `true`, and returns `true` if the
284 // operation overflowed.
285 fn conditionalSubWithOverflow(x: *Self, on: bool, y: Self) u1 {
286 var borrow: u1 = 0;
287 for (x.limbs(), y.limbsConst()) |*x_limb, y_limb| {
288 const res = x_limb.* -% y_limb -% borrow;
289 x_limb.* = ct.select(on, @as(TLimb, @truncate(res)), x_limb.*);
290 borrow = @truncate(res >> t_bits);
291 }
292 return borrow;
293 }
294 };
295}
296
297/// A field element.
298fn Fe_(comptime bits: comptime_int) type {
299 return struct {
300 const Self = @This();
301
302 const FeUint = Uint(bits);
303
304 /// The element value as a `Uint`.
305 v: FeUint,
306
307 /// `true` if the element is in Montgomery form.
308 montgomery: bool = false,
309
310 /// The maximum number of bytes required to encode a field element.
311 pub const encoded_bytes = FeUint.encoded_bytes;
312
313 // The number of active limbs to represent the field element.
314 fn limbs_count(self: Self) usize {
315 return self.v.limbs_len;
316 }
317
318 /// Creates a field element from a primitive.
319 /// This function may not run in constant time.
320 pub fn fromPrimitive(comptime T: type, m: Modulus(bits), x: T) (OverflowError || FieldElementError)!Self {
321 comptime assert(@bitSizeOf(T) <= bits); // Primitive type is larger than the modulus type.
322 const v = try FeUint.fromPrimitive(T, x);
323 var fe = Self{ .v = v };
324 try m.shrink(&fe);
325 try m.rejectNonCanonical(fe);
326 return fe;
327 }
328
329 /// Converts the field element to a primitive.
330 /// This function may not run in constant time.
331 /// Returns an error if the element is in Montgomery form.
332 pub fn toPrimitive(self: Self, comptime T: type) (OverflowError || RepresentationError)!T {
333 if (self.montgomery) {
334 return error.UnexpectedRepresentation;
335 }
336 return self.v.toPrimitive(T);
337 }
338
339 /// Creates a field element from a byte string.
340 pub fn fromBytes(m: Modulus(bits), bytes: []const u8, comptime endian: Endian) (OverflowError || FieldElementError)!Self {
341 const v = try FeUint.fromBytes(bytes, endian);
342 var fe = Self{ .v = v };
343 try m.shrink(&fe);
344 try m.rejectNonCanonical(fe);
345 return fe;
346 }
347
348 /// Converts the field element to a byte string.
349 /// Returns an error if the element is in Montgomery form.
350 pub fn toBytes(self: Self, bytes: []u8, comptime endian: Endian) (OverflowError || RepresentationError)!void {
351 if (self.montgomery) {
352 return error.UnexpectedRepresentation;
353 }
354 return self.v.toBytes(bytes, endian);
355 }
356
357 /// Returns `true` if the field elements are equal, in constant time.
358 pub fn eql(x: Self, y: Self) bool {
359 return x.v.eql(y.v);
360 }
361
362 /// Compares two field elements in constant time.
363 pub fn compare(x: Self, y: Self) math.Order {
364 return x.v.compare(y.v);
365 }
366
367 /// Returns `true` if the element is zero.
368 pub fn isZero(self: Self) bool {
369 return self.v.isZero();
370 }
371
372 /// Returns `true` is the element is odd.
373 pub fn isOdd(self: Self) bool {
374 return self.v.isOdd();
375 }
376 };
377}
378
379/// A modulus, defining a finite field.
380/// All operations within the field are performed modulo this modulus, without heap allocations.
381/// `max_bits` represents the number of bits in the maximum value the modulus can be set to.
382pub fn Modulus(comptime max_bits: comptime_int) type {
383 return struct {
384 const Self = @This();
385
386 /// A field element, representing a value within the field defined by this modulus.
387 pub const Fe = Fe_(max_bits);
388
389 const FeUint = Fe.FeUint;
390
391 /// The neutral element.
392 zero: Fe,
393
394 /// The modulus value.
395 v: FeUint,
396
397 /// R^2 for the Montgomery representation.
398 rr: Fe,
399 /// Inverse of the first limb
400 m0inv: Limb,
401 /// Number of leading zero bits in the modulus.
402 leading: usize,
403
404 // Number of active limbs in the modulus.
405 fn limbs_count(self: Self) usize {
406 return self.v.limbs_len;
407 }
408
409 /// Actual size of the modulus, in bits.
410 pub fn bits(self: Self) usize {
411 return self.limbs_count() * t_bits - self.leading;
412 }
413
414 /// Returns the element `1`.
415 pub fn one(self: Self) Fe {
416 var fe = self.zero;
417 fe.v.limbs()[0] = 1;
418 return fe;
419 }
420
421 /// Creates a new modulus from a `Uint` value.
422 /// The modulus must be odd and larger than 2.
423 pub fn fromUint(v_: FeUint) InvalidModulusError!Self {
424 if (!v_.isOdd()) return error.EvenModulus;
425
426 var v = v_.normalize();
427 const hi = v.limbsConst()[v.limbs_len - 1];
428 const lo = v.limbsConst()[0];
429
430 if (v.limbs_len < 2 and lo < 3) {
431 return error.ModulusTooSmall;
432 }
433
434 const leading = @clz(hi) - carry_bits;
435
436 var y = lo;
437
438 inline for (0..comptime math.log2_int(usize, t_bits)) |_| {
439 y = y *% (2 -% lo *% y);
440 }
441 const m0inv = (@as(Limb, 1) << t_bits) - (@as(TLimb, @truncate(y)));
442
443 const zero = Fe{ .v = FeUint.zero };
444
445 var m = Self{
446 .zero = zero,
447 .v = v,
448 .leading = leading,
449 .m0inv = m0inv,
450 .rr = undefined, // will be computed right after
451 };
452 m.shrink(&m.zero) catch unreachable;
453 computeRR(&m);
454
455 return m;
456 }
457
458 /// Creates a new modulus from a primitive value.
459 /// The modulus must be odd and larger than 2.
460 pub fn fromPrimitive(comptime T: type, x: T) (InvalidModulusError || OverflowError)!Self {
461 comptime assert(@bitSizeOf(T) <= max_bits); // Primitive type is larger than the modulus type.
462 const v = try FeUint.fromPrimitive(T, x);
463 return try Self.fromUint(v);
464 }
465
466 /// Creates a new modulus from a byte string.
467 pub fn fromBytes(bytes: []const u8, comptime endian: Endian) (InvalidModulusError || OverflowError)!Self {
468 const v = try FeUint.fromBytes(bytes, endian);
469 return try Self.fromUint(v);
470 }
471
472 /// Serializes the modulus to a byte string.
473 pub fn toBytes(self: Self, bytes: []u8, comptime endian: Endian) OverflowError!void {
474 return self.v.toBytes(bytes, endian);
475 }
476
477 /// Rejects field elements that are not in the canonical form.
478 pub fn rejectNonCanonical(self: Self, fe: Fe) error{NonCanonical}!void {
479 if (fe.limbs_count() != self.limbs_count() or ct.limbsCmpGeq(fe.v, self.v)) {
480 return error.NonCanonical;
481 }
482 }
483
484 // Makes the number of active limbs in a field element match the one of the modulus.
485 fn shrink(self: Self, fe: *Fe) OverflowError!void {
486 const new_len = self.limbs_count();
487 if (fe.limbs_count() < new_len) return error.Overflow;
488 var acc: Limb = 0;
489 for (fe.v.limbsConst()[new_len..]) |limb| {
490 acc |= limb;
491 }
492 if (acc != 0) return error.Overflow;
493 if (new_len > fe.v.limbs_buffer.len) return error.Overflow;
494 fe.v.limbs_len = new_len;
495 }
496
497 // Computes R^2 for the Montgomery representation.
498 fn computeRR(self: *Self) void {
499 self.rr = self.zero;
500 const n = self.rr.limbs_count();
501 self.rr.v.limbs()[n - 1] = 1;
502 for ((n - 1)..(2 * n)) |_| {
503 self.shiftIn(&self.rr, 0);
504 }
505 self.shrink(&self.rr) catch unreachable;
506 }
507
508 /// Computes x << t_bits + y (mod m)
509 fn shiftIn(self: Self, x: *Fe, y: Limb) void {
510 var d = self.zero;
511 const x_limbs = x.v.limbs();
512 const d_limbs = d.v.limbs();
513 const m_limbs = self.v.limbsConst();
514
515 var need_sub = false;
516 var i: usize = t_bits - 1;
517 while (true) : (i -= 1) {
518 var carry: u1 = @truncate(math.shr(Limb, y, i));
519 var borrow: u1 = 0;
520 for (0..self.limbs_count()) |j| {
521 const l = ct.select(need_sub, d_limbs[j], x_limbs[j]);
522 var res = (l << 1) + carry;
523 x_limbs[j] = @as(TLimb, @truncate(res));
524 carry = @truncate(res >> t_bits);
525
526 res = x_limbs[j] -% m_limbs[j] -% borrow;
527 d_limbs[j] = @as(TLimb, @truncate(res));
528
529 borrow = @truncate(res >> t_bits);
530 }
531 need_sub = ct.eql(carry, borrow);
532 if (i == 0) break;
533 }
534 x.v.cmov(need_sub, d.v);
535 }
536
537 /// Adds two field elements (mod m).
538 pub fn add(self: Self, x: Fe, y: Fe) Fe {
539 var out = x;
540 if (x.montgomery == y.montgomery) {
541 @branchHint(.likely);
542 const overflow = out.v.addWithOverflow(y.v);
543 const underflow: u1 = @bitCast(ct.limbsCmpLt(out.v, self.v));
544 const need_sub = ct.eql(overflow, underflow);
545 _ = out.v.conditionalSubWithOverflow(need_sub, self.v);
546 return out;
547 } else {
548 var y_ = y;
549 if (y.montgomery) {
550 self.fromMontgomery(&y_) catch unreachable;
551 } else {
552 self.toMontgomery(&y_) catch unreachable;
553 }
554 const overflow = out.v.addWithOverflow(y_.v);
555 const underflow: u1 = @bitCast(ct.limbsCmpLt(out.v, self.v));
556 const need_sub = ct.eql(overflow, underflow);
557 _ = out.v.conditionalSubWithOverflow(need_sub, self.v);
558 return out;
559 }
560 }
561
562 /// Subtracts two field elements (mod m).
563 pub fn sub(self: Self, x: Fe, y: Fe) Fe {
564 var out = x;
565 if (x.montgomery == y.montgomery) {
566 const underflow: bool = @bitCast(out.v.subWithOverflow(y.v));
567 _ = out.v.conditionalAddWithOverflow(underflow, self.v);
568 return out;
569 } else {
570 var y_ = y;
571 if (y.montgomery) {
572 self.fromMontgomery(&y_) catch unreachable;
573 } else {
574 self.toMontgomery(&y_) catch unreachable;
575 }
576 const underflow: bool = @bitCast(out.v.subWithOverflow(y_.v));
577 _ = out.v.conditionalAddWithOverflow(underflow, self.v);
578 return out;
579 }
580 }
581
582 /// Converts a field element to the Montgomery form.
583 pub fn toMontgomery(self: Self, x: *Fe) RepresentationError!void {
584 if (x.montgomery) {
585 return error.UnexpectedRepresentation;
586 }
587 self.shrink(x) catch unreachable;
588 x.* = self.montgomeryMul(x.*, self.rr);
589 x.montgomery = true;
590 }
591
592 /// Takes a field element out of the Montgomery form.
593 pub fn fromMontgomery(self: Self, x: *Fe) RepresentationError!void {
594 if (!x.montgomery) {
595 return error.UnexpectedRepresentation;
596 }
597 self.shrink(x) catch unreachable;
598 x.* = self.montgomeryMul(x.*, self.one());
599 x.montgomery = false;
600 }
601
602 /// Reduces an arbitrary `Uint`, converting it to a field element.
603 pub fn reduce(self: Self, x: anytype) Fe {
604 var out = self.zero;
605 var i = x.limbs_len - 1;
606 if (self.limbs_count() >= 2) {
607 const start = @min(i, self.limbs_count() - 2);
608 var j = start;
609 while (true) : (j -= 1) {
610 out.v.limbs()[j] = x.limbsConst()[i];
611 i -= 1;
612 if (j == 0) break;
613 }
614 }
615 while (true) : (i -= 1) {
616 self.shiftIn(&out, x.limbsConst()[i]);
617 if (i == 0) break;
618 }
619 return out;
620 }
621
622 fn montgomeryLoop(self: Self, d: *Fe, x: Fe, y: Fe) u1 {
623 assert(d.limbs_count() == x.limbs_count());
624 assert(d.limbs_count() == y.limbs_count());
625 assert(d.limbs_count() == self.limbs_count());
626
627 const a_limbs = x.v.limbsConst();
628 const b_limbs = y.v.limbsConst();
629 const d_limbs = d.v.limbs();
630 const m_limbs = self.v.limbsConst();
631
632 var overflow: u1 = 0;
633 for (0..self.limbs_count()) |i| {
634 var carry: Limb = 0;
635
636 var wide = ct.mulWide(a_limbs[i], b_limbs[0]);
637 var z_lo = @addWithOverflow(d_limbs[0], wide.lo);
638 const f = @as(TLimb, @truncate(z_lo[0] *% self.m0inv));
639 var z_hi = wide.hi +% z_lo[1];
640 wide = ct.mulWide(f, m_limbs[0]);
641 z_lo = @addWithOverflow(z_lo[0], wide.lo);
642 z_hi +%= z_lo[1];
643 z_hi +%= wide.hi;
644 carry = (z_hi << 1) | (z_lo[0] >> t_bits);
645
646 for (1..self.limbs_count()) |j| {
647 wide = ct.mulWide(a_limbs[i], b_limbs[j]);
648 z_lo = @addWithOverflow(d_limbs[j], wide.lo);
649 z_hi = wide.hi +% z_lo[1];
650 wide = ct.mulWide(f, m_limbs[j]);
651 z_lo = @addWithOverflow(z_lo[0], wide.lo);
652 z_hi +%= z_lo[1];
653 z_hi +%= wide.hi;
654 z_lo = @addWithOverflow(z_lo[0], carry);
655 z_hi +%= z_lo[1];
656 if (j > 0) {
657 d_limbs[j - 1] = @as(TLimb, @truncate(z_lo[0]));
658 }
659 carry = (z_hi << 1) | (z_lo[0] >> t_bits);
660 }
661 const z = overflow + carry;
662 d_limbs[self.limbs_count() - 1] = @as(TLimb, @truncate(z));
663 overflow = @as(u1, @truncate(z >> t_bits));
664 }
665 return overflow;
666 }
667
668 // Montgomery multiplication.
669 fn montgomeryMul(self: Self, x: Fe, y: Fe) Fe {
670 var d = self.zero;
671 assert(x.limbs_count() == self.limbs_count());
672 assert(y.limbs_count() == self.limbs_count());
673 const overflow = self.montgomeryLoop(&d, x, y);
674 const underflow = 1 -% @intFromBool(ct.limbsCmpGeq(d.v, self.v));
675 const need_sub = ct.eql(overflow, underflow);
676 _ = d.v.conditionalSubWithOverflow(need_sub, self.v);
677 d.montgomery = x.montgomery == y.montgomery;
678 return d;
679 }
680
681 // Montgomery squaring.
682 fn montgomerySq(self: Self, x: Fe) Fe {
683 var d = self.zero;
684 assert(x.limbs_count() == self.limbs_count());
685 const overflow = self.montgomeryLoop(&d, x, x);
686 const underflow = 1 -% @intFromBool(ct.limbsCmpGeq(d.v, self.v));
687 const need_sub = ct.eql(overflow, underflow);
688 _ = d.v.conditionalSubWithOverflow(need_sub, self.v);
689 d.montgomery = true;
690 return d;
691 }
692
693 // Returns x^e (mod m), with the exponent provided as a byte string.
694 // `public` must be set to `false` if the exponent it secret.
695 fn powWithEncodedExponentInternal(self: Self, x: Fe, e: []const u8, endian: Endian, comptime public: bool) NullExponentError!Fe {
696 var acc: u8 = 0;
697 for (e) |b| acc |= b;
698 if (acc == 0) return error.NullExponent;
699
700 const was_montgomery = x.montgomery;
701
702 var out = self.one();
703 self.toMontgomery(&out) catch unreachable;
704
705 if (public and
706 (e.len < 3 or (e.len == 3 and e[if (endian == .big) 0 else 2] <= 0b1111)))
707 {
708 // Do not use a precomputation table for short, public exponents
709 var x_m = x;
710 if (!x.montgomery) {
711 self.toMontgomery(&x_m) catch unreachable;
712 }
713 var s = switch (endian) {
714 .big => 0,
715 .little => e.len - 1,
716 };
717 while (true) {
718 const b = e[s];
719 var j: u3 = 7;
720 while (true) : (j -= 1) {
721 out = self.montgomerySq(out);
722 const k: u1 = @truncate(b >> j);
723 if (k != 0) {
724 const t = self.montgomeryMul(out, x_m);
725 @memcpy(out.v.limbs(), t.v.limbsConst());
726 }
727 if (j == 0) break;
728 }
729 switch (endian) {
730 .big => {
731 s += 1;
732 if (s == e.len) break;
733 },
734 .little => {
735 if (s == 0) break;
736 s -= 1;
737 },
738 }
739 }
740 } else {
741 // Use a precomputation table for large exponents
742 var pc: [15]Fe = [1]Fe{x} ++ @as([14]Fe, @splat(self.zero));
743 if (!x.montgomery) {
744 self.toMontgomery(&pc[0]) catch unreachable;
745 }
746 for (1..pc.len) |i| {
747 pc[i] = self.montgomeryMul(pc[i - 1], pc[0]);
748 }
749 var t0 = self.zero;
750 var s = switch (endian) {
751 .big => 0,
752 .little => e.len - 1,
753 };
754 while (true) {
755 const b = e[s];
756 for ([_]u3{ 4, 0 }) |j| {
757 for (0..4) |_| {
758 out = self.montgomerySq(out);
759 }
760 const k = (b >> j) & 0b1111;
761 if (public or std.options.side_channels_mitigations == .none) {
762 if (k == 0) continue;
763 t0 = pc[k - 1];
764 } else {
765 for (pc, 0..) |t, i| {
766 t0.v.cmov(ct.eql(k, @as(u8, @truncate(i + 1))), t.v);
767 }
768 }
769 const t1 = self.montgomeryMul(out, t0);
770 if (public) {
771 @memcpy(out.v.limbs(), t1.v.limbsConst());
772 } else {
773 out.v.cmov(!ct.eql(k, 0), t1.v);
774 }
775 }
776 switch (endian) {
777 .big => {
778 s += 1;
779 if (s == e.len) break;
780 },
781 .little => {
782 if (s == 0) break;
783 s -= 1;
784 },
785 }
786 }
787 }
788 if (!was_montgomery) {
789 self.fromMontgomery(&out) catch unreachable;
790 }
791 return out;
792 }
793
794 /// Multiplies two field elements.
795 /// Result preserves the first operand's form.
796 pub fn mul(self: Self, x: Fe, y: Fe) Fe {
797 if (x.montgomery) {
798 const y_ = if (!y.montgomery) blk: {
799 var yy = y;
800 self.toMontgomery(&yy) catch unreachable;
801 break :blk yy;
802 } else y;
803 return self.montgomeryMul(x, y_);
804 } else {
805 var x_m = x;
806 var y_m = if (y.montgomery) blk: {
807 var yy = y;
808 self.fromMontgomery(&yy) catch unreachable;
809 break :blk yy;
810 } else y;
811 self.toMontgomery(&x_m) catch unreachable;
812 self.toMontgomery(&y_m) catch unreachable;
813 var out = self.montgomeryMul(x_m, y_m);
814 self.fromMontgomery(&out) catch unreachable;
815 return out;
816 }
817 }
818
819 /// Squares a field element.
820 pub fn sq(self: Self, x: Fe) Fe {
821 if (x.montgomery) {
822 return self.montgomerySq(x);
823 } else {
824 var out = x;
825 self.toMontgomery(&out) catch unreachable;
826 out = self.montgomerySq(out);
827 self.fromMontgomery(&out) catch unreachable;
828 return out;
829 }
830 }
831
832 /// Returns x^e (mod m) in constant time.
833 pub fn pow(self: Self, x: Fe, e: Fe) (NullExponentError || RepresentationError)!Fe {
834 if (e.montgomery) {
835 return error.UnexpectedRepresentation;
836 }
837 var buf: [Fe.encoded_bytes]u8 = undefined;
838 e.toBytes(&buf, native_endian) catch unreachable;
839 return self.powWithEncodedExponent(x, &buf, native_endian);
840 }
841
842 /// Returns x^e (mod m), assuming that the exponent is public.
843 /// The function remains constant time with respect to `x`.
844 pub fn powPublic(self: Self, x: Fe, e: Fe) (NullExponentError || RepresentationError)!Fe {
845 if (e.montgomery) {
846 return error.UnexpectedRepresentation;
847 }
848 var e_normalized = Fe{ .v = e.v.normalize() };
849 var buf_: [Fe.encoded_bytes]u8 = undefined;
850 var buf = buf_[0..@divCeil(e_normalized.v.limbs_len * t_bits, 8)];
851 e_normalized.toBytes(buf, .little) catch unreachable;
852 const leading = @clz(e_normalized.v.limbsConst()[e_normalized.v.limbs_len - carry_bits]);
853 buf = buf[0 .. buf.len - leading / 8];
854 return self.powWithEncodedPublicExponent(x, buf, .little);
855 }
856
857 /// Returns x^e (mod m), with the exponent provided as a byte string.
858 /// Exponents are usually small, so this function is faster than `powPublic` as a field element
859 /// doesn't have to be created if a serialized representation is already available.
860 ///
861 /// If the exponent is public, `powWithEncodedPublicExponent()` can be used instead for a slight speedup.
862 pub fn powWithEncodedExponent(self: Self, x: Fe, e: []const u8, endian: Endian) NullExponentError!Fe {
863 return self.powWithEncodedExponentInternal(x, e, endian, false);
864 }
865
866 /// Returns x^e (mod m), the exponent being public and provided as a byte string.
867 /// Exponents are usually small, so this function is faster than `powPublic` as a field element
868 /// doesn't have to be created if a serialized representation is already available.
869 ///
870 /// If the exponent is secret, `powWithEncodedExponent` must be used instead.
871 pub fn powWithEncodedPublicExponent(self: Self, x: Fe, e: []const u8, endian: Endian) NullExponentError!Fe {
872 return self.powWithEncodedExponentInternal(x, e, endian, true);
873 }
874 };
875}
876
877const ct = if (std.options.side_channels_mitigations == .none) ct_unprotected else ct_protected;
878
879const ct_protected = struct {
880 // Returns x if on is true, otherwise y.
881 fn select(on: bool, x: Limb, y: Limb) Limb {
882 const mask = @as(Limb, 0) -% @intFromBool(on);
883 return y ^ (mask & (y ^ x));
884 }
885
886 // Compares two values in constant time.
887 fn eql(x: anytype, y: @TypeOf(x)) bool {
888 const c1 = @subWithOverflow(x, y)[1];
889 const c2 = @subWithOverflow(y, x)[1];
890 return @as(bool, @bitCast(1 - (c1 | c2)));
891 }
892
893 // Compares two big integers in constant time, returning true if x < y.
894 fn limbsCmpLt(x: anytype, y: @TypeOf(x)) bool {
895 var c: u1 = 0;
896 for (x.limbsConst(), y.limbsConst()) |x_limb, y_limb| {
897 c = @truncate((x_limb -% y_limb -% c) >> t_bits);
898 }
899 return c != 0;
900 }
901
902 // Compares two big integers in constant time, returning true if x >= y.
903 fn limbsCmpGeq(x: anytype, y: @TypeOf(x)) bool {
904 return !limbsCmpLt(x, y);
905 }
906
907 // Multiplies two limbs and returns the result as a wide limb.
908 fn mulWide(x: Limb, y: Limb) WideLimb {
909 const half_bits = @typeInfo(Limb).int.bits / 2;
910 const Half = @Int(.unsigned, half_bits);
911 const x0 = @as(Half, @truncate(x));
912 const x1 = @as(Half, @truncate(x >> half_bits));
913 const y0 = @as(Half, @truncate(y));
914 const y1 = @as(Half, @truncate(y >> half_bits));
915 const w0 = math.mulWide(Half, x0, y0);
916 const t = math.mulWide(Half, x1, y0) + (w0 >> half_bits);
917 var w1: Limb = @as(Half, @truncate(t));
918 const w2 = @as(Half, @truncate(t >> half_bits));
919 w1 += math.mulWide(Half, x0, y1);
920 const hi = math.mulWide(Half, x1, y1) + w2 + (w1 >> half_bits);
921 const lo = x *% y;
922 return .{ .hi = hi, .lo = lo };
923 }
924};
925
926const ct_unprotected = struct {
927 // Returns x if on is true, otherwise y.
928 fn select(on: bool, x: Limb, y: Limb) Limb {
929 return if (on) x else y;
930 }
931
932 // Compares two values in constant time.
933 fn eql(x: anytype, y: @TypeOf(x)) bool {
934 return x == y;
935 }
936
937 // Compares two big integers in constant time, returning true if x < y.
938 fn limbsCmpLt(x: anytype, y: @TypeOf(x)) bool {
939 const x_limbs = x.limbsConst();
940 const y_limbs = y.limbsConst();
941 assert(x_limbs.len == y_limbs.len);
942
943 var i = x_limbs.len;
944 while (i != 0) {
945 i -= 1;
946 if (x_limbs[i] != y_limbs[i]) {
947 return x_limbs[i] < y_limbs[i];
948 }
949 }
950 return false;
951 }
952
953 // Compares two big integers in constant time, returning true if x >= y.
954 fn limbsCmpGeq(x: anytype, y: @TypeOf(x)) bool {
955 return !limbsCmpLt(x, y);
956 }
957
958 // Multiplies two limbs and returns the result as a wide limb.
959 fn mulWide(x: Limb, y: Limb) WideLimb {
960 const wide = math.mulWide(Limb, x, y);
961 return .{
962 .hi = @as(Limb, @truncate(wide >> @typeInfo(Limb).int.bits)),
963 .lo = @as(Limb, @truncate(wide)),
964 };
965 }
966};
967
968test "finite field arithmetic" {
969 const M = Modulus(256);
970 const m = try M.fromPrimitive(u256, 3429938563481314093726330772853735541133072814650493833233);
971 var x = try M.Fe.fromPrimitive(u256, m, 80169837251094269539116136208111827396136208141182357733);
972 var y = try M.Fe.fromPrimitive(u256, m, 24620149608466364616251608466389896540098571);
973
974 const x_ = try x.toPrimitive(u256);
975 try testing.expect((try M.Fe.fromPrimitive(@TypeOf(x_), m, x_)).eql(x));
976 try testing.expectError(error.Overflow, x.toPrimitive(u50));
977
978 const bits = m.bits();
979 try testing.expectEqual(bits, 192);
980
981 var x_y = m.mul(x, y);
982 try testing.expectEqual(x_y.toPrimitive(u256), 1666576607955767413750776202132407807424848069716933450241);
983
984 try m.toMontgomery(&x);
985 x_y = m.mul(x, y);
986 try testing.expect(x_y.montgomery); // result preserves first operand's form
987 try m.fromMontgomery(&x_y);
988 try testing.expectEqual(x_y.toPrimitive(u256), 1666576607955767413750776202132407807424848069716933450241);
989 try m.fromMontgomery(&x);
990
991 x = m.add(x, y);
992 try testing.expectEqual(x.toPrimitive(u256), 80169837251118889688724602572728079004602598037722456304);
993 x = m.sub(x, y);
994 try testing.expectEqual(x.toPrimitive(u256), 80169837251094269539116136208111827396136208141182357733);
995
996 const big = try Uint(512).fromPrimitive(u495, 77285373554113307281465049383342993856348131409372633077285373554113307281465049383323332333429938563481314093726330772853735541133072814650493833233);
997 const reduced = m.reduce(big);
998 try testing.expectEqual(reduced.toPrimitive(u495), 858047099884257670294681641776170038885500210968322054970);
999
1000 const x_pow_y = try m.powPublic(x, y);
1001 try testing.expectEqual(x_pow_y.toPrimitive(u256), 1631933139300737762906024873185789093007782131928298618473);
1002 try testing.expect(!x_pow_y.montgomery);
1003 try m.toMontgomery(&x);
1004 var x_pow_y2 = try m.powPublic(x, y);
1005 try testing.expect(x_pow_y2.montgomery);
1006 try m.fromMontgomery(&x_pow_y2);
1007 try m.fromMontgomery(&x);
1008 try testing.expect(x_pow_y2.eql(x_pow_y));
1009 try testing.expectError(error.NullExponent, m.powPublic(x, m.zero));
1010
1011 try testing.expect(!x.isZero());
1012 try testing.expect(!y.isZero());
1013 try testing.expect(m.v.isOdd());
1014
1015 const x_sq = m.sq(x);
1016 const x_sq2 = m.mul(x, x);
1017 try testing.expect(!x_sq.montgomery);
1018 try testing.expect(!x_sq2.montgomery);
1019 try testing.expect(x_sq.eql(x_sq2));
1020 try m.toMontgomery(&x);
1021 var x_sq3 = m.sq(x);
1022 var x_sq4 = m.mul(x, x);
1023 try testing.expect(x_sq3.montgomery);
1024 try testing.expect(x_sq4.montgomery);
1025 try m.fromMontgomery(&x_sq3);
1026 try m.fromMontgomery(&x_sq4);
1027 try testing.expect(x_sq.eql(x_sq3));
1028 try testing.expect(x_sq3.eql(x_sq4));
1029 try m.fromMontgomery(&x);
1030
1031 var x_mont = x;
1032 try m.toMontgomery(&x_mont);
1033
1034 // Non-montgomery + montgomery
1035 const add_nm_m = m.add(x, x_mont);
1036 try testing.expect(!add_nm_m.montgomery);
1037 var add_m_nm = m.add(x_mont, x);
1038 try testing.expect(add_m_nm.montgomery);
1039 try m.fromMontgomery(&add_m_nm);
1040 try testing.expect(add_nm_m.eql(add_m_nm));
1041
1042 // Non-montgomery - montgomery
1043 const sub_nm_m = m.sub(x, y);
1044 try testing.expect(!sub_nm_m.montgomery);
1045 var y_mont = y;
1046 try m.toMontgomery(&y_mont);
1047 var sub_m_nm = m.sub(x_mont, y);
1048 try testing.expect(sub_m_nm.montgomery);
1049 try m.fromMontgomery(&sub_m_nm);
1050 try testing.expect(sub_nm_m.eql(sub_m_nm));
1051
1052 // mul: preserves first operand's form
1053 const mul_nm_m = m.mul(x, x_mont);
1054 try testing.expect(!mul_nm_m.montgomery);
1055 const mul_nm_nm = m.mul(x, x);
1056 try testing.expect(mul_nm_m.eql(mul_nm_nm));
1057 var mul_m_nm = m.mul(x_mont, x);
1058 try testing.expect(mul_m_nm.montgomery);
1059 try m.fromMontgomery(&mul_m_nm);
1060 try testing.expect(mul_m_nm.eql(mul_nm_nm));
1061
1062 try testing.expectEqual(x.toPrimitive(u256), 80169837251094269539116136208111827396136208141182357733);
1063 try testing.expectError(error.UnexpectedRepresentation, x_mont.toPrimitive(u256));
1064}
1065
1066fn testCt(ct_: anytype) !void {
1067 const l0: Limb = 0;
1068 const l1: Limb = 1;
1069 try testing.expectEqual(l1, ct_.select(true, l1, l0));
1070 try testing.expectEqual(l0, ct_.select(false, l1, l0));
1071 try testing.expectEqual(false, ct_.eql(l1, l0));
1072 try testing.expectEqual(true, ct_.eql(l1, l1));
1073
1074 const M = Modulus(256);
1075 const m = try M.fromPrimitive(u256, 3429938563481314093726330772853735541133072814650493833233);
1076 const x = try M.Fe.fromPrimitive(u256, m, 80169837251094269539116136208111827396136208141182357733);
1077 const y = try M.Fe.fromPrimitive(u256, m, 24620149608466364616251608466389896540098571);
1078 try testing.expectEqual(false, ct_.limbsCmpLt(x.v, y.v));
1079 try testing.expectEqual(true, ct_.limbsCmpGeq(x.v, y.v));
1080
1081 try testing.expectEqual(WideLimb{ .hi = 0, .lo = 0x88 }, ct_.mulWide(1 << 3, (1 << 4) + 1));
1082}
1083
1084test ct {
1085 try testCt(ct_protected);
1086 try testCt(ct_unprotected);
1087}