diff --git a/Cargo.toml b/Cargo.toml index 7b901d7..d634479 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "curve25519-dalek" -version = "0.12.1" +version = "0.13.0" authors = ["Isis Lovecruft ", "Henry de Valence "] readme = "README.md" @@ -15,6 +15,9 @@ exclude = [ ".gitignore" ] +[package.metadata.docs.rs] +rustdoc-args = ["--html-in-header rustdoc-include-katex-header.html"] + [badges] travis-ci = { repository = "isislovecruft/curve25519-dalek", branch = "master"} @@ -51,59 +54,8 @@ nightly = ["radix_51", "subtle/nightly"] default = ["std"] std = ["rand", "subtle/std"] alloc = [] +# This isn't used at the moment, but keep it around for future yolocrypto features. yolocrypto = [] bench = [] # Radix-51 arithmetic using u128 radix_51 = [] - -# The development profile, used for `cargo build`. -[profile.dev] -opt-level = 0 # controls the `--opt-level` the compiler builds with -debug = true # controls whether the compiler passes `-g` -rpath = false # controls whether the compiler passes `-C rpath` -lto = false # controls `-C lto` for binaries and staticlibs -debug-assertions = true # controls whether debug assertions are enabled -codegen-units = 1 # controls whether the compiler passes `-C codegen-units` - # `codegen-units` is ignored when `lto = true` -panic = 'unwind' # panic strategy (`-C panic=...`), can also be 'abort' - -# The release profile, used for `cargo build --release`. -[profile.release] -opt-level = 3 -debug = false -rpath = false -lto = false -debug-assertions = false -codegen-units = 1 -panic = 'unwind' - -# The testing profile, used for `cargo test`. -[profile.test] -opt-level = 0 -debug = true -rpath = false -lto = false -debug-assertions = true -codegen-units = 1 -panic = 'unwind' -required-features = ['yolocrypto'] - -# The benchmarking profile, used for `cargo bench`. -[profile.bench] -opt-level = 3 -debug = false -rpath = false -lto = false -debug-assertions = false -codegen-units = 1 -panic = 'unwind' - -# The documentation profile, used for `cargo doc`. -[profile.doc] -opt-level = 0 -debug = true -rpath = false -lto = false -debug-assertions = true -codegen-units = 1 -panic = 'unwind' diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..eb4f935 --- /dev/null +++ b/Makefile @@ -0,0 +1,3 @@ + +doc: + cargo rustdoc --features "nightly yolocrypto" -- --html-in-header katex-header.html diff --git a/README.md b/README.md index 77420cd..0e27d0d 100644 --- a/README.md +++ b/README.md @@ -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.12" + curve25519-dalek = "^0.13" Then, in your library or executable source, add: @@ -60,3 +60,11 @@ We intend to stabilise the following before curve25519-dalek-1.0.0: * Implement hashing to a point on the curve (Elligator). * Finish Ristretto (Decaf for curve25519) implementation. + +## Contributing + +Please see +[CONTRIBUTING.md](https://github.com/isislovecruft/curve25519-dalek/blob/master/CONTRIBUTING.md). + +Patches and pull requests should be make against the `develop` +branch, **not** `master`. diff --git a/fuzz/Cargo.toml b/fuzz/Cargo.toml index 5245cbe..ceb067b 100644 --- a/fuzz/Cargo.toml +++ b/fuzz/Cargo.toml @@ -10,7 +10,11 @@ cargo-fuzz = true [dependencies.curve25519-dalek] path = ".." -features = ["yolocrypto"] + +[features] +yolocrypto = ["curve25519-dalek/yolocrypto"] +nightly = ["curve25519-dalek/nightly"] +radix_51 = ["curve25519-dalek/radix_51"] [dependencies.libfuzzer-sys] git = "https://github.com/rust-fuzz/libfuzzer-sys.git" @@ -21,4 +25,8 @@ members = ["."] [[bin]] name = "decaf" -path = "fuzzers/decaf.rs" +path = "fuzz_targets/decaf.rs" + +[[bin]] +name = "scalar_constructor_accepts_256bit_values" +path = "fuzz_targets/scalar_constructor_accepts_256bit_values.rs" diff --git a/fuzz/fuzzers/decaf.rs b/fuzz/fuzz_targets/decaf.rs similarity index 100% rename from fuzz/fuzzers/decaf.rs rename to fuzz/fuzz_targets/decaf.rs diff --git a/fuzz/fuzz_targets/scalar_constructor_accepts_256bit_values.rs b/fuzz/fuzz_targets/scalar_constructor_accepts_256bit_values.rs new file mode 100644 index 0000000..38bba88 --- /dev/null +++ b/fuzz/fuzz_targets/scalar_constructor_accepts_256bit_values.rs @@ -0,0 +1,36 @@ +#![no_main] +#[macro_use] extern crate libfuzzer_sys; +extern crate curve25519_dalek; + +use curve25519_dalek::scalar::Scalar; + +/// Check that the Scalar constructor accepts 256-bit input values and +/// behaves correctly on them. +/// +/// Specifically, we take 256-bit values `a` and `b` from the fuzzer +/// input data and check that `(a mod l) * (b mod l) == (a * b) mod l`. +fuzz_target!(|data: &[u8]| { + if data.len() != 64 { + return; + } + let mut a_bytes = [0u8; 32]; + let mut b_bytes = [0u8; 32]; + + // Set a, b to be random 256-bit integers + a_bytes.copy_from_slice(&data[ 0..32]); + b_bytes.copy_from_slice(&data[32..64]); + + // Compute c = a*b (mod l) + let c1 = &Scalar(a_bytes) * &Scalar(b_bytes); + + // Compute c = (a mod l) * (b mod l) + let mut tmp = [0u8; 64]; + tmp[0..32].copy_from_slice(&a_bytes[..]); + let a_mod_l = Scalar::reduce(&tmp); + tmp[0..32].copy_from_slice(&b_bytes[..]); + let b_mod_l = Scalar::reduce(&tmp); + + let c2 = &a_mod_l * &b_mod_l; + + assert_eq!(c1, c2); +}); diff --git a/rustdoc-include-katex-header.html b/rustdoc-include-katex-header.html new file mode 100644 index 0000000..695c735 --- /dev/null +++ b/rustdoc-include-katex-header.html @@ -0,0 +1,10 @@ + + + + + diff --git a/src/constants.rs b/src/constants.rs index 4c08102..63648e5 100644 --- a/src/constants.rs +++ b/src/constants.rs @@ -11,16 +11,26 @@ //! This module contains various constants (such as curve parameters //! and useful field elements like `sqrt(-1)`), as well as //! lookup tables of pre-computed points. +//! +//! Most of the constants are given with +//! `LONG_DESCRIPTIVE_UPPER_CASE_NAMES`, but they can be brought into +//! scope using a `let` binding: +//! +//! ``` +//! use curve25519_dalek::constants; +//! use curve25519_dalek::edwards::IsIdentity; +//! +//! let B = &constants::RISTRETTO_BASEPOINT_TABLE; +//! let l = &constants::BASEPOINT_ORDER; +//! +//! let A = l * B; +//! assert!(A.is_identity()); +//! ``` -#![allow(dead_code)] -#![allow(non_snake_case)] -#![allow(non_upper_case_globals)] -#![allow(missing_docs)] #![allow(non_snake_case)] use edwards::CompressedEdwardsY; -#[cfg(feature = "yolocrypto")] -use decaf::{DecafPoint, DecafBasepointTable}; +use ristretto::{RistrettoPoint, RistrettoBasepointTable}; use montgomery::CompressedMontgomeryU; use scalar::Scalar; @@ -29,20 +39,6 @@ pub use constants_64bit::*; #[cfg(not(feature="radix_51"))] pub use constants_32bit::*; -/// (p-1)/2, in little-endian bytes. -pub const HALF_P_MINUS_1_BYTES: [u8; 32] = - [0xf6, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, - 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, - 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, - 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x3f]; - -/// `HALF_Q_MINUS_1_BYTES` is (2^255-20)/2 expressed in little endian form. -pub const HALF_Q_MINUS_1_BYTES: [u8; 32] = [ // halfQMinus1Bytes - 0xf6, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, - 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, - 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, - 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x3f, ]; - /// Basepoint has y = 4/5. /// /// Generated with Sage: these are the bytes of 4/5 in 𝔽_p. The @@ -61,36 +57,38 @@ pub const BASE_COMPRESSED_MONTGOMERY: CompressedMontgomeryU = 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00]); -/// The Ed25519 basepoint, as a `DecafPoint`. This is called `_POINT` to distinguish it from +/// The Ed25519 basepoint, as a `RistrettoPoint`. This is called `_POINT` to distinguish it from /// `_TABLE`, which provides fast scalar multiplication. -#[cfg(feature = "yolocrypto")] pub const DECAF_ED25519_BASEPOINT_POINT: DecafPoint = -DecafPoint(ED25519_BASEPOINT_POINT); +pub const RISTRETTO_BASEPOINT_POINT: RistrettoPoint = RistrettoPoint(ED25519_BASEPOINT_POINT); -/// `l` is the order of base point, i.e. 2^252 + -/// 27742317777372353535851937790883648493, in little-endian form -pub const l: Scalar = Scalar([ 0xed, 0xd3, 0xf5, 0x5c, 0x1a, 0x63, 0x12, 0x58, - 0xd6, 0x9c, 0xf7, 0xa2, 0xde, 0xf9, 0xde, 0x14, - 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, - 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x10 ]); +/// `BASEPOINT_ORDER` is the order of base point, i.e. `l = 2^252 + +/// 27742317777372353535851937790883648493`, in little-endian bytes. +pub const BASEPOINT_ORDER: Scalar = Scalar([ + 0xed, 0xd3, 0xf5, 0x5c, 0x1a, 0x63, 0x12, 0x58, + 0xd6, 0x9c, 0xf7, 0xa2, 0xde, 0xf9, 0xde, 0x14, + 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, + 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x10, +]); -/// `l_minus_1` is the order of base point minus one, i.e. 2^252 + -/// 27742317777372353535851937790883648493 - 1, in little-endian form -pub const l_minus_1: Scalar = Scalar([ 0xec, 0xd3, 0xf5, 0x5c, 0x1a, 0x63, 0x12, 0x58, - 0xd6, 0x9c, 0xf7, 0xa2, 0xde, 0xf9, 0xde, 0x14, - 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, - 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x10 ]); +/// `BASEPOINT_ORDER_MINUS_1` is the order of base point minus one, i.e. `l-1`, in little-endian bytes. +pub const BASEPOINT_ORDER_MINUS_1: Scalar = Scalar([ + 0xec, 0xd3, 0xf5, 0x5c, 0x1a, 0x63, 0x12, 0x58, + 0xd6, 0x9c, 0xf7, 0xa2, 0xde, 0xf9, 0xde, 0x14, + 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, + 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x10, +]); -/// `lminus1` is the order of base point minus two, i.e. 2^252 + -/// 27742317777372353535851937790883648493 - 2, in little-endian form -pub const l_minus_2: Scalar = Scalar([ 0xeb, 0xd3, 0xf5, 0x5c, 0x1a, 0x63, 0x12, 0x58, - 0xd6, 0x9c, 0xf7, 0xa2, 0xde, 0xf9, 0xde, 0x14, - 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, - 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x10 ]); +/// `BASEPOINT_ORDER_MINUS_2` is the order of base point minus two, i.e. `l-2`, in little-endian bytes. +pub const BASEPOINT_ORDER_MINUS_2: Scalar = Scalar([ + 0xeb, 0xd3, 0xf5, 0x5c, 0x1a, 0x63, 0x12, 0x58, + 0xd6, 0x9c, 0xf7, 0xa2, 0xde, 0xf9, 0xde, 0x14, + 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, + 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x10, +]); -#[cfg(feature = "yolocrypto")] -/// The Ed25519 basepoint -pub const DECAF_ED25519_BASEPOINT_TABLE: DecafBasepointTable - = DecafBasepointTable(ED25519_BASEPOINT_TABLE); +/// The Ed25519 basepoint, as a RistrettoPoint +pub const RISTRETTO_BASEPOINT_TABLE: RistrettoBasepointTable + = RistrettoBasepointTable(ED25519_BASEPOINT_TABLE); #[cfg(test)] mod test { @@ -126,13 +124,6 @@ mod test { } } - #[test] - fn test_half() { - let one = FieldElement::one(); - let two = &one + &one; - assert_eq!(one, &two * &constants::HALF); - } - /// Test that the constant for sqrt(-486664) really is a square /// root of -486664. #[test] @@ -158,27 +149,22 @@ mod test { } #[test] - /// Test that SQRT_M1 and MSQRT_M1 are square roots of -1 + /// Test that SQRT_M1 is a square root of -1 fn test_sqrt_minus_one() { let minus_one = FieldElement::minus_one(); let sqrt_m1_sq = &constants::SQRT_M1 * &constants::SQRT_M1; - let msqrt_m1_sq = &constants::MSQRT_M1 * &constants::MSQRT_M1; assert_eq!(minus_one, sqrt_m1_sq); - assert_eq!(minus_one, msqrt_m1_sq); } #[test] fn test_sqrt_constants_sign() { - let one = FieldElement::one(); let minus_one = FieldElement::minus_one(); let (was_nonzero_square, invsqrt_m1) = minus_one.invsqrt(); assert_eq!(was_nonzero_square, 1u8); let sign_test_sqrt = &invsqrt_m1 * &constants::SQRT_M1; - let sign_test_msqrt = &invsqrt_m1 * &constants::MSQRT_M1; // XXX it seems we have flipped the sign relative to // the invsqrt function? assert_eq!(sign_test_sqrt, minus_one); - assert_eq!(sign_test_msqrt, one); } /// Test that d = -121665/121666 @@ -190,8 +176,8 @@ mod test { let b = FieldElement32([ 121666,0,0,0,0,0,0,0,0,0]); let d = &a * &b.invert(); let d2 = &d + &d; - assert_eq!(d, constants::d); - assert_eq!(d2, constants::d2); + assert_eq!(d, constants::EDWARDS_D); + assert_eq!(d2, constants::EDWARDS_D2); } /// Test that d = -121665/121666 @@ -203,27 +189,16 @@ mod test { let b = FieldElement64([121666,0,0,0,0]); let d = &a * &b.invert(); let d2 = &d + &d; - assert_eq!(d, constants::d); - assert_eq!(d2, constants::d2); + assert_eq!(d, constants::EDWARDS_D); + assert_eq!(d2, constants::EDWARDS_D2); } #[test] - fn test_d4() { - let mut four = FieldElement::zero(); - // XXX should have a way to create small field elements - four.0[0] = 4; - assert_eq!(&constants::d * &four, constants::d4); - } - - #[test] - fn test_a_minus_d() { + fn test_sqrt_ad_minus_one() { let a = FieldElement::minus_one(); - let a_minus_d = &a - &constants::d; - assert_eq!(a_minus_d, constants::a_minus_d); - let (_, invsqrt_a_minus_d) = constants::a_minus_d.invsqrt(); - assert_eq!(invsqrt_a_minus_d, constants::invsqrt_a_minus_d); - let inv_a_minus_d = invsqrt_a_minus_d.square(); - assert_eq!(inv_a_minus_d, constants::inv_a_minus_d); - assert_eq!(&inv_a_minus_d * &a_minus_d, FieldElement::one()); + let ad_minus_one = &(&a * &constants::EDWARDS_D) + &a; + let should_be_ad_minus_one = constants::SQRT_AD_MINUS_ONE.square(); + assert_eq!(should_be_ad_minus_one, ad_minus_one); } + } diff --git a/src/constants_32bit.rs b/src/constants_32bit.rs index bf46fe5..d5a0dd7 100644 --- a/src/constants_32bit.rs +++ b/src/constants_32bit.rs @@ -12,82 +12,49 @@ //! and useful field elements like `sqrt(-1)`), as well as //! lookup tables of pre-computed points. -#![allow(dead_code)] -#![allow(non_snake_case)] -#![allow(non_upper_case_globals)] -#![allow(missing_docs)] #![allow(non_snake_case)] use field_32bit::FieldElement32; +use scalar_32bit::Scalar32; use edwards::ExtendedPoint; use edwards::AffineNielsPoint; use edwards::EdwardsBasepointTable; -pub const d: FieldElement32 = FieldElement32([ +/// Edwards `d` value, equal to `-121665/121666 mod p`. +pub(crate) const EDWARDS_D: FieldElement32 = FieldElement32([ -10913610, 13857413, -15372611, 6949391, 114729, -8787816, -6275908, -3247719, -18696448, -12055116, ]); -pub const d2: FieldElement32 = FieldElement32([ +/// Edwards `2*d` value, equal to `2*(-121665/121666) mod p`. +pub(crate) const EDWARDS_D2: FieldElement32 = FieldElement32([ -21827239, -5839606, -30745221, 13898782, 229458, 15978800, -12551817, -6495438, 29715968, 9444199, ]); -pub const d4: FieldElement32 = FieldElement32([ - 23454405, -11679213, 5618422, -5756869, 458917, - -1596832, -25103633, -12990876, -7676928, -14666033 ]); +/// `= sqrt(a*d - 1)`, where `a = -1 (mod p)`, `d` are the Edwards curve parameters. +pub(crate) const SQRT_AD_MINUS_ONE: FieldElement32 = FieldElement32([ + 24849947, -153582, -23613485, 6347715, -21072328, -667138, -25271143, -15367704, -870347, 14525639 +]); -pub const a_minus_d: FieldElement32 = FieldElement32([ - 10913609, -13857413, 15372611, -6949391, -114729, - 8787816, 6275908, 3247719, 18696448, 12055116, ]); - -pub const invsqrt_a_minus_d: FieldElement32 = FieldElement32([ +/// `= 1/sqrt(a-d)`, where `a = -1 (mod p)`, `d` are the Edwards curve parameters. +pub(crate) const INVSQRT_A_MINUS_D: FieldElement32 = FieldElement32([ 6111485, 4156064, -27798727, 12243468, -25904040, 120897, 20826367, -7060776, 6093568, -1986012 ]); -#[cfg(not(feature="radix_51"))] -pub const inv_a_minus_d: FieldElement32 = FieldElement32([ - -121666, 0, 0, 0, 0, 0, 0, 0, 0, 0 -]); - -/// (p-1)/2, in little-endian bytes. -pub const HALF_P_MINUS_1_BYTES: [u8; 32] = - [0xf6, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, - 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, - 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, - 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x3f]; - /// Precomputed value of one of the square roots of -1 (mod p) -pub const SQRT_M1: FieldElement32 = FieldElement32([ +pub(crate) const SQRT_M1: FieldElement32 = FieldElement32([ -32595792, -7943725, 9377950, 3500415, 12389472, -272473, -25146209, -2005654, 326686, 11406482, ]); -/// Precomputed value of the other square root of -1 (mod p), -/// i.e., `MSQRT_M1 = -SQRT_M1`. -pub const MSQRT_M1: FieldElement32 = FieldElement32([ - 32595792, 7943725, -9377950, -3500415, -12389472, - 272473, 25146209, 2005654, -326686, -11406482, ]); - -/// Precomputed value of 1/2 (mod p). -pub const HALF: FieldElement32 = FieldElement32([ - 10, 0, 0, 0, 0, 0, 0, 0, 0, -16777216, ]); - /// In Montgomery form y² = x³+Ax²+x, Curve25519 has A=486662. -pub const A: FieldElement32 = FieldElement32([ +pub(crate) const MONTGOMERY_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) -// instead...? - hdevalence -pub const SQRT_MINUS_A: FieldElement32 = FieldElement32([ // sqrtMinusA - 12222970, 8312128, 11511410, -9067497, 15300785, - 241793, -25456130, -14121551, 12187136, -3972024, ]); +pub(crate) const APLUS2_OVER_FOUR: FieldElement32 = FieldElement32([121666, 0, 0, 0, 0, 0, 0, 0, 0, 0]); /// `SQRT_MINUS_APLUS2` is sqrt(-486664) -pub const SQRT_MINUS_APLUS2: FieldElement32 = FieldElement32([ +pub(crate) const SQRT_MINUS_APLUS2: FieldElement32 = FieldElement32([ -12222970, -8312128, -11511410, 9067497, -15300785, -241793, 25456130, 14121551, -12187136, 3972024]); @@ -96,8 +63,28 @@ pub const SQRT_MINUS_HALF: FieldElement32 = FieldElement32([ // sqrtMinusHalf -17256545, 3971863, 28865457, -1750208, 27359696, -16640980, 12573105, 1002827, -163343, 11073975, ]); -/// 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). +/// `L` is the order of base point, i.e. 2^252 + +/// 27742317777372353535851937790883648493 +pub(crate) const L: Scalar32 = Scalar32([ 0x1cf5d3ed, 0x009318d2, 0x1de73596, 0x1df3bd45, + 0x0000014d, 0x00000000, 0x00000000, 0x00000000, + 0x00100000 ]); + +/// `L` * `LFACTOR` = -1 (mod 2^29) +pub(crate) const LFACTOR: u32 = 0x12547e1b; + +/// `R` = R % L where R = 2^261 +pub(crate) const R: Scalar32 = Scalar32([ 0x114df9ed, 0x1a617303, 0x0f7c098c, 0x16793167, + 0x1ffd656e, 0x1fffffff, 0x1fffffff, 0x1fffffff, + 0x000fffff ]); + +/// `RR` = (R^2) % L where R = 2^261 +pub(crate) const RR: Scalar32 = Scalar32([ 0x0b5f9d12, 0x1e141b17, 0x158d7f3d, 0x143f3757, + 0x1972d781, 0x042feb7c, 0x1ceec73d, 0x1e184d1e, + 0x0005046d ]); + +/// The Ed25519 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{ X: FieldElement32([-14297830, -7645148, 16144683, -16471763, 27570974, -2696100, -26142465, 8378389, 20764389, 8758491]), Y: FieldElement32([-26843541, -6710886, 13421773, -13421773, 26843546, 6710886, -13421773, 13421773, -26843546, -6710886]), @@ -163,7 +150,8 @@ pub const EIGHT_TORSION: [ExtendedPoint; 8] = [ }, ]; -pub const bi: [AffineNielsPoint; 8] = [ +/// Odd multiples of the basepoint `[B, 3B, 5B, 7B, 9B, 11B, 13B, 15B]`. +pub(crate) const AFFINE_ODD_MULTIPLES_OF_BASEPOINT: [AffineNielsPoint; 8] = [ AffineNielsPoint{ y_plus_x: FieldElement32([25967493, -14356035, 29566456, 3660896, -12694345, 4014787, 27544626, -11754271, -6079156, 2047605]), y_minus_x: FieldElement32([-12545711, 934262, -2722910, 3049990, -727428, 9406986, 12720692, 5043384, 19500929, -15469378]), diff --git a/src/constants_64bit.rs b/src/constants_64bit.rs index ba7b2d9..e4e8fb4 100644 --- a/src/constants_64bit.rs +++ b/src/constants_64bit.rs @@ -12,65 +12,57 @@ //! and useful field elements like `sqrt(-1)`), as well as //! lookup tables of pre-computed points. -#![allow(dead_code)] -#![allow(non_snake_case)] -#![allow(non_upper_case_globals)] -#![allow(missing_docs)] #![allow(non_snake_case)] use field_64bit::FieldElement64; +use scalar_64bit::Scalar64; use edwards::ExtendedPoint; use edwards::AffineNielsPoint; use edwards::EdwardsBasepointTable; -pub const p: FieldElement64 = FieldElement64([2251799813685229, 2251799813685247, 2251799813685247, 2251799813685247, 2251799813685247]); +/// Edwards `d` value, equal to `-121665/121666 mod p`. +pub(crate) const EDWARDS_D: FieldElement64 = FieldElement64([929955233495203, 466365720129213, 1662059464998953, 2033849074728123, 1442794654840575]); -pub const d: FieldElement64 = FieldElement64([929955233495203, 466365720129213, 1662059464998953, 2033849074728123, 1442794654840575]); +/// Edwards `2*d` value, equal to `2*(-121665/121666) mod p`. +pub(crate) const EDWARDS_D2: FieldElement64 = FieldElement64([1859910466990425, 932731440258426, 1072319116312658, 1815898335770999, 633789495995903]); -pub const d2: FieldElement64 = FieldElement64([1859910466990425, 932731440258426, 1072319116312658, 1815898335770999, 633789495995903]); +/// `= sqrt(a*d - 1)`, where `a = -1 (mod p)`, `d` are the Edwards curve parameters. +pub(crate) const SQRT_AD_MINUS_ONE: FieldElement64 = FieldElement64([ + 2241493124984347, 425987919032274, 2207028919301688, 1220490630685848, 974799131293748 +]); -pub const d4: FieldElement64 = FieldElement64([1468021120295602, 1865462880516853, 2144638232625316, 1379996857856750, 1267578991991807]); - -pub const a_minus_d: FieldElement64 = FieldElement64([1321844580190025, 1785434093556034, 589740348686294, 217950738957124, 809005158844672]); - -pub const invsqrt_a_minus_d: FieldElement64 = FieldElement64([ +/// `= 1/sqrt(a-d)`, where `a = -1 (mod p)`, `d` are the Edwards curve parameters. +pub(crate) const INVSQRT_A_MINUS_D: FieldElement64 = FieldElement64([ 278908739862762, 821645201101625, 8113234426968, 1777959178193151, 2118520810568447 ]); -pub const inv_a_minus_d: FieldElement64 = FieldElement64([ - 2251799813563563, 2251799813685247, 2251799813685247, 2251799813685247, 2251799813685247 -]); - /// Precomputed value of one of the square roots of -1 (mod p) -pub const SQRT_M1: FieldElement64 = FieldElement64([1718705420411056, 234908883556509, 2233514472574048, 2117202627021982, 765476049583133]); - -/// Precomputed value of the other square root of -1 (mod p), -/// i.e., `MSQRT_M1 = -SQRT_M1`. -pub const MSQRT_M1: FieldElement64 = FieldElement64([533094393274173, 2016890930128738, 18285341111199, 134597186663265, 1486323764102114]); - -/// Precomputed value of 1/2 (mod p). -pub const HALF: FieldElement64 = FieldElement64([2251799813685239, 2251799813685247, 2251799813685247, 2251799813685247, 1125899906842623]); +pub(crate) const SQRT_M1: FieldElement64 = FieldElement64([1718705420411056, 234908883556509, 2233514472574048, 2117202627021982, 765476049583133]); /// In Montgomery form y² = x³+Ax²+x, Curve25519 has A=486662. -pub const A: FieldElement64 = FieldElement64([486662, 0, 0, 0, 0]); +pub(crate) const MONTGOMERY_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) -// instead...? - hdevalence -pub const SQRT_MINUS_A: FieldElement64 = FieldElement64([557817479725543, 1643290402203250, 16226468853936, 1304118542701054, 1985241807451647]); +pub(crate) const APLUS2_OVER_FOUR: FieldElement64 = FieldElement64([121666, 0, 0, 0, 0]); /// `SQRT_MINUS_APLUS2` is sqrt(-486664) -pub const SQRT_MINUS_APLUS2: FieldElement64 = FieldElement64([1693982333959686, 608509411481997, 2235573344831311, 947681270984193, 266558006233600]); +pub(crate) const SQRT_MINUS_APLUS2: FieldElement64 = FieldElement64([1693982333959686, 608509411481997, 2235573344831311, 947681270984193, 266558006233600]); -/// `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(crate) const L: Scalar64 = Scalar64([ 0x0002631a5cf5d3ed, 0x000dea2f79cd6581, 0x000000000014def9, 0x0000000000000000, 0x0000100000000000 ]); -/// 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). +/// `L` * `LFACTOR` = -1 (mod 2^51) +pub(crate) const LFACTOR: u64 = 0x51da312547e1b; + +/// `R` = R % L where R = 2^260 +pub(crate) const R: Scalar64 = Scalar64([ 0x000f48bd6721e6ed, 0x0003bab5ac67e45a, 0x000fffffeb35e51b, 0x000fffffffffffff, 0x00000fffffffffff ]); + +/// `RR` = (R^2) % L where R = 2^260 +pub(crate) const RR: Scalar64 = Scalar64([ 0x0009d265e952d13b, 0x000d63c715bea69f, 0x0005be65cb687604, 0x0003dceec73d217f, 0x000009411b7c309a ]); + +/// The Ed25519 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{ X: FieldElement64([1738742601995546, 1146398526822698, 2070867633025821, 562264141797630, 587772402128613]), Y: FieldElement64([1801439850948184, 1351079888211148, 450359962737049, 900719925474099, 1801439850948198]), @@ -143,7 +135,8 @@ pub const EIGHT_TORSION: [ExtendedPoint; 8] = [ } ]; -pub const bi: [AffineNielsPoint; 8] = [ +/// Odd multiples of the basepoint `[B, 3B, 5B, 7B, 9B, 11B, 13B, 15B]`. +pub(crate) const AFFINE_ODD_MULTIPLES_OF_BASEPOINT: [AffineNielsPoint; 8] = [ AffineNielsPoint { y_plus_x: FieldElement64([1288382639258501, 245678601348599, 269427782077623, 1462984067271730, 137412439391563]), y_minus_x: FieldElement64([62697248952638, 204681361388450, 631292143396476, 338455783676468, 1213667448819585]), diff --git a/src/decaf.rs b/src/decaf.rs deleted file mode 100644 index 9bc2c8a..0000000 --- a/src/decaf.rs +++ /dev/null @@ -1,871 +0,0 @@ -// -*- mode: rust; -*- -// -// This file is part of curve25519-dalek. -// Copyright (c) 2016-2017 Isis Lovecruft, Henry de Valence -// See LICENSE for licensing information. -// -// Authors: -// - Isis Agora Lovecruft -// - Henry de Valence - -//! An implementation of Mike Hamburg's Decaf cofactor-eliminating -//! point-compression scheme, providing a prime-order group on top of -//! a non-prime-order elliptic curve. -//! -//! Note: this code is currently feature-gated with the `yolocrypto` -//! feature flag, because our implementation is still unfinished. - -// We allow non snake_case names because coordinates in projective space are -// traditionally denoted by the capitalisation of their respective -// counterparts in affine space. Yeah, you heard me, rustc, I'm gonna have my -// affine and projective cakes and eat both of them too. -#![allow(non_snake_case)] - -use core::fmt::Debug; - -#[cfg(feature = "std")] -use rand::Rng; - -use digest::Digest; -use generic_array::typenum::U32; - -use constants; -use field::FieldElement; - -use core::ops::{Add, Sub, Neg}; -use core::ops::{AddAssign, SubAssign}; -use core::ops::{Mul, MulAssign}; - -use edwards; -use edwards::ExtendedPoint; -use edwards::CompletedPoint; -use edwards::EdwardsBasepointTable; -use edwards::Identity; -use scalar::Scalar; - -use subtle::ConditionallyAssignable; -use subtle::ConditionallyNegatable; - -// ------------------------------------------------------------------------ -// Compressed points -// ------------------------------------------------------------------------ - -/// A point serialized using Mike Hamburg's Decaf scheme. -/// -/// XXX think about how this API should work -#[derive(Copy, Clone, Eq, PartialEq)] -pub struct CompressedDecaf(pub [u8; 32]); - -/// The result of compressing a `DecafPoint`. -impl CompressedDecaf { - /// Convert this `CompressedDecaf` to an array of bytes. - pub fn to_bytes(&self) -> [u8; 32] { - self.0 - } - - /// View this `CompressedDecaf` as an array of bytes. - pub fn as_bytes<'a>(&'a self) -> &'a [u8; 32] { - &self.0 - } - - /// Attempt to decompress to an `DecafPoint`. - pub fn decompress(&self) -> Option { - // XXX should decoding be CT ? - // XXX need to check that xy is nonnegative and reject otherwise - let s = FieldElement::from_bytes(self.as_bytes()); - - // Check that s = |s| and reject otherwise. - let mut abs_s = s; - let neg = abs_s.is_negative_decaf(); - abs_s.conditional_negate(neg); - if abs_s != s { return None; } - - let ss = s.square(); - let X = &s + &s; // X = 2s - let Z = &FieldElement::one() - &ss; // Z = 1+as^2 - let ZZ = Z.square(); - let u = &ZZ- &(&constants::d4 * &ss); // u = Z^2 - 4ds^2 - let uss = &u * &ss; - let ussZZ = &uss * &ZZ; - - if Z.is_zero() == 1u8 { return None; } - - // Batch inversion: set b = 1/sqrt(us^2 Z^2) - let (ussZZ_is_nonzero_square, b) = ussZZ.invsqrt(); - if (ussZZ_is_nonzero_square | uss.is_zero()) == 0u8 { - return None; // us^2 is nonzero nonsquare - } - - let mut v = &b * &Z; // now v = 1/sqrt(us^2) - let Zinv = &b * &(&v * &uss); // now Zinv = b^2 Z us^2 = 1/Z - - // Now v = 1/sqrt(us^2) if us^2 is a nonzero square, 0 if us^2 is zero. - let uv = &v * &u; - if uv.is_negative_decaf() == 1u8 { - v.negate(); - } - let mut two_minus_Z = -&Z; two_minus_Z.0[0] += 2; - let mut w = &v * &(&s * &two_minus_Z); - w.conditional_assign(&FieldElement::one(), s.is_zero()); - let Y = &w * &Z; - let T = &w * &X; - - // "To decode the point, one must decode it to affine form - // instead of projective, and check that xy is non-negative." - - // Use the value of 1/Z previously computed in the batch inversion - let xy = &T * &Zinv; - if (Y.is_nonzero() & xy.is_nonnegative_decaf()) == 1u8 { - Some(DecafPoint(ExtendedPoint{ X: X, Y: Y, Z: Z, T: T })) - } else { - None - } - } -} - -impl Identity for CompressedDecaf { - fn identity() -> CompressedDecaf { - CompressedDecaf([0u8; 32]) - } -} - -// ------------------------------------------------------------------------ -// Serde support -// ------------------------------------------------------------------------ -// Serializes to and from `DecafPoint` directly, doing compression -// and decompression internally. This means that users can create -// structs containing `DecafPoint`s and use Serde's derived -// serializers to serialize those structures. - -#[cfg(feature = "serde")] -use serde::{self, Serialize, Deserialize, Serializer, Deserializer}; -#[cfg(feature = "serde")] -use serde::de::Visitor; - -#[cfg(feature = "serde")] -impl Serialize for DecafPoint { - fn serialize(&self, serializer: S) -> Result - where S: Serializer - { - serializer.serialize_bytes(self.compress().as_bytes()) - } -} - -#[cfg(feature = "serde")] -impl<'de> Deserialize<'de> for DecafPoint { - fn deserialize(deserializer: D) -> Result - where D: Deserializer<'de> - { - struct DecafPointVisitor; - - impl<'de> Visitor<'de> for DecafPointVisitor { - type Value = DecafPoint; - - fn expecting(&self, formatter: &mut ::core::fmt::Formatter) -> ::core::fmt::Result { - formatter.write_str("a valid point in Decaf format") - } - - fn visit_bytes(self, v: &[u8]) -> Result - where E: serde::de::Error - { - if v.len() == 32 { - let arr32 = array_ref!(v, 0, 32); // &[u8;32] from &[u8] - CompressedDecaf(*arr32) - .decompress() - .ok_or(serde::de::Error::custom("decompression failed")) - } else { - Err(serde::de::Error::invalid_length(v.len(), &self)) - } - } - } - - deserializer.deserialize_bytes(DecafPointVisitor) - } -} - -// ------------------------------------------------------------------------ -// Internal point representations -// ------------------------------------------------------------------------ - -/// A point in a prime-order group. -/// -// XXX think about how this API should work -#[derive(Copy, Clone)] -pub struct DecafPoint(pub ExtendedPoint); - -impl DecafPoint { - /// Compress in Decaf format. - pub fn compress(&self) -> CompressedDecaf { - // Q: Do we want to encode twisted or untwisted? - // - // Notes: - // Recall that the twisted Edwards curve E_{a,d} is of the form - // - // ax^2 + y^2 = 1 + dx^2y^2. - // - // Internally, we operate on the curve with a = -1, d = - // -121665/121666, a.k.a., the twist. But maybe we would like - // to use Decaf on the untwisted curve with a = 1, d = - // 121665/121666. (why? interop?) - // - // Fix i, a square root of -1 (mod p). - // - // The map x -> ix is an isomorphism from E_{a,d} to E_{-a,-d}. - // Its inverse is x -> -ix. - // let untwisted_X = &self.X * &constants::MSQRT_M1; - // etc. - // - // Step 0: pre-rotation, needed for Decaf with E[8] = Z/8. - // - // We want to select a point (x,y) in the coset P + E[4] with - // y nonzero and xy nonnegative. The naive approach is as - // follows. First, compute xy = T/Z and check that Y is - // nonzero and xy is nonnegative. If not, then "rotate" the - // original point by adding (x,y) + (i,0) = (iy,ix) = (x',y'): - // this rotated point has x'y' = -xy. Then perform the normal - // Decaf encoding, as described in Appendix A.1 of the Decaf - // paper, using the rotated point (x',y'). - // - // This is straightforward but requires an extra inversion. - // We would like to batch the inversion in xy = T/Z with the - // inverse square root in the computation of - // - // r = invsqrt((a-d)*(Z+Y)*(Z-Y)) - // = invsqrt(a-d)*invsqrt(Z^2-Y^2), - // - // but the X and Y we are trying to decode depend on whether - // we rotated the coset representative! - // - // However, it is possible to batch these inversions. Credit: - // the following explanation (and trick) is adapted from an - // email from Mike Hamburg, but of course any errors are ours. - // - // Let the initial point be ( X_0 : Y_0 : Z_0 : T_0). - // The rotated point is then (iY_0 : iX_0 : Z_0 : -T_0). - // - // We want to relate the computation of: - // - // invsqrt(Z^2 - Y^2) = invsqrt(Z_0^2 - Y_0^2) [non-rotated] - // invsqrt(Z^2 - Y^2) = invsqrt(Z_0^2 + X_0^2) [rotated] - // - // The curve equation in extended coordinates is - // - // 0 = (-X^2 + Y^2)*Z^2 - Z^4 - d*X^2*Y^2, - // - // so - // 0 = (-X^2 + Y^2)*Z^2 - Z^4 - d*T^2*Z^2 since XY=TZ - // = (-X^2 + Y^2 - Z^2 - d*T^2)*Z^2 - // = ( X^2 - Y^2 + Z^2 + d*T^2)*Z^2 mult by -1 - // -T^2*Z^2 = (X^2 - Y^2 + Z^2 + d*T^2)*Z^2 - T^2*Z^2 sub T^2*Z^2 - // -T^2*Z^2 = (X^2 - Y^2 + Z^2 - T^2)*Z^2 + d*T^2*Z^2 - // (-1-d)*T^2*Z^2 = (X^2 - Y^2 + Z^2 - T^2)*Z^2 - // - // for any point (X:Y:Z:T) in extended coordinates. Therefore, - // - // (Z^2 - Y^2)*(Z^2 + X^2) = Z^4 + Z^2*X^2 - Y^2*Z^2 - Y^2*X^2 - // = Z^4 + Z^2*X^2 - Y^2*Z^2 - T^2*Z^2 since XY=TZ - // = Z^2*(X^2 - Y^2 + Z^2 - T^2) - // = (-1-d)*T^2*Z^2. - // - // Taking square roots of both sides and rearranging, we get - // - // invsqrt(Z^2 - Y^2) = invsqrt(-1-d)*(Z^2+X^2)*(1/TZ)*invsqrt(Z^2+X^2) - // `-----------' `---------------------' - // curve constant batchable - // - // for any point (X:Y:Z:T) in extended coordinates. - // - // Therefore, we can do the computation with only one inverse - // square root like so: - // - // W <--- invsqrt((T_0 * Z_0)^2 * (Z_0^2+X_0^2)) - // = 1/(T_0 * Z_0 * sqrt(Z_0^2 + X_0^2)) - // - // xy <--- T_0^2 * W^2 * (T_0 * Z_0) * (Z_0^2 + X_0^2) - // = T_0 / Z_0 = xy - // - // if Y_0 nonzero and xy nonnegative: - // (X : Y : Z : T) <--- (X_0 : Y_0 : Z_0 : T_0) - // r <--- (1/(-1-d)) * (Z_0^2 + X_0^2) * W - // = invsqrt(a-d) * invsqrt(Z_0^2 - Y_0^2) since a = -1 - // = invsqrt(a-d) * invsqrt(Z^2 - Y^2) - // otherwise: - // (X : Y : Z : T) <--- (i*Y_0 : i*X_0 : Z_0 : -T_0) - // r <--- invsqrt(a-d) * (T_0 * Z_0) * W - // = invsqrt(a-d) * invsqrt(Z_0^2 + X_0^2) - // = invsqrt(a-d) * invsqrt(Z^2 - Y^2) - // - // The rest of the compression follows the steps in the - // appendix of the Decaf paper. - - let mut X = self.0.X; - let mut Y = self.0.Y; - let mut T = self.0.T; - let Z = &self.0.Z; - - let TZ = &T * Z; - let ZZ_plus_XX = &Z.square() + &X.square(); - let tmp = &TZ.square() * &ZZ_plus_XX; - let (tmp_is_nonzero_square, W) = tmp.invsqrt(); - // tmp should always be a square (why? related to being in the - // image of the isogeny?) - debug_assert_eq!(tmp_is_nonzero_square | tmp.is_zero(), 1u8); - - let xy = &T.square() * &(&W.square() * &(&TZ * &ZZ_plus_XX)); - let rotate = 1u8 & !(Y.is_nonzero() & xy.is_nonnegative_decaf()); - - let mut r = &W * &(&ZZ_plus_XX * &constants::inv_a_minus_d); - let r_rot = &W * &(&TZ * &constants::invsqrt_a_minus_d); - - let iX = &X * &constants::SQRT_M1; - let iY = &Y * &constants::SQRT_M1; - - r.conditional_assign(&r_rot, rotate); - X.conditional_assign(&iY, rotate); - Y.conditional_assign(&iX, rotate); - T.conditional_negate(rotate); - - // Step 2: Compute u = (a-d)r - let u = &constants::a_minus_d * &r; - - // Step 3: Negate r if -2uZ is negative. - let uZ = &u * Z; - let m2uZ = -&(&uZ + &uZ); - r.conditional_negate(m2uZ.is_negative_decaf()); - - // Step 4: Compute s = | u(r(aZX - dYT)+Y)/a| - // = |u(r(-ZX - dYT)+Y)| since a = -1 - let minus_ZX = -&(Z * &X); - let dYT = &constants::d * &(&Y * &T); - // Compute s = u(r(aZX - dYT)+Y) and cnegate for abs - let mut s = &u * &(&(&r * &(&minus_ZX - &dYT)) + &Y); - let neg = s.is_negative_decaf(); - s.conditional_negate(neg); - CompressedDecaf(s.to_bytes()) - } - - /// Return the coset self + E[4], for debugging. - fn coset4(&self) -> [ExtendedPoint; 4] { - [ self.0 - , &self.0 + &constants::EIGHT_TORSION[2] - , &self.0 + &constants::EIGHT_TORSION[4] - , &self.0 + &constants::EIGHT_TORSION[6] - ] - } - - /// Computes the Elligator map as described in the Decaf paper. - /// - /// # Note - /// - /// This method is not public because it's just used for hashing - /// to a point -- proper elligator support is deferred for now. - pub fn elligator_decaf_flavour(r_0: &FieldElement) -> DecafPoint { - // Follows Appendix C of the Decaf paper. - // Use n = 2 as the quadratic nonresidue so that n*x = x + x. - let minus_one = -&FieldElement::one(); - - // 1. Compute r <--- nr_0^2. - let r_0_squared = r_0.square(); - let r = &r_0_squared + &r_0_squared; - - // 2. Compute D <--- (dr + (a-d)) * (dr - (d + ar)) - let dr = &constants::d * &r; - // D = (dr + (a-d)) * (dr - (d + ar)) - // = (dr + (a-d)) * (dr - (d-r)) since a=-1 - // writing as - // = (dr + (a-d)) * dr - (dr + (a-d)) * (d - r) - // avoids two consecutive additions (could cause overflow) - let dr_plus_amd = &dr + &constants::a_minus_d; - let D = &(&dr_plus_amd * &dr) - &(&dr_plus_amd * &(&constants::d - &r)); - - // 3. Compute N <--- (r+1) * (a-2d) - let N = &(&r + &FieldElement::one()) * &(&minus_one - &constants::d2); - - // 4. Compute - // / +1, 1 / sqrt(ND) if ND is square - // c, e <--- | +1, 0 if N or D = 0 - // \ -1, nr_0 / sqrt(nND) otherwise - let ND = &N * &D; - let nND = &ND + &ND; - let mut c = FieldElement::one(); - let mut e = FieldElement::zero(); - let (ND_is_nonzero_square, ND_invsqrt) = ND.invsqrt(); - e.conditional_assign(&ND_invsqrt, ND_is_nonzero_square); - let (nND_is_nonzero_square, nND_invsqrt) = nND.invsqrt(); - let nr_0_nND_invsqrt = &nND_invsqrt * &(r_0 + r_0); - c.conditional_assign(&minus_one, nND_is_nonzero_square); - e.conditional_assign(&nr_0_nND_invsqrt, nND_is_nonzero_square); - - // 5. Compute s <--- c*|N*e| - let mut s = &N * &e; - let neg = s.is_negative_decaf(); - s.conditional_negate(neg); - s *= &c; - - // 6. Compute t <--- -c*N*(r-1)* ((a-2d)*e)^2 -1 - let a_minus_2d_e_sq = (&(&minus_one - &constants::d2) * &e).square(); - let c_N_r_minus_1 = &c * &(&N * &(&r + &minus_one)); - let t = &minus_one - &(&c_N_r_minus_1 * &a_minus_2d_e_sq); - - // 7. Apply the isogeny: - // (x,y) = ((2s)/(1+as^2), (1-as^2)/(t)) - let as_sq = &minus_one * &s.square(); - let P = CompletedPoint{ - X: &s + &s, - Z: &FieldElement::one() + &as_sq, - Y: &FieldElement::one() - &as_sq, - T: t, - }; - - // Convert to extended and return. - DecafPoint(P.to_extended()) - } - - /// Return a `DecafPoint` chosen uniformly at random using a user-provided RNG. - /// - /// # Inputs - /// - /// * `rng`: any RNG which implements the `rand::Rng` interface. - /// - /// # Returns - /// - /// A random element of the Decaf group. - /// - /// # Implementation - /// - /// Uses the Decaf-flavoured Elligator 2 map, so that the discrete log of the - /// output point with respect to any other point should be unknown. - #[cfg(feature = "std")] - pub fn random(rng: &mut T) -> Self { - let mut field_bytes = [0u8; 32]; - rng.fill_bytes(&mut field_bytes); - let r_0 = FieldElement::from_bytes(&field_bytes); - DecafPoint::elligator_decaf_flavour(&r_0) - } - - /// Hash a slice of bytes into a `DecafPoint`. - /// - /// Takes a type parameter `D`, which is any `Digest` producing 32 - /// bytes (256 bits) of output. - /// - /// Convenience wrapper around `from_hash`. - /// - /// # Implementation - /// - /// Uses the Decaf-flavoured Elligator 2 map, so that the discrete log of the - /// output point with respect to any other point should be unknown. - /// - /// # Example - /// - /// ``` - /// # extern crate curve25519_dalek; - /// # use curve25519_dalek::decaf::DecafPoint; - /// extern crate sha2; - /// use sha2::Sha256; - /// - /// # // Need fn main() here in comment so the doctest compiles - /// # // See https://doc.rust-lang.org/book/documentation.html#documentation-as-tests - /// # fn main() { - /// let msg = "To really appreciate architecture, you may even need to commit a murder"; - /// let P = DecafPoint::hash_from_bytes::(msg.as_bytes()); - /// # } - /// ``` - /// - pub fn hash_from_bytes(input: &[u8]) -> DecafPoint - where D: Digest + Default - { - let mut hash = D::default(); - hash.input(input); - DecafPoint::from_hash(hash) - } - - /// Construct a `DecafPoint` from an existing `Digest` instance. - /// - /// Use this instead of `hash_from_bytes` if it is more convenient - /// to stream data into the `Digest` than to pass a single byte - /// slice. - pub fn from_hash(hash: D) -> DecafPoint - where D: Digest + Default - { - // XXX this seems clumsy - let mut output = [0u8; 32]; - output.copy_from_slice(hash.result().as_slice()); - let r_0 = FieldElement::from_bytes(&output); - DecafPoint::elligator_decaf_flavour(&r_0) - } -} - -impl Identity for DecafPoint { - fn identity() -> DecafPoint { - DecafPoint(ExtendedPoint::identity()) - } -} - -// ------------------------------------------------------------------------ -// Equality -// ------------------------------------------------------------------------ - -/// XXX check whether there's a simple way to do equality checking -/// with cofactor 8, not just cofactor 4, and add a CT equality function? -impl PartialEq for DecafPoint { - fn eq(&self, other: &DecafPoint) -> bool { - let self_compressed = self.compress(); - let other_compressed = other.compress(); - self_compressed == other_compressed - } -} - -impl Eq for DecafPoint {} - -// ------------------------------------------------------------------------ -// Arithmetic -// ------------------------------------------------------------------------ - -impl<'a, 'b> Add<&'b DecafPoint> for &'a DecafPoint { - type Output = DecafPoint; - - fn add(self, other: &'b DecafPoint) -> DecafPoint { - DecafPoint(&self.0 + &other.0) - } -} - -impl<'b> AddAssign<&'b DecafPoint> for DecafPoint { - fn add_assign(&mut self, _rhs: &DecafPoint) { - *self = (self as &DecafPoint) + _rhs; - } -} - -impl<'a, 'b> Sub<&'b DecafPoint> for &'a DecafPoint { - type Output = DecafPoint; - - fn sub(self, other: &'b DecafPoint) -> DecafPoint { - DecafPoint(&self.0 - &other.0) - } -} - -impl<'b> SubAssign<&'b DecafPoint> for DecafPoint { - fn sub_assign(&mut self, _rhs: &DecafPoint) { - *self = (self as &DecafPoint) - _rhs; - } -} - -impl<'a> Neg for &'a DecafPoint { - type Output = DecafPoint; - - fn neg(self) -> DecafPoint { - DecafPoint(-&self.0) - } -} - -impl<'b> MulAssign<&'b Scalar> for DecafPoint { - fn mul_assign(&mut self, scalar: &'b Scalar) { - let result = (self as &DecafPoint) * scalar; - *self = result; - } -} - -impl<'a, 'b> Mul<&'b Scalar> for &'a DecafPoint { - type Output = DecafPoint; - /// Scalar multiplication: compute `scalar * self`. - fn mul(self, scalar: &'b Scalar) -> DecafPoint { - DecafPoint(&self.0 * scalar) - } -} - -impl<'a, 'b> Mul<&'b DecafPoint> for &'a Scalar { - type Output = DecafPoint; - - /// Scalar multiplication: compute `self * scalar`. - fn mul(self, point: &'b DecafPoint) -> DecafPoint { - DecafPoint(self * &point.0) - } -} - -/// Given a vector of (possibly secret) scalars and a vector of -/// (possibly secret) points, compute `c_1 P_1 + ... + c_n P_n`. -/// -/// This function has the same behaviour as -/// `vartime::multiscalar_mult` but is constant-time. -/// -/// # Input -/// -/// A vector of `Scalar`s and a vector of `DecafPoints`. It is an -/// error to call this function with two vectors of different lengths. -#[cfg(any(feature = "alloc", feature = "std"))] -pub fn multiscalar_mult<'a, 'b, I, J>(scalars: I, points: J) -> DecafPoint - where I: IntoIterator, - J: IntoIterator, -{ - let extended_points = points.into_iter().map(|P| &P.0); - DecafPoint(edwards::multiscalar_mult(scalars, extended_points)) -} - -/// Precomputation -#[derive(Clone)] -pub struct DecafBasepointTable(pub EdwardsBasepointTable); - -impl<'a, 'b> Mul<&'b Scalar> for &'a DecafBasepointTable { - type Output = DecafPoint; - - fn mul(self, scalar: &'b Scalar) -> DecafPoint { - DecafPoint(&self.0 * scalar) - } -} - -impl<'a, 'b> Mul<&'a DecafBasepointTable> for &'b Scalar { - type Output = DecafPoint; - - fn mul(self, basepoint_table: &'a DecafBasepointTable) -> DecafPoint { - DecafPoint(self * &basepoint_table.0) - } -} - -impl DecafBasepointTable { - /// Create a precomputed table of multiples of the given `basepoint`. - pub fn create(basepoint: &DecafPoint) -> DecafBasepointTable { - DecafBasepointTable(EdwardsBasepointTable::create(&basepoint.0)) - } - - /// Get the basepoint for this table as a `DecafPoint`. - pub fn basepoint(&self) -> DecafPoint { - DecafPoint(self.0.basepoint()) - } -} - -// ------------------------------------------------------------------------ -// Constant-time conditional assignment -// ------------------------------------------------------------------------ - -impl ConditionallyAssignable for DecafPoint { - /// Conditionally assign `other` to `self`, if `choice == 1u8`. - /// - /// # Example - /// - /// ``` - /// # extern crate subtle; - /// # extern crate curve25519_dalek; - /// # - /// # use subtle::ConditionallyAssignable; - /// # - /// # use curve25519_dalek::edwards::Identity; - /// # use curve25519_dalek::decaf::DecafPoint; - /// # use curve25519_dalek::constants; - /// # fn main() { - /// let A = DecafPoint::identity(); - /// let B = constants::DECAF_ED25519_BASEPOINT_POINT; - /// - /// let mut P = A; - /// - /// P.conditional_assign(&B, 0u8); - /// assert!(P == A); - /// P.conditional_assign(&B, 1u8); - /// assert!(P == B); - /// # } - /// ``` - fn conditional_assign(&mut self, other: &DecafPoint, choice: u8) { - self.0.X.conditional_assign(&other.0.X, choice); - self.0.Y.conditional_assign(&other.0.Y, choice); - self.0.Z.conditional_assign(&other.0.Z, choice); - self.0.T.conditional_assign(&other.0.T, choice); - } -} - -// ------------------------------------------------------------------------ -// Debug traits -// ------------------------------------------------------------------------ - -impl Debug for CompressedDecaf { - fn fmt(&self, f: &mut ::core::fmt::Formatter) -> ::core::fmt::Result { - write!(f, "CompressedDecaf: {:?}", self.as_bytes()) - } -} - -impl Debug for DecafPoint { - fn fmt(&self, f: &mut ::core::fmt::Formatter) -> ::core::fmt::Result { - let coset = self.coset4(); - write!(f, "DecafPoint: coset \n{:?}\n{:?}\n{:?}\n{:?}", - coset[0], coset[1], coset[2], coset[3]) - } -} - -// ------------------------------------------------------------------------ -// Variable-time functions -// ------------------------------------------------------------------------ - -pub mod vartime { - //! Variable-time operations on decaf points, useful for non-secret data. - use super::*; - - /// Given a vector of public scalars and a vector of (possibly secret) - /// points, compute - /// - /// c_1 P_1 + ... + c_n P_n. - /// - /// # Input - /// - /// A vector of `Scalar`s and a vector of `ExtendedPoints`. It is an - /// error to call this function with two vectors of different lengths. - pub fn multiscalar_mult<'a, 'b, I, J>(scalars: I, points: J) -> DecafPoint - where I: IntoIterator, - J: IntoIterator - { - let extended_points = points.into_iter().map(|P| &P.0); - DecafPoint(edwards::vartime::multiscalar_mult(scalars, extended_points)) - } -} - -// ------------------------------------------------------------------------ -// Tests -// ------------------------------------------------------------------------ - -#[cfg(test)] -mod test { - use rand::OsRng; - - use scalar::Scalar; - use constants; - use edwards::CompressedEdwardsY; - use edwards::Identity; - use edwards::ValidityCheck; - use super::*; - - #[cfg(feature = "serde")] - use serde_cbor; - - #[test] - #[cfg(feature = "serde")] - fn serde_cbor_basepoint_roundtrip() { - let output = serde_cbor::to_vec(&constants::DECAF_ED25519_BASEPOINT_POINT).unwrap(); - let parsed: DecafPoint = serde_cbor::from_slice(&output).unwrap(); - assert_eq!(parsed, constants::DECAF_ED25519_BASEPOINT_POINT); - } - - - #[test] - fn decaf_decompress_negative_s_fails() { - // constants::d is neg, so decompression should fail as |d| != d. - let bad_compressed = CompressedDecaf(constants::d.to_bytes()); - assert!(bad_compressed.decompress().is_none()); - } - - #[test] - fn decaf_decompress_id() { - let compressed_id = CompressedDecaf::identity(); - let id = compressed_id.decompress().unwrap(); - assert_eq!(id.0.compress(), CompressedEdwardsY::identity()); - } - - #[test] - fn decaf_compress_id() { - let id = DecafPoint::identity(); - assert_eq!(id.compress(), CompressedDecaf::identity()); - } - - #[test] - fn decaf_basepoint_roundtrip() { - let bp_compressed_decaf = constants::DECAF_ED25519_BASEPOINT_POINT.compress(); - let bp_recaf = bp_compressed_decaf.decompress().unwrap().0; - // 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(), CompressedEdwardsY::identity()); - } - - #[test] - fn encodings_of_small_multiples_of_basepoint() { - // Table of encodings of (1+i)*basepoint - // Generated using the previous naive implementation. - let compressed = [ - CompressedDecaf([141, 190, 226, 107, 177, 201, 35, 118, 14, 55, 160, 165, 242, 207, 121, 161, 177, 80, 8, 132, 205, 254, 101, 169, 233, 65, 124, 96, 255, 182, 249, 40]), - CompressedDecaf([131, 57, 148, 16, 8, 196, 141, 82, 144, 220, 105, 112, 66, 33, 48, 16, 182, 198, 173, 35, 248, 181, 92, 231, 222, 35, 85, 56, 5, 252, 91, 40]), - CompressedDecaf([199, 132, 32, 144, 156, 143, 81, 170, 240, 56, 232, 6, 178, 37, 118, 190, 110, 201, 26, 173, 156, 97, 59, 162, 240, 247, 226, 107, 197, 111, 107, 26]), - CompressedDecaf([210, 120, 34, 214, 175, 27, 61, 6, 229, 181, 216, 36, 11, 245, 146, 232, 130, 215, 77, 29, 210, 30, 54, 155, 191, 81, 59, 124, 174, 3, 135, 36]), - CompressedDecaf([155, 52, 159, 52, 189, 27, 181, 0, 245, 131, 0, 197, 79, 208, 252, 122, 104, 161, 245, 143, 67, 94, 13, 129, 153, 173, 129, 179, 118, 231, 90, 52]), - CompressedDecaf([42, 117, 252, 118, 8, 1, 72, 25, 111, 246, 247, 103, 236, 86, 235, 29, 100, 156, 186, 209, 159, 21, 61, 26, 249, 25, 137, 228, 84, 23, 10, 27]), - CompressedDecaf([21, 126, 181, 117, 58, 90, 216, 28, 184, 57, 9, 23, 158, 68, 159, 171, 109, 150, 232, 140, 144, 73, 139, 122, 124, 105, 125, 160, 94, 185, 150, 52]), - CompressedDecaf([232, 167, 112, 233, 126, 33, 105, 63, 151, 6, 88, 225, 181, 17, 223, 12, 116, 138, 203, 47, 243, 225, 50, 171, 21, 220, 186, 179, 132, 20, 48, 6]), - CompressedDecaf([99, 44, 97, 48, 242, 174, 78, 198, 112, 154, 146, 36, 239, 34, 94, 4, 0, 244, 175, 34, 46, 0, 83, 187, 5, 163, 225, 63, 51, 237, 234, 22]), - CompressedDecaf([2, 33, 89, 176, 178, 123, 159, 75, 235, 172, 251, 11, 137, 177, 90, 122, 149, 186, 52, 243, 153, 190, 185, 202, 59, 137, 204, 160, 150, 152, 148, 55]), - CompressedDecaf([245, 79, 78, 226, 114, 69, 247, 112, 18, 54, 90, 225, 176, 77, 231, 235, 196, 123, 49, 221, 34, 205, 151, 228, 244, 112, 82, 58, 30, 31, 58, 12]), - CompressedDecaf([135, 53, 175, 167, 13, 94, 62, 31, 29, 248, 13, 132, 29, 69, 7, 188, 145, 49, 62, 55, 181, 109, 214, 11, 248, 162, 70, 15, 236, 126, 100, 60]), - CompressedDecaf([98, 150, 69, 229, 144, 122, 237, 107, 127, 177, 33, 64, 59, 173, 210, 102, 74, 34, 23, 16, 252, 117, 14, 97, 231, 178, 63, 193, 157, 28, 178, 17]), - CompressedDecaf([222, 104, 6, 1, 72, 12, 72, 178, 204, 238, 128, 70, 41, 150, 235, 96, 153, 150, 18, 4, 141, 206, 0, 38, 122, 112, 249, 51, 94, 251, 20, 57]), - CompressedDecaf([7, 221, 140, 57, 13, 146, 248, 27, 56, 4, 128, 23, 145, 120, 126, 4, 158, 173, 52, 213, 164, 250, 26, 55, 89, 96, 187, 111, 211, 18, 63, 19]), - CompressedDecaf([91, 213, 193, 10, 102, 92, 199, 124, 61, 176, 1, 47, 111, 59, 183, 91, 79, 56, 208, 109, 172, 209, 17, 167, 229, 216, 3, 236, 200, 208, 15, 20]), - ]; - let mut bp = constants::DECAF_ED25519_BASEPOINT_POINT; - for i in 0..16 { - assert_eq!(bp.compress(), compressed[i]); - bp = &bp + &constants::DECAF_ED25519_BASEPOINT_POINT; - } - } - - #[test] - fn decaf_four_torsion_basepoint() { - let bp = constants::DECAF_ED25519_BASEPOINT_POINT; - let bp_coset = bp.coset4(); - for i in 0..4 { - assert_eq!(bp, DecafPoint(bp_coset[i])); - } - } - - #[test] - fn decaf_four_torsion_random() { - let mut rng = OsRng::new().unwrap(); - let B = &constants::DECAF_ED25519_BASEPOINT_TABLE; - let P = B * &Scalar::random(&mut rng); - let P_coset = P.coset4(); - for i in 0..4 { - assert_eq!(P, DecafPoint(P_coset[i])); - } - } - - #[test] - fn decaf_random_roundtrip() { - let mut rng = OsRng::new().unwrap(); - let B = &constants::DECAF_ED25519_BASEPOINT_TABLE; - for _ in 0..100 { - let P = B * &Scalar::random(&mut rng); - let compressed_P = P.compress(); - let Q = compressed_P.decompress().unwrap(); - assert_eq!(P, Q); - } - } - - #[test] - fn decaf_random_is_valid() { - let mut rng = OsRng::new().unwrap(); - for _ in 0..100 { - let P = DecafPoint::random(&mut rng); - // Check that P is on the curve - assert!(P.0.is_valid()); - // Check that P is in the image of the decaf map - P.compress(); - } - } -} - -#[cfg(all(test, feature = "bench"))] -mod bench { - use rand::OsRng; - use test::Bencher; - - use super::*; - - #[bench] - fn decompression(b: &mut Bencher) { - let mut rng = OsRng::new().unwrap(); - let B = &constants::DECAF_ED25519_BASEPOINT_TABLE; - let P = B * &Scalar::random(&mut rng); - let P_compressed = P.compress(); - b.iter(|| P_compressed.decompress().unwrap()); - } - - #[bench] - fn compression(b: &mut Bencher) { - let mut rng = OsRng::new().unwrap(); - let B = &constants::DECAF_ED25519_BASEPOINT_TABLE; - let P = B * &Scalar::random(&mut rng); - b.iter(|| P.compress()); - } -} diff --git a/src/edwards.rs b/src/edwards.rs index 0b8b333..8fad123 100644 --- a/src/edwards.rs +++ b/src/edwards.rs @@ -135,15 +135,15 @@ impl CompressedEdwardsY { let Y = FieldElement::from_bytes(self.as_bytes()); let Z = FieldElement::one(); let YY = Y.square(); - let u = &YY - &Z; // u = y²-1 - let v = &(&YY * &constants::d) + &Z; // v = dy²+1 + let u = &YY - &Z; // u = y²-1 + let v = &(&YY * &constants::EDWARDS_D) + &Z; // v = dy²+1 let (is_nonzero_square, mut X) = FieldElement::sqrt_ratio(&u, &v); if is_nonzero_square != 1u8 { return None; } // Flip the sign of X if it's not correct let compressed_sign_bit = self.as_bytes()[31] >> 7; - let current_sign_bit = X.is_negative_ed25519(); + let current_sign_bit = X.is_negative(); X.conditional_negate(current_sign_bit ^ compressed_sign_bit); Some(ExtendedPoint{ X: X, Y: Y, Z: Z, T: &X * &Y }) @@ -342,7 +342,7 @@ impl ValidityCheck for ProjectivePoint { let ZZ = self.Z.square(); let ZZZZ = ZZ.square(); let lhs = &(&YY - &XX) * &ZZ; - let rhs = &ZZZZ + &(&constants::d * &(&XX * &YY)); + let rhs = &ZZZZ + &(&constants::EDWARDS_D * &(&XX * &YY)); lhs == rhs } @@ -442,7 +442,7 @@ impl ProjectivePoint { let mut s: [u8; 32]; s = y.to_bytes(); - s[31] ^= (x.is_negative_ed25519() << 7) as u8; + s[31] ^= (x.is_negative() << 7) as u8; CompressedEdwardsY(s) } @@ -518,7 +518,7 @@ impl ExtendedPoint { Y_plus_X: &self.Y + &self.X, Y_minus_X: &self.Y - &self.X, Z: self.Z, - T2d: &self.T * &constants::d2, + T2d: &self.T * &constants::EDWARDS_D2, } } @@ -541,7 +541,7 @@ impl ExtendedPoint { let recip = self.Z.invert(); let x = &self.X * &recip; let y = &self.Y * &recip; - let xy2d = &(&x * &y) * &constants::d2; + let xy2d = &(&x * &y) * &constants::EDWARDS_D2; AffineNielsPoint{ y_plus_x: &y + &x, y_minus_x: &y - &x, @@ -1233,6 +1233,7 @@ pub mod vartime { } let odd_multiples_of_A = OddMultiples::create(A); + let odd_multiples_of_B = &constants::AFFINE_ODD_MULTIPLES_OF_BASEPOINT; let mut r = ProjectivePoint::identity(); loop { @@ -1245,9 +1246,9 @@ pub mod vartime { } if b_naf[i] > 0 { - t = &t.to_extended() + &constants::bi[( b_naf[i]/2) as usize]; + t = &t.to_extended() + &odd_multiples_of_B[( b_naf[i]/2) as usize]; } else if b_naf[i] < 0 { - t = &t.to_extended() - &constants::bi[(-b_naf[i]/2) as usize]; + t = &t.to_extended() - &odd_multiples_of_B[(-b_naf[i]/2) as usize]; } r = t.to_projective(); @@ -1269,8 +1270,6 @@ pub mod vartime { #[cfg(test)] mod test { - #[cfg(feature = "yolocrypto")] - use decaf::DecafPoint; use field::FieldElement; use scalar::Scalar; use subtle::ConditionallyAssignable; @@ -1433,7 +1432,7 @@ mod test { #[test] fn basepoint_mult_by_basepoint_order() { let B = &constants::ED25519_BASEPOINT_TABLE; - let should_be_id = B * &constants::l; + let should_be_id = B * &constants::BASEPOINT_ORDER; assert!(should_be_id.is_identity()); } @@ -1550,18 +1549,6 @@ mod test { assert!(P1.compress().to_bytes() == P2.compress().to_bytes()); } - #[test] - #[cfg(feature = "yolocrypto")] - fn scalarmult_decafpoint_works_both_ways() { - let P: DecafPoint = DecafPoint(constants::ED25519_BASEPOINT_POINT); - let s: Scalar = A_SCALAR; - - let P1 = &P * &s; - let P2 = &s * &P; - - assert!(P1.compress().as_bytes() == P2.compress().as_bytes()); - } - mod vartime { use super::super::*; use super::{A_SCALAR, B_SCALAR, A_TIMES_BASEPOINT, DOUBLE_SCALAR_MULT_RESULT}; diff --git a/src/field.rs b/src/field.rs index 730c3fa..840355d 100644 --- a/src/field.rs +++ b/src/field.rs @@ -23,6 +23,7 @@ use core::cmp::{Eq, PartialEq}; use subtle::slices_equal; use subtle::byte_is_nonzero; use subtle::ConditionallyAssignable; +use subtle::ConditionallyNegatable; use subtle::Equal; use constants; @@ -81,44 +82,11 @@ impl FieldElement { /// # Return /// /// If negative, return `1u8`. Otherwise, return `0u8`. - pub fn is_negative_ed25519(&self) -> u8 { //FeIsNegative + pub fn is_negative(&self) -> u8 { let bytes = self.to_bytes(); (bytes[0] & 1) as u8 } - /// Determine if this `FieldElement` is negative, in the - /// sense used by Decaf: `x` is nonnegative if the least - /// absolute residue for `x` lies in `[0, (p-1)/2]`, and - /// is negative otherwise. - /// - /// # Return - /// - /// Returns `1u8` if negative, `0u8` if nonnegative. - /// - /// # Implementation - /// - /// Uses a trick borrowed from Mike Hamburg's code. Let `x \in - /// F_p` and let `y \in Z` be the least absolute residue for `x`. - /// Suppose `y ≤ (p-1)/2`. Then `2y < p` so `2y = 2y mod p` and - /// `2y mod p` is even. On the other hand, if `y > (p-1)/2` then - /// `2y ≥ p`; since `y < p`, `2y \in [p, 2p)`, so `2y mod p = - /// 2y-p`, which is odd. - /// - /// Thus we can test whether `y ≤ (p-1)/2` by checking whether `2y - /// mod p` is even. - pub fn is_negative_decaf(&self) -> u8 { - let y = self + self; - (y.to_bytes()[0] & 1) as u8 - } - - /// Determine if this `FieldElement` is nonnegative, in the - /// sense used by Decaf: `x` is nonnegative if the least - /// absolute residue for `x` lies in `[0, (p-1)/2]`, and - /// is negative otherwise. - pub fn is_nonnegative_decaf(&self) -> u8 { - 1u8 & (!self.is_negative_decaf()) - } - /// Determine if this `FieldElement` is zero. /// /// # Return @@ -232,6 +200,8 @@ impl FieldElement { /// Given `FieldElements` `u` and `v`, attempt to compute /// `sqrt(u/v)` in constant time. /// + /// This function always returns the nonnegative square root, if it exists. + /// /// It would be much better to use an `Option` type here, but /// doing so forces the caller to branch, which we don't want to /// do. This seems like the least bad solution. @@ -278,6 +248,10 @@ impl FieldElement { let r_prime = &constants::SQRT_M1 * &r; r.conditional_assign(&r_prime, flipped_sign_sqrt); + // Choose the nonnegative square root. + let r_is_negative = r.is_negative(); + r.conditional_negate(r_is_negative); + let was_nonzero_square = correct_sign_sqrt | flipped_sign_sqrt; (was_nonzero_square, r) @@ -465,6 +439,20 @@ mod test { x.conditional_negate(1u8); assert_eq!(x, one); } + + #[test] + fn encoding_is_canonical() { + // Encode 1 wrongly as 1 + (2^255 - 19) = 2^255 - 18 + let one_encoded_wrongly_bytes: [u8;32] = [0xee, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x7f]; + // Decode to a field element + let one = FieldElement::from_bytes(&one_encoded_wrongly_bytes); + // .. then check that the encoding is correct + let one_bytes = one.to_bytes(); + assert_eq!(one_bytes[0], 1); + for i in 1..32 { + assert_eq!(one_bytes[i], 0); + } + } } #[cfg(all(test, feature = "bench"))] diff --git a/src/lib.rs b/src/lib.rs index 31c18ff..524cd4f 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -38,9 +38,8 @@ #[cfg(all(test, feature = "bench"))] 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,12 +70,15 @@ 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; -// Feature gate decaf while our implementation is unfinished and probably incorrect. -#[cfg(feature = "yolocrypto")] -pub mod decaf; +pub mod ristretto; // Other miscelaneous utilities. diff --git a/src/montgomery.rs b/src/montgomery.rs index 9ba71ed..45bd8d5 100644 --- a/src/montgomery.rs +++ b/src/montgomery.rs @@ -146,8 +146,9 @@ impl CompressedMontgomeryU { /// actually a square and `0` if otherwise, along with a `FieldElement`: the /// Montgomery `v` corresponding to this `u`. pub fn to_montgomery_v(u: &FieldElement) -> (u8, FieldElement) { + let A = &constants::MONTGOMERY_A; let one: FieldElement = FieldElement::one(); - let v_squared: FieldElement = u * &(&u.square() + &(&(&constants::A * u) + &one)); + let v_squared: FieldElement = u * &(&u.square() + &(&(A * u) + &one)); let (okay, v_inv) = v_squared.invsqrt(); let v = &v_inv * &v_squared; @@ -183,7 +184,7 @@ impl CompressedMontgomeryU { pub fn to_edwards_x(u: &FieldElement, v: &FieldElement, sign: &u8) -> FieldElement { let mut x: FieldElement = &(u * &v.invert()) * &constants::SQRT_MINUS_APLUS2; let neg_x: FieldElement = -(&x); - let current_sign: u8 = x.is_negative_ed25519(); + let current_sign: u8 = x.is_negative(); // Negate x to match the sign: x.conditional_assign(&neg_x, current_sign ^ sign); diff --git a/src/ristretto.rs b/src/ristretto.rs new file mode 100644 index 0000000..f537b25 --- /dev/null +++ b/src/ristretto.rs @@ -0,0 +1,1245 @@ +// -*- mode: rust; -*- +// +// This file is part of curve25519-dalek. +// Copyright (c) 2016-2017 Isis Lovecruft, Henry de Valence +// See LICENSE for licensing information. +// +// Authors: +// - Isis Agora Lovecruft +// - Henry de Valence + +//! An implementation of Ristretto, which provides a prime-order group. +//! +//! Ristretto is a modification of Mike Hamburg's [Decaf +//! cofactor-eliminating point-compression +//! scheme](https://eprint.iacr.org/2015/673.pdf) to work on top of the +//! Curve25519 group. +//! +//! Below are some notes on Ristretto, which are *NOT* a full writeup and which may have errors. +//! +//! # Notes on Ristretto +//! +//! ## Decaf +//! +//! The introduction of the Decaf paper, [_Decaf: Eliminating cofactors +//! through point compression_](https://eprint.iacr.org/2015/673.pdf) +//! notes that while most cryptographic systems require a group of prime +//! order, most concrete implementations using elliptic curve groups +//! fall short -- they either provide a group of prime order, but with +//! incomplete or variable-time addition formulae (for instance, most +//! Weierstrass models), or else they provide a fast and safe +//! implementation of a group whose order is not quite a prime \\(q\\), +//! but \\(hq\\) for a small cofactor \\(h\\) (for instance, Edwards +//! curves, which have cofactor at least \\(4\\)). +//! +//! This abstraction mismatch requires ad-hoc protocol modifications to +//! ensure security; these modifications require careful analysis and +//! are a recurring source of vulnerabilities. +//! +//! The Decaf suggestion is to use a quotient group, such as \\(\mathcal +//! E / \mathcal E[4]\\) or \\(2 \mathcal E / \mathcal E[2] \\), to +//! implement a prime-order group. +//! +//! This requires only changing +//! +//! 1. the function for equality checking (so that two representatives +//! of the same coset are considered equal); +//! 2. the function for encoding (so that two representatives of the +//! same coset are encoded as identical bitstrings); +//! 3. the function for decoding (so that only the canonical encoding of +//! a coset is accepted). +//! +//! Internally, each coset is represented by a curve point; two points +//! may represent the same coset in the same way that two points with +//! different \\(X,Y,Z\\) coordinates may represent the same point. The +//! group operations are carried out using the fast, safe Edwards +//! formulas. +//! +//! The Decaf paper suggests implementing the compression and +//! decompression routines using an isogeny from a Jacobi quartic; for +//! curves of cofactor \\(4\\), this eliminates the cofactor, and +//! explains the name: Decaf is named "after the procedure which divides +//! the effect of coffee by \\(4\\)". However, Curve25519 has a +//! cofactor of \\(8\\). To eliminate its cofactor, we tweak Decaf to +//! restrict further. This gives the +//! [Ristretto](https://en.wikipedia.org/wiki/Ristretto) encoding. +//! +//! ## The Jacobi Quartic +//! +//! The Jacobi quartic is parameterized by \\(e, A\\), and is of the +//! form $$ \mathcal J\_{e,A} : t\^2 = es\^4 + 2As\^2 + 1, $$ with +//! identity point \\((0,1)\\). For more details on the Jacobi quartic, +//! see the [Decaf paper](https://eprint.iacr.org/2015/673.pdf) or +//! [_Jacobi Quartic Curves +//! Revisited_](https://eprint.iacr.org/2009/312.pdf) by Hisil, Wong, +//! Carter, and Dawson). +//! +//! When \\(e = a\^2\\), \\(\mathcal J\_{e,A}\\) has full +//! \\(2\\)-torsion (i.e., \\(\mathcal J[2] \cong \mathbb Z /2 \times +//! \mathbb Z/2\\)), and +//! we can write the \\(\mathcal J[2]\\)-coset of a point \\(P = +//! (s,t)\\) as +//! $$ +//! P + \mathcal J[2] = \left\\{ +//! (s,t), +//! (-s,-t), +//! (1/as, -t/as\^2), +//! (-1/as, t/as\^2) \right\\}. +//! $$ +//! Notice that replacing \\(a\\) by \\(-a\\) just swaps the last two +//! points, so this set does not depend on the choice of \\(a\\). In +//! what follows we require \\(a = \pm 1\\). +//! +//! ## Encoding \\(\mathcal J / \mathcal J[2]\\) +//! +//! To encode points on \\(\mathcal J\\) modulo \\(\mathcal J[2]\\), +//! we need to choose a canonical representative of the above coset. +//! To do this, it's sufficient to make two independent sign choices: +//! the Decaf paper suggests choosing \\((s,t)\\) with \\(s\\) +//! non-negative and finite, and \\(t/s\\) non-negative or infinite. +//! +//! The encoding is then the (canonical byte encoding of the) +//! \\(s\\)-value of the canonical representative. +//! +//! ## The Edwards Curve +//! +//! Our primary internal model for Curve25519 points are the [_Extended +//! Twisted Edwards Coordinates_](https://eprint.iacr.org/2008/522.pdf) +//! of Hisil, Wong, Carter, and Dawson. +//! These correspond to the affine model +//! +//! $$\mathcal E\_{a,d} : ax\^2 + y\^2 = 1 + dx\^2y\^2.$$ +//! +//! In projective coordinates, we represent a point as \\((X:Y:Z:T)\\) +//! with $$XY = ZT, \quad aX\^2 + Y\^2 = Z\^2 + dT\^2.$$ (For more +//! details on this model, see the documentation for the `edwards` +//! module). The case \\(a = 1\\) is the _untwisted_ case; we only +//! consider \\(a = \pm 1\\), and in particular we focus on the twisted +//! Edwards form of Curve25519, which has \\(a = -1, d = +//! -121665/121666\\). When not otherwise specified, we write +//! \\(\mathcal E\\) for \\(\mathcal E\_{-1, -121665/121666}\\). +//! +//! When both \\(d\\) and \\(ad\\) are nonsquare (which forces \\(a\\) +//! to be square), the curve is *complete*. In this case the +//! four-torsion subgroup is cyclic, and we +//! can write it explicitly as +//! $$ +//! \mathcal E\_{a,d}[4] = \\{ (0,1),\; (1/\sqrt a, 0),\; (0, -1),\; (-1/\sqrt{a}, 0)\\}. +//! $$ +//! These are the only points with \\(xy = 0\\); the points with \\( y +//! \neq 0 \\) are \\(2\\)-torsion. The \\(\mathcal +//! E\_{a,d}[4]\\)-coset of \\(P = (x,y)\\) is then +//! $$ +//! P + \mathcal E\_{a,d}[4] = \\{ (x,y),\; (y/\sqrt a, -x\sqrt a),\; (-x, -y),\; (-y/\sqrt a, x\sqrt a)\\}. +//! $$ +//! Notice that if \\(xy \neq 0 \\), then exactly two of +//! these points have \\( xy \\) non-negative, and they differ by the +//! \\(2\\)-torsion point \\( (0,-1) \\). This means that we can select +//! a representative modulo \\(\mathcal +//! E\_{a,d}[2] \\) by requiring \\(xy\\) nonnegative and \\(y \neq +//! 0\\), and we can ensure this condition by conditionally adding a +//! \\(4\\)-torsion point if \\(xy\\) is negative or \\(y = 0\\). +//! +//! This procedure gives a canonical lift from \\(\mathcal E / \mathcal +//! E[4]\\) to \\(\mathcal E / \mathcal E[2]\\). Since it involves a +//! conditional rotation, we refer to it as *torquing* the point. +//! +//! The structure of the Curve25519 group is \\( \mathcal E(\mathbb +//! F\_p) \cong \mathbb Z / 8 \times \mathbb Z / \ell\\), where \\( \ell +//! = 2\^{252} + \cdots \\) is a large prime. Because \\(\mathcal E[8] +//! \cong \mathbb Z / 8\\), we have \\(\[2\](\mathcal E[8]) = \mathcal +//! E[4]\\), \\(\mathcal E[4] \cong \mathbb Z / 4 +//! \\) and \\( \mathcal E[2] \cong \mathbb Z / 2\\). In particular +//! this tells us that the group +//! $$ +//! \frac{\[2\](\mathcal E)}{\mathcal E[4]} +//! $$ +//! is well-defined and has prime order \\( (8\ell / 2) / 4 = \ell \\). +//! This is the group we will construct using Ristretto. +//! +//! ## The Isogeny +//! +//! For \\(a = \pm 1\\), we have a \\(2\\)-isogeny +//! $$ +//! \theta\_{a,d} : \mathcal J\_{a\^2, -a(a+d)/(a-d)} \longrightarrow \mathcal E\_{a,d} +//! $$ +//! (or simply \\(\theta\\)) defined by +//! $$ +//! \theta\_{a,d} : (s,t) \mapsto \left( \frac{1}{\sqrt{ad-1}} \cdot \frac{2s}{t},\quad \frac{1+as\^2}{1-as\^2} \right). +//! $$ +//! +//! XXX Its dual is ... ? +//! +//! The kernel of the isogeny is \\( \{(0, \pm 1)\} \\). +//! The image of the isogeny is \\(\[2\](\mathcal E)\\). To see this, +//! first note that because \\( \theta \circ \hat{\theta} = [2] \\), we +//! know that \\( \[2\](\mathcal E) \subseteq \theta(\mathcal J)\\); then, to see that +//! \\(\theta(\mathcal J)\\) is exactly \\(\[2\](\mathcal E)\\), +//! recall that isogenous elliptic curves over a finite field have the +//! same number of points (exercise 5.4 of Silverman), so that +//! $$ +//! \\# \theta(\mathcal J) = \frac {\\# \mathcal J} {\\# \ker \theta} +//! = \frac {\\# \mathcal E}{2} = \\# \[2\](\mathcal E). +//! $$ +//! +//! To determine the image \\(\theta(\mathcal J[2])\\) of the +//! \\(2\\)-torsion, we consider the image of the coset \\(\theta((s,t) +//! + \mathcal J[2])\\). Let \\((x,y) = \theta(s,t)\\); then +//! \\(\theta(-s,-t) = (x,y)\\) and \\(\theta(1/as, -t/as\^2) = (-x, +//! -y)\\), so that \\(\theta(\mathcal J[2]) = \mathcal E[2]\\). +//! +//! The Decaf paper recalls that, for a group \\( G \\) with normal +//! subgroup \\(G' \leq G\\), a group homomorphism \\( \phi : G +//! \rightarrow H \\) induces a homomorphism +//! $$ +//! \bar{\phi} : \frac G {G'} \longrightarrow \frac {\phi(G)}{\phi(G')} \leq \frac {H} {\phi(G')}, +//! $$ +//! and that the induced homomorphism \\(\bar{\phi}\\) is injective if +//! \\( \ker \phi \leq G' \\). In our context, the kernel of +//! \\(\theta\\) is \\( \\{(0, \pm 1)\\} \leq \mathcal J[2] \\), +//! so \\(\theta\\) gives an isomorphism +//! $$ +//! \frac {\mathcal J} {\mathcal J[2]} +//! \cong +//! \frac {\theta(\mathcal J)} {\theta(\mathcal J[2])} +//! \cong +//! \frac {\[2\](\mathcal E)} {\mathcal E[2]}. +//! $$ +//! +//! We can use the isomorphism to transfer the encoding of \\(\mathcal +//! J / \mathcal J[2] \\) defined above to \\(\[2\](\mathcal E)/\mathcal +//! E[2]\\), by encoding the Edwards point \\((x,y)\\) using the Jacobi +//! quartic encoding of \\(\theta\^{-1}(x,y)\\). +//! +//! Since \\(\\# (\[2\](\mathcal E) / \mathcal E[2]) = (\\#\mathcal +//! E)/4\\), if \\(\mathcal E\\) has cofactor \\(4\\), we're done. +//! Otherwise, if \\(\mathcal E\\) has cofactor \\(8\\), as in the +//! Curve25519 case, we use the torquing procedure to lift \\(\mathcal E +//! / \mathcal E[4]\\) to \\(\mathcal E / \mathcal E[2]\\), and then +//! apply the encoding for \\( \[2\](\mathcal E) / \mathcal E[2] \\). +//! +//! ## The Ristretto Encoding +//! +//! We can write the above encoding/decoding procedure concretely (in affine +//! coordinates) as follows: +//! +//! ### Encoding +//! +//! On input \\( (x,y) \in \[2\](\mathcal E)\\), a representative for a +//! coset in \\( \[2\](\mathcal E) / \mathcal E[4] \\): +//! +//! 1. Check if \\( xy \\) is negative or \\( x = 0 \\); if so, torque +//! the point by setting \\( (x,y) \gets (x,y) + P_4 \\), where +//! \\(P_4\\) is a \\(4\\)-torsion point. +//! +//! 2. Check if \\(x\\) is negative or \\( y = -1 \\); if so, set +//! \\( (x,y) \gets (x,y) + (0,-1) = (-x, -y) \\). +//! +//! 3. Compute $$ s = +\sqrt {(-a) \frac {1 - y} {1 + y} }, $$ choosing +//! the positive square root. +//! +//! The output is then the (canonical) byte-encoding of \\(s\\). +//! +//! If \\(\mathcal E\\) has cofactor \\(4\\), we skip the first step, +//! since our input already represents a coset in +//! \\( \[2\](\mathcal E) / \mathcal E[2] \\). +//! +//! To see that this corresponds to the encoding procedure above, notice +//! that the first step lifts from \\( \mathcal E / \mathcal E[4] \\) to +//! \\(\mathcal E / \mathcal E[2]\\). To understand steps 2 and 3, +//! notice that the \\(y\\)-coordinate of \\(\theta(s,t)\\) is +//! $$ +//! y = \frac {1 + as\^2}{1 - as\^2}, +//! $$ +//! so that the \\(s\\)-coordinate of \\(\theta\^{-1}(x,y)\\) has +//! $$ +//! s\^2 = (-a)\frac {1-y}{1+y}. +//! $$ +//! Since +//! $$ +//! x = \frac 1 {\sqrt {ad - 1}} \frac {2s} {t}, +//! $$ +//! we also have +//! $$ +//! \frac s t = x \frac {\sqrt {ad-1}} 2, +//! $$ +//! so that the sign of \\(s/t\\) is determined by the sign of \\(x\\). +//! +//! Recall that to choose a canonical representative of \\( (s,t) + +//! \mathcal J[2] \\), it's sufficient to make two sign choices: the +//! sign of \\(s\\) and the sign of \\(s/t\\). Step 2 determines the +//! sign of \\(s/t\\), while step 3 computes \\(s\\) and determines its +//! sign (by choosing the positive square root). Finally, the check +//! that \\(y \neq -1\\) prevents division-by-zero when encoding the +//! identity; it falls out of the optimized formulas below. +//! +//! ### Decoding +//! +//! On input `s_bytes`, decoding proceeds as follows: +//! +//! 1. Decode `s_bytes` to \\(s\\); reject if `s_bytes` is not the +//! canonical encoding of \\(s\\). +//! +//! 2. Check whether \\(s\\) is negative; if so, reject. +//! +//! 3. Compute +//! $$ +//! y \gets \frac {1 + as\^2}{1 - as\^2}. +//! $$ +//! +//! 4. Compute +//! $$ +//! x \gets +\sqrt{ \frac{4s\^2} {ad(1+as\^2)\^2 - (1-as\^2)\^2}}, +//! $$ +//! choosing the positive square root, or reject if the square root does +//! not exist. +//! +//! 5. Check whether \\(xy\\) is negative or \\(y = 0\\); if so, reject. +//! +//! ## Encoding in Extended Coordinates +//! +//! The formulas above are given in affine coordinates, but the usual +//! internal representation is extended twisted Edwards coordinates \\( +//! (X:Y:Z:T) \\) with \\( x = X/Z \\), \\(y = Y/Z\\), \\(xy = T/Z \\). +//! Selecting the distinguished representative of the coset +//! requires the affine coordinates \\( (x,y) \\), and computing \\( s +//! \\) requires an inverse square root. +//! As inversions are expensive, we'd like to be able to do this +//! whole computation with only one inverse square root, by batching +//! together the inversion and the inverse square root. +//! +//! However, it is not obvious how to do this, since the inverse square +//! root computation depends on the affine coordinates (which select the +//! distinguished representative). +//! +//! In what follows we consider only the case +//! \\(a = -1\\); a similar argument applies to the case \\( a = 1\\). +//! +//! Since \\(y = Y/Z\\), in extended coordinates the formula for \\(s\\) becomes +//! $$ +//! s = \sqrt{ \frac{ 1 - Y/Z}{1+Y/Z}} = \sqrt{\frac{Z - Y}{Z+Y}} +//! = \frac {Z - Y} {\sqrt{Z\^2 - Y\^2}}. +//! $$ +//! +//! Here \\( (X:Y:Z:T) \\) are the coordinates of the distinguished +//! representative of the coset. +//! Write \\( (X\_0 : Y\_0 : Z\_0 : T\_0) \\) +//! for the coordinates of the initial representative. Then the +//! torquing procedure in step 1 replaces \\( (X\_0 : Y\_0 : Z\_0 : +//! T\_0) \\) by \\( (iY\_0 : iX\_0 : Z\_0 : -T\_0) \\). This means we +//! want to obtain either +//! $$ +//! \frac {1} { \sqrt{Z\_0\^2 - Y\_0\^2}} +//! \quad \text{or} \quad +//! \frac {1} { \sqrt{Z\_0\^2 + X\_0\^2}}. +//! $$ +//! +//! We can relate these using the identity +//! $$ +//! (a-d)X\^2Y\^2 = (Z\^2 - aX\^2)(Z\^2 - Y\^2), +//! $$ +//! which is valid for all curve points. To see this, recall from the curve equation that +//! $$ +//! -dX\^2Y\^2 = Z\^4 - aZ\^2X\^2 - Z\^2Y\^2, +//! $$ +//! so that +//! $$ +//! (a-d)X\^2Y\^2 = Z\^4 - aZ\^2X\^2 - Z\^2Y\^2 + aX\^2Y\^2 = (Z\^2 - Y\^2)(Z\^2 + X\^2). +//! $$ +//! +//! The encoding procedure is as follows: +//! +//! 1. \\(u\_1 \gets (Z\_0 + Y\_0)(Z\_0 - Y\_0) = Z\_0\^2 - Y\_0\^2 \\) +//! 2. \\(u\_2 \gets X\_0 Y\_0 \\) +//! 3. \\(I \gets \mathrm{invsqrt}(u\_1 u\_2\^2) = 1/\sqrt{X\_0\^2 Y\_0\^2 (Z\_0\^2 - Y\_0\^2)} \\) +//! 4. \\(D\_1 \gets u\_1 I = \sqrt{(Z\_0\^2 - Y\_0\^2)/(X\_0\^2 Y\_0\^2)} \\) +//! 5. \\(D\_2 \gets u\_2 I = \pm \sqrt{1/(Z\_0\^2 - Y\_0\^2)} \\) +//! 6. \\(Z\_{inv} \gets D\_1 D\_2 T\_0 = (u\_1 u\_2)/(u\_1 u\_2\^2) T\_0 = T\_0 / X\_0 Y\_0 = 1/Z\_0 \\) +//! 7. If \\( T\_0 Z\_{inv} = x\_0 y\_0 \\) is negative: +//! 1. \\( X \gets iY\_0 \\) +//! 2. \\( Y \gets iX\_0 \\) +//! 3. \\( D \gets D\_1 / \sqrt{a-d} = 1/\sqrt{Z\_0\^2 + X\_0\^2} \\) +//! 8. Otherwise: +//! 1. \\( X \gets X\_0 \\) +//! 2. \\( Y \gets Y\_0 \\) +//! 3. \\( D \gets D\_2 = \pm \sqrt{1/(Z\_0\^2 - Y\_0\^2)} \\) +//! 9. If \\( X Z\_{inv} = x \\) is negative, set \\( Y \gets - Y\\) +//! 10. Compute \\( s \gets (Z - Y) D = (Z - Y) / \sqrt{Z\^2 - Y\^2} \\) and return. +//! +//! ## Decoding to Extended Coordinates +//! +//! ## Equality Testing +//! +//! ## Elligator +//! +//! ## The Double-Ristretto Encoding +//! +//! It's possible to do batch encoding of \\( [2]P \\) using the dual +//! isogeny \\(\hat{\theta}\\). Defer this for now. +//! +//! ## ??? + +// We allow non snake_case names because coordinates in projective space are +// traditionally denoted by the capitalisation of their respective +// counterparts in affine space. Yeah, you heard me, rustc, I'm gonna have my +// affine and projective cakes and eat both of them too. +#![allow(non_snake_case)] + +use core::fmt::Debug; + +#[cfg(feature = "std")] +use rand::Rng; + +use digest::Digest; +use generic_array::typenum::U32; + +use constants; +use field::FieldElement; + +use core::ops::{Add, Sub, Neg}; +use core::ops::{AddAssign, SubAssign}; +use core::ops::{Mul, MulAssign}; + +use edwards; +use edwards::ExtendedPoint; +use edwards::CompletedPoint; +use edwards::EdwardsBasepointTable; +use edwards::Identity; +use scalar::Scalar; + +use subtle; +use subtle::ConditionallyAssignable; +use subtle::ConditionallyNegatable; +use subtle::Equal; + +// ------------------------------------------------------------------------ +// Compressed points +// ------------------------------------------------------------------------ + +/// A point serialized using Mike Hamburg's Ristretto scheme. +/// +/// XXX think about how this API should work +#[derive(Copy, Clone, Eq, PartialEq)] +pub struct CompressedRistretto(pub [u8; 32]); + +/// The result of compressing a `RistrettoPoint`. +impl CompressedRistretto { + /// View this `CompressedRistretto` as an array of bytes. + pub fn as_bytes<'a>(&'a self) -> &'a [u8; 32] { + &self.0 + } + + /// Attempt to decompress to an `RistrettoPoint`. + /// + /// This function executes in constant time for all valid inputs. + /// Inputs which do not decode to a RistrettoPoint may return + /// early. + pub fn decompress(&self) -> Option { + // Step 1. Check s for validity: + // 1.a) s must be 32 bytes (we get this from the type system) + // 1.b) s < p + // 1.c) s is nonnegative + // + // Our decoding routine ignores the high bit, so the only + // possible failure for 1.b) is if someone encodes s in 0..18 + // as s+p in 2^255-19..2^255-1. We can check this by + // converting back to bytes, and checking that we get the + // original input, since our encoding routine is canonical. + + let s = FieldElement::from_bytes(self.as_bytes()); + let s_bytes_check = s.to_bytes(); + let s_encoding_is_canonical = + subtle::slices_equal(&s_bytes_check[..], self.as_bytes()); + let s_is_negative = s.is_negative(); + + if s_encoding_is_canonical == 0u8 || s_is_negative == 1u8 { + return None; + } + + // Step 2. The rest. (XXX write comments) + let one = FieldElement::one(); + let ss = s.square(); + let yden = &one + &ss; // 1 - a*s^2 + let ynum = &one - &ss; // 1 + a*s^2 + let yden_sqr = yden.square(); + let xden_sqr = &(&(-&constants::EDWARDS_D) * &ynum.square()) - &yden_sqr; + + let (ok, invsqrt) = (&xden_sqr * &yden_sqr).invsqrt(); + + let xden_inv = &invsqrt * &yden; + let yden_inv = &invsqrt * &(&xden_inv * &xden_sqr); + + let mut x = &(&s + &s) * &xden_inv; // 2*s*xden_inv + let x_is_negative = x.is_negative(); + x.conditional_negate(x_is_negative); + let y = &ynum * &yden_inv; + + let t = &x * &y; + + if ok == 0u8 || t.is_negative() == 1u8 || y.is_zero() == 1u8 { + return None; + } else { + return Some(RistrettoPoint(ExtendedPoint{X: x, Y: y, Z: one, T: t})); + } + } +} + +impl Identity for CompressedRistretto { + fn identity() -> CompressedRistretto { + CompressedRistretto([0u8; 32]) + } +} + +// ------------------------------------------------------------------------ +// Serde support +// ------------------------------------------------------------------------ +// Serializes to and from `RistrettoPoint` directly, doing compression +// and decompression internally. This means that users can create +// structs containing `RistrettoPoint`s and use Serde's derived +// serializers to serialize those structures. + +#[cfg(feature = "serde")] +use serde::{self, Serialize, Deserialize, Serializer, Deserializer}; +#[cfg(feature = "serde")] +use serde::de::Visitor; + +#[cfg(feature = "serde")] +impl Serialize for RistrettoPoint { + fn serialize(&self, serializer: S) -> Result + where S: Serializer + { + serializer.serialize_bytes(self.compress().as_bytes()) + } +} + +#[cfg(feature = "serde")] +impl<'de> Deserialize<'de> for RistrettoPoint { + fn deserialize(deserializer: D) -> Result + where D: Deserializer<'de> + { + struct RistrettoPointVisitor; + + impl<'de> Visitor<'de> for RistrettoPointVisitor { + type Value = RistrettoPoint; + + fn expecting(&self, formatter: &mut ::core::fmt::Formatter) -> ::core::fmt::Result { + formatter.write_str("a valid point in Ristretto format") + } + + fn visit_bytes(self, v: &[u8]) -> Result + where E: serde::de::Error + { + if v.len() == 32 { + let arr32 = array_ref!(v, 0, 32); // &[u8;32] from &[u8] + CompressedRistretto(*arr32) + .decompress() + .ok_or(serde::de::Error::custom("decompression failed")) + } else { + Err(serde::de::Error::invalid_length(v.len(), &self)) + } + } + } + + deserializer.deserialize_bytes(RistrettoPointVisitor) + } +} + +// ------------------------------------------------------------------------ +// Internal point representations +// ------------------------------------------------------------------------ + +/// A `RistrettoPoint` represents a point in the Ristretto group for +/// Curve25519. Ristretto, a variant of Decaf, constructs a +/// prime-order group as a quotient group of a subgroup of (the +/// Edwards form of) Curve25519. +/// +/// Internally, a `RistrettoPoint` is a wrapper type around +/// `ExtendedPoint`, with custom equality, compression, and +/// decompression routines to account for the quotient. +#[derive(Copy, Clone)] +pub struct RistrettoPoint(pub ExtendedPoint); + +impl RistrettoPoint { + /// Compress in Ristretto format. + /// + /// # Implementation Notes + /// + /// The Ristretto encoding is as follows, on input in affine coordinates `(x,y)`: + /// + /// 1. If `xy` is negative or `x = 0`, "rotate" the point by + /// setting `(x,y) = (iy, ix)`. + /// 2. If `x` is negative, set `(x,y) = (-x, -y)`. + /// 3. Compute `s = +sqrt((1-y)/(1+y))`. + /// 4. Return the little-endian 32-byte encoding of `s`. + /// + /// However, our input is in extended twisted Edwards coordinates + /// `(X:Y:Z:T)` with `x = X/Z`, `y = Y/Z`, `xy = T/Z` (see the + /// module-level documentation on curve representations for more + /// details). Since inversions are expensive, we'd like to be + /// able to do this whole computation with only one inversion. + /// + /// Since `y = Y/Z`, in extended coordinates the formula for `s` becomes + /// + ///     s = sqrt((1 - Y/Z)/(1 + Y/Z)) = sqrt((Z-Y)/(Z+Y)). (1) + /// + /// We can compute this as + /// + ///     s = (Z - Y) / sqrt((Z-Y)(Z+Y)). (1) + /// + /// The denominator is + /// + ///     invsqrt((Z-Y)(Z+Y)) = invsqrt(Z² - Y²). (1) + /// + /// Write the input point as `(X₀:Y₀:Z₀:T₀)`. The rotation in + /// step 1 of the encoding procedure replaces `(X₀:Y₀:Z₀:T₀)` by + /// `(iY₀:iX₀:Z₀:-T₀)`. We therefore wish to relate the + /// computation of + /// + ///     invsqrt(Z² - Y²) = invsqrt(Z₀² - Y₀²) [non-rotated case] + /// + /// with the computation of + /// + ///     invsqrt(Z² - Y²) = invsqrt(Z₀² + X₀²). [rotated case] + /// + /// Recall the curve equation (in the 𝗣² model): + /// + ///     (-X² + Y²)Z² = Z⁴ + dX²Y². (1) + /// + /// This means that, for any point `(X:Y:Z:T)` in extended coordinates, we have + /// + ///     -dX²Y² = Z⁴ + Z²X² - Z²Y², (2) + /// + /// so that + /// + ///     (-1-d)X²Y² = Z⁴ + Z²X² - Z²Y² - X²Y², (3) + /// + /// and hence + /// + ///     (-1-d)X²Y² = (Z² - Y²)(Z² + X²). (4) + /// + /// Taking inverse square roots gives + /// + ///     invsqrt(Z² + X²) = invsqrt(-1-d) sqrt((Z² - Y²)/(X²Y²)). (4) + /// + /// + pub fn compress(&self) -> CompressedRistretto { + let mut X = self.0.X; + let mut Y = self.0.Y; + let Z = &self.0.Z; + let T = &self.0.T; + + let u1 = &(Z + &Y) * &(Z - &Y); + let u2 = &X * &Y; + // Ignore return value since this is always square + let (_, invsqrt) = (&u1 * &u2.square()).invsqrt(); + let i1 = &invsqrt * &u1; + let i2 = &invsqrt * &u2; + let z_inv = &i1 * &(&i2 * T); + let mut den_inv = i2; + + let iX = &X * &constants::SQRT_M1; + let iY = &Y * &constants::SQRT_M1; + let ristretto_magic = &constants::INVSQRT_A_MINUS_D; + let enchanted_denominator = &i1 * ristretto_magic; + + let rotate = (T * &z_inv).is_negative(); + + X.conditional_assign(&iY, rotate); + Y.conditional_assign(&iX, rotate); + den_inv.conditional_assign(&enchanted_denominator, rotate); + + Y.conditional_negate((&X * &z_inv).is_negative()); + + let mut s = &den_inv * &(Z - &Y); + let s_is_negative = s.is_negative(); + s.conditional_negate(s_is_negative); + + CompressedRistretto(s.to_bytes()) + } + + /// Return the coset self + E[4], for debugging. + fn coset4(&self) -> [ExtendedPoint; 4] { + [ self.0 + , &self.0 + &constants::EIGHT_TORSION[2] + , &self.0 + &constants::EIGHT_TORSION[4] + , &self.0 + &constants::EIGHT_TORSION[6] + ] + } + + /// Computes the Ristretto Elligator map. + /// + /// # Note + /// + /// This method is not public because it's just used for hashing + /// to a point -- proper elligator support is deferred for now. + pub fn elligator_ristretto_flavour(r_0: &FieldElement) -> RistrettoPoint { + let (i, d) = (&constants::SQRT_M1, &constants::EDWARDS_D); + let one = FieldElement::one(); + + let r = i * &r_0.square(); + + // D = (dr -a)(ar-d) = -(dr+1)(r+d) + let D = -&( &(&(d * &r) + &one) * &(&r + d) ); + // N = a(d-a)(d+a)(r+1) = -(r+1)(d^2 -1) + let d_sq = d.square(); + let N = -&( &(&d_sq - &one) * &(&r + &one) ); + + let mut s = FieldElement::zero(); + let mut c = -&one; + + let (N_over_D_is_square, maybe_s) = FieldElement::sqrt_ratio(&N, &D); + // s = sqrt(N/D) if N/D is square + s.conditional_assign(&maybe_s, N_over_D_is_square); + + // XXX how do we reuse the computation of sqrt(N/D) to find sqrt(rN/D) ? + let (rN_over_D_is_square, mut maybe_s) = FieldElement::sqrt_ratio(&(&r*&N), &D); + maybe_s.negate(); + + // s = -sqrt(rN/D) if rN/D is square (should happen exactly when N/D is nonsquare) + debug_assert_eq!(N_over_D_is_square ^ rN_over_D_is_square, 1u8); + s.conditional_assign(&maybe_s, rN_over_D_is_square); + c.conditional_assign(&r, rN_over_D_is_square); + + // T = (c * (r - one) * (d-one).square()) - D; + let T = &(&c * &(&(&r - &one) * &((d - &one).square()))) - &D; + + let s_sq = s.square(); + let P = CompletedPoint{ + X: &(&s + &s) * &D, + Z: &T * &constants::SQRT_AD_MINUS_ONE, + Y: &FieldElement::one() - &s_sq, + T: &FieldElement::one() + &s_sq, + }; + + // Convert to extended and return. + RistrettoPoint(P.to_extended()) + } + + /// Return a `RistrettoPoint` chosen uniformly at random using a user-provided RNG. + /// + /// # Inputs + /// + /// * `rng`: any RNG which implements the `rand::Rng` interface. + /// + /// # Returns + /// + /// A random element of the Ristretto group. + /// + /// # Implementation + /// + /// Uses the Ristretto-flavoured Elligator 2 map, so that the discrete log of the + /// output point with respect to any other point should be unknown. + #[cfg(feature = "std")] + pub fn random(rng: &mut T) -> Self { + let mut field_bytes = [0u8; 32]; + rng.fill_bytes(&mut field_bytes); + let r_0 = FieldElement::from_bytes(&field_bytes); + RistrettoPoint::elligator_ristretto_flavour(&r_0) + } + + /// Hash a slice of bytes into a `RistrettoPoint`. + /// + /// Takes a type parameter `D`, which is any `Digest` producing 32 + /// bytes (256 bits) of output. + /// + /// Convenience wrapper around `from_hash`. + /// + /// # Implementation + /// + /// Uses the Ristretto-flavoured Elligator 2 map, so that the discrete log of the + /// output point with respect to any other point should be unknown. + /// + /// # Example + /// + /// ``` + /// # extern crate curve25519_dalek; + /// # use curve25519_dalek::ristretto::RistrettoPoint; + /// extern crate sha2; + /// use sha2::Sha256; + /// + /// # // Need fn main() here in comment so the doctest compiles + /// # // See https://doc.rust-lang.org/book/documentation.html#documentation-as-tests + /// # fn main() { + /// let msg = "To really appreciate architecture, you may even need to commit a murder"; + /// let P = RistrettoPoint::hash_from_bytes::(msg.as_bytes()); + /// # } + /// ``` + /// + pub fn hash_from_bytes(input: &[u8]) -> RistrettoPoint + where D: Digest + Default + { + let mut hash = D::default(); + hash.input(input); + RistrettoPoint::from_hash(hash) + } + + /// Construct a `RistrettoPoint` from an existing `Digest` instance. + /// + /// Use this instead of `hash_from_bytes` if it is more convenient + /// to stream data into the `Digest` than to pass a single byte + /// slice. + pub fn from_hash(hash: D) -> RistrettoPoint + where D: Digest + Default + { + // XXX this seems clumsy + let mut output = [0u8; 32]; + output.copy_from_slice(hash.result().as_slice()); + let r_0 = FieldElement::from_bytes(&output); + RistrettoPoint::elligator_ristretto_flavour(&r_0) + } +} + +impl Identity for RistrettoPoint { + fn identity() -> RistrettoPoint { + RistrettoPoint(ExtendedPoint::identity()) + } +} + +// ------------------------------------------------------------------------ +// Equality +// ------------------------------------------------------------------------ + +impl PartialEq for RistrettoPoint { + fn eq(&self, other: &RistrettoPoint) -> bool { + self.ct_eq(other) == 1u8 + } +} + +impl Equal for RistrettoPoint { + /// Test equality between two `RistrettoPoint`s. + /// + /// # Returns + /// + /// `1u8` if the two `RistrettoPoint`s are equal, and `0u8` otherwise. + fn ct_eq(&self, other: &RistrettoPoint) -> u8 { + let X1Y2 = &self.0.X * &other.0.Y; + let Y1X2 = &self.0.Y * &other.0.X; + let X1X2 = &self.0.X * &other.0.X; + let Y1Y2 = &self.0.Y * &other.0.Y; + + X1Y2.ct_eq(&Y1X2) | X1X2.ct_eq(&Y1Y2) + } +} + +impl Eq for RistrettoPoint {} + +// ------------------------------------------------------------------------ +// Arithmetic +// ------------------------------------------------------------------------ + +impl<'a, 'b> Add<&'b RistrettoPoint> for &'a RistrettoPoint { + type Output = RistrettoPoint; + + fn add(self, other: &'b RistrettoPoint) -> RistrettoPoint { + RistrettoPoint(&self.0 + &other.0) + } +} + +impl<'b> AddAssign<&'b RistrettoPoint> for RistrettoPoint { + fn add_assign(&mut self, _rhs: &RistrettoPoint) { + *self = (self as &RistrettoPoint) + _rhs; + } +} + +impl<'a, 'b> Sub<&'b RistrettoPoint> for &'a RistrettoPoint { + type Output = RistrettoPoint; + + fn sub(self, other: &'b RistrettoPoint) -> RistrettoPoint { + RistrettoPoint(&self.0 - &other.0) + } +} + +impl<'b> SubAssign<&'b RistrettoPoint> for RistrettoPoint { + fn sub_assign(&mut self, _rhs: &RistrettoPoint) { + *self = (self as &RistrettoPoint) - _rhs; + } +} + +impl<'a> Neg for &'a RistrettoPoint { + type Output = RistrettoPoint; + + fn neg(self) -> RistrettoPoint { + RistrettoPoint(-&self.0) + } +} + +impl<'b> MulAssign<&'b Scalar> for RistrettoPoint { + fn mul_assign(&mut self, scalar: &'b Scalar) { + let result = (self as &RistrettoPoint) * scalar; + *self = result; + } +} + +impl<'a, 'b> Mul<&'b Scalar> for &'a RistrettoPoint { + type Output = RistrettoPoint; + /// Scalar multiplication: compute `scalar * self`. + fn mul(self, scalar: &'b Scalar) -> RistrettoPoint { + RistrettoPoint(&self.0 * scalar) + } +} + +impl<'a, 'b> Mul<&'b RistrettoPoint> for &'a Scalar { + type Output = RistrettoPoint; + + /// Scalar multiplication: compute `self * scalar`. + fn mul(self, point: &'b RistrettoPoint) -> RistrettoPoint { + RistrettoPoint(self * &point.0) + } +} + +/// Given a vector of (possibly secret) scalars and a vector of +/// (possibly secret) points, compute `c_1 P_1 + ... + c_n P_n`. +/// +/// This function has the same behaviour as +/// `vartime::multiscalar_mult` but is constant-time. +/// +/// # Input +/// +/// An iterable of `Scalar`s and a iterable of `DecafPoints`. It is an +/// error to call this function with two iterators of different lengths. +#[cfg(any(feature = "alloc", feature = "std"))] +pub fn multiscalar_mult<'a, 'b, I, J>(scalars: I, points: J) -> RistrettoPoint + where I: IntoIterator, + J: IntoIterator, +{ + let extended_points = points.into_iter().map(|P| &P.0); + RistrettoPoint(edwards::multiscalar_mult(scalars, extended_points)) +} + +/// Precomputation +#[derive(Clone)] +pub struct RistrettoBasepointTable(pub EdwardsBasepointTable); + +impl<'a, 'b> Mul<&'b Scalar> for &'a RistrettoBasepointTable { + type Output = RistrettoPoint; + + fn mul(self, scalar: &'b Scalar) -> RistrettoPoint { + RistrettoPoint(&self.0 * scalar) + } +} + +impl<'a, 'b> Mul<&'a RistrettoBasepointTable> for &'b Scalar { + type Output = RistrettoPoint; + + fn mul(self, basepoint_table: &'a RistrettoBasepointTable) -> RistrettoPoint { + RistrettoPoint(self * &basepoint_table.0) + } +} + +impl RistrettoBasepointTable { + /// Create a precomputed table of multiples of the given `basepoint`. + pub fn create(basepoint: &RistrettoPoint) -> RistrettoBasepointTable { + RistrettoBasepointTable(EdwardsBasepointTable::create(&basepoint.0)) + } + + /// Get the basepoint for this table as a `RistrettoPoint`. + pub fn basepoint(&self) -> RistrettoPoint { + RistrettoPoint(self.0.basepoint()) + } +} + +// ------------------------------------------------------------------------ +// Constant-time conditional assignment +// ------------------------------------------------------------------------ + +impl ConditionallyAssignable for RistrettoPoint { + /// Conditionally assign `other` to `self`, if `choice == 1u8`. + /// + /// # Example + /// + /// ``` + /// # extern crate subtle; + /// # extern crate curve25519_dalek; + /// # + /// # use subtle::ConditionallyAssignable; + /// # + /// # use curve25519_dalek::edwards::Identity; + /// # use curve25519_dalek::ristretto::RistrettoPoint; + /// # use curve25519_dalek::constants; + /// # fn main() { + /// let A = RistrettoPoint::identity(); + /// let B = constants::RISTRETTO_BASEPOINT_POINT; + /// + /// let mut P = A; + /// + /// P.conditional_assign(&B, 0u8); + /// assert!(P == A); + /// P.conditional_assign(&B, 1u8); + /// assert!(P == B); + /// # } + /// ``` + fn conditional_assign(&mut self, other: &RistrettoPoint, choice: u8) { + self.0.X.conditional_assign(&other.0.X, choice); + self.0.Y.conditional_assign(&other.0.Y, choice); + self.0.Z.conditional_assign(&other.0.Z, choice); + self.0.T.conditional_assign(&other.0.T, choice); + } +} + +// ------------------------------------------------------------------------ +// Debug traits +// ------------------------------------------------------------------------ + +impl Debug for CompressedRistretto { + fn fmt(&self, f: &mut ::core::fmt::Formatter) -> ::core::fmt::Result { + write!(f, "CompressedRistretto: {:?}", self.as_bytes()) + } +} + +impl Debug for RistrettoPoint { + fn fmt(&self, f: &mut ::core::fmt::Formatter) -> ::core::fmt::Result { + let coset = self.coset4(); + write!(f, "RistrettoPoint: coset \n{:?}\n{:?}\n{:?}\n{:?}", + coset[0], coset[1], coset[2], coset[3]) + } +} + +// ------------------------------------------------------------------------ +// Variable-time functions +// ------------------------------------------------------------------------ + +pub mod vartime { + //! Variable-time operations on ristretto points, useful for non-secret data. + use super::*; + + /// Given a vector of public scalars and a vector of (possibly secret) + /// points, compute + /// + /// c_1 P_1 + ... + c_n P_n. + /// + /// # Input + /// + /// A vector of `Scalar`s and a vector of `RistrettoPoints`. It is an + /// error to call this function with two vectors of different lengths. + #[cfg(any(feature = "alloc", feature = "std"))] + pub fn multiscalar_mult<'a, 'b, I, J>(scalars: I, points: J) -> RistrettoPoint + where I: IntoIterator, + J: IntoIterator + { + let extended_points = points.into_iter().map(|P| &P.0); + RistrettoPoint(edwards::vartime::multiscalar_mult(scalars, extended_points)) + } +} + +// ------------------------------------------------------------------------ +// Tests +// ------------------------------------------------------------------------ + +#[cfg(test)] +mod test { + use rand::OsRng; + + use scalar::Scalar; + use constants; + use edwards::CompressedEdwardsY; + use edwards::Identity; + use edwards::ValidityCheck; + use super::*; + + #[cfg(feature = "serde")] + use serde_cbor; + + #[test] + #[cfg(feature = "serde")] + fn serde_cbor_basepoint_roundtrip() { + let output = serde_cbor::to_vec(&constants::RISTRETTO_BASEPOINT_POINT).unwrap(); + let parsed: RistrettoPoint = serde_cbor::from_slice(&output).unwrap(); + assert_eq!(parsed, constants::RISTRETTO_BASEPOINT_POINT); + } + + #[test] + fn scalarmult_ristrettopoint_works_both_ways() { + let P = constants::RISTRETTO_BASEPOINT_POINT; + let s = Scalar::from_u64(999); + + let P1 = &P * &s; + let P2 = &s * &P; + + assert!(P1.compress().as_bytes() == P2.compress().as_bytes()); + } + + #[test] + fn decompress_negative_s_fails() { + // constants::d is neg, so decompression should fail as |d| != d. + let bad_compressed = CompressedRistretto(constants::EDWARDS_D.to_bytes()); + assert!(bad_compressed.decompress().is_none()); + } + + #[test] + fn decompress_id() { + let compressed_id = CompressedRistretto::identity(); + let id = compressed_id.decompress().unwrap(); + let mut identity_in_coset = false; + for P in &id.coset4() { + if P.compress() == CompressedEdwardsY::identity() { + identity_in_coset = true; + } + } + assert!(identity_in_coset); + } + + #[test] + fn compress_id() { + let id = RistrettoPoint::identity(); + assert_eq!(id.compress(), CompressedRistretto::identity()); + } + + #[test] + fn basepoint_roundtrip() { + let bp_compressed_ristretto = constants::RISTRETTO_BASEPOINT_POINT.compress(); + let bp_recaf = bp_compressed_ristretto.decompress().unwrap().0; + // Check that bp_recaf differs from bp by a point of order 4 + let diff = &constants::RISTRETTO_BASEPOINT_POINT.0 - &bp_recaf; + let diff4 = diff.mult_by_pow_2(2); + assert_eq!(diff4.compress(), CompressedEdwardsY::identity()); + } + + #[test] + fn encodings_of_small_multiples_of_basepoint() { + // Table of encodings of i*basepoint + // Generated using ristretto.sage + let compressed = [ + CompressedRistretto([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]), + CompressedRistretto([226, 242, 174, 10, 106, 188, 78, 113, 168, 132, 169, 97, 197, 0, 81, 95, 88, 227, 11, 106, 165, 130, 221, 141, 182, 166, 89, 69, 224, 141, 45, 118]), + CompressedRistretto([106, 73, 50, 16, 247, 73, 156, 209, 127, 236, 181, 16, 174, 12, 234, 35, 161, 16, 232, 213, 185, 1, 248, 172, 173, 211, 9, 92, 115, 163, 185, 25]), + CompressedRistretto([148, 116, 31, 93, 93, 82, 117, 94, 206, 79, 35, 240, 68, 238, 39, 213, 209, 234, 30, 43, 209, 150, 180, 98, 22, 107, 22, 21, 42, 157, 2, 89]), + CompressedRistretto([218, 128, 134, 39, 115, 53, 139, 70, 111, 250, 223, 224, 179, 41, 58, 179, 217, 253, 83, 197, 234, 108, 149, 83, 88, 245, 104, 50, 45, 175, 106, 87]), + CompressedRistretto([232, 130, 177, 49, 1, 107, 82, 193, 211, 51, 112, 128, 24, 124, 247, 104, 66, 62, 252, 203, 181, 23, 187, 73, 90, 184, 18, 196, 22, 15, 244, 78]), + CompressedRistretto([246, 71, 70, 211, 201, 43, 19, 5, 14, 216, 216, 2, 54, 167, 240, 0, 124, 59, 63, 150, 47, 91, 167, 147, 209, 154, 96, 30, 187, 29, 244, 3]), + CompressedRistretto([68, 245, 53, 32, 146, 110, 200, 31, 189, 90, 56, 120, 69, 190, 183, 223, 133, 169, 106, 36, 236, 225, 135, 56, 189, 207, 166, 167, 130, 42, 23, 109]), + CompressedRistretto([144, 50, 147, 216, 242, 40, 126, 190, 16, 226, 55, 77, 193, 165, 62, 11, 200, 135, 229, 146, 105, 159, 2, 208, 119, 213, 38, 60, 221, 85, 96, 28]), + CompressedRistretto([2, 98, 42, 206, 143, 115, 3, 163, 28, 175, 198, 63, 143, 196, 143, 220, 22, 225, 200, 200, 210, 52, 178, 240, 214, 104, 82, 130, 169, 7, 96, 49]), + CompressedRistretto([32, 112, 111, 215, 136, 178, 114, 10, 30, 210, 165, 218, 212, 149, 43, 1, 244, 19, 188, 240, 231, 86, 77, 232, 205, 200, 22, 104, 158, 45, 185, 95]), + CompressedRistretto([188, 232, 63, 139, 165, 221, 47, 165, 114, 134, 76, 36, 186, 24, 16, 249, 82, 43, 198, 0, 74, 254, 149, 135, 122, 199, 50, 65, 202, 253, 171, 66]), + CompressedRistretto([228, 84, 158, 225, 107, 154, 160, 48, 153, 202, 32, 140, 103, 173, 175, 202, 250, 76, 63, 62, 78, 83, 3, 222, 96, 38, 227, 202, 143, 248, 68, 96]), + CompressedRistretto([170, 82, 224, 0, 223, 46, 22, 245, 95, 177, 3, 47, 195, 59, 196, 39, 66, 218, 214, 189, 90, 143, 192, 190, 1, 103, 67, 108, 89, 72, 80, 31]), + CompressedRistretto([70, 55, 107, 128, 244, 9, 178, 157, 194, 181, 246, 240, 197, 37, 145, 153, 8, 150, 229, 113, 111, 65, 71, 124, 211, 0, 133, 171, 127, 16, 48, 30]), + CompressedRistretto([224, 196, 24, 247, 200, 217, 196, 205, 215, 57, 91, 147, 234, 18, 79, 58, 217, 144, 33, 187, 104, 29, 252, 51, 2, 169, 217, 154, 46, 83, 230, 78]), + ]; + let mut bp = RistrettoPoint::identity(); + for i in 0..16 { + assert_eq!(bp.compress(), compressed[i]); + bp = &bp + &constants::RISTRETTO_BASEPOINT_POINT; + } + } + + #[test] + fn four_torsion_basepoint() { + let bp = constants::RISTRETTO_BASEPOINT_POINT; + let bp_coset = bp.coset4(); + for i in 0..4 { + assert_eq!(bp, RistrettoPoint(bp_coset[i])); + } + } + + #[test] + fn four_torsion_random() { + let mut rng = OsRng::new().unwrap(); + let B = &constants::RISTRETTO_BASEPOINT_TABLE; + let P = B * &Scalar::random(&mut rng); + let P_coset = P.coset4(); + for i in 0..4 { + assert_eq!(P, RistrettoPoint(P_coset[i])); + } + } + + #[test] + fn elligator_vs_ristretto_sage() { + // Test vectors extracted from ristretto.sage. + // + // Notice that all of the byte sequences have bit 255 set to 0; this is because + // ristretto.sage does not mask the high bit of a field element. When the high bit is set, + // the ristretto.sage elligator implementation gives different results, since it takes a + // different field element as input. + let bytes: [[u8;32]; 16] = [ + [184, 249, 135, 49, 253, 123, 89, 113, 67, 160, 6, 239, 7, 105, 211, 41, 192, 249, 185, 57, 9, 102, 70, 198, 15, 127, 7, 26, 160, 102, 134, 71], + [229, 14, 241, 227, 75, 9, 118, 60, 128, 153, 226, 21, 183, 217, 91, 136, 98, 0, 231, 156, 124, 77, 82, 139, 142, 134, 164, 169, 169, 62, 250, 52], + [115, 109, 36, 220, 180, 223, 99, 6, 204, 169, 19, 29, 169, 68, 84, 23, 21, 109, 189, 149, 127, 205, 91, 102, 172, 35, 112, 35, 134, 69, 186, 34], + [16, 49, 96, 107, 171, 199, 164, 9, 129, 16, 64, 62, 241, 63, 132, 173, 209, 160, 112, 215, 105, 50, 157, 81, 253, 105, 1, 154, 229, 25, 120, 83], + [156, 131, 161, 162, 236, 251, 5, 187, 167, 171, 17, 178, 148, 210, 90, 207, 86, 21, 79, 161, 167, 215, 234, 1, 136, 242, 182, 248, 38, 85, 79, 86], + [251, 177, 124, 54, 18, 101, 75, 235, 245, 186, 19, 46, 133, 157, 229, 64, 10, 136, 181, 185, 78, 144, 254, 167, 137, 49, 107, 10, 61, 10, 21, 25], + [232, 193, 20, 68, 240, 77, 186, 77, 183, 40, 44, 86, 150, 31, 198, 212, 76, 81, 3, 217, 197, 8, 126, 128, 126, 152, 164, 208, 153, 44, 189, 77], + [173, 229, 149, 177, 37, 230, 30, 69, 61, 56, 172, 190, 219, 115, 167, 194, 71, 134, 59, 75, 28, 244, 118, 26, 162, 97, 64, 16, 15, 189, 30, 64], + [106, 71, 61, 107, 250, 117, 42, 151, 91, 202, 212, 100, 52, 188, 190, 21, 125, 218, 31, 18, 253, 241, 160, 133, 57, 242, 3, 164, 189, 68, 111, 75], + [112, 204, 182, 90, 220, 198, 120, 73, 173, 107, 193, 17, 227, 40, 162, 36, 150, 141, 235, 55, 172, 183, 12, 39, 194, 136, 43, 153, 244, 118, 91, 89], + [111, 24, 203, 123, 254, 189, 11, 162, 51, 196, 163, 136, 204, 143, 10, 222, 33, 112, 81, 205, 34, 35, 8, 66, 90, 6, 164, 58, 170, 177, 34, 25], + [225, 183, 30, 52, 236, 82, 6, 183, 109, 25, 227, 181, 25, 82, 41, 193, 80, 77, 161, 80, 242, 203, 79, 204, 136, 245, 131, 110, 237, 106, 3, 58], + [207, 246, 38, 56, 30, 86, 176, 90, 27, 200, 61, 42, 221, 27, 56, 210, 79, 178, 189, 120, 68, 193, 120, 167, 77, 185, 53, 197, 124, 128, 191, 126], + [1, 136, 215, 80, 240, 46, 63, 147, 16, 244, 230, 207, 82, 189, 74, 50, 106, 169, 138, 86, 30, 131, 214, 202, 166, 125, 251, 228, 98, 24, 36, 21], + [210, 207, 228, 56, 155, 116, 207, 54, 84, 195, 251, 215, 249, 199, 116, 75, 109, 239, 196, 251, 194, 246, 252, 228, 70, 146, 156, 35, 25, 39, 241, 4], + [34, 116, 123, 9, 8, 40, 93, 189, 9, 103, 57, 103, 66, 227, 3, 2, 157, 107, 134, 219, 202, 74, 230, 154, 78, 107, 219, 195, 214, 14, 84, 80], + ]; + let encoded_images: [CompressedRistretto; 16] = [ + CompressedRistretto([176, 157, 237, 97, 66, 29, 140, 166, 168, 94, 26, 157, 212, 216, 229, 160, 195, 246, 232, 239, 169, 112, 63, 193, 64, 32, 152, 69, 11, 190, 246, 86]), + CompressedRistretto([234, 141, 77, 203, 181, 225, 250, 74, 171, 62, 15, 118, 78, 212, 150, 19, 131, 14, 188, 238, 194, 244, 141, 138, 166, 162, 83, 122, 228, 201, 19, 26]), + CompressedRistretto([232, 231, 51, 92, 5, 168, 80, 36, 173, 179, 104, 68, 186, 149, 68, 40, 140, 170, 27, 103, 99, 140, 21, 242, 43, 62, 250, 134, 208, 255, 61, 89]), + CompressedRistretto([208, 120, 140, 129, 177, 179, 237, 159, 252, 160, 28, 13, 206, 5, 211, 241, 192, 218, 1, 97, 130, 241, 20, 169, 119, 46, 246, 29, 79, 80, 77, 84]), + CompressedRistretto([202, 11, 236, 145, 58, 12, 181, 157, 209, 6, 213, 88, 75, 147, 11, 119, 191, 139, 47, 142, 33, 36, 153, 193, 223, 183, 178, 8, 205, 120, 248, 110]), + CompressedRistretto([26, 66, 231, 67, 203, 175, 116, 130, 32, 136, 62, 253, 215, 46, 5, 214, 166, 248, 108, 237, 216, 71, 244, 173, 72, 133, 82, 6, 143, 240, 104, 41]), + CompressedRistretto([40, 157, 102, 96, 201, 223, 200, 197, 150, 181, 106, 83, 103, 126, 143, 33, 145, 230, 78, 6, 171, 146, 210, 143, 112, 5, 245, 23, 183, 138, 18, 120]), + CompressedRistretto([220, 37, 27, 203, 239, 196, 176, 131, 37, 66, 188, 243, 185, 250, 113, 23, 167, 211, 154, 243, 168, 215, 54, 171, 159, 36, 195, 81, 13, 150, 43, 43]), + CompressedRistretto([232, 121, 176, 222, 183, 196, 159, 90, 238, 193, 105, 52, 101, 167, 244, 170, 121, 114, 196, 6, 67, 152, 80, 185, 221, 7, 83, 105, 176, 208, 224, 121]), + CompressedRistretto([226, 181, 183, 52, 241, 163, 61, 179, 221, 207, 220, 73, 245, 242, 25, 236, 67, 84, 179, 222, 167, 62, 167, 182, 32, 9, 92, 30, 165, 127, 204, 68]), + CompressedRistretto([226, 119, 16, 242, 200, 139, 240, 87, 11, 222, 92, 146, 156, 243, 46, 119, 65, 59, 1, 248, 92, 183, 50, 175, 87, 40, 206, 53, 208, 220, 148, 13]), + CompressedRistretto([70, 240, 79, 112, 54, 157, 228, 146, 74, 122, 216, 88, 232, 62, 158, 13, 14, 146, 115, 117, 176, 222, 90, 225, 244, 23, 94, 190, 150, 7, 136, 96]), + CompressedRistretto([22, 71, 241, 103, 45, 193, 195, 144, 183, 101, 154, 50, 39, 68, 49, 110, 51, 44, 62, 0, 229, 113, 72, 81, 168, 29, 73, 106, 102, 40, 132, 24]), + CompressedRistretto([196, 133, 107, 11, 130, 105, 74, 33, 204, 171, 133, 221, 174, 193, 241, 36, 38, 179, 196, 107, 219, 185, 181, 253, 228, 47, 155, 42, 231, 73, 41, 78]), + CompressedRistretto([58, 255, 225, 197, 115, 208, 160, 143, 39, 197, 82, 69, 143, 235, 92, 170, 74, 40, 57, 11, 171, 227, 26, 185, 217, 207, 90, 185, 197, 190, 35, 60]), + CompressedRistretto([88, 43, 92, 118, 223, 136, 105, 145, 238, 186, 115, 8, 214, 112, 153, 253, 38, 108, 205, 230, 157, 130, 11, 66, 101, 85, 253, 110, 110, 14, 148, 112]), + ]; + for i in 0..16 { + let r_0 = FieldElement::from_bytes(&bytes[i]); + let Q = RistrettoPoint::elligator_ristretto_flavour(&r_0); + assert_eq!(Q.compress(), encoded_images[i]); + } + } + + #[test] + fn random_roundtrip() { + let mut rng = OsRng::new().unwrap(); + let B = &constants::RISTRETTO_BASEPOINT_TABLE; + for _ in 0..100 { + let P = B * &Scalar::random(&mut rng); + let compressed_P = P.compress(); + let Q = compressed_P.decompress().unwrap(); + assert_eq!(P, Q); + } + } + + #[test] + fn random_is_valid() { + let mut rng = OsRng::new().unwrap(); + for _ in 0..100 { + let P = RistrettoPoint::random(&mut rng); + // Check that P is on the curve + assert!(P.0.is_valid()); + // Check that P is in the image of the ristretto map + P.compress(); + } + } +} + +#[cfg(all(test, feature = "bench"))] +mod bench { + use rand::OsRng; + use test::Bencher; + + use super::*; + + #[bench] + fn decompression(b: &mut Bencher) { + let mut rng = OsRng::new().unwrap(); + let B = &constants::RISTRETTO_BASEPOINT_TABLE; + let P = B * &Scalar::random(&mut rng); + let P_compressed = P.compress(); + b.iter(|| P_compressed.decompress().unwrap()); + } + + #[bench] + fn compression(b: &mut Bencher) { + let mut rng = OsRng::new().unwrap(); + let B = &constants::RISTRETTO_BASEPOINT_TABLE; + let P = B * &Scalar::random(&mut rng); + b.iter(|| P.compress()); + } +} diff --git a/src/scalar.rs b/src/scalar.rs index 00124c1..c793cc3 100644 --- a/src/scalar.rs +++ b/src/scalar.rs @@ -22,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; @@ -43,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; @@ -108,50 +104,47 @@ impl IndexMut 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) } } @@ -234,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. /// @@ -384,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 /// @@ -442,127 +427,41 @@ 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 for UnpackedScalar { - type Output = i64; - - fn index(&self, _index: usize) -> &i64 { - &(self.0[_index]) - } -} - -impl IndexMut 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. @@ -571,25 +470,25 @@ impl UnpackedScalar { // https://briansmith.org/ecc-inversion-addition-chains-01#curve25519_scalar_inversion // as it was published on 2017-09-03. - let _1 = *self; - let _10 = _1.square(); - let _100 = _10.square(); - let _11 = UnpackedScalar::multiply_add(&_10, &_1, &UnpackedScalar::zero()); - let _101 = UnpackedScalar::multiply_add(&_10, &_11, &UnpackedScalar::zero()); - let _111 = UnpackedScalar::multiply_add(&_10, &_101, &UnpackedScalar::zero()); - let _1001 = UnpackedScalar::multiply_add(&_10, &_111, &UnpackedScalar::zero()); - let _1011 = UnpackedScalar::multiply_add(&_10, &_1001, &UnpackedScalar::zero()); - let _1111 = UnpackedScalar::multiply_add(&_100, &_1011, &UnpackedScalar::zero()); + 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::multiply_add(&_1111, &_1, &UnpackedScalar::zero()); + 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.square(); + *y = y.montgomery_square(); } - *y = UnpackedScalar::multiply_add(y, x, &UnpackedScalar::zero()); + *y = UnpackedScalar::montgomery_mul(y, x); } square_multiply(&mut y, 123 + 3, &_101); @@ -620,194 +519,14 @@ impl UnpackedScalar { square_multiply(&mut y, 3, &_101); square_multiply(&mut y, 1 + 2, &_11); - y - } - - /// Compute `a^2 (mod l)`. - pub fn square(&self) -> UnpackedScalar { - let a = self.0; - let mut result = [0i64; 24]; - - result[0] = a[0]*a[0]; - result[1] = 2i64 * a[0]*a[1]; - result[2] = 2i64 * (a[0]*a[2]) + a[1]*a[1]; - result[3] = 2i64 * (a[0]*a[3] + a[1]*a[2]); - result[4] = 2i64 * (a[0]*a[4] + a[1]*a[3]) + a[2]*a[2]; - result[5] = 2i64 * (a[0]*a[5] + a[1]*a[4] + a[2]*a[3]); - result[6] = 2i64 * (a[0]*a[6] + a[1]*a[5] + a[2]*a[4]) + a[3]*a[3]; - result[7] = 2i64 * (a[0]*a[7] + a[1]*a[6] + a[2]*a[5] + a[3]*a[4]); - result[8] = 2i64 * (a[0]*a[8] + a[1]*a[7] + a[2]*a[6] + a[3]*a[5]) + a[4]*a[4]; - result[9] = 2i64 * (a[0]*a[9] + a[1]*a[8] + a[2]*a[7] + a[3]*a[6] + a[4]*a[5]); - result[10] = 2i64 * (a[0]*a[10] + a[1]*a[9] + a[2]*a[8] + a[3]*a[7] + a[4]*a[6]) + a[5]*a[5]; - result[11] = 2i64 * (a[0]*a[11] + a[1]*a[10] + a[2]*a[9] + a[3]*a[8] + a[4]*a[7] + a[5]*a[6]); - result[12] = 2i64 * (a[1]*a[11] + a[2]*a[10] + a[3]*a[9] + a[4]*a[8] + a[5]*a[7]) + a[6]*a[6]; - result[13] = 2i64 * (a[2]*a[11] + a[3]*a[10] + a[4]*a[9] + a[5]*a[8] + a[6]*a[7]); - result[14] = 2i64 * (a[3]*a[11] + a[4]*a[10] + a[5]*a[9] + a[6]*a[8]) + a[7]*a[7]; - result[15] = 2i64 * (a[4]*a[11] + a[5]*a[10] + a[6]*a[9] + a[7]*a[8]); - result[16] = 2i64 * (a[5]*a[11] + a[6]*a[10] + a[7]*a[9]) + a[8]*a[8]; - result[17] = 2i64 * (a[6]*a[11] + a[7]*a[10] + a[8]*a[9]); - result[18] = 2i64 * (a[7]*a[11] + a[8]*a[10]) + a[9]*a[9]; - result[19] = 2i64 * (a[8]*a[11] + a[9]*a[10]); - result[20] = 2i64 * (a[9]*a[11]) + a[10]*a[10]; - result[21] = 2i64 * (a[10]*a[11]); - result[22] = a[11]*a[11]; - result[23] = 0i64; - - // Reduce limbs - UnpackedScalar::reduce_limbs(&mut result) - } - - /// 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; - } - - 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); - } - - 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( @@ -815,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, @@ -857,6 +582,36 @@ mod test { 0,0,0,11,0,0,0,0,0,15,0,0,0,0,0,-9,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,7, 0,0,0,0,0,-15,0,0,0,0,0,15,0,0,0,0,15,0,0,0,0,15,0,0,0,0,0,1,0,0,0,0]; + #[test] + fn fuzzer_testcase_reduction() { + // LE bytes of 24519928653854221733733552434404946937899825954937634815 + let a_bytes = [255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 0, 0, 0, 0, 0, 0, 0, 0, 0]; + // LE bytes of 4975441334397345751130612518500927154628011511324180036903450236863266160640 + let b_bytes = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 255, 210, 210, 210, 255, 255, 255, 255, 10]; + // LE bytes of 6432735165214683820902750800207468552549813371247423777071615116673864412038 + let c_bytes = [134, 171, 119, 216, 180, 128, 178, 62, 171, 132, 32, 62, 34, 119, 104, 193, 47, 215, 181, 250, 14, 207, 172, 93, 75, 207, 211, 103, 144, 204, 56, 14]; + + let a = Scalar(a_bytes); + let b = Scalar(b_bytes); + let c = Scalar(c_bytes); + + let mut tmp = [0u8; 64]; + + // also_a = (a mod l) + tmp[0..32].copy_from_slice(&a_bytes[..]); + let also_a = Scalar::reduce(&tmp); + + // also_b = (b mod l) + tmp[0..32].copy_from_slice(&b_bytes[..]); + let also_b = Scalar::reduce(&tmp); + + let expected_c = &a * &b; + let also_expected_c = &also_a * &also_b; + + assert_eq!(c, expected_c); + assert_eq!(c, also_expected_c); + } + #[test] fn non_adjacent_form() { let naf = A_SCALAR.non_adjacent_form(); @@ -900,7 +655,7 @@ mod test { #[test] fn impl_sub() { - let should_be_one = &constants::l - &constants::l_minus_1; + let should_be_one = &constants::BASEPOINT_ORDER - &constants::BASEPOINT_ORDER_MINUS_1; assert_eq!(should_be_one, Scalar::one()); } @@ -952,6 +707,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()); } @@ -966,6 +722,47 @@ mod test { assert_eq!(should_be_X, X); } + #[test] + fn to_bytes_from_bytes_roundtrips() { + let unpacked = X.unpack(); + let bytes = unpacked.to_bytes(); + let should_be_unpacked = UnpackedScalar::from_bytes(&bytes); + + assert_eq!(should_be_unpacked.0, unpacked.0); + } + + + #[test] + fn montgomery_reduce_matches_reduce() { + let mut bignum = [0u8; 64]; + + // set bignum = x + 2^256x + for i in 0..32 { + bignum[ i] = X[i]; + bignum[32+i] = X[i]; + } + // x + 2^256x (mod l) + // = 3958878930004874126169954872055634648693766179881526445624823978500314864344 + let expected = Scalar([216, 154, 179, 139, 210, 121, 2, 71, + 69, 99, 158, 216, 23, 173, 63, 100, + 204, 0, 91, 50, 219, 153, 57, 249, + 28, 82, 31, 197, 100, 165, 192, 8]); + let reduced = Scalar::reduce(&bignum); + + // The reduced scalar should match the expected + assert_eq!(reduced.0, expected.0); + + // (x + 2^256x) * R + let interim = UnpackedScalar::mul_internal(&UnpackedScalar::from_bytes_wide(&bignum), + &constants::R); + // ((x + 2^256x) * R) / R (mod l) + let montgomery_reduced = UnpackedScalar::montgomery_reduce(&interim); + + // The Montgomery reduced scalar should match the reduced one, as well as the expected + assert_eq!(montgomery_reduced.0, reduced.unpack().0); + assert_eq!(montgomery_reduced.0, expected.unpack().0) + } + #[cfg(feature = "serde")] use serde_cbor; @@ -984,7 +781,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) { @@ -993,28 +790,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 square(b: &mut Bencher) { - let x = X.unpack(); - b.iter(|| x.square()); - } - - #[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)); - } } diff --git a/src/scalar_32bit.rs b/src/scalar_32bit.rs new file mode 100644 index 0000000..238ee09 --- /dev/null +++ b/src/scalar_32bit.rs @@ -0,0 +1,555 @@ +//! 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. + +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 for Scalar32 { + type Output = u32; + fn index(&self, _index: usize) -> &u32 { + &(self.0[_index]) + } +} + +impl IndexMut 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 / 256 bit scalar into 9 29-bit limbs. + 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 << 24) - 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 ; + + 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)] + pub (crate) 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)] + pub (crate) 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)); + } +} diff --git a/src/scalar_64bit.rs b/src/scalar_64bit.rs new file mode 100644 index 0000000..a6d14ad --- /dev/null +++ b/src/scalar_64bit.rs @@ -0,0 +1,473 @@ +//! 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). + +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 for Scalar64 { + type Output = u64; + fn index(&self, _index: usize) -> &u64 { + &(self.0[_index]) + } +} + +impl IndexMut 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. + pub fn from_bytes(bytes: &[u8; 32]) -> Scalar64 { + let mut words = [0u64; 4]; + 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 << 48) - 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 = [0u64; 8]; + 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)] + pub (crate) fn mul_internal(a: &Scalar64, b: &Scalar64) -> [u128; 9] { + let mut z = [0u128; 9]; + + z[0] = m(a[0],b[0]); + z[1] = m(a[0],b[1]) + m(a[1],b[0]); + z[2] = m(a[0],b[2]) + m(a[1],b[1]) + m(a[2],b[0]); + z[3] = m(a[0],b[3]) + m(a[1],b[2]) + m(a[2],b[1]) + m(a[3],b[0]); + 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]); + z[5] = m(a[1],b[4]) + m(a[2],b[3]) + m(a[3],b[2]) + m(a[4],b[1]); + z[6] = m(a[2],b[4]) + m(a[3],b[3]) + m(a[4],b[2]); + z[7] = m(a[3],b[4]) + m(a[4],b[3]); + z[8] = m(a[4],b[4]); + + z + } + + /// 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)] + pub (crate) 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)); + } +} diff --git a/vendor/ristretto.sage b/vendor/ristretto.sage new file mode 100644 index 0000000..f40bebd --- /dev/null +++ b/vendor/ristretto.sage @@ -0,0 +1,611 @@ +import binascii +class InvalidEncodingException(Exception): pass +class NotOnCurveException(Exception): pass +class SpecException(Exception): pass + +def lobit(x): return int(x) & 1 +def hibit(x): return lobit(2*x) +def negative(x): return lobit(x) +def enc_le(x,n): return bytearray([int(x)>>(8*i) & 0xFF for i in xrange(n)]) +def dec_le(x): return sum(b<<(8*i) for i,b in enumerate(x)) +def randombytes(n): return bytearray([randint(0,255) for _ in range(n)]) + +def optimized_version_of(spec): + """Decorator: This function is an optimized version of some specification""" + def decorator(f): + def wrapper(self,*args,**kwargs): + def pr(x): + if isinstance(x,bytearray): return binascii.hexlify(x) + else: return str(x) + try: spec_ans = getattr(self,spec,spec)(*args,**kwargs),None + except Exception as e: spec_ans = None,e + try: opt_ans = f(self,*args,**kwargs),None + except Exception as e: opt_ans = None,e + if spec_ans[1] is None and opt_ans[1] is not None: + raise + #raise SpecException("Mismatch in %s: spec returned %s but opt threw %s" + # % (f.__name__,str(spec_ans[0]),str(opt_ans[1]))) + if spec_ans[1] is not None and opt_ans[1] is None: + raise + #raise SpecException("Mismatch in %s: spec threw %s but opt returned %s" + # % (f.__name__,str(spec_ans[1]),str(opt_ans[0]))) + if spec_ans[0] != opt_ans[0]: + raise SpecException("Mismatch in %s: %s != %s" + % (f.__name__,pr(spec_ans[0]),pr(opt_ans[0]))) + if opt_ans[1] is not None: raise + else: return opt_ans[0] + wrapper.__name__ = f.__name__ + return wrapper + return decorator + +def xsqrt(x,exn=InvalidEncodingException("Not on curve")): + """Return sqrt(x)""" + if not is_square(x): raise exn + s = sqrt(x) + if negative(s): s=-s + return s + +def isqrt(x,exn=InvalidEncodingException("Not on curve")): + """Return 1/sqrt(x)""" + if x==0: return 0 + if not is_square(x): raise exn + return 1/sqrt(x) + +def isqrt_i(x): + """Return 1/sqrt(x) or 1/sqrt(zeta * x)""" + if x==0: return True,0 + gen = x.parent(-1) + while is_square(gen): gen = sqrt(gen) + if is_square(x): return True,1/sqrt(x) + else: return False,1/sqrt(x*gen) + +class QuotientEdwardsPoint(object): + """Abstract class for point an a quotiented Edwards curve; needs F,a,d,cofactor to work""" + def __init__(self,x=0,y=1): + x = self.x = self.F(x) + y = self.y = self.F(y) + if y^2 + self.a*x^2 != 1 + self.d*x^2*y^2: + raise NotOnCurveException(str(self)) + + def __repr__(self): + return "%s(0x%x,0x%x)" % (self.__class__.__name__, self.x, self.y) + + def __iter__(self): + yield self.x + yield self.y + + def __add__(self,other): + x,y = self + X,Y = other + a,d = self.a,self.d + return self.__class__( + (x*Y+y*X)/(1+d*x*y*X*Y), + (y*Y-a*x*X)/(1-d*x*y*X*Y) + ) + + def __neg__(self): return self.__class__(-self.x,self.y) + def __sub__(self,other): return self + (-other) + def __rmul__(self,other): return self*other + def __eq__(self,other): + """NB: this is the only method that is different from the usual one""" + x,y = self + X,Y = other + return x*Y == X*y or (self.cofactor==8 and -self.a*x*X == y*Y) + def __ne__(self,other): return not (self==other) + + def __mul__(self,exp): + exp = int(exp) + if exp < 0: exp,self = -exp,-self + total = self.__class__() + work = self + while exp != 0: + if exp & 1: total += work + work += work + exp >>= 1 + return total + + def xyzt(self): + x,y = self + z = self.F.random_element() + return x*z,y*z,z,x*y*z + + def torque(self): + """Apply cofactor group, except keeping the point even""" + if self.cofactor == 8: + if self.a == -1: return self.__class__(self.y*self.i, self.x*self.i) + if self.a == 1: return self.__class__(-self.y, self.x) + else: + return self.__class__(-self.x, -self.y) + + + # Utility functions + @classmethod + def bytesToGf(cls,bytes,mustBeProper=True,mustBePositive=False): + """Convert little-endian bytes to field element, sanity check length""" + if len(bytes) != cls.encLen: + raise InvalidEncodingException("wrong length %d" % len(bytes)) + s = dec_le(bytes) + if mustBeProper and s >= cls.F.modulus(): + raise InvalidEncodingException("%d out of range!" % s) + s = cls.F(s) + if mustBePositive and negative(s): + raise InvalidEncodingException("%d is negative!" % s) + return s + + @classmethod + def gfToBytes(cls,x,mustBePositive=False): + """Convert little-endian bytes to field element, sanity check length""" + if negative(x) and mustBePositive: x = -x + return enc_le(x,cls.encLen) + +class RistrettoPoint(QuotientEdwardsPoint): + """The new Ristretto group""" + def encodeSpec(self): + """Unoptimized specification for encoding""" + x,y = self + if self.cofactor==8 and (negative(x*y) or y==0): (x,y) = self.torque() + if y == -1: y = 1 # Avoid divide by 0; doesn't affect impl + + if negative(x): x,y = -x,-y + s = xsqrt(self.mneg*(1-y)/(1+y),exn=Exception("Unimplemented: point is odd: " + str(self))) + return self.gfToBytes(s) + + @classmethod + def decodeSpec(cls,s): + """Unoptimized specification for decoding""" + s = cls.bytesToGf(s,mustBePositive=True) + + a,d = cls.a,cls.d + x = xsqrt(4*s^2 / (a*d*(1+a*s^2)^2 - (1-a*s^2)^2)) + y = (1+a*s^2) / (1-a*s^2) + + if cls.cofactor==8 and (negative(x*y) or y==0): + raise InvalidEncodingException("x*y has high bit") + + return cls(x,y) + + @optimized_version_of("encodeSpec") + def encode(self): + """Encode, optimized version""" + a,d,mneg = self.a,self.d,self.mneg + x,y,z,t = self.xyzt() + + if self.cofactor==8: + u1 = mneg*(z+y)*(z-y) + u2 = x*y # = t*z + isr = isqrt(u1*u2^2) + i1 = isr*u1 # sqrt(mneg*(z+y)*(z-y))/(x*y) + i2 = isr*u2 # 1/sqrt(a*(y+z)*(y-z)) + z_inv = i1*i2*t # 1/z + + if negative(t*z_inv): + if a==-1: + x,y = y*self.i,x*self.i + den_inv = self.magic * i1 + else: + x,y = -y,x + den_inv = self.i * self.magic * i1 + + else: + den_inv = i2 + + if negative(x*z_inv): y = -y + s = (z-y) * den_inv + else: + num = mneg*(z+y)*(z-y) + isr = isqrt(num*y^2) + if negative(isr^2*num*y*t): y = -y + s = isr*y*(z-y) + + + return self.gfToBytes(s,mustBePositive=True) + + @classmethod + @optimized_version_of("decodeSpec") + def decode(cls,s): + """Decode, optimized version""" + s = cls.bytesToGf(s,mustBePositive=True) + + a,d = cls.a,cls.d + yden = 1-a*s^2 + ynum = 1+a*s^2 + yden_sqr = yden^2 + xden_sqr = a*d*ynum^2 - yden_sqr + + isr = isqrt(xden_sqr * yden_sqr) + + xden_inv = isr * yden + yden_inv = xden_inv * isr * xden_sqr + + x = 2*s*xden_inv + if negative(x): x = -x + y = ynum * yden_inv + + if cls.cofactor==8 and (negative(x*y) or y==0): + raise InvalidEncodingException("x*y is invalid: %d, %d" % (x,y)) + + return cls(x,y) + + @classmethod + def fromJacobiQuartic(cls,s,t,sgn=1): + """Convert point from its Jacobi Quartic representation""" + a,d = cls.a,cls.d + assert s^4 - 2*cls.a*(1-2*d/(d-a))*s^2 + 1 == t^2 + x = 2*s*cls.magic / t + y = (1+a*s^2) / (1-a*s^2) + return cls(sgn*x,y) + + @classmethod + def elligatorSpec(cls,r0): + a,d = cls.a,cls.d + r = cls.qnr * cls.bytesToGf(r0)^2 + den = (d*r-a)*(a*r-d) + n1 = cls.a*(r+1)*(a+d)*(d-a)/den + n2 = r*n1 + if is_square(n1): + sgn,s,t = 1, xsqrt(n1), -(r-1)*(a+d)^2 / den - 1 + else: + sgn,s,t = -1,-xsqrt(n2), r*(r-1)*(a+d)^2 / den - 1 + + return cls.fromJacobiQuartic(s,t) + + @classmethod + @optimized_version_of("elligatorSpec") + def elligator(cls,r0): + a,d = cls.a,cls.d + r0 = cls.bytesToGf(r0) + r = cls.qnr * r0^2 + den = (d*r-a)*(a*r-d) + num = cls.a*(r+1)*(a+d)*(d-a) + + iss,isri = isqrt_i(num*den) + if iss: sgn,twiddle = 1,1 + else: sgn,twiddle = -1,r0*cls.qnr + isri *= twiddle + s = isri*num + t = -sgn*isri*s*(r-1)*(d+a)^2 - 1 + if negative(s) == iss: s = -s + return cls.fromJacobiQuartic(s,t) + + +class Decaf_1_1_Point(QuotientEdwardsPoint): + """Like current decaf but tweaked for simplicity""" + def encodeSpec(self): + """Unoptimized specification for encoding""" + a,d = self.a,self.d + x,y = self + if x==0 or y==0: return(self.gfToBytes(0)) + + if self.cofactor==8 and negative(x*y*self.isoMagic): + x,y = self.torque() + + isr2 = isqrt(a*(y^2-1)) * sqrt(a*d-1) + + sr = xsqrt(1-a*x^2) + assert sr in [isr2*x*y,-isr2*x*y] + + altx = 1/isr2*self.isoMagic + if negative(altx): s = (1+x*y*isr2)/(a*x) + else: s = (1-x*y*isr2)/(a*x) + + return self.gfToBytes(s,mustBePositive=True) + + @classmethod + def decodeSpec(cls,s): + """Unoptimized specification for decoding""" + a,d = cls.a,cls.d + s = cls.bytesToGf(s,mustBePositive=True) + + if s==0: return cls() + isr = isqrt(s^4 + 2*(a-2*d)*s^2 + 1) + altx = 2*s*isr*cls.isoMagic + if negative(altx): isr = -isr + x = 2*s / (1+a*s^2) + y = (1-a*s^2) * isr + + if cls.cofactor==8 and (negative(x*y*cls.isoMagic) or y==0): + raise InvalidEncodingException("x*y is invalid: %d, %d" % (x,y)) + + return cls(x,y) + + @optimized_version_of("encodeSpec") + def encode(self): + """Encode, optimized version""" + a,d = self.a,self.d + x,y,z,t = self.xyzt() + + if self.cofactor == 8: + # Cofactor 8 version + num = (z+y)*(z-y) + den = x*y + tmp = isqrt(num*(a-d)*den^2) + + if negative(tmp^2*den*num*(a-d)*t^2*self.isoMagic): + den,num = num,den + tmp *= sqrt(a-d) # witness that cofactor is 8 + yisr = x*sqrt(a) + toggle = (a==1) + else: + yisr = y*(a*d-1) + toggle = False + + tiisr = tmp*num + altx = tiisr*t*self.isoMagic + if negative(altx) != toggle: tiisr =- tiisr + s = tmp*den*yisr*(tiisr*z - 1) + + else: + # Much simpler cofactor 4 version + num = (x+t)*(x-t) + isr = isqrt(num*(a-d)*x^2) + ratio = isr*num + if negative(ratio*self.isoMagic): ratio=-ratio + s = (a-d)*isr*x*(ratio*z - t) + + return self.gfToBytes(s,mustBePositive=True) + + @classmethod + @optimized_version_of("decodeSpec") + def decode(cls,s): + """Decode, optimized version""" + a,d = cls.a,cls.d + s = cls.bytesToGf(s,mustBePositive=True) + + if s==0: return cls() + s2 = s^2 + den = 1+a*s2 + num = den^2 - 4*d*s2 + isr = isqrt(num*den^2) + altx = 2*s*isr*den*cls.isoMagic + if negative(altx): isr = -isr + x = 2*s *isr^2*den*num + y = (1-a*s^2) * isr*den + + if cls.cofactor==8 and (negative(x*y*cls.isoMagic) or y==0): + raise InvalidEncodingException("x*y is invalid: %d, %d" % (x,y)) + + return cls(x,y) + + @classmethod + def fromJacobiQuartic(cls,s,t,sgn=1): + """Convert point from its Jacobi Quartic representation""" + a,d = cls.a,cls.d + if s==0: return cls() + x = 2*s / (1+a*s^2) + y = (1-a*s^2) / t + return cls(x,sgn*y) + + @classmethod + def elligatorSpec(cls,r0): + a,d = cls.a,cls.d + r = cls.qnr * cls.bytesToGf(r0)^2 + + den = (d*r-(d-a))*((d-a)*r-d) + n1 = (r+1)*(a-2*d)/den + n2 = r*n1 + if is_square(n1): + sgn,s,t = 1, xsqrt(n1), -(r-1)*(a-2*d)^2 / den - 1 + else: + sgn,s,t = -1, -xsqrt(n2), r*(r-1)*(a-2*d)^2 / den - 1 + + return cls.fromJacobiQuartic(s,t) + + @classmethod + @optimized_version_of("elligatorSpec") + def elligator(cls,r0): + a,d = cls.a,cls.d + r0 = cls.bytesToGf(r0) + r = cls.qnr * r0^2 + den = (d*r-(d-a))*((d-a)*r-d) + num = (r+1)*(a-2*d) + + iss,isri = isqrt_i(num*den) + if iss: sgn,twiddle = 1,1 + else: sgn,twiddle = -1,r0*cls.qnr + isri *= twiddle + s = isri*num + t = -sgn*isri*s*(r-1)*(a-2*d)^2 - 1 + if negative(s) == iss: s = -s + return cls.fromJacobiQuartic(s,t) + +class Ed25519Point(RistrettoPoint): + F = GF(2^255-19) + d = F(-121665/121666) + a = F(-1) + i = sqrt(F(-1)) + mneg = F(1) + qnr = i + magic = isqrt(a*d-1) + cofactor = 8 + encLen = 32 + + @classmethod + def base(cls): + return cls( 15112221349535400772501151409588531511454012693041857206046113283949847762202, 46316835694926478169428394003475163141307993866256225615783033603165251855960 + ) + +class NegEd25519Point(RistrettoPoint): + F = GF(2^255-19) + d = F(121665/121666) + a = F(1) + i = sqrt(F(-1)) + mneg = F(-1) # TODO checkme vs 1-ad or whatever + qnr = i + magic = isqrt(a*d-1) + cofactor = 8 + encLen = 32 + + @classmethod + def base(cls): + y = cls.F(4/5) + x = sqrt((y^2-1)/(cls.d*y^2-cls.a)) + if negative(x): x = -x + return cls(x,y) + +class IsoEd448Point(RistrettoPoint): + F = GF(2^448-2^224-1) + d = F(39082/39081) + a = F(1) + mneg = F(-1) + qnr = -1 + magic = isqrt(a*d-1) + cofactor = 4 + encLen = 56 + + @classmethod + def base(cls): + return cls( # RFC has it wrong + -345397493039729516374008604150537410266655260075183290216406970281645695073672344430481787759340633221708391583424041788924124567700732, + -363419362147803445274661903944002267176820680343659030140745099590306164083365386343198191849338272965044442230921818680526749009182718 + ) + +class TwistedEd448GoldilocksPoint(Decaf_1_1_Point): + F = GF(2^448-2^224-1) + d = F(-39082) + a = F(-1) + qnr = -1 + magic = isqrt(a*d-1) + cofactor = 4 + encLen = 56 + isoMagic = IsoEd448Point.magic + + @classmethod + def base(cls): + return cls.decodeSpec(Ed448GoldilocksPoint.base().encodeSpec()) + +class Ed448GoldilocksPoint(Decaf_1_1_Point): + F = GF(2^448-2^224-1) + d = F(-39081) + a = F(1) + qnr = -1 + magic = isqrt(a*d-1) + cofactor = 4 + encLen = 56 + isoMagic = IsoEd448Point.magic + + @classmethod + def base(cls): + return -2*cls( # FIXME: make not negative + 224580040295924300187604334099896036246789641632564134246125461686950415467406032909029192869357953282578032075146446173674602635247710, 298819210078481492676017930443930673437544040154080242095928241372331506189835876003536878655418784733982303233503462500531545062832660 + ) + +class IsoEd25519Point(Decaf_1_1_Point): + # TODO: twisted iso too! + # TODO: twisted iso might have to IMAGINE_TWIST or whatever + F = GF(2^255-19) + d = F(-121665) + a = F(1) + i = sqrt(F(-1)) + qnr = i + magic = isqrt(a*d-1) + cofactor = 8 + encLen = 32 + isoMagic = Ed25519Point.magic + isoA = Ed25519Point.a + + @classmethod + def base(cls): + return cls.decodeSpec(Ed25519Point.base().encode()) + +class TestFailedException(Exception): pass + +def test(cls,n): + print "Testing curve %s" % cls.__name__ + + specials = [1] + ii = cls.F(-1) + while is_square(ii): + specials.append(ii) + ii = sqrt(ii) + specials.append(ii) + for i in specials: + if negative(cls.F(i)): i = -i + i = enc_le(i,cls.encLen) + try: + Q = cls.decode(i) + QE = Q.encode() + if QE != i: + raise TestFailedException("Round trip special %s != %s" % + (binascii.hexlify(QE),binascii.hexlify(i))) + except NotOnCurveException: pass + except InvalidEncodingException: pass + + + P = cls.base() + print "base", list(P.encode()) + for i in xrange(16): + Q = P*i + print i, list(Q.encode()) + + Q = cls() + for i in xrange(n): + #print i, binascii.hexlify(Q.encode()) + QQ = cls.decode(Q.encode()) + if QQ != Q: raise TestFailedException("Round trip %s != %s" % (str(QQ),str(Q))) + + QT = Q + QE = Q.encode() + for h in xrange(cls.cofactor): + QT = QT.torque() + if QT.encode() != QE: + raise TestFailedException("Can't torque %s,%d" % (str(Q),h+1)) + + Q0 = Q + P + if Q0 == Q: raise TestFailedException("Addition doesn't work") + if Q0-P != Q: raise TestFailedException("Subtraction doesn't work") + + r = randint(1,1000) + Q1 = Q0*r + Q2 = Q0*(r+1) + if Q1 + Q0 != Q2: raise TestFailedException("Scalarmul doesn't work") + Q = Q1 + +test(Ed25519Point,100) +#test(NegEd25519Point,100) +#test(IsoEd25519Point,100) +#test(IsoEd448Point,100) +#test(TwistedEd448GoldilocksPoint,100) +#test(Ed448GoldilocksPoint,100) + + +def testElligator(cls,n): + print "Testing elligator on %s" % cls.__name__ + for i in xrange(n): + r = randombytes(cls.encLen) + Q = cls.elligator(r) + print list(r), list(Q.encode()) + +testElligator(Ed25519Point,100) +#testElligator(NegEd25519Point,100) +#testElligator(IsoEd448Point,100) +#testElligator(Ed448GoldilocksPoint,100) +#testElligator(TwistedEd448GoldilocksPoint,100) + +def gangtest(classes,n): + print "Gang test",[cls.__name__ for cls in classes] + specials = [1] + ii = classes[0].F(-1) + while is_square(ii): + specials.append(ii) + ii = sqrt(ii) + specials.append(ii) + + for i in xrange(n): + rets = [bytes((cls.base()*i).encode()) for cls in classes] + if len(set(rets)) != 1: + print "Divergence in encode at %d" % i + for c,ret in zip(classes,rets): + print c,binascii.hexlify(ret) + print + + if i < len(specials): r0 = enc_le(specials[i],classes[0].encLen) + else: r0 = randombytes(classes[0].encLen) + + rets = [bytes((cls.elligator(r0)*i).encode()) for cls in classes] + if len(set(rets)) != 1: + print "Divergence in elligator at %d" % i + for c,ret in zip(classes,rets): + print c,binascii.hexlify(ret) + print +gangtest([IsoEd448Point,TwistedEd448GoldilocksPoint,Ed448GoldilocksPoint],100) +gangtest([Ed25519Point,IsoEd25519Point],100)