# Numeric error in trivial compilation

I came across this strange problem:

dd = Compile[{}, Abs[(10^-200 + 10^-200 I)/(10^-200 + I 10^-200)]^2,
CompilationTarget -> "C"]

dd[]


CompiledFunction::cfne: Numerical error encountered; proceeding with uncompiled evaluation. >>

1


It is same with 10^-100, but these exponents cannot be out of bound with 64bit Visual Studio. What's the problem?

Here is the full code, with the corrected power, which still fails:

f =
Compile[{{σt, _Real, 1}, {sρ, _Complex, 1}, {dt, _Real, 1},
{x, _Real}, {nnt, _Integer}},
Module[{kz = 1/2 x, σ = σt, d = dt, i, kp, c, β, r, k, ρ = sρ * 10^-6, nn = nnt},
kp = kz + 0. I;
k = Sqrt[kz^2 - 4 π (ρ[] - ρ[])];
r = (kp - k)/(kp + k) E^(-2 kp k σ[]^2);
c = {{1 + 0. I, r }, {r , 1 + 0.}};
kp = k;
For[i = 2, i < 2 nn + 2, i++,
If[EvenQ[i],
k = Sqrt[kz^2 - 4 π (ρ[] - ρ[])];
r = (kp - k)/(kp + k) E^(-2 kp k σ[]^2);
β = I kp d[],
k = Sqrt[kz^2 - 4 π (ρ[] - ρ[])];
r = (kp - k)/(kp + k) E^(-2 kp k σ[]^2);
β = I kp d[]
];
c = c.{{E^β, r E^β}, {r E^-β, E^-β}};
(*previous k*)
kp = k
];
Abs[(c[[2, 1]] + 1*^-100 + I 1*^-200 )/(c[[1, 1]] + 1*^-100 + I 1*^-200)]^2]
, CompilationTarget -> "C"]

f[{0, 0, 0}, {0, 24, 5}, {6, 4}, 0.025, 1000]


CompiledFunction::cfne: Numerical error encountered; proceeding with uncompiled evaluation. >>

1.00000000000000


Update: I think what is happening is that the matrix c starts having components too large and at this point it fails. Giving a small imaginary component to \rho helps and the failure comes only at higher values of nnt. I assume this can't be fixed unless I put conditionals to handle the big values of c but this would defeat the purpose of compiling.

 In:= AbsoluteTiming@
f[{0, 0, 0}, {0, 24 + 47 I, 5 + 47 I}, {6, 4}, 0.025, 1000]

Out= {0., 0.296593}


Compile gives two orders of magintude increase in speed!

• Read the documentation. Try, e.g. Compile[{}, Module[{x = 10^10}, Floor[x]]][], (or 10^whatever for your machine arch.) then understand why...
– ciao
Aug 27 '15 at 5:59
• Try using 1.*^-200 instead. Aug 27 '15 at 6:27
• @Guesswhoitis. that appears to compile the result, i.e. CompilePrint[] just has one instruction, namely Return... Aug 27 '15 at 6:56
• @blochwave That's compile-golf Aug 27 '15 at 6:57
• @bloch, as I thought. The numbers involved are after all within the range of $MinMachineNumber and $MaxMachineNumber Aug 27 '15 at 7:16

The reason for the error message is machine integer overflow. The largest machine integer on a 64-bit platform is 2^63 - 1

Developer$MaxMachineInteger == 2^63 - 1 (* True *)  Compare the following examples where the second one overflows but the first one doesn't cf1 = Compile[{}, Developer$MaxMachineInteger];
cf1[]

(* 9223372036854775807 *)

cf2 = Compile[{}, Developer$MaxMachineInteger + 1]; cf2[] (* CompiledFunction::cfn: Numerical error encountered at instruction 2; proceeding with uncompiled evaluation. >> *) (* 9223372036854775808 *)  Turning off integer overflow checking will eliminate the error message but also give an unexpected result cf3 = Compile[{}, Developer$MaxMachineInteger + 1,
RuntimeOptions -> {"CatchMachineIntegerOverflow" -> False}];
cf3[]

(* 0 *)


Inspecting the compiled code for the dd function shows where to expect an integer overflow

Needs["CompiledFunctionTools"]

dd = Compile[{}, Abs[(10^-200 + 10^-200 I)/(10^-200 + I 10^-200)]^2 ];
CompilePrint[dd]

(*  No argument
3 Integer registers
4 Real registers
4 Complex registers
Underflow checking off
Overflow checking off
Integer overflow checking on
RuntimeAttributes -> {}

I1 = 200
I0 = 10
C0 = 0. + 1. I
R2 = 0.
Result = R3

1   I2 = Power[ I0, I1]
2   R0 = I2
3   R1 = Reciprocal[ R0]
4   C1 = R1 + R2 I
5   C1 = C1 * C0
6   C2 = R1 + R2 I
7   C2 = C2 + C1
8   C1 = R1 + R2 I
9   C3 = C0 * C1
10  C1 = R1 + R2 I
11  C1 = C1 + C3
12  C3 = Reciprocal[ C1]
13  C2 = C2 * C3
14  R0 = Abs[ C2]
15  R3 = Square[ R0]
16  Return *)


so the very first instruction is trying to compute 10^200 which exceeds 2^63 - 1.

The following version which uses machine floating point numbers instead of integers does not have an overflow problem, as mentioned in the comments, it directly returns 1.

ddd = Compile[{}, Abs[(10.^-200 + 10.^-200 I)/(10.^-200 + I 10.^-200)]^2 ];
ddd[]

(* 1. *)


See also the documentation for the CompiledFunction::cfn message and for RuntimeOptions`.