Basically, something like a**(1.0/3.0) lose precision when a is large because 1.0/3.0 doesn’t have a finite representation in floating point.
Now, Python has implemented a cbrt function that doesn’t incur in the loosing of the precision, and have been suggested to implement a function like root(a, n) that may calculating the n-root without this problem.
Going back to Fortran,
I’m suggesting to implement in Fortran a new intrinsic that we may call powerroot(value, m, n) that may calculate value**(real(m)/real(n)) without loosing precision do to incorrectness of real(m)/real(n) and may avoid overflows.
I’m sorry I don’t know which algorithm one may use, so it is just an half proposal.
Hi, this is why to compute nth roots I use either root() from the mpmath package, or root() from the gmpy2 package. I’m pretty sure that these algorithms are available in some C libraries that can be wrapped.
Fortran stdlib appears the right place to try out such things where the experimental attempts can eventually turn into proposals for the language. So please look there first and consider trying out options there
Expressions like this are evaluated as exp(log(a)*1.0/3.0) or maybe exp(log(a)/3.0)with some optimization. Unless there is something wrong with the exp() and log() evaluations, from where would the error arise?
Of course, efficiency is a separate issue. There may be faster ways to compute specific values, including cube roots, but however they are evaluated they should all agree to the last bit or two, which is all that you can expect from floating point arithmetic. Are there data anywhere that shows otherwise?
See this code, that use “quadruple precision” as a reference
program main
USE ISO_FORTRAN_ENV, only : qp => real128
real(qp) :: a, b
real :: c
do i = 0, 30
a = 10.0_qp ** i
b = a**(1.0_qp/3.0_qp)
c = real(a)**(1.0/3.0)
write(*,*) i, spacing(real(b)), real(c-b)
end do
end program
As a result, the deviation between the single and quadruple precision computation gets larger than the single precision spacing() for values of a larger than 1e3 (compiled with gfortran)
program root
use, intrinsic :: iso_fortran_env, only: real128
real :: a, b, c, db, dc, qb
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)
db = (b - qb) / spacing(b)
dc = (c - qb) / spacing(c)
write(*,'(i2,*(1x,f5.1))') i, db, dc
enddo
end program root
It is interesting that the ** expression differs from the exp(log) expression. Most of the errors are within 2 or 3 bits. The largest error is for 10.0**27, which for some reason has 4 and 5 bits in error.
Just a matter of style, but would it possible to define a binary operator for the Nth root, just like the Nth power is? // would have been the perfect twin of **, but it’s already used for the string concatenation. What about */ ? (x */ n would be the approximation of x ** (1.0/n))
Interesting idea. If we discuss the paper I will mention it. // could be used, as the interface would not conflict with it’s existing use. I’m not sure how much I like it, but we can think about it.
I agree that notation would be misleading for the root operator. In fact, I’ve wondered why ^ has not been introduced into fortran as a substitute for ** decades ago. I’ve always assumed it is being held back and saved for something better than exponentiation.
As for an operator notation for the integer root, something like n\/x could also work, although the combination of backslash and forward slash looks like an uppercase V with some fonts.
You could change + to have documented behavior of multiplication if the second number is 63.1 but that doesn’t make that a good idea. Saying “It’s only misleading if you don’t bother to learn the language” ignores the fact that there are a lot of ways to design a language badly.
If introducing a new notation could be easily manageable by the compilers, why not simply ‘^/’ ? This looks natural and might even pave the way for getting ‘^’ to actually work as a replacement for ‘**’ ?
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?
Ok, two decades then. In any case, this would be like allowing [ and ] to denote arrays in addition to (/ and /), it would be backwards compatible with prior source code.
Here is one way to avoid division, using one Newton-Raphson step to find the reciprocal, and using that reciprocal in another Newton-Raphson step for the cube root.
! Given N, find the cube root, x, and the reciprocal, y
! without using any division ops
!
program cubrt
implicit none
real x,y,N,dx,dy
integer i
N = 27.0
x = 2.0 ! approx cub.rt
y = 0.3 ! approx reciprocal of x
print '(i1,2es13.6)',0,x,y
do i = 1, 10
dy = y * (1.0-x*y)
y = y + dy
dx = (N*y*y - x)*0.333333333
x = x + dx
print '(i1,2es13.6)',i,x,y
if (abs(dx) < 1e-7 .and. abs(dy) < 1e-7)exit
end do
end
A good working solution for the cube root problem would probably require a good way of generating sufficiently close first approximations without resorting to calling exp or log.
Similar issues occur in converting a binary floating point number to a decimal string. A good approximation for the base-10 log of the input number is needed.
Interesting idea. If we discuss the paper I will mention it.