The issue is due to rounding during double -> float and float->double conversions. This algorithm modifies slightly its input seismic line at each iteration, causing the computed normal to be slightly different each time because of this rounding issue. Converting the normalized normal to float first increases a bit the precision of the algorithm. But the best work around so far, is to store the first normal and use it at each iteration.
Cheers, Priya