1
$\begingroup$

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 π (ρ[[2]] - ρ[[1]])];
    r = (kp - k)/(kp + k) E^(-2 kp k σ[[1]]^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 π (ρ[[3]] - ρ[[1]])];
        r = (kp - k)/(kp + k) E^(-2 kp k σ[[2]]^2); 
        β = I kp d[[1]],
        k = Sqrt[kz^2 - 4 π (ρ[[2]] - ρ[[1]])];
        r = (kp - k)/(kp + k) E^(-2 kp k σ[[3]]^2); 
        β = I kp d[[2]] 
      ];
      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[67]:= AbsoluteTiming@
 f[{0, 0, 0}, {0, 24 + 47 I, 5 + 47 I}, {6, 4}, 0.025, 1000]

Out[67]= {0., 0.296593}

Compile gives two orders of magintude increase in speed!

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

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.