From dd1ad9f1148ebf4457bb36d414e5f32c01630408 Mon Sep 17 00:00:00 2001 From: Sean Bowe Date: Sat, 22 Aug 2020 14:15:39 -0600 Subject: [PATCH] Initial commit --- .gitignore | 3 + CHANGELOG.md | 0 COPYING | 5 + Cargo.toml | 27 ++ Contributor_Agreement | 58 +++ LICENSE-TGPPL | 182 ++++++++ README.md | 19 + katex-header.html | 15 + rust-toolchain | 1 + src/arithmetic.rs | 348 +++++++++++++++ src/arithmetic/curves.rs | 860 ++++++++++++++++++++++++++++++++++++ src/arithmetic/fields.rs | 296 +++++++++++++ src/arithmetic/fields/fp.rs | 642 +++++++++++++++++++++++++++ src/arithmetic/fields/fq.rs | 656 +++++++++++++++++++++++++++ src/arithmetic/macros.rs | 153 +++++++ src/lib.rs | 19 + src/plonk.rs | 145 ++++++ src/plonk/circuit.rs | 66 +++ src/plonk/domain.rs | 202 +++++++++ src/plonk/prover.rs | 266 +++++++++++ src/plonk/srs.rs | 105 +++++ src/plonk/verifier.rs | 108 +++++ src/polycommit.rs | 600 +++++++++++++++++++++++++ src/transcript.rs | 45 ++ 24 files changed, 4821 insertions(+) create mode 100644 .gitignore create mode 100644 CHANGELOG.md create mode 100644 COPYING create mode 100644 Cargo.toml create mode 100644 Contributor_Agreement create mode 100644 LICENSE-TGPPL create mode 100644 README.md create mode 100644 katex-header.html create mode 100644 rust-toolchain create mode 100644 src/arithmetic.rs create mode 100644 src/arithmetic/curves.rs create mode 100644 src/arithmetic/fields.rs create mode 100644 src/arithmetic/fields/fp.rs create mode 100644 src/arithmetic/fields/fq.rs create mode 100644 src/arithmetic/macros.rs create mode 100644 src/lib.rs create mode 100644 src/plonk.rs create mode 100644 src/plonk/circuit.rs create mode 100644 src/plonk/domain.rs create mode 100644 src/plonk/prover.rs create mode 100644 src/plonk/srs.rs create mode 100644 src/plonk/verifier.rs create mode 100644 src/polycommit.rs create mode 100644 src/transcript.rs diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..6936990 --- /dev/null +++ b/.gitignore @@ -0,0 +1,3 @@ +/target +**/*.rs.bk +Cargo.lock diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 0000000..e69de29 diff --git a/COPYING b/COPYING new file mode 100644 index 0000000..ec6f958 --- /dev/null +++ b/COPYING @@ -0,0 +1,5 @@ +Copyright 2020 The Electric Coin Company + +You may use this package under the Transitive Grace Period Public Licence, +version 1.0, or at your option, any later version. See the file ./LICENSE-TGPPL +for the terms of the Transitive Grace Period Public Licence, version 1.0. diff --git a/Cargo.toml b/Cargo.toml new file mode 100644 index 0000000..37df7c5 --- /dev/null +++ b/Cargo.toml @@ -0,0 +1,27 @@ +[package] +name = "halo2" +version = "0.0.1" +authors = [ + "Sean Bowe ", + "Ying Tong Lai ", +] +edition = "2018" +description = """ +Fast proof-carrying data implementation with no trusted setup +""" +license = "TGPPL" +repository = "https://github.com/Electric-Coin-Company/halo2" +documentation = "https://docs.rs/halo2" +readme = "README.md" + +# We are not publishing this yet. +publish = false + +[package.metadata.docs.rs] +rustdoc-args = [ "--html-in-header", "katex-header.html" ] + +[dependencies] +subtle = "2.2.1" +crossbeam-utils = "0.7" +num_cpus = "1.13" +rand = "0.7" diff --git a/Contributor_Agreement b/Contributor_Agreement new file mode 100644 index 0000000..d296a6e --- /dev/null +++ b/Contributor_Agreement @@ -0,0 +1,58 @@ +ECC Contributor Agreement + +This Contributor Agreement ("Agreement") applies to your contribution of any +software, documentation, or other materials to any project ("project") +established or managed by ZeroCoin Electric Coin Company LLC (dba Electric Coin +Company) ("ECC", "we", "us" or "our") and establishes the intellectual property +rights you grant to ECC in your contribution. + +Please read this Agreement carefully. If you agree to this Agreement, please +complete and sign the Agreement below. We suggest you also keep a copy of this +Agreement for your records. + +Please provide your information below. If you are making a contribution +individually, the term "you" refers to the individual identified below. If you +are making a contribution on behalf of a company or other entity, "you" will +also mean the company or entity you identify below. + +Your name: + +Company/entity name (if applicable): + +Your mailing address: + +Your telephone: + +Your preferred email: + +1. The term "contribution" means any software, documentation, materials, or +other work of authorship, including any modifications or additions to any +existing software, documentation, materials, or other work, that you submit to +ECC in any form, including electronic, verbal, or written. + +2. You grant to ECC, and recipients of any contributions distributed by ECC, the +following rights and licenses: + +3. A perpetual, worldwide, non-exclusive, no-charge, royalty-free, fully +sublicensable, irrevocable right and license to reproduce, prepare derivative +works of, display, perform, sublicense, and distribute and make available any +contributions and such derivative works. + +4. A perpetual, worldwide, non-exclusive, no-charge, royalty-free, fully +sublicensable, irrevocable[ (except as stated in this section)] right and +license to make, have made, use, offer to sell, sell, import, transfer and +otherwise practice any contributions, in whole or in part, alone or in +combination with or included in any other contributions or any other software, +documentation, materials, or works. [If any entity institutes patent litigation +against you or any other entity (including a cross-claim or counterclaim in a +lawsuit) alleging that one of your contributions, or any work to which you have +contributed, constitutes direct or contributory patent infringement, then any +patent licenses granted to that entity under this Agreement for that +contribution or work shall terminate as of the date such litigation is filed.] + +You also agree that you will not assert any moral rights in any contribution +against us, our licensees or transferees. The rights in this Section are +effective on the date a contribution is first submitted to us, even if the +submission took place before the date you agree to this Agreement. Except as +stated in this Section, as between you and us, you keep all right, title, and +interest in your contribution. diff --git a/LICENSE-TGPPL b/LICENSE-TGPPL new file mode 100644 index 0000000..ce788ad --- /dev/null +++ b/LICENSE-TGPPL @@ -0,0 +1,182 @@ + +======================================================= +Transitive Grace Period Public Licence ("TGPPL") v. 1.0 +======================================================= + +This Transitive Grace Period Public Licence (the "License") applies to any +original work of authorship (the "Original Work") whose owner (the +"Licensor") has placed the following licensing notice adjacent to the +copyright notice for the Original Work: + + *Licensed under the Transitive Grace Period Public Licence version 1.0* + +1. **Grant of Copyright License.** Licensor grants You a worldwide, + royalty-free, non-exclusive, sublicensable license, for the duration of + the copyright, to do the following: + + a. to reproduce the Original Work in copies, either alone or as part of a + collective work; + + b. to translate, adapt, alter, transform, modify, or arrange the Original + Work, thereby creating derivative works ("Derivative Works") based upon + the Original Work; + + c. to distribute or communicate copies of the Original Work and Derivative + Works to the public, with the proviso that copies of Original Work or + Derivative Works that You distribute or communicate shall be licensed + under this Transitive Grace Period Public Licence no later than 12 + months after You distributed or communicated said copies; + + d. to perform the Original Work publicly; and + + e. to display the Original Work publicly. + +2. **Grant of Patent License.** Licensor grants You a worldwide, + royalty-free, non-exclusive, sublicensable license, under patent claims + owned or controlled by the Licensor that are embodied in the Original + Work as furnished by the Licensor, for the duration of the patents, to + make, use, sell, offer for sale, have made, and import the Original Work + and Derivative Works. + +3. **Grant of Source Code License.** The term "Source Code" means the + preferred form of the Original Work for making modifications to it and + all available documentation describing how to modify the Original + Work. Licensor agrees to provide a machine-readable copy of the Source + Code of the Original Work along with each copy of the Original Work that + Licensor distributes. Licensor reserves the right to satisfy this + obligation by placing a machine-readable copy of the Source Code in an + information repository reasonably calculated to permit inexpensive and + convenient access by You for as long as Licensor continues to distribute + the Original Work. + +4. **Exclusions From License Grant.** Neither the names of Licensor, nor the + names of any contributors to the Original Work, nor any of their + trademarks or service marks, may be used to endorse or promote products + derived from this Original Work without express prior permission of the + Licensor. Except as expressly stated herein, nothing in this License + grants any license to Licensor's trademarks, copyrights, patents, trade + secrets or any other intellectual property. No patent license is granted + to make, use, sell, offer for sale, have made, or import embodiments of + any patent claims other than the licensed claims defined in Section 2. No + license is granted to the trademarks of Licensor even if such marks are + included in the Original Work. Nothing in this License shall be + interpreted to prohibit Licensor from licensing under terms different + from this License any Original Work that Licensor otherwise would have a + right to license. + +5. **External Deployment.** The term "External Deployment" means the use, + distribution, or communication of the Original Work or Derivative Works + in any way such that the Original Work or Derivative Works may be used by + anyone other than You, whether those works are distributed or + communicated to those persons or made available as an application + intended for use over a network. As an express condition for the grants + of license hereunder, You must treat any External Deployment by You of + the Original Work or a Derivative Work as a distribution under section + 1(c). + +6. **Attribution Rights.** You must retain, in the Source Code of any + Derivative Works that You create, all copyright, patent, or trademark + notices from the Source Code of the Original Work, as well as any notices + of licensing and any descriptive text identified therein as an + "Attribution Notice." You must cause the Source Code for any Derivative + Works that You create to carry a prominent Attribution Notice reasonably + calculated to inform recipients that You have modified the Original Work. + +7. **Warranty of Provenance and Disclaimer of Warranty.** Licensor warrants + that the copyright in and to the Original Work and the patent rights + granted herein by Licensor are owned by the Licensor or are sublicensed + to You under the terms of this License with the permission of the + contributor(s) of those copyrights and patent rights. Except as expressly + stated in the immediately preceding sentence, the Original Work is + provided under this License on an "AS IS" BASIS and WITHOUT WARRANTY, + either express or implied, including, without limitation, the warranties + of non-infringement, merchantability or fitness for a particular + purpose. THE ENTIRE RISK AS TO THE QUALITY OF THE ORIGINAL WORK IS WITH + YOU. This DISCLAIMER OF WARRANTY constitutes an essential part of this + License. No license to the Original Work is granted by this License + except under this disclaimer. + +8. **Limitation of Liability.** Under no circumstances and under no legal + theory, whether in tort (including negligence), contract, or otherwise, + shall the Licensor be liable to anyone for any indirect, special, + incidental, or consequential damages of any character arising as a result + of this License or the use of the Original Work including, without + limitation, damages for loss of goodwill, work stoppage, computer failure + or malfunction, or any and all other commercial damages or losses. This + limitation of liability shall not apply to the extent applicable law + prohibits such limitation. + +9. **Acceptance and Termination.** If, at any time, You expressly assented + to this License, that assent indicates your clear and irrevocable + acceptance of this License and all of its terms and conditions. If You + distribute or communicate copies of the Original Work or a Derivative + Work, You must make a reasonable effort under the circumstances to obtain + the express assent of recipients to the terms of this License. This + License conditions your rights to undertake the activities listed in + Section 1, including your right to create Derivative Works based upon the + Original Work, and doing so without honoring these terms and conditions + is prohibited by copyright law and international treaty. Nothing in this + License is intended to affect copyright exceptions and limitations + (including 'fair use' or 'fair dealing'). This License shall terminate + immediately and You may no longer exercise any of the rights granted to + You by this License upon your failure to honor the conditions in Section + 1(c). + +10. **Termination for Patent Action.** This License shall terminate + automatically and You may no longer exercise any of the rights granted to + You by this License as of the date You commence an action, including a + cross-claim or counterclaim, against Licensor or any licensee alleging + that the Original Work infringes a patent. This termination provision + shall not apply for an action alleging patent infringement by + combinations of the Original Work with other software or hardware. + +11. **Jurisdiction, Venue and Governing Law.** Any action or suit relating to + this License may be brought only in the courts of a jurisdiction wherein + the Licensor resides or in which Licensor conducts its primary business, + and under the laws of that jurisdiction excluding its conflict-of-law + provisions. The application of the United Nations Convention on Contracts + for the International Sale of Goods is expressly excluded. Any use of the + Original Work outside the scope of this License or after its termination + shall be subject to the requirements and penalties of copyright or patent + law in the appropriate jurisdiction. This section shall survive the + termination of this License. + +12. **Attorneys' Fees.** In any action to enforce the terms of this License + or seeking damages relating thereto, the prevailing party shall be + entitled to recover its costs and expenses, including, without + limitation, reasonable attorneys' fees and costs incurred in connection + with such action, including any appeal of such action. This section shall + survive the termination of this License. + +13. **Miscellaneous.** If any provision of this License is held to be + unenforceable, such provision shall be reformed only to the extent + necessary to make it enforceable. + +14. **Definition of "You" in This License.** "You" throughout this License, + whether in upper or lower case, means an individual or a legal entity + exercising rights under, and complying with all of the terms of, this + License. For legal entities, "You" includes any entity that controls, is + controlled by, or is under common control with you. For purposes of this + definition, "control" means (i) the power, direct or indirect, to cause + the direction or management of such entity, whether by contract or + otherwise, or (ii) ownership of fifty percent (50%) or more of the + outstanding shares, or (iii) beneficial ownership of such entity. + +15. **Right to Use.** You may use the Original Work in all ways not otherwise + restricted or conditioned by this License or by law, and Licensor + promises not to interfere with or be responsible for such uses by You. + +16. **Modification of This License.** This License is Copyright © 2007 Zooko + Wilcox-O'Hearn. Permission is granted to copy, distribute, or communicate + this License without modification. Nothing in this License permits You to + modify this License as applied to the Original Work or to Derivative + Works. However, You may modify the text of this License and copy, + distribute or communicate your modified version (the "Modified License") + and apply it to other original works of authorship subject to the + following conditions: (i) You may not indicate in any way that your + Modified License is the "Transitive Grace Period Public Licence" or + "TGPPL" and you may not use those names in the name of your Modified + License; and (ii) You must replace the notice specified in the first + paragraph above with the notice "Licensed under " or with a notice of your own that is not confusingly similar to + the notice in this License. \ No newline at end of file diff --git a/README.md b/README.md new file mode 100644 index 0000000..85bf23c --- /dev/null +++ b/README.md @@ -0,0 +1,19 @@ +# halo2 [![Crates.io](https://img.shields.io/crates/v/halo2.svg)](https://crates.io/crates/halo2) # + +**IMPORTANT**: This library is being actively developed and should not be used in production software. + +## [Documentation](https://docs.rs/halo2) + +## License + +Copyright 2020 The Electric Coin Company. + +You may use this package under the Transitive Grace Period Public Licence, +version 1.0, or at your option, any later version. See the file +[`LICENSE-TGPPL`](LICENSE-TGPPL) for the terms of the Transitive Grace Period +Public Licence, version 1.0. + +The purpose of the TGPPL is to allow commercial improvements to the package +while ensuring that all improvements are eventually open source. See +[here](https://tahoe-lafs.org/~zooko/tgppl.pdf) for why the TGPPL exists, +graphically illustrated on three slides. diff --git a/katex-header.html b/katex-header.html new file mode 100644 index 0000000..98e8590 --- /dev/null +++ b/katex-header.html @@ -0,0 +1,15 @@ + + + + \ No newline at end of file diff --git a/rust-toolchain b/rust-toolchain new file mode 100644 index 0000000..9206354 --- /dev/null +++ b/rust-toolchain @@ -0,0 +1 @@ +1.45 \ No newline at end of file diff --git a/src/arithmetic.rs b/src/arithmetic.rs new file mode 100644 index 0000000..39d3920 --- /dev/null +++ b/src/arithmetic.rs @@ -0,0 +1,348 @@ +//! This module provides common utilities, traits and structures for group, +//! field and polynomial arithmetic. + +use crossbeam_utils::thread; + +#[macro_use] +mod macros; +mod curves; +mod fields; + +pub use curves::*; +pub use fields::*; + +/// This represents an element of a group with basic operations that can be +/// performed. This allows an FFT implementation (for example) to operate +/// generically over either a field or elliptic curve group. +pub trait Group: Copy + Clone + Send + Sync + 'static { + /// The group is assumed to be of prime order $p$. `Scalar` is the + /// associated scalar field of size $p$. + type Scalar: Field; + + /// Returns the additive identity of the group. + fn group_zero() -> Self; + + /// Adds `rhs` to this group element. + fn group_add(&mut self, rhs: &Self); + + /// Subtracts `rhs` from this group element. + fn group_sub(&mut self, rhs: &Self); + + /// Scales this group element by a scalar. + fn group_scale(&mut self, by: &Self::Scalar); +} + +/// This is a 128-bit verifier challenge. +#[derive(Copy, Clone, Debug)] +pub struct Challenge(pub(crate) u128); + +/// This algorithm applies the mapping of Algorithm 1 from the +/// [Halo](https://eprint.iacr.org/2019/1021) paper. +pub fn get_challenge_scalar(challenge: Challenge) -> F { + let mut acc = (F::ZETA + F::one()).double(); + + for i in (0..64).rev() { + let should_negate = ((challenge.0 >> ((i << 1) + 1)) & 1) == 1; + let should_endo = ((challenge.0 >> (i << 1)) & 1) == 1; + + let q = if should_negate { -F::one() } else { F::one() }; + let q = if should_endo { q * F::ZETA } else { q }; + acc = acc + q + acc; + } + + acc +} + +fn multiexp_serial(coeffs: &[C::Scalar], bases: &[C], acc: &mut C::Projective) { + let coeffs: Vec<[u8; 32]> = coeffs.iter().map(|a| a.to_bytes()).collect(); + + let c = if bases.len() < 4 { + 1 + } else if bases.len() < 32 { + 3 + } else { + (f64::from(bases.len() as u32)).ln().ceil() as usize + }; + + fn get_at(segment: usize, c: usize, bytes: &[u8; 32]) -> usize { + let skip_bits = segment * c; + let skip_bytes = skip_bits / 8; + + if skip_bytes >= 32 { + return 0; + } + + let mut v = [0; 8]; + for (v, o) in v.iter_mut().zip(bytes[skip_bytes..].iter()) { + *v = *o; + } + + let mut tmp = u64::from_le_bytes(v); + tmp >>= skip_bits - (skip_bytes * 8); + tmp = tmp % (1 << c); + + tmp as usize + } + + let segments = (256 / c) + 1; + + for current_segment in (0..segments).rev() { + for _ in 0..c { + *acc = acc.double(); + } + + #[derive(Clone, Copy)] + enum Bucket { + None, + Affine(C), + Projective(C::Projective), + } + + impl Bucket { + fn add_assign(&mut self, other: &C) { + *self = match *self { + Bucket::None => Bucket::Affine(*other), + Bucket::Affine(a) => Bucket::Projective(a + *other), + Bucket::Projective(mut a) => { + a += *other; + Bucket::Projective(a) + } + } + } + + fn add(self, mut other: C::Projective) -> C::Projective { + match self { + Bucket::None => other, + Bucket::Affine(a) => { + other += a; + other + } + Bucket::Projective(a) => other + &a, + } + } + } + + let mut buckets: Vec> = vec![Bucket::None; (1 << c) - 1]; + + for (coeff, base) in coeffs.iter().zip(bases.iter()) { + let coeff = get_at(current_segment, c, coeff); + if coeff != 0 { + buckets[coeff - 1].add_assign(base); + } + } + + // Summation by parts + // e.g. 3a + 2b + 1c = a + + // (a) + b + + // ((a) + b) + c + let mut running_sum = C::Projective::zero(); + for exp in buckets.into_iter().rev() { + running_sum = exp.add(running_sum); + *acc = *acc + &running_sum; + } + } +} + +/// Performs a multi-exponentiation operation. +/// +/// This function will panic if coeffs and bases have a different length. +/// +/// This will use multithreading if beneficial. +pub fn best_multiexp(coeffs: &[C::Scalar], bases: &[C]) -> C::Projective { + assert_eq!(coeffs.len(), bases.len()); + + let num_cpus = num_cpus::get(); + if coeffs.len() > num_cpus { + let chunk = coeffs.len() / num_cpus; + let num_chunks = coeffs.chunks(chunk).len(); + let mut results = vec![C::Projective::zero(); num_chunks]; + thread::scope(|scope| { + let chunk = coeffs.len() / num_cpus; + + for ((coeffs, bases), acc) in coeffs + .chunks(chunk) + .zip(bases.chunks(chunk)) + .zip(results.iter_mut()) + { + scope.spawn(move |_| { + multiexp_serial(coeffs, bases, acc); + }); + } + }) + .unwrap(); + results.iter().fold(C::Projective::zero(), |a, b| a + b) + } else { + let mut acc = C::Projective::zero(); + multiexp_serial(coeffs, bases, &mut acc); + acc + } +} + +/// Performs a radix-$2$ Fast-Fourier Transformation (FFT) on a vector of size +/// $n = 2^k$, when provided `log_n` = $k$ and an element of multiplicative +/// order $n$ called `omega` ($\omega$). The result is that the vector `a`, when +/// interpreted as the coefficients of a polynomial of degree $n - 1$, is +/// transformed into the evaluations of this polynomial at each of the $n$ +/// distinct powers of $\omega$. This transformation is invertible by providing +/// $\omega^{-1}$ in place of $\omega$ and dividing each resulting field element +/// by $n$. +/// +/// This will use multithreading if beneficial. +pub fn best_fft(a: &mut [G], omega: G::Scalar, log_n: u32) { + let cpus = num_cpus::get(); + let log_cpus = log2_floor(cpus); + + if log_n <= log_cpus { + serial_fft(a, omega, log_n); + } else { + parallel_fft(a, omega, log_n, log_cpus); + } +} + +fn serial_fft(a: &mut [G], omega: G::Scalar, log_n: u32) { + fn bitreverse(mut n: u32, l: u32) -> u32 { + let mut r = 0; + for _ in 0..l { + r = (r << 1) | (n & 1); + n >>= 1; + } + r + } + + let n = a.len() as u32; + assert_eq!(n, 1 << log_n); + + for k in 0..n { + let rk = bitreverse(k, log_n); + if k < rk { + a.swap(rk as usize, k as usize); + } + } + + let mut m = 1; + for _ in 0..log_n { + let w_m = omega.pow(&[u64::from(n / (2 * m)), 0, 0, 0]); + + let mut k = 0; + while k < n { + let mut w = G::Scalar::one(); + for j in 0..m { + let mut t = a[(k + j + m) as usize]; + t.group_scale(&w); + a[(k + j + m) as usize] = a[(k + j) as usize]; + a[(k + j + m) as usize].group_sub(&t); + a[(k + j) as usize].group_add(&t); + w *= &w_m; + } + + k += 2 * m; + } + + m *= 2; + } +} + +fn parallel_fft(a: &mut [G], omega: G::Scalar, log_n: u32, log_cpus: u32) { + assert!(log_n >= log_cpus); + + let num_cpus = 1 << log_cpus; + let log_new_n = log_n - log_cpus; + let mut tmp = vec![vec![G::group_zero(); 1 << log_new_n]; num_cpus]; + let new_omega = omega.pow(&[num_cpus as u64, 0, 0, 0]); + + thread::scope(|scope| { + let a = &*a; + + for (j, tmp) in tmp.iter_mut().enumerate() { + scope.spawn(move |_| { + // Shuffle into a sub-FFT + let omega_j = omega.pow(&[j as u64, 0, 0, 0]); + let omega_step = omega.pow(&[(j as u64) << log_new_n, 0, 0, 0]); + + let mut elt = G::Scalar::one(); + + for (i, tmp) in tmp.iter_mut().enumerate() { + for s in 0..num_cpus { + let idx = (i + (s << log_new_n)) % (1 << log_n); + let mut t = a[idx]; + t.group_scale(&elt); + tmp.group_add(&t); + elt *= &omega_step; + } + elt *= &omega_j; + } + + // Perform sub-FFT + serial_fft(tmp, new_omega, log_new_n); + }); + } + }) + .unwrap(); + + // Unshuffle + let mask = (1 << log_cpus) - 1; + for (idx, a) in a.iter_mut().enumerate() { + *a = tmp[idx & mask][idx >> log_cpus]; + } +} + +/// This evaluates a provided polynomial (in coefficient form) at `point`. +pub fn eval_polynomial(poly: &[F], point: F) -> F { + // TODO: parallelize? + let mut acc = F::zero(); + let mut cur = F::one(); + for coeff in poly { + acc += &(cur * coeff); + cur *= &point; + } + acc +} + +/// This computes the inner product of two vectors `a` and `b`. +/// +/// This function will panic if the two vectors are not the same size. +pub fn compute_inner_product(a: &[F], b: &[F]) -> F { + // TODO: parallelize? + assert_eq!(a.len(), b.len()); + + let mut acc = F::zero(); + for (a, b) in a.iter().zip(b.iter()) { + acc += (*a) * (*b); + } + + acc +} + +/// This simple utility function will parallelize an operation that is to be +/// performed over a mutable slice. +pub fn parallelize(v: &mut [T], f: F) { + let n = v.len(); + let num_cpus = num_cpus::get(); + let mut chunk = (n as usize) / num_cpus; + if chunk < num_cpus { + chunk = n as usize; + } + + thread::scope(|scope| { + for (chunk_num, v) in v.chunks_mut(chunk).enumerate() { + let f = f.clone(); + scope.spawn(move |_| { + let start = chunk_num * chunk; + f(v, start); + }); + } + }) + .unwrap(); +} + +fn log2_floor(num: usize) -> u32 { + assert!(num > 0); + + let mut pow = 0; + + while (1 << (pow + 1)) <= num { + pow += 1; + } + + pow +} diff --git a/src/arithmetic/curves.rs b/src/arithmetic/curves.rs new file mode 100644 index 0000000..7bb2850 --- /dev/null +++ b/src/arithmetic/curves.rs @@ -0,0 +1,860 @@ +//! This module contains implementations for the Tweedledum and Tweedledee +//! elliptic curve groups. The `Curve`/`CurveAffine` abstractions allow us to +//! write code that generalizes over these two groups. + +use core::cmp; +use core::fmt::Debug; +use core::ops::{Add, AddAssign, Mul, MulAssign, Neg, Sub, SubAssign}; +use subtle::{Choice, ConditionallySelectable, ConstantTimeEq, CtOption}; + +use super::{Field, Fp, Fq, Group}; + +/// This trait is a common interface for dealing with elements of an elliptic +/// curve group in the "projective" form, where that arithmetic is usually more +/// efficient. +pub trait Curve: + Sized + + Default + + Copy + + Clone + + Send + + Sync + + 'static + + Debug + + Add + + Sub + + Mul<::Scalar, Output = Self> + + Neg + + for<'a> Add<&'a Self, Output = Self> + + for<'a> Sub<&'a Self, Output = Self> + + MulAssign<::Scalar> + + AddAssign + + SubAssign + + for<'a> AddAssign<&'a Self> + + for<'a> SubAssign<&'a Self> + + AddAssign<::Affine> + + SubAssign<::Affine> + + PartialEq + + cmp::Eq + + ConditionallySelectable + + ConstantTimeEq + + From<::Affine> + + Group::Scalar> +{ + /// The representation of a point on this curve in the affine coordinate space. + type Affine: CurveAffine< + Projective = Self, + Scalar = ::Scalar, + Base = ::Base, + > + Add + + Sub + + Mul<::Scalar, Output = Self> + + Neg::Affine> + + From; + /// The scalar field of this elliptic curve. + type Scalar: Field; + /// The base field over which this elliptic curve is constructed. + type Base: Field; + + /// Obtains the additive identity. + fn zero() -> Self; + + /// Obtains the base point of the curve. + fn one() -> Self; + + /// Doubles this element. + fn double(&self) -> Self; + + /// Returns whether or not this element is the identity. + fn is_zero(&self) -> Choice; + + /// Apply the curve endomorphism by multiplying the x-coordinate + /// by an element of multiplicative order 3. + fn endo(&self) -> Self; + + /// Converts this element into its affine form. + fn to_affine(&self) -> Self::Affine; + + /// Returns whether or not this element is on the curve; should + /// always be true unless an "unchecked" API was used. + fn is_on_curve(&self) -> Choice; + + /// Converts many elements into their affine form. Panics if the + /// sizes of the slices are different. + fn batch_to_affine(v: &[Self], target: &mut [Self::Affine]); + + /// Returns the curve constant b + fn b() -> Self::Base; +} + +/// This trait is the affine counterpart to `Curve` and is used for +/// serialization, storage in memory, and inspection of $x$ and $y$ coordinates. +pub trait CurveAffine: + Sized + + Default + + Copy + + Clone + + Send + + Sync + + 'static + + Debug + + Add::Projective> + + Sub::Projective> + + Mul<::Scalar, Output = ::Projective> + + Neg + + PartialEq + + cmp::Eq + + ConditionallySelectable + + ConstantTimeEq + + From<::Projective> +{ + /// The representation of a point on this curve in the projective coordinate space. + type Projective: Curve< + Affine = Self, + Scalar = ::Scalar, + Base = ::Base, + > + Mul<::Scalar, Output = ::Projective> + + MulAssign<::Scalar> + + AddAssign + + SubAssign + + From; + /// The scalar field of this elliptic curve. + type Scalar: Field; + /// The base field over which this elliptic curve is constructed. + type Base: Field; + + /// Obtains the additive identity. + fn zero() -> Self; + + /// Obtains the base point of the curve. + fn one() -> Self; + + /// Returns whether or not this element is the identity. + fn is_zero(&self) -> Choice; + + /// Converts this element into its projective form. + fn to_projective(&self) -> Self::Projective; + + /// Gets the $(x, y)$ coordinates of this point. + fn get_xy(&self) -> CtOption<(Self::Base, Self::Base)>; + + /// Obtains a point given $(x, y)$, failing if it is not on the + /// curve. + fn from_xy(x: Self::Base, y: Self::Base) -> CtOption; + + /// Returns whether or not this element is on the curve; should + /// always be true unless an "unchecked" API was used. + fn is_on_curve(&self) -> Choice; + + /// Attempts to obtain a group element from its compressed 32-byte little + /// endian representation. + fn from_bytes(bytes: &[u8; 32]) -> CtOption; + + /// Obtains the compressed, 32-byte little endian representation of this + /// element. + fn to_bytes(&self) -> [u8; 32]; + + /// Attempts to obtain a group element from its uncompressed 64-byte little + /// endian representation. + fn from_bytes_wide(bytes: &[u8; 64]) -> CtOption; + + /// Obtains the uncompressed, 64-byte little endian representation of this + /// element. + fn to_bytes_wide(&self) -> [u8; 64]; + + /// Returns the curve constant $b$ + fn b() -> Self::Base; +} + +macro_rules! new_curve_impl { + ($name:ident, $name_affine:ident, $base:ident, $scalar:ident) => { + /// Represents a point in the projective coordinate space. + #[derive(Copy, Clone, Debug)] + pub struct $name { + x: $base, + y: $base, + z: $base, + } + + impl $name { + const fn curve_constant_b() -> $base { + // NOTE: this is specific to b = 5 + $base::from_raw([5, 0, 0, 0]) + } + } + + /// Represents a point in the affine coordinate space (or the point at + /// infinity). + #[derive(Copy, Clone, Debug)] + pub struct $name_affine { + x: $base, + y: $base, + infinity: Choice, + } + + impl Curve for $name { + type Affine = $name_affine; + type Scalar = $scalar; + type Base = $base; + + fn zero() -> Self { + Self { + x: $base::zero(), + y: $base::zero(), + z: $base::zero(), + } + } + + fn one() -> Self { + // NOTE: This is specific to b = 5 + + const NEGATIVE_ONE: $base = $base::neg(&$base::one()); + const TWO: $base = $base::from_raw([2, 0, 0, 0]); + + Self { + x: NEGATIVE_ONE, + y: TWO, + z: $base::one(), + } + } + + fn is_zero(&self) -> Choice { + self.z.is_zero() + } + + fn to_affine(&self) -> Self::Affine { + let zinv = self.z.invert().unwrap_or($base::zero()); + let zinv2 = zinv.square(); + let x = self.x * zinv2; + let zinv3 = zinv2 * zinv; + let y = self.y * zinv3; + + let tmp = $name_affine { + x, + y, + infinity: Choice::from(0u8), + }; + + $name_affine::conditional_select(&tmp, &$name_affine::zero(), zinv.is_zero()) + } + + fn double(&self) -> Self { + // http://www.hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-0.html#doubling-dbl-2009-l + // + // There are no points of order 2. + + let a = self.x.square(); + let b = self.y.square(); + let c = b.square(); + let d = self.x + b; + let d = d.square(); + let d = d - a - c; + let d = d + d; + let e = a + a + a; + let f = e.square(); + let z3 = self.z * self.y; + let z3 = z3 + z3; + let x3 = f - (d + d); + let c = c + c; + let c = c + c; + let c = c + c; + let y3 = e * (d - x3) - c; + + let tmp = $name { + x: x3, + y: y3, + z: z3, + }; + + $name::conditional_select(&tmp, &$name::zero(), self.is_zero()) + } + + /// Apply the curve endomorphism by multiplying the x-coordinate + /// by an element of multiplicative order 3. + fn endo(&self) -> Self { + $name { + x: self.x * $base::ZETA, + y: self.y, + z: self.z, + } + } + + fn b() -> Self::Base { + $name::curve_constant_b() + } + + fn is_on_curve(&self) -> Choice { + // Y^2 - X^3 = 5(Z^6) + + (self.y.square() - (self.x.square() * self.x)) + .ct_eq(&((self.z.square() * self.z).square() * $name::curve_constant_b())) + | self.z.is_zero() + } + + fn batch_to_affine(p: &[Self], q: &mut [Self::Affine]) { + assert_eq!(p.len(), q.len()); + + let mut acc = $base::one(); + for (p, q) in p.iter().zip(q.iter_mut()) { + // We use the `x` field of $name_affine to store the product + // of previous z-coordinates seen. + q.x = acc; + + // We will end up skipping all identities in p + acc = $base::conditional_select(&(acc * p.z), &acc, p.is_zero()); + } + + // This is the inverse, as all z-coordinates are nonzero and the ones + // that are not are skipped. + acc = acc.invert().unwrap(); + + for (p, q) in p.iter().rev().zip(q.iter_mut().rev()) { + let skip = p.is_zero(); + + // Compute tmp = 1/z + let tmp = q.x * acc; + + // Cancel out z-coordinate in denominator of `acc` + acc = $base::conditional_select(&(acc * p.z), &acc, skip); + + // Set the coordinates to the correct value + let tmp2 = tmp.square(); + let tmp3 = tmp2 * tmp; + + q.x = p.x * tmp2; + q.y = p.y * tmp3; + q.infinity = Choice::from(0u8); + + *q = $name_affine::conditional_select(&q, &$name_affine::zero(), skip); + } + } + } + + impl<'a> From<&'a $name_affine> for $name { + fn from(p: &'a $name_affine) -> $name { + p.to_projective() + } + } + + impl From<$name_affine> for $name { + fn from(p: $name_affine) -> $name { + p.to_projective() + } + } + + impl Default for $name { + fn default() -> $name { + $name::zero() + } + } + + impl ConstantTimeEq for $name { + fn ct_eq(&self, other: &Self) -> Choice { + // Is (xz^2, yz^3, z) equal to (x'z'^2, yz'^3, z') when converted to affine? + + let z = other.z.square(); + let x1 = self.x * z; + let z = z * other.z; + let y1 = self.y * z; + let z = self.z.square(); + let x2 = other.x * z; + let z = z * self.z; + let y2 = other.y * z; + + let self_is_zero = self.is_zero(); + let other_is_zero = other.is_zero(); + + (self_is_zero & other_is_zero) // Both point at infinity + | ((!self_is_zero) & (!other_is_zero) & x1.ct_eq(&x2) & y1.ct_eq(&y2)) + // Neither point at infinity, coordinates are the same + } + } + + impl PartialEq for $name { + fn eq(&self, other: &Self) -> bool { + self.ct_eq(other).into() + } + } + + impl cmp::Eq for $name {} + + impl ConditionallySelectable for $name { + fn conditional_select(a: &Self, b: &Self, choice: Choice) -> Self { + $name { + x: $base::conditional_select(&a.x, &b.x, choice), + y: $base::conditional_select(&a.y, &b.y, choice), + z: $base::conditional_select(&a.z, &b.z, choice), + } + } + } + + impl<'a> Neg for &'a $name { + type Output = $name; + + fn neg(self) -> $name { + $name { + x: self.x, + y: -self.y, + z: self.z, + } + } + } + + impl Neg for $name { + type Output = $name; + + fn neg(self) -> $name { + -&self + } + } + + impl<'a, 'b> Add<&'a $name> for &'b $name { + type Output = $name; + + fn add(self, rhs: &'a $name) -> $name { + if bool::from(self.is_zero()) { + *rhs + } else if bool::from(rhs.is_zero()) { + *self + } else { + let z1z1 = self.z.square(); + let z2z2 = rhs.z.square(); + let u1 = self.x * z2z2; + let u2 = rhs.x * z1z1; + let s1 = self.y * z2z2 * rhs.z; + let s2 = rhs.y * z1z1 * self.z; + + if u1 == u2 { + if s1 == s2 { + self.double() + } else { + $name::zero() + } + } else { + let h = u2 - u1; + let i = (h + h).square(); + let j = h * i; + let r = s2 - s1; + let r = r + r; + let v = u1 * i; + let x3 = r.square() - j - v - v; + let s1 = s1 * j; + let s1 = s1 + s1; + let y3 = r * (v - x3) - s1; + let z3 = (self.z + rhs.z).square() - z1z1 - z2z2; + let z3 = z3 * h; + + $name { + x: x3, y: y3, z: z3 + } + } + } + } + } + + impl<'a, 'b> Add<&'a $name_affine> for &'b $name { + type Output = $name; + + fn add(self, rhs: &'a $name_affine) -> $name { + if bool::from(self.is_zero()) { + rhs.to_projective() + } else if bool::from(rhs.is_zero()) { + *self + } else { + let z1z1 = self.z.square(); + let u2 = rhs.x * z1z1; + let s2 = rhs.y * z1z1 * self.z; + + if self.x == u2 { + if self.y == s2 { + self.double() + } else { + $name::zero() + } + } else { + let h = u2 - self.x; + let hh = h.square(); + let i = hh + hh; + let i = i + i; + let j = h * i; + let r = s2 - self.y; + let r = r + r; + let v = self.x * i; + let x3 = r.square() - j - v - v; + let j = self.y * j; + let j = j + j; + let y3 = r * (v - x3) - j; + let z3 = (self.z + h).square() - z1z1 - hh; + + $name { + x: x3, y: y3, z: z3 + } + } + } + } + } + + impl<'a, 'b> Sub<&'a $name> for &'b $name { + type Output = $name; + + fn sub(self, other: &'a $name) -> $name { + self + (-other) + } + } + + impl<'a, 'b> Sub<&'a $name_affine> for &'b $name { + type Output = $name; + + fn sub(self, other: &'a $name_affine) -> $name { + self + (-other) + } + } + + impl<'a, 'b> Mul<&'b $scalar> for &'a $name { + type Output = $name; + + fn mul(self, other: &'b $scalar) -> Self::Output { + // TODO: make this faster + + let mut acc = $name::zero(); + + // This is a simple double-and-add implementation of point + // multiplication, moving from most significant to least + // significant bit of the scalar. + // + // NOTE: We skip the leading bit because it's always unset. + for bit in other + .to_bytes() + .iter() + .rev() + .flat_map(|byte| (0..8).rev().map(move |i| Choice::from((byte >> i) & 1u8))) + .skip(1) + { + acc = acc.double(); + acc = $name::conditional_select(&acc, &(acc + self), bit); + } + + acc + } + } + + impl<'a> Neg for &'a $name_affine { + type Output = $name_affine; + + fn neg(self) -> $name_affine { + $name_affine { + x: self.x, + y: -self.y, + infinity: self.infinity, + } + } + } + + impl Neg for $name_affine { + type Output = $name_affine; + + fn neg(self) -> $name_affine { + -&self + } + } + + impl<'a, 'b> Add<&'a $name> for &'b $name_affine { + type Output = $name; + + fn add(self, rhs: &'a $name) -> $name { + rhs + self + } + } + + impl<'a, 'b> Add<&'a $name_affine> for &'b $name_affine { + type Output = $name; + + fn add(self, rhs: &'a $name_affine) -> $name { + if bool::from(self.is_zero()) { + rhs.to_projective() + } else if bool::from(rhs.is_zero()) { + self.to_projective() + } else { + if self.x == rhs.x { + if self.y == rhs.y { + self.to_projective().double() + } else { + $name::zero() + } + } else { + let h = rhs.x - self.x; + let hh = h.square(); + let i = hh + hh; + let i = i + i; + let j = h * i; + let r = rhs.y - self.y; + let r = r + r; + let v = self.x * i; + let x3 = r.square() - j - v - v; + let j = self.y * j; + let j = j + j; + let y3 = r * (v - x3) - j; + let z3 = h + h; + + $name { + x: x3, y: y3, z: z3 + } + } + } + } + } + + impl<'a, 'b> Sub<&'a $name_affine> for &'b $name_affine { + type Output = $name; + + fn sub(self, other: &'a $name_affine) -> $name { + self + (-other) + } + } + + impl<'a, 'b> Sub<&'a $name> for &'b $name_affine { + type Output = $name; + + fn sub(self, other: &'a $name) -> $name { + self + (-other) + } + } + + impl<'a, 'b> Mul<&'b $scalar> for &'a $name_affine { + type Output = $name; + + fn mul(self, other: &'b $scalar) -> Self::Output { + // TODO: make this faster + + let mut acc = $name::zero(); + + // This is a simple double-and-add implementation of point + // multiplication, moving from most significant to least + // significant bit of the scalar. + // + // NOTE: We skip the leading bit because it's always unset. + for bit in other + .to_bytes() + .iter() + .rev() + .flat_map(|byte| (0..8).rev().map(move |i| Choice::from((byte >> i) & 1u8))) + .skip(1) + { + acc = acc.double(); + acc = $name::conditional_select(&acc, &(acc + self), bit); + } + + acc + } + } + + impl CurveAffine for $name_affine { + type Projective = $name; + type Scalar = $scalar; + type Base = $base; + + fn zero() -> Self { + Self { + x: $base::zero(), + y: $base::zero(), + infinity: Choice::from(1u8), + } + } + + fn one() -> Self { + // NOTE: This is specific to b = 5 + + const NEGATIVE_ONE: $base = $base::neg(&$base::from_raw([1, 0, 0, 0])); + const TWO: $base = $base::from_raw([2, 0, 0, 0]); + + Self { + x: NEGATIVE_ONE, + y: TWO, + infinity: Choice::from(0u8), + } + } + + fn is_zero(&self) -> Choice { + self.infinity + } + + fn is_on_curve(&self) -> Choice { + // y^2 - x^3 ?= b + (self.y.square() - (self.x.square() * self.x)).ct_eq(&$name::curve_constant_b()) + | self.infinity + } + + fn to_projective(&self) -> Self::Projective { + $name { + x: self.x, + y: self.y, + z: $base::conditional_select(&$base::one(), &$base::zero(), self.infinity), + } + } + + fn get_xy(&self) -> CtOption<(Self::Base, Self::Base)> { + CtOption::new((self.x, self.y), !self.is_zero()) + } + + fn from_xy(x: Self::Base, y: Self::Base) -> CtOption { + let p = $name_affine { + x, y, infinity: 0u8.into() + }; + CtOption::new(p, p.is_on_curve()) + } + + fn from_bytes(bytes: &[u8; 32]) -> CtOption { + let mut tmp = *bytes; + let ysign = Choice::from(tmp[31] >> 7); + tmp[31] &= 0b0111_1111; + + $base::from_bytes(&tmp).and_then(|x| { + CtOption::new(Self::zero(), x.is_zero() & (!ysign)).or_else(|| { + let x3 = x.square() * x; + (x3 + $name::curve_constant_b()).sqrt().and_then(|y| { + let sign = Choice::from(y.to_bytes()[0] & 1); + + let y = $base::conditional_select(&y, &-y, ysign ^ sign); + + CtOption::new( + $name_affine { + x, + y, + infinity: Choice::from(0u8), + }, + Choice::from(1u8), + ) + }) + }) + }) + } + + fn to_bytes(&self) -> [u8; 32] { + // TODO: not constant time + if bool::from(self.is_zero()) { + [0; 32] + } else { + let (x, y) = (self.x, self.y); + let sign = (y.to_bytes()[0] & 1) << 7; + let mut xbytes = x.to_bytes(); + xbytes[31] |= sign; + xbytes + } + } + + fn from_bytes_wide(bytes: &[u8; 64]) -> CtOption { + let mut xbytes = [0u8; 32]; + let mut ybytes = [0u8; 32]; + xbytes.copy_from_slice(&bytes[0..32]); + ybytes.copy_from_slice(&bytes[32..64]); + + $base::from_bytes(&xbytes).and_then(|x| { + $base::from_bytes(&ybytes).and_then(|y| { + CtOption::new(Self::zero(), x.is_zero() & y.is_zero()).or_else(|| { + let on_curve = + (x * x.square() + $name::curve_constant_b()).ct_eq(&y.square()); + + CtOption::new( + $name_affine { + x, + y, + infinity: Choice::from(0u8), + }, + Choice::from(on_curve), + ) + }) + }) + }) + } + + fn to_bytes_wide(&self) -> [u8; 64] { + // TODO: not constant time + if bool::from(self.is_zero()) { + [0; 64] + } else { + let mut out = [0u8; 64]; + (&mut out[0..32]).copy_from_slice(&self.x.to_bytes()); + (&mut out[32..64]).copy_from_slice(&self.y.to_bytes()); + + out + } + } + + fn b() -> Self::Base { + $name::curve_constant_b() + } + } + + impl Default for $name_affine { + fn default() -> $name_affine { + $name_affine::zero() + } + } + + impl<'a> From<&'a $name> for $name_affine { + fn from(p: &'a $name) -> $name_affine { + p.to_affine() + } + } + + impl From<$name> for $name_affine { + fn from(p: $name) -> $name_affine { + p.to_affine() + } + } + + impl ConstantTimeEq for $name_affine { + fn ct_eq(&self, other: &Self) -> Choice { + let z1 = self.infinity; + let z2 = other.infinity; + + (z1 & z2) | ((!z1) & (!z2) & (self.x.ct_eq(&other.x)) & (self.y.ct_eq(&other.y))) + } + } + + impl PartialEq for $name_affine { + fn eq(&self, other: &Self) -> bool { + self.ct_eq(other).into() + } + } + + impl cmp::Eq for $name_affine {} + + impl ConditionallySelectable for $name_affine { + fn conditional_select(a: &Self, b: &Self, choice: Choice) -> Self { + $name_affine { + x: $base::conditional_select(&a.x, &b.x, choice), + y: $base::conditional_select(&a.y, &b.y, choice), + infinity: Choice::conditional_select(&a.infinity, &b.infinity, choice), + } + } + } + + impl_binops_additive!($name, $name); + impl_binops_additive!($name, $name_affine); + impl_binops_additive_specify_output!($name_affine, $name_affine, $name); + impl_binops_additive_specify_output!($name_affine, $name, $name); + impl_binops_multiplicative!($name, $scalar); + impl_binops_multiplicative_mixed!($name_affine, $scalar, $name); + + impl Group for $name { + type Scalar = $scalar; + + fn group_zero() -> Self { + Self::zero() + } + fn group_add(&mut self, rhs: &Self) { + *self = *self + *rhs; + } + fn group_sub(&mut self, rhs: &Self) { + *self = *self - *rhs; + } + fn group_scale(&mut self, by: &Self::Scalar) { + *self = *self * (*by); + } + } + }; +} + +new_curve_impl!(Ep, EpAffine, Fp, Fq); +new_curve_impl!(Eq, EqAffine, Fq, Fp); diff --git a/src/arithmetic/fields.rs b/src/arithmetic/fields.rs new file mode 100644 index 0000000..bc75289 --- /dev/null +++ b/src/arithmetic/fields.rs @@ -0,0 +1,296 @@ +//! This module contains implementations for the two finite fields of the +//! Tweedledum and Tweedledee curves. The `Field` abstraction allows us to write +//! code that generalizes over these two fields. + +use std::fmt::Debug; +use std::ops::{Add, AddAssign, Mul, MulAssign, Neg, Sub, SubAssign}; +use subtle::{Choice, ConditionallySelectable, ConstantTimeEq, CtOption}; + +use super::Group; + +/// This trait is a common interface for dealing with elements of a finite +/// field. +pub trait Field: + Sized + + Default + + Copy + + Clone + + Send + + Sync + + 'static + + Debug + + From + + From + + Add + + Sub + + Mul + + Neg + + for<'a> Add<&'a Self, Output = Self> + + for<'a> Mul<&'a Self, Output = Self> + + for<'a> Sub<&'a Self, Output = Self> + + MulAssign + + AddAssign + + SubAssign + + for<'a> MulAssign<&'a Self> + + for<'a> AddAssign<&'a Self> + + for<'a> SubAssign<&'a Self> + + PartialEq + + Eq + + ConditionallySelectable + + ConstantTimeEq + + Group +{ + /// How many bits needed to express the modulus $p$? + const NUM_BITS: u32; + + /// How many bits of information can be stored reliably? + const CAPACITY: u32; + + /// Represents $S$ where $p - 1 = 2^S \cdot t$ with $t$ odd. + const S: u32; + + /// Generator of the $2^S$ multiplicative subgroup + const ROOT_OF_UNITY: Self; + + /// Inverse of `ROOT_OF_UNITY` + const ROOT_OF_UNITY_INV: Self; + + /// The value $(2^S)^{-1} \mod t$. + const UNROLL_T_EXPONENT: [u64; 4]; + + /// Represents $t$ where $2^S \cdot t = p - 1$ with $t$ odd. + const T_EXPONENT: [u64; 4]; + + /// The value $t^{-1} \mod 2^S$. + const UNROLL_S_EXPONENT: u64; + + /// Inverse of $2$ in the field. + const TWO_INV: Self; + + /// Ideally the smallest prime $\alpha$ such that gcd($p - 1$, $\alpha$) = $1$ + const RESCUE_ALPHA: u64; + + /// $RESCUE_INVALPHA \cdot RESCUE_ALPHA = 1 \mod p - 1$ such that + /// `(a^RESCUE_ALPHA)^RESCUE_INVALPHA = a`. + const RESCUE_INVALPHA: [u64; 4]; + + /// Element of multiplicative order $3$. + const ZETA: Self; + + /// This computes a random element of the field using system randomness. + fn random() -> Self { + use rand::{thread_rng, RngCore}; + + let mut random_bytes = [0; 64]; + thread_rng().fill_bytes(&mut random_bytes[..]); + + Self::from_bytes_wide(&random_bytes) + } + + /// Returns whether or not this element is zero. + fn is_zero(&self) -> Choice; + + /// Doubles this element in the field. + fn double(&self) -> Self; + + /// Obtains a field element congruent to the integer `v`. + fn from_u64(v: u64) -> Self; + + /// Obtains a field element congruent to the integer `v`. + fn from_u128(v: u128) -> Self; + + /// Attempts to obtain the square root of this field element. + fn sqrt(&self) -> CtOption; + + /// Attempts to find the multiplicative inverse of this field element. + fn invert(&self) -> CtOption; + + /// Returns zero, the additive identity. + fn zero() -> Self; + + /// Returns one, the multiplicative identity. + fn one() -> Self; + + /// Squares this element in the field. + fn square(&self) -> Self; + + /// Converts this field element to its normalized, little endian byte + /// representation. + fn to_bytes(&self) -> [u8; 32]; + + /// Attempts to obtain a field element from its normalized, little endian + /// byte representation. + fn from_bytes(bytes: &[u8; 32]) -> CtOption; + + /// Obtains a field element that is congruent to the provided little endian + /// byte representation of an integer. + fn from_bytes_wide(bytes: &[u8; 64]) -> Self; + + /// Returns a square root of this element, if it exists and this element is + /// nonzero. Always returns the same square root, and it is efficient to + /// check that it has done so using `extract_radix2_vartime`. + fn deterministic_sqrt(&self) -> Option { + let sqrt = self.sqrt(); + if bool::from(sqrt.is_none()) { + return None; + } + let sqrt = sqrt.unwrap(); + let extracted = sqrt.extract_radix2_vartime()?; + + if extracted.1 >> (Self::S - 1) == 1 { + Some(-sqrt) + } else { + Some(sqrt) + } + } + + /// Returns an element $a$ of multiplicative order $t$ together with an + /// integer `s` such that `self` is the square of $a \cdot \omega^{s}$ if + /// indeed `self` is a square. + fn extract_radix2_vartime(&self) -> Option<(Self, u64)> { + if bool::from(self.is_zero()) { + return None; + } + + // TODO: these can probably be simplified + let t = self.pow_vartime(&[1 << Self::S, 0, 0, 0]); + let t = t.pow_vartime(&Self::UNROLL_T_EXPONENT); + let t = t.pow_vartime(&Self::UNROLL_T_EXPONENT); + let s = self.pow_vartime(&Self::T_EXPONENT); + let mut s = s.pow_vartime(&[Self::UNROLL_S_EXPONENT, 0, 0, 0]); + + let mut m = Self::S; + let mut c = Self::ROOT_OF_UNITY_INV; + + let mut extract: u64 = 0; + + let mut cur = 1; + while s != Self::one() { + let mut i = 1; + { + let mut s2i = s; + s2i = s2i.square(); + while s2i != Self::one() { + i += 1; + s2i = s2i.square(); + } + } + + for _ in 0..(m - i) { + c = c.square(); + cur <<= 1; + } + extract |= cur; + s *= c; + m = i; + } + + Some((t, extract)) + } + + /// Exponentiates `self` by `by`, where `by` is a little-endian order + /// integer exponent. + fn pow(&self, by: &[u64; 4]) -> Self { + let mut res = Self::one(); + for e in by.iter().rev() { + for i in (0..64).rev() { + res = res.square(); + let mut tmp = res; + tmp *= self; + res.conditional_assign(&tmp, (((*e >> i) & 0x1) as u8).into()); + } + } + res + } + + /// Exponentiates `self` by `by`, where `by` is a little-endian order + /// integer exponent. + /// + /// **This operation is variable time with respect to the exponent.** If the + /// exponent is fixed, this operation is effectively constant time. + fn pow_vartime(&self, by: &[u64; 4]) -> Self { + let mut res = Self::one(); + let mut found_one = false; + for e in by.iter().rev() { + for i in (0..64).rev() { + if found_one { + res = res.square(); + } + + if ((*e >> i) & 1) == 1 { + found_one = true; + res.mul_assign(self); + } + } + } + res + } + + /// Gets the lower 128 bits of this field element when expressed + /// canonically. + fn get_lower_128(&self) -> u128; + + /// Performs a batch inversion using Montgomery's trick, returns the product + /// of every inverse. Zero inputs are ignored. + fn batch_invert(v: &mut [Self]) -> Self { + let mut tmp = Vec::with_capacity(v.len()); + + let mut acc = Self::one(); + for p in v.iter() { + tmp.push(acc); + acc = Self::conditional_select(&(acc * p), &acc, p.is_zero()); + } + + acc = acc.invert().unwrap(); + let allinv = acc; + + for (p, tmp) in v.iter_mut().rev().zip(tmp.into_iter().rev()) { + let skip = p.is_zero(); + + let tmp = tmp * acc; + acc = Self::conditional_select(&(acc * *p), &acc, skip); + *p = Self::conditional_select(&tmp, p, skip); + } + + allinv + } +} + +mod fp; +mod fq; + +/// Compute a + b + carry, returning the result and the new carry over. +#[inline(always)] +const fn adc(a: u64, b: u64, carry: u64) -> (u64, u64) { + let ret = (a as u128) + (b as u128) + (carry as u128); + (ret as u64, (ret >> 64) as u64) +} + +/// Compute a - (b + borrow), returning the result and the new borrow. +#[inline(always)] +const fn sbb(a: u64, b: u64, borrow: u64) -> (u64, u64) { + let ret = (a as u128).wrapping_sub((b as u128) + ((borrow >> 63) as u128)); + (ret as u64, (ret >> 64) as u64) +} + +/// Compute a + (b * c) + carry, returning the result and the new carry over. +#[inline(always)] +const fn mac(a: u64, b: u64, c: u64, carry: u64) -> (u64, u64) { + let ret = (a as u128) + ((b as u128) * (c as u128)) + (carry as u128); + (ret as u64, (ret >> 64) as u64) +} + +pub use fp::*; +pub use fq::*; + +#[test] +fn test_extract() { + let a = Fq::random(); + let a = a.square(); + let (t, s) = a.extract_radix2_vartime().unwrap(); + assert_eq!( + t.pow_vartime(&[1 << Fq::S, 0, 0, 0]) * Fq::ROOT_OF_UNITY.pow_vartime(&[s, 0, 0, 0]), + a + ); + assert_eq!(a.deterministic_sqrt().unwrap().square(), a); +} diff --git a/src/arithmetic/fields/fp.rs b/src/arithmetic/fields/fp.rs new file mode 100644 index 0000000..f79fb07 --- /dev/null +++ b/src/arithmetic/fields/fp.rs @@ -0,0 +1,642 @@ +use super::{Field, Group}; + +use core::convert::TryInto; +use core::fmt; +use core::ops::{Add, AddAssign, Mul, MulAssign, Neg, Sub, SubAssign}; + +use subtle::{Choice, ConditionallySelectable, ConstantTimeEq, CtOption}; + +use super::{adc, mac, sbb}; + +/// This represents an element of $\mathbb{F}_p$ where +/// +/// `p = 0x40000000000000000000000000000000038aa1276c3f59b9a14064e200000001` +/// +/// is the base field of the Tweedledum curve. +// The internal representation of this type is four 64-bit unsigned +// integers in little-endian order. `Fp` values are always in +// Montgomery form; i.e., Fp(a) = aR mod p, with R = 2^256. +#[derive(Clone, Copy, Eq)] +pub struct Fp(pub(crate) [u64; 4]); + +impl fmt::Debug for Fp { + fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { + let tmp = self.to_bytes(); + write!(f, "0x")?; + for &b in tmp.iter().rev() { + write!(f, "{:02x}", b)?; + } + Ok(()) + } +} + +impl From for Fp { + fn from(bit: bool) -> Fp { + if bit { + Fp::one() + } else { + Fp::zero() + } + } +} + +impl From for Fp { + fn from(val: u64) -> Fp { + Fp([val, 0, 0, 0]) * R2 + } +} + +impl ConstantTimeEq for Fp { + fn ct_eq(&self, other: &Self) -> Choice { + self.0[0].ct_eq(&other.0[0]) + & self.0[1].ct_eq(&other.0[1]) + & self.0[2].ct_eq(&other.0[2]) + & self.0[3].ct_eq(&other.0[3]) + } +} + +impl PartialEq for Fp { + #[inline] + fn eq(&self, other: &Self) -> bool { + self.ct_eq(other).unwrap_u8() == 1 + } +} + +impl ConditionallySelectable for Fp { + fn conditional_select(a: &Self, b: &Self, choice: Choice) -> Self { + Fp([ + u64::conditional_select(&a.0[0], &b.0[0], choice), + u64::conditional_select(&a.0[1], &b.0[1], choice), + u64::conditional_select(&a.0[2], &b.0[2], choice), + u64::conditional_select(&a.0[3], &b.0[3], choice), + ]) + } +} + +/// Constant representing the modulus +/// p = 0x40000000000000000000000000000000038aa1276c3f59b9a14064e200000001 +const MODULUS: Fp = Fp([ + 0xa14064e200000001, + 0x38aa1276c3f59b9, + 0x0, + 0x4000000000000000, +]); + +impl<'a> Neg for &'a Fp { + type Output = Fp; + + #[inline] + fn neg(self) -> Fp { + self.neg() + } +} + +impl Neg for Fp { + type Output = Fp; + + #[inline] + fn neg(self) -> Fp { + -&self + } +} + +impl<'a, 'b> Sub<&'b Fp> for &'a Fp { + type Output = Fp; + + #[inline] + fn sub(self, rhs: &'b Fp) -> Fp { + self.sub(rhs) + } +} + +impl<'a, 'b> Add<&'b Fp> for &'a Fp { + type Output = Fp; + + #[inline] + fn add(self, rhs: &'b Fp) -> Fp { + self.add(rhs) + } +} + +impl<'a, 'b> Mul<&'b Fp> for &'a Fp { + type Output = Fp; + + #[inline] + fn mul(self, rhs: &'b Fp) -> Fp { + self.mul(rhs) + } +} + +impl_binops_additive!(Fp, Fp); +impl_binops_multiplicative!(Fp, Fp); + +/// INV = -(p^{-1} mod 2^64) mod 2^64 +const INV: u64 = 0xa14064e1ffffffff; + +/// R = 2^256 mod p +const R: Fp = Fp([ + 0x1c3ed159fffffffd, + 0xf5601c89bb41f2d3, + 0xffffffffffffffff, + 0x3fffffffffffffff, +]); + +/// R^2 = 2^512 mod p +const R2: Fp = Fp([ + 0x280c9c4000000010, + 0x91a4409b5400af74, + 0xdd7b28e19094c659, + 0xc8ad9107ccca0e, +]); + +/// R^3 = 2^768 mod p +const R3: Fp = Fp([ + 0x98fb3d144380a737, + 0xf9fdbeb55b7eb87c, + 0x63f75cb999eafa89, + 0x217cb214ebb8fc72, +]); + +const S: u32 = 33; + +/// GENERATOR^t where t * 2^s + 1 = p +/// with t odd. In other words, this +/// is a 2^s root of unity. +/// +/// `GENERATOR = 5 mod p` is a generator +/// of the p - 1 order multiplicative +/// subgroup. +const ROOT_OF_UNITY: Fp = Fp::from_raw([ + 0x53de9f31b88837ce, + 0xff46e8f3f3ea99d6, + 0xf624f2eaaf8c2d57, + 0x2ae45117890ee2fc, +]); + +impl Default for Fp { + #[inline] + fn default() -> Self { + Self::zero() + } +} + +impl Fp { + /// Returns zero, the additive identity. + #[inline] + pub const fn zero() -> Fp { + Fp([0, 0, 0, 0]) + } + + /// Returns one, the multiplicative identity. + #[inline] + pub const fn one() -> Fp { + R + } + + /// Doubles this field element. + #[inline] + pub const fn double(&self) -> Fp { + // TODO: This can be achieved more efficiently with a bitshift. + self.add(self) + } + + fn from_u512(limbs: [u64; 8]) -> Fp { + // We reduce an arbitrary 512-bit number by decomposing it into two 256-bit digits + // with the higher bits multiplied by 2^256. Thus, we perform two reductions + // + // 1. the lower bits are multiplied by R^2, as normal + // 2. the upper bits are multiplied by R^2 * 2^256 = R^3 + // + // and computing their sum in the field. It remains to see that arbitrary 256-bit + // numbers can be placed into Montgomery form safely using the reduction. The + // reduction works so long as the product is less than R=2^256 multipled by + // the modulus. This holds because for any `c` smaller than the modulus, we have + // that (2^256 - 1)*c is an acceptable product for the reduction. Therefore, the + // reduction always works so long as `c` is in the field; in this case it is either the + // constant `R2` or `R3`. + let d0 = Fp([limbs[0], limbs[1], limbs[2], limbs[3]]); + let d1 = Fp([limbs[4], limbs[5], limbs[6], limbs[7]]); + // Convert to Montgomery form + d0 * R2 + d1 * R3 + } + + /// Converts from an integer represented in little endian + /// into its (congruent) `Fp` representation. + pub const fn from_raw(val: [u64; 4]) -> Self { + (&Fp(val)).mul(&R2) + } + + /// Squares this element. + #[inline] + pub const fn square(&self) -> Fp { + let (r1, carry) = mac(0, self.0[0], self.0[1], 0); + let (r2, carry) = mac(0, self.0[0], self.0[2], carry); + let (r3, r4) = mac(0, self.0[0], self.0[3], carry); + + let (r3, carry) = mac(r3, self.0[1], self.0[2], 0); + let (r4, r5) = mac(r4, self.0[1], self.0[3], carry); + + let (r5, r6) = mac(r5, self.0[2], self.0[3], 0); + + let r7 = r6 >> 63; + let r6 = (r6 << 1) | (r5 >> 63); + let r5 = (r5 << 1) | (r4 >> 63); + let r4 = (r4 << 1) | (r3 >> 63); + let r3 = (r3 << 1) | (r2 >> 63); + let r2 = (r2 << 1) | (r1 >> 63); + let r1 = r1 << 1; + + let (r0, carry) = mac(0, self.0[0], self.0[0], 0); + let (r1, carry) = adc(0, r1, carry); + let (r2, carry) = mac(r2, self.0[1], self.0[1], carry); + let (r3, carry) = adc(0, r3, carry); + let (r4, carry) = mac(r4, self.0[2], self.0[2], carry); + let (r5, carry) = adc(0, r5, carry); + let (r6, carry) = mac(r6, self.0[3], self.0[3], carry); + let (r7, _) = adc(0, r7, carry); + + Fp::montgomery_reduce(r0, r1, r2, r3, r4, r5, r6, r7) + } + + #[inline(always)] + const fn montgomery_reduce( + r0: u64, + r1: u64, + r2: u64, + r3: u64, + r4: u64, + r5: u64, + r6: u64, + r7: u64, + ) -> Self { + // The Montgomery reduction here is based on Algorithm 14.32 in + // Handbook of Applied Cryptography + // . + + let k = r0.wrapping_mul(INV); + let (_, carry) = mac(r0, k, MODULUS.0[0], 0); + let (r1, carry) = mac(r1, k, MODULUS.0[1], carry); + let (r2, carry) = mac(r2, k, MODULUS.0[2], carry); + let (r3, carry) = mac(r3, k, MODULUS.0[3], carry); + let (r4, carry2) = adc(r4, 0, carry); + + let k = r1.wrapping_mul(INV); + let (_, carry) = mac(r1, k, MODULUS.0[0], 0); + let (r2, carry) = mac(r2, k, MODULUS.0[1], carry); + let (r3, carry) = mac(r3, k, MODULUS.0[2], carry); + let (r4, carry) = mac(r4, k, MODULUS.0[3], carry); + let (r5, carry2) = adc(r5, carry2, carry); + + let k = r2.wrapping_mul(INV); + let (_, carry) = mac(r2, k, MODULUS.0[0], 0); + let (r3, carry) = mac(r3, k, MODULUS.0[1], carry); + let (r4, carry) = mac(r4, k, MODULUS.0[2], carry); + let (r5, carry) = mac(r5, k, MODULUS.0[3], carry); + let (r6, carry2) = adc(r6, carry2, carry); + + let k = r3.wrapping_mul(INV); + let (_, carry) = mac(r3, k, MODULUS.0[0], 0); + let (r4, carry) = mac(r4, k, MODULUS.0[1], carry); + let (r5, carry) = mac(r5, k, MODULUS.0[2], carry); + let (r6, carry) = mac(r6, k, MODULUS.0[3], carry); + let (r7, _) = adc(r7, carry2, carry); + + // Result may be within MODULUS of the correct value + (&Fp([r4, r5, r6, r7])).sub(&MODULUS) + } + + /// Multiplies `rhs` by `self`, returning the result. + #[inline] + pub const fn mul(&self, rhs: &Self) -> Self { + // Schoolbook multiplication + + let (r0, carry) = mac(0, self.0[0], rhs.0[0], 0); + let (r1, carry) = mac(0, self.0[0], rhs.0[1], carry); + let (r2, carry) = mac(0, self.0[0], rhs.0[2], carry); + let (r3, r4) = mac(0, self.0[0], rhs.0[3], carry); + + let (r1, carry) = mac(r1, self.0[1], rhs.0[0], 0); + let (r2, carry) = mac(r2, self.0[1], rhs.0[1], carry); + let (r3, carry) = mac(r3, self.0[1], rhs.0[2], carry); + let (r4, r5) = mac(r4, self.0[1], rhs.0[3], carry); + + let (r2, carry) = mac(r2, self.0[2], rhs.0[0], 0); + let (r3, carry) = mac(r3, self.0[2], rhs.0[1], carry); + let (r4, carry) = mac(r4, self.0[2], rhs.0[2], carry); + let (r5, r6) = mac(r5, self.0[2], rhs.0[3], carry); + + let (r3, carry) = mac(r3, self.0[3], rhs.0[0], 0); + let (r4, carry) = mac(r4, self.0[3], rhs.0[1], carry); + let (r5, carry) = mac(r5, self.0[3], rhs.0[2], carry); + let (r6, r7) = mac(r6, self.0[3], rhs.0[3], carry); + + Fp::montgomery_reduce(r0, r1, r2, r3, r4, r5, r6, r7) + } + + /// Subtracts `rhs` from `self`, returning the result. + #[inline] + pub const fn sub(&self, rhs: &Self) -> Self { + let (d0, borrow) = sbb(self.0[0], rhs.0[0], 0); + let (d1, borrow) = sbb(self.0[1], rhs.0[1], borrow); + let (d2, borrow) = sbb(self.0[2], rhs.0[2], borrow); + let (d3, borrow) = sbb(self.0[3], rhs.0[3], borrow); + + // If underflow occurred on the final limb, borrow = 0xfff...fff, otherwise + // borrow = 0x000...000. Thus, we use it as a mask to conditionally add the modulus. + let (d0, carry) = adc(d0, MODULUS.0[0] & borrow, 0); + let (d1, carry) = adc(d1, MODULUS.0[1] & borrow, carry); + let (d2, carry) = adc(d2, MODULUS.0[2] & borrow, carry); + let (d3, _) = adc(d3, MODULUS.0[3] & borrow, carry); + + Fp([d0, d1, d2, d3]) + } + + /// Adds `rhs` to `self`, returning the result. + #[inline] + pub const fn add(&self, rhs: &Self) -> Self { + let (d0, carry) = adc(self.0[0], rhs.0[0], 0); + let (d1, carry) = adc(self.0[1], rhs.0[1], carry); + let (d2, carry) = adc(self.0[2], rhs.0[2], carry); + let (d3, _) = adc(self.0[3], rhs.0[3], carry); + + // Attempt to subtract the modulus, to ensure the value + // is smaller than the modulus. + (&Fp([d0, d1, d2, d3])).sub(&MODULUS) + } + + /// Negates `self`. + #[inline] + pub const fn neg(&self) -> Self { + // Subtract `self` from `MODULUS` to negate. Ignore the final + // borrow because it cannot underflow; self is guaranteed to + // be in the field. + let (d0, borrow) = sbb(MODULUS.0[0], self.0[0], 0); + let (d1, borrow) = sbb(MODULUS.0[1], self.0[1], borrow); + let (d2, borrow) = sbb(MODULUS.0[2], self.0[2], borrow); + let (d3, _) = sbb(MODULUS.0[3], self.0[3], borrow); + + // `tmp` could be `MODULUS` if `self` was zero. Create a mask that is + // zero if `self` was zero, and `u64::max_value()` if self was nonzero. + let mask = (((self.0[0] | self.0[1] | self.0[2] | self.0[3]) == 0) as u64).wrapping_sub(1); + + Fp([d0 & mask, d1 & mask, d2 & mask, d3 & mask]) + } +} + +impl<'a> From<&'a Fp> for [u8; 32] { + fn from(value: &'a Fp) -> [u8; 32] { + value.to_bytes() + } +} + +impl Group for Fp { + type Scalar = Fp; + + fn group_zero() -> Self { + Self::zero() + } + fn group_add(&mut self, rhs: &Self) { + *self = *self + *rhs; + } + fn group_sub(&mut self, rhs: &Self) { + *self = *self - *rhs; + } + fn group_scale(&mut self, by: &Self::Scalar) { + *self = *self * (*by); + } +} + +impl Field for Fp { + const NUM_BITS: u32 = 255; + const CAPACITY: u32 = 254; + const S: u32 = S; + const ROOT_OF_UNITY: Self = ROOT_OF_UNITY; + const ROOT_OF_UNITY_INV: Self = Fp::from_raw([ + 0x9246674078fa45bb, + 0xd822ebd60888c5ea, + 0x56d579133a11731f, + 0x1c88fa9e942120bb, + ]); + const UNROLL_T_EXPONENT: [u64; 4] = [ + 0x3b3a6633d1897d83, + 0x0000000000c93d5b, + 0xf000000000000000, + 0xe34ab16, + ]; + const T_EXPONENT: [u64; 4] = [ + 0xb61facdcd0a03271, + 0x0000000001c55093, + 0x0000000000000000, + 0x20000000, + ]; + const UNROLL_S_EXPONENT: u64 = 0x11cb54e91; + const TWO_INV: Self = Fp::from_raw([ + 0xd0a0327100000001, + 0x01c55093b61facdc, + 0x0000000000000000, + 0x2000000000000000, + ]); + const RESCUE_ALPHA: u64 = 5; + const RESCUE_INVALPHA: [u64; 4] = [ + 0x810050b4cccccccd, + 0x360880ec56991494, + 0x3333333333333333, + 0x3333333333333333, + ]; + const ZETA: Self = Fp::from_raw([ + 0x8598abb3a410c9c8, + 0x7881fb239ba41a26, + 0x9bebc9146ef83d9a, + 0x1508415ab5e97c94, + ]); + + fn is_zero(&self) -> Choice { + self.ct_eq(&Self::zero()) + } + + fn zero() -> Self { + Self::zero() + } + + fn one() -> Self { + Self::one() + } + + fn from_u64(v: u64) -> Self { + Fp::from_raw([v as u64, 0, 0, 0]) + } + + fn from_u128(v: u128) -> Self { + Fp::from_raw([v as u64, (v >> 64) as u64, 0, 0]) + } + + fn double(&self) -> Self { + self.double() + } + + #[inline(always)] + fn square(&self) -> Self { + self.square() + } + + /// Computes the square root of this element, if it exists. + fn sqrt(&self) -> CtOption { + // Tonelli-Shank's algorithm for p mod 16 = 1 + // https://eprint.iacr.org/2012/685.pdf (page 12, algorithm 5) + + // w = self^((t - 1) // 2) + let w = self.pow_vartime(&[0xdb0fd66e68501938, 0xe2a849, 0x0, 0x10000000]); + + let mut v = S; + let mut x = self * w; + let mut b = x * w; + + // Initialize z as the 2^S root of unity. + let mut z = ROOT_OF_UNITY; + + for max_v in (1..=S).rev() { + let mut k = 1; + let mut tmp = b.square(); + let mut j_less_than_v: Choice = 1.into(); + + for j in 2..max_v { + let tmp_is_one = tmp.ct_eq(&Fp::one()); + let squared = Fp::conditional_select(&tmp, &z, tmp_is_one).square(); + tmp = Fp::conditional_select(&squared, &tmp, tmp_is_one); + let new_z = Fp::conditional_select(&z, &squared, tmp_is_one); + j_less_than_v &= !j.ct_eq(&v); + k = u32::conditional_select(&j, &k, tmp_is_one); + z = Fp::conditional_select(&z, &new_z, j_less_than_v); + } + + let result = x * z; + x = Fp::conditional_select(&result, &x, b.ct_eq(&Fp::one())); + z = z.square(); + b *= z; + v = k; + } + + CtOption::new( + x, + (x * x).ct_eq(self), // Only return Some if it's the square root. + ) + } + + /// Computes the multiplicative inverse of this element, + /// failing if the element is zero. + fn invert(&self) -> CtOption { + let tmp = self.pow_vartime(&[ + 0xa14064e1ffffffff, + 0x38aa1276c3f59b9, + 0x0, + 0x4000000000000000, + ]); + + CtOption::new(tmp, !self.ct_eq(&Self::zero())) + } + + /// Attempts to convert a little-endian byte representation of + /// a scalar into a `Fp`, failing if the input is not canonical. + fn from_bytes(bytes: &[u8; 32]) -> CtOption { + let mut tmp = Fp([0, 0, 0, 0]); + + tmp.0[0] = u64::from_le_bytes(bytes[0..8].try_into().unwrap()); + tmp.0[1] = u64::from_le_bytes(bytes[8..16].try_into().unwrap()); + tmp.0[2] = u64::from_le_bytes(bytes[16..24].try_into().unwrap()); + tmp.0[3] = u64::from_le_bytes(bytes[24..32].try_into().unwrap()); + + // Try to subtract the modulus + let (_, borrow) = sbb(tmp.0[0], MODULUS.0[0], 0); + let (_, borrow) = sbb(tmp.0[1], MODULUS.0[1], borrow); + let (_, borrow) = sbb(tmp.0[2], MODULUS.0[2], borrow); + let (_, borrow) = sbb(tmp.0[3], MODULUS.0[3], borrow); + + // If the element is smaller than MODULUS then the + // subtraction will underflow, producing a borrow value + // of 0xffff...ffff. Otherwise, it'll be zero. + let is_some = (borrow as u8) & 1; + + // Convert to Montgomery form by computing + // (a.R^0 * R^2) / R = a.R + tmp *= &R2; + + CtOption::new(tmp, Choice::from(is_some)) + } + + /// Converts an element of `Fp` into a byte representation in + /// little-endian byte order. + fn to_bytes(&self) -> [u8; 32] { + // Turn into canonical form by computing + // (a.R) / R = a + let tmp = Fp::montgomery_reduce(self.0[0], self.0[1], self.0[2], self.0[3], 0, 0, 0, 0); + + let mut res = [0; 32]; + res[0..8].copy_from_slice(&tmp.0[0].to_le_bytes()); + res[8..16].copy_from_slice(&tmp.0[1].to_le_bytes()); + res[16..24].copy_from_slice(&tmp.0[2].to_le_bytes()); + res[24..32].copy_from_slice(&tmp.0[3].to_le_bytes()); + + res + } + + /// Converts a 512-bit little endian integer into + /// a `Fp` by reducing by the modulus. + fn from_bytes_wide(bytes: &[u8; 64]) -> Fp { + Fp::from_u512([ + u64::from_le_bytes(bytes[0..8].try_into().unwrap()), + u64::from_le_bytes(bytes[8..16].try_into().unwrap()), + u64::from_le_bytes(bytes[16..24].try_into().unwrap()), + u64::from_le_bytes(bytes[24..32].try_into().unwrap()), + u64::from_le_bytes(bytes[32..40].try_into().unwrap()), + u64::from_le_bytes(bytes[40..48].try_into().unwrap()), + u64::from_le_bytes(bytes[48..56].try_into().unwrap()), + u64::from_le_bytes(bytes[56..64].try_into().unwrap()), + ]) + } + + fn get_lower_128(&self) -> u128 { + let tmp = Fp::montgomery_reduce(self.0[0], self.0[1], self.0[2], self.0[3], 0, 0, 0, 0); + + u128::from(tmp.0[0]) | (u128::from(tmp.0[1]) << 64) + } +} + +#[test] +fn test_inv() { + // Compute -(r^{-1} mod 2^64) mod 2^64 by exponentiating + // by totient(2**64) - 1 + + let mut inv = 1u64; + for _ in 0..63 { + inv = inv.wrapping_mul(inv); + inv = inv.wrapping_mul(MODULUS.0[0]); + } + inv = inv.wrapping_neg(); + + assert_eq!(inv, INV); +} + +#[test] +fn test_zeta() { + assert_eq!( + format!("{:?}", Fp::ZETA), + "0x1508415ab5e97c949bebc9146ef83d9a7881fb239ba41a268598abb3a410c9c8" + ); + + let a = Fp::ZETA; + assert!(bool::from(a != Fp::one())); + let b = a * a; + assert!(bool::from(b != Fp::one())); + let c = b * a; + assert!(bool::from(c == Fp::one())); +} + +#[test] +fn test_inv_root_of_unity() { + assert_eq!(Fp::ROOT_OF_UNITY_INV, Fp::ROOT_OF_UNITY.invert().unwrap()); +} + +#[test] +fn test_inv_2() { + assert_eq!(Fp::TWO_INV, Fp::from(2).invert().unwrap()); +} diff --git a/src/arithmetic/fields/fq.rs b/src/arithmetic/fields/fq.rs new file mode 100644 index 0000000..24358e3 --- /dev/null +++ b/src/arithmetic/fields/fq.rs @@ -0,0 +1,656 @@ +use super::{Field, Group}; + +use core::convert::TryInto; +use core::fmt; +use core::ops::{Add, AddAssign, Mul, MulAssign, Neg, Sub, SubAssign}; + +use subtle::{Choice, ConditionallySelectable, ConstantTimeEq, CtOption}; + +use super::{adc, mac, sbb}; + +/// This represents an element of $\mathbb{F}_q$ where +/// +/// `q = 0x40000000000000000000000000000000038aa127696286c9842cafd400000001` +/// +/// is the base field of the Tweedledee curve. +// The internal representation of this type is four 64-bit unsigned +// integers in little-endian order. `Fq` values are always in +// Montgomery form; i.e., Fq(a) = aR mod q, with R = 2^256. +#[derive(Clone, Copy, Eq)] +pub struct Fq(pub(crate) [u64; 4]); + +impl fmt::Debug for Fq { + fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { + let tmp = self.to_bytes(); + write!(f, "0x")?; + for &b in tmp.iter().rev() { + write!(f, "{:02x}", b)?; + } + Ok(()) + } +} + +impl From for Fq { + fn from(bit: bool) -> Fq { + if bit { + Fq::one() + } else { + Fq::zero() + } + } +} + +impl From for Fq { + fn from(val: u64) -> Fq { + Fq([val, 0, 0, 0]) * R2 + } +} + +impl ConstantTimeEq for Fq { + fn ct_eq(&self, other: &Self) -> Choice { + self.0[0].ct_eq(&other.0[0]) + & self.0[1].ct_eq(&other.0[1]) + & self.0[2].ct_eq(&other.0[2]) + & self.0[3].ct_eq(&other.0[3]) + } +} + +impl PartialEq for Fq { + #[inline] + fn eq(&self, other: &Self) -> bool { + self.ct_eq(other).unwrap_u8() == 1 + } +} + +impl ConditionallySelectable for Fq { + fn conditional_select(a: &Self, b: &Self, choice: Choice) -> Self { + Fq([ + u64::conditional_select(&a.0[0], &b.0[0], choice), + u64::conditional_select(&a.0[1], &b.0[1], choice), + u64::conditional_select(&a.0[2], &b.0[2], choice), + u64::conditional_select(&a.0[3], &b.0[3], choice), + ]) + } +} + +/// Constant representing the modulus +/// q = 0x40000000000000000000000000000000038aa127696286c9842cafd400000001 +const MODULUS: Fq = Fq([ + 0x842cafd400000001, + 0x38aa127696286c9, + 0x0, + 0x4000000000000000, +]); + +impl<'a> Neg for &'a Fq { + type Output = Fq; + + #[inline] + fn neg(self) -> Fq { + self.neg() + } +} + +impl Neg for Fq { + type Output = Fq; + + #[inline] + fn neg(self) -> Fq { + -&self + } +} + +impl<'a, 'b> Sub<&'b Fq> for &'a Fq { + type Output = Fq; + + #[inline] + fn sub(self, rhs: &'b Fq) -> Fq { + self.sub(rhs) + } +} + +impl<'a, 'b> Add<&'b Fq> for &'a Fq { + type Output = Fq; + + #[inline] + fn add(self, rhs: &'b Fq) -> Fq { + self.add(rhs) + } +} + +impl<'a, 'b> Mul<&'b Fq> for &'a Fq { + type Output = Fq; + + #[inline] + fn mul(self, rhs: &'b Fq) -> Fq { + self.mul(rhs) + } +} + +impl_binops_additive!(Fq, Fq); +impl_binops_multiplicative!(Fq, Fq); + +/// INV = -(q^{-1} mod 2^64) mod 2^64 +const INV: u64 = 0x842cafd3ffffffff; + +/// R = 2^256 mod q +const R: Fq = Fq([ + 0x7379f083fffffffd, + 0xf5601c89c3d86ba3, + 0xffffffffffffffff, + 0x3fffffffffffffff, +]); + +/// R^2 = 2^512 mod q +const R2: Fq = Fq([ + 0x8595fa8000000010, + 0x7e16a565c6895230, + 0xf4c0e6fcb03aa0a2, + 0xc8ad9106886013, +]); + +/// R^3 = 2^768 mod q +const R3: Fq = Fq([ + 0xa624f338075cdb5e, + 0x57964eacb8fe21f2, + 0xcb266d18c0413bc2, + 0xa42cdf95f959577, +]); + +const S: u32 = 34; + +/// GENERATOR^t where t * 2^s + 1 = q +/// with t odd. In other words, this +/// is a 2^s root of unity. +/// +/// `GENERATOR = 5 mod q` is a generator +/// of the q - 1 order multiplicative +/// subgroup. +const ROOT_OF_UNITY: Fq = Fq::from_raw([ + 0x1cbd3234869d57ec, + 0xa287dd1b8084fbf, + 0xf1dbcb645a987293, + 0x113efc510dc03c0b, +]); + +impl Default for Fq { + #[inline] + fn default() -> Self { + Self::zero() + } +} + +impl Fq { + /// Returns zero, the additive identity. + #[inline] + pub const fn zero() -> Fq { + Fq([0, 0, 0, 0]) + } + + /// Returns one, the multiplicative identity. + #[inline] + pub const fn one() -> Fq { + R + } + + /// Doubles this field element. + #[inline] + pub const fn double(&self) -> Fq { + // TODO: This can be achieved more efficiently with a bitshift. + self.add(self) + } + + /// Converts a 512-bit little endian integer into + /// a `Fq` by reducing by the modulus. + pub fn from_bytes_wide(bytes: &[u8; 64]) -> Fq { + Fq::from_u512([ + u64::from_le_bytes(bytes[0..8].try_into().unwrap()), + u64::from_le_bytes(bytes[8..16].try_into().unwrap()), + u64::from_le_bytes(bytes[16..24].try_into().unwrap()), + u64::from_le_bytes(bytes[24..32].try_into().unwrap()), + u64::from_le_bytes(bytes[32..40].try_into().unwrap()), + u64::from_le_bytes(bytes[40..48].try_into().unwrap()), + u64::from_le_bytes(bytes[48..56].try_into().unwrap()), + u64::from_le_bytes(bytes[56..64].try_into().unwrap()), + ]) + } + + fn from_u512(limbs: [u64; 8]) -> Fq { + // We reduce an arbitrary 512-bit number by decomposing it into two 256-bit digits + // with the higher bits multiplied by 2^256. Thus, we perform two reductions + // + // 1. the lower bits are multiplied by R^2, as normal + // 2. the upper bits are multiplied by R^2 * 2^256 = R^3 + // + // and computing their sum in the field. It remains to see that arbitrary 256-bit + // numbers can be placed into Montgomery form safely using the reduction. The + // reduction works so long as the product is less than R=2^256 multipled by + // the modulus. This holds because for any `c` smaller than the modulus, we have + // that (2^256 - 1)*c is an acceptable product for the reduction. Therefore, the + // reduction always works so long as `c` is in the field; in this case it is either the + // constant `R2` or `R3`. + let d0 = Fq([limbs[0], limbs[1], limbs[2], limbs[3]]); + let d1 = Fq([limbs[4], limbs[5], limbs[6], limbs[7]]); + // Convert to Montgomery form + d0 * R2 + d1 * R3 + } + + /// Converts from an integer represented in little endian + /// into its (congruent) `Fq` representation. + pub const fn from_raw(val: [u64; 4]) -> Self { + (&Fq(val)).mul(&R2) + } + + /// Squares this element. + #[inline] + pub const fn square(&self) -> Fq { + let (r1, carry) = mac(0, self.0[0], self.0[1], 0); + let (r2, carry) = mac(0, self.0[0], self.0[2], carry); + let (r3, r4) = mac(0, self.0[0], self.0[3], carry); + + let (r3, carry) = mac(r3, self.0[1], self.0[2], 0); + let (r4, r5) = mac(r4, self.0[1], self.0[3], carry); + + let (r5, r6) = mac(r5, self.0[2], self.0[3], 0); + + let r7 = r6 >> 63; + let r6 = (r6 << 1) | (r5 >> 63); + let r5 = (r5 << 1) | (r4 >> 63); + let r4 = (r4 << 1) | (r3 >> 63); + let r3 = (r3 << 1) | (r2 >> 63); + let r2 = (r2 << 1) | (r1 >> 63); + let r1 = r1 << 1; + + let (r0, carry) = mac(0, self.0[0], self.0[0], 0); + let (r1, carry) = adc(0, r1, carry); + let (r2, carry) = mac(r2, self.0[1], self.0[1], carry); + let (r3, carry) = adc(0, r3, carry); + let (r4, carry) = mac(r4, self.0[2], self.0[2], carry); + let (r5, carry) = adc(0, r5, carry); + let (r6, carry) = mac(r6, self.0[3], self.0[3], carry); + let (r7, _) = adc(0, r7, carry); + + Fq::montgomery_reduce(r0, r1, r2, r3, r4, r5, r6, r7) + } + + #[inline(always)] + const fn montgomery_reduce( + r0: u64, + r1: u64, + r2: u64, + r3: u64, + r4: u64, + r5: u64, + r6: u64, + r7: u64, + ) -> Self { + // The Montgomery reduction here is based on Algorithm 14.32 in + // Handbook of Applied Cryptography + // . + + let k = r0.wrapping_mul(INV); + let (_, carry) = mac(r0, k, MODULUS.0[0], 0); + let (r1, carry) = mac(r1, k, MODULUS.0[1], carry); + let (r2, carry) = mac(r2, k, MODULUS.0[2], carry); + let (r3, carry) = mac(r3, k, MODULUS.0[3], carry); + let (r4, carry2) = adc(r4, 0, carry); + + let k = r1.wrapping_mul(INV); + let (_, carry) = mac(r1, k, MODULUS.0[0], 0); + let (r2, carry) = mac(r2, k, MODULUS.0[1], carry); + let (r3, carry) = mac(r3, k, MODULUS.0[2], carry); + let (r4, carry) = mac(r4, k, MODULUS.0[3], carry); + let (r5, carry2) = adc(r5, carry2, carry); + + let k = r2.wrapping_mul(INV); + let (_, carry) = mac(r2, k, MODULUS.0[0], 0); + let (r3, carry) = mac(r3, k, MODULUS.0[1], carry); + let (r4, carry) = mac(r4, k, MODULUS.0[2], carry); + let (r5, carry) = mac(r5, k, MODULUS.0[3], carry); + let (r6, carry2) = adc(r6, carry2, carry); + + let k = r3.wrapping_mul(INV); + let (_, carry) = mac(r3, k, MODULUS.0[0], 0); + let (r4, carry) = mac(r4, k, MODULUS.0[1], carry); + let (r5, carry) = mac(r5, k, MODULUS.0[2], carry); + let (r6, carry) = mac(r6, k, MODULUS.0[3], carry); + let (r7, _) = adc(r7, carry2, carry); + + // Result may be within MODULUS of the correct value + (&Fq([r4, r5, r6, r7])).sub(&MODULUS) + } + + /// Multiplies `rhs` by `self`, returning the result. + #[inline] + pub const fn mul(&self, rhs: &Self) -> Self { + // Schoolbook multiplication + + let (r0, carry) = mac(0, self.0[0], rhs.0[0], 0); + let (r1, carry) = mac(0, self.0[0], rhs.0[1], carry); + let (r2, carry) = mac(0, self.0[0], rhs.0[2], carry); + let (r3, r4) = mac(0, self.0[0], rhs.0[3], carry); + + let (r1, carry) = mac(r1, self.0[1], rhs.0[0], 0); + let (r2, carry) = mac(r2, self.0[1], rhs.0[1], carry); + let (r3, carry) = mac(r3, self.0[1], rhs.0[2], carry); + let (r4, r5) = mac(r4, self.0[1], rhs.0[3], carry); + + let (r2, carry) = mac(r2, self.0[2], rhs.0[0], 0); + let (r3, carry) = mac(r3, self.0[2], rhs.0[1], carry); + let (r4, carry) = mac(r4, self.0[2], rhs.0[2], carry); + let (r5, r6) = mac(r5, self.0[2], rhs.0[3], carry); + + let (r3, carry) = mac(r3, self.0[3], rhs.0[0], 0); + let (r4, carry) = mac(r4, self.0[3], rhs.0[1], carry); + let (r5, carry) = mac(r5, self.0[3], rhs.0[2], carry); + let (r6, r7) = mac(r6, self.0[3], rhs.0[3], carry); + + Fq::montgomery_reduce(r0, r1, r2, r3, r4, r5, r6, r7) + } + + /// Subtracts `rhs` from `self`, returning the result. + #[inline] + pub const fn sub(&self, rhs: &Self) -> Self { + let (d0, borrow) = sbb(self.0[0], rhs.0[0], 0); + let (d1, borrow) = sbb(self.0[1], rhs.0[1], borrow); + let (d2, borrow) = sbb(self.0[2], rhs.0[2], borrow); + let (d3, borrow) = sbb(self.0[3], rhs.0[3], borrow); + + // If underflow occurred on the final limb, borrow = 0xfff...fff, otherwise + // borrow = 0x000...000. Thus, we use it as a mask to conditionally add the modulus. + let (d0, carry) = adc(d0, MODULUS.0[0] & borrow, 0); + let (d1, carry) = adc(d1, MODULUS.0[1] & borrow, carry); + let (d2, carry) = adc(d2, MODULUS.0[2] & borrow, carry); + let (d3, _) = adc(d3, MODULUS.0[3] & borrow, carry); + + Fq([d0, d1, d2, d3]) + } + + /// Adds `rhs` to `self`, returning the result. + #[inline] + pub const fn add(&self, rhs: &Self) -> Self { + let (d0, carry) = adc(self.0[0], rhs.0[0], 0); + let (d1, carry) = adc(self.0[1], rhs.0[1], carry); + let (d2, carry) = adc(self.0[2], rhs.0[2], carry); + let (d3, _) = adc(self.0[3], rhs.0[3], carry); + + // Attempt to subtract the modulus, to ensure the value + // is smaller than the modulus. + (&Fq([d0, d1, d2, d3])).sub(&MODULUS) + } + + /// Negates `self`. + #[inline] + pub const fn neg(&self) -> Self { + // Subtract `self` from `MODULUS` to negate. Ignore the final + // borrow because it cannot underflow; self is guaranteed to + // be in the field. + let (d0, borrow) = sbb(MODULUS.0[0], self.0[0], 0); + let (d1, borrow) = sbb(MODULUS.0[1], self.0[1], borrow); + let (d2, borrow) = sbb(MODULUS.0[2], self.0[2], borrow); + let (d3, _) = sbb(MODULUS.0[3], self.0[3], borrow); + + // `tmp` could be `MODULUS` if `self` was zero. Create a mask that is + // zero if `self` was zero, and `u64::max_value()` if self was nonzero. + let mask = (((self.0[0] | self.0[1] | self.0[2] | self.0[3]) == 0) as u64).wrapping_sub(1); + + Fq([d0 & mask, d1 & mask, d2 & mask, d3 & mask]) + } +} + +impl<'a> From<&'a Fq> for [u8; 32] { + fn from(value: &'a Fq) -> [u8; 32] { + value.to_bytes() + } +} + +impl Group for Fq { + type Scalar = Fq; + + fn group_zero() -> Self { + Self::zero() + } + fn group_add(&mut self, rhs: &Self) { + *self = *self + *rhs; + } + fn group_sub(&mut self, rhs: &Self) { + *self = *self - *rhs; + } + fn group_scale(&mut self, by: &Self::Scalar) { + *self = *self * (*by); + } +} + +impl Field for Fq { + const NUM_BITS: u32 = 255; + const CAPACITY: u32 = 254; + const S: u32 = S; + const ROOT_OF_UNITY: Self = ROOT_OF_UNITY; + const ROOT_OF_UNITY_INV: Self = Fq::from_raw([ + 0x824860d3eb30de02, + 0xad9f0afd0ea63acc, + 0xd250318c11a16fe1, + 0x12a9f5e1dd62dabc, + ]); + const UNROLL_T_EXPONENT: [u64; 4] = [ + 0x9b71de17e6d2d5a0, + 0x0000000000296ee0, + 0x8c00000000000000, + 0x2ecc05e, + ]; + const T_EXPONENT: [u64; 4] = [ + 0xda58a1b2610b2bf5, + 0x0000000000e2a849, + 0x0000000000000000, + 0x10000000, + ]; + const UNROLL_S_EXPONENT: u64 = 0x344cfe85d; + const TWO_INV: Self = Fq::from_raw([ + 0xc21657ea00000001, + 0x01c55093b4b14364, + 0x0000000000000000, + 0x2000000000000000, + ]); + const RESCUE_ALPHA: u64 = 5; + const RESCUE_INVALPHA: [u64; 4] = [ + 0xd023bfdccccccccd, + 0x360880ec544ed23a, + 0x3333333333333333, + 0x3333333333333333, + ]; + const ZETA: Self = Fq::from_raw([ + 0x4394c2bd148fa4fd, + 0x69cf8de720e52ec1, + 0x87ad8b5ff9731ffe, + 0x36c66d3a1e049a58, + ]); + + fn is_zero(&self) -> Choice { + self.ct_eq(&Self::zero()) + } + + fn zero() -> Self { + Self::zero() + } + + fn one() -> Self { + Self::one() + } + + fn from_u64(v: u64) -> Self { + Fq::from_raw([v as u64, 0, 0, 0]) + } + + fn from_u128(v: u128) -> Self { + Fq::from_raw([v as u64, (v >> 64) as u64, 0, 0]) + } + + fn double(&self) -> Self { + self.double() + } + + #[inline(always)] + fn square(&self) -> Self { + self.square() + } + + /// Computes the square root of this element, if it exists. + fn sqrt(&self) -> CtOption { + // Tonelli-Shank's algorithm for q mod 16 = 1 + // https://eprint.iacr.org/2012/685.pdf (page 12, algorithm 5) + + // w = self^((t - 1) // 2) + let w = self.pow_vartime(&[0xed2c50d9308595fa, 0x715424, 0x0, 0x8000000]); + + let mut v = S; + let mut x = self * w; + let mut b = x * w; + + // Initialize z as the 2^S root of unity. + let mut z = ROOT_OF_UNITY; + + for max_v in (1..=S).rev() { + let mut k = 1; + let mut tmp = b.square(); + let mut j_less_than_v: Choice = 1.into(); + + for j in 2..max_v { + let tmp_is_one = tmp.ct_eq(&Fq::one()); + let squared = Fq::conditional_select(&tmp, &z, tmp_is_one).square(); + tmp = Fq::conditional_select(&squared, &tmp, tmp_is_one); + let new_z = Fq::conditional_select(&z, &squared, tmp_is_one); + j_less_than_v &= !j.ct_eq(&v); + k = u32::conditional_select(&j, &k, tmp_is_one); + z = Fq::conditional_select(&z, &new_z, j_less_than_v); + } + + let result = x * z; + x = Fq::conditional_select(&result, &x, b.ct_eq(&Fq::one())); + z = z.square(); + b *= z; + v = k; + } + + CtOption::new( + x, + (x * x).ct_eq(self), // Only return Some if it's the square root. + ) + } + + /// Computes the multiplicative inverse of this element, + /// failing if the element is zero. + fn invert(&self) -> CtOption { + let tmp = self.pow_vartime(&[ + 0x842cafd3ffffffff, + 0x38aa127696286c9, + 0x0, + 0x4000000000000000, + ]); + + CtOption::new(tmp, !self.ct_eq(&Self::zero())) + } + + /// Attempts to convert a little-endian byte representation of + /// a scalar into a `Fq`, failing if the input is not canonical. + fn from_bytes(bytes: &[u8; 32]) -> CtOption { + let mut tmp = Fq([0, 0, 0, 0]); + + tmp.0[0] = u64::from_le_bytes(bytes[0..8].try_into().unwrap()); + tmp.0[1] = u64::from_le_bytes(bytes[8..16].try_into().unwrap()); + tmp.0[2] = u64::from_le_bytes(bytes[16..24].try_into().unwrap()); + tmp.0[3] = u64::from_le_bytes(bytes[24..32].try_into().unwrap()); + + // Try to subtract the modulus + let (_, borrow) = sbb(tmp.0[0], MODULUS.0[0], 0); + let (_, borrow) = sbb(tmp.0[1], MODULUS.0[1], borrow); + let (_, borrow) = sbb(tmp.0[2], MODULUS.0[2], borrow); + let (_, borrow) = sbb(tmp.0[3], MODULUS.0[3], borrow); + + // If the element is smaller than MODULUS then the + // subtraction will underflow, producing a borrow value + // of 0xffff...ffff. Otherwise, it'll be zero. + let is_some = (borrow as u8) & 1; + + // Convert to Montgomery form by computing + // (a.R^0 * R^2) / R = a.R + tmp *= &R2; + + CtOption::new(tmp, Choice::from(is_some)) + } + + /// Converts an element of `Fq` into a byte representation in + /// little-endian byte order. + fn to_bytes(&self) -> [u8; 32] { + // Turn into canonical form by computing + // (a.R) / R = a + let tmp = Fq::montgomery_reduce(self.0[0], self.0[1], self.0[2], self.0[3], 0, 0, 0, 0); + + let mut res = [0; 32]; + res[0..8].copy_from_slice(&tmp.0[0].to_le_bytes()); + res[8..16].copy_from_slice(&tmp.0[1].to_le_bytes()); + res[16..24].copy_from_slice(&tmp.0[2].to_le_bytes()); + res[24..32].copy_from_slice(&tmp.0[3].to_le_bytes()); + + res + } + + /// Converts a 512-bit little endian integer into + /// a `Fq` by reducing by the modulus. + fn from_bytes_wide(bytes: &[u8; 64]) -> Fq { + Fq::from_u512([ + u64::from_le_bytes(bytes[0..8].try_into().unwrap()), + u64::from_le_bytes(bytes[8..16].try_into().unwrap()), + u64::from_le_bytes(bytes[16..24].try_into().unwrap()), + u64::from_le_bytes(bytes[24..32].try_into().unwrap()), + u64::from_le_bytes(bytes[32..40].try_into().unwrap()), + u64::from_le_bytes(bytes[40..48].try_into().unwrap()), + u64::from_le_bytes(bytes[48..56].try_into().unwrap()), + u64::from_le_bytes(bytes[56..64].try_into().unwrap()), + ]) + } + + fn get_lower_128(&self) -> u128 { + let tmp = Fq::montgomery_reduce(self.0[0], self.0[1], self.0[2], self.0[3], 0, 0, 0, 0); + + u128::from(tmp.0[0]) | (u128::from(tmp.0[1]) << 64) + } +} + +#[test] +fn test_inv() { + // Compute -(r^{-1} mod 2^64) mod 2^64 by exponentiating + // by totient(2**64) - 1 + + let mut inv = 1u64; + for _ in 0..63 { + inv = inv.wrapping_mul(inv); + inv = inv.wrapping_mul(MODULUS.0[0]); + } + inv = inv.wrapping_neg(); + + assert_eq!(inv, INV); +} + +#[test] +fn test_zeta() { + assert_eq!( + format!("{:?}", Fq::ZETA), + "0x36c66d3a1e049a5887ad8b5ff9731ffe69cf8de720e52ec14394c2bd148fa4fd" + ); + let a = Fq::ZETA; + assert!(bool::from(a != Fq::one())); + let b = a * a; + assert!(bool::from(b != Fq::one())); + let c = b * a; + assert!(bool::from(c == Fq::one())); +} + +#[test] +fn test_inv_root_of_unity() { + assert_eq!(Fq::ROOT_OF_UNITY_INV, Fq::ROOT_OF_UNITY.invert().unwrap()); +} + +#[test] +fn test_inv_2() { + assert_eq!(Fq::TWO_INV, Fq::from(2).invert().unwrap()); +} diff --git a/src/arithmetic/macros.rs b/src/arithmetic/macros.rs new file mode 100644 index 0000000..8145c99 --- /dev/null +++ b/src/arithmetic/macros.rs @@ -0,0 +1,153 @@ +macro_rules! impl_add_binop_specify_output { + ($lhs:ident, $rhs:ident, $output:ident) => { + impl<'b> Add<&'b $rhs> for $lhs { + type Output = $output; + + #[inline] + fn add(self, rhs: &'b $rhs) -> $output { + &self + rhs + } + } + + impl<'a> Add<$rhs> for &'a $lhs { + type Output = $output; + + #[inline] + fn add(self, rhs: $rhs) -> $output { + self + &rhs + } + } + + impl Add<$rhs> for $lhs { + type Output = $output; + + #[inline] + fn add(self, rhs: $rhs) -> $output { + &self + &rhs + } + } + }; +} + +macro_rules! impl_sub_binop_specify_output { + ($lhs:ident, $rhs:ident, $output:ident) => { + impl<'b> Sub<&'b $rhs> for $lhs { + type Output = $output; + + #[inline] + fn sub(self, rhs: &'b $rhs) -> $output { + &self - rhs + } + } + + impl<'a> Sub<$rhs> for &'a $lhs { + type Output = $output; + + #[inline] + fn sub(self, rhs: $rhs) -> $output { + self - &rhs + } + } + + impl Sub<$rhs> for $lhs { + type Output = $output; + + #[inline] + fn sub(self, rhs: $rhs) -> $output { + &self - &rhs + } + } + }; +} + +macro_rules! impl_binops_additive_specify_output { + ($lhs:ident, $rhs:ident, $output:ident) => { + impl_add_binop_specify_output!($lhs, $rhs, $output); + impl_sub_binop_specify_output!($lhs, $rhs, $output); + }; +} + +macro_rules! impl_binops_multiplicative_mixed { + ($lhs:ident, $rhs:ident, $output:ident) => { + impl<'b> Mul<&'b $rhs> for $lhs { + type Output = $output; + + #[inline] + fn mul(self, rhs: &'b $rhs) -> $output { + &self * rhs + } + } + + impl<'a> Mul<$rhs> for &'a $lhs { + type Output = $output; + + #[inline] + fn mul(self, rhs: $rhs) -> $output { + self * &rhs + } + } + + impl Mul<$rhs> for $lhs { + type Output = $output; + + #[inline] + fn mul(self, rhs: $rhs) -> $output { + &self * &rhs + } + } + }; +} + +macro_rules! impl_binops_additive { + ($lhs:ident, $rhs:ident) => { + impl_binops_additive_specify_output!($lhs, $rhs, $lhs); + + impl SubAssign<$rhs> for $lhs { + #[inline] + fn sub_assign(&mut self, rhs: $rhs) { + *self = &*self - &rhs; + } + } + + impl AddAssign<$rhs> for $lhs { + #[inline] + fn add_assign(&mut self, rhs: $rhs) { + *self = &*self + &rhs; + } + } + + impl<'b> SubAssign<&'b $rhs> for $lhs { + #[inline] + fn sub_assign(&mut self, rhs: &'b $rhs) { + *self = &*self - rhs; + } + } + + impl<'b> AddAssign<&'b $rhs> for $lhs { + #[inline] + fn add_assign(&mut self, rhs: &'b $rhs) { + *self = &*self + rhs; + } + } + }; +} + +macro_rules! impl_binops_multiplicative { + ($lhs:ident, $rhs:ident) => { + impl_binops_multiplicative_mixed!($lhs, $rhs, $lhs); + + impl MulAssign<$rhs> for $lhs { + #[inline] + fn mul_assign(&mut self, rhs: $rhs) { + *self = &*self * &rhs; + } + } + + impl<'b> MulAssign<&'b $rhs> for $lhs { + #[inline] + fn mul_assign(&mut self, rhs: &'b $rhs) { + *self = &*self * rhs; + } + } + }; +} diff --git a/src/lib.rs b/src/lib.rs new file mode 100644 index 0000000..78f2364 --- /dev/null +++ b/src/lib.rs @@ -0,0 +1,19 @@ +//! # halo2 + +#![allow( + clippy::op_ref, + clippy::assign_op_pattern, + clippy::too_many_arguments, + clippy::suspicious_arithmetic_impl, + clippy::many_single_char_names, + clippy::same_item_push +)] +#![deny(intra_doc_link_resolution_failure)] +#![deny(missing_debug_implementations)] +#![deny(missing_docs)] +#![deny(unsafe_code)] + +pub mod arithmetic; +pub mod plonk; +pub mod polycommit; +pub mod transcript; diff --git a/src/plonk.rs b/src/plonk.rs new file mode 100644 index 0000000..66f3e32 --- /dev/null +++ b/src/plonk.rs @@ -0,0 +1,145 @@ +//! This module provides an implementation of a variant of (Turbo)[PLONK][plonk] +//! that is designed specifically for the polynomial commitment scheme described +//! in the [Halo][halo] paper. +//! +//! [halo]: https://eprint.iacr.org/2019/1021 +//! [plonk]: https://eprint.iacr.org/2019/953 + +use crate::arithmetic::CurveAffine; +use crate::polycommit::OpeningProof; +use crate::transcript::Hasher; + +#[macro_use] +mod circuit; +mod domain; +mod prover; +mod srs; +mod verifier; + +pub use circuit::*; +pub use prover::*; +pub use srs::*; +pub use verifier::*; + +use domain::EvaluationDomain; + +// TODO: remove this +const GATE_DEGREE: u32 = 3; + +/// This is a structured reference string (SRS) that is (deterministically) +/// computed from a specific circuit and parameters for the polynomial +/// commitment scheme. +#[derive(Debug)] +pub struct SRS { + sa: (Vec, Vec), + sb: (Vec, Vec), + sc: (Vec, Vec), + sd: (Vec, Vec), + sm: (Vec, Vec), + sa_commitment: C, + sb_commitment: C, + sc_commitment: C, + sd_commitment: C, + sm_commitment: C, + domain: EvaluationDomain, +} + +/// This is an object which represents a (Turbo)PLONK proof. +#[derive(Debug, Clone)] +pub struct Proof { + a_commitment: C, + b_commitment: C, + c_commitment: C, + d_commitment: C, + h_commitments: Vec, + a_eval_x: C::Scalar, + b_eval_x: C::Scalar, + c_eval_x: C::Scalar, + d_eval_x: C::Scalar, + sa_eval_x: C::Scalar, + sb_eval_x: C::Scalar, + sc_eval_x: C::Scalar, + sd_eval_x: C::Scalar, + sm_eval_x: C::Scalar, + h_evals_x: Vec, + opening: OpeningProof, +} + +/// This is an error that could occur during proving or circuit synthesis. +// TODO: these errors need to be cleaned up +#[derive(Debug)] +pub enum Error { + /// This is an error that can occur during synthesis of the circuit, for + /// example, when the witness is not present. + SynthesisError, + /// The structured reference string or the parameters are not compatible + /// with the circuit being synthesized. + IncompatibleParams, + /// The constraint system is not satisfied. + ConstraintSystemFailure, +} + +fn hash_point>( + transcript: &mut H, + point: &C, +) -> Result<(), Error> { + let tmp = point.get_xy(); + if bool::from(tmp.is_none()) { + return Err(Error::SynthesisError); + }; + let tmp = tmp.unwrap(); + transcript.absorb(tmp.0); + transcript.absorb(tmp.1); + Ok(()) +} + +#[test] +fn test_proving() { + use crate::arithmetic::{EqAffine, Field, Fp, Fq}; + use crate::polycommit::Params; + use crate::transcript::DummyHash; + const K: u32 = 5; + + // Initialize the polynomial commitment parameters + let params: Params = Params::new::>(K); + + struct MyCircuit { + a: Option, + } + + impl Circuit for MyCircuit { + fn synthesize(&self, cs: &mut impl ConstraintSystem) -> Result<(), Error> { + for _ in 0..10 { + let (_, _, _, _) = cs.multiply(|| { + let a = self.a.ok_or(Error::SynthesisError)?; + let a2 = a.square(); + Ok((a, a, a2)) + })?; + //cs.copy(a, b); + let (_, _, _, _) = cs.add(|| { + let a = self.a.ok_or(Error::SynthesisError)?; + let a2 = a.square(); + let a3 = a + a2; + Ok((a, a2, a3)) + })?; + //cs.copy(a, d); + //cs.copy(c, e); + } + + Ok(()) + } + } + + let circuit: MyCircuit = MyCircuit { + a: Some((-Fp::from_u64(2) + Fp::ROOT_OF_UNITY).pow(&[100, 0, 0, 0])), + }; + + // Initialize the SRS + let srs = SRS::generate(¶ms, &circuit).expect("SRS generation should not fail"); + + // Create a proof + let proof = Proof::create::, DummyHash, _>(¶ms, &srs, &circuit) + .expect("proof generation should not fail"); + + assert!(proof.verify::, DummyHash>(¶ms, &srs)); +} diff --git a/src/plonk/circuit.rs b/src/plonk/circuit.rs new file mode 100644 index 0000000..c0168d2 --- /dev/null +++ b/src/plonk/circuit.rs @@ -0,0 +1,66 @@ +use super::Error; + +use crate::arithmetic::Field; + +/// This represents a PLONK wire, which could be a fixed (selector) wire or an +/// advice wire. +#[derive(Debug)] +pub enum Wire { + /// A wires + A(usize), + /// B wires + B(usize), + /// C wires + C(usize), + /// D wires + D(usize), +} + +/// This trait allows a [`Circuit`] to direct some backend to assign a witness +/// for a constraint system. +pub trait ConstraintSystem { + /// Creates a gate. + fn create_gate( + &mut self, + sa: F, + sb: F, + sc: F, + sd: F, + sm: F, + f: impl Fn() -> Result<(F, F, F, F), Error>, + ) -> Result<(Wire, Wire, Wire, Wire), Error>; + + /// a * b - c = 0 + fn multiply( + &mut self, + f: impl Fn() -> Result<(F, F, F), Error>, + ) -> Result<(Wire, Wire, Wire, Wire), Error> { + self.create_gate(F::zero(), F::zero(), F::one(), F::zero(), F::one(), || { + let (a, b, c) = f()?; + Ok((a, b, c, F::zero())) + }) + } + + /// a + b - c = 0 + fn add( + &mut self, + f: impl Fn() -> Result<(F, F, F), Error>, + ) -> Result<(Wire, Wire, Wire, Wire), Error> { + self.create_gate(F::one(), F::one(), F::one(), F::zero(), F::zero(), || { + let (a, b, c) = f()?; + Ok((a, b, c, F::zero())) + }) + } + + // fn copy(&mut self, left: Wire, right: Wire); +} + +/// This is a trait that circuits provide implementations for so that the +/// backend prover can ask the circuit to synthesize using some given +/// [`ConstraintSystem`] implementation. +pub trait Circuit { + /// Given the provided `cs`, synthesize the circuit. The concrete type of + /// the caller will be different depending on the context, and they may or + /// may not expect to have a witness present. + fn synthesize(&self, cs: &mut impl ConstraintSystem) -> Result<(), Error>; +} diff --git a/src/plonk/domain.rs b/src/plonk/domain.rs new file mode 100644 index 0000000..1aa70fa --- /dev/null +++ b/src/plonk/domain.rs @@ -0,0 +1,202 @@ +use crate::arithmetic::{best_fft, parallelize, Field, Group}; + +/// This structure contains precomputed constants and other details needed for +/// performing operations on an evaluation domain of size $2^k$ in the context +/// of PLONK. +#[derive(Debug)] +pub struct EvaluationDomain { + n: u64, + k: u32, + extended_k: u32, + omega_inv: G::Scalar, + extended_omega: G::Scalar, + extended_omega_inv: G::Scalar, + g_coset: G::Scalar, + g_coset_inv: G::Scalar, + quotient_poly_degree: u64, + ifft_divisor: G::Scalar, + extended_ifft_divisor: G::Scalar, + t_evaluations: Vec, +} + +impl EvaluationDomain { + /// This constructs a new evaluation domain object (containing precomputed + /// constants) for operating on an evaluation domain of size $2^k$ and for + /// some operations over an extended domain of size $2^{k + j}$ where $j$ is + /// sufficiently large to describe the quotient polynomial depending on the + /// maximum degree of all PLONK gates. + pub fn new(gate_degree: u32, k: u32) -> Self { + // quotient_poly_degree * params.n - 1 is the degree of the quotient polynomial + let quotient_poly_degree = (gate_degree - 1) as u64; + + let n = 1u64 << k; + + // We need to work within an extended domain, not params.k but params.k + j + // such that 2^(params.k + j) is sufficiently large to describe the quotient + // polynomial. + let mut extended_k = k; + while (1 << extended_k) < (n * quotient_poly_degree) { + extended_k += 1; + } + + let mut extended_omega = G::Scalar::ROOT_OF_UNITY; + for _ in extended_k..G::Scalar::S { + extended_omega = extended_omega.square(); + } + let extended_omega = extended_omega; // 2^{j+k}'th root of unity + let extended_omega_inv = extended_omega.invert().unwrap(); + + let mut omega = extended_omega; + for _ in k..extended_k { + omega = omega.square(); + } + let omega = omega; // 2^{k}'th root of unity + let omega_inv = omega.invert().unwrap(); + + // We use zeta here because we know it generates a coset, and it's available + // already. + let g_coset = G::Scalar::ZETA; + let g_coset_inv = g_coset.square(); + + // TODO: merge these inversions together with t_evaluations batch inversion? + let ifft_divisor = G::Scalar::from_u64(1 << k).invert().unwrap(); + let extended_ifft_divisor = G::Scalar::from_u64(1 << extended_k).invert().unwrap(); + + let mut t_evaluations = Vec::with_capacity(1 << (extended_k - k)); + { + // Compute the evaluations of t(X) in the coset evaluation domain. + // We don't have to compute all of them, because it will repeat. + let orig = G::Scalar::ZETA.pow_vartime(&[n as u64, 0, 0, 0]); + let step = extended_omega.pow_vartime(&[n as u64, 0, 0, 0]); + let mut cur = orig; + loop { + t_evaluations.push(cur); + cur *= &step; + if cur == orig { + break; + } + } + assert_eq!(t_evaluations.len(), 1 << (extended_k - k)); + + // Subtract 1 from each to give us t_evaluations[i] = t(zeta * extended_omega^i) + for coeff in &mut t_evaluations { + *coeff -= &G::Scalar::one(); + } + + // Invert, because we're dividing by this polynomial. + G::Scalar::batch_invert(&mut t_evaluations); + } + + EvaluationDomain { + n, + k, + extended_k, + omega_inv, + extended_omega, + extended_omega_inv, + g_coset, + g_coset_inv, + quotient_poly_degree, + ifft_divisor, + extended_ifft_divisor, + t_evaluations, + } + } + + /// This takes us from an n-length vector into the coset evaluation domain. + /// Also returns the polynomial. + /// + /// This function will panic if the provided vector is not the correct + /// length. + pub fn obtain_coset(&self, mut a: Vec) -> (Vec, Vec) { + assert_eq!(a.len(), 1 << self.k); + + // Perform inverse FFT to obtain the polynomial in coefficient form + Self::ifft(&mut a, self.omega_inv, self.k, self.ifft_divisor); + + // Keep this polynomial around; we'll need to evaluate it at arbitrary + // points later. + let old = a.clone(); + + // Distributes powers so that an FFT will move us into the coset + // evaluation domain. + Self::distribute_powers(&mut a, self.g_coset); + + // Resize to account for the quotient polynomial's size + a.resize(1 << self.extended_k, G::group_zero()); + + // Move into coset evaluation domain + best_fft(&mut a, self.extended_omega, self.extended_k); + + (a, old) + } + + /// This takes us from the coset evaluation domain and gets us the quotient + /// polynomial coefficients. + /// + /// This function will panic if the provided vector is not the correct + /// length. + pub fn from_coset(&self, mut a: Vec) -> Vec { + assert_eq!(a.len(), 1 << self.extended_k); + + // Inverse FFT + Self::ifft( + &mut a, + self.extended_omega_inv, + self.extended_k, + self.extended_ifft_divisor, + ); + + // Distribute powers to move from coset; opposite from the + // transformation we performed earlier. + Self::distribute_powers(&mut a, self.g_coset_inv); + + // Truncate it to match the size of the quotient polynomial; the + // evaluation domain might be slightly larger than necessary because + // it always lies on a power-of-two boundary. + a.truncate((&self.n * self.quotient_poly_degree) as usize); + + a + } + + /// This divides the polynomial (in the coset domain) by the vanishing + /// polynomial. + pub fn divide_by_vanishing_poly(&self, mut h_poly: Vec) -> Vec { + assert_eq!(h_poly.len(), 1 << self.extended_k); + + // Divide to obtain the quotient polynomial in the coset evaluation + // domain. + parallelize(&mut h_poly, |h, mut index| { + for h in h { + h.group_scale(&self.t_evaluations[index % self.t_evaluations.len()]); + index += 1; + } + }); + + h_poly + } + + fn distribute_powers(mut a: &mut [G], g: G::Scalar) { + let coset_powers = [g, g.square()]; + parallelize(&mut a, |a, mut index| { + for a in a { + // Distribute powers to move into coset + let i = index % (coset_powers.len() + 1); + if i != 0 { + a.group_scale(&coset_powers[i - 1]); + } + index += 1; + } + }); + } + + fn ifft(a: &mut [G], omega_inv: G::Scalar, log_n: u32, divisor: G::Scalar) { + best_fft(a, omega_inv, log_n); + parallelize(a, |a, _| { + for a in a { + // Finish iFFT + a.group_scale(&divisor); + } + }); + } +} diff --git a/src/plonk/prover.rs b/src/plonk/prover.rs new file mode 100644 index 0000000..d06119d --- /dev/null +++ b/src/plonk/prover.rs @@ -0,0 +1,266 @@ +use super::{ + circuit::{Circuit, ConstraintSystem, Wire}, + hash_point, Error, Proof, SRS, +}; +use crate::arithmetic::{ + eval_polynomial, get_challenge_scalar, Challenge, Curve, CurveAffine, Field, +}; +use crate::polycommit::Params; +use crate::transcript::Hasher; + +impl Proof { + /// This creates a proof for the provided `circuit` when given the public + /// parameters `params` and the structured reference string `srs` that was + /// previously computed for the same circuit. + pub fn create< + HBase: Hasher, + HScalar: Hasher, + ConcreteCircuit: Circuit, + >( + params: &Params, + srs: &SRS, + circuit: &ConcreteCircuit, + ) -> Result { + struct WitnessCollection { + a: Vec, + b: Vec, + c: Vec, + d: Vec, + sa: Vec, + sb: Vec, + sc: Vec, + sd: Vec, + sm: Vec, + } + + impl ConstraintSystem for WitnessCollection { + fn create_gate( + &mut self, + sa: F, + sb: F, + sc: F, + sd: F, + sm: F, + f: impl Fn() -> Result<(F, F, F, F), Error>, + ) -> Result<(Wire, Wire, Wire, Wire), Error> { + let (a, b, c, d) = f()?; + let tmp = Ok(( + Wire::A(self.a.len()), + Wire::B(self.a.len()), + Wire::C(self.a.len()), + Wire::D(self.a.len()), + )); + self.a.push(a); + self.b.push(b); + self.c.push(c); + self.d.push(d); + self.sa.push(sa); + self.sb.push(sb); + self.sc.push(sc); + self.sd.push(sd); + self.sm.push(sm); + tmp + } + // fn copy(&mut self, left: Wire, right: Wire) { + // unimplemented!() + // } + } + + let mut witness = WitnessCollection { + a: vec![], + b: vec![], + c: vec![], + d: vec![], + sa: vec![], + sb: vec![], + sc: vec![], + sd: vec![], + sm: vec![], + }; + + // Synthesize the circuit to obtain the witness and other information. + circuit.synthesize(&mut witness)?; + + // Create a transcript for obtaining Fiat-Shamir challenges. + let mut transcript = HBase::init(C::Base::one()); + + if witness.a.len() > params.n as usize { + // The polynomial commitment does not support a high enough degree + // polynomial to commit to our wires because this circuit has too + // many gates. + return Err(Error::IncompatibleParams); + } + + witness.a.resize(params.n as usize, C::Scalar::zero()); + witness.b.resize(params.n as usize, C::Scalar::zero()); + witness.c.resize(params.n as usize, C::Scalar::zero()); + witness.d.resize(params.n as usize, C::Scalar::zero()); + witness.sa.resize(params.n as usize, C::Scalar::zero()); + witness.sb.resize(params.n as usize, C::Scalar::zero()); + witness.sc.resize(params.n as usize, C::Scalar::zero()); + witness.sd.resize(params.n as usize, C::Scalar::zero()); + witness.sm.resize(params.n as usize, C::Scalar::zero()); + + // Compute commitments to the various wire values + let a_blind = C::Scalar::one(); // TODO: not random + let b_blind = C::Scalar::one(); // TODO: not random + let c_blind = C::Scalar::one(); // TODO: not random + let d_blind = C::Scalar::one(); // TODO: not random + let a_commitment = params.commit_lagrange(&witness.a, a_blind).to_affine(); + let b_commitment = params.commit_lagrange(&witness.b, b_blind).to_affine(); + let c_commitment = params.commit_lagrange(&witness.c, c_blind).to_affine(); + let d_commitment = params.commit_lagrange(&witness.d, d_blind).to_affine(); + + hash_point(&mut transcript, &a_commitment)?; + hash_point(&mut transcript, &b_commitment)?; + hash_point(&mut transcript, &c_commitment)?; + hash_point(&mut transcript, &d_commitment)?; + + let domain = &srs.domain; + + let (a_coset, a_poly) = domain.obtain_coset(witness.a); + let (b_coset, b_poly) = domain.obtain_coset(witness.b); + let (c_coset, c_poly) = domain.obtain_coset(witness.c); + let (d_coset, d_poly) = domain.obtain_coset(witness.d); + + // (a * sa) + (b * sb) + (a * sm * b) + (d * sd) - (c * sc) + let mut h_poly = Vec::with_capacity(a_coset.len()); + for ((((((((a, b), c), d), sa), sb), sc), sd), sm) in a_coset + .iter() + .zip(b_coset.iter()) + .zip(c_coset.iter()) + .zip(d_coset.iter()) + .zip(srs.sa.0.iter()) + .zip(srs.sb.0.iter()) + .zip(srs.sc.0.iter()) + .zip(srs.sd.0.iter()) + .zip(srs.sm.0.iter()) + { + h_poly.push((*a) * sa + &((*b) * sb) + &((*a) * sm * b) + &((*d) * sd) - &((*c) * sc)); + } + + // Divide by t(X) = X^{params.n} - 1. + let h_poly = domain.divide_by_vanishing_poly(h_poly); + + // Obtain final h(X) polynomial + let h_poly = domain.from_coset(h_poly); + + // Split h(X) up into pieces + let h_pieces = h_poly + .chunks_exact(params.n as usize) + .map(|v| v.to_vec()) + .collect::>(); + drop(h_poly); + let h_blinds = vec![C::Scalar::one(); h_pieces.len()]; // TODO: not random + + // Compute commitments to each h(X) piece + let h_commitments: Vec<_> = h_pieces + .iter() + .zip(h_blinds.iter()) + .map(|(h_piece, blind)| params.commit(&h_piece, *blind).to_affine()) + .collect(); + + // Hash each h(X) piece + for c in h_commitments.iter() { + hash_point(&mut transcript, c)?; + } + + let x: C::Scalar = get_challenge_scalar(Challenge(transcript.squeeze().get_lower_128())); + + // Evaluate polynomials at x + let a_eval_x = eval_polynomial(&a_poly, x); + let b_eval_x = eval_polynomial(&b_poly, x); + let c_eval_x = eval_polynomial(&c_poly, x); + let d_eval_x = eval_polynomial(&d_poly, x); + let sa_eval_x = eval_polynomial(&srs.sa.1, x); + let sb_eval_x = eval_polynomial(&srs.sb.1, x); + let sc_eval_x = eval_polynomial(&srs.sc.1, x); + let sd_eval_x = eval_polynomial(&srs.sd.1, x); + let sm_eval_x = eval_polynomial(&srs.sm.1, x); + + let h_evals_x: Vec<_> = h_pieces + .iter() + .map(|poly| eval_polynomial(poly, x)) + .collect(); + + // We set up a second transcript on the scalar field to hash in openings of + // our polynomial commitments. + let mut transcript_scalar = HScalar::init(C::Scalar::one()); + transcript_scalar.absorb(a_eval_x); + transcript_scalar.absorb(b_eval_x); + transcript_scalar.absorb(c_eval_x); + transcript_scalar.absorb(d_eval_x); + transcript_scalar.absorb(sa_eval_x); + transcript_scalar.absorb(sb_eval_x); + transcript_scalar.absorb(sc_eval_x); + transcript_scalar.absorb(sd_eval_x); + transcript_scalar.absorb(sm_eval_x); + + // Hash each h(x) piece + for eval in h_evals_x.iter() { + transcript_scalar.absorb(*eval); + } + + let transcript_scalar_point = + C::Base::from_bytes(&(transcript_scalar.squeeze()).to_bytes()).unwrap(); + transcript.absorb(transcript_scalar_point); + + let y: C::Scalar = get_challenge_scalar(Challenge(transcript.squeeze().get_lower_128())); + + let mut q_commitment = h_commitments[0].clone().to_projective(); + let mut q_poly = h_pieces[0].clone(); + let mut q_blind = h_blinds[0]; + { + let mut accumulate = |poly: &[_], blind: &C::Scalar, commitment: C| { + for (a, q) in poly.iter().zip(q_poly.iter_mut()) { + *q = (*q * &y) + a; + } + q_commitment = (q_commitment * y) + &commitment.to_projective(); + q_blind = (q_blind * &y) + blind; + }; + + for ((poly, blind), commitment) in h_pieces + .iter() + .zip(h_blinds.iter()) + .zip(h_commitments.iter()) + .skip(1) + { + accumulate(&poly, blind, *commitment); + } + + accumulate(&a_poly, &a_blind, a_commitment); + accumulate(&b_poly, &b_blind, b_commitment); + accumulate(&c_poly, &c_blind, c_commitment); + accumulate(&d_poly, &d_blind, d_commitment); + accumulate(&srs.sa.1, &Field::one(), srs.sa_commitment); + accumulate(&srs.sb.1, &Field::one(), srs.sb_commitment); + accumulate(&srs.sc.1, &Field::one(), srs.sc_commitment); + accumulate(&srs.sd.1, &Field::one(), srs.sd_commitment); + accumulate(&srs.sm.1, &Field::one(), srs.sm_commitment); + } + + // Let's prove that the q_commitment opens at x to the expected value. + let opening = params + .create_proof(&mut transcript, &q_poly, q_blind, x) + .map_err(|_| Error::ConstraintSystemFailure)?; + + Ok(Proof { + a_commitment, + b_commitment, + c_commitment, + d_commitment, + h_commitments, + a_eval_x, + b_eval_x, + c_eval_x, + d_eval_x, + sa_eval_x, + sb_eval_x, + sc_eval_x, + sd_eval_x, + sm_eval_x, + h_evals_x, + opening, + }) + } +} diff --git a/src/plonk/srs.rs b/src/plonk/srs.rs new file mode 100644 index 0000000..0a9acdb --- /dev/null +++ b/src/plonk/srs.rs @@ -0,0 +1,105 @@ +use super::{ + circuit::{Circuit, ConstraintSystem, Wire}, + domain::EvaluationDomain, + Error, GATE_DEGREE, SRS, +}; +use crate::arithmetic::{Curve, CurveAffine, Field}; +use crate::polycommit::Params; + +impl SRS { + /// This generates a structured reference string for the provided `circuit` + /// and `params`. + pub fn generate>( + params: &Params, + circuit: &ConcreteCircuit, + ) -> Result { + struct Assembly { + sa: Vec, + sb: Vec, + sc: Vec, + sd: Vec, + sm: Vec, + } + + impl ConstraintSystem for Assembly { + fn create_gate( + &mut self, + sa: F, + sb: F, + sc: F, + sd: F, + sm: F, + _: impl Fn() -> Result<(F, F, F, F), Error>, + ) -> Result<(Wire, Wire, Wire, Wire), Error> { + let tmp = Ok(( + Wire::A(self.sa.len()), + Wire::B(self.sa.len()), + Wire::C(self.sa.len()), + Wire::D(self.sa.len()), + )); + self.sa.push(sa); + self.sb.push(sb); + self.sc.push(sc); + self.sd.push(sd); + self.sm.push(sm); + tmp + } + } + + let mut assembly: Assembly = Assembly { + sa: vec![], + sb: vec![], + sc: vec![], + sd: vec![], + sm: vec![], + }; + + // Synthesize the circuit to obtain SRS + circuit.synthesize(&mut assembly)?; + + assembly.sa.resize(params.n as usize, C::Scalar::zero()); + assembly.sb.resize(params.n as usize, C::Scalar::zero()); + assembly.sc.resize(params.n as usize, C::Scalar::zero()); + assembly.sd.resize(params.n as usize, C::Scalar::zero()); + assembly.sm.resize(params.n as usize, C::Scalar::zero()); + + // Compute commitments to the fixed wire values + let sa_commitment = params + .commit_lagrange(&assembly.sa, C::Scalar::one()) + .to_affine(); + let sb_commitment = params + .commit_lagrange(&assembly.sb, C::Scalar::one()) + .to_affine(); + let sc_commitment = params + .commit_lagrange(&assembly.sc, C::Scalar::one()) + .to_affine(); + let sd_commitment = params + .commit_lagrange(&assembly.sd, C::Scalar::one()) + .to_affine(); + let sm_commitment = params + .commit_lagrange(&assembly.sm, C::Scalar::one()) + .to_affine(); + + let domain = EvaluationDomain::new(GATE_DEGREE, params.k); + + let sa = domain.obtain_coset(assembly.sa); + let sb = domain.obtain_coset(assembly.sb); + let sc = domain.obtain_coset(assembly.sc); + let sd = domain.obtain_coset(assembly.sd); + let sm = domain.obtain_coset(assembly.sm); + + Ok(SRS { + sa, + sb, + sc, + sd, + sm, + sa_commitment, + sb_commitment, + sc_commitment, + sd_commitment, + sm_commitment, + domain, + }) + } +} diff --git a/src/plonk/verifier.rs b/src/plonk/verifier.rs new file mode 100644 index 0000000..4759fd4 --- /dev/null +++ b/src/plonk/verifier.rs @@ -0,0 +1,108 @@ +use super::{hash_point, Proof, SRS}; +use crate::arithmetic::{get_challenge_scalar, Challenge, Curve, CurveAffine, Field}; +use crate::polycommit::Params; +use crate::transcript::Hasher; + +impl Proof { + /// Returns + pub fn verify, HScalar: Hasher>( + &self, + params: &Params, + srs: &SRS, + ) -> bool { + // Create a transcript for obtaining Fiat-Shamir challenges. + let mut transcript = HBase::init(C::Base::one()); + + hash_point(&mut transcript, &self.a_commitment) + .expect("proof cannot contain points at infinity"); + hash_point(&mut transcript, &self.b_commitment) + .expect("proof cannot contain points at infinity"); + hash_point(&mut transcript, &self.c_commitment) + .expect("proof cannot contain points at infinity"); + hash_point(&mut transcript, &self.d_commitment) + .expect("proof cannot contain points at infinity"); + + for c in &self.h_commitments { + hash_point(&mut transcript, c).expect("proof cannot contain points at infinity"); + } + + let x: C::Scalar = get_challenge_scalar(Challenge(transcript.squeeze().get_lower_128())); + + // We set up a second transcript on the scalar field to hash in openings of + // our polynomial commitments. + let mut transcript_scalar = HScalar::init(C::Scalar::one()); + transcript_scalar.absorb(self.a_eval_x); + transcript_scalar.absorb(self.b_eval_x); + transcript_scalar.absorb(self.c_eval_x); + transcript_scalar.absorb(self.d_eval_x); + transcript_scalar.absorb(self.sa_eval_x); + transcript_scalar.absorb(self.sb_eval_x); + transcript_scalar.absorb(self.sc_eval_x); + transcript_scalar.absorb(self.sd_eval_x); + transcript_scalar.absorb(self.sm_eval_x); + + for eval in &self.h_evals_x { + transcript_scalar.absorb(*eval); + } + + let transcript_scalar_point = + C::Base::from_bytes(&(transcript_scalar.squeeze()).to_bytes()).unwrap(); + transcript.absorb(transcript_scalar_point); + + let y: C::Scalar = get_challenge_scalar(Challenge(transcript.squeeze().get_lower_128())); + + let mut q_commitment = self.h_commitments[0].clone().to_projective(); + let mut expected_opening = self.h_evals_x[0]; + { + let mut accumulate = |commitment: C, opening: C::Scalar| { + q_commitment = commitment.to_projective() + &(q_commitment * y); + expected_opening = opening + &(expected_opening * &y); + }; + + for (commitment, eval) in self.h_commitments.iter().zip(self.h_evals_x.iter()).skip(1) { + accumulate(*commitment, *eval); + } + + accumulate(self.a_commitment, self.a_eval_x); + accumulate(self.b_commitment, self.b_eval_x); + accumulate(self.c_commitment, self.c_eval_x); + accumulate(self.d_commitment, self.d_eval_x); + accumulate(srs.sa_commitment, self.sa_eval_x); + accumulate(srs.sb_commitment, self.sb_eval_x); + accumulate(srs.sc_commitment, self.sc_eval_x); + accumulate(srs.sd_commitment, self.sd_eval_x); + accumulate(srs.sm_commitment, self.sm_eval_x); + } + let q_commitment = q_commitment.to_affine(); + + let xn = x.pow(&[params.n as u64, 0, 0, 0]); + + // Compute the expected h(x) value + let mut h_eval_x = C::Scalar::zero(); + let mut cur = C::Scalar::one(); + for eval in &self.h_evals_x { + h_eval_x += &(cur * eval); + cur *= &xn; + } + + // Check that the circuit is satisfied. + // (a * sa) + (b * sb) + (a * sm * b) + (d * sd) - (c * sc) + if self.a_eval_x * &self.sa_eval_x + + &(self.b_eval_x * &self.sb_eval_x) + + &(self.a_eval_x * &self.sm_eval_x * &self.b_eval_x) + + &(self.d_eval_x * &self.sd_eval_x) + - &(self.c_eval_x * &self.sc_eval_x) + != h_eval_x * &(xn - &C::Scalar::one()) + { + return false; + } + + params.verify_proof( + &self.opening, + &mut transcript, + x, + &q_commitment, + expected_opening, + ) + } +} diff --git a/src/polycommit.rs b/src/polycommit.rs new file mode 100644 index 0000000..29f6c5b --- /dev/null +++ b/src/polycommit.rs @@ -0,0 +1,600 @@ +//! This module contains an implementation of the polynomial commitment scheme +//! described in the [Halo][halo] paper. +//! +//! [halo]: https://eprint.iacr.org/2019/1021 + +use crate::arithmetic::{ + best_fft, best_multiexp, compute_inner_product, get_challenge_scalar, parallelize, Challenge, + Curve, CurveAffine, Field, +}; +use crate::transcript::Hasher; + +/// This is a proof object for the polynomial commitment scheme opening. +#[derive(Debug, Clone)] +pub struct OpeningProof { + fork: u8, + rounds: Vec<(C, C)>, + delta: C, + z1: C::Scalar, + z2: C::Scalar, +} + +/// These are the public parameters for the polynomial commitment scheme. +#[derive(Debug)] +pub struct Params { + pub(crate) k: u32, + pub(crate) n: u64, + pub(crate) g: Vec, + pub(crate) g_lagrange: Vec, + pub(crate) h: C, +} + +impl Params { + /// Initializes parameters for the curve, given a random oracle to draw + /// points from. + pub fn new>(k: u32) -> Self { + // This is usually a limitation on the curve, but we also want 32-bit + // architectures to be supported. + assert!(k < 32); + // No goofy hardware please. + assert!(core::mem::size_of::() >= 4); + + let n: u64 = 1 << k; + + let g = { + let hasher = &H::init(C::Base::zero()); + + let mut g = Vec::with_capacity(n as usize); + g.resize(n as usize, C::zero()); + + parallelize(&mut g, move |g, start| { + let mut cur_value = C::Base::from(start as u64); + for g in g.iter_mut() { + let mut hasher = hasher.clone(); + hasher.absorb(cur_value); + cur_value += &C::Base::one(); + loop { + let x = hasher.squeeze().to_bytes(); + let p = C::from_bytes(&x); + if bool::from(p.is_some()) { + *g = p.unwrap(); + break; + } + } + } + }); + + g + }; + + // Let's evaluate all of the Lagrange basis polynomials + // using an inverse FFT. + let mut alpha_inv = C::Scalar::ROOT_OF_UNITY_INV; + for _ in k..C::Scalar::S { + alpha_inv = alpha_inv.square(); + } + let mut g_lagrange_projective = g.iter().map(|g| g.to_projective()).collect::>(); + best_fft(&mut g_lagrange_projective, alpha_inv, k); + let minv = C::Scalar::TWO_INV.pow_vartime(&[k as u64, 0, 0, 0]); + parallelize(&mut g_lagrange_projective, |g, _| { + for g in g.iter_mut() { + *g *= minv; + } + }); + + let g_lagrange = { + let mut g_lagrange = vec![C::zero(); n as usize]; + parallelize(&mut g_lagrange, |g_lagrange, starts| { + C::Projective::batch_to_affine( + &g_lagrange_projective[starts..(starts + g_lagrange.len())], + g_lagrange, + ); + }); + drop(g_lagrange_projective); + g_lagrange + }; + + let h = { + let mut hasher = H::init(C::Base::zero()); + hasher.absorb(-C::Base::one()); + let x = hasher.squeeze().to_bytes(); + let p = C::from_bytes(&x); + p.unwrap() + }; + + Params { + k, + n, + g, + g_lagrange, + h, + } + } + + /// This computes a commitment to a polynomial described by the provided + /// slice of coefficients. The commitment will be blinded by the blinding + /// factor `r`. + pub fn commit(&self, poly: &[C::Scalar], r: C::Scalar) -> C::Projective { + let mut tmp_scalars = Vec::with_capacity(poly.len() + 1); + let mut tmp_bases = Vec::with_capacity(poly.len() + 1); + + tmp_scalars.extend(poly.iter()); + tmp_scalars.push(r); + + tmp_bases.extend(self.g.iter()); + tmp_bases.push(self.h); + + best_multiexp::(&tmp_scalars, &tmp_bases) + } + + /// This commits to a polynomial using its evaluations over the $2^k$ size + /// evaluation domain. The commitment will be blinded by the blinding factor + /// `r`. + pub fn commit_lagrange(&self, poly: &[C::Scalar], r: C::Scalar) -> C::Projective { + let mut tmp_scalars = Vec::with_capacity(poly.len() + 1); + let mut tmp_bases = Vec::with_capacity(poly.len() + 1); + + tmp_scalars.extend(poly.iter()); + tmp_scalars.push(r); + + tmp_bases.extend(self.g_lagrange.iter()); + tmp_bases.push(self.h); + + best_multiexp::(&tmp_scalars, &tmp_bases) + } + + /// Create a polynomial commitment opening proof for the polynomial defined + /// by the coefficients `px`, the blinding factor `blind` used for the + /// polynomial commitment, and the point `x` that the polynomial is + /// evaluated at. + /// + /// This function will panic if the provided polynomial is too large with + /// respect to the polynomial commitment parameters. + /// + /// **Important:** This function assumes that the provided `transcript` has + /// already seen the common inputs: the polynomial commitment P, the claimed + /// opening v, and the point x. It's probably also nice for the transcript + /// to have seen the elliptic curve description and the SRS, if you want to + /// be rigorous. + pub fn create_proof>( + &self, + transcript: &mut H, + px: &[C::Scalar], + mut blind: C::Scalar, + x: C::Scalar, + ) -> Result, ()> { + // We're limited to polynomials of degree n - 1. + assert!(px.len() <= self.n as usize); + + let mut fork = 0; + + // TODO: remove this hack and force the caller to deal with it + loop { + let mut transcript = transcript.clone(); + transcript.absorb(C::Base::from_u64(fork as u64)); + let u_x = transcript.squeeze(); + // y^2 = x^3 + B + let u_y2 = u_x.square() * &u_x + &C::b(); + let u_y = u_y2.deterministic_sqrt(); + + if u_y.is_none() { + fork += 1; + } else { + break; + } + } + + transcript.absorb(C::Base::from_u64(fork as u64)); + + // Compute U + let u = { + let u_x = transcript.squeeze(); + // y^2 = x^3 + B + let u_y2 = u_x.square() * &u_x + &C::b(); + let u_y = u_y2.deterministic_sqrt().unwrap(); + + C::from_xy(u_x, u_y).unwrap() + }; + + // Initialize the vector `a` as the coefficients of the polynomial, + // rounding up to the parameters. + let mut a = px.to_vec(); + a.resize(self.n as usize, C::Scalar::zero()); + + // Initialize the vector `b` as the powers of `x`. The inner product of + // `a` and `b` is the evaluation of the polynomial at `x`. + let mut b = Vec::with_capacity(1 << self.k); + { + let mut cur = C::Scalar::one(); + for _ in 0..(1 << self.k) { + b.push(cur); + cur *= &x; + } + } + + // Initialize the vector `G` from the SRS. We'll be progressively + // collapsing this vector into smaller and smaller vectors until it is + // of length 1. + let mut g = self.g.clone(); + + // Perform the inner product argument, round by round. + let mut rounds = Vec::with_capacity(self.k as usize); + for k in (1..=self.k).rev() { + let half = 1 << (k - 1); // half the length of `a`, `b`, `G` + + // Compute L, R + // + // TODO: If we modify multiexp to take "extra" bases, we could speed + // this piece up a bit by combining the multiexps. + let l = best_multiexp(&a[0..half], &g[half..]); + let r = best_multiexp(&a[half..], &g[0..half]); + let value_l = compute_inner_product(&a[0..half], &b[half..]); + let value_r = compute_inner_product(&a[half..], &b[0..half]); + let mut l_randomness = C::Scalar::random(); + let r_randomness = C::Scalar::random(); + let l = l + &best_multiexp(&[value_l, l_randomness], &[u, self.h]); + let r = r + &best_multiexp(&[value_r, r_randomness], &[u, self.h]); + let mut l = l.to_affine(); + let r = r.to_affine(); + + let challenge = loop { + // We'll fork the transcript and adjust our randomness + // until the challenge is a square. + let mut transcript = transcript.clone(); + + // We expect these to not be points at infinity due to the randomness. + let (l_x, l_y) = l.get_xy().unwrap(); + let (r_x, r_y) = r.get_xy().unwrap(); + + // Feed L and R into the cloned transcript... + transcript.absorb(l_x); + transcript.absorb(l_y); + transcript.absorb(r_x); + transcript.absorb(r_y); + + // ... and get the squared challenge. + let challenge_sq_packed = transcript.squeeze().get_lower_128(); + let challenge_sq: C::Scalar = get_challenge_scalar(Challenge(challenge_sq_packed)); + + // There might be no square root, in which case we'll fork the + // transcript. + let challenge = challenge_sq.deterministic_sqrt(); + if let Some(challenge) = challenge { + break challenge; + } else { + // Try again, with slightly different randomness + l = (l + self.h).to_affine(); + l_randomness += &C::Scalar::one(); + } + }; + + // Challenge is unlikely to be zero. + let challenge_inv = challenge.invert().unwrap(); + let challenge_sq_inv = challenge_inv.square(); + let challenge_sq = challenge.square(); + + // Feed L and R into the real transcript + let (l_x, l_y) = l.get_xy().unwrap(); + let (r_x, r_y) = r.get_xy().unwrap(); + transcript.absorb(l_x); + transcript.absorb(l_y); + transcript.absorb(r_x); + transcript.absorb(r_y); + + // And obtain the challenge, even though we already have it, since + // squeezing affects the transcript. + { + let challenge_sq_packed = transcript.squeeze().get_lower_128(); + let challenge_sq_expected = get_challenge_scalar(Challenge(challenge_sq_packed)); + assert_eq!(challenge_sq, challenge_sq_expected); + } + + // Done with this round. + rounds.push((l, r)); + + // Collapse `a` and `b`. + // TODO: parallelize + for i in 0..half { + a[i] = (a[i] * &challenge) + &(a[i + half] * &challenge_inv); + b[i] = (b[i] * &challenge_inv) + &(b[i + half] * &challenge); + } + a.truncate(half); + b.truncate(half); + + // Collapse `G` + parallel_generator_collapse(&mut g, challenge, challenge_inv); + g.truncate(half); + + // Update randomness (the synthetic blinding factor at the end) + blind += &(l_randomness * &challenge_sq); + blind += &(r_randomness * &challenge_sq_inv); + } + + // We have fully colapsed `a`, `b`, `G` + assert_eq!(a.len(), 1); + let a = a[0]; + assert_eq!(b.len(), 1); + let b = b[0]; + assert_eq!(g.len(), 1); + let g = g[0]; + + // Random nonces for the zero-knowledge opening + let d = C::Scalar::random(); + let s = C::Scalar::random(); + + let delta = best_multiexp(&[d, d * &b, s], &[g, u, self.h]).to_affine(); + + let (delta_x, delta_y) = delta.get_xy().unwrap(); + + // Feed delta into the transcript + transcript.absorb(delta_x); + transcript.absorb(delta_y); + + // Obtain the challenge c. + let c_packed = transcript.squeeze().get_lower_128(); + let c: C::Scalar = get_challenge_scalar(Challenge(c_packed)); + + // Compute z1 and z2 as described in the Halo paper. + let z1 = a * &c + &d; + let z2 = c * &blind + &s; + + Ok(OpeningProof { + fork, + rounds, + delta, + z1, + z2, + }) + } + + /// Checks to see if an [`OpeningProof`] is valid given the current + /// `transcript`, and a point `x` that the polynomial commitment `p` opens + /// purportedly to the value `v`. + pub fn verify_proof>( + &self, + proof: &OpeningProof, + transcript: &mut H, + x: C::Scalar, + p: &C, + v: C::Scalar, + ) -> bool { + // Check for well-formedness + if proof.rounds.len() != self.k as usize { + return false; + } + + transcript.absorb(C::Base::from_u64(proof.fork as u64)); + + // Compute U + let u = { + let u_x = transcript.squeeze(); + // y^2 = x^3 + B + let u_y2 = u_x.square() * &u_x + &C::b(); + let u_y = u_y2.deterministic_sqrt(); + if u_y.is_none() { + return false; + } + let u_y = u_y.unwrap(); + + C::from_xy(u_x, u_y).unwrap() + }; + + let mut extra_scalars = Vec::with_capacity(proof.rounds.len() * 2 + 4 + self.n as usize); + let mut extra_bases = Vec::with_capacity(proof.rounds.len() * 2 + 4 + self.n as usize); + + // Data about the challenges from each of the rounds. + let mut challenges = Vec::with_capacity(proof.rounds.len()); + let mut challenges_inv = Vec::with_capacity(proof.rounds.len()); + let mut challenges_sq = Vec::with_capacity(proof.rounds.len()); + let mut allinv = Field::one(); + + for round in &proof.rounds { + // Feed L and R into the transcript. + let l = round.0.get_xy(); + let r = round.1.get_xy(); + if bool::from(l.is_none() | r.is_none()) { + return false; + } + let l = l.unwrap(); + let r = r.unwrap(); + transcript.absorb(l.0); + transcript.absorb(l.1); + transcript.absorb(r.0); + transcript.absorb(r.1); + let challenge_sq_packed = transcript.squeeze().get_lower_128(); + let challenge_sq: C::Scalar = get_challenge_scalar(Challenge(challenge_sq_packed)); + + let challenge = challenge_sq.deterministic_sqrt(); + if challenge.is_none() { + // We didn't sample a square. + return false; + } + let challenge = challenge.unwrap(); + + let challenge_inv = challenge.invert(); + if bool::from(challenge_inv.is_none()) { + // We sampled zero for some reason, unlikely to happen by + // chance. + return false; + } + let challenge_inv = challenge_inv.unwrap(); + allinv *= challenge_inv; + + let challenge_sq_inv = challenge_inv.square(); + + extra_scalars.push(challenge_sq); + extra_bases.push(round.0); + extra_scalars.push(challenge_sq_inv); + extra_bases.push(round.1); + + challenges.push(challenge); + challenges_inv.push(challenge_inv); + challenges_sq.push(challenge_sq); + } + + let delta = proof.delta.get_xy(); + if bool::from(delta.is_none()) { + return false; + } + let delta = delta.unwrap(); + + // Feed delta into the transcript + transcript.absorb(delta.0); + transcript.absorb(delta.1); + + // Get the challenge `c` + let c_packed = transcript.squeeze().get_lower_128(); + let c: C::Scalar = get_challenge_scalar(Challenge(c_packed)); + + // Check + // [c] P + [c * v] U + [c] sum(L_i * u_i^2) + [c] sum(R_i * u_i^-2) + delta - [z1] G - [z1 * b] U - [z2] H + // = 0 + + for scalar in &mut extra_scalars { + *scalar *= &c; + } + + let b = compute_b(x, &challenges, &challenges_inv); + + let neg_z1 = -proof.z1; + + // [c] P + extra_bases.push(*p); + extra_scalars.push(c); + + // [c * v] U - [z1 * b] U + extra_bases.push(u); + extra_scalars.push((c * &v) + &(neg_z1 * &b)); + + // delta + extra_bases.push(proof.delta); + extra_scalars.push(Field::one()); + + // - [z2] H + extra_bases.push(self.h); + extra_scalars.push(-proof.z2); + + // - [z1] G + extra_bases.extend(&self.g); + let mut s = compute_s(&challenges_sq, allinv); + // TODO: parallelize + for s in &mut s { + *s *= &neg_z1; + } + extra_scalars.extend(s); + + bool::from(best_multiexp(&extra_scalars, &extra_bases).is_zero()) + } +} + +fn compute_b(x: F, challenges: &[F], challenges_inv: &[F]) -> F { + assert!(!challenges.is_empty()); + assert_eq!(challenges.len(), challenges_inv.len()); + if challenges.len() == 1 { + *challenges_inv.last().unwrap() + *challenges.last().unwrap() * x + } else { + (*challenges_inv.last().unwrap() + *challenges.last().unwrap() * x) + * compute_b( + x.square(), + &challenges[0..(challenges.len() - 1)], + &challenges_inv[0..(challenges.len() - 1)], + ) + } +} + +// TODO: parallelize +fn compute_s(challenges_sq: &[F], allinv: F) -> Vec { + let lg_n = challenges_sq.len(); + let n = 1 << lg_n; + + let mut s = Vec::with_capacity(n); + s.push(allinv); + for i in 1..n { + let lg_i = (32 - 1 - (i as u32).leading_zeros()) as usize; + let k = 1 << lg_i; + let u_lg_i_sq = challenges_sq[(lg_n - 1) - lg_i]; + s.push(s[i - k] * u_lg_i_sq); + } + + s +} + +fn parallel_generator_collapse( + g: &mut [C], + challenge: C::Scalar, + challenge_inv: C::Scalar, +) { + let len = g.len() / 2; + let (mut g_lo, g_hi) = g.split_at_mut(len); + + parallelize(&mut g_lo, |g_lo, start| { + let g_hi = &g_hi[start..]; + let mut tmp = Vec::with_capacity(g_lo.len()); + for (g_lo, g_hi) in g_lo.iter().zip(g_hi.iter()) { + // TODO: could use multiexp + tmp.push(((*g_lo) * challenge_inv) + &((*g_hi) * challenge)); + } + C::Projective::batch_to_affine(&tmp, g_lo); + }); +} + +#[test] +fn test_commit_lagrange() { + const K: u32 = 6; + + use crate::arithmetic::{EpAffine, Fp, Fq}; + use crate::transcript::DummyHash; + let params = Params::::new::>(K); + + let a = (0..(1 << K)).map(|l| Fq::from(l)).collect::>(); + let mut b = a.clone(); + let mut alpha = Fq::ROOT_OF_UNITY; + for _ in K..Fq::S { + alpha = alpha.square(); + } + best_fft(&mut b, alpha, K); + + assert_eq!(params.commit(&a, alpha), params.commit_lagrange(&b, alpha)); +} + +#[test] +fn test_opening_proof() { + const K: u32 = 6; + + use crate::arithmetic::{eval_polynomial, EpAffine, Fp, Fq}; + use crate::transcript::DummyHash; + let params = Params::::new::>(K); + + let px = (0..(1 << K)) + .map(|l| Fq::from(l + 1) * Fq::ZETA) + .collect::>(); + let blind = Fq::random(); + + let p = params.commit(&px, blind).to_affine(); + + let mut transcript = DummyHash::init(Field::one()); + let (p_x, p_y) = p.get_xy().unwrap(); + transcript.absorb(p_x); + transcript.absorb(p_y); + let x_packed = transcript.squeeze().get_lower_128(); + let x: Fq = get_challenge_scalar(Challenge(x_packed)); + + // Evaluate the polynomial + let v = eval_polynomial(&px, x); + + transcript.absorb(Fp::from_bytes(&v.to_bytes()).unwrap()); // unlikely to fail since p ~ q + + loop { + let mut transcript_dup = transcript.clone(); + + let opening_proof = params.create_proof(&mut transcript, &px, blind, x); + if opening_proof.is_err() { + transcript = transcript_dup; + transcript.absorb(Field::one()); + } else { + let opening_proof = opening_proof.unwrap(); + assert!(params.verify_proof(&opening_proof, &mut transcript_dup, x, &p, v)); + break; + } + } +} diff --git a/src/transcript.rs b/src/transcript.rs new file mode 100644 index 0000000..8d4fde0 --- /dev/null +++ b/src/transcript.rs @@ -0,0 +1,45 @@ +//! This module contains utilities and traits for dealing with Fiat-Shamir +//! transcripts. + +use crate::arithmetic::Field; + +/// This is a generic interface for a sponge function that can be used for +/// Fiat-Shamir transformations. +pub trait Hasher: Clone + Send + Sync + 'static { + /// Initialize the sponge with some key. + fn init(key: F) -> Self; + /// Absorb a field element into the sponge. + fn absorb(&mut self, value: F); + /// Square a field element out of the sponge. + fn squeeze(&mut self) -> F; +} + +/// This is just a simple (and completely broken) hash function, standing in for +/// some algebraic hash function that we'll switch to later. +#[derive(Debug, Clone)] +pub struct DummyHash { + power: F, + state: F, +} + +impl Hasher for DummyHash { + fn init(value: F) -> Self { + DummyHash { + power: F::ZETA + F::one() + value, + state: F::ZETA, + } + } + fn absorb(&mut self, value: F) { + for _ in 0..10 { + self.state += value; + self.state *= self.power; + self.power += self.power.square(); + self.state += self.power; + } + } + fn squeeze(&mut self) -> F { + let tmp = self.state; + self.absorb(tmp); + tmp + } +}