I was simulating results of the function in LibreOffice Calc, and what I discovered, that function is doing : Output = Input x coefficient (1.00273790935) - Input.
With Input 86398 , Outputs = 235. This result was received using rounddown and mod sub-function,
That is pretty close, same time not exact, correct result: 236.5498920213
My suggestion, instead of multistage calculation do simple multiply - shift, which is basically always work as substitute for floating point math. First step would be rearrange formula:
Output = Input x coefficient (1.00273790935) - Input
to Output = Input x coefficient (0.00273790935).
Next, get value 0.00273790935 as long integer, that is easy, just multiply 2^32 :
0.00273790935 x 2^32 = 11759231
Here you are, one line code: Output = ( 11759231UL * Input ) >> 32.
Simulation: 86398 * 11759231 >> 32 = 236, that is better than original math function.
*Not tested, could be casting issue.
Eddited:
There is a mistake in original code comments, valid values for z is not 0<=z<=86 399, but
0<=z<=86 163, so my simulation with 86398 isn't correct.
Plus, after thinking on 32-bit overflow prevention, I came up with another "scaled" down formula with coefficient 0.00273790935 x 2^22 = 11 484
Output = Input + ((Input x 11 484) >> 22).
Check in LibreOffice:
Input Output dz cz bz dz2 cz2 bz1 bz2 n result
86163 86398.906483324 80562405 635796777 86335326 8056 63580 8906 8639 86399 86399
Approximation:
Output = 86163UL + ((86 163UL x 11 484) >> 22) = 86 398
Well, it's "off" with error = 1 ( strictly, error = 0.90......) . Same time it's fastest. To get more accurate results, little bit slower, we need to do "round":
Output = 86163UL + ((86163UL x 11 484 + 2097152UL) >> 22) = 86 399 ,
where 2 097 152 is 2^22 / 2, or "0.5" in float arithmetics.
Finally: Output = Input + ((Input x 11 484 + 2097152) >> 22).