Mathematica CForm / FortranForm bug?

I am trying to output a very complicated matrix (~1.3MB in plain text) from Mathematica to use in a Fortran program. When I do this (via Splice) the resulting matrix is off by ~2% when numerical values are given to the variables. This is a problem since there needs to be an eigenvalue that is exactly zero, and the composition of the eigenvectors needs to be exactly correct.

I've done all the usual due diligence regarding precision, correct variables, proper diagonalization code, etc. and It came down to either Fortran itself being unable to cope with such a large matrix or Mathematica messing up the FortranForm output.

So I made Mathematica give me the CForm of the matrix and tried that. It was also ~2% off from what it should be, more astoundingly, it was the same (within machine precision) as the FortranForm matrix!

Has anyone come into contact with this kind of problem? Do you have any idea what might cause it? I dread the need to go through 25000 lines of Mathematica formatted Fortran code to figure this one out.

EDIT: The matrix in question is complicated, not large. It is only 6x6 but each element is individually algebraically very messy including trigonometric functions, logarithms, and various roots and powers.

The Plaintext of the (1,1) element of our matrix, the C code, and the Fortran code. Sane parameter开发者_StackOverflow社区 values are: 0 < lambda, kappa, Y*** < 1; all others between 100 and 1000.

Looking at the expressions for the first element, this is almost certainty a problem because of finite precision floating point arithmetic. It seems to me there are two ways to try and approach this problem. The first is to use more precision in the C code. Perhaps you could use the macro pre-processor (or just find and replace) to change all of the declarations in the C code from double to long double. Depending on what compiler you are using that will get you from 64 bits to 80 or 128 (unless you use Microsoft's compiler which treats long doubles as standard doubles). If you post what C compiler you are using I could help look at what changes you would need to make to go that route. At the very end you can cast down to a double if you need to pass it to Mathematica or another C program that is just using doubles.

Also, C compilers can have different floating point modes. Might be worth checking whether the mode is set for precise, strict or fast. If it is on fast or fast-math then that might be contributing to your rounding problems.

The other route is to use Mathematica to rearrange the terms of each element to optimize for floating point accuracy. While in general (a*b)*c=(c*b)*a in the land of floating point arithmetic this isn't true. The best heuristic I have heard of is:

1. Addition: add the terms from smallest to largest to minimize rounding errors
2. Subtraction: avoid or delay subtracting values that are very close to each other in value so you maintain as many significant digits as possible
3. Multiplication/Division: avoid/delay multiplying by very small numbers or dividing by very large numbers for reasons of significant digits again.

In practical terms this means factoring these expressions is your best bet, especially problematic terms. if a1~=a2 then the expression (a1-a2)(b-d) should have more accuracy/significant digits than (a1-a2)*b-(a1-a2)*d because you have one less multiply operation.

Could it be the following simple thing: Mathematica is quite clever about numerical evaluation, so for example these two expressions

N[(10^18+100)-10^18]
N[10^18+100]-N[10^18]

resul in quite different output

Out[3]= 100.
Out[4]= 128. 

With C or Fortran you don't get any babysitting and the expression would be evaluated as in the second line.
Maybe you could look at a single matrix entry and try to evaluate it in Mathematica as C/Fortran would, and see if this matches what you see?

EDIT:
It seems that a simple test for problems with intermediate precision is to use Compile which drops all precision checks, so that

In[10]:= Compile[{a,b},a+100-b][10^18,10^18]
Out[10]= 128.

Could you check whether a compiled version of your matrix agrees with C/Fortran?