Function contract_site_tensors 

pub fn contract_site_tensors<T: SVDScalar + EinsumScalar>(
    a: &Tensor4<T>,
    b: &Tensor4<T>,
) -> Result<Tensor4<T>>
where
    <T as ComplexFloat>::Real: Into<f64>,

Contract two 4D site tensors over their shared physical index

Given two 4D tensors:

• A: (left_a, s1_a, s2_a, right_a) where s2_a is the shared index
• B: (left_b, s1_b, s2_b, right_b) where s1_b is the shared index

Contracts over s2_a = s1_b to produce:

• C: (left_a * left_b, s1_a, s2_b, right_a * right_b)

This is the Rust equivalent of _contractsitetensors from Julia.