How to code the expression with least roundoff error

[bsmile]

The expression to code is: (z complex, m some positive integer)

z^m ( ln(z/(z-1)) - sum( 1/(k z^k), {k=1,m} ) )

The way I code is (in fortran) (in case z<2)

Z_1=1.0_q/Z

ZV1=LOG(Z/(Z-1))
ZK=1.0_q
DO K=1,M; ZK=ZK*Z_1
ZV1=ZV1-ZK/K
ENDDO
ZV1=ZV1/ZK

But it gives quite big error by this bruteforce way of coding. Can anybody help me out to suggest a better way to code it?

Thanks,

 

[DrClaude]

I would do the sum first, then subtract it from the logarithm. Another possibility is to calculate
$$
z^m \ln \frac{z}{z-1} - \sum_{k=1}^m \frac{z^{m-k}}{k}
$$
Depending on the value of z and m, one approach might be better than the other.

 

[rcgldr]

Summing smaller numbers first will reduce truncation error, so this may be better. I assume the common factor z^m can be placed outside the sum, even if z is complex?

$z^m ( \ln \frac{z}{z-1} - \sum_{k=m}^1 \frac{z^{-k}}{k})$

 

[bsmile]

Thanks for your help. I wonder whether the issue comes from Log[z/(z-1)]. To test it, I evaluated two equivalent expressions with Mathematica making use of the following identity

ln(z/(z-1))=∫^{1}_{0} 1/(z-x) dx

Ideally, the difference between left and right hand sides should be zero, but Mathematica tells me otherwise, as shown from the attached figure.

Can you explain to me why as usually we would assume the elementary functions have been extremely well carried out as built-in functions and should attain the highest machine precision?

(how to upload a figure?)

 

[DrClaude]

(how to upload a figure?)

Below the text box, to the right, you should find a botton marked "UPLOAD".

 

[bsmile]

Below the text box, to the right, you should find a botton marked "UPLOAD".

Thanks, I had been clicking on the figure icon on top of the textbox, and it really didn't do what I wanted. Here is the figure showing the two otherwise equivalent expressions have error bigger than machine precision and has a strong dependence on the input parameters.

I understand the precision of the output from the built-in functions will have a dependence on the input parameters, but they should do at least what we ordinary people can achieve. Obviously if |z| is small, we should be able to make it more accurate by using taylor expansion than what gfortran

returns us.

 

[bsmile]

One thing I might mention is that the lack of accuracy can be in NIntegrate as well as it might fail to reach prescribed precision goal I set. I have written to Mathematica to see whether there is way to detect the condition satisfied or not.

 

[pbuk]

One thing I might mention is that the lack of accuracy can be in NIntegrate as well as it might fail to reach prescribed precision goal I set. I have written to Mathematica to see whether there is way to detect the condition satisfied or not.

Have you set WorkingPrecision as well as PrecisionGoal?

 

[pbuk]

The expression to code is: (z complex, m some positive integer)

z^m ( ln(z/(z-1)) - sum( 1/(k z^k), {k=1,m} ) )

The way I code is (in fortran) (in case z<2)

Z_1=1.0_q/Z

ZV1=LOG(Z/(Z-1))
ZK=1.0_q
DO K=1,M; ZK=ZK*Z_1
ZV1=ZV1-ZK/K
ENDDO
ZV1=ZV1/ZK

But it gives quite big error by this bruteforce way of coding. Can anybody help me out to suggest a better way to code it?

Thanks,

I have forgotten all the FORTRAN I once knew (good riddance!) but I do remember some hard-to-track-down bugs due to not forcing double precision - particularly with complex variables.