Implement bce_with_logits(logits, targets) that computes the mean binary cross-entropy loss from raw logits in a numerically stable way (log-sum-exp / softplus form), without applying a sigmoid first.

Use the stable per-element formula:

$\ell = \max(x, 0) - x \cdot y + \log(1 + e^{-|x|})$

where $x$ is a logit and $y \in \{0, 1\}$ is the corresponding target. Average over all elements, then return the result as a Python float rounded to 4 decimal places.

logits and targets are 1-D torch.Tensor values of the same shape (CPU, float).