This is a short code to demonstrate some features of fortran using complex arithmetic.
There is a familiar online puzzle to find solutions to the equation 1^x=4 (or some other right-hand-side value). There are of course no real solutions to this problem. In the online comments, you will invariably see declarations like “1 raised to any power is always 1, so there are no solutions.” However, there are complex solutions because the inverse of complex exponentiation is multivalued, and if you work on different branches you can get different solutions. It turns out that you can find numerical approximations to these solutions in fortran, but you must be careful exactly how the expressions are evaluated. Here is a code that demonstrates this.
program xxx
! print some complex solutions to the equation: 1**x = rhs.
use, intrinsic :: iso_fortran_env, only: wp=>real64
implicit none
complex(wp) :: z1, arg, x
real(wp), parameter :: rhs = 4.0_wp, twopi = 8 * atan(1.0_wp), lnrhs = log(rhs)
integer :: k
character(*), parameter :: cfmth='("Complex solutions to 1**x=",g0.4/a2,*(4x,a11,5x))', &
& cfmt='(i2,*(" (",es0.2,",",es0.2,")":))'
write(*,cfmth) rhs, 'k', 'arg=i2Pik', 'exp(arg)', 'x', 'exp(arg)**x', 'exp(arg*x)'
do k = 1, 10
arg = cmplx( 0.0_wp, twopi * k, kind=wp )
z1 = exp( arg )
x = cmplx( 0.0_wp, -lnrhs/(twopi * k), kind=wp )
write(*,cfmt) k, arg, z1, x, z1**x, exp(arg*x)
enddo
end program xxx
$ flang xxx.f90 && a.out
Complex solutions to 1**x=4.000
k arg=i2Pik exp(arg) x exp(arg)**x exp(arg*x)
1 (0.00E+00,6.28E+00) (1.00E+00,-2.45E-16) (0.00E+00,-2.21E-01) (1.00E+00,-6.62E-33) (4.00E+00,0.00E+00)
2 (0.00E+00,1.26E+01) (1.00E+00,-4.90E-16) (0.00E+00,-1.10E-01) (1.00E+00,-1.32E-32) (4.00E+00,0.00E+00)
3 (0.00E+00,1.88E+01) (1.00E+00,-7.35E-16) (0.00E+00,-7.35E-02) (1.00E+00,-1.99E-32) (4.00E+00,0.00E+00)
4 (0.00E+00,2.51E+01) (1.00E+00,-9.80E-16) (0.00E+00,-5.52E-02) (1.00E+00,-2.65E-32) (4.00E+00,0.00E+00)
5 (0.00E+00,3.14E+01) (1.00E+00,-1.22E-15) (0.00E+00,-4.41E-02) (1.00E+00,-3.31E-32) (4.00E+00,0.00E+00)
6 (0.00E+00,3.77E+01) (1.00E+00,-1.47E-15) (0.00E+00,-3.68E-02) (1.00E+00,-3.97E-32) (4.00E+00,0.00E+00)
7 (0.00E+00,4.40E+01) (1.00E+00,-1.71E-15) (0.00E+00,-3.15E-02) (1.00E+00,-4.63E-32) (4.00E+00,0.00E+00)
8 (0.00E+00,5.03E+01) (1.00E+00,-1.96E-15) (0.00E+00,-2.76E-02) (1.00E+00,-5.29E-32) (4.00E+00,0.00E+00)
9 (0.00E+00,5.65E+01) (1.00E+00,-2.20E-15) (0.00E+00,-2.45E-02) (1.00E+00,-5.96E-32) (4.00E+00,0.00E+00)
10 (0.00E+00,6.28E+01) (1.00E+00,-2.45E-15) (0.00E+00,-2.21E-02) (1.00E+00,-6.62E-32) (4.00E+00,0.00E+00)
end program xxx
The code is written so that a student can easily change the precision (the wp
parameter) and the right hand side value.
First note that the cmplx()
intrinsic needs to use the kind=
optional argument to ensure consistent precision. Without that argument the result will always be returned with the default real kind. This quirk is due to fortran’s history when it originally supported multiple real floating point kinds but only a single complex kind (using modern fortran terminology here).
The mathematical trick here is that the number 1
can be represented as e^{i 2 \pi k} for any integer k. There are an infinite number of integers, so there are an infinite number of ways to represent the number 1
. For the solution to this problem, only nonzero k values are appropriate. If you let k be a real number, then you can vary k in the neighborhood of an integer value and the result is well defined and smooth, so the inverse function ln() is also well defined and smooth in that neighborhood. The different values of k define the different branches of the ln() function.
if you substitute that expression into the equation to be solved, take the natural log of both sides, and solve for x
you get x=-i ln(4)/(2 \pi k) for k≠0. These are the infinite number of solutions to the problem. The above code loops over some values of the integer k
and evaluates the lhs of the expression in two different ways. One way as (e^{i 2 \pi k})^x and the other way as e^{i 2 \pi k x}. Although these two expressions appear to be the same mathematically, they differ when evaluated in fortran. This is because the first expression loses it branch identity too soon, and there is no way to evaluate the exponential and to stay on the correct branch. This is the z1
value that is printed in the code above, it is always close to 1, or to the complex number (1,0). Only the latter way works to give the desired rhs value of 4. This works because the k
in the numerator correctly cancels the k
in the denominator within x
, keeping both the numerator and the denominator consistent and on the same branch.
So yes, there are solutions to the equation 1^x=4 and you can find them in fortran and verify that they are correct.