# How to code the expression with least roundoff error

Tags:
1. Jan 29, 2016

### 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,

2. Jan 29, 2016

### Staff: Mentor

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.

3. Jan 29, 2016

### 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})$

4. Jan 29, 2016

### 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?

5. Jan 30, 2016

### Staff: Mentor

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

6. Jan 30, 2016

### bsmile

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.

7. Feb 1, 2016

### 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.

8. Feb 1, 2016

### MrAnchovy

Have you set WorkingPrecision as well as PrecisionGoal?

9. Feb 1, 2016

### MrAnchovy

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.