Integer root of floating point value

23

I added another value to these results where I did one Newton step in order to refine the values. This also reduces the errors down to at most one bit in the same way that ifort does by default.

program root
   use, intrinsic :: iso_fortran_env, only: real128
   real :: a, b, c, e, db, dc, qb, de
   do i = 0, 30
      a = 10.0 ** i
      b = a**(1.0/3.0)
      qb = real((10.0_real128 ** i) ** (1.0_real128 / 3.0_real128))
      c = exp(log(a)/3.0)
      e = (2.0*c+a/c**2)/3.0   ! refine with a newton step.
      db = (b - qb) / spacing(b)
      dc = (c - qb) / spacing(c)
      de = (e - qb) / spacing(e)
      write(*,'(i2,*(1x,f5.1))') i, db, dc, de
   enddo
 end program root

gfortran root.f90 && a.out
 0   0.0   0.0   0.0
 1   0.0   0.0   0.0
 2   1.0   0.0   0.0
 3   1.0   0.0   0.0
 4   1.0  -1.0   0.0
 5   1.0  -2.0   0.0
 6   2.0   1.0   0.0
 7   2.0  -3.0   1.0
 8   3.0  -1.0   1.0
 9   3.0   2.0   0.0
10   3.0  -1.0   0.0
11   2.0   0.0  -1.0
12   3.0   1.0   0.0
13   3.0   3.0  -1.0
14   4.0  -5.0   0.0
15   4.0  -4.0   0.0
16   5.0  -2.0  -1.0
17   6.0   1.0   0.0
18   7.0   3.0   0.0
19   3.0  -5.0   0.0
20   4.0  -4.0  -1.0
21   5.0  -3.0   0.0
22   6.0   0.0   0.0
23   6.0   1.0   0.0
24   7.0   3.0   0.0
25   8.0   6.0   0.0
26   9.0   9.0   0.0
27  10.0 -18.0   0.0
28   5.0  -9.0   0.0
29   6.0  10.0   0.0
30   7.0  -6.0   0.0

That refinement step requires a division. Maybe there are other refinement expressions that don’t require that. But on the other hand, a division is not too expensive after computing both a log() and an exp() for the general case.

I’m still puzzled by the unusually large errors in the 10.0**27 case. You can do that one in your head, so why does the computer struggle with it?