Loss of precision in exponentiation - questions

[Goicochea__Javier]

Dear Fortran Masters, 

I have two simple questions.

1 - In the operation x**b, where b is always an integer. Is there any loss of precision if b is declared as real(8) or integer?

2 - In the operation x**b, where b is always real positive number. Is there any loss of precision if b is declared as real(8) or real(4)?

Why or why not? Kindly advise. Thank you, Javier

program powx

	implicit none

	real(8)	:: x
	real(8)	:: a
	real(4)	:: b
	integer	:: i
	
	x = 1.23456789d0
	i = 4
	a = 3d0
	b = 1.23456789
	
	! Question 1
	write (*,*) x**(i-1), x**dble(i-1), x**a, x**3.0
	
	! Question 2
	write (*,*) x**1.23456789, x**dble(1.23456789), x**x, x**b
	
end program powx

[Arjen_Markus]

I do not want to be pedantic but the correct term is accuracy, not precision. Precision relates to the number of bits in the value as stored in the computer, accuracy relates to the possible deviation of the answer from the mathematically exact one. Assuming you mean accuracy, here are my answers (with emphasis on "my").

ad 1.

The compiler will very probably use different methods for the calculation. I am no expert on how this is actually arranged, but one way for the calculation of x**b (b real) is: exp( b * log(x) ), whereas with x**n (n integer) you can do repeated multiplications. For small values of the exponent that may be faster and more accurate than the exp/log method. But it all depends. And which one is more accurate

ad 2.

If I understand it correctly, then it makes no difference if you do x**b (x double precision, b single precision) or x**double(b), since b has to be promoted to double precision anyway. What does make a difference is: Define b as double precision and set its value to b = 1.23456789d0. If you leave out the "d0" (or anything equivalent, like "_double", with "double" the kind for double precision) , then the value on the right-hand side is single precision, no matter how many decimals you specify, and will be converted into double precision. That may be slightly different from the double-precision value you get with "d0".

Another question is: does this influence the outcome of your program in a significant way? If so, you may need to examine the algorithms used more closely as they are very sensitive to the numerical accuracy.

Regards,

Arjen

[mecej4]

The following program may be instructive in understanding the nice points made by Arjen Markus. The value 1.23456789 cannot be represented exactly in a 4-byte IEEE real.

program tst
real :: b_sp=1.23456789
double precision :: b_dp=1.23456789d0
integer div(2),civ(2)
!
div=transfer(b_dp,0,2)         ! double precision value -> hex
civ=transfer(dble(b_sp),0,2)   ! single precision -> double -> hex
!
write(*,'(F15.10)')b_sp        ! print out more digits to reveal roundoff
write(*,'(2Z8,/,2Z8)')div(2),div(1),civ(2),civ(1)
end

The first line of output from the program reveals that when the value 1.23456789 is stored into a 32-bit real and converted back to decimal the output value is not exactly the same as the input value. The subsequent lines illustrate the points that Arjen made.

   1.2345678806
3FF3C0CA4283DE1B
3FF3C0CA40000000

[TimP]

A brief summary of accuracy expected from Intel math functions is in

https://software.intel.com/en-us/articles/differences-in-floating-point-arithmetic-between-intel-xeon-processors-and-the-intel-xeon

You will note the higher expectations for scalar math functions than for the vector ones.

Scalar exponentiation typically takes precautions such that integral real exponents should produce the same results as integer exponents (at possibly value-dependent cost in time).  These precautions aren't conducive to vectorization, which might involve the log and exponent scheme mentioned previously.

If you like to experiment, you might compare x**(i+.5) vs. x**i*sqrt(x) and the like in scalar and vector modes.

You may be interested also in checking whether the compiler treats a constant integral real exponent as an integer (does it throw a compile or run-time error on negative base, does it call an integer function or expand with multiplication in line?).  Compilers were pushed in this direction by popular benchmarks which inexplicably used real exponents (maybe written by C programmers).

[LRaim]

But the Intel White Paper does not include any information about the method used to evaluate basic math functions (exp,log,log10,sin,cos, etc) and the precision expected. I remember the documentation of an IBM  Fortran compiler  describing  for each math function which series expansion was used and the precision.

Is such a document available from Intel ?

Regards

[Steven_L_Intel1]

The method used will vary depending on various compiler options and also on what the processor type is. In some cases, CPU instructions are used for these basic operations, but there can be variations using the Short Vector Math Library, etc. The setting of options such as /fp can also affect this.

[Goicochea__Javier]

Thank you for your replies.

Arjan, very often we have legacy code that uses x**b where x is real(8) and b real(4) or integer. So I was not sure about the impact of this. It is also difficult to anticipate or identify if the change in the results of legacy code might be produced due to the declaration used for b or for instance the change in the compiler options (i.e. /fp). In these codes there are many places where real(4) and real(8) are mixed. 

Mecej and Prince, I will keep in mind your suggestions. 

Luigi, yes documentation will be beneficial as always! 

Thank you, Javier