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$
8
  • 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, 2015 at 5:59
  • $\begingroup$ Try using 1.*^-200 instead. $\endgroup$ Aug 27, 2015 at 6:27
  • 1
    $\begingroup$ @Guesswhoitis. that appears to compile the result, i.e. CompilePrint[] just has one instruction, namely Return... $\endgroup$ Aug 27, 2015 at 6:56
  • 2
    $\begingroup$ @blochwave That's compile-golf $\endgroup$ Aug 27, 2015 at 6:57
  • $\begingroup$ @bloch, as I thought. The numbers involved are after all within the range of $MinMachineNumber and $MaxMachineNumber $\endgroup$ Aug 27, 2015 at 7:16

1 Answer 1

5
$\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 and acknowledge you have read our privacy policy.

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