Re: Help on compiling a function
- To: mathgroup at smc.vnet.net
- Subject: [mg115628] Re: Help on compiling a function
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Mon, 17 Jan 2011 05:36:11 -0500 (EST)
----- Original Message -----
> From: "Ramiro" <ramiro.barrantes at gmail.com>
> To: mathgroup at smc.vnet.net
> Sent: Sunday, January 16, 2011 4:54:33 AM
> Subject: [mg115614] Help on compiling a function
> Hi everyone,
>
> I had written about this but wanted to revisit the problem in light of
> the latest Mathematica options for compiling (C and parallelization),
> which I was hoping to try to use on this function, which is the main
> arithmetic and time-consuming part of my simulation.
>
> I have the following code:
>
> example1 =
> Compile[{{n, _Real, 1}, {a, _Real}, {b, _Real}, {t, _Real, 1}},
> Gamma[Total[n] + a]/(Times @@ (Gamma[n + 1])*Gamma[a])*b^a*
> Times @@ (t^n)/(Total[t] + b)^(Total[n] + a)];
>
> but under input such as the following it breaks:
>
> example1[{97.6203, 8.4788, 21.4204, 46.1755}, 1, 1, {39.9342, 7.5820,
> 5.8656, 10.0553}]
> CompiledFunction::cfse: Compiled expression
> 1.33128164105722332870399207`12.920368310128136*^315 should be a
> machine-size real number. >>
> CompiledFunction::cfex: Could not complete external evaluation at
> instruction 4; proceeding with uncompiled evaluation. >>
>
> the problem is just that in the calculation Gamma ends up being really
> big, larger than $MaxMachineNumber so it complains.
>
> Anybody see a way around this? This is the workhorse function of my
> simulation, so any gain in speed helps a lot
>
> Thanks in advance,
> Ramiro
>
> Longer explanation (if interested):
> [...]
You have products/quotients with intermediate overflow because some factors get too big before division brings them back down. The easiest way to avoid this is to use logarithms, sum, then exponentiate. This is workable because your factors are all explicit powers of positives or Gamma function evaluations, so you can use the power-of-log rule and LogGamma.
In[285]:= example1[{97.6203, 8.4788, 21.4204,
46.1755}, 1, 1, {39.9342, 7.5820, 5.8656, 10.0553}]
During evaluation of In[285]:= CompiledFunction::cfex: Could not complete external evaluation at instruction 3; proceeding with uncompiled evaluation. >>
Out[285]= 1.062645560684518*10^-11
In[286]:= example2 =
Compile[{{n, _Real, 1}, {a, _Real}, {b, _Real}, {t, _Real, 1}},
Module[{log1, log2, log3, log4, log5},
log1 = LogGamma[Total[n] + a];
log2 = Plus @@ LogGamma[n + 1];
log3 = LogGamma[a];
log4 = Plus @@ (n*Log[t]);
log5 = (Total[n] + a)*Log[Total[t] + b];
Exp[log1 - (log2 + log3) + log4 - log5]*b^a]]
In[287]:= example2[{97.6203, 8.4788, 21.4204,
46.1755}, 1, 1, {39.9342, 7.5820, 5.8656, 10.0553}]
Out[287]= 1.062645560684574*10^-11
Daniel Lichtblau
Wolfram Research
- Follow-Ups:
- CUDADot bug?
- From: Ulrich Arndt <ulrich.arndt@data2knowledge.de>
- CUDADot bug?