The following is not exactly right:
part1a = (M - sqrt(M.^2))/2
part1b = part1a./(part1a - delta)
part2a = (M + sqrt(M.^2))/2
part2b = part2a./(part2a + delta)
cost = part1b * costa + part2b * costb
Mathematically it is wrong at exactly two points, M = -eps(realmin) and M = +eps(realmin) . For those two points, the output should be costa and costb respectively, but instead the formula mathematically gives costa/2 and costb/2 at those two points instead.
In practice, though, for values sufficiently close to +/- realmin, numeric evaluation might return costa+costb and for values sufficiently clost to eps(realmin) numeric evaluation might return 0 instead of costa or costb .
The exact result in a range close to +/- realmin is going to depend on the exact order of evaluation, which is not something that we have control over; optimization could potentially re-arrange the evaluation.
This code has been constructed so that it should never return NaN.
How important is it for your purposes that the values must be correct near +/- realmin, given that it is designed to return 0 for 0 exactly?