Function carlson_rj 

pub fn carlson_rj<P, E, V>(x: V, y: V, z: V, p: V) -> V
where
    P: Policy,
    E: FloatElement + EllipticConsts,
    V: SpecializedSpecialMath<E>,

Carlson symmetric integral of the third kind, R_J(x, y, z, p), via duplication. Handles p < 0 (a Cauchy principal value) through Carlson’s transform to a positive parameter. Each step accumulates an R_C term, so this is the most expensive Carlson primitive.

Accuracy note: when p coincides with one of x, y, z (so (p-x)(p-y)(p-z) -> 0), the per-step R_C(1, b) term has b -> 1 every iteration. That used to lose ~7 digits, but carlson_rc now switches to its small-argument series there, so this case holds full precision without the R_D special-case or extended precision Boost resorts to.