saymrwulf/betrusted-curve25519-dalek-source
1
0
Fork 0
mirror of https://github.com/saymrwulf/betrusted-curve25519-dalek-source.git synced 2026-09-05 20:30:54 +00:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge branch 'release/0.12.0'

Browse source
This commit is contained in:
Isis Lovecruft 2017-10-05 06:31:58 +00:00
commit a6aaa6fae6
Failed to extract signature
15 changed files with 1729 additions and 405 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

4
Cargo.toml
View file

@ -1,6 +1,6 @@
[package]
name = "curve25519-dalek"
version = "0.11.0"
version = "0.12.0"
authors = ["Isis Lovecruft <isis@patternsinthevoid.net>",
"Henry de Valence <hdevalence@hdevalence.ca>"]
readme = "README.md"
@ -33,7 +33,7 @@ version = "0.3"
version = "0.6"
[dependencies.subtle]
version = "^0.2"
version = "^0.3"
default-features = false
[dependencies.generic-array]

2
README.md
View file

@ -44,7 +44,7 @@ Extensive documentation is available [here](https://docs.rs/curve25519-dalek).
To install, add the following to the dependencies section of your project's
`Cargo.toml`:
curve25519-dalek = "^0.11"
curve25519-dalek = "^0.12"
Then, in your library or executable source, add:

9
src/constants.rs
View file

@ -21,6 +21,7 @@
use edwards::CompressedEdwardsY;
#[cfg(feature = "yolocrypto")]
use decaf::{DecafPoint, DecafBasepointTable};
use montgomery::CompressedMontgomeryU;
use scalar::Scalar;
#[cfg(feature="radix_51")]
@ -52,6 +53,14 @@ pub const BASE_CMPRSSD: CompressedEdwardsY =
0x66, 0x66, 0x66, 0x66, 0x66, 0x66, 0x66, 0x66,
0x66, 0x66, 0x66, 0x66, 0x66, 0x66, 0x66, 0x66]);
/// The X25519 basepoint, in compressed Montgomery form.
pub const BASE_COMPRESSED_MONTGOMERY: CompressedMontgomeryU =
CompressedMontgomeryU([0x09, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00]);
/// The Ed25519 basepoint, as a `DecafPoint`. This is called `_POINT` to distinguish it from
/// `_TABLE`, which provides fast scalar multiplication.
#[cfg(feature = "yolocrypto")] pub const DECAF_ED25519_BASEPOINT_POINT: DecafPoint =

24
src/constants_32bit.rs
View file

@ -19,6 +19,7 @@
#![allow(non_snake_case)]
use field_32bit::FieldElement32;
use scalar_32bit::Scalar32;
use edwards::ExtendedPoint;
use edwards::AffineNielsPoint;
use edwards::EdwardsBasepointTable;
@ -75,6 +76,9 @@ pub const HALF: FieldElement32 = FieldElement32([
pub const A: FieldElement32 = FieldElement32([
486662, 0, 0, 0, 0, 0, 0, 0, 0, 0, ]);
/// `APLUS2_OVER_FOUR` is (A+2)/4. (This is used internally within the Montgomery ladder.)
pub const APLUS2_OVER_FOUR: FieldElement32 = FieldElement32([121666, 0, 0, 0, 0, 0, 0, 0, 0, 0]);
/// `SQRT_MINUS_A` is sqrt(-486662)
// XXX I think that this was used in Adam's code for his elligator
// implementation, but that should maybe be using sqrt(-486664)
@ -93,6 +97,26 @@ pub const SQRT_MINUS_HALF: FieldElement32 = FieldElement32([ // sqrtMinusHalf
-17256545, 3971863, 28865457, -1750208, 27359696,
-16640980, 12573105, 1002827, -163343, 11073975, ]);
/// `L` is the order of base point, i.e. 2^252 +
/// 27742317777372353535851937790883648493
pub const L: Scalar32 = Scalar32([ 0x1cf5d3ed, 0x009318d2, 0x1de73596, 0x1df3bd45,
0x0000014d, 0x00000000, 0x00000000, 0x00000000,
0x00100000 ]);
/// `L` * `LFACTOR` = -1 (mod 2^29)
pub const LFACTOR: u32 = 0x12547e1b;
/// `R` = R % L where R = 2^261
pub const R: Scalar32 = Scalar32([ 0x114df9ed, 0x1a617303, 0x0f7c098c, 0x16793167,
0x1ffd656e, 0x1fffffff, 0x1fffffff, 0x1fffffff,
0x000fffff ]);
/// `RR` = (R^2) % L where R = 2^261
pub const RR: Scalar32 = Scalar32([ 0x0b5f9d12, 0x1e141b17, 0x158d7f3d, 0x143f3757,
0x1972d781, 0x042feb7c, 0x1ceec73d, 0x1e184d1e,
0x0005046d ]);
/// Basepoint has y = 4/5. This is called `_POINT` to distinguish it from `_TABLE`, which should
/// be used for scalar multiplication (it's much faster).
pub const ED25519_BASEPOINT_POINT: ExtendedPoint = ExtendedPoint{

16
src/constants_64bit.rs
View file

@ -19,6 +19,7 @@
#![allow(non_snake_case)]
use field_64bit::FieldElement64;
use scalar_64bit::Scalar64;
use edwards::ExtendedPoint;
use edwards::AffineNielsPoint;
use edwards::EdwardsBasepointTable;
@ -54,6 +55,9 @@ pub const HALF: FieldElement64 = FieldElement64([2251799813685239, 2251799813685
/// In Montgomery form y² = x³+Ax²+x, Curve25519 has A=486662.
pub const A: FieldElement64 = FieldElement64([486662, 0, 0, 0, 0]);
/// `APLUS2_OVER_FOUR` is (A+2)/4. (This is used internally within the Montgomery ladder.)
pub const APLUS2_OVER_FOUR: FieldElement64 = FieldElement64([121666, 0, 0, 0, 0]);
/// `SQRT_MINUS_A` is sqrt(-486662)
// XXX I think that this was used in Adam's code for his elligator
// implementation, but that should maybe be using sqrt(-486664)
@ -66,6 +70,18 @@ pub const SQRT_MINUS_APLUS2: FieldElement64 = FieldElement64([1693982333959686,
/// `SQRT_MINUS_HALF` is sqrt(-1/2)
pub const SQRT_MINUS_HALF: FieldElement64 = FieldElement64([266547196637087, 2134345371906993, 1135042577398223, 67298593331632, 743161882051057]);
/// `L` is the order of base point, i.e. 2^252 + 27742317777372353535851937790883648493
pub const L: Scalar64 = Scalar64([ 0x0002631a5cf5d3ed, 0x000dea2f79cd6581, 0x000000000014def9, 0x0000000000000000, 0x0000100000000000 ]);
/// `L` * `LFACTOR` = -1 (mod 2^51)
pub const LFACTOR: u64 = 0x51da312547e1b;
/// `R` = R % L where R = 2^260
pub const R: Scalar64 = Scalar64([ 0x000f48bd6721e6ed, 0x0003bab5ac67e45a, 0x000fffffeb35e51b, 0x000fffffffffffff, 0x00000fffffffffff ]);
/// `RR` = (R^2) % L where R = 2^260
pub const RR: Scalar64 = Scalar64([ 0x0009d265e952d13b, 0x000d63c715bea69f, 0x0005be65cb687604, 0x0003dceec73d217f, 0x000009411b7c309a ]);
/// Basepoint has y = 4/5. This is called `_POINT` to distinguish it from `_TABLE`, which should
/// be used for scalar multiplication (it's much faster).
pub const ED25519_BASEPOINT_POINT: ExtendedPoint = ExtendedPoint{

6
src/decaf.rs
View file

@ -189,7 +189,7 @@ impl<'de> Deserialize<'de> for DecafPoint {
/// A point in a prime-order group.
///
/// XXX think about how this API should work
// XXX think about how this API should work
#[derive(Copy, Clone)]
pub struct DecafPoint(pub ExtendedPoint);
@ -752,7 +752,7 @@ mod test {
fn decaf_decompress_id() {
let compressed_id = CompressedDecaf::identity();
let id = compressed_id.decompress().unwrap();
assert_eq!(id.0.compress_edwards(), CompressedEdwardsY::identity());
assert_eq!(id.0.compress(), CompressedEdwardsY::identity());
}
#[test]
@ -768,7 +768,7 @@ mod test {
// Check that bp_recaf differs from bp by a point of order 4
let diff = &constants::ED25519_BASEPOINT_POINT - &bp_recaf;
let diff4 = diff.mult_by_pow_2(4); // XXX this is wrong
assert_eq!(diff4.compress_edwards(), CompressedEdwardsY::identity());
assert_eq!(diff4.compress(), CompressedEdwardsY::identity());
}
#[test]

157
src/edwards.rs
View file

@ -89,7 +89,7 @@ use core::ops::Index;
use constants;
use field::FieldElement;
use scalar::Scalar;
use montgomery::CompressedMontgomeryU;
use montgomery::MontgomeryPoint;
use subtle::slices_equal;
use subtle::bytes_equal;
@ -123,7 +123,6 @@ impl CompressedEdwardsY {
}
/// Copy this `CompressedEdwardsY` to an array of bytes.
/// XXX is this useful?
pub fn to_bytes(&self) -> [u8; 32] {
self.0
}
@ -169,7 +168,7 @@ impl Serialize for ExtendedPoint {
fn serialize<S>(&self, serializer: S) -> Result<S::Ok, S::Error>
where S: Serializer
{
serializer.serialize_bytes(self.compress_edwards().as_bytes())
serializer.serialize_bytes(self.compress().as_bytes())
}
}
@ -393,8 +392,8 @@ impl ConditionallyAssignable for ExtendedPoint {
impl Equal for ExtendedPoint {
fn ct_eq(&self, other: &ExtendedPoint) -> u8 {
slices_equal(self.compress_edwards().as_bytes(),
other.compress_edwards().as_bytes())
slices_equal(self.compress().as_bytes(),
other.compress().as_bytes())
}
}
@ -436,7 +435,7 @@ impl ProjectivePoint {
}
/// Convert this point to a `CompressedEdwardsY`
pub fn compress_edwards(&self) -> CompressedEdwardsY {
pub fn compress(&self) -> CompressedEdwardsY {
let recip = self.Z.invert();
let x = &self.X * &recip;
let y = &self.Y * &recip;
@ -447,32 +446,67 @@ impl ProjectivePoint {
CompressedEdwardsY(s)
}
/// Convert this point to a `CompressedMontgomeryU`.
/// Note that this discards the sign.
/// Convert this projective point in the Edwards model to its equivalent
/// projective point on the Montgomery form of the curve.
///
/// # Return
/// - `None` if `self` is the identity point;
/// - `Some(CompressedMontgomeryU)` otherwise.
/// Taking the Montgomery curve equation in affine coordinates:
///
pub fn compress_montgomery(&self) -> Option<CompressedMontgomeryU> {
// u = (1 + y) / (1 - y)
// v = sqrt(-486664) * u / x
//
// since y = Y/Z, x = X/Z,
//
// u = (1 + Y/Z) / (1 - Y/Z);
// = (Z + Y) / (Z - Y);
//
// exceptional points:
// y = 1 <=> Y/Z = 1 <=> Z - Y = 0
let Z_plus_Y = &self.Z + &self.Y;
let Z_minus_Y = &self.Z - &self.Y;
let u = &Z_plus_Y * &Z_minus_Y.invert();
if Z_minus_Y.is_zero() == 0u8 {
Some(CompressedMontgomeryU(u.to_bytes()))
} else {
None
///     E_(A,B) = Bv² = u³ + Au² + u   <span style="float: right">(1)</span>
///
/// and given its relations to the coordinates of the Edwards model:
///
///     u = (1+y)/(1-y)        <span style="float: right">(2)</span>
///     v = (λu)/(x)
///
/// Converting from affine to projective coordinates in the Montgomery
/// model, we arrive at:
///
///     u = (Z+Y)/(Z-Y)        <span style="float: right">(3)</span>
///     v = λ * ((Z+Y)/(Z-Y)) * (Z/X)
///
/// The transition between affine and projective is given by
///
///     u → U/W        <span style="float: right">(4)</span>
///     v → V/W
///
/// thus the Montgomery curve equation (1) becomes
///
///     E_(A,B) : BV²W = U³ + AU²W + UW² ⊆ 𝗣^2  <span style="float: right">(5)</span>
///
/// Here, again, to differentiate from points in the twisted Edwards model, we
/// call the point `(x,y)` in affine coordinates `(u,v)` and similarly in projective
/// space we use `(U:V:W)`. However, since (as per Montgomery's original work) the
/// v-coordinate is superfluous to the definition of the group law, we merely
/// use `(U:W)`.
///
/// Therefore, the direct translation between projective Montgomery points
/// and projective twisted Edwards points is
///
///     (U:W) = (Z+Y:Z-Y) <span style="float: right">(6)</span>
///
/// Note, however, that there appears to be an exception where `Z=Y`,
/// since—from equation 2—this would imply that `y=1` (thus causing the
/// denominator to be zero). If this is the case, then it follows from the
/// twisted Edwards curve equation
///
///     -x² + y² = 1 + dx²y² <span style="float: right">(7)</span>
///
/// that
///
///     -x² + 1 = 1 + dx²
///
/// and, assuming that `d ≠ -1`,
///
///     -x² = x²
/// x = 0
///
/// Therefore, the only valid point with `y=1` is the twisted Edwards
/// identity point, which correctly becomes `(1:0)`, that is, the identity,
/// in the Montgomery model.
pub fn to_montgomery(&self) -> MontgomeryPoint {
MontgomeryPoint{
U: &self.Z + &self.Y,
W: &self.Z - &self.Y,
}
}
}
@ -515,20 +549,15 @@ impl ExtendedPoint {
}
}
/// Compress this point to `CompressedEdwardsY` format.
pub fn compress_edwards(&self) -> CompressedEdwardsY {
self.to_projective().compress_edwards()
/// Convert this point to its equivalent on the Montgomery form of the
/// curve.
pub fn to_montgomery(&self) -> MontgomeryPoint {
self.to_projective().to_montgomery()
}
/// Convert this point to a `CompressedMontgomeryU`.
/// Note that this discards the sign.
///
/// # Return
/// - `None` if `self` is the identity point;
/// - `Some(CompressedMontgomeryU)` otherwise.
///
pub fn compress_montgomery(&self) -> Option<CompressedMontgomeryU> {
self.to_projective().compress_montgomery()
/// Compress this point to `CompressedEdwardsY` format.
pub fn compress(&self) -> CompressedEdwardsY {
self.to_projective().compress()
}
}
@ -1305,7 +1334,7 @@ mod test {
assert!(bp.is_valid());
// Check that decompression actually gives the correct X coordinate
assert_eq!(base_X, bp.X);
assert_eq!(bp.compress_edwards(), constants::BASE_CMPRSSD);
assert_eq!(bp.compress(), constants::BASE_CMPRSSD);
}
/// Test sign handling in decompression
@ -1328,7 +1357,7 @@ mod test {
#[test]
fn basepoint_mult_one_vs_basepoint() {
let bp = &constants::ED25519_BASEPOINT_TABLE * &Scalar::one();
let compressed = bp.compress_edwards();
let compressed = bp.compress();
assert_eq!(compressed, constants::BASE_CMPRSSD);
}
@ -1336,7 +1365,7 @@ mod test {
#[test]
fn basepoint_table_basepoint_function_correct() {
let bp = constants::ED25519_BASEPOINT_TABLE.basepoint();
assert_eq!(bp.compress_edwards(), constants::BASE_CMPRSSD);
assert_eq!(bp.compress(), constants::BASE_CMPRSSD);
}
/// Test `impl Add<ExtendedPoint> for ExtendedPoint`
@ -1345,7 +1374,7 @@ mod test {
fn basepoint_plus_basepoint_vs_basepoint2() {
let bp = constants::ED25519_BASEPOINT_POINT;
let bp_added = &bp + &bp;
assert_eq!(bp_added.compress_edwards(), BASE2_CMPRSSD);
assert_eq!(bp_added.compress(), BASE2_CMPRSSD);
}
/// Test `impl Add<ProjectiveNielsPoint> for ExtendedPoint`
@ -1354,7 +1383,7 @@ mod test {
fn basepoint_plus_basepoint_projective_niels_vs_basepoint2() {
let bp = constants::ED25519_BASEPOINT_POINT;
let bp_added = (&bp + &bp.to_projective_niels()).to_extended();
assert_eq!(bp_added.compress_edwards(), BASE2_CMPRSSD);
assert_eq!(bp_added.compress(), BASE2_CMPRSSD);
}
/// Test `impl Add<AffineNielsPoint> for ExtendedPoint`
@ -1364,7 +1393,7 @@ mod test {
let bp = constants::ED25519_BASEPOINT_POINT;
let bp_affine_niels = bp.to_affine_niels();
let bp_added = (&bp + &bp_affine_niels).to_extended();
assert_eq!(bp_added.compress_edwards(), BASE2_CMPRSSD);
assert_eq!(bp_added.compress(), BASE2_CMPRSSD);
}
/// Check that equality of `ExtendedPoints` handles projective
@ -1389,15 +1418,15 @@ mod test {
let aB = &constants::ED25519_BASEPOINT_TABLE * &A_SCALAR;
let aB_affine_niels = aB.to_affine_niels();
let also_aB = (&ExtendedPoint::identity() + &aB_affine_niels).to_extended();
assert_eq!( aB.compress_edwards(),
also_aB.compress_edwards());
assert_eq!( aB.compress(),
also_aB.compress());
}
/// Test basepoint_mult versus a known scalar multiple from ed25519.py
#[test]
fn basepoint_mult_vs_ed25519py() {
let aB = &constants::ED25519_BASEPOINT_TABLE * &A_SCALAR;
assert_eq!(aB.compress_edwards(), A_TIMES_BASEPOINT);
assert_eq!(aB.compress(), A_TIMES_BASEPOINT);
}
/// Test that multiplication by the basepoint order kills the basepoint
@ -1415,20 +1444,20 @@ mod test {
let table = EdwardsBasepointTable::create(&constants::ED25519_BASEPOINT_POINT);
let aB_1 = &constants::ED25519_BASEPOINT_TABLE * &A_SCALAR;
let aB_2 = &table * &A_SCALAR;
assert_eq!(aB_1.compress_edwards(), aB_2.compress_edwards());
assert_eq!(aB_1.compress(), aB_2.compress());
}
/// Test scalar_mult versus a known scalar multiple from ed25519.py
#[test]
fn scalar_mult_vs_ed25519py() {
let aB = &constants::ED25519_BASEPOINT_POINT * &A_SCALAR;
assert_eq!(aB.compress_edwards(), A_TIMES_BASEPOINT);
assert_eq!(aB.compress(), A_TIMES_BASEPOINT);
}
/// Test basepoint.double() versus the 2*basepoint constant.
#[test]
fn basepoint_double_vs_basepoint2() {
assert_eq!(constants::ED25519_BASEPOINT_POINT.double().compress_edwards(),
assert_eq!(constants::ED25519_BASEPOINT_POINT.double().compress(),
BASE2_CMPRSSD);
}
@ -1437,14 +1466,14 @@ mod test {
fn basepoint_mult_two_vs_basepoint2() {
let mut two_bytes = [0u8; 32]; two_bytes[0] = 2;
let bp2 = &constants::ED25519_BASEPOINT_TABLE * &Scalar(two_bytes);
assert_eq!(bp2.compress_edwards(), BASE2_CMPRSSD);
assert_eq!(bp2.compress(), BASE2_CMPRSSD);
}
/// Check that converting to projective and then back to extended round-trips.
#[test]
fn basepoint_projective_extended_round_trip() {
assert_eq!(constants::ED25519_BASEPOINT_POINT
.to_projective().to_extended().compress_edwards(),
.to_projective().to_extended().compress(),
constants::BASE_CMPRSSD);
}
@ -1452,7 +1481,7 @@ mod test {
#[test]
fn basepoint16_vs_mult_by_pow_2_4() {
let bp16 = constants::ED25519_BASEPOINT_POINT.mult_by_pow_2(4);
assert_eq!(bp16.compress_edwards(), BASE16_CMPRSSD);
assert_eq!(bp16.compress(), BASE16_CMPRSSD);
}
/// Test that the conditional assignment trait works for AffineNielsPoints.
@ -1480,7 +1509,7 @@ mod test {
#[test]
fn compressed_identity() {
assert_eq!(ExtendedPoint::identity().compress_edwards(),
assert_eq!(ExtendedPoint::identity().compress(),
CompressedEdwardsY::identity());
}
@ -1518,7 +1547,7 @@ mod test {
let P1 = &G * &s;
let P2 = &s * &G;
assert!(P1.compress_edwards().to_bytes() == P2.compress_edwards().to_bytes());
assert!(P1.compress().to_bytes() == P2.compress().to_bytes());
}
#[test]
@ -1542,7 +1571,7 @@ mod test {
fn double_scalar_mult_basepoint_vs_ed25519py() {
let A = A_TIMES_BASEPOINT.decompress().unwrap();
let result = vartime::double_scalar_mult_basepoint(&A_SCALAR, &A, &B_SCALAR);
assert_eq!(result.compress_edwards(), DOUBLE_SCALAR_MULT_RESULT);
assert_eq!(result.compress(), DOUBLE_SCALAR_MULT_RESULT);
}
#[test]
@ -1552,7 +1581,7 @@ mod test {
&[A_SCALAR, B_SCALAR],
&[A, constants::ED25519_BASEPOINT_POINT]
);
assert_eq!(result.compress_edwards(), DOUBLE_SCALAR_MULT_RESULT);
assert_eq!(result.compress(), DOUBLE_SCALAR_MULT_RESULT);
}
#[test]
@ -1567,7 +1596,7 @@ mod test {
&[A, constants::ED25519_BASEPOINT_POINT]
);
assert_eq!(result_vartime.compress_edwards(), result_consttime.compress_edwards());
assert_eq!(result_vartime.compress(), result_consttime.compress());
}
}
@ -1579,7 +1608,7 @@ mod test {
fn serde_cbor_basepoint_roundtrip() {
let output = serde_cbor::to_vec(&constants::ED25519_BASEPOINT_POINT).unwrap();
let parsed: ExtendedPoint = serde_cbor::from_slice(&output).unwrap();
assert_eq!(parsed.compress_edwards(), constants::BASE_CMPRSSD);
assert_eq!(parsed.compress(), constants::BASE_CMPRSSD);
}
#[test]
@ -1615,7 +1644,7 @@ mod bench {
#[bench]
fn edwards_compress(b: &mut Bencher) {
let B = &constants::ED25519_BASEPOINT_POINT;
b.iter(|| B.compress_edwards());
b.iter(|| B.compress());
}
#[bench]

8
src/field.rs
View file

@ -197,11 +197,15 @@ impl FieldElement {
}
/// Given a nonzero field element, compute its inverse.
///
/// The inverse is computed as self^(p-2), since
/// x^(p-2)x = x^(p-1) = 1 (mod p).
///
/// XXX should we add a debug_assert that self is nonzero?
//
// XXX do we want the debug assertion to check for zero? it breaks behaviour
// such as that such as in curve25519_dalek::montgomery::test::identity_to_monty.
pub fn invert(&self) -> FieldElement {
// debug_assert!(*self != FieldElement::zero());
// The bits of p-2 = 2^255 -19 -2 are 11010111111...11.
//
// nonzero bits of exponent

2
src/field_32bit.rs
View file

@ -56,7 +56,7 @@ use utils::{load3, load4};
/// faster. However, the `FieldElement64` implementation requires Rust's
/// `u128`, which is not yet stable.
#[derive(Copy, Clone)]
pub struct FieldElement32(pub [i32; 10]);
pub struct FieldElement32(pub (crate) [i32; 10]);
impl Debug for FieldElement32 {
fn fmt(&self, f: &mut ::core::fmt::Formatter) -> ::core::fmt::Result {

2
src/field_64bit.rs
View file

@ -51,7 +51,7 @@ pub type Limb = u64;
/// which gives even better performance. This implementation requires Rust's
/// `u128`, which is not yet stable.
#[derive(Copy, Clone)]
pub struct FieldElement64(pub [u64; 5]);
pub struct FieldElement64(pub (crate) [u64; 5]);
impl Debug for FieldElement64 {
fn fmt(&self, f: &mut ::core::fmt::Formatter) -> ::core::fmt::Result {

7
src/lib.rs
View file

@ -41,6 +41,8 @@ extern crate test;
#[cfg(test)]
extern crate sha2;
// this appears to only be used for serde support right now?
#[cfg(feature = "serde")]
#[macro_use]
extern crate arrayref;
@ -71,6 +73,11 @@ mod field_32bit;
mod field_64bit;
pub mod scalar;
#[cfg(not(feature="radix_51"))]
mod scalar_32bit;
#[cfg(feature="radix_51")]
mod scalar_64bit;
pub mod edwards;
pub mod montgomery;

434
src/montgomery.rs
View file

@ -8,7 +8,19 @@
// - Isis Agora Lovecruft <isis@patternsinthevoid.net>
// - Henry de Valence <hdevalence@hdevalence.ca>
//! Montgomery arithmetic prototype, subject to revision.
//! Montgomery arithmetic.
//!
//! Apart from the compressed point implementation
//! (i.e. `CompressedMontgomeryU`), this module is a "clean room" implementation
//! of the Montgomery arithmetic described in the following papers:
//!
//! * Costello, Craig, and Benjamin Smith. "Montgomery curves and their
//! arithmetic." Journal of Cryptographic Engineering (2017): 1-14.
//! [PDF](http://eprint.iacr.org/2017/212.pdf)
//!
//! * Montgomery, Peter L. "Speeding the Pollard and elliptic curve methods of
//! factorization." Mathematics of computation 48.177 (1987): 243-264.
//! [PDF](http://www.ams.org/mcom/1987-48-177/S0025-5718-1987-0866113-7/)
// We allow non snake_case names because coordinates in projective space are
// traditionally denoted by the capitalisation of their respective
@ -16,12 +28,22 @@
// affine and projective cakes and eat both of them too.
#![allow(non_snake_case)]
use core::ops::{Mul, MulAssign};
use constants;
use field::FieldElement;
use edwards::{ExtendedPoint, CompressedEdwardsY};
use scalar::Scalar;
// XXX Move these to a common "group" module? At the same time, we should
// XXX probably make a `trait Group` once const generics are implemented in
// XXX Rust. —isis
use edwards::{Identity, ValidityCheck};
use subtle::ConditionallyAssignable;
use subtle::ConditionallySwappable;
use subtle::Equal;
use subtle::Mask;
/// In "Montgomery u" format, as used in X25519, a point `(u,v)` on
/// the Montgomery curve
@ -33,13 +55,16 @@ use subtle::ConditionallyAssignable;
/// coordinates. For Montgomery curves, it is possible to compute the
/// `u`-coordinate of `n(u,v)` just from `n` and `u`, so it is not
/// necessary to use `v` for a Diffie-Hellman key exchange.
///
/// XXX add note on monty, twist security, edwards impl of x25519, rfc7748
#[derive(Copy, Clone, Debug, PartialEq, Eq)]
pub struct CompressedMontgomeryU(pub [u8; 32]);
impl CompressedMontgomeryU {
/// View this `CompressedMontgomeryU` as an array of bytes.
pub fn as_bytes<'a>(&'a self) -> &'a [u8; 32] {
&self.0
}
/// Convert this `CompressedMontgomeryU` to an array of bytes.
pub fn to_bytes(&self) -> [u8; 32] {
self.0
}
@ -65,7 +90,7 @@ impl CompressedMontgomeryU {
/// * `v` is not square.
//
// XXX any other exceptional points for the birational map?
pub fn decompress(&self) -> Option<ExtendedPoint> {
pub fn decompress_edwards(&self) -> Option<ExtendedPoint> {
let u: FieldElement = FieldElement::from_bytes(&self.0);
// If u = -1, then v^2 = u*(u^2+486662*u+1) = 486660.
@ -84,6 +109,22 @@ impl CompressedMontgomeryU {
CompressedEdwardsY(y.to_bytes()).decompress()
}
/// Decompress this `CompressedMontgomeryU` to a `MontgomeryPoint`.
///
/// Going from affine to projective coordinates, we have:
///
///     u → U/W
///
/// # Returns
///
/// A projective `MontgomeryPoint` corresponding to this compressed point.
pub fn decompress(&self) -> MontgomeryPoint {
MontgomeryPoint{
U: FieldElement::from_bytes(&self.0),
W: FieldElement::one(),
}
}
/// Given a Montgomery `u` coordinate, compute an Edwards `y` via
/// `y = (u-1)/(u+1)`.
///
@ -150,34 +191,259 @@ impl CompressedMontgomeryU {
}
}
/// A point on the Montgomery form of the curve, in projective 𝗣^2 coordinates.
///
/// The transition between affine and projective is given by
///
///     u → U/W
///     v → V/W
///
/// thus the Montgomery curve equation
///
///     E_(A,B) : Bv² = u(u² + Au + 1)
///
/// becomes
///
///     E_(A,B) : BV²W = U(U² + AUW + W²) ⊆ 𝗣^2
///
/// Here, again, to differentiate from points in the twisted Edwards model, we
/// call the point `(x,y)` in affine coordinates `(u,v)` and similarly in projective
/// space we use `(U:V:W)`. However, since (as per Montgomery's original work) the
/// v-coordinate is superfluous for the purposes of scalar multiplication, we merely
/// use `(U:W)`.
#[derive(Copy, Clone, Debug)]
#[allow(missing_docs)]
pub struct MontgomeryPoint{
pub U: FieldElement,
pub W: FieldElement,
}
/// The identity point is a unique point (the only where `W = 0`) on the curve.
///
/// In projective coordinates, the quotient map `x : E (A,B) → E/<⦵> = 𝗣¹` is
///
///     ⎧ (x_P:1) if P = (x_P:y_P:1) ,
///     x : P ↦ ⎨
///     ⎩ (1:0) if P = O = (0:1:0) .
///
/// We emphasize that the formula `x((U: V : W)) = (U : W)` only holds on the
/// open subset of `E_(A,B)` where `W ≠ 0`; it does not extend to the point
/// `O = (0:1:0)` at infinity, because `(0:0)` is not a projective point.
///
/// # Returns
///
/// The (exceptional) point at infinity in the Montgomery model.
impl Identity for MontgomeryPoint {
fn identity() -> MontgomeryPoint {
MontgomeryPoint {
U: FieldElement::one(),
W: FieldElement::zero(),
}
}
}
/// Determine if two `MontgomeryPoint`s are equal, in constant time.
///
/// # Note
///
/// Because a compressed point on the Montgomery form of the curve doesn't
/// include the sign bit, there's two points here (if translated from the
/// Edwards form) which will equate.
///
/// # Returns
///
/// `1` if the points are equal, and `0` otherwise.
impl Equal for MontgomeryPoint {
fn ct_eq(&self, that: &MontgomeryPoint) -> u8 {
// (U_P:W_P) = (U_Q:W_Q) iff U_P * W_Q == U_Q * W_P,
// since U_P/W_P == U_Q/W_Q.
(&self.U * &that.W).ct_eq(&(&self.W * &that.U))
}
}
/// Determine if this `MontgomeryPoint` is valid.
///
/// # Note
///
/// All projective points, except for `(X:W) = (0:0)`, are valid, since the
/// projective model is linear through the origin and is comprised by all `X` in
/// ℤ/(2²⁵⁵-19), thus `(0:0)` is the only element in Fₚ² which is not a
/// projective point.
///
/// # Returns
///
/// `true` if it is valid, and `false` otherwise.
impl ValidityCheck for MontgomeryPoint {
fn is_valid(&self) -> bool {
let zero = FieldElement::zero();
if (self.U.ct_eq(&zero) & self.W.ct_eq(&zero)) == 1 {
return true;
}
false
}
}
/// Conditionally assign another `MontgomeryPoint` to this point, in constant time.
///
/// If `choice == 1`, assign `that` to `self`. Otherwise, leave `self`
/// unchanged.
impl ConditionallyAssignable for MontgomeryPoint {
fn conditional_assign(&mut self, that: &MontgomeryPoint, choice: Mask) {
self.U.conditional_assign(&that.U, choice);
self.W.conditional_assign(&that.W, choice);
}
}
impl MontgomeryPoint {
/// Compress this point to only its u-coordinate (note: affine).
///
/// # Returns
///
/// A `CompressedMontgomeryU`.
pub fn compress(&self) -> CompressedMontgomeryU {
let u_affine: FieldElement = &self.U * &self.W.invert();
CompressedMontgomeryU(u_affine.to_bytes())
}
/// Differential addition for single-coordinate Montgomery points.
///
/// Montgomery coordinates in projective 𝗣¹ space are odd in that 𝗣¹
/// inherits none of the group structure from E_(A,B). Hence, the mapping
/// of the group operation, `⊕`, is undefined for the pair `(x(P), x(Q))`;
/// that is, given `x(P)` and `x(Q)`, we cannot derive `x(P ⊕ Q)`. This is
/// due to the fact that, in Montgomery coordinates, `x(P)` determines `P`
/// only up to a sign, and thus we cannot differentiate `x(P ⊕ Q)` from
/// `x(P ⊖ Q)`. However, via differential addition, any three of the values
/// `{x(P), x(Q), x(P ⊕ Q), x(P ⊖ Q)}` determines the forth, so we can
/// define *pseudo-addition* for a singular coordinate.
///
/// # Warning
///
/// If the `difference` is the identity point, or a two torsion point, the
/// results of this method are not correct, but instead result in `(0:0)`
/// (an invalid projective point in the Montgomery model).
///
/// The doubling case is degenerate, in that `P ⦵ Q ∉ {O,T}`, where `T` is
/// the two torsion point.
fn differential_add(&self, that: &MontgomeryPoint,
difference: &MontgomeryPoint) -> MontgomeryPoint {
// XXX Do we want these debug assertions? We would need to implement
// XXX is_two_torsion_point(). —isis
// debug_assert!(!difference.is_identity()); // P ⦵ Q ∉ {O,T}
// debug_assert!(!difference.is_two_torsion_point());
let v1: FieldElement = &(&self.U + &self.W) * &(&that.U - &that.W);
let v2: FieldElement = &(&self.U - &self.W) * &(&that.U + &that.W);
MontgomeryPoint {
U: &difference.W * &(&v1 + &v2).square(), // does reduction on square()
W: &difference.U * &(&v1 - &v2).square(), // does reduction on square()
}
}
/// Pseudo-doubling for single-coordinate Montgomery points.
///
/// Given a Montgomery U-coordinate of a point `P`, compute the
/// U-coordinate given by
///
///     differential_double: x(P) ⟼ x([2]P)
///
/// # Returns
///
/// A Montgomery point equal to doubling this one.
///
// XXX It seems possible that combining the differential_add() and
// XXX differential_double() methods would save a non-trivial amount of
// XXX computation in the ladder. —isis
fn differential_double(&self) -> MontgomeryPoint {
let mut v1: FieldElement;
let v2: FieldElement;
let v3: FieldElement;
v1 = (&self.U + &self.W).square();
v2 = (&self.U - &self.W).square();
let U: FieldElement = &v1 * &v2;
v1 -= &v2;
v3 = &(&constants::APLUS2_OVER_FOUR * &v1) + &v2;
let W: FieldElement = &v1 * &v3;
MontgomeryPoint{ U: U, W: W }
}
}
/// Multiply this `MontgomeryPoint` by a `Scalar`.
///
/// The reader is refered to §5.3 of ["Montgomery Curves and Their Arithmetic"
/// by Craig Costello and Benjamin Smith](https://eprint.iacr.org/2017/212.pdf)
/// for an overview of side-channel-free Montgomery laddering algorithms.
impl<'a, 'b> Mul<&'b Scalar> for &'a MontgomeryPoint {
type Output = MontgomeryPoint;
fn mul(self, scalar: &'b Scalar) -> MontgomeryPoint {
let mut x0: MontgomeryPoint = MontgomeryPoint::identity();
let mut x1: MontgomeryPoint = *self;
let bits: [i8; 256] = scalar.bits();
for i in (0..255).rev() {
let mask: u8 = (bits[i+1] ^ bits[i]) as u8;
debug_assert!(mask == 0 || mask == 1);
x0.conditional_swap(&mut x1, mask);
x1 = x0.differential_add(&x1, &self);
x0 = x0.differential_double();
}
x0.conditional_swap(&mut x1, bits[0] as u8);
x0
}
}
impl<'b> MulAssign<&'b Scalar> for MontgomeryPoint {
fn mul_assign(&mut self, scalar: &'b Scalar) {
let result = (self as &MontgomeryPoint) * scalar;
*self = result;
}
}
impl<'a, 'b> Mul<&'b MontgomeryPoint> for &'a Scalar {
type Output = MontgomeryPoint;
fn mul(self, point: &'b MontgomeryPoint) -> MontgomeryPoint {
point * &self
}
}
// ------------------------------------------------------------------------
// Tests
// ------------------------------------------------------------------------
#[cfg(test)]
mod test {
use constants::ED25519_BASEPOINT_TABLE;
use constants::BASE_COMPRESSED_MONTGOMERY;
use edwards::Identity;
use super::*;
/// The X25519 basepoint, in compressed Montgomery form.
static BASE_CMPRSSD_MONTY: CompressedMontgomeryU =
CompressedMontgomeryU([0x09, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00]);
use rand::OsRng;
/// Test Montgomery conversion against the X25519 basepoint.
#[test]
fn basepoint_to_montgomery() {
assert_eq!(constants::ED25519_BASEPOINT_POINT.compress_montgomery().unwrap(),
BASE_CMPRSSD_MONTY);
assert_eq!(constants::ED25519_BASEPOINT_POINT.to_montgomery().compress(),
BASE_COMPRESSED_MONTGOMERY);
}
/// Test Montgomery conversion against the X25519 basepoint.
#[test]
fn basepoint_from_montgomery() {
assert_eq!(BASE_CMPRSSD_MONTY.decompress().unwrap().compress_edwards(),
constants::BASE_CMPRSSD);
assert_eq!(BASE_COMPRESSED_MONTGOMERY,
constants::BASE_CMPRSSD.decompress().unwrap().to_montgomery().compress());
}
/// If u = -1, then v^2 = u*(u^2+486662*u+1) = 486660.
@ -189,14 +455,146 @@ mod test {
let minus_one = FieldElement::minus_one();
let minus_one_bytes = minus_one.to_bytes();
let div_by_zero_u = CompressedMontgomeryU(minus_one_bytes);
assert!(div_by_zero_u.decompress().is_none());
assert!(div_by_zero_u.decompress_edwards().is_none());
}
/// Montgomery compression of the identity point should
/// fail (it's sent to infinity).
/// Montgomery compression of the identity point should not fail (since the
/// mapping in `ProjectivePoint.to_montgomery()` should be valid for the
/// identity.
#[test]
fn identity_to_monty() {
let id = ExtendedPoint::identity();
assert!(id.compress_montgomery().is_none());
assert_eq!(id.to_montgomery().compress(), MontgomeryPoint::identity().compress());
}
#[test]
fn projective_to_affine_roundtrips() {
assert_eq!(BASE_COMPRESSED_MONTGOMERY.decompress().compress(),
BASE_COMPRESSED_MONTGOMERY);
}
#[test]
fn differential_double_matches_double() {
let p: ExtendedPoint = constants::ED25519_BASEPOINT_POINT.double();
let q: MontgomeryPoint = BASE_COMPRESSED_MONTGOMERY.decompress().differential_double();
assert_eq!(p.to_montgomery().compress(), q.compress());
}
#[test]
fn montgomery_ct_eq_ne() {
let mut csprng: OsRng = OsRng::new().unwrap();
let s1: Scalar = Scalar::random(&mut csprng);
let s2: Scalar = Scalar::random(&mut csprng);
let p1: MontgomeryPoint = (&s1 * &ED25519_BASEPOINT_TABLE).to_montgomery();
let p2: MontgomeryPoint = (&s2 * &ED25519_BASEPOINT_TABLE).to_montgomery();
assert_eq!(p1.ct_eq(&p2), 0);
}
#[test]
fn montgomery_ct_eq_eq() {
let mut csprng: OsRng = OsRng::new().unwrap();
let s1: Scalar = Scalar::random(&mut csprng);
let p1: MontgomeryPoint = (&s1 * &ED25519_BASEPOINT_TABLE).to_montgomery();
assert_eq!(p1.ct_eq(&p1), 1);
}
#[test]
fn differential_add_matches_edwards_model() {
let mut csprng: OsRng = OsRng::new().unwrap();
let s1: Scalar = Scalar::random(&mut csprng);
let s2: Scalar = Scalar::random(&mut csprng);
let p1: ExtendedPoint = &constants::ED25519_BASEPOINT_TABLE * &s1;
let p2: ExtendedPoint = &constants::ED25519_BASEPOINT_TABLE * &s2;
let diff: ExtendedPoint = &p1 - &p2;
let p1m: MontgomeryPoint = p1.to_montgomery();
let p2m: MontgomeryPoint = p2.to_montgomery();
let diffm: MontgomeryPoint = diff.to_montgomery();
let result = p1m.differential_add(&p2m, &diffm);
assert_eq!(result.compress(), (&p1 + &p2).to_montgomery().compress());
}
#[test]
fn ladder_matches_scalarmult() {
let mut csprng: OsRng = OsRng::new().unwrap();
let s: Scalar = Scalar::random(&mut csprng);
let p_edwards: ExtendedPoint = &constants::ED25519_BASEPOINT_TABLE * &s;
let p_montgomery: MontgomeryPoint = p_edwards.to_montgomery();
let expected = &s * &p_edwards;
let result = &s * &p_montgomery;
assert_eq!(result.compress(), expected.to_montgomery().compress())
}
#[test]
fn ladder_basepoint_times_two_matches_double() {
let two: Scalar = Scalar::from_u64(2u64);
let result: MontgomeryPoint = &BASE_COMPRESSED_MONTGOMERY.decompress() * &two;
let expected: ExtendedPoint = constants::ED25519_BASEPOINT_POINT.double();
assert_eq!(result.compress(), expected.to_montgomery().compress());
}
#[test]
#[should_panic(expected = "assertion failed: self[31] <= 127")]
fn ladder_matches_scalarmult_with_scalar_high_bit_set() {
let mut s: Scalar = Scalar::one();
s[31] = 255;
let result: MontgomeryPoint = &BASE_COMPRESSED_MONTGOMERY.decompress() * &s;
let expected: ExtendedPoint = &constants::ED25519_BASEPOINT_TABLE * &s;
assert_eq!(result.compress(), expected.to_montgomery().compress())
}
}
#[cfg(all(test, feature = "bench"))]
mod bench {
use rand::OsRng;
use constants::ED25519_BASEPOINT_TABLE;
use constants::BASE_COMPRESSED_MONTGOMERY;
use test::Bencher;
use super::*;
#[bench]
fn montgomery_ct_eq(b: &mut Bencher) {
let mut csprng: OsRng = OsRng::new().unwrap();
let s1: Scalar = Scalar::random(&mut csprng);
let s2: Scalar = Scalar::random(&mut csprng);
let p1: MontgomeryPoint = (&s1 * &ED25519_BASEPOINT_TABLE).to_montgomery();
let p2: MontgomeryPoint = (&s2 * &ED25519_BASEPOINT_TABLE).to_montgomery();
b.iter(| | p1.ct_eq(&p2))
}
#[bench]
fn montgomery_decompress(b: &mut Bencher) {
b.iter(| | BASE_COMPRESSED_MONTGOMERY.decompress());
}
#[bench]
fn montgomery_compress(b: &mut Bencher) {
let p: MontgomeryPoint = BASE_COMPRESSED_MONTGOMERY.decompress();
b.iter(| | p.compress());
}
#[bench]
fn montgomery_ladder(b: &mut Bencher) {
let mut csprng: OsRng = OsRng::new().unwrap();
let s: Scalar = Scalar::random(&mut csprng);
let p: MontgomeryPoint = (&Scalar::random(&mut csprng) * &ED25519_BASEPOINT_TABLE).to_montgomery();
b.iter(| | &s * &p);
}
}

429
src/scalar.rs
View file

@ -2,11 +2,13 @@
//
// This file is part of curve25519-dalek.
// Copyright (c) 2016-2017 Isis Lovecruft, Henry de Valence
// Portions Copyright 2017 Brian Smith
// See LICENSE for licensing information.
//
// Authors:
// - Isis Agora Lovecruft <isis@patternsinthevoid.net>
// - Henry de Valence <hdevalence@hdevalence.ca>
// - Brian Smith <brian@briansmith.org>
//! Arithmetic for scalar multiplication.
//!
@ -20,12 +22,11 @@
//!
//! The `Scalar` struct represents an element in ℤ/lℤ.
//!
//! Arithmetic operations on `Scalar`s are done using 12 21-bit limbs.
//! However, in contrast to `FieldElement`s, `Scalar`s are stored in
//! In contrast to `FieldElement`s, `Scalar`s are stored in
//! memory as bytes, allowing easy access to the bits of the `Scalar`
//! when multiplying a point by a scalar. For efficient arithmetic
//! between two scalars, the `UnpackedScalar` struct is stored as
//! limbs.
//! between two scalars, the `UnpackedScalar` struct (internally
//! either `Scalar32` or `Scalar64`) is stored as limbs.
use core::fmt::Debug;
use core::ops::Neg;
@ -41,9 +42,6 @@ use rand::Rng;
use digest::Digest;
use generic_array::typenum::U64;
use constants;
use utils::{load3, load4};
use subtle::slices_equal;
use subtle::ConditionallyAssignable;
use subtle::Equal;
@ -106,50 +104,47 @@ impl IndexMut<usize> for Scalar {
impl<'b> MulAssign<&'b Scalar> for Scalar {
fn mul_assign(&mut self, _rhs: &'b Scalar) {
let result = (self as &Scalar) * _rhs;
self.0 = result.0;
*self = Scalar::mul(self, _rhs)
}
}
impl<'a, 'b> Mul<&'b Scalar> for &'a Scalar {
type Output = Scalar;
fn mul(self, _rhs: &'b Scalar) -> Scalar {
Scalar::multiply_add(self, _rhs, &Scalar::zero())
Scalar::mul(self, _rhs)
}
}
impl<'b> AddAssign<&'b Scalar> for Scalar {
fn add_assign(&mut self, _rhs: &'b Scalar) {
*self = Scalar::multiply_add(&Scalar::one(), self, _rhs);
*self = Scalar::add(self, _rhs);
}
}
impl<'a, 'b> Add<&'b Scalar> for &'a Scalar {
type Output = Scalar;
fn add(self, _rhs: &'b Scalar) -> Scalar {
Scalar::multiply_add(&Scalar::one(), self, _rhs)
Scalar::add(self, _rhs)
}
}
impl<'b> SubAssign<&'b Scalar> for Scalar {
fn sub_assign(&mut self, _rhs: &'b Scalar) {
// (l-1)*_rhs + self = self - _rhs
*self = Scalar::multiply_add(&constants::l_minus_1, _rhs, self);
*self = Scalar::sub(self, _rhs);
}
}
impl<'a, 'b> Sub<&'b Scalar> for &'a Scalar {
type Output = Scalar;
fn sub(self, _rhs: &'b Scalar) -> Scalar {
// (l-1)*_rhs + self = self - _rhs
Scalar::multiply_add(&constants::l_minus_1, _rhs, self)
Scalar::sub(self, _rhs)
}
}
impl<'a> Neg for &'a Scalar {
type Output = Scalar;
fn neg(self) -> Scalar {
self * &constants::l_minus_1
Scalar::sub(&Scalar::zero(), self)
}
}
@ -232,6 +227,18 @@ impl<'de> Deserialize<'de> for Scalar {
}
}
/// An `UnpackedScalar` represents an element of the field GF(l), optimized for speed.
#[cfg(feature="radix_51")]
type UnpackedScalar = Scalar64;
#[cfg(feature="radix_51")]
use scalar_64bit::*;
/// An `UnpackedScalar` represents an element of the field GF(l), optimized for speed.
#[cfg(not(feature="radix_51"))]
type UnpackedScalar = Scalar32;
#[cfg(not(feature="radix_51"))]
use scalar_32bit::*;
impl Scalar {
/// Return a `Scalar` chosen uniformly at random using a user-provided RNG.
///
@ -382,26 +389,6 @@ impl Scalar {
naf
}
// Unpack a scalar into 12 21-bit limbs.
fn unpack(&self) -> UnpackedScalar {
let mask_21bits: i64 = (1 << 21) - 1;
let mut a = UnpackedScalar([0i64; 12]);
a[ 0] = mask_21bits & load3(&self.0[ 0..]) ;
a[ 1] = mask_21bits & (load4(&self.0[ 2..]) >> 5);
a[ 2] = mask_21bits & (load3(&self.0[ 5..]) >> 2);
a[ 3] = mask_21bits & (load4(&self.0[ 7..]) >> 7);
a[ 4] = mask_21bits & (load4(&self.0[10..]) >> 4);
a[ 5] = mask_21bits & (load3(&self.0[13..]) >> 1);
a[ 6] = mask_21bits & (load4(&self.0[15..]) >> 6);
a[ 7] = mask_21bits & (load3(&self.0[18..]) >> 3);
a[ 8] = mask_21bits & load3(&self.0[21..]) ;
a[ 9] = mask_21bits & (load4(&self.0[23..]) >> 5);
a[10] = mask_21bits & (load3(&self.0[26..]) >> 2);
a[11] = load4(&self.0[28..]) >> 7 ;
a
}
/// Write this scalar in radix 16, with coefficients in `[-8,8)`,
/// i.e., compute `a_i` such that
///
@ -440,293 +427,106 @@ impl Scalar {
output
}
/// Compute `ab+c (mod l)`.
/// XXX should this exist, or should we just have Mul, Add etc impls
/// that unpack and then call UnpackedScalar::multiply_add ?
pub fn multiply_add(a: &Scalar, b: &Scalar, c: &Scalar) -> Scalar {
// Unpack scalars into limbs
let al = a.unpack();
let bl = b.unpack();
let cl = c.unpack();
/// Unpack this `Scalar` to an `UnpackedScalar`
pub fn unpack(&self) -> UnpackedScalar {
UnpackedScalar::from_bytes(&self.0)
}
// Multiply and repack
UnpackedScalar::multiply_add(&al, &bl, &cl).pack()
/// Compute `a + b` (mod l)
pub fn add(a: &Scalar, b: &Scalar) -> Scalar {
UnpackedScalar::add(&a.unpack(), &b.unpack()).pack()
}
/// Compute `a - b` (mod l).
pub fn sub(a: &Scalar, b: &Scalar) -> Scalar {
UnpackedScalar::sub(&a.unpack(), &b.unpack()).pack()
}
/// Compute `a * b` (mod l).
pub fn mul(a: &Scalar, b: &Scalar) -> Scalar {
UnpackedScalar::mul(&a.unpack(), &b.unpack()).pack()
}
/// Compute `(a * b) + c` (mod l).
pub fn multiply_add(a: &Scalar, b: &Scalar, c: &Scalar) -> Scalar {
UnpackedScalar::add(&UnpackedScalar::mul(&a.unpack(), &b.unpack()), &c.unpack()).pack()
}
/// Reduce a 512-bit little endian number mod l
pub fn reduce(input: &[u8; 64]) -> Scalar {
let mut s = [0i64; 24];
// XXX express this as two unpack_limbs
// some issues re: masking with the top byte of the 32byte input
let mask_21bits: i64 = (1 << 21) -1;
s[0] = mask_21bits & load3(&input[ 0..]) ;
s[1] = mask_21bits & (load4(&input[ 2..]) >> 5);
s[2] = mask_21bits & (load3(&input[ 5..]) >> 2);
s[3] = mask_21bits & (load4(&input[ 7..]) >> 7);
s[4] = mask_21bits & (load4(&input[10..]) >> 4);
s[5] = mask_21bits & (load3(&input[13..]) >> 1);
s[6] = mask_21bits & (load4(&input[15..]) >> 6);
s[7] = mask_21bits & (load3(&input[18..]) >> 3);
s[8] = mask_21bits & load3(&input[21..]) ;
s[9] = mask_21bits & (load4(&input[23..]) >> 5);
s[10] = mask_21bits & (load3(&input[26..]) >> 2);
s[11] = mask_21bits & (load4(&input[28..]) >> 7);
s[12] = mask_21bits & (load4(&input[31..]) >> 4);
s[13] = mask_21bits & (load3(&input[34..]) >> 1);
s[14] = mask_21bits & (load4(&input[36..]) >> 6);
s[15] = mask_21bits & (load3(&input[39..]) >> 3);
s[16] = mask_21bits & load3(&input[42..]) ;
s[17] = mask_21bits & (load4(&input[44..]) >> 5);
s[18] = mask_21bits & (load3(&input[47..]) >> 2);
s[19] = mask_21bits & (load4(&input[49..]) >> 7);
s[20] = mask_21bits & (load4(&input[52..]) >> 4);
s[21] = mask_21bits & (load3(&input[55..]) >> 1);
s[22] = mask_21bits & (load4(&input[57..]) >> 6);
s[23] = load4(&input[60..]) >> 3 ;
// XXX replacing the previous code in this function with the
// call to reduce_limbs adds two extra carry passes (the ones
// at the top of the reduce_limbs function). Otherwise they
// are identical. The test seems to work OK but it would be
// good to check that this really is OK to add.
UnpackedScalar::reduce_limbs(&mut s).pack()
}
}
/// The `UnpackedScalar` struct represents an element in ℤ/lℤ as 12
/// 21-bit limbs.
#[derive(Copy,Clone)]
pub struct UnpackedScalar(pub [i64; 12]);
impl Index<usize> for UnpackedScalar {
type Output = i64;
fn index(&self, _index: usize) -> &i64 {
&(self.0[_index])
}
}
impl IndexMut<usize> for UnpackedScalar {
fn index_mut(&mut self, _index: usize) -> &mut i64 {
&mut (self.0[_index])
UnpackedScalar::from_bytes_wide(input).pack()
}
}
impl UnpackedScalar {
/// Pack the limbs of this `UnpackedScalar` into a `Scalar`.
fn pack(&self) -> Scalar {
let mut s = Scalar::zero();
s[0] = (self.0[ 0] >> 0) as u8;
s[1] = (self.0[ 0] >> 8) as u8;
s[2] = ((self.0[ 0] >> 16) | (self.0[ 1] << 5)) as u8;
s[3] = (self.0[ 1] >> 3) as u8;
s[4] = (self.0[ 1] >> 11) as u8;
s[5] = ((self.0[ 1] >> 19) | (self.0[ 2] << 2)) as u8;
s[6] = (self.0[ 2] >> 6) as u8;
s[7] = ((self.0[ 2] >> 14) | (self.0[ 3] << 7)) as u8;
s[8] = (self.0[ 3] >> 1) as u8;
s[9] = (self.0[ 3] >> 9) as u8;
s[10] = ((self.0[ 3] >> 17) | (self.0[ 4] << 4)) as u8;
s[11] = (self.0[ 4] >> 4) as u8;
s[12] = (self.0[ 4] >> 12) as u8;
s[13] = ((self.0[ 4] >> 20) | (self.0[ 5] << 1)) as u8;
s[14] = (self.0[ 5] >> 7) as u8;
s[15] = ((self.0[ 5] >> 15) | (self.0[ 6] << 6)) as u8;
s[16] = (self.0[ 6] >> 2) as u8;
s[17] = (self.0[ 6] >> 10) as u8;
s[18] = ((self.0[ 6] >> 18) | (self.0[ 7] << 3)) as u8;
s[19] = (self.0[ 7] >> 5) as u8;
s[20] = (self.0[ 7] >> 13) as u8;
s[21] = (self.0[ 8] >> 0) as u8;
s[22] = (self.0[ 8] >> 8) as u8;
s[23] = ((self.0[ 8] >> 16) | (self.0[ 9] << 5)) as u8;
s[24] = (self.0[ 9] >> 3) as u8;
s[25] = (self.0[ 9] >> 11) as u8;
s[26] = ((self.0[ 9] >> 19) | (self.0[10] << 2)) as u8;
s[27] = (self.0[10] >> 6) as u8;
s[28] = ((self.0[10] >> 14) | (self.0[11] << 7)) as u8;
s[29] = (self.0[11] >> 1) as u8;
s[30] = (self.0[11] >> 9) as u8;
s[31] = (self.0[11] >> 17) as u8;
s
}
/// Return the zero scalar.
pub fn zero() -> UnpackedScalar {
UnpackedScalar([0,0,0,0,0,0,0,0,0,0,0,0])
}
/// Return the one scalar.
pub fn one() -> UnpackedScalar {
UnpackedScalar([1,0,0,0,0,0,0,0,0,0,0,0])
Scalar(self.to_bytes())
}
/// Compute the multiplicative inverse of this scalar.
pub fn invert(&self) -> UnpackedScalar {
let mut y = UnpackedScalar::one();
// Run through bits of l-2 from highest to least
for bit in constants::l_minus_2.bits().iter().rev() {
y = UnpackedScalar::multiply_add(&y, &y, &UnpackedScalar::zero());
if *bit == 1 {
y = UnpackedScalar::multiply_add(&y, self, &UnpackedScalar::zero());
// This is a direct transliteration of the addition chain from
// https://briansmith.org/ecc-inversion-addition-chains-01#curve25519_scalar_inversion
// as it was published on 2017-09-03.
let _1 = self.to_montgomery();
let _10 = _1.montgomery_square();
let _100 = _10.montgomery_square();
let _11 = UnpackedScalar::montgomery_mul(&_10, &_1);
let _101 = UnpackedScalar::montgomery_mul(&_10, &_11);
let _111 = UnpackedScalar::montgomery_mul(&_10, &_101);
let _1001 = UnpackedScalar::montgomery_mul(&_10, &_111);
let _1011 = UnpackedScalar::montgomery_mul(&_10, &_1001);
let _1111 = UnpackedScalar::montgomery_mul(&_100, &_1011);
// _10000
let mut y = UnpackedScalar::montgomery_mul(&_1111, &_1);
#[inline]
fn square_multiply(y: &mut UnpackedScalar, squarings: usize, x: &UnpackedScalar) {
for _ in 0..squarings {
*y = y.montgomery_square();
}
}
y
}
/// Compute `ab+c (mod l)`.
pub fn multiply_add(a: &UnpackedScalar,
b: &UnpackedScalar,
c: &UnpackedScalar) -> UnpackedScalar {
let mut result = [0i64; 24];
// Multiply a and b, and add c
result[0] = c[0] + a[0]*b[0];
result[1] = c[1] + a[0]*b[1] + a[1]*b[0];
result[2] = c[2] + a[0]*b[2] + a[1]*b[1] + a[2]*b[0];
result[3] = c[3] + a[0]*b[3] + a[1]*b[2] + a[2]*b[1] + a[3]*b[0];
result[4] = c[4] + a[0]*b[4] + a[1]*b[3] + a[2]*b[2] + a[3]*b[1] + a[4]*b[0];
result[5] = c[5] + a[0]*b[5] + a[1]*b[4] + a[2]*b[3] + a[3]*b[2] + a[4]*b[1] + a[5]*b[0];
result[6] = c[6] + a[0]*b[6] + a[1]*b[5] + a[2]*b[4] + a[3]*b[3] + a[4]*b[2] + a[5]*b[1] + a[6]*b[0];
result[7] = c[7] + a[0]*b[7] + a[1]*b[6] + a[2]*b[5] + a[3]*b[4] + a[4]*b[3] + a[5]*b[2] + a[6]*b[1] + a[7]*b[0];
result[8] = c[8] + a[0]*b[8] + a[1]*b[7] + a[2]*b[6] + a[3]*b[5] + a[4]*b[4] + a[5]*b[3] + a[6]*b[2] + a[7]*b[1] + a[8]*b[0];
result[9] = c[9] + a[0]*b[9] + a[1]*b[8] + a[2]*b[7] + a[3]*b[6] + a[4]*b[5] + a[5]*b[4] + a[6]*b[3] + a[7]*b[2] + a[8]*b[1] + a[9]*b[0];
result[10] = c[10] + a[0]*b[10] + a[1]*b[9] + a[2]*b[8] + a[3]*b[7] + a[4]*b[6] + a[5]*b[5] + a[6]*b[4] + a[7]*b[3] + a[8]*b[2] + a[9]*b[1] + a[10]*b[0];
result[11] = c[11] + a[0]*b[11] + a[1]*b[10] + a[2]*b[9] + a[3]*b[8] + a[4]*b[7] + a[5]*b[6] + a[6]*b[5] + a[7]*b[4] + a[8]*b[3] + a[9]*b[2] + a[10]*b[1] + a[11]*b[0];
result[12] = a[1]*b[11] + a[2]*b[10] + a[3]*b[9] + a[4]*b[8] + a[5]*b[7] + a[6]*b[6] + a[7]*b[5] + a[8]*b[4] + a[9]*b[3] + a[10]*b[2] + a[11]*b[1];
result[13] = a[2]*b[11] + a[3]*b[10] + a[4]*b[9] + a[5]*b[8] + a[6]*b[7] + a[7]*b[6] + a[8]*b[5] + a[9]*b[4] + a[10]*b[3] + a[11]*b[2];
result[14] = a[3]*b[11] + a[4]*b[10] + a[5]*b[9] + a[6]*b[8] + a[7]*b[7] + a[8]*b[6] + a[9]*b[5] + a[10]*b[4] + a[11]*b[3];
result[15] = a[4]*b[11] + a[5]*b[10] + a[6]*b[9] + a[7]*b[8] + a[8]*b[7] + a[9]*b[6] + a[10]*b[5] + a[11]*b[4];
result[16] = a[5]*b[11] + a[6]*b[10] + a[7]*b[9] + a[8]*b[8] + a[9]*b[7] + a[10]*b[6] + a[11]*b[5];
result[17] = a[6]*b[11] + a[7]*b[10] + a[8]*b[9] + a[9]*b[8] + a[10]*b[7] + a[11]*b[6];
result[18] = a[7]*b[11] + a[8]*b[10] + a[9]*b[9] + a[10]*b[8] + a[11]*b[7];
result[19] = a[8]*b[11] + a[9]*b[10] + a[10]*b[9] + a[11]*b[8];
result[20] = a[9]*b[11] + a[10]*b[10] + a[11]*b[9];
result[21] = a[10]*b[11] + a[11]*b[10];
result[22] = a[11]*b[11];
result[23] = 0i64;
// Reduce limbs
UnpackedScalar::reduce_limbs(&mut result)
}
/// Reduce 24 limbs to 12, consuming the input. Reduction is mod
///
/// l = 2^252 + 27742317777372353535851937790883648493,
///
/// so
///
/// 2^252 = -27742317777372353535851937790883648493 (mod l).
///
/// We can write the right-hand side in 21-bit limbs as
///
/// rhs = 666643 * 2^0
/// + 470296 * 2^21
/// + 654183 * 2^42
/// - 997805 * 2^63
/// + 136657 * 2^84
/// - 683901 * 2^105
///
/// The (12+k)-th limb of `limbs` is the coefficient of
///
/// 2^(252 + 21*k)
///
/// since 12*21 = 252. By the above, we have that
///
/// c * 2^(252 + 21*k) = c * 666643 * 2^(21*k)
/// + c * 470296 * 2^(42*k) + ...
///
/// so we can eliminate it by adding those values to the lower
/// limbs. Reduction mod l amounts to eliminating all of the
/// high limbs while carrying as appropriate to prevent
/// overflows in the lower limbs.
fn reduce_limbs(mut limbs: &mut [i64; 24]) -> UnpackedScalar {
#[inline]
#[allow(dead_code)]
fn do_reduction(limbs: &mut [i64; 24], i: usize) {
limbs[i - 12] += limbs[i] * 666643;
limbs[i - 11] += limbs[i] * 470296;
limbs[i - 10] += limbs[i] * 654183;
limbs[i - 9] -= limbs[i] * 997805;
limbs[i - 8] += limbs[i] * 136657;
limbs[i - 7] -= limbs[i] * 683901;
limbs[i] = 0;
}
/// Carry excess from the `i`-th limb into the `(i+1)`-th limb.
/// Postcondition: `0 <= limbs[i] < 2^21`.
#[inline]
#[allow(dead_code)]
fn do_carry_uncentered(limbs: &mut [i64; 24], i: usize) {
let carry: i64 = limbs[i] >> 21;
limbs[i+1] += carry;
limbs[i ] -= carry << 21;
}
#[inline]
#[allow(dead_code)]
/// Carry excess from the `i`-th limb into the `(i+1)`-th limb.
/// Postcondition: `-2^20 <= limbs[i] < 2^20`.
fn do_carry_centered(limbs: &mut [i64; 24], i: usize) {
let carry: i64 = (limbs[i] + (1<<20)) >> 21;
limbs[i+1] += carry;
limbs[i ] -= carry << 21;
*y = UnpackedScalar::montgomery_mul(y, x);
}
for i in 0..23 {
do_carry_centered(&mut limbs, i);
}
for i in (0..23).filter(|x| x % 2 == 1) {
do_carry_centered(&mut limbs, i);
}
square_multiply(&mut y, 123 + 3, &_101);
square_multiply(&mut y, 2 + 2, &_11);
square_multiply(&mut y, 1 + 4, &_1111);
square_multiply(&mut y, 1 + 4, &_1111);
square_multiply(&mut y, 4, &_1001);
square_multiply(&mut y, 2, &_11);
square_multiply(&mut y, 1 + 4, &_1111);
square_multiply(&mut y, 1 + 3, &_101);
square_multiply(&mut y, 3 + 3, &_101);
square_multiply(&mut y, 3, &_111);
square_multiply(&mut y, 1 + 4, &_1111);
square_multiply(&mut y, 2 + 3, &_111);
square_multiply(&mut y, 2 + 2, &_11);
square_multiply(&mut y, 1 + 4, &_1011);
square_multiply(&mut y, 2 + 4, &_1011);
square_multiply(&mut y, 6 + 4, &_1001);
square_multiply(&mut y, 2 + 2, &_11);
square_multiply(&mut y, 3 + 2, &_11);
square_multiply(&mut y, 3 + 2, &_11);
square_multiply(&mut y, 1 + 4, &_1001);
square_multiply(&mut y, 1 + 3, &_111);
square_multiply(&mut y, 2 + 4, &_1111);
square_multiply(&mut y, 1 + 4, &_1011);
square_multiply(&mut y, 3, &_101);
square_multiply(&mut y, 2 + 4, &_1111);
square_multiply(&mut y, 3, &_101);
square_multiply(&mut y, 1 + 2, &_11);
do_reduction(&mut limbs, 23);
do_reduction(&mut limbs, 22);
do_reduction(&mut limbs, 21);
do_reduction(&mut limbs, 20);
do_reduction(&mut limbs, 19);
do_reduction(&mut limbs, 18);
for i in (6..18).filter(|x| x % 2 == 0) {
do_carry_centered(&mut limbs, i);
}
for i in (6..16).filter(|x| x % 2 == 1) {
do_carry_centered(&mut limbs, i);
}
do_reduction(&mut limbs, 17);
do_reduction(&mut limbs, 16);
do_reduction(&mut limbs, 15);
do_reduction(&mut limbs, 14);
do_reduction(&mut limbs, 13);
do_reduction(&mut limbs, 12);
for i in (0..12).filter(|x| x % 2 == 0) {
do_carry_centered(&mut limbs, i);
}
for i in (0..12).filter(|x| x % 2 == 1) {
do_carry_centered(&mut limbs, i);
}
do_reduction(&mut limbs, 12);
for i in 0..12 {
do_carry_uncentered(&mut limbs, i);
}
do_reduction(&mut limbs, 12);
for i in 0..11 {
do_carry_uncentered(&mut limbs, i);
}
UnpackedScalar(*array_ref!(limbs, 0, 12))
y.from_montgomery()
}
}
#[cfg(test)]
mod test {
use super::*;
use constants;
/// x = 2238329342913194256032495932344128051776374960164957527413114840482143558222
pub static X: Scalar = Scalar(
@ -734,6 +534,12 @@ mod test {
0x59, 0x13, 0xb4, 0x64, 0x1b, 0xc2, 0x7d, 0x52,
0x52, 0xa5, 0x85, 0x10, 0x1b, 0xcc, 0x42, 0x44,
0xd4, 0x49, 0xf4, 0xa8, 0x79, 0xd9, 0xf2, 0x04]);
/// 1/x = 6859937278830797291664592131120606308688036382723378951768035303146619657244
pub static XINV: Scalar = Scalar(
[0x1c, 0xdc, 0x17, 0xfc, 0xe0, 0xe9, 0xa5, 0xbb,
0xd9, 0x24, 0x7e, 0x56, 0xbb, 0x01, 0x63, 0x47,
0xbb, 0xba, 0x31, 0xed, 0xd5, 0xa9, 0xbb, 0x96,
0xd5, 0x0b, 0xcd, 0x7a, 0x3f, 0x96, 0x2a, 0x0f]);
/// y = 2592331292931086675770238855846338635550719849568364935475441891787804997264
pub static Y: Scalar = Scalar(
[0x90, 0x76, 0x33, 0xfe, 0x1c, 0x4b, 0x66, 0xa4,
@ -838,6 +644,15 @@ mod test {
}
}
#[test]
fn square() {
let expected = Scalar::multiply_add(&X, &X, &Scalar::zero());
let actual = X.unpack().square().pack();
for i in 0..32 {
assert!(expected[i] == actual[i]);
}
}
#[test]
fn scalar_reduce() {
let mut bignum = [0u8; 64];
@ -862,6 +677,7 @@ mod test {
#[test]
fn invert() {
let inv_X = X.invert();
assert_eq!(inv_X, XINV);
let should_be_one = &inv_X * &X;
assert_eq!(should_be_one, Scalar::one());
}
@ -894,7 +710,7 @@ mod bench {
use test::Bencher;
use super::*;
use super::test::{X, Y, Z};
use super::test::{X};
#[bench]
fn scalar_random(b: &mut Bencher) {
@ -903,22 +719,9 @@ mod bench {
b.iter(|| Scalar::random(&mut csprng));
}
#[bench]
fn scalar_multiply_add(b: &mut Bencher) {
b.iter(|| Scalar::multiply_add(&X, &Y, &Z));
}
#[bench]
fn invert(b: &mut Bencher) {
let x = X.unpack();
b.iter(|| x.invert());
}
#[bench]
fn scalar_unpacked_multiply_add(b: &mut Bencher) {
let x = X.unpack();
let y = Y.unpack();
let z = Z.unpack();
b.iter(|| UnpackedScalar::multiply_add(&x, &y, &z));
}
}

559
src/scalar_32bit.rs Normal file
View file

@ -0,0 +1,559 @@
//! Arithmetic mod 2^252 + 27742317777372353535851937790883648493
//! with 9 29-bit unsigned limbs
//!
//! To see that this is safe for intermediate results, note that
//! the largest limb in a 9 by 9 product of 29-bit limbs will be
//! (0x1fffffff^2) * 9 = 0x23fffffdc0000009 (62 bits).
//!
//! For a one level Karatsuba decomposition, the specific ranges
//! depend on how the limbs are combined, but will stay within
//! -0x1ffffffe00000008 (62 bits with sign bit) to
//! 0x43fffffbc0000011 (63 bits), which is still safe.
//!
//! (the 9th limb will never exceed 21 bits, so the actual
//! ranges are slightly smaller)
use core::fmt::Debug;
use core::ops::{Index, IndexMut};
use constants;
/// The `Scalar32` struct represents an element in ℤ/lℤ as 9 29-bit limbs
#[derive(Copy,Clone)]
pub struct Scalar32(pub [u32; 9]);
impl Debug for Scalar32 {
fn fmt(&self, f: &mut ::core::fmt::Formatter) -> ::core::fmt::Result {
write!(f, "Scalar32: {:?}", &self.0[..])
}
}
impl Index<usize> for Scalar32 {
type Output = u32;
fn index(&self, _index: usize) -> &u32 {
&(self.0[_index])
}
}
impl IndexMut<usize> for Scalar32 {
fn index_mut(&mut self, _index: usize) -> &mut u32 {
&mut (self.0[_index])
}
}
/// u32 * u32 = u64 multiply helper
#[inline(always)]
fn m(x: u32, y: u32) -> u64 {
(x as u64) * (y as u64)
}
impl Scalar32 {
/// Return the zero scalar.
pub fn zero() -> Scalar32 {
Scalar32([0,0,0,0,0,0,0,0,0])
}
/// Unpack a 32 byte / 512 bit scalar into 9 29-bit limbs, ignoring the upper 3 bits.
pub fn from_bytes(bytes: &[u8; 32]) -> Scalar32 {
let mut words = [0u32; 8];
for i in 0..8 {
for j in 0..4 {
words[i] |= (bytes[(i * 4) + j] as u32) << (j * 8);
}
}
let mask = (1u32 << 29) - 1;
let top_mask = (1u32 << 21) - 1;
let mut s = Scalar32::zero();
s[ 0] = words[0] & mask;
s[ 1] = ((words[0] >> 29) | (words[1] << 3)) & mask;
s[ 2] = ((words[1] >> 26) | (words[2] << 6)) & mask;
s[ 3] = ((words[2] >> 23) | (words[3] << 9)) & mask;
s[ 4] = ((words[3] >> 20) | (words[4] << 12)) & mask;
s[ 5] = ((words[4] >> 17) | (words[5] << 15)) & mask;
s[ 6] = ((words[5] >> 14) | (words[6] << 18)) & mask;
s[ 7] = ((words[6] >> 11) | (words[7] << 21)) & mask;
s[ 8] = (words[7] >> 8) & top_mask;
s
}
/// Reduce a 64 byte / 512 bit scalar mod l.
pub fn from_bytes_wide(bytes: &[u8; 64]) -> Scalar32 {
let mut words = [0u32; 16];
for i in 0..16 {
for j in 0..4 {
words[i] |= (bytes[(i * 4) + j] as u32) << (j * 8);
}
}
let mask = (1u32 << 29) - 1;
let mut lo = Scalar32::zero();
let mut hi = Scalar32::zero();
lo[0] = words[ 0] & mask;
lo[1] = ((words[ 0] >> 29) | (words[ 1] << 3)) & mask;
lo[2] = ((words[ 1] >> 26) | (words[ 2] << 6)) & mask;
lo[3] = ((words[ 2] >> 23) | (words[ 3] << 9)) & mask;
lo[4] = ((words[ 3] >> 20) | (words[ 4] << 12)) & mask;
lo[5] = ((words[ 4] >> 17) | (words[ 5] << 15)) & mask;
lo[6] = ((words[ 5] >> 14) | (words[ 6] << 18)) & mask;
lo[7] = ((words[ 6] >> 11) | (words[ 7] << 21)) & mask;
lo[8] = ((words[ 7] >> 8) | (words[ 8] << 24)) & mask;
hi[0] = ((words[ 8] >> 5) | (words[ 9] << 27)) & mask;
hi[1] = (words[ 9] >> 2) & mask;
hi[2] = ((words[ 9] >> 31) | (words[10] << 1)) & mask;
hi[3] = ((words[10] >> 28) | (words[11] << 4)) & mask;
hi[4] = ((words[11] >> 25) | (words[12] << 7)) & mask;
hi[5] = ((words[12] >> 22) | (words[13] << 10)) & mask;
hi[6] = ((words[13] >> 19) | (words[14] << 13)) & mask;
hi[7] = ((words[14] >> 16) | (words[15] << 16)) & mask;
hi[8] = (words[15] >> 13) & mask;
lo = Scalar32::montgomery_mul(&lo, &constants::R); // (lo * R) / R = lo
hi = Scalar32::montgomery_mul(&hi, &constants::RR); // (hi * R^2) / R = hi * R
Scalar32::add(&hi, &lo) // (hi * R) + lo
}
/// Pack the limbs of this `Scalar32` into 32 bytes.
pub fn to_bytes(&self) -> [u8; 32] {
let mut s = [0u8; 32];
s[0] = (self.0[ 0] >> 0) as u8;
s[1] = (self.0[ 0] >> 8) as u8;
s[2] = (self.0[ 0] >> 16) as u8;
s[3] = ((self.0[ 0] >> 24) | (self.0[ 1] << 5)) as u8;
s[4] = (self.0[ 1] >> 3) as u8;
s[5] = (self.0[ 1] >> 11) as u8;
s[6] = (self.0[ 1] >> 19) as u8;
s[7] = ((self.0[ 1] >> 27) | (self.0[ 2] << 2)) as u8;
s[8] = (self.0[ 2] >> 6) as u8;
s[9] = (self.0[ 2] >> 14) as u8;
s[10] = ((self.0[ 2] >> 22) | (self.0[ 3] << 7)) as u8;
s[11] = (self.0[ 3] >> 1) as u8;
s[12] = (self.0[ 3] >> 9) as u8;
s[13] = (self.0[ 3] >> 17) as u8;
s[14] = ((self.0[ 3] >> 25) | (self.0[ 4] << 4)) as u8;
s[15] = (self.0[ 4] >> 4) as u8;
s[16] = (self.0[ 4] >> 12) as u8;
s[17] = (self.0[ 4] >> 20) as u8;
s[18] = ((self.0[ 4] >> 28) | (self.0[ 5] << 1)) as u8;
s[19] = (self.0[ 5] >> 7) as u8;
s[20] = (self.0[ 5] >> 15) as u8;
s[21] = ((self.0[ 5] >> 23) | (self.0[ 6] << 6)) as u8;
s[22] = (self.0[ 6] >> 2) as u8;
s[23] = (self.0[ 6] >> 10) as u8;
s[24] = (self.0[ 6] >> 18) as u8;
s[25] = ((self.0[ 6] >> 26) | (self.0[ 7] << 3)) as u8;
s[26] = (self.0[ 7] >> 5) as u8;
s[27] = (self.0[ 7] >> 13) as u8;
s[28] = (self.0[ 7] >> 21) as u8;
s[29] = (self.0[ 8] >> 0) as u8;
s[30] = (self.0[ 8] >> 8) as u8;
s[31] = (self.0[ 8] >> 16) as u8;
s
}
/// Compute `a + b` (mod l).
pub fn add(a: &Scalar32, b: &Scalar32) -> Scalar32 {
let mut sum = Scalar32::zero();
let mask = (1u32 << 29) - 1;
// a + b
let mut carry: u32 = 0;
for i in 0..9 {
carry = a[i] + b[i] + (carry >> 29);
sum[i] = carry & mask;
}
// subtract l if the sum is >= l
Scalar32::sub(&sum, &constants::L)
}
/// Compute `a - b` (mod l).
pub fn sub(a: &Scalar32, b: &Scalar32) -> Scalar32 {
let mut difference = Scalar32::zero();
let mask = (1u32 << 29) - 1;
// a - b
let mut borrow: u32 = 0;
for i in 0..9 {
borrow = a[i].wrapping_sub(b[i] + (borrow >> 31));
difference[i] = borrow & mask;
}
// conditionally add l if the difference is negative
let underflow_mask = ((borrow >> 31) ^ 1).wrapping_sub(1);
let mut carry: u32 = 0;
for i in 0..9 {
carry = (carry >> 29) + difference[i] + (constants::L[i] & underflow_mask);
difference[i] = carry & mask;
}
difference
}
/// Compute `a * b`.
///
/// This is implemented with a one-level refined Karatsuba decomposition
#[inline(always)]
fn mul_internal(a: &Scalar32, b: &Scalar32) -> [u64; 17] {
let mut z = [0u64; 17];
z[0] = m(a[0],b[0]); // c00
z[1] = m(a[0],b[1]) + m(a[1],b[0]); // c01
z[2] = m(a[0],b[2]) + m(a[1],b[1]) + m(a[2],b[0]); // c02
z[3] = m(a[0],b[3]) + m(a[1],b[2]) + m(a[2],b[1]) + m(a[3],b[0]); // c03
z[4] = m(a[0],b[4]) + m(a[1],b[3]) + m(a[2],b[2]) + m(a[3],b[1]) + m(a[4],b[0]); // c04
z[5] = m(a[1],b[4]) + m(a[2],b[3]) + m(a[3],b[2]) + m(a[4],b[1]); // c05
z[6] = m(a[2],b[4]) + m(a[3],b[3]) + m(a[4],b[2]); // c06
z[7] = m(a[3],b[4]) + m(a[4],b[3]); // c07
z[8] = (m(a[4],b[4])).wrapping_sub(z[3]); // c08 - c03
z[10] = z[5].wrapping_sub(m(a[5],b[5])); // c05mc10
z[11] = z[6].wrapping_sub(m(a[5],b[6]) + m(a[6],b[5])); // c06mc11
z[12] = z[7].wrapping_sub(m(a[5],b[7]) + m(a[6],b[6]) + m(a[7],b[5])); // c07mc12
z[13] = m(a[5],b[8]) + m(a[6],b[7]) + m(a[7],b[6]) + m(a[8],b[5]); // c13
z[14] = m(a[6],b[8]) + m(a[7],b[7]) + m(a[8],b[6]); // c14
z[15] = m(a[7],b[8]) + m(a[8],b[7]); // c15
z[16] = m(a[8],b[8]); // c16
z[ 5] = z[10].wrapping_sub(z[ 0]); // c05mc10 - c00
z[ 6] = z[11].wrapping_sub(z[ 1]); // c06mc11 - c01
z[ 7] = z[12].wrapping_sub(z[ 2]); // c07mc12 - c02
z[ 8] = z[ 8].wrapping_sub(z[13]); // c08mc13 - c03
z[ 9] = z[14].wrapping_add(z[ 4]); // c14 + c04
z[10] = z[15].wrapping_add(z[10]); // c15 + c05mc10
z[11] = z[16].wrapping_add(z[11]); // c16 + c06mc11
let aa = [
a[0]+a[5],
a[1]+a[6],
a[2]+a[7],
a[3]+a[8]
];
let bb = [
b[0]+b[5],
b[1]+b[6],
b[2]+b[7],
b[3]+b[8]
];
z[ 5] = (m(aa[0],bb[0])) .wrapping_add(z[ 5]); // c20 + c05mc10 - c00
z[ 6] = (m(aa[0],bb[1]) + m(aa[1],bb[0])) .wrapping_add(z[ 6]); // c21 + c06mc11 - c01
z[ 7] = (m(aa[0],bb[2]) + m(aa[1],bb[1]) + m(aa[2],bb[0])) .wrapping_add(z[ 7]); // c22 + c07mc12 - c02
z[ 8] = (m(aa[0],bb[3]) + m(aa[1],bb[2]) + m(aa[2],bb[1]) + m(aa[3],bb[0])) .wrapping_add(z[ 8]); // c23 + c08mc13 - c03
z[ 9] = (m(aa[0], b[4]) + m(aa[1],bb[3]) + m(aa[2],bb[2]) + m(aa[3],bb[1]) + m(a[4],bb[0])).wrapping_sub(z[ 9]); // c24 - c14 - c04
z[10] = ( m(aa[1], b[4]) + m(aa[2],bb[3]) + m(aa[3],bb[2]) + m(a[4],bb[1])).wrapping_sub(z[10]); // c25 - c15 - c05mc10
z[11] = ( m(aa[2], b[4]) + m(aa[3],bb[3]) + m(a[4],bb[2])).wrapping_sub(z[11]); // c26 - c16 - c06mc11
z[12] = ( m(aa[3], b[4]) + m(a[4],bb[3])).wrapping_sub(z[12]); // c27 - c07mc12
z
}
/// Compute `a^2`.
#[inline(always)]
fn square_internal(a: &Scalar32) -> [u64; 17] {
let aa = [
a[0]*2,
a[1]*2,
a[2]*2,
a[3]*2,
a[4]*2,
a[5]*2,
a[6]*2,
a[7]*2
];
[
m( a[0],a[0]),
m(aa[0],a[1]),
m(aa[0],a[2]) + m( a[1],a[1]),
m(aa[0],a[3]) + m(aa[1],a[2]),
m(aa[0],a[4]) + m(aa[1],a[3]) + m( a[2],a[2]),
m(aa[0],a[5]) + m(aa[1],a[4]) + m(aa[2],a[3]),
m(aa[0],a[6]) + m(aa[1],a[5]) + m(aa[2],a[4]) + m( a[3],a[3]),
m(aa[0],a[7]) + m(aa[1],a[6]) + m(aa[2],a[5]) + m(aa[3],a[4]),
m(aa[0],a[8]) + m(aa[1],a[7]) + m(aa[2],a[6]) + m(aa[3],a[5]) + m( a[4],a[4]),
m(aa[1],a[8]) + m(aa[2],a[7]) + m(aa[3],a[6]) + m(aa[4],a[5]),
m(aa[2],a[8]) + m(aa[3],a[7]) + m(aa[4],a[6]) + m( a[5],a[5]),
m(aa[3],a[8]) + m(aa[4],a[7]) + m(aa[5],a[6]),
m(aa[4],a[8]) + m(aa[5],a[7]) + m( a[6],a[6]),
m(aa[5],a[8]) + m(aa[6],a[7]),
m(aa[6],a[8]) + m( a[7],a[7]),
m(aa[7],a[8]),
m( a[8],a[8]),
]
}
/// Compute `limbs/R` (mod l), where R is the Montgomery modulus 2^261
#[inline(always)]
fn montgomery_reduce(limbs: &[u64; 17]) -> Scalar32 {
#[inline(always)]
fn part1(sum: u64) -> (u64, u32) {
let p = (sum as u32).wrapping_mul(constants::LFACTOR) & ((1u32 << 29) - 1);
((sum + m(p,constants::L[0])) >> 29, p)
}
#[inline(always)]
fn part2(sum: u64) -> (u64, u32) {
let w = (sum as u32) & ((1u32 << 29) - 1);
(sum >> 29, w)
}
// note: l5,l6,l7 are zero, so their multiplies can be skipped
let l = &constants::L;
// the first half computes the Montgomery adjustment factor n, and begins adding n*l to make limbs divisible by R
let (carry, n0) = part1( limbs[ 0]);
let (carry, n1) = part1(carry + limbs[ 1] + m(n0,l[1]));
let (carry, n2) = part1(carry + limbs[ 2] + m(n0,l[2]) + m(n1,l[1]));
let (carry, n3) = part1(carry + limbs[ 3] + m(n0,l[3]) + m(n1,l[2]) + m(n2,l[1]));
let (carry, n4) = part1(carry + limbs[ 4] + m(n0,l[4]) + m(n1,l[3]) + m(n2,l[2]) + m(n3,l[1]));
let (carry, n5) = part1(carry + limbs[ 5] + m(n1,l[4]) + m(n2,l[3]) + m(n3,l[2]) + m(n4,l[1]));
let (carry, n6) = part1(carry + limbs[ 6] + m(n2,l[4]) + m(n3,l[3]) + m(n4,l[2]) + m(n5,l[1]));
let (carry, n7) = part1(carry + limbs[ 7] + m(n3,l[4]) + m(n4,l[3]) + m(n5,l[2]) + m(n6,l[1]));
let (carry, n8) = part1(carry + limbs[ 8] + m(n0,l[8]) + m(n4,l[4]) + m(n5,l[3]) + m(n6,l[2]) + m(n7,l[1]));
// limbs is divisible by R now, so we can divide by R by simply storing the upper half as the result
let (carry, r0) = part2(carry + limbs[ 9] + m(n1,l[8]) + m(n5,l[4]) + m(n6,l[3]) + m(n7,l[2]) + m(n8,l[1]));
let (carry, r1) = part2(carry + limbs[10] + m(n2,l[8]) + m(n6,l[4]) + m(n7,l[3]) + m(n8,l[2]));
let (carry, r2) = part2(carry + limbs[11] + m(n3,l[8]) + m(n7,l[4]) + m(n8,l[3]));
let (carry, r3) = part2(carry + limbs[12] + m(n4,l[8]) + m(n8,l[4]));
let (carry, r4) = part2(carry + limbs[13] + m(n5,l[8]) );
let (carry, r5) = part2(carry + limbs[14] + m(n6,l[8]) );
let (carry, r6) = part2(carry + limbs[15] + m(n7,l[8]) );
let (carry, r7) = part2(carry + limbs[16] + m(n8,l[8]));
let r8 = carry as u32;
// result may be >= l, so attempt to subtract l
Scalar32::sub(&Scalar32([r0,r1,r2,r3,r4,r5,r6,r7,r8]), l)
}
/// Compute `a * b` (mod l).
#[inline(never)]
pub fn mul(a: &Scalar32, b: &Scalar32) -> Scalar32 {
let ab = Scalar32::montgomery_reduce(&Scalar32::mul_internal(a, b));
Scalar32::montgomery_reduce(&Scalar32::mul_internal(&ab, &constants::RR))
}
/// Compute `a^2` (mod l).
#[inline(never)]
pub fn square(&self) -> Scalar32 {
let aa = Scalar32::montgomery_reduce(&Scalar32::square_internal(self));
Scalar32::montgomery_reduce(&Scalar32::mul_internal(&aa, &constants::RR))
}
/// Compute `(a * b) / R` (mod l), where R is the Montgomery modulus 2^261
#[inline(never)]
pub fn montgomery_mul(a: &Scalar32, b: &Scalar32) -> Scalar32 {
Scalar32::montgomery_reduce(&Scalar32::mul_internal(a, b))
}
/// Compute `(a^2) / R` (mod l) in Montgomery form, where R is the Montgomery modulus 2^261
#[inline(never)]
pub fn montgomery_square(&self) -> Scalar32 {
Scalar32::montgomery_reduce(&Scalar32::square_internal(self))
}
/// Puts a Scalar32 in to Montgomery form, i.e. computes `a*R (mod l)`
#[inline(never)]
pub fn to_montgomery(&self) -> Scalar32 {
Scalar32::montgomery_mul(self, &constants::RR)
}
/// Takes a Scalar32 out of Montgomery form, i.e. computes `a/R (mod l)`
pub fn from_montgomery(&self) -> Scalar32 {
let mut limbs = [0u64; 17];
for i in 0..9 {
limbs[i] = self[i] as u64;
}
Scalar32::montgomery_reduce(&limbs)
}
}
#[cfg(test)]
mod test {
use super::*;
/// Note: x is 2^253-1 which is slightly larger than the largest scalar produced by
/// this implementation (l-1), and should verify there are no overflows for valid scalars
///
/// x = 2^253-1 = 14474011154664524427946373126085988481658748083205070504932198000989141204991
/// x = 7237005577332262213973186563042994240801631723825162898930247062703686954002 mod l
/// x = 5147078182513738803124273553712992179887200054963030844803268920753008712037*R mod l in Montgomery form
pub static X: Scalar32 = Scalar32(
[0x1fffffff, 0x1fffffff, 0x1fffffff, 0x1fffffff,
0x1fffffff, 0x1fffffff, 0x1fffffff, 0x1fffffff,
0x001fffff]);
/// x^2 = 3078544782642840487852506753550082162405942681916160040940637093560259278169 mod l
pub static XX: Scalar32 = Scalar32(
[0x00217559, 0x000b3401, 0x103ff43b, 0x1462a62c,
0x1d6f9f38, 0x18e7a42f, 0x09a3dcee, 0x008dbe18,
0x0006ce65]);
/// x^2 = 2912514428060642753613814151688322857484807845836623976981729207238463947987*R mod l in Montgomery form
pub static XX_MONT: Scalar32 = Scalar32(
[0x152b4d2e, 0x0571d53b, 0x1da6d964, 0x188663b6,
0x1d1b5f92, 0x19d50e3f, 0x12306c29, 0x0c6f26fe,
0x00030edb]);
/// y = 6145104759870991071742105800796537629880401874866217824609283457819451087098
pub static Y: Scalar32 = Scalar32(
[0x1e1458fa, 0x165ba838, 0x1d787b36, 0x0e577f3a,
0x1d2baf06, 0x1d689a19, 0x1fff3047, 0x117704ab,
0x000d9601]);
/// x*y = 36752150652102274958925982391442301741
pub static XY: Scalar32 = Scalar32(
[0x0ba7632d, 0x017736bb, 0x15c76138, 0x0c69daa1,
0x000001ba, 0x00000000, 0x00000000, 0x00000000,
0x00000000]);
/// x*y = 3783114862749659543382438697751927473898937741870308063443170013240655651591*R mod l in Montgomery form
pub static XY_MONT: Scalar32 = Scalar32(
[0x077b51e1, 0x1c64e119, 0x02a19ef5, 0x18d2129e,
0x00de0430, 0x045a7bc8, 0x04cfc7c9, 0x1c002681,
0x000bdc1c]);
/// a = 2351415481556538453565687241199399922945659411799870114962672658845158063753
pub static A: Scalar32 = Scalar32(
[0x07b3be89, 0x02291b60, 0x14a99f03, 0x07dc3787,
0x0a782aae, 0x16262525, 0x0cfdb93f, 0x13f5718d,
0x000532da]);
/// b = 4885590095775723760407499321843594317911456947580037491039278279440296187236
pub static B: Scalar32 = Scalar32(
[0x15421564, 0x1e69fd72, 0x093d9692, 0x161785be,
0x1587d69f, 0x09d9dada, 0x130246c0, 0x0c0a8e72,
0x000acd25]);
/// a+b = 0
/// a-b = 4702830963113076907131374482398799845891318823599740229925345317690316127506
pub static AB: Scalar32 = Scalar32(
[0x0f677d12, 0x045236c0, 0x09533e06, 0x0fb86f0f,
0x14f0555c, 0x0c4c4a4a, 0x19fb727f, 0x07eae31a,
0x000a65b5]);
// c = (2^512 - 1) % l = 1627715501170711445284395025044413883736156588369414752970002579683115011840
pub static C: Scalar32 = Scalar32(
[0x049c0f00, 0x00308f1a, 0x0164d1e9, 0x1c374ed1,
0x1be65d00, 0x19e90bfa, 0x08f73bb1, 0x036f8613,
0x00039941]);
#[test]
fn mul_max() {
let res = Scalar32::mul(&X, &X);
for i in 0..9 {
assert!(res[i] == XX[i]);
}
}
#[test]
fn square_max() {
let res = X.square();
for i in 0..9 {
assert!(res[i] == XX[i]);
}
}
#[test]
fn montgomery_mul_max() {
let res = Scalar32::montgomery_mul(&X, &X);
for i in 0..9 {
assert!(res[i] == XX_MONT[i]);
}
}
#[test]
fn montgomery_square_max() {
let res = X.montgomery_square();
for i in 0..9 {
assert!(res[i] == XX_MONT[i]);
}
}
#[test]
fn mul() {
let res = Scalar32::mul(&X, &Y);
for i in 0..9 {
assert!(res[i] == XY[i]);
}
}
#[test]
fn montgomery_mul() {
let res = Scalar32::montgomery_mul(&X, &Y);
for i in 0..9 {
assert!(res[i] == XY_MONT[i]);
}
}
#[test]
fn add() {
let res = Scalar32::add(&A, &B);
let zero = Scalar32::zero();
for i in 0..9 {
assert!(res[i] == zero[i]);
}
}
#[test]
fn sub() {
let res = Scalar32::sub(&A, &B);
for i in 0..9 {
assert!(res[i] == AB[i]);
}
}
#[test]
fn from_bytes_wide() {
let bignum = [255u8; 64]; // 2^512 - 1
let reduced = Scalar32::from_bytes_wide(&bignum);
for i in 0..9 {
assert!(reduced[i] == C[i]);
}
}
}
#[cfg(all(test, feature = "bench"))]
mod bench {
use test::Bencher;
use super::*;
use super::test::{X, Y};
#[bench]
fn square(b: &mut Bencher) {
b.iter(|| X.square());
}
#[bench]
fn mul(b: &mut Bencher) {
b.iter(|| Scalar32::mul(&X, &Y));
}
#[bench]
fn montgomery_square(b: &mut Bencher) {
b.iter(|| X.montgomery_square());
}
#[bench]
fn montgomery_mul(b: &mut Bencher) {
b.iter(|| Scalar32::montgomery_mul(&X, &Y));
}
#[bench]
fn from_bytes_wide(b: &mut Bencher) {
let bignum = [255u8; 64]; // 2^512 - 1
b.iter(|| Scalar32::from_bytes_wide(&bignum));
}
}

475
src/scalar_64bit.rs Normal file
View file

@ -0,0 +1,475 @@
//! Arithmetic mod 2^252 + 27742317777372353535851937790883648493
//! with 5 52-bit unsigned limbs. 51-bit limbs would cover the
//! desired bit range (253 bits), but isn't large enough to reduce
//! a 512 bit number with Montgomery multiplication, so 52 bits is
//! used instead
//!
//! To see that this is safe for intermediate results, note that
//! the largest limb in a 5 by 5 product of 52-bit limbs will be
//! (0xfffffffffffff^2) * 5 = 0x4ffffffffffff60000000000005 (107 bits).
//!
//! (the 5th limb will never exceed 45 bits, so the actual
//! ranges are slightly smaller)
use core::fmt::Debug;
use core::ops::{Index, IndexMut};
use constants;
/// The `Scalar64` struct represents an element in ℤ/lℤ as 5 52-bit limbs
#[derive(Copy,Clone)]
pub struct Scalar64(pub [u64; 5]);
impl Debug for Scalar64 {
fn fmt(&self, f: &mut ::core::fmt::Formatter) -> ::core::fmt::Result {
write!(f, "Scalar64: {:?}", &self.0[..])
}
}
impl Index<usize> for Scalar64 {
type Output = u64;
fn index(&self, _index: usize) -> &u64 {
&(self.0[_index])
}
}
impl IndexMut<usize> for Scalar64 {
fn index_mut(&mut self, _index: usize) -> &mut u64 {
&mut (self.0[_index])
}
}
/// u64 * u64 = u128 multiply helper
#[inline(always)]
fn m(x: u64, y: u64) -> u128 {
(x as u128) * (y as u128)
}
impl Scalar64 {
/// Return the zero scalar
pub fn zero() -> Scalar64 {
Scalar64([0,0,0,0,0])
}
/// Unpack a 32 byte / 256 bit scalar into 5 52-bit limbs, ignoring the upper 3 bits
pub fn from_bytes(bytes: &[u8; 32]) -> Scalar64 {
let mut words = [0u64; 8];
for i in 0..4 {
for j in 0..8 {
words[i] |= (bytes[(i * 8) + j] as u64) << (j * 8);
}
}
let mask = (1u64 << 52) - 1;
let top_mask = (1u64 << 45) - 1;
let mut s = Scalar64::zero();
s[ 0] = words[0] & mask;
s[ 1] = ((words[0] >> 52) | (words[1] << 12)) & mask;
s[ 2] = ((words[1] >> 40) | (words[2] << 24)) & mask;
s[ 3] = ((words[2] >> 28) | (words[3] << 36)) & mask;
s[ 4] = (words[3] >> 16) & top_mask;
s
}
/// Reduce a 64 byte / 512 bit scalar mod l
pub fn from_bytes_wide(bytes: &[u8; 64]) -> Scalar64 {
let mut words = [064; 16];
for i in 0..8 {
for j in 0..8 {
words[i] |= (bytes[(i * 8) + j] as u64) << (j * 8);
}
}
let mask = (1u64 << 52) - 1;
let mut lo = Scalar64::zero();
let mut hi = Scalar64::zero();
lo[0] = words[ 0] & mask;
lo[1] = ((words[ 0] >> 52) | (words[ 1] << 12)) & mask;
lo[2] = ((words[ 1] >> 40) | (words[ 2] << 24)) & mask;
lo[3] = ((words[ 2] >> 28) | (words[ 3] << 36)) & mask;
lo[4] = ((words[ 3] >> 16) | (words[ 4] << 48)) & mask;
hi[0] = (words[ 4] >> 4) & mask;
hi[1] = ((words[ 4] >> 56) | (words[ 5] << 8)) & mask;
hi[2] = ((words[ 5] >> 44) | (words[ 6] << 20)) & mask;
hi[3] = ((words[ 6] >> 32) | (words[ 7] << 32)) & mask;
hi[4] = words[ 7] >> 20 ;
lo = Scalar64::montgomery_mul(&lo, &constants::R); // (lo * R) / R = lo
hi = Scalar64::montgomery_mul(&hi, &constants::RR); // (hi * R^2) / R = hi * R
Scalar64::add(&hi, &lo)
}
/// Pack the limbs of this `Scalar64` into 32 bytes
pub fn to_bytes(&self) -> [u8; 32] {
let mut s = [0u8; 32];
s[0] = (self.0[ 0] >> 0) as u8;
s[1] = (self.0[ 0] >> 8) as u8;
s[2] = (self.0[ 0] >> 16) as u8;
s[3] = (self.0[ 0] >> 24) as u8;
s[4] = (self.0[ 0] >> 32) as u8;
s[5] = (self.0[ 0] >> 40) as u8;
s[6] = ((self.0[ 0] >> 48) | (self.0[ 1] << 4)) as u8;
s[7] = (self.0[ 1] >> 4) as u8;
s[8] = (self.0[ 1] >> 12) as u8;
s[9] = (self.0[ 1] >> 20) as u8;
s[10] = (self.0[ 1] >> 28) as u8;
s[11] = (self.0[ 1] >> 36) as u8;
s[12] = (self.0[ 1] >> 44) as u8;
s[13] = (self.0[ 2] >> 0) as u8;
s[14] = (self.0[ 2] >> 8) as u8;
s[15] = (self.0[ 2] >> 16) as u8;
s[16] = (self.0[ 2] >> 24) as u8;
s[17] = (self.0[ 2] >> 32) as u8;
s[18] = (self.0[ 2] >> 40) as u8;
s[19] = ((self.0[ 2] >> 48) | (self.0[ 3] << 4)) as u8;
s[20] = (self.0[ 3] >> 4) as u8;
s[21] = (self.0[ 3] >> 12) as u8;
s[22] = (self.0[ 3] >> 20) as u8;
s[23] = (self.0[ 3] >> 28) as u8;
s[24] = (self.0[ 3] >> 36) as u8;
s[25] = (self.0[ 3] >> 44) as u8;
s[26] = (self.0[ 4] >> 0) as u8;
s[27] = (self.0[ 4] >> 8) as u8;
s[28] = (self.0[ 4] >> 16) as u8;
s[29] = (self.0[ 4] >> 24) as u8;
s[30] = (self.0[ 4] >> 32) as u8;
s[31] = (self.0[ 4] >> 40) as u8;
s
}
/// Compute `a + b` (mod l)
pub fn add(a: &Scalar64, b: &Scalar64) -> Scalar64 {
let mut sum = Scalar64::zero();
let mask = (1u64 << 52) - 1;
// a + b
let mut carry: u64 = 0;
for i in 0..5 {
carry = a[i] + b[i] + (carry >> 52);
sum[i] = carry & mask;
}
// subtract l if the sum is >= l
Scalar64::sub(&sum, &constants::L)
}
/// Compute `a - b` (mod l)
pub fn sub(a: &Scalar64, b: &Scalar64) -> Scalar64 {
let mut difference = Scalar64::zero();
let mask = (1u64 << 52) - 1;
// a - b
let mut borrow: u64 = 0;
for i in 0..5 {
borrow = a[i].wrapping_sub(b[i] + (borrow >> 63));
difference[i] = borrow & mask;
}
// conditionally add l if the difference is negative
let underflow_mask = ((borrow >> 63) ^ 1).wrapping_sub(1);
let mut carry: u64 = 0;
for i in 0..5 {
carry = (carry >> 52) + difference[i] + (constants::L[i] & underflow_mask);
difference[i] = carry & mask;
}
difference
}
/// Compute `a * b`
#[inline(always)]
fn mul_internal(a: &Scalar64, b: &Scalar64) -> [u128; 9] {
[
m(a[0],b[0]),
m(a[0],b[1]) + m(a[1],b[0]),
m(a[0],b[2]) + m(a[1],b[1]) + m(a[2],b[0]),
m(a[0],b[3]) + m(a[1],b[2]) + m(a[2],b[1]) + m(a[3],b[0]),
m(a[0],b[4]) + m(a[1],b[3]) + m(a[2],b[2]) + m(a[3],b[1]) + m(a[4],b[0]),
m(a[1],b[4]) + m(a[2],b[3]) + m(a[3],b[2]) + m(a[4],b[1]),
m(a[2],b[4]) + m(a[3],b[3]) + m(a[4],b[2]),
m(a[3],b[4]) + m(a[4],b[3]),
m(a[4],b[4])
]
}
/// Compute `a^2`
#[inline(always)]
fn square_internal(a: &Scalar64) -> [u128; 9] {
let aa = [
a[0]*2,
a[1]*2,
a[2]*2,
a[3]*2,
];
[
m( a[0],a[0]),
m(aa[0],a[1]),
m(aa[0],a[2]) + m( a[1],a[1]),
m(aa[0],a[3]) + m(aa[1],a[2]),
m(aa[0],a[4]) + m(aa[1],a[3]) + m( a[2],a[2]),
m(aa[1],a[4]) + m(aa[2],a[3]),
m(aa[2],a[4]) + m( a[3],a[3]),
m(aa[3],a[4]),
m(a[4],a[4])
]
}
/// Compute `limbs/R` (mod l), where R is the Montgomery modulus 2^260
#[inline(always)]
fn montgomery_reduce(limbs: &[u128; 9]) -> Scalar64 {
#[inline(always)]
fn part1(sum: u128) -> (u128, u64) {
let p = (sum as u64).wrapping_mul(constants::LFACTOR) & ((1u64 << 52) - 1);
((sum + m(p,constants::L[0])) >> 52, p)
}
#[inline(always)]
fn part2(sum: u128) -> (u128, u64) {
let w = (sum as u64) & ((1u64 << 52) - 1);
(sum >> 52, w)
}
// note: l3 is zero, so its multiplies can be skipped
let l = &constants::L;
// the first half computes the Montgomery adjustment factor n, and begins adding n*l to make limbs divisible by R
let (carry, n0) = part1( limbs[0]);
let (carry, n1) = part1(carry + limbs[1] + m(n0,l[1]));
let (carry, n2) = part1(carry + limbs[2] + m(n0,l[2]) + m(n1,l[1]));
let (carry, n3) = part1(carry + limbs[3] + m(n1,l[2]) + m(n2,l[1]));
let (carry, n4) = part1(carry + limbs[4] + m(n0,l[4]) + m(n2,l[2]) + m(n3,l[1]));
// limbs is divisible by R now, so we can divide by R by simply storing the upper half as the result
let (carry, r0) = part2(carry + limbs[5] + m(n1,l[4]) + m(n3,l[2]) + m(n4,l[1]));
let (carry, r1) = part2(carry + limbs[6] + m(n2,l[4]) + m(n4,l[2]));
let (carry, r2) = part2(carry + limbs[7] + m(n3,l[4]) );
let (carry, r3) = part2(carry + limbs[8] + m(n4,l[4]));
let r4 = carry as u64;
// result may be >= l, so attempt to subtract l
Scalar64::sub(&Scalar64([r0,r1,r2,r3,r4]), l)
}
/// Compute `a * b` (mod l)
#[inline(never)]
pub fn mul(a: &Scalar64, b: &Scalar64) -> Scalar64 {
let ab = Scalar64::montgomery_reduce(&Scalar64::mul_internal(a, b));
Scalar64::montgomery_reduce(&Scalar64::mul_internal(&ab, &constants::RR))
}
/// Compute `a^2` (mod l)
#[inline(never)]
pub fn square(&self) -> Scalar64 {
let aa = Scalar64::montgomery_reduce(&Scalar64::square_internal(self));
Scalar64::montgomery_reduce(&Scalar64::mul_internal(&aa, &constants::RR))
}
/// Compute `(a * b) / R` (mod l), where R is the Montgomery modulus 2^260
#[inline(never)]
pub fn montgomery_mul(a: &Scalar64, b: &Scalar64) -> Scalar64 {
Scalar64::montgomery_reduce(&Scalar64::mul_internal(a, b))
}
/// Compute `(a^2) / R` (mod l) in Montgomery form, where R is the Montgomery modulus 2^260
#[inline(never)]
pub fn montgomery_square(&self) -> Scalar64 {
Scalar64::montgomery_reduce(&Scalar64::square_internal(self))
}
/// Puts a Scalar64 in to Montgomery form, i.e. computes `a*R (mod l)`
#[inline(never)]
pub fn to_montgomery(&self) -> Scalar64 {
Scalar64::montgomery_mul(self, &constants::RR)
}
/// Takes a Scalar64 out of Montgomery form, i.e. computes `a/R (mod l)`
#[inline(never)]
pub fn from_montgomery(&self) -> Scalar64 {
let mut limbs = [0u128; 9];
for i in 0..5 {
limbs[i] = self[i] as u128;
}
Scalar64::montgomery_reduce(&limbs)
}
}
#[cfg(test)]
mod test {
use super::*;
/// Note: x is 2^253-1 which is slightly larger than the largest scalar produced by
/// this implementation (l-1), and should show there are no overflows for valid scalars
///
/// x = 14474011154664524427946373126085988481658748083205070504932198000989141204991
/// x = 7237005577332262213973186563042994240801631723825162898930247062703686954002 mod l
/// x = 3057150787695215392275360544382990118917283750546154083604586903220563173085*R mod l in Montgomery form
pub static X: Scalar64 = Scalar64(
[0x000fffffffffffff, 0x000fffffffffffff, 0x000fffffffffffff, 0x000fffffffffffff,
0x00001fffffffffff]);
/// x^2 = 3078544782642840487852506753550082162405942681916160040940637093560259278169 mod l
pub static XX: Scalar64 = Scalar64(
[0x0001668020217559, 0x000531640ffd0ec0, 0x00085fd6f9f38a31, 0x000c268f73bb1cf4,
0x000006ce65046df0]);
/// x^2 = 4413052134910308800482070043710297189082115023966588301924965890668401540959*R mod l in Montgomery form
pub static XX_MONT: Scalar64 = Scalar64(
[0x000c754eea569a5c, 0x00063b6ed36cb215, 0x0008ffa36bf25886, 0x000e9183614e7543,
0x0000061db6c6f26f]);
/// y = 6145104759870991071742105800796537629880401874866217824609283457819451087098
pub static Y: Scalar64 = Scalar64(
[0x000b75071e1458fa, 0x000bf9d75e1ecdac, 0x000433d2baf0672b, 0x0005fffcc11fad13,
0x00000d96018bb825]);
/// x*y = 36752150652102274958925982391442301741 mod l
pub static XY: Scalar64 = Scalar64(
[0x000ee6d76ba7632d, 0x000ed50d71d84e02, 0x00000000001ba634, 0x0000000000000000,
0x0000000000000000]);
/// x*y = 658448296334113745583381664921721413881518248721417041768778176391714104386*R mod l in Montgomery form
pub static XY_MONT: Scalar64 = Scalar64(
[0x0006d52bf200cfd5, 0x00033fb1d7021570, 0x000f201bc07139d8, 0x0001267e3e49169e,
0x000007b839c00268]);
/// a = 2351415481556538453565687241199399922945659411799870114962672658845158063753
pub static A: Scalar64 = Scalar64(
[0x0005236c07b3be89, 0x0001bc3d2a67c0c4, 0x000a4aa782aae3ee, 0x0006b3f6e4fec4c4,
0x00000532da9fab8c]);
/// b = 4885590095775723760407499321843594317911456947580037491039278279440296187236
pub static B: Scalar64 = Scalar64(
[0x000d3fae55421564, 0x000c2df24f65a4bc, 0x0005b5587d69fb0b, 0x00094c091b013b3b,
0x00000acd25605473]);
/// a+b = 0
/// a-b = 4702830963113076907131374482398799845891318823599740229925345317690316127506
pub static AB: Scalar64 = Scalar64(
[0x000a46d80f677d12, 0x0003787a54cf8188, 0x0004954f0555c7dc, 0x000d67edc9fd8989,
0x00000a65b53f5718]);
// c = (2^512 - 1) % l = 1627715501170711445284395025044413883736156588369414752970002579683115011840
pub static C: Scalar64 = Scalar64(
[0x000611e3449c0f00, 0x000a768859347a40, 0x0007f5be65d00e1b, 0x0009a3dceec73d21,
0x00000399411b7c30]);
#[test]
fn mul_max() {
let res = Scalar64::mul(&X, &X);
for i in 0..5 {
assert!(res[i] == XX[i]);
}
}
#[test]
fn square_max() {
let res = X.square();
for i in 0..5 {
assert!(res[i] == XX[i]);
}
}
#[test]
fn montgomery_mul_max() {
let res = Scalar64::montgomery_mul(&X, &X);
for i in 0..5 {
assert!(res[i] == XX_MONT[i]);
}
}
#[test]
fn montgomery_square_max() {
let res = X.montgomery_square();
for i in 0..5 {
assert!(res[i] == XX_MONT[i]);
}
}
#[test]
fn mul() {
let res = Scalar64::mul(&X, &Y);
for i in 0..5 {
assert!(res[i] == XY[i]);
}
}
#[test]
fn montgomery_mul() {
let res = Scalar64::montgomery_mul(&X, &Y);
for i in 0..5 {
assert!(res[i] == XY_MONT[i]);
}
}
#[test]
fn add() {
let res = Scalar64::add(&A, &B);
let zero = Scalar64::zero();
for i in 0..5 {
assert!(res[i] == zero[i]);
}
}
#[test]
fn sub() {
let res = Scalar64::sub(&A, &B);
for i in 0..5 {
assert!(res[i] == AB[i]);
}
}
#[test]
fn from_bytes_wide() {
let bignum = [255u8; 64]; // 2^512 - 1
let reduced = Scalar64::from_bytes_wide(&bignum);
println!("{:?}", reduced);
for i in 0..5 {
assert!(reduced[i] == C[i]);
}
}
}
#[cfg(all(test, feature = "bench"))]
mod bench {
use test::Bencher;
use super::*;
use super::test::{X, Y};
#[bench]
fn square(b: &mut Bencher) {
b.iter(|| X.square());
}
#[bench]
fn mul(b: &mut Bencher) {
b.iter(|| Scalar64::mul(&X, &Y));
}
#[bench]
fn montgomery_square(b: &mut Bencher) {
b.iter(|| X.montgomery_square());
}
#[bench]
fn montgomery_mul(b: &mut Bencher) {
b.iter(|| Scalar64::montgomery_mul(&X, &Y));
}
#[bench]
fn from_bytes_wide(b: &mut Bencher) {
let bignum = [255u8; 64]; // 2^512 - 1
b.iter(|| Scalar64::from_bytes_wide(&bignum));
}
}