Yeah, there's a similar formula for all homogeneous linear recurrences if I'm not mistaken, but you can't evaluate it naively with floating point arithmetic. We can view (1 + sqrt(5)) and (1 - sqrt(5)) as elements of Z[sqrt(5)] though. It's probably extremely annoying to write because you have to implement the arithmetic operations in that field, but it wouldn't surprise me if it had faster runtime.
Probably Claude can do it? Maybe I'll give it a try later :)
Good luck calculating (1+sqrt(5))^n in Z[sqrt(5)].