A minor mistake by me, sorry
. The bound depends on rounding mode and should be
| n * phi^(n-1) * epsilon / sqrt(5) | < 1/2 // assume round to nearest.
For rounding up, down and towards zero, the bound is 1 but we need to be careful about the sign of epsilon. For example, we must guarantee that the error in phi >= 0 if we are to rounding down or towards 0.
And yes, I would do 2). If we compute sqrt() using Newton's method, we can fold this check into sqrt.