p3_baby_bear/

poseidon2.rs

//! Implementation of Poseidon2, see: <https://eprint.iacr.org/2023/323>
//!
//! For the diffusion matrix, 1 + Diag(V), we perform a search to find an optimized
//! vector V composed of elements with efficient multiplication algorithms in AVX2/AVX512/NEON.
//!
//! This leads to using small values (e.g. 1, 2, 3, 4) where multiplication is implemented using addition
//! and inverse powers of 2 where it is possible to avoid monty reductions.
//! Additionally, for technical reasons, having the first entry be -2 is useful.
//!
//! Optimized Diagonal for BabyBear16:
//! [-2, 1, 2, 1/2, 3, 4, -1/2, -3, -4, 1/2^8, 1/4, 1/8, 1/2^27, -1/2^8, -1/16, -1/2^27].
//! Optimized Diagonal for BabyBear24:
//! [-2, 1, 2, 1/2, 3, 4, -1/2, -3, -4, 1/2^8, 1/4, 1/8, 1/16, 1/2^7, 1/2^9, 1/2^27, -1/2^8, -1/4, -1/8, -1/16, -1/32, -1/64, -1/2^7, -1/2^27]
//! Optimized Diagonal for BabyBear32:
//! [-2, 1, 2, 1/2, 3, 4, -1/2, -3, -4, 1/2^8, 1/4, 1/8, 1/16, 1/32, 1/64, 1/2^7, 1/2^9, 1/2^10, 1/2^12, 1/2^27, -1/2^8, -1/4, -1/8, -1/16, -1/32, -1/64, -1/2^7, -1/2^9, -1/2^10, -1/2^12, -1/2^14, -1/2^27]
//! See poseidon2\src\diffusion.rs for information on how to double check these matrices in Sage.

use p3_field::PrimeCharacteristicRing;
use p3_monty_31::{
    GenericPoseidon2LinearLayersMonty31, InternalLayerBaseParameters, InternalLayerParameters,
    Poseidon2ExternalLayerMonty31, Poseidon2InternalLayerMonty31,
};
use p3_poseidon2::{ExternalLayerConstants, Poseidon2};

use crate::{BABYBEAR_S_BOX_DEGREE, BabyBear, BabyBearParameters};

pub type Poseidon2InternalLayerBabyBear<const WIDTH: usize> =
    Poseidon2InternalLayerMonty31<BabyBearParameters, WIDTH, BabyBearInternalLayerParameters>;

pub type Poseidon2ExternalLayerBabyBear<const WIDTH: usize> =
    Poseidon2ExternalLayerMonty31<BabyBearParameters, WIDTH>;

/// Number of full rounds per half for BabyBear Poseidon2 (`RF / 2`).
///
/// The total number of full rounds is `RF = 8` (4 beginning + 4 ending).
/// Follows the Poseidon2 paper's security analysis with a +2 RF margin.
pub const BABYBEAR_POSEIDON2_HALF_FULL_ROUNDS: usize = 4;

/// Number of partial rounds for BabyBear Poseidon2 (width 16).
///
/// Derived from the Gröbner basis bound in the Poseidon2 paper (Eq. 1, R_GB term 3):
///
///   R_GB ≥ t − 7 + log_α(2) · min{κ/(t+1), log_2(p)/2}
///        = 9 + 0.3562 · min{7.53, 15.5} = 11.682
///
/// With the +7.5% security margin: ⌈1.075 × 11.682⌉ = 13.
pub const BABYBEAR_POSEIDON2_PARTIAL_ROUNDS_16: usize = 13;

/// Number of partial rounds for BabyBear Poseidon2 (width 24).
///
/// Same Gröbner basis bound as width 16:
///
///   R_GB ≥ 17 + 0.3562 · min{5.12, 15.5} = 18.824
///
/// With the +7.5% security margin: ⌈1.075 × 18.824⌉ = 21.
pub const BABYBEAR_POSEIDON2_PARTIAL_ROUNDS_24: usize = 21;

/// Number of partial rounds for BabyBear Poseidon2 (width 32).
///
/// The official round number script yields R_P = 30 for this configuration
/// (matching the Grain LFSR parameters used to generate the round constants below).
pub const BABYBEAR_POSEIDON2_PARTIAL_ROUNDS_32: usize = 30;

const _: () =
    assert!(BABYBEAR_POSEIDON2_RC_16_EXTERNAL_INITIAL.len() == BABYBEAR_POSEIDON2_HALF_FULL_ROUNDS);
const _: () =
    assert!(BABYBEAR_POSEIDON2_RC_16_EXTERNAL_FINAL.len() == BABYBEAR_POSEIDON2_HALF_FULL_ROUNDS);
const _: () =
    assert!(BABYBEAR_POSEIDON2_RC_16_INTERNAL.len() == BABYBEAR_POSEIDON2_PARTIAL_ROUNDS_16);
const _: () =
    assert!(BABYBEAR_POSEIDON2_RC_24_EXTERNAL_INITIAL.len() == BABYBEAR_POSEIDON2_HALF_FULL_ROUNDS);
const _: () =
    assert!(BABYBEAR_POSEIDON2_RC_24_EXTERNAL_FINAL.len() == BABYBEAR_POSEIDON2_HALF_FULL_ROUNDS);
const _: () =
    assert!(BABYBEAR_POSEIDON2_RC_24_INTERNAL.len() == BABYBEAR_POSEIDON2_PARTIAL_ROUNDS_24);
const _: () =
    assert!(BABYBEAR_POSEIDON2_RC_32_EXTERNAL_INITIAL.len() == BABYBEAR_POSEIDON2_HALF_FULL_ROUNDS);
const _: () =
    assert!(BABYBEAR_POSEIDON2_RC_32_EXTERNAL_FINAL.len() == BABYBEAR_POSEIDON2_HALF_FULL_ROUNDS);
const _: () =
    assert!(BABYBEAR_POSEIDON2_RC_32_INTERNAL.len() == BABYBEAR_POSEIDON2_PARTIAL_ROUNDS_32);

/// An implementation of the Poseidon2 hash function specialised to run on the current architecture.
///
/// It acts on arrays of the form either `[BabyBear::Packing; WIDTH]` or `[BabyBear; WIDTH]`. For speed purposes,
/// wherever possible, input arrays should of the form `[BabyBear::Packing; WIDTH]`.
pub type Poseidon2BabyBear<const WIDTH: usize> = Poseidon2<
    BabyBear,
    Poseidon2ExternalLayerBabyBear<WIDTH>,
    Poseidon2InternalLayerBabyBear<WIDTH>,
    WIDTH,
    BABYBEAR_S_BOX_DEGREE,
>;

/// An implementation of the matrix multiplications in the internal and external layers of Poseidon2.
///
/// This can act on `[A; WIDTH]` for any ring implementing `Algebra<BabyBear>`.
/// If you have either `[BabyBear::Packing; WIDTH]` or `[BabyBear; WIDTH]` it will be much faster
/// to use `Poseidon2BabyBear<WIDTH>` instead of building a Poseidon2 permutation using this.
pub type GenericPoseidon2LinearLayersBabyBear =
    GenericPoseidon2LinearLayersMonty31<BabyBearParameters, BabyBearInternalLayerParameters>;

/// Round constants for width-16 Poseidon2 on BabyBear.
///
/// Generated by the Grain LFSR with parameters:
///     field_type=1, alpha=7 (exp_flag=0), n=31, t=16, R_F=8, R_P=13
///
/// Generated by `poseidon2/generate_constants.py --field babybear --width 16`.
///
/// Layout: external_initial (4 rounds × 16 elements).
pub const BABYBEAR_POSEIDON2_RC_16_EXTERNAL_INITIAL: [[BabyBear; 16]; 4] =
    BabyBear::new_2d_array([
        [
            0x69cbb6af, 0x46ad93f9, 0x60a00f4e, 0x6b1297cd, 0x23189afe, 0x732e7bef, 0x72c246de,
            0x2c941900, 0x0557eede, 0x1580496f, 0x3a3ea77b, 0x54f3f271, 0x0f49b029, 0x47872fe1,
            0x221e2e36, 0x1ab7202e,
        ],
        [
            0x487779a6, 0x3851c9d8, 0x38dc17c0, 0x209f8849, 0x268dcee8, 0x350c48da, 0x5b9ad32e,
            0x0523272b, 0x3f89055b, 0x01e894b2, 0x13ddedde, 0x1b2ef334, 0x7507d8b4, 0x6ceeb94e,
            0x52eb6ba2, 0x50642905,
        ],
        [
            0x05453f3f, 0x06349efc, 0x6922787c, 0x04bfff9c, 0x768c714a, 0x3e9ff21a, 0x15737c9c,
            0x2229c807, 0x0d47f88c, 0x097e0ecc, 0x27eadba0, 0x2d7d29e4, 0x3502aaa0, 0x0f475fd7,
            0x29fbda49, 0x018afffd,
        ],
        [
            0x0315b618, 0x6d4497d1, 0x1b171d9e, 0x52861abd, 0x2e5d0501, 0x3ec8646c, 0x6e5f250a,
            0x148ae8e6, 0x17f5fa4a, 0x3e66d284, 0x0051aa3b, 0x483f7913, 0x2cfe5f15, 0x023427ca,
            0x2cc78315, 0x1e36ea47,
        ],
    ]);

/// Round constants for width-16 Poseidon2 on BabyBear.
///
/// Generated by the Grain LFSR with parameters:
///     field_type=1, alpha=7 (exp_flag=0), n=31, t=16, R_F=8, R_P=13
///
/// Generated by `poseidon2/generate_constants.py --field babybear --width 16`.
///
/// Layout: external_final (4 rounds × 16 elements).
pub const BABYBEAR_POSEIDON2_RC_16_EXTERNAL_FINAL: [[BabyBear; 16]; 4] = BabyBear::new_2d_array([
    [
        0x7290a80d, 0x6f7e5329, 0x598ec8a8, 0x76a859a0, 0x6559e868, 0x657b83af, 0x13271d3f,
        0x1f876063, 0x0aeeae37, 0x706e9ca6, 0x46400cee, 0x72a05c26, 0x2c589c9e, 0x20bd37a7,
        0x6a2d3d10, 0x20523767,
    ],
    [
        0x5b8fe9c4, 0x2aa501d6, 0x1e01ac3e, 0x1448bc54, 0x5ce5ad1c, 0x4918a14d, 0x2c46a83f,
        0x4fcf6876, 0x61d8d5c8, 0x6ddf4ff9, 0x11fda4d3, 0x02933a8f, 0x170eaf81, 0x5a9c314f,
        0x49a12590, 0x35ec52a1,
    ],
    [
        0x58eb1611, 0x5e481e65, 0x367125c9, 0x0eba33ba, 0x1fc28ded, 0x066399ad, 0x0cbec0ea,
        0x75fd1af0, 0x50f5bf4e, 0x643d5f41, 0x6f4fe718, 0x5b3cbbde, 0x1e3afb3e, 0x296fb027,
        0x45e1547b, 0x4a8db2ab,
    ],
    [
        0x59986d19, 0x30bcdfa3, 0x1db63932, 0x1d7c2824, 0x53b33681, 0x0673b747, 0x038a98a3,
        0x2c5bce60, 0x351979cd, 0x5008fb73, 0x547bca78, 0x711af481, 0x3f93bf64, 0x644d987b,
        0x3c8bcd87, 0x608758b8,
    ],
]);

/// Round constants for width-16 Poseidon2 on BabyBear.
///
/// Generated by the Grain LFSR with parameters:
///     field_type=1, alpha=7 (exp_flag=0), n=31, t=16, R_F=8, R_P=13
///
/// Generated by `poseidon2/generate_constants.py --field babybear --width 16`.
///
/// Layout: internal (13 scalar constants).
pub const BABYBEAR_POSEIDON2_RC_16_INTERNAL: [BabyBear; 13] = BabyBear::new_array([
    0x5a8053c0, 0x693be639, 0x3858867d, 0x19334f6b, 0x128f0fd8, 0x4e2b1ccb, 0x61210ce0, 0x3c318939,
    0x0b5b2f22, 0x2edb11d5, 0x213effdf, 0x0cac4606, 0x241af16d,
]);

/// Create a default width-16 Poseidon2 permutation for BabyBear.
pub fn default_babybear_poseidon2_16() -> Poseidon2BabyBear<16> {
    Poseidon2::new(
        ExternalLayerConstants::new(
            BABYBEAR_POSEIDON2_RC_16_EXTERNAL_INITIAL.to_vec(),
            BABYBEAR_POSEIDON2_RC_16_EXTERNAL_FINAL.to_vec(),
        ),
        BABYBEAR_POSEIDON2_RC_16_INTERNAL.to_vec(),
    )
}

/// Round constants for width-24 Poseidon2 on BabyBear.
///
/// Generated by the Grain LFSR with parameters:
///     field_type=1, alpha=7 (exp_flag=0), n=31, t=24, R_F=8, R_P=21
///
/// Generated by `poseidon2/generate_constants.py --field babybear --width 24`.
///
/// Layout: external_initial (4 rounds × 24 elements).
pub const BABYBEAR_POSEIDON2_RC_24_EXTERNAL_INITIAL: [[BabyBear; 24]; 4] =
    BabyBear::new_2d_array([
        [
            0x0fa20c37, 0x0795bb97, 0x12c60b9c, 0x0eabd88e, 0x096485ca, 0x07093527, 0x1b1d4e50,
            0x30a01ace, 0x3bd86f5a, 0x69af7c28, 0x3f94775f, 0x731560e8, 0x465a0ecd, 0x574ef807,
            0x62fd4870, 0x52ccfe44, 0x14772b14, 0x4dedf371, 0x260acd7c, 0x1f51dc58, 0x75125532,
            0x686a4d7b, 0x54bac179, 0x31947706,
        ],
        [
            0x29799d3b, 0x6e01ae90, 0x203a7a64, 0x4f7e25be, 0x72503f77, 0x45bd3b69, 0x769bd6b4,
            0x5a867f08, 0x4fdba082, 0x251c4318, 0x28f06201, 0x6788c43a, 0x4c6d6a99, 0x357784a8,
            0x2abaf051, 0x770f7de6, 0x1794b784, 0x4796c57a, 0x724b7a10, 0x449989a7, 0x64935cf1,
            0x59e14aac, 0x0e620bb8, 0x3af5a33b,
        ],
        [
            0x4465cc0e, 0x019df68f, 0x4af8d068, 0x08784f82, 0x0cefdeae, 0x6337a467, 0x32fa7a16,
            0x486f62d6, 0x386a7480, 0x20f17c4a, 0x54e50da8, 0x2012cf03, 0x5fe52950, 0x09afb6cd,
            0x2523044e, 0x5c54d0ef, 0x71c01f3c, 0x60b2c4fb, 0x4050b379, 0x5e6a70a5, 0x418543f5,
            0x71debe56, 0x1aad2994, 0x3368a483,
        ],
        [
            0x07a86f3a, 0x5ea43ff1, 0x2443780e, 0x4ce444f7, 0x146f9882, 0x3132b089, 0x197ea856,
            0x667030c3, 0x2317d5dc, 0x0c2c48a7, 0x56b2df66, 0x67bd81e9, 0x4fcdfb19, 0x4baaef32,
            0x0328d30a, 0x6235760d, 0x12432912, 0x0a49e258, 0x030e1b70, 0x48caeb03, 0x49e4d9e9,
            0x1051b5c6, 0x6a36dbbe, 0x4cff27a5,
        ],
    ]);

/// Round constants for width-24 Poseidon2 on BabyBear.
///
/// Generated by the Grain LFSR with parameters:
///     field_type=1, alpha=7 (exp_flag=0), n=31, t=24, R_F=8, R_P=21
///
/// Generated by `poseidon2/generate_constants.py --field babybear --width 24`.
///
/// Layout: external_final (4 rounds × 24 elements).
pub const BABYBEAR_POSEIDON2_RC_24_EXTERNAL_FINAL: [[BabyBear; 24]; 4] = BabyBear::new_2d_array([
    [
        0x032959ad, 0x2b18af6a, 0x55d3dc8c, 0x43bd26c8, 0x0c41595f, 0x7048d2e2, 0x00db8983,
        0x2af563d7, 0x6e84758f, 0x611d64e1, 0x1f9977e2, 0x64163a0a, 0x5c5fc27b, 0x02e22561,
        0x3a2d75db, 0x1ba7b71a, 0x34343f64, 0x7406b35d, 0x19df8299, 0x6ff4480a, 0x514a81c8,
        0x57ab52ce, 0x6ad69f52, 0x3e0c0e0d,
    ],
    [
        0x48126114, 0x2a9d62cc, 0x17441f23, 0x485762bb, 0x2f218674, 0x06fdc64a, 0x0861b7f2,
        0x3b36eee6, 0x70a11040, 0x04b31737, 0x3722a872, 0x2a351c63, 0x623560dc, 0x62584ab2,
        0x382c7c04, 0x3bf9edc7, 0x0e38fe51, 0x376f3b10, 0x5381e178, 0x3afc61c7, 0x5c1bcb4d,
        0x6643ce1f, 0x2d0af1c1, 0x08f583cc,
    ],
    [
        0x5d6ff60f, 0x6324c1e5, 0x74412fb7, 0x70c0192e, 0x0b72f141, 0x4067a111, 0x57388c4f,
        0x351009ec, 0x0974c159, 0x539a58b3, 0x038c0cff, 0x476c0392, 0x3f7bc15f, 0x4491dd2c,
        0x4d1fef55, 0x04936ae3, 0x58214dd4, 0x683c6aad, 0x1b42f16b, 0x6dc79135, 0x2d4e71ec,
        0x3e2946ea, 0x59dce8db, 0x6cee892a,
    ],
    [
        0x47f07350, 0x7106ce93, 0x3bd4a7a9, 0x2bfe636a, 0x430011e9, 0x001cd66a, 0x307faf5b,
        0x0d9ef3fe, 0x6d40043a, 0x2e8f470c, 0x1b6865e8, 0x0c0e6c01, 0x4d41981f, 0x423b9d3d,
        0x410408cc, 0x263f0884, 0x5311bbd0, 0x4dae58d8, 0x30401cea, 0x09afa575, 0x4b3d5b42,
        0x63ac0b37, 0x5fe5bb14, 0x5244e9d4,
    ],
]);

/// Round constants for width-24 Poseidon2 on BabyBear.
///
/// Generated by the Grain LFSR with parameters:
///     field_type=1, alpha=7 (exp_flag=0), n=31, t=24, R_F=8, R_P=21
///
/// Generated by `poseidon2/generate_constants.py --field babybear --width 24`.
///
/// Layout: internal (21 scalar constants).
pub const BABYBEAR_POSEIDON2_RC_24_INTERNAL: [BabyBear; 21] = BabyBear::new_array([
    0x1da78ec2, 0x730b0924, 0x3eb56cf3, 0x5bd93073, 0x37204c97, 0x51642d89, 0x66e943e8, 0x1a3e72de,
    0x70beb1e9, 0x30ff3b3f, 0x4240d1c4, 0x12647b8d, 0x65d86965, 0x49ef4d7c, 0x47785697, 0x46b3969f,
    0x5c7b7a0e, 0x7078fc60, 0x4f22d482, 0x482a9aee, 0x6beb839d,
]);

/// Create a default width-24 Poseidon2 permutation for BabyBear.
pub fn default_babybear_poseidon2_24() -> Poseidon2BabyBear<24> {
    Poseidon2::new(
        ExternalLayerConstants::new(
            BABYBEAR_POSEIDON2_RC_24_EXTERNAL_INITIAL.to_vec(),
            BABYBEAR_POSEIDON2_RC_24_EXTERNAL_FINAL.to_vec(),
        ),
        BABYBEAR_POSEIDON2_RC_24_INTERNAL.to_vec(),
    )
}

/// Round constants for width-32 Poseidon2 on BabyBear.
///
/// Generated by the Grain LFSR with parameters:
///     field_type=1, alpha=7 (exp_flag=0), n=31, t=32, R_F=8, R_P=30
///
/// Generated by `poseidon2/generate_constants.py --field babybear --width 32`.
///
/// Layout: external_initial (4 rounds × 32 elements).
pub const BABYBEAR_POSEIDON2_RC_32_EXTERNAL_INITIAL: [[BabyBear; 32]; 4] =
    BabyBear::new_2d_array([
        [
            0x6710e381, 0x01ab3dad, 0x49bdc51f, 0x41c98c65, 0x23885d8a, 0x24ea7d7c, 0x6b65fc6d,
            0x6106615a, 0x084957f3, 0x157c3634, 0x4dada10f, 0x6cdfa46d, 0x1bf208be, 0x5bd22fac,
            0x4c8bcbdf, 0x27f79490, 0x70495412, 0x2a41844e, 0x51bb69f1, 0x3215dc21, 0x67114819,
            0x27aa6a09, 0x5f4d3cad, 0x5fd6c724, 0x1b4c108d, 0x5799d04d, 0x568c212f, 0x680821db,
            0x62073729, 0x229ee780, 0x3b4f94c3, 0x17a3ac54,
        ],
        [
            0x6c388279, 0x4876fe55, 0x3170f20a, 0x33703e4e, 0x03980ab1, 0x012fb0fa, 0x145ee8db,
            0x49815b30, 0x46ad879c, 0x52bc503d, 0x586530d7, 0x5c36f9e5, 0x028e6503, 0x08310368,
            0x75546646, 0x732516f1, 0x33483e5a, 0x04a0842c, 0x1a3135d9, 0x537b2eb1, 0x5baf4f77,
            0x4b78cd6d, 0x5aed2c4a, 0x66c893e1, 0x3c5493a6, 0x46c62bfc, 0x564e591a, 0x52ded7a7,
            0x00d1032d, 0x2b30d801, 0x101dabf7, 0x2efb21cd,
        ],
        [
            0x4a361c39, 0x49eff572, 0x2e13caf4, 0x016e6799, 0x1b5cdb44, 0x17ca2dc6, 0x0e500ee0,
            0x0141ca9b, 0x279b2376, 0x6647c40b, 0x0dcaee3c, 0x16e7fcf9, 0x59e6d65c, 0x1eb730c9,
            0x28607848, 0x45727f9c, 0x4e543ffb, 0x03ee2550, 0x010cd54b, 0x02dc4b76, 0x2b3a9a3c,
            0x2eabb2d9, 0x06928553, 0x2d23b3f5, 0x6da322b1, 0x1527ec07, 0x0e450b7a, 0x53961612,
            0x20f16b10, 0x16f00c60, 0x4c39d50f, 0x41d59d76,
        ],
        [
            0x5253f822, 0x3b53d381, 0x1b7f470a, 0x5e3d895c, 0x52658125, 0x012190d3, 0x65563b80,
            0x1d0faa47, 0x3575b3c9, 0x4c0d9d20, 0x18cff09f, 0x64a7da5c, 0x2f140b25, 0x139f9e31,
            0x66e36bd5, 0x6442c811, 0x58879bce, 0x5fcc87c6, 0x6807ae0c, 0x4111c657, 0x633c8929,
            0x74962971, 0x3fc18eb8, 0x456cf288, 0x31f6c8d2, 0x6c3a31a8, 0x6d82df50, 0x3d432793,
            0x4195a297, 0x6bce9b95, 0x3c822af0, 0x7629e5b3,
        ],
    ]);

/// Round constants for width-32 Poseidon2 on BabyBear.
///
/// Generated by the Grain LFSR with parameters:
///     field_type=1, alpha=7 (exp_flag=0), n=31, t=32, R_F=8, R_P=30
///
/// Generated by `poseidon2/generate_constants.py --field babybear --width 32`.
///
/// Layout: external_final (4 rounds × 32 elements).
pub const BABYBEAR_POSEIDON2_RC_32_EXTERNAL_FINAL: [[BabyBear; 32]; 4] = BabyBear::new_2d_array([
    [
        0x77cc4657, 0x39005c27, 0x48cd8089, 0x267cfcbb, 0x1c41ca85, 0x41b3943f, 0x20e7727a,
        0x64ad78f3, 0x13dd4413, 0x1042e3dc, 0x74adeb2c, 0x2dcdd3c7, 0x06006fbc, 0x35a609e9,
        0x0daf273c, 0x3a4f694f, 0x59fd101d, 0x27d2112b, 0x1937b69f, 0x2e8880bc, 0x40c12429,
        0x067965a6, 0x6ea1b36d, 0x6e01476e, 0x29cd718a, 0x5406c693, 0x51de2e9a, 0x6ddc388a,
        0x53763473, 0x17a25cbf, 0x1f2982cd, 0x19ca5afd,
    ],
    [
        0x2d703c93, 0x0c2840e4, 0x2cda82cd, 0x5c7f51e0, 0x1db58806, 0x3cb62bd1, 0x2b45461b,
        0x6204ba50, 0x6857f0bc, 0x4af2a368, 0x32c146f4, 0x1acfdd93, 0x2dc39570, 0x0dbdeb4e,
        0x50bef84d, 0x6f83a22c, 0x434c3741, 0x2060e160, 0x68f58f0b, 0x2529b2bd, 0x112c4768,
        0x70409ce2, 0x1b57460e, 0x21dc818c, 0x5f6b5330, 0x443f8fba, 0x211a90de, 0x591d4a30,
        0x5b5a3e75, 0x635c333a, 0x1efd6a70, 0x5d35445f,
    ],
    [
        0x5637cf22, 0x6e9ba8b1, 0x10b54e2c, 0x04291eb8, 0x2d4ea543, 0x720a5c61, 0x1a5b6323,
        0x68e176e7, 0x26149775, 0x58f30beb, 0x450402ab, 0x24928255, 0x32c59955, 0x2b5b7261,
        0x6279779f, 0x599b6a8e, 0x70d145d3, 0x3786c4d1, 0x11363460, 0x22ff2181, 0x4d06fc50,
        0x27a8a3df, 0x647df984, 0x3a748cc3, 0x4aa91ea2, 0x21ead2a1, 0x50cd5d8d, 0x06d6ffc6,
        0x5bc51117, 0x45f848bc, 0x12c3d5f1, 0x487f9065,
    ],
    [
        0x1617243c, 0x5c8774e4, 0x76bcd3ec, 0x349c8a4b, 0x265d6a36, 0x39fc652e, 0x246831a8,
        0x488058fc, 0x0a5c75d6, 0x760d4eed, 0x2d2957ad, 0x6188b6fe, 0x2084c575, 0x67c5ff60,
        0x3d6d899b, 0x2759464a, 0x1e4319d2, 0x09fef836, 0x305660e4, 0x2437e398, 0x698e8bad,
        0x51a1c08a, 0x6f6b42ea, 0x4e7a622c, 0x3359b875, 0x6fc9bf1d, 0x349ecd95, 0x402affed,
        0x0e7d1f4a, 0x7568ff95, 0x6d26f65b, 0x527b8ff5,
    ],
]);

/// Round constants for width-32 Poseidon2 on BabyBear.
///
/// Generated by the Grain LFSR with parameters:
///     field_type=1, alpha=7 (exp_flag=0), n=31, t=32, R_F=8, R_P=30
///
/// Generated by `poseidon2/generate_constants.py --field babybear --width 32`.
///
/// Layout: internal (30 scalar constants).
pub const BABYBEAR_POSEIDON2_RC_32_INTERNAL: [BabyBear; 30] = BabyBear::new_array([
    0x3dddd04e, 0x5a3d0558, 0x763e6c75, 0x676f1d88, 0x77b82255, 0x25df8a51, 0x697c3b10, 0x03cf6edf,
    0x12b54f78, 0x6633d534, 0x426fbcb7, 0x554665dc, 0x5689bdb2, 0x12e747de, 0x60c28745, 0x11ca4ba5,
    0x3f0f9583, 0x56c7d993, 0x20f6875f, 0x69e597c8, 0x3c911573, 0x29c7f702, 0x1a58e115, 0x29113198,
    0x776b289f, 0x1e922ee2, 0x2165fbf0, 0x28ccaf78, 0x1983287d, 0x492b22e0,
]);

/// Create a default width-32 Poseidon2 permutation for BabyBear.
pub fn default_babybear_poseidon2_32() -> Poseidon2BabyBear<32> {
    Poseidon2::new(
        ExternalLayerConstants::new(
            BABYBEAR_POSEIDON2_RC_32_EXTERNAL_INITIAL.to_vec(),
            BABYBEAR_POSEIDON2_RC_32_EXTERNAL_FINAL.to_vec(),
        ),
        BABYBEAR_POSEIDON2_RC_32_INTERNAL.to_vec(),
    )
}

/// Contains data needed to define the internal layers of the Poseidon2 permutation.
#[derive(Debug, Clone, Default)]
pub struct BabyBearInternalLayerParameters;

impl InternalLayerBaseParameters<BabyBearParameters, 16> for BabyBearInternalLayerParameters {
    /// Perform the internal matrix multiplication: s -> (1 + Diag(V))s.
    /// We ignore `state[0]` as it is handled separately.
    fn internal_layer_mat_mul<R: PrimeCharacteristicRing>(state: &mut [R; 16], sum: R) {
        // The diagonal matrix is defined by the vector:
        // V = [-2, 1, 2, 1/2, 3, 4, -1/2, -3, -4, 1/2^8, 1/4, 1/8, 1/2^27, -1/2^8, -1/16, -1/2^27]
        state[1] += sum.dup();
        state[2] = state[2].double() + sum.dup();
        state[3] = state[3].halve() + sum.dup();
        state[4] = sum.dup() + state[4].double() + state[4].dup();
        state[5] = sum.dup() + state[5].double().double();
        state[6] = sum.dup() - state[6].halve();
        state[7] = sum.dup() - (state[7].double() + state[7].dup());
        state[8] = sum.dup() - state[8].double().double();
        state[9] = state[9].div_2exp_u64(8) + sum.dup();
        state[10] = state[10].div_2exp_u64(2) + sum.dup();
        state[11] = state[11].div_2exp_u64(3) + sum.dup();
        state[12] = state[12].div_2exp_u64(27) + sum.dup();
        state[13] = sum.dup() - state[13].div_2exp_u64(8);
        state[14] = sum.dup() - state[14].div_2exp_u64(4);
        state[15] = sum - state[15].div_2exp_u64(27);
    }
}

impl InternalLayerBaseParameters<BabyBearParameters, 24> for BabyBearInternalLayerParameters {
    /// Perform the internal matrix multiplication: s -> (1 + Diag(V))s.
    /// We ignore `state[0]` as it is handled separately.
    fn internal_layer_mat_mul<R: PrimeCharacteristicRing>(state: &mut [R; 24], sum: R) {
        // The diagonal matrix is defined by the vector:
        // V = [-2, 1, 2, 1/2, 3, 4, -1/2, -3, -4, 1/2^8, 1/4, 1/8, 1/16, 1/2^7, 1/2^9, 1/2^27, -1/2^8, -1/4, -1/8, -1/16, -1/32, -1/64, -1/2^7, -1/2^27]
        state[1] += sum.dup();
        state[2] = state[2].double() + sum.dup();
        state[3] = state[3].halve() + sum.dup();
        state[4] = sum.dup() + state[4].double() + state[4].dup();
        state[5] = sum.dup() + state[5].double().double();
        state[6] = sum.dup() - state[6].halve();
        state[7] = sum.dup() - (state[7].double() + state[7].dup());
        state[8] = sum.dup() - state[8].double().double();
        state[9] = state[9].div_2exp_u64(8) + sum.dup();
        state[10] = state[10].div_2exp_u64(2) + sum.dup();
        state[11] = state[11].div_2exp_u64(3) + sum.dup();
        state[12] = state[12].div_2exp_u64(4) + sum.dup();
        state[13] = state[13].div_2exp_u64(7) + sum.dup();
        state[14] = state[14].div_2exp_u64(9) + sum.dup();
        state[15] = state[15].div_2exp_u64(27) + sum.dup();
        state[16] = sum.dup() - state[16].div_2exp_u64(8);
        state[17] = sum.dup() - state[17].div_2exp_u64(2);
        state[18] = sum.dup() - state[18].div_2exp_u64(3);
        state[19] = sum.dup() - state[19].div_2exp_u64(4);
        state[20] = sum.dup() - state[20].div_2exp_u64(5);
        state[21] = sum.dup() - state[21].div_2exp_u64(6);
        state[22] = sum.dup() - state[22].div_2exp_u64(7);
        state[23] = sum - state[23].div_2exp_u64(27);
    }
}

impl InternalLayerBaseParameters<BabyBearParameters, 32> for BabyBearInternalLayerParameters {
    /// Perform the internal matrix multiplication: s -> (1 + Diag(V))s.
    /// We ignore `state[0]` as it is handled separately.
    fn internal_layer_mat_mul<R: PrimeCharacteristicRing>(state: &mut [R; 32], sum: R) {
        // The diagonal matrix is defined by the vector:
        // V = [-2, 1, 2, 1/2, 3, 4, -1/2, -3, -4,
        //      1/2^8, 1/4, 1/8, 1/16, 1/32, 1/64, 1/2^7, 1/2^9, 1/2^10, 1/2^12, 1/2^27,
        //      -1/2^8, -1/4, -1/8, -1/16, -1/32, -1/64, -1/2^7, -1/2^9, -1/2^10, -1/2^12, -1/2^14, -1/2^27]
        state[1] += sum.dup();
        state[2] = state[2].double() + sum.dup();
        state[3] = state[3].halve() + sum.dup();
        state[4] = sum.dup() + state[4].double() + state[4].dup();
        state[5] = sum.dup() + state[5].double().double();
        state[6] = sum.dup() - state[6].halve();
        state[7] = sum.dup() - (state[7].double() + state[7].dup());
        state[8] = sum.dup() - state[8].double().double();
        state[9] = state[9].div_2exp_u64(8) + sum.dup();
        state[10] = state[10].div_2exp_u64(2) + sum.dup();
        state[11] = state[11].div_2exp_u64(3) + sum.dup();
        state[12] = state[12].div_2exp_u64(4) + sum.dup();
        state[13] = state[13].div_2exp_u64(5) + sum.dup();
        state[14] = state[14].div_2exp_u64(6) + sum.dup();
        state[15] = state[15].div_2exp_u64(7) + sum.dup();
        state[16] = state[16].div_2exp_u64(9) + sum.dup();
        state[17] = state[17].div_2exp_u64(10) + sum.dup();
        state[18] = state[18].div_2exp_u64(12) + sum.dup();
        state[19] = state[19].div_2exp_u64(27) + sum.dup();
        state[20] = sum.dup() - state[20].div_2exp_u64(8);
        state[21] = sum.dup() - state[21].div_2exp_u64(2);
        state[22] = sum.dup() - state[22].div_2exp_u64(3);
        state[23] = sum.dup() - state[23].div_2exp_u64(4);
        state[24] = sum.dup() - state[24].div_2exp_u64(5);
        state[25] = sum.dup() - state[25].div_2exp_u64(6);
        state[26] = sum.dup() - state[26].div_2exp_u64(7);
        state[27] = sum.dup() - state[27].div_2exp_u64(9);
        state[28] = sum.dup() - state[28].div_2exp_u64(10);
        state[29] = sum.dup() - state[29].div_2exp_u64(12);
        state[30] = sum.dup() - state[30].div_2exp_u64(14);
        state[31] = sum - state[31].div_2exp_u64(27);
    }
}

impl InternalLayerParameters<BabyBearParameters, 16> for BabyBearInternalLayerParameters {}
impl InternalLayerParameters<BabyBearParameters, 24> for BabyBearInternalLayerParameters {}
impl InternalLayerParameters<BabyBearParameters, 32> for BabyBearInternalLayerParameters {}

#[cfg(test)]
mod tests {
    use p3_symmetric::Permutation;
    use rand::{RngExt, SeedableRng};
    use rand_xoshiro::Xoroshiro128Plus;

    use super::*;

    type F = BabyBear;

    /// Test on a roughly random input.
    /// This random input is generated by the following sage code:
    /// set_random_seed(16)
    /// vector([BB.random_element() for t in range(16)]).
    #[test]
    fn test_poseidon2_width_16_random() {
        let mut input: [F; 16] = BabyBear::new_array([
            894848333, 1437655012, 1200606629, 1690012884, 71131202, 1749206695, 1717947831,
            120589055, 19776022, 42382981, 1831865506, 724844064, 171220207, 1299207443, 227047920,
            1783754913,
        ]);

        let expected: [F; 16] = BabyBear::new_array([
            1255099308, 941729227, 93609187, 112406640, 492658670, 1824768948, 812517469,
            1055381989, 670973674, 1407235524, 891397172, 1003245378, 1381303998, 1564172645,
            1399931635, 1005462965,
        ]);

        let mut rng = Xoroshiro128Plus::seed_from_u64(1);
        let perm = Poseidon2BabyBear::new_from_rng_128(&mut rng);

        perm.permute_mut(&mut input);
        assert_eq!(input, expected);
    }

    /// Test on a roughly random input.
    /// This random input is generated by the following sage code:
    /// set_random_seed(24)
    /// vector([BB.random_element() for t in range(24)]).
    #[test]
    fn test_poseidon2_width_24_random() {
        let mut input: [F; 24] = BabyBear::new_array([
            886409618, 1327899896, 1902407911, 591953491, 648428576, 1844789031, 1198336108,
            355597330, 1799586834, 59617783, 790334801, 1968791836, 559272107, 31054313,
            1042221543, 474748436, 135686258, 263665994, 1962340735, 1741539604, 449439011,
            1131357108, 50869465, 1589724894,
        ]);

        let expected: [F; 24] = BabyBear::new_array([
            249424342, 562262148, 757431114, 354243402, 57767055, 976981973, 1393169022,
            1774550827, 1527742125, 1019514605, 1776327602, 266236737, 1412355182, 1070239213,
            426390978, 1775539440, 1527732214, 1101406020, 1417710778, 1699632661, 413672313,
            820348291, 1067197851, 1669055675,
        ]);

        let mut rng = Xoroshiro128Plus::seed_from_u64(1);
        let perm = Poseidon2BabyBear::new_from_rng_128(&mut rng);

        perm.permute_mut(&mut input);

        assert_eq!(input, expected);
    }

    /// Test the generic internal layer against the optimized internal layer
    /// for a random input of width 16.
    #[test]
    fn test_generic_internal_linear_layer_16() {
        let mut rng = Xoroshiro128Plus::seed_from_u64(1);
        let mut input1: [F; 16] = rng.random();
        let mut input2 = input1;

        let part_sum: F = input1[1..].iter().copied().sum();
        let full_sum = part_sum + input1[0];

        input1[0] = part_sum - input1[0];

        BabyBearInternalLayerParameters::internal_layer_mat_mul(&mut input1, full_sum);
        BabyBearInternalLayerParameters::generic_internal_linear_layer(&mut input2);

        assert_eq!(input1, input2);
    }

    /// Test the generic internal layer against the optimized internal layer
    /// for a random input of width 24.
    #[test]
    fn test_generic_internal_linear_layer_24() {
        let mut rng = Xoroshiro128Plus::seed_from_u64(1);
        let mut input1: [F; 24] = rng.random();
        let mut input2 = input1;

        let part_sum: F = input1[1..].iter().copied().sum();
        let full_sum = part_sum + input1[0];

        input1[0] = part_sum - input1[0];

        BabyBearInternalLayerParameters::internal_layer_mat_mul(&mut input1, full_sum);
        BabyBearInternalLayerParameters::generic_internal_linear_layer(&mut input2);

        assert_eq!(input1, input2);
    }

    #[test]
    fn test_default_babybear_poseidon2_width_16() {
        let mut input: [F; 16] = BabyBear::new_array([
            894848333, 1437655012, 1200606629, 1690012884, 71131202, 1749206695, 1717947831,
            120589055, 19776022, 42382981, 1831865506, 724844064, 171220207, 1299207443, 227047920,
            1783754913,
        ]);

        let expected: [F; 16] = BabyBear::new_array([
            516096821, 90309867, 1101817252, 1660784290, 360715097, 1789519026, 1788910906,
            563338433, 319524748, 1741414159, 1650859320, 894311162, 1121347488, 1692793758,
            1052633829, 1344246938,
        ]);

        let perm = default_babybear_poseidon2_16();
        perm.permute_mut(&mut input);

        assert_eq!(input, expected);
    }

    #[test]
    fn test_default_babybear_poseidon2_width_24() {
        let mut input: [F; 24] = BabyBear::new_array([
            886409618, 1327899896, 1902407911, 591953491, 648428576, 1844789031, 1198336108,
            355597330, 1799586834, 59617783, 790334801, 1968791836, 559272107, 31054313,
            1042221543, 474748436, 135686258, 263665994, 1962340735, 1741539604, 2026927696,
            449439011, 1131357108, 50869465,
        ]);

        let expected: [F; 24] = BabyBear::new_array([
            882297297, 1264077610, 512812497, 782602970, 867738552, 1251075457, 309180082,
            340784773, 524041877, 351272188, 404451680, 15001466, 322926653, 1773004150,
            1718440818, 674682955, 1154713225, 1719133502, 324232301, 1005243141, 443371079,
            268735940, 770060019, 718377682,
        ]);

        let perm = default_babybear_poseidon2_24();
        perm.permute_mut(&mut input);

        assert_eq!(input, expected);
    }

    #[test]
    fn test_default_babybear_poseidon2_width_32() {
        let mut input: [F; 32] = BabyBear::new_array([
            377682961, 1957793603, 980981814, 6565119, 1583211709, 176593168, 1672635515,
            226854190, 1096634172, 1317773742, 1472230830, 1621534427, 559807320, 1484241910,
            1847825942, 3491998, 950152427, 1935451636, 275759400, 227625951, 1271142011,
            1492341973, 1502961189, 147694103, 1939834518, 1449972249, 1822424048, 1518111482,
            714203295, 383863563, 411489861, 1253612091,
        ]);

        let expected: [F; 32] = BabyBear::new_array([
            303440672, 985419733, 780962554, 1395263823, 188752116, 1348917749, 677984158,
            667170017, 97281439, 178741618, 1770541242, 1894441262, 847173187, 1374453653,
            1242356754, 1485142795, 1019698843, 334329175, 540395852, 918117757, 1288401072,
            508687761, 996827321, 1660764537, 546969284, 1848510002, 334793951, 736596659,
            1928951999, 1444080616, 55017699, 1832626373,
        ]);

        let perm = default_babybear_poseidon2_32();
        perm.permute_mut(&mut input);

        assert_eq!(input, expected);
    }

    /// Test the generic internal layer against the optimized internal layer
    /// for a random input of width 32.
    #[test]
    fn test_generic_internal_linear_layer_32() {
        let mut rng = Xoroshiro128Plus::seed_from_u64(1);
        let mut input1: [F; 32] = rng.random();
        let mut input2 = input1;

        let part_sum: F = input1[1..].iter().copied().sum();
        let full_sum = part_sum + input1[0];

        input1[0] = part_sum - input1[0];

        BabyBearInternalLayerParameters::internal_layer_mat_mul(&mut input1, full_sum);
        BabyBearInternalLayerParameters::generic_internal_linear_layer(&mut input2);

        assert_eq!(input1, input2);
    }
}