russell_tensor 3.2.1

use crate::{ADD, SQRT_2};
use crate::{Tensor2, Tensor4};

/// Performs the underbar dyadic product between two Tensor2 resulting in a (general) Tensor4
///
/// Computes:
///
/// ```text
/// ADD: D += s A ⊗ B  or  SET: D = s A ⊗ B
///               ‾                     ‾
/// ```
///
/// With Cartesian components (example with SET):
///
/// ```text
/// Dᵢⱼₖₗ = s Aᵢₗ Bⱼₖ
/// ```
///
/// **Important:** The result is **not** necessarily minor-symmetric; therefore `D` is general.
///
/// # Output
///
/// * `dd` -- the tensor `D`
///
/// # Input
///
/// * `op` -- operation: ADD or SET
/// * `s` -- the multiplier
/// * `aa` -- first tensor
/// * `bb` -- second tensor
pub fn t2_udyad_t2<const N: usize>(dd: &mut Tensor4<9>, op: u8, s: f64, aa: &Tensor2<N>, bb: &Tensor2<N>) {
    t2_udyad_t2_vec::<N>(dd, op, s, aa.as_vec(), bb.as_vec());
}

/// Internal (unrolled) underbar dyadic product on raw Kelvin-Mandel vectors.
#[rustfmt::skip]
#[inline]
pub(crate) fn t2_udyad_t2_vec<const N:usize>(dd: &mut Tensor4<9>, op: u8, s: f64, a: &[f64; N], b: &[f64; N]) {
    let tsq2 = 2.0 * SQRT_2;
    if op == ADD {
        if N == 4 {
            dd.add(0, 0, s*a[0]*b[0]);
            dd.add(0, 1, s*(a[3]*b[3])/2.0);
            dd.add(0, 2, 0.0);
            dd.add(0, 3, s*(a[3]*b[0] + a[0]*b[3])/2.0);
            dd.add(0, 4, 0.0);
            dd.add(0, 5, 0.0);
            dd.add(0, 6, s*(a[3]*b[0] - a[0]*b[3])/2.0);
            dd.add(0, 7, 0.0);
            dd.add(0, 8, 0.0);

            dd.add(1, 0, s*(a[3]*b[3])/2.0);
            dd.add(1, 1, s*a[1]*b[1]);
            dd.add(1, 2, 0.0);
            dd.add(1, 3, s*(a[3]*b[1] + a[1]*b[3])/2.0);
            dd.add(1, 4, 0.0);
            dd.add(1, 5, 0.0);
            dd.add(1, 6, s*(-(a[3]*b[1]) + a[1]*b[3])/2.0);
            dd.add(1, 7, 0.0);
            dd.add(1, 8, 0.0);

            dd.add(2, 0, 0.0);
            dd.add(2, 1, 0.0);
            dd.add(2, 2, s*a[2]*b[2]);
            dd.add(2, 3, 0.0);
            dd.add(2, 4, 0.0);
            dd.add(2, 5, 0.0);
            dd.add(2, 6, 0.0);
            dd.add(2, 7, 0.0);
            dd.add(2, 8, 0.0);

            dd.add(3, 0, s*(a[3]*b[0] + a[0]*b[3])/2.0);
            dd.add(3, 1, s*(a[3]*b[1] + a[1]*b[3])/2.0);
            dd.add(3, 2, 0.0);
            dd.add(3, 3, s*(a[1]*b[0] + a[0]*b[1] + a[3]*b[3])/2.0);
            dd.add(3, 4, 0.0);
            dd.add(3, 5, 0.0);
            dd.add(3, 6, s*(a[1]*b[0] - a[0]*b[1])/2.0);
            dd.add(3, 7, 0.0);
            dd.add(3, 8, 0.0);

            dd.add(4, 0, 0.0);
            dd.add(4, 1, 0.0);
            dd.add(4, 2, 0.0);
            dd.add(4, 3, 0.0);
            dd.add(4, 4, s*(a[2]*b[1] + a[1]*b[2])/2.0);
            dd.add(4, 5, s*(a[3]*b[2] + a[2]*b[3])/tsq2);
            dd.add(4, 6, 0.0);
            dd.add(4, 7, s*(a[2]*b[1] - a[1]*b[2])/2.0);
            dd.add(4, 8, s*(-(a[3]*b[2]) + a[2]*b[3])/tsq2);

            dd.add(5, 0, 0.0);
            dd.add(5, 1, 0.0);
            dd.add(5, 2, 0.0);
            dd.add(5, 3, 0.0);
            dd.add(5, 4, s*(a[3]*b[2] + a[2]*b[3])/tsq2);
            dd.add(5, 5, s*(a[2]*b[0] + a[0]*b[2])/2.0);
            dd.add(5, 6, 0.0);
            dd.add(5, 7, s*(-(a[3]*b[2]) + a[2]*b[3])/tsq2);
            dd.add(5, 8, s*(a[2]*b[0] - a[0]*b[2])/2.0);

            dd.add(6, 0, s*(-(a[3]*b[0]) + a[0]*b[3])/2.0);
            dd.add(6, 1, s*(a[3]*b[1] - a[1]*b[3])/2.0);
            dd.add(6, 2, 0.0);
            dd.add(6, 3, s*(-(a[1]*b[0]) + a[0]*b[1])/2.0);
            dd.add(6, 4, 0.0);
            dd.add(6, 5, 0.0);
            dd.add(6, 6, s*(-(a[1]*b[0]) - a[0]*b[1] + a[3]*b[3])/2.0);
            dd.add(6, 7, 0.0);
            dd.add(6, 8, 0.0);

            dd.add(7, 0, 0.0);
            dd.add(7, 1, 0.0);
            dd.add(7, 2, 0.0);
            dd.add(7, 3, 0.0);
            dd.add(7, 4, s*(-(a[2]*b[1]) + a[1]*b[2])/2.0);
            dd.add(7, 5, s*(a[3]*b[2] - a[2]*b[3])/tsq2);
            dd.add(7, 6, 0.0);
            dd.add(7, 7, s*(-(a[2]*b[1]) - a[1]*b[2])/2.0);
            dd.add(7, 8, s*(-(a[3]*b[2] + a[2]*b[3])/tsq2));

            dd.add(8, 0, 0.0);
            dd.add(8, 1, 0.0);
            dd.add(8, 2, 0.0);
            dd.add(8, 3, 0.0);
            dd.add(8, 4, s*(a[3]*b[2] - a[2]*b[3])/tsq2);
            dd.add(8, 5, s*(-(a[2]*b[0]) + a[0]*b[2])/2.0);
            dd.add(8, 6, 0.0);
            dd.add(8, 7, s*(-(a[3]*b[2] + a[2]*b[3])/tsq2));
            dd.add(8, 8, s*(-(a[2]*b[0]) - a[0]*b[2])/2.0);
        } else if N == 6 {
            dd.add(0, 0, s*a[0]*b[0]);
            dd.add(0, 1, s*(a[3]*b[3])/2.0);
            dd.add(0, 2, s*(a[5]*b[5])/2.0);
            dd.add(0, 3, s*(a[3]*b[0] + a[0]*b[3])/2.0);
            dd.add(0, 4, s*(a[5]*b[3] + a[3]*b[5])/tsq2);
            dd.add(0, 5, s*(a[5]*b[0] + a[0]*b[5])/2.0);
            dd.add(0, 6, s*(a[3]*b[0] - a[0]*b[3])/2.0);
            dd.add(0, 7, s*(a[5]*b[3] - a[3]*b[5])/tsq2);
            dd.add(0, 8, s*(a[5]*b[0] - a[0]*b[5])/2.0);

            dd.add(1, 0, s*(a[3]*b[3])/2.0);
            dd.add(1, 1, s*a[1]*b[1]);
            dd.add(1, 2, s*(a[4]*b[4])/2.0);
            dd.add(1, 3, s*(a[3]*b[1] + a[1]*b[3])/2.0);
            dd.add(1, 4, s*(a[4]*b[1] + a[1]*b[4])/2.0);
            dd.add(1, 5, s*(a[4]*b[3] + a[3]*b[4])/tsq2);
            dd.add(1, 6, s*(-(a[3]*b[1]) + a[1]*b[3])/2.0);
            dd.add(1, 7, s*(a[4]*b[1] - a[1]*b[4])/2.0);
            dd.add(1, 8, s*(a[4]*b[3] - a[3]*b[4])/tsq2);

            dd.add(2, 0, s*(a[5]*b[5])/2.0);
            dd.add(2, 1, s*(a[4]*b[4])/2.0);
            dd.add(2, 2, s*a[2]*b[2]);
            dd.add(2, 3, s*(a[ 5]*b[4] + a[4]*b[5])/tsq2);
            dd.add(2, 4, s*(a[4]*b[2] + a[2]*b[4])/2.0);
            dd.add(2, 5, s*(a[5]*b[2] + a[2]*b[5])/2.0);
            dd.add(2, 6, s*(-(a[5]*b[4]) + a[4]*b[5])/tsq2);
            dd.add(2, 7, s*(-(a[4]*b[2]) + a[2]*b[4])/2.0);
            dd.add(2, 8, s*(-(a[5]*b[2]) + a[2]*b[5])/2.0);

            dd.add(3, 0, s*(a[3]*b[0] + a[0]*b[3])/2.0);
            dd.add(3, 1, s*(a[3]*b[1] + a[1]*b[3])/2.0);
            dd.add(3, 2, s*(a[5]*b[4] + a[4]*b[5])/tsq2);
            dd.add(3, 3, s*(a[1]*b[0] + a[0]*b[1] + a[3]*b[3])/2.0);
            dd.add(3, 4, s*(SQRT_2*a[5]*b[1] + a[4]*b[3] + a[3]*b[4] + SQRT_2*a[1]*b[5])/4.0);
            dd.add(3, 5, s*(SQRT_2*a[4]*b[0] + a[5]*b[3] + SQRT_2*a[0]*b[4] + a[3]*b[5])/4.0);
            dd.add(3, 6, s*(a[1]*b[0] - a[0]*b[1])/2.0);
            dd.add(3, 7, s*(SQRT_2*a[5]*b[1] + a[4]*b[3] - a[3]*b[4] - SQRT_2*a[1]*b[5])/4.0);
            dd.add(3, 8, s*(SQRT_2*a[4]*b[0] + a[5]*b[3] - SQRT_2*a[0]*b[4] - a[3]*b[5])/4.0);

            dd.add(4, 0, s*(a[5]*b[3] + a[3]*b[5])/tsq2);
            dd.add(4, 1, s*(a[4]*b[1] + a[1]*b[4])/2.0);
            dd.add(4, 2, s*(a[4]*b[2] + a[2]*b[4])/2.0);
            dd.add(4, 3, s*(SQRT_2*a[5]*b[1] + a[4]*b[3] + a[3]*b[4] + SQRT_2*a[1]*b[5])/4.0);
            dd.add(4, 4, s*(a[2]*b[1] + a[1]*b[2] + a[4]*b[4])/2.0);
            dd.add(4, 5, s*(SQRT_2*a[3]*b[2] + SQRT_2*a[2]*b[3] + a[5]*b[4] + a[4]*b[5])/4.0);
            dd.add(4, 6, s*(-(SQRT_2*a[5]*b[1]) + a[4]*b[3] - a[3]*b[4] + SQRT_2*a[1]*b[5])/4.0);
            dd.add(4, 7, s*(a[2]*b[1] - a[1]*b[2])/2.0);
            dd.add(4, 8, s*(-(SQRT_2*a[3]*b[2]) + SQRT_2*a[2]*b[3] - a[5]*b[4] + a[4]*b[5])/4.0);

            dd.add(5, 0, s*(a[5]*b[0] + a[0]*b[5])/2.0);
            dd.add(5, 1, s*(a[4]*b[3] + a[3]*b[4])/tsq2);
            dd.add(5, 2, s*(a[5]*b[2] + a[2]*b[5])/2.0);
            dd.add(5, 3, s*(SQRT_2*a[4]*b[0] + a[5]*b[3] + SQRT_2*a[0]*b[4] + a[3]*b[5])/4.0);
            dd.add(5, 4, s*(SQRT_2*a[3]*b[2] + SQRT_2*a[2]*b[3] + a[5]*b[4] + a[4]*b[5])/4.0);
            dd.add(5, 5, s*(a[2]*b[0] + a[0]*b[2] + a[5]*b[5])/2.0);
            dd.add(5, 6, s*(SQRT_2*a[4]*b[0] - a[5]*b[3] - SQRT_2*a[0]*b[4] + a[3]*b[5])/4.0);
            dd.add(5, 7, s*(-(SQRT_2*a[3]*b[2]) + SQRT_2*a[2]*b[3] + a[5]*b[4] - a[4]*b[5])/4.0);
            dd.add(5, 8, s*(a[2]*b[0] - a[0]*b[2])/2.0);

            dd.add(6, 0, s*(-(a[3]*b[0]) + a[0]*b[3])/2.0);
            dd.add(6, 1, s*(a[3]*b[1] - a[1]*b[3])/2.0);
            dd.add(6, 2, s*(a[5]*b[4] - a[4]*b[5])/tsq2);
            dd.add(6, 3, s*(-(a[1]*b[0]) + a[0]*b[1])/2.0);
            dd.add(6, 4, s*(SQRT_2*a[5]*b[1] - a[4]*b[3] + a[3]*b[4] - SQRT_2*a[1]*b[5])/4.0);
            dd.add(6, 5, s*(-(SQRT_2*a[4]*b[0]) + a[5]*b[3] + SQRT_2*a[0]*b[4] - a[3]*b[5])/4.0);
            dd.add(6, 6, s*(-(a[1]*b[0]) - a[0]*b[1] + a[3]*b[3])/2.0);
            dd.add(6, 7, s*(SQRT_2*a[5]*b[1] - a[4]*b[3] - a[3]*b[4] + SQRT_2*a[1]*b[5])/4.0);
            dd.add(6, 8, s*(-(SQRT_2*a[4]*b[0]) + a[5]*b[3] - SQRT_2*a[0]*b[4] + a[3]*b[5])/4.0);

            dd.add(7, 0, s*(-(a[5]*b[3]) + a[3]*b[5])/tsq2);
            dd.add(7, 1, s*(-(a[4]*b[1]) + a[1]*b[4])/2.0);
            dd.add(7, 2, s*(a[4]*b[2] - a[2]*b[4])/2.0);
            dd.add(7, 3, s*(-(SQRT_2*a[5]*b[1]) - a[4]*b[3] + a[3]*b[4] + SQRT_2*a[1]*b[5])/4.0);
            dd.add(7, 4, s*(-(a[2]*b[1]) + a[1]*b[2])/2.0);
            dd.add(7, 5, s*(SQRT_2*a[3]*b[2] - SQRT_2*a[2]*b[3] - a[5]*b[4] + a[4]*b[5])/4.0);
            dd.add(7, 6, s*(SQRT_2*a[5]*b[1] - a[4]*b[3] - a[3]*b[4] + SQRT_2*a[1]*b[5])/4.0);
            dd.add(7, 7, s*(-(a[2]*b[1]) - a[1]*b[2] + a[4]*b[4])/2.0);
            dd.add(7, 8, s*(-(SQRT_2*a[3]*b[2]) - SQRT_2*a[2]*b[3] + a[5]*b[4] + a[4]*b[5])/4.0);

            dd.add(8, 0, s*(-(a[5]*b[0]) + a[0]*b[5])/2.0);
            dd.add(8, 1, s*(-(a[4]*b[3]) + a[3]*b[4])/tsq2);
            dd.add(8, 2, s*(a[5]*b[2] - a[2]*b[5])/2.0);
            dd.add(8, 3, s*(-(SQRT_2*a[4]*b[0]) - a[5]*b[3] + SQRT_2*a[0]*b[4] + a[3]*b[5])/4.0);
            dd.add(8, 4, s*(SQRT_2*a[3]*b[2] - SQRT_2*a[2]*b[3] + a[5]*b[4] - a[4]*b[5])/4.0);
            dd.add(8, 5, s*(-(a[2]*b[0]) + a[0]*b[2])/2.0);
            dd.add(8, 6, s*(-(SQRT_2*a[4]*b[0]) + a[5]*b[3] - SQRT_2*a[0]*b[4] + a[3]*b[5])/4.0);
            dd.add(8, 7, s*(-(SQRT_2*a[3]*b[2]) - SQRT_2*a[2]*b[3] + a[5]*b[4] + a[4]*b[5])/4.0);
            dd.add(8, 8, s*(-(a[2]*b[0]) - a[0]*b[2] + a[5]*b[5])/2.0);
        } else {
            debug_assert!(N == 9);

            dd.add(0, 0, s*a[0]*b[0]);
            dd.add(0, 1, s*((a[3] + a[6])*(b[3] + b[6]))/2.0);
            dd.add(0, 2, s*((a[5] + a[8])*(b[5] + b[8]))/2.0);
            dd.add(0, 3, s*(a[3]*b[0] + a[6]*b[0] + a[0]*(b[3] + b[6]))/2.0);
            dd.add(0, 4, s*((a[5] + a[8])*(b[3] + b[6]) + (a[3] + a[6])*(b[5] + b[8]))/tsq2);
            dd.add(0, 5, s*(a[5]*b[0] + a[8]*b[0] + a[0]*(b[5] + b[8]))/2.0);
            dd.add(0, 6, s*(a[3]*b[0] + a[6]*b[0] - a[0]*(b[3] + b[6]))/2.0);
            dd.add(0, 7, s*((a[5] + a[8])*(b[3] + b[6]) - (a[3] + a[6])*(b[5] + b[8]))/tsq2);
            dd.add(0, 8, s*(a[5]*b[0] + a[8]*b[0] - a[0]*(b[5] + b[8]))/2.0);

            dd.add(1, 0, s*((a[3] - a[6])*(b[3] - b[6]))/2.0);
            dd.add(1, 1, s*a[1]*b[1]);
            dd.add(1, 2, s*((a[4] + a[7])*(b[4] + b[7]))/2.0);
            dd.add(1, 3, s*(a[3]*b[1] - a[6]*b[1] + a[1]*(b[3] - b[6]))/2.0);
            dd.add(1, 4, s*(a[4]*b[1] + a[7]*b[1] + a[1]*(b[4] + b[7]))/2.0);
            dd.add(1, 5, s*((a[4] + a[7])*(b[3] - b[6]) + (a[3] - a[6])*(b[4] + b[7]))/tsq2);
            dd.add(1, 6, s*(-(a[3]*b[1]) + a[6]*b[1] + a[1]*(b[3] - b[6]))/2.0);
            dd.add(1, 7, s*(a[4]*b[1] + a[7]*b[1] - a[1]*(b[4] + b[7]))/2.0);
            dd.add(1, 8, s*((a[4] + a[7])*(b[3] - b[6]) - (a[3] - a[6])*(b[4] + b[7]))/tsq2);

            dd.add(2, 0, s*((a[5] - a[8])*(b[5] - b[8]))/2.0);
            dd.add(2, 1, s*((a[4] - a[7])*(b[4] - b[7]))/2.0);
            dd.add(2, 2, s*a[2]*b[2]);
            dd.add(2, 3, s*((a[5] - a[8])*(b[4] - b[7]) + (a[4] - a[7])*(b[5] - b[8]))/tsq2);
            dd.add(2, 4, s*(a[4]*b[2] - a[7]*b[2] + a[2]*(b[4] - b[7]))/2.0);
            dd.add(2, 5, s*(a[5]*b[2] - a[8]*b[2] + a[2]*(b[5] - b[8]))/2.0);
            dd.add(2, 6, s*(-((a[5] - a[8])*(b[4] - b[7])) + (a[4] - a[7])*(b[5] - b[8]))/tsq2);
            dd.add(2, 7, s*(-(a[4]*b[2]) + a[7]*b[2] + a[2]*(b[4] - b[7]))/2.0);
            dd.add(2, 8, s*(-(a[5]*b[2]) + a[8]*b[2] + a[2]*(b[5] - b[8]))/2.0);

            dd.add(3, 0, s*(a[3]*b[0] - a[6]*b[0] + a[0]*(b[3] - b[6]))/2.0);
            dd.add(3, 1, s*(a[3]*b[1] + a[6]*b[1] + a[1]*(b[3] + b[6]))/2.0);
            dd.add(3, 2, s*((a[5] + a[8])*(b[4] + b[7]) + (a[4] + a[7])*(b[5] + b[8]))/tsq2);
            dd.add(3, 3, s*(a[1]*b[0] + a[0]*b[1] + a[3]*b[3] - a[6]*b[6])/2.0);
            dd.add(3, 4, s*(SQRT_2*(a[5] + a[8])*b[1] + (a[4] + a[7])*(b[3] + b[6]) + (a[3] + a[6])*(b[4] + b[7]) + SQRT_2*a[1]*(b[5] + b[8]))/4.0);
            dd.add(3, 5, s*(SQRT_2*(a[4] + a[7])*b[0] + (a[5] + a[8])*(b[3] - b[6]) + SQRT_2*a[0]*(b[4] + b[7]) + (a[3] - a[6])*(b[5] + b[8]))/4.0);
            dd.add(3, 6, s*(a[1]*b[0] - a[0]*b[1] + a[6]*b[3] - a[3]*b[6])/2.0);
            dd.add(3, 7, s*(SQRT_2*(a[5] + a[8])*b[1] + (a[4] + a[7])*(b[3] + b[6]) - (a[3] + a[6])*(b[4] + b[7]) - SQRT_2*a[1]*(b[5] + b[8]))/4.0);
            dd.add(3, 8, s*(SQRT_2*(a[4] + a[7])*b[0] + (a[5] + a[8])*(b[3] - b[6]) - SQRT_2*a[0]*(b[4] + b[7]) - (a[3] - a[6])*(b[5] + b[8]))/4.0);

            dd.add(4, 0, s*((a[5] - a[8])*(b[3] - b[6]) + (a[3] - a[6])*(b[5] - b[8]))/tsq2);
            dd.add(4, 1, s*(a[4]*b[1] - a[7]*b[1] + a[1]*(b[4] - b[7]))/2.0);
            dd.add(4, 2, s*(a[4]*b[2] + a[7]*b[2] + a[2]*(b[4] + b[7]))/2.0);
            dd.add(4, 3, s*(SQRT_2*(a[5] - a[8])*b[1] + (a[4] - a[7])*(b[3] - b[6]) + (a[3] - a[6])*(b[4] - b[7]) + SQRT_2*a[1]*(b[5] - b[8]))/4.0);
            dd.add(4, 4, s*(a[2]*b[1] + a[1]*b[2] + a[4]*b[4] - a[7]*b[7])/2.0);
            dd.add(4, 5, s*(SQRT_2*(a[3] - a[6])*b[2] + SQRT_2*a[2]*(b[3] - b[6]) + (a[5] - a[8])*(b[4] + b[7]) + (a[4] + a[7])*(b[5] - b[8]))/4.0);
            dd.add(4, 6, s*(-(SQRT_2*(a[5] - a[8])*b[1]) + (a[4] - a[7])*(b[3] - b[6]) - (a[3] - a[6])*(b[4] - b[7]) + SQRT_2*a[1]*(b[5] - b[8]))/4.0);
            dd.add(4, 7, s*(a[2]*b[1] - a[1]*b[2] + a[7]*b[4] - a[4]*b[7])/2.0);
            dd.add(4, 8, s*(-(SQRT_2*(a[3] - a[6])*b[2]) + SQRT_2*a[2]*(b[3] - b[6]) - (a[5] - a[8])*(b[4] + b[7]) + (a[4] + a[7])*(b[5] - b[8]))/4.0);

            dd.add(5, 0, s*(a[5]*b[0] - a[8]*b[0] + a[0]*(b[5] - b[8]))/2.0);
            dd.add(5, 1, s*((a[4] - a[7])*(b[3] + b[6]) + (a[3] + a[6])*(b[4] - b[7]))/tsq2);
            dd.add(5, 2, s*(a[5]*b[2] + a[8]*b[2] + a[2]*(b[5] + b[8]))/2.0);
            dd.add(5, 3, s*(SQRT_2*(a[4] - a[7])*b[0] + (a[5] - a[8])*(b[3] + b[6]) + SQRT_2*a[0]*(b[4] - b[7]) + (a[3] + a[6])*(b[5] - b[8]))/4.0);
            dd.add(5, 4, s*(SQRT_2*(a[3] + a[6])*b[2] + SQRT_2*a[2]*(b[3] + b[6]) + (a[5] + a[8])*(b[4] - b[7]) + (a[4] - a[7])*(b[5] + b[8]))/4.0);
            dd.add(5, 5, s*(a[2]*b[0] + a[0]*b[2] + a[5]*b[5] - a[8]*b[8])/2.0);
            dd.add(5, 6, s*(SQRT_2*(a[4] - a[7])*b[0] - (a[5] - a[8])*(b[3] + b[6]) - SQRT_2*a[0]*(b[4] - b[7]) + (a[3] + a[6])*(b[5] - b[8]))/4.0);
            dd.add(5, 7, s*(-(SQRT_2*(a[3] + a[6])*b[2]) + SQRT_2*a[2]*(b[3] + b[6]) + (a[5] + a[8])*(b[4] - b[7]) - (a[4] - a[7])*(b[5] + b[8]))/4.0);
            dd.add(5, 8, s*(a[2]*b[0] - a[0]*b[2] + a[8]*b[5] - a[5]*b[8])/2.0);

            dd.add(6, 0, s*(-(a[3]*b[0]) + a[6]*b[0] + a[0]*(b[3] - b[6]))/2.0);
            dd.add(6, 1, s*(a[3]*b[1] + a[6]*b[1] - a[1]*(b[3] + b[6]))/2.0);
            dd.add(6, 2, s*((a[5] + a[8])*(b[4] + b[7]) - (a[4] + a[7])*(b[5] + b[8]))/tsq2);
            dd.add(6, 3, s*(-(a[1]*b[0]) + a[0]*b[1] + a[6]*b[3] - a[3]*b[6])/2.0);
            dd.add(6, 4, s*(SQRT_2*(a[5] + a[8])*b[1] - (a[4] + a[7])*(b[3] + b[6]) + (a[3] + a[6])*(b[4] + b[7]) - SQRT_2*a[1]*(b[5] + b[8]))/4.0);
            dd.add(6, 5, s*(-(SQRT_2*(a[4] + a[7])*b[0]) + (a[5] + a[8])*(b[3] - b[6]) + SQRT_2*a[0]*(b[4] + b[7]) - (a[3] - a[6])*(b[5] + b[8]))/4.0);
            dd.add(6, 6, s*(-(a[1]*b[0]) - a[0]*b[1] + a[3]*b[3] - a[6]*b[6])/2.0);
            dd.add(6, 7, s*(SQRT_2*(a[5] + a[8])*b[1] - (a[4] + a[7])*(b[3] + b[6]) - (a[3] + a[6])*(b[4] + b[7]) + SQRT_2*a[1]*(b[5] + b[8]))/4.0);
            dd.add(6, 8, s*(-(SQRT_2*(a[4] + a[7])*b[0]) + (a[5] + a[8])*(b[3] - b[6]) - SQRT_2*a[0]*(b[4] + b[7]) + (a[3] - a[6])*(b[5] + b[8]))/4.0);

            dd.add(7, 0, s*(-((a[5] - a[8])*(b[3] - b[6])) + (a[3] - a[6])*(b[5] - b[8]))/tsq2);
            dd.add(7, 1, s*(-(a[4]*b[1]) + a[7]*b[1] + a[1]*(b[4] - b[7]))/2.0);
            dd.add(7, 2, s*(a[4]*b[2] + a[7]*b[2] - a[2]*(b[4] + b[7]))/2.0);
            dd.add(7, 3, s*(-(SQRT_2*(a[5] - a[8])*b[1]) - (a[4] - a[7])*(b[3] - b[6]) + (a[3] - a[6])*(b[4] - b[7]) + SQRT_2*a[1]*(b[5] - b[8]))/4.0);
            dd.add(7, 4, s*(-(a[2]*b[1]) + a[1]*b[2] + a[7]*b[4] - a[4]*b[7])/2.0);
            dd.add(7, 5, s*(SQRT_2*(a[3] - a[6])*b[2] - SQRT_2*a[2]*(b[3] - b[6]) - (a[5] - a[8])*(b[4] + b[7]) + (a[4] + a[7])*(b[5] - b[8]))/4.0);
            dd.add(7, 6, s*(SQRT_2*(a[5] - a[8])*b[1] - (a[4] - a[7])*(b[3] - b[6]) - (a[3] - a[6])*(b[4] - b[7]) + SQRT_2*a[1]*(b[5] - b[8]))/4.0);
            dd.add(7, 7, s*(-(a[2]*b[1]) - a[1]*b[2] + a[4]*b[4] - a[7]*b[7])/2.0);
            dd.add(7, 8, s*(-(SQRT_2*(a[3] - a[6])*b[2]) - SQRT_2*a[2]*(b[3] - b[6]) + (a[5] - a[8])*(b[4] + b[7]) + (a[4] + a[7])*(b[5] - b[8]))/4.0);

            dd.add(8, 0, s*(-(a[5]*b[0]) + a[8]*b[0] + a[0]*(b[5] - b[8]))/2.0);
            dd.add(8, 1, s*(-((a[4] - a[7])*(b[3] + b[6])) + (a[3] + a[6])*(b[4] - b[7]))/tsq2);
            dd.add(8, 2, s*(a[5]*b[2] + a[8]*b[2] - a[2]*(b[5] + b[8]))/2.0);
            dd.add(8, 3, s*(-(SQRT_2*(a[4] - a[7])*b[0]) - (a[5] - a[8])*(b[3] + b[6]) + SQRT_2*a[0]*(b[4] - b[7]) + (a[3] + a[6])*(b[5] - b[8]))/4.0);
            dd.add(8, 4, s*(SQRT_2*(a[3] + a[6])*b[2] - SQRT_2*a[2]*(b[3] + b[6]) + (a[5] + a[8])*(b[4] - b[7]) - (a[4] - a[7])*(b[5] + b[8]))/4.0);
            dd.add(8, 5, s*(-(a[2]*b[0]) + a[0]*b[2] + a[8]*b[5] - a[5]*b[8])/2.0);
            dd.add(8, 6, s*(-(SQRT_2*(a[4] - a[7])*b[0]) + (a[5] - a[8])*(b[3] + b[6]) - SQRT_2*a[0]*(b[4] - b[7]) + (a[3] + a[6])*(b[5] - b[8]))/4.0);
            dd.add(8, 7, s*(-(SQRT_2*(a[3] + a[6])*b[2]) - SQRT_2*a[2]*(b[3] + b[6]) + (a[5] + a[8])*(b[4] - b[7]) + (a[4] - a[7])*(b[5] + b[8]))/4.0);
            dd.add(8, 8, s*(-(a[2]*b[0]) - a[0]*b[2] + a[5]*b[5] - a[8]*b[8])/2.0);
        }
    } else {
        if N == 4 {
            dd.set(0, 0, s*a[0]*b[0]);
            dd.set(0, 1, s*(a[3]*b[3])/2.0);
            dd.set(0, 2, 0.0);
            dd.set(0, 3, s*(a[3]*b[0] + a[0]*b[3])/2.0);
            dd.set(0, 4, 0.0);
            dd.set(0, 5, 0.0);
            dd.set(0, 6, s*(a[3]*b[0] - a[0]*b[3])/2.0);
            dd.set(0, 7, 0.0);
            dd.set(0, 8, 0.0);

            dd.set(1, 0, s*(a[3]*b[3])/2.0);
            dd.set(1, 1, s*a[1]*b[1]);
            dd.set(1, 2, 0.0);
            dd.set(1, 3, s*(a[3]*b[1] + a[1]*b[3])/2.0);
            dd.set(1, 4, 0.0);
            dd.set(1, 5, 0.0);
            dd.set(1, 6, s*(-(a[3]*b[1]) + a[1]*b[3])/2.0);
            dd.set(1, 7, 0.0);
            dd.set(1, 8, 0.0);

            dd.set(2, 0, 0.0);
            dd.set(2, 1, 0.0);
            dd.set(2, 2, s*a[2]*b[2]);
            dd.set(2, 3, 0.0);
            dd.set(2, 4, 0.0);
            dd.set(2, 5, 0.0);
            dd.set(2, 6, 0.0);
            dd.set(2, 7, 0.0);
            dd.set(2, 8, 0.0);

            dd.set(3, 0, s*(a[3]*b[0] + a[0]*b[3])/2.0);
            dd.set(3, 1, s*(a[3]*b[1] + a[1]*b[3])/2.0);
            dd.set(3, 2, 0.0);
            dd.set(3, 3, s*(a[1]*b[0] + a[0]*b[1] + a[3]*b[3])/2.0);
            dd.set(3, 4, 0.0);
            dd.set(3, 5, 0.0);
            dd.set(3, 6, s*(a[1]*b[0] - a[0]*b[1])/2.0);
            dd.set(3, 7, 0.0);
            dd.set(3, 8, 0.0);

            dd.set(4, 0, 0.0);
            dd.set(4, 1, 0.0);
            dd.set(4, 2, 0.0);
            dd.set(4, 3, 0.0);
            dd.set(4, 4, s*(a[2]*b[1] + a[1]*b[2])/2.0);
            dd.set(4, 5, s*(a[3]*b[2] + a[2]*b[3])/tsq2);
            dd.set(4, 6, 0.0);
            dd.set(4, 7, s*(a[2]*b[1] - a[1]*b[2])/2.0);
            dd.set(4, 8, s*(-(a[3]*b[2]) + a[2]*b[3])/tsq2);

            dd.set(5, 0, 0.0);
            dd.set(5, 1, 0.0);
            dd.set(5, 2, 0.0);
            dd.set(5, 3, 0.0);
            dd.set(5, 4, s*(a[3]*b[2] + a[2]*b[3])/tsq2);
            dd.set(5, 5, s*(a[2]*b[0] + a[0]*b[2])/2.0);
            dd.set(5, 6, 0.0);
            dd.set(5, 7, s*(-(a[3]*b[2]) + a[2]*b[3])/tsq2);
            dd.set(5, 8, s*(a[2]*b[0] - a[0]*b[2])/2.0);

            dd.set(6, 0, s*(-(a[3]*b[0]) + a[0]*b[3])/2.0);
            dd.set(6, 1, s*(a[3]*b[1] - a[1]*b[3])/2.0);
            dd.set(6, 2, 0.0);
            dd.set(6, 3, s*(-(a[1]*b[0]) + a[0]*b[1])/2.0);
            dd.set(6, 4, 0.0);
            dd.set(6, 5, 0.0);
            dd.set(6, 6, s*(-(a[1]*b[0]) - a[0]*b[1] + a[3]*b[3])/2.0);
            dd.set(6, 7, 0.0);
            dd.set(6, 8, 0.0);

            dd.set(7, 0, 0.0);
            dd.set(7, 1, 0.0);
            dd.set(7, 2, 0.0);
            dd.set(7, 3, 0.0);
            dd.set(7, 4, s*(-(a[2]*b[1]) + a[1]*b[2])/2.0);
            dd.set(7, 5, s*(a[3]*b[2] - a[2]*b[3])/tsq2);
            dd.set(7, 6, 0.0);
            dd.set(7, 7, s*(-(a[2]*b[1]) - a[1]*b[2])/2.0);
            dd.set(7, 8, s*(-(a[3]*b[2] + a[2]*b[3])/tsq2));

            dd.set(8, 0, 0.0);
            dd.set(8, 1, 0.0);
            dd.set(8, 2, 0.0);
            dd.set(8, 3, 0.0);
            dd.set(8, 4, s*(a[3]*b[2] - a[2]*b[3])/tsq2);
            dd.set(8, 5, s*(-(a[2]*b[0]) + a[0]*b[2])/2.0);
            dd.set(8, 6, 0.0);
            dd.set(8, 7, s*(-(a[3]*b[2] + a[2]*b[3])/tsq2));
            dd.set(8, 8, s*(-(a[2]*b[0]) - a[0]*b[2])/2.0);
        } else if N == 6 {
            dd.set(0, 0, s*a[0]*b[0]);
            dd.set(0, 1, s*(a[3]*b[3])/2.0);
            dd.set(0, 2, s*(a[5]*b[5])/2.0);
            dd.set(0, 3, s*(a[3]*b[0] + a[0]*b[3])/2.0);
            dd.set(0, 4, s*(a[5]*b[3] + a[3]*b[5])/tsq2);
            dd.set(0, 5, s*(a[5]*b[0] + a[0]*b[5])/2.0);
            dd.set(0, 6, s*(a[3]*b[0] - a[0]*b[3])/2.0);
            dd.set(0, 7, s*(a[5]*b[3] - a[3]*b[5])/tsq2);
            dd.set(0, 8, s*(a[5]*b[0] - a[0]*b[5])/2.0);

            dd.set(1, 0, s*(a[3]*b[3])/2.0);
            dd.set(1, 1, s*a[1]*b[1]);
            dd.set(1, 2, s*(a[4]*b[4])/2.0);
            dd.set(1, 3, s*(a[3]*b[1] + a[1]*b[3])/2.0);
            dd.set(1, 4, s*(a[4]*b[1] + a[1]*b[4])/2.0);
            dd.set(1, 5, s*(a[4]*b[3] + a[3]*b[4])/tsq2);
            dd.set(1, 6, s*(-(a[3]*b[1]) + a[1]*b[3])/2.0);
            dd.set(1, 7, s*(a[4]*b[1] - a[1]*b[4])/2.0);
            dd.set(1, 8, s*(a[4]*b[3] - a[3]*b[4])/tsq2);

            dd.set(2, 0, s*(a[5]*b[5])/2.0);
            dd.set(2, 1, s*(a[4]*b[4])/2.0);
            dd.set(2, 2, s*a[2]*b[2]);
            dd.set(2, 3, s*(a[ 5]*b[4] + a[4]*b[5])/tsq2);
            dd.set(2, 4, s*(a[4]*b[2] + a[2]*b[4])/2.0);
            dd.set(2, 5, s*(a[5]*b[2] + a[2]*b[5])/2.0);
            dd.set(2, 6, s*(-(a[5]*b[4]) + a[4]*b[5])/tsq2);
            dd.set(2, 7, s*(-(a[4]*b[2]) + a[2]*b[4])/2.0);
            dd.set(2, 8, s*(-(a[5]*b[2]) + a[2]*b[5])/2.0);

            dd.set(3, 0, s*(a[3]*b[0] + a[0]*b[3])/2.0);
            dd.set(3, 1, s*(a[3]*b[1] + a[1]*b[3])/2.0);
            dd.set(3, 2, s*(a[5]*b[4] + a[4]*b[5])/tsq2);
            dd.set(3, 3, s*(a[1]*b[0] + a[0]*b[1] + a[3]*b[3])/2.0);
            dd.set(3, 4, s*(SQRT_2*a[5]*b[1] + a[4]*b[3] + a[3]*b[4] + SQRT_2*a[1]*b[5])/4.0);
            dd.set(3, 5, s*(SQRT_2*a[4]*b[0] + a[5]*b[3] + SQRT_2*a[0]*b[4] + a[3]*b[5])/4.0);
            dd.set(3, 6, s*(a[1]*b[0] - a[0]*b[1])/2.0);
            dd.set(3, 7, s*(SQRT_2*a[5]*b[1] + a[4]*b[3] - a[3]*b[4] - SQRT_2*a[1]*b[5])/4.0);
            dd.set(3, 8, s*(SQRT_2*a[4]*b[0] + a[5]*b[3] - SQRT_2*a[0]*b[4] - a[3]*b[5])/4.0);

            dd.set(4, 0, s*(a[5]*b[3] + a[3]*b[5])/tsq2);
            dd.set(4, 1, s*(a[4]*b[1] + a[1]*b[4])/2.0);
            dd.set(4, 2, s*(a[4]*b[2] + a[2]*b[4])/2.0);
            dd.set(4, 3, s*(SQRT_2*a[5]*b[1] + a[4]*b[3] + a[3]*b[4] + SQRT_2*a[1]*b[5])/4.0);
            dd.set(4, 4, s*(a[2]*b[1] + a[1]*b[2] + a[4]*b[4])/2.0);
            dd.set(4, 5, s*(SQRT_2*a[3]*b[2] + SQRT_2*a[2]*b[3] + a[5]*b[4] + a[4]*b[5])/4.0);
            dd.set(4, 6, s*(-(SQRT_2*a[5]*b[1]) + a[4]*b[3] - a[3]*b[4] + SQRT_2*a[1]*b[5])/4.0);
            dd.set(4, 7, s*(a[2]*b[1] - a[1]*b[2])/2.0);
            dd.set(4, 8, s*(-(SQRT_2*a[3]*b[2]) + SQRT_2*a[2]*b[3] - a[5]*b[4] + a[4]*b[5])/4.0);

            dd.set(5, 0, s*(a[5]*b[0] + a[0]*b[5])/2.0);
            dd.set(5, 1, s*(a[4]*b[3] + a[3]*b[4])/tsq2);
            dd.set(5, 2, s*(a[5]*b[2] + a[2]*b[5])/2.0);
            dd.set(5, 3, s*(SQRT_2*a[4]*b[0] + a[5]*b[3] + SQRT_2*a[0]*b[4] + a[3]*b[5])/4.0);
            dd.set(5, 4, s*(SQRT_2*a[3]*b[2] + SQRT_2*a[2]*b[3] + a[5]*b[4] + a[4]*b[5])/4.0);
            dd.set(5, 5, s*(a[2]*b[0] + a[0]*b[2] + a[5]*b[5])/2.0);
            dd.set(5, 6, s*(SQRT_2*a[4]*b[0] - a[5]*b[3] - SQRT_2*a[0]*b[4] + a[3]*b[5])/4.0);
            dd.set(5, 7, s*(-(SQRT_2*a[3]*b[2]) + SQRT_2*a[2]*b[3] + a[5]*b[4] - a[4]*b[5])/4.0);
            dd.set(5, 8, s*(a[2]*b[0] - a[0]*b[2])/2.0);

            dd.set(6, 0, s*(-(a[3]*b[0]) + a[0]*b[3])/2.0);
            dd.set(6, 1, s*(a[3]*b[1] - a[1]*b[3])/2.0);
            dd.set(6, 2, s*(a[5]*b[4] - a[4]*b[5])/tsq2);
            dd.set(6, 3, s*(-(a[1]*b[0]) + a[0]*b[1])/2.0);
            dd.set(6, 4, s*(SQRT_2*a[5]*b[1] - a[4]*b[3] + a[3]*b[4] - SQRT_2*a[1]*b[5])/4.0);
            dd.set(6, 5, s*(-(SQRT_2*a[4]*b[0]) + a[5]*b[3] + SQRT_2*a[0]*b[4] - a[3]*b[5])/4.0);
            dd.set(6, 6, s*(-(a[1]*b[0]) - a[0]*b[1] + a[3]*b[3])/2.0);
            dd.set(6, 7, s*(SQRT_2*a[5]*b[1] - a[4]*b[3] - a[3]*b[4] + SQRT_2*a[1]*b[5])/4.0);
            dd.set(6, 8, s*(-(SQRT_2*a[4]*b[0]) + a[5]*b[3] - SQRT_2*a[0]*b[4] + a[3]*b[5])/4.0);

            dd.set(7, 0, s*(-(a[5]*b[3]) + a[3]*b[5])/tsq2);
            dd.set(7, 1, s*(-(a[4]*b[1]) + a[1]*b[4])/2.0);
            dd.set(7, 2, s*(a[4]*b[2] - a[2]*b[4])/2.0);
            dd.set(7, 3, s*(-(SQRT_2*a[5]*b[1]) - a[4]*b[3] + a[3]*b[4] + SQRT_2*a[1]*b[5])/4.0);
            dd.set(7, 4, s*(-(a[2]*b[1]) + a[1]*b[2])/2.0);
            dd.set(7, 5, s*(SQRT_2*a[3]*b[2] - SQRT_2*a[2]*b[3] - a[5]*b[4] + a[4]*b[5])/4.0);
            dd.set(7, 6, s*(SQRT_2*a[5]*b[1] - a[4]*b[3] - a[3]*b[4] + SQRT_2*a[1]*b[5])/4.0);
            dd.set(7, 7, s*(-(a[2]*b[1]) - a[1]*b[2] + a[4]*b[4])/2.0);
            dd.set(7, 8, s*(-(SQRT_2*a[3]*b[2]) - SQRT_2*a[2]*b[3] + a[5]*b[4] + a[4]*b[5])/4.0);

            dd.set(8, 0, s*(-(a[5]*b[0]) + a[0]*b[5])/2.0);
            dd.set(8, 1, s*(-(a[4]*b[3]) + a[3]*b[4])/tsq2);
            dd.set(8, 2, s*(a[5]*b[2] - a[2]*b[5])/2.0);
            dd.set(8, 3, s*(-(SQRT_2*a[4]*b[0]) - a[5]*b[3] + SQRT_2*a[0]*b[4] + a[3]*b[5])/4.0);
            dd.set(8, 4, s*(SQRT_2*a[3]*b[2] - SQRT_2*a[2]*b[3] + a[5]*b[4] - a[4]*b[5])/4.0);
            dd.set(8, 5, s*(-(a[2]*b[0]) + a[0]*b[2])/2.0);
            dd.set(8, 6, s*(-(SQRT_2*a[4]*b[0]) + a[5]*b[3] - SQRT_2*a[0]*b[4] + a[3]*b[5])/4.0);
            dd.set(8, 7, s*(-(SQRT_2*a[3]*b[2]) - SQRT_2*a[2]*b[3] + a[5]*b[4] + a[4]*b[5])/4.0);
            dd.set(8, 8, s*(-(a[2]*b[0]) - a[0]*b[2] + a[5]*b[5])/2.0);
        } else {
            debug_assert!(N == 9);

            dd.set(0, 0, s*a[0]*b[0]);
            dd.set(0, 1, s*((a[3] + a[6])*(b[3] + b[6]))/2.0);
            dd.set(0, 2, s*((a[5] + a[8])*(b[5] + b[8]))/2.0);
            dd.set(0, 3, s*(a[3]*b[0] + a[6]*b[0] + a[0]*(b[3] + b[6]))/2.0);
            dd.set(0, 4, s*((a[5] + a[8])*(b[3] + b[6]) + (a[3] + a[6])*(b[5] + b[8]))/tsq2);
            dd.set(0, 5, s*(a[5]*b[0] + a[8]*b[0] + a[0]*(b[5] + b[8]))/2.0);
            dd.set(0, 6, s*(a[3]*b[0] + a[6]*b[0] - a[0]*(b[3] + b[6]))/2.0);
            dd.set(0, 7, s*((a[5] + a[8])*(b[3] + b[6]) - (a[3] + a[6])*(b[5] + b[8]))/tsq2);
            dd.set(0, 8, s*(a[5]*b[0] + a[8]*b[0] - a[0]*(b[5] + b[8]))/2.0);

            dd.set(1, 0, s*((a[3] - a[6])*(b[3] - b[6]))/2.0);
            dd.set(1, 1, s*a[1]*b[1]);
            dd.set(1, 2, s*((a[4] + a[7])*(b[4] + b[7]))/2.0);
            dd.set(1, 3, s*(a[3]*b[1] - a[6]*b[1] + a[1]*(b[3] - b[6]))/2.0);
            dd.set(1, 4, s*(a[4]*b[1] + a[7]*b[1] + a[1]*(b[4] + b[7]))/2.0);
            dd.set(1, 5, s*((a[4] + a[7])*(b[3] - b[6]) + (a[3] - a[6])*(b[4] + b[7]))/tsq2);
            dd.set(1, 6, s*(-(a[3]*b[1]) + a[6]*b[1] + a[1]*(b[3] - b[6]))/2.0);
            dd.set(1, 7, s*(a[4]*b[1] + a[7]*b[1] - a[1]*(b[4] + b[7]))/2.0);
            dd.set(1, 8, s*((a[4] + a[7])*(b[3] - b[6]) - (a[3] - a[6])*(b[4] + b[7]))/tsq2);

            dd.set(2, 0, s*((a[5] - a[8])*(b[5] - b[8]))/2.0);
            dd.set(2, 1, s*((a[4] - a[7])*(b[4] - b[7]))/2.0);
            dd.set(2, 2, s*a[2]*b[2]);
            dd.set(2, 3, s*((a[5] - a[8])*(b[4] - b[7]) + (a[4] - a[7])*(b[5] - b[8]))/tsq2);
            dd.set(2, 4, s*(a[4]*b[2] - a[7]*b[2] + a[2]*(b[4] - b[7]))/2.0);
            dd.set(2, 5, s*(a[5]*b[2] - a[8]*b[2] + a[2]*(b[5] - b[8]))/2.0);
            dd.set(2, 6, s*(-((a[5] - a[8])*(b[4] - b[7])) + (a[4] - a[7])*(b[5] - b[8]))/tsq2);
            dd.set(2, 7, s*(-(a[4]*b[2]) + a[7]*b[2] + a[2]*(b[4] - b[7]))/2.0);
            dd.set(2, 8, s*(-(a[5]*b[2]) + a[8]*b[2] + a[2]*(b[5] - b[8]))/2.0);

            dd.set(3, 0, s*(a[3]*b[0] - a[6]*b[0] + a[0]*(b[3] - b[6]))/2.0);
            dd.set(3, 1, s*(a[3]*b[1] + a[6]*b[1] + a[1]*(b[3] + b[6]))/2.0);
            dd.set(3, 2, s*((a[5] + a[8])*(b[4] + b[7]) + (a[4] + a[7])*(b[5] + b[8]))/tsq2);
            dd.set(3, 3, s*(a[1]*b[0] + a[0]*b[1] + a[3]*b[3] - a[6]*b[6])/2.0);
            dd.set(3, 4, s*(SQRT_2*(a[5] + a[8])*b[1] + (a[4] + a[7])*(b[3] + b[6]) + (a[3] + a[6])*(b[4] + b[7]) + SQRT_2*a[1]*(b[5] + b[8]))/4.0);
            dd.set(3, 5, s*(SQRT_2*(a[4] + a[7])*b[0] + (a[5] + a[8])*(b[3] - b[6]) + SQRT_2*a[0]*(b[4] + b[7]) + (a[3] - a[6])*(b[5] + b[8]))/4.0);
            dd.set(3, 6, s*(a[1]*b[0] - a[0]*b[1] + a[6]*b[3] - a[3]*b[6])/2.0);
            dd.set(3, 7, s*(SQRT_2*(a[5] + a[8])*b[1] + (a[4] + a[7])*(b[3] + b[6]) - (a[3] + a[6])*(b[4] + b[7]) - SQRT_2*a[1]*(b[5] + b[8]))/4.0);
            dd.set(3, 8, s*(SQRT_2*(a[4] + a[7])*b[0] + (a[5] + a[8])*(b[3] - b[6]) - SQRT_2*a[0]*(b[4] + b[7]) - (a[3] - a[6])*(b[5] + b[8]))/4.0);

            dd.set(4, 0, s*((a[5] - a[8])*(b[3] - b[6]) + (a[3] - a[6])*(b[5] - b[8]))/tsq2);
            dd.set(4, 1, s*(a[4]*b[1] - a[7]*b[1] + a[1]*(b[4] - b[7]))/2.0);
            dd.set(4, 2, s*(a[4]*b[2] + a[7]*b[2] + a[2]*(b[4] + b[7]))/2.0);
            dd.set(4, 3, s*(SQRT_2*(a[5] - a[8])*b[1] + (a[4] - a[7])*(b[3] - b[6]) + (a[3] - a[6])*(b[4] - b[7]) + SQRT_2*a[1]*(b[5] - b[8]))/4.0);
            dd.set(4, 4, s*(a[2]*b[1] + a[1]*b[2] + a[4]*b[4] - a[7]*b[7])/2.0);
            dd.set(4, 5, s*(SQRT_2*(a[3] - a[6])*b[2] + SQRT_2*a[2]*(b[3] - b[6]) + (a[5] - a[8])*(b[4] + b[7]) + (a[4] + a[7])*(b[5] - b[8]))/4.0);
            dd.set(4, 6, s*(-(SQRT_2*(a[5] - a[8])*b[1]) + (a[4] - a[7])*(b[3] - b[6]) - (a[3] - a[6])*(b[4] - b[7]) + SQRT_2*a[1]*(b[5] - b[8]))/4.0);
            dd.set(4, 7, s*(a[2]*b[1] - a[1]*b[2] + a[7]*b[4] - a[4]*b[7])/2.0);
            dd.set(4, 8, s*(-(SQRT_2*(a[3] - a[6])*b[2]) + SQRT_2*a[2]*(b[3] - b[6]) - (a[5] - a[8])*(b[4] + b[7]) + (a[4] + a[7])*(b[5] - b[8]))/4.0);

            dd.set(5, 0, s*(a[5]*b[0] - a[8]*b[0] + a[0]*(b[5] - b[8]))/2.0);
            dd.set(5, 1, s*((a[4] - a[7])*(b[3] + b[6]) + (a[3] + a[6])*(b[4] - b[7]))/tsq2);
            dd.set(5, 2, s*(a[5]*b[2] + a[8]*b[2] + a[2]*(b[5] + b[8]))/2.0);
            dd.set(5, 3, s*(SQRT_2*(a[4] - a[7])*b[0] + (a[5] - a[8])*(b[3] + b[6]) + SQRT_2*a[0]*(b[4] - b[7]) + (a[3] + a[6])*(b[5] - b[8]))/4.0);
            dd.set(5, 4, s*(SQRT_2*(a[3] + a[6])*b[2] + SQRT_2*a[2]*(b[3] + b[6]) + (a[5] + a[8])*(b[4] - b[7]) + (a[4] - a[7])*(b[5] + b[8]))/4.0);
            dd.set(5, 5, s*(a[2]*b[0] + a[0]*b[2] + a[5]*b[5] - a[8]*b[8])/2.0);
            dd.set(5, 6, s*(SQRT_2*(a[4] - a[7])*b[0] - (a[5] - a[8])*(b[3] + b[6]) - SQRT_2*a[0]*(b[4] - b[7]) + (a[3] + a[6])*(b[5] - b[8]))/4.0);
            dd.set(5, 7, s*(-(SQRT_2*(a[3] + a[6])*b[2]) + SQRT_2*a[2]*(b[3] + b[6]) + (a[5] + a[8])*(b[4] - b[7]) - (a[4] - a[7])*(b[5] + b[8]))/4.0);
            dd.set(5, 8, s*(a[2]*b[0] - a[0]*b[2] + a[8]*b[5] - a[5]*b[8])/2.0);

            dd.set(6, 0, s*(-(a[3]*b[0]) + a[6]*b[0] + a[0]*(b[3] - b[6]))/2.0);
            dd.set(6, 1, s*(a[3]*b[1] + a[6]*b[1] - a[1]*(b[3] + b[6]))/2.0);
            dd.set(6, 2, s*((a[5] + a[8])*(b[4] + b[7]) - (a[4] + a[7])*(b[5] + b[8]))/tsq2);
            dd.set(6, 3, s*(-(a[1]*b[0]) + a[0]*b[1] + a[6]*b[3] - a[3]*b[6])/2.0);
            dd.set(6, 4, s*(SQRT_2*(a[5] + a[8])*b[1] - (a[4] + a[7])*(b[3] + b[6]) + (a[3] + a[6])*(b[4] + b[7]) - SQRT_2*a[1]*(b[5] + b[8]))/4.0);
            dd.set(6, 5, s*(-(SQRT_2*(a[4] + a[7])*b[0]) + (a[5] + a[8])*(b[3] - b[6]) + SQRT_2*a[0]*(b[4] + b[7]) - (a[3] - a[6])*(b[5] + b[8]))/4.0);
            dd.set(6, 6, s*(-(a[1]*b[0]) - a[0]*b[1] + a[3]*b[3] - a[6]*b[6])/2.0);
            dd.set(6, 7, s*(SQRT_2*(a[5] + a[8])*b[1] - (a[4] + a[7])*(b[3] + b[6]) - (a[3] + a[6])*(b[4] + b[7]) + SQRT_2*a[1]*(b[5] + b[8]))/4.0);
            dd.set(6, 8, s*(-(SQRT_2*(a[4] + a[7])*b[0]) + (a[5] + a[8])*(b[3] - b[6]) - SQRT_2*a[0]*(b[4] + b[7]) + (a[3] - a[6])*(b[5] + b[8]))/4.0);

            dd.set(7, 0, s*(-((a[5] - a[8])*(b[3] - b[6])) + (a[3] - a[6])*(b[5] - b[8]))/tsq2);
            dd.set(7, 1, s*(-(a[4]*b[1]) + a[7]*b[1] + a[1]*(b[4] - b[7]))/2.0);
            dd.set(7, 2, s*(a[4]*b[2] + a[7]*b[2] - a[2]*(b[4] + b[7]))/2.0);
            dd.set(7, 3, s*(-(SQRT_2*(a[5] - a[8])*b[1]) - (a[4] - a[7])*(b[3] - b[6]) + (a[3] - a[6])*(b[4] - b[7]) + SQRT_2*a[1]*(b[5] - b[8]))/4.0);
            dd.set(7, 4, s*(-(a[2]*b[1]) + a[1]*b[2] + a[7]*b[4] - a[4]*b[7])/2.0);
            dd.set(7, 5, s*(SQRT_2*(a[3] - a[6])*b[2] - SQRT_2*a[2]*(b[3] - b[6]) - (a[5] - a[8])*(b[4] + b[7]) + (a[4] + a[7])*(b[5] - b[8]))/4.0);
            dd.set(7, 6, s*(SQRT_2*(a[5] - a[8])*b[1] - (a[4] - a[7])*(b[3] - b[6]) - (a[3] - a[6])*(b[4] - b[7]) + SQRT_2*a[1]*(b[5] - b[8]))/4.0);
            dd.set(7, 7, s*(-(a[2]*b[1]) - a[1]*b[2] + a[4]*b[4] - a[7]*b[7])/2.0);
            dd.set(7, 8, s*(-(SQRT_2*(a[3] - a[6])*b[2]) - SQRT_2*a[2]*(b[3] - b[6]) + (a[5] - a[8])*(b[4] + b[7]) + (a[4] + a[7])*(b[5] - b[8]))/4.0);

            dd.set(8, 0, s*(-(a[5]*b[0]) + a[8]*b[0] + a[0]*(b[5] - b[8]))/2.0);
            dd.set(8, 1, s*(-((a[4] - a[7])*(b[3] + b[6])) + (a[3] + a[6])*(b[4] - b[7]))/tsq2);
            dd.set(8, 2, s*(a[5]*b[2] + a[8]*b[2] - a[2]*(b[5] + b[8]))/2.0);
            dd.set(8, 3, s*(-(SQRT_2*(a[4] - a[7])*b[0]) - (a[5] - a[8])*(b[3] + b[6]) + SQRT_2*a[0]*(b[4] - b[7]) + (a[3] + a[6])*(b[5] - b[8]))/4.0);
            dd.set(8, 4, s*(SQRT_2*(a[3] + a[6])*b[2] - SQRT_2*a[2]*(b[3] + b[6]) + (a[5] + a[8])*(b[4] - b[7]) - (a[4] - a[7])*(b[5] + b[8]))/4.0);
            dd.set(8, 5, s*(-(a[2]*b[0]) + a[0]*b[2] + a[8]*b[5] - a[5]*b[8])/2.0);
            dd.set(8, 6, s*(-(SQRT_2*(a[4] - a[7])*b[0]) + (a[5] - a[8])*(b[3] + b[6]) - SQRT_2*a[0]*(b[4] - b[7]) + (a[3] + a[6])*(b[5] - b[8]))/4.0);
            dd.set(8, 7, s*(-(SQRT_2*(a[3] + a[6])*b[2]) - SQRT_2*a[2]*(b[3] + b[6]) + (a[5] + a[8])*(b[4] - b[7]) + (a[4] - a[7])*(b[5] + b[8]))/4.0);
            dd.set(8, 8, s*(-(a[2]*b[0]) - a[0]*b[2] + a[5]*b[5] - a[8]*b[8])/2.0);
        }
    }
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

#[cfg(test)]
mod tests {
    use super::t2_udyad_t2;
    use crate::{ADD, MN_TO_IJKL, SET, SQRT_2};
    use crate::{Tensor2, Tensor4};
    use russell_lab::{Matrix, mat_approx_eq};

    fn kelvin_matrix<const N: usize>(dd: &Tensor4<N>) -> Matrix {
        let mut m = Matrix::new(N, N);
        for i in 0..N {
            for j in 0..N {
                m.set(i, j, dd.get(i, j));
            }
        }
        m
    }

    fn check_udyad<const N: usize>(s: f64, a_ten: &Tensor2<N>, b_ten: &Tensor2<N>, dd_ten: &Tensor4<9>, tol: f64) {
        let a = a_ten.as_std_matrix();
        let b = b_ten.as_std_matrix();
        let dd = dd_ten.as_std_matrix();
        let mut correct = Matrix::new(9, 9); // Use 9 here due to the conversion to "STD"
        for m in 0..9 {
            for n in 0..9 {
                let (i, j, k, l) = MN_TO_IJKL[m][n];
                correct.set(m, n, s * a.get(i, l) * b.get(j, k));
            }
        }
        mat_approx_eq(&dd, &correct, tol);
    }

    #[test]
    fn t2_udyad_t2_works() {
        // general udyad general
        #[rustfmt::skip]
        let a = Tensor2::<9>::from_std_matrix(&[
            [1.0, 2.0, 3.0],
            [4.0, 5.0, 6.0],
            [7.0, 8.0, 9.0],
        ]).unwrap();
        #[rustfmt::skip]
        let b = Tensor2::<9>::from_std_matrix(&[
            [9.0, 8.0, 7.0],
            [6.0, 5.0, 4.0],
            [3.0, 2.0, 1.0],
        ]).unwrap();
        let mut dd = Tensor4::<9>::new();
        t2_udyad_t2(&mut dd, SET, 2.0, &a, &b);
        let mat = dd.as_std_matrix();
        let correct = Matrix::from(&[
            [18.0, 32.0, 42.0, 36.0, 48.0, 54.0, 16.0, 28.0, 14.0],
            [48.0, 50.0, 48.0, 60.0, 60.0, 72.0, 40.0, 40.0, 32.0],
            [42.0, 32.0, 18.0, 48.0, 36.0, 54.0, 28.0, 16.0, 14.0],
            [12.0, 20.0, 24.0, 24.0, 30.0, 36.0, 10.0, 16.0, 8.0],
            [24.0, 20.0, 12.0, 30.0, 24.0, 36.0, 16.0, 10.0, 8.0],
            [6.0, 8.0, 6.0, 12.0, 12.0, 18.0, 4.0, 4.0, 2.0],
            [72.0, 80.0, 84.0, 90.0, 96.0, 108.0, 64.0, 70.0, 56.0],
            [84.0, 80.0, 72.0, 96.0, 90.0, 108.0, 70.0, 64.0, 56.0],
            [126.0, 128.0, 126.0, 144.0, 144.0, 162.0, 112.0, 112.0, 98.0],
        ]);
        mat_approx_eq(&mat, &correct, 1e-13);
        check_udyad(2.0, &a, &b, &dd, 1e-13);

        // symmetric udyad symmetric
        #[rustfmt::skip]
        let a = Tensor2::<6>::from_std_matrix(&[
            [1.0, 4.0, 6.0],
            [4.0, 2.0, 5.0],
            [6.0, 5.0, 3.0],
        ]).unwrap();
        #[rustfmt::skip]
        let b = Tensor2::<6>::from_std_matrix(&[
            [3.0, 5.0, 6.0],
            [5.0, 2.0, 4.0],
            [6.0, 4.0, 1.0],
        ]).unwrap();
        let mut dd = Tensor4::<9>::new();
        t2_udyad_t2(&mut dd, SET, 2.0, &a, &b);
        let mat = dd.as_std_matrix();
        let correct = Matrix::from(&[
            [6.0, 40.0, 72.0, 24.0, 60.0, 36.0, 10.0, 48.0, 12.0],
            [40.0, 8.0, 40.0, 20.0, 20.0, 50.0, 16.0, 16.0, 32.0],
            [72.0, 40.0, 6.0, 60.0, 24.0, 36.0, 48.0, 10.0, 12.0],
            [10.0, 16.0, 48.0, 40.0, 24.0, 60.0, 4.0, 32.0, 8.0],
            [48.0, 16.0, 10.0, 24.0, 40.0, 60.0, 32.0, 4.0, 8.0],
            [12.0, 32.0, 12.0, 48.0, 48.0, 72.0, 8.0, 8.0, 2.0],
            [24.0, 20.0, 60.0, 12.0, 50.0, 30.0, 40.0, 24.0, 48.0],
            [60.0, 20.0, 24.0, 50.0, 12.0, 30.0, 24.0, 40.0, 48.0],
            [36.0, 50.0, 36.0, 30.0, 30.0, 18.0, 60.0, 60.0, 72.0],
        ]);
        mat_approx_eq(&mat, &correct, 1e-13);
        check_udyad(2.0, &a, &b, &dd, 1e-13);

        // symmetric generalized plane udyad symmetric generalized plane
        #[rustfmt::skip]
        let a = Tensor2::<4>::from_std_matrix(&[
            [1.0, 4.0, 0.0],
            [4.0, 2.0, 0.0],
            [0.0, 0.0, 3.0],
        ]).unwrap();
        #[rustfmt::skip]
        let b = Tensor2::<4>::from_std_matrix(&[
            [3.0, 4.0, 0.0],
            [4.0, 2.0, 0.0],
            [0.0, 0.0, 1.0],
        ]).unwrap();
        let mut dd = Tensor4::<9>::new();
        t2_udyad_t2(&mut dd, SET, 2.0, &a, &b);
        let kelvin_mat = Matrix::from(&[
            [6.0, 32.0, 0.0, 16.0 * SQRT_2, 0.0, 0.0, 8.0 * SQRT_2, 0.0, 0.0],
            [32.0, 8.0, 0.0, 16.0 * SQRT_2, 0.0, 0.0, 0.0, 0.0, 0.0],
            [0.0, 0.0, 6.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
            [16.0 * SQRT_2, 16.0 * SQRT_2, 0.0, 40.0, 0.0, 0.0, 4.0, 0.0, 0.0],
            [0.0, 0.0, 0.0, 0.0, 8.0, 16.0, 0.0, 4.0, 8.0],
            [0.0, 0.0, 0.0, 0.0, 16.0, 10.0, 0.0, 8.0, 8.0],
            [-8.0 * SQRT_2, 0.0, 0.0, -4.0, 0.0, 0.0, 24.0, 0.0, 0.0],
            [0.0, 0.0, 0.0, 0.0, -4.0, -8.0, 0.0, -8.0, -16.0],
            [0.0, 0.0, 0.0, 0.0, -8.0, -8.0, 0.0, -16.0, -10.0],
        ]);
        mat_approx_eq(&kelvin_matrix(&dd), &kelvin_mat, 1e-14);
        let mat = dd.as_std_matrix();
        let correct = Matrix::from(&[
            [6.0, 32.0, 0.0, 24.0, 0.0, 0.0, 8.0, 0.0, 0.0],
            [32.0, 8.0, 0.0, 16.0, 0.0, 0.0, 16.0, 0.0, 0.0],
            [0.0, 0.0, 6.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
            [8.0, 16.0, 0.0, 32.0, 0.0, 0.0, 4.0, 0.0, 0.0],
            [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 8.0],
            [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 8.0, 2.0],
            [24.0, 16.0, 0.0, 12.0, 0.0, 0.0, 32.0, 0.0, 0.0],
            [0.0, 0.0, 0.0, 0.0, 12.0, 24.0, 0.0, 0.0, 0.0],
            [0.0, 0.0, 0.0, 0.0, 24.0, 18.0, 0.0, 0.0, 0.0],
        ]);
        mat_approx_eq(&mat, &correct, 1e-14);
        check_udyad(2.0, &a, &b, &dd, 1e-15);
    }

    #[test]
    fn t2_udyad_t2_add_works() {
        // general udyad general
        #[rustfmt::skip]
        let a = Tensor2::<9>::from_std_matrix(&[
            [1.0, 2.0, 3.0],
            [4.0, 5.0, 6.0],
            [7.0, 8.0, 9.0],
        ]).unwrap();
        #[rustfmt::skip]
        let b = Tensor2::<9>::from_std_matrix(&[
            [9.0, 8.0, 7.0],
            [6.0, 5.0, 4.0],
            [3.0, 2.0, 1.0],
        ]).unwrap();
        let mut dd = Tensor4::<9>::new();
        t2_udyad_t2(&mut dd, SET, 2.0, &a, &b);
        t2_udyad_t2(&mut dd, ADD, 3.0, &a, &b);
        check_udyad(5.0, &a, &b, &dd, 1e-12);

        // symmetric udyad symmetric
        #[rustfmt::skip]
        let a = Tensor2::<6>::from_std_matrix(&[
            [1.0, 4.0, 6.0],
            [4.0, 2.0, 5.0],
            [6.0, 5.0, 3.0],
        ]).unwrap();
        #[rustfmt::skip]
        let b = Tensor2::<6>::from_std_matrix(&[
            [3.0, 5.0, 6.0],
            [5.0, 2.0, 4.0],
            [6.0, 4.0, 1.0],
        ]).unwrap();
        let mut dd = Tensor4::<9>::new();
        t2_udyad_t2(&mut dd, SET, 2.0, &a, &b);
        t2_udyad_t2(&mut dd, ADD, 3.0, &a, &b);
        check_udyad(5.0, &a, &b, &dd, 1e-12);

        // symmetric generalized plane udyad symmetric generalized plane
        #[rustfmt::skip]
        let a = Tensor2::<4>::from_std_matrix(&[
            [1.0, 4.0, 0.0],
            [4.0, 2.0, 0.0],
            [0.0, 0.0, 3.0],
        ]).unwrap();
        #[rustfmt::skip]
        let b = Tensor2::<4>::from_std_matrix(&[
            [3.0, 4.0, 0.0],
            [4.0, 2.0, 0.0],
            [0.0, 0.0, 1.0],
        ]).unwrap();
        let mut dd = Tensor4::<9>::new();
        t2_udyad_t2(&mut dd, SET, 2.0, &a, &b);
        t2_udyad_t2(&mut dd, ADD, 3.0, &a, &b);
        check_udyad(5.0, &a, &b, &dd, 1e-12);
    }
}