Function compensated_threshold_distance 

fn compensated_threshold_distance(q: [Complex; 3], k: Complex) -> Complex

q . q - k^2 of a diffraction order to a few ulps of itself, from exact products and compensated additions: next to a threshold the rounded q . q / k^2 - 1 loses about eps k^2 / |k_q^2| of it. Only the zeroth order of an unreduced Bloch vector takes it (within NEAR_THRESHOLD): every other q carries the rounding of its reciprocal vector, as large as the rounding this removes.