By using this site, you agree to our updated Privacy Policy and our Terms of Use. Manage your Cookies Settings.
440,665 Members | 1,900 Online
Bytes IT Community
+ Ask a Question
Need help? Post your question and get tips & solutions from a community of 440,665 IT Pros & Developers. It's quick & easy.

Math libarary

P: n/a
I wrote a routine to replace Math's Exp method but it turns out to be almost
2x slower ;/ (well, actually its about 1.5x in release)

I'm essentially using a lookup table and interpolate between integer
values(but even just looking up its still as slow).

public static double Exp(double x)
{
int n = (int)Math.Floor(x);
//if ((n 2351) || (n < -2351)) return Math.Sign(x) *
double.PositiveInfinity;

double f = x - n;
double ff = 1;
double e = 1;

if (f != 0)
for (int i = 1; i < 10; i++)
{
ff = ff * f;
e += invnfac[i] * ff;
}

return ExpIp[n + 2351]*e;
}

What I did was precompute n! and ExpIp is exp(k) for integer k. exp(+-2351)
is the maximum that a double can hold so no reason to go beyond that.

If I comment out the loop I get almost there speed but still about 1%
slower.

How the hell are they computing so fast? Normally these Math library
routines are very slow but either my code/method sucks or Math .NET is
fairly optimized?

Anyone know whats going on here?

Thanks,
Jon
Sep 12 '07 #1
Share this Question
Share on Google+
5 Replies


P: n/a
On Sep 12, 5:00 am, "Jon Slaughter" <Jon_Slaugh...@Hotmail.comwrote:

<snip>
How the hell are they computing so fast? Normally these Math library
routines are very slow but either my code/method sucks or Math .NET is
fairly optimized?
I'd expect the Math library to be very heavily optimised for things
like Exp - it's a natural place to put such optimisation, after all.
Math.Exp is actually an extern method, not managed code.

Jon

Sep 12 '07 #2

P: n/a
"Jon Slaughter" <Jo***********@Hotmail.comwrote in message
news:pv*****************@newssvr25.news.prodigy.ne t...
>I wrote a routine to replace Math's Exp method but it turns out to be
almost 2x slower ;/ (well, actually its about 1.5x in release)

I'm essentially using a lookup table and interpolate between integer
values(but even just looking up its still as slow).

public static double Exp(double x)
{
int n = (int)Math.Floor(x);
//if ((n 2351) || (n < -2351)) return Math.Sign(x) *
double.PositiveInfinity;

double f = x - n;
double ff = 1;
double e = 1;

if (f != 0)
for (int i = 1; i < 10; i++)
{
ff = ff * f;
e += invnfac[i] * ff;
}

return ExpIp[n + 2351]*e;
}

What I did was precompute n! and ExpIp is exp(k) for integer k.
exp(+-2351) is the maximum that a double can hold so no reason to go
beyond that.

If I comment out the loop I get almost there speed but still about 1%
slower.

How the hell are they computing so fast? Normally these Math library
routines are very slow but either my code/method sucks or Math .NET is
fairly optimized?

Anyone know whats going on here?

Thanks,
Jon

Adding to Jon's ( Skeet) reply, Math.Exp is finally calling the CRT library
function Exp, this is unmanaged highly optimized code that uses the FPU for
this.

Willy.

Sep 12 '07 #3

P: n/a

"Jon Slaughter" <Jo***********@Hotmail.comwrote in message
news:pv*****************@newssvr25.news.prodigy.ne t...
>I wrote a routine to replace Math's Exp method but it turns out to be
almost 2x slower ;/ (well, actually its about 1.5x in release)
Jon, when I suggested a lookup table for your particular problem, I didn't
mean for you to make a general-purpose Exp routine implemented with a lookup
table. I meant to have a lookup table mapping the 256 possible values that
the integer variable appearing in your expression could take, to the value
of the whole expression.
>
I'm essentially using a lookup table and interpolate between integer
values(but even just looking up its still as slow).

public static double Exp(double x)
{
int n = (int)Math.Floor(x);
//if ((n 2351) || (n < -2351)) return Math.Sign(x) *
double.PositiveInfinity;

double f = x - n;
double ff = 1;
double e = 1;

if (f != 0)
for (int i = 1; i < 10; i++)
{
ff = ff * f;
e += invnfac[i] * ff;
}

return ExpIp[n + 2351]*e;
}

What I did was precompute n! and ExpIp is exp(k) for integer k.
exp(+-2351) is the maximum that a double can hold so no reason to go
beyond that.

If I comment out the loop I get almost there speed but still about 1%
slower.

How the hell are they computing so fast? Normally these Math library
routines are very slow but either my code/method sucks or Math .NET is
fairly optimized?

Anyone know whats going on here?

Thanks,
Jon

Sep 12 '07 #4

P: n/a
Thanks guys.
..NET just gets better and better ;)

Sep 12 '07 #5

P: n/a

"Ben Voigt [C++ MVP]" <rb*@nospam.nospamwrote in message
news:%2****************@TK2MSFTNGP04.phx.gbl...
>
"Jon Slaughter" <Jo***********@Hotmail.comwrote in message
news:pv*****************@newssvr25.news.prodigy.ne t...
>>I wrote a routine to replace Math's Exp method but it turns out to be
almost 2x slower ;/ (well, actually its about 1.5x in release)

Jon, when I suggested a lookup table for your particular problem, I didn't
mean for you to make a general-purpose Exp routine implemented with a
lookup table. I meant to have a lookup table mapping the 256 possible
values that the integer variable appearing in your expression could take,
to the value of the whole expression.
Hehe, no. I was writing it because I'm kinda putting together a small vector
library so I can do some problems that I've been wanting and are to slow in
matlab. I figured that I might as well put together a simple numerical
library that uses lookup tables while I'm at it. So I threw together what I
thought would have been a bette routine than Math.Exp and I was wrong ;/
Actually all I'll probably never use Exp inside a critical loop as its just
for initialization but who knows. Its nice to see that Microsoft finally
did there homework with .NET.

My real problem is that for a 2D physics problem I have, it requires 4
nested loops just to caculate one period of time. (so technically 5 if I
want to evolve the system)

The issue is that to do the simulation I have to take into account boundary
values. For example, Suppose I want to take the gradient of a scalar field
S. This requires I compute (S[i+1, j] - S[i,j])/dx or equivilently (S[i,
j] - S[i-1,j])/dx or (S[i+1, j] - S[i-1,j])/dx/2.

So if I'm at the boundary then I have a problem. Doing checks and applying
the correct formula would be inefficient.

What I do is when I want to specifically compute the something like this
that as to deal with the boundary is break it up into 9 seconds and deal
with each one seperately. But how to codify this in general where there is
ultimately only 1 main loop(which might contian sub loops but never
duplicating loops).

What I mean suppose I want to compute the gradient then the laplacian.

Both of these depend on S only but have to take into account the boundaries.

Right now I have two functions that essentially loop over S seperately.

So if I call

Grad(S);
Laplacian(S);

I have 2 loops when in reality I could have "consolidated" them into one.
Now this isn't a big deal until you consider that I have to break such a
loop into 9 parts to handle the boundaries. (I can wrap the boundaries but I
do not want to do that for this problem... I will implement that later as an
option... which will be much easier).

The only real option is to codify whatever I'm doing by hand but it gets
messy real quick. I do have routines that can compute the gradient and
laplacian at only a point instead of computing it on the entire field but it
has to check if the point is on the boundary. What this means is that if the
user called those function it would be slow because of the checks(abit
easier to use though which is why I included it).

The following is the code I use to compute the laplacian of a scalar field.
(I do have another function that actually returns a new scalar field but I
want to try to avoid allocating new fields because its slow and could get
very mess quick if one starts doing some algebra on the field)
The general scheme of the laplacian is applied to every manipulation of a
scalar field. Basically just handle corners, sides, and then inner part. In
this case I'm calculating the second derivatives but in reality it could be
much more complicated.

The real meat of the function is this which is simple... the problem is the
boundaries as they require special techniques.
// Inner Matrix

for (int j = 1; j < SF.Ny - 1; j++)

for (int i = 1; i < SF.Nx - 1; i++)

{

x = (SF[i + 1, j] + SF[i - 1, j] - 2 * SF[i, j]) / dx2;

y = (SF[i, j + 1] + SF[i, j - 1] - 2 * SF[i, j]) / dy2;

SF2[i, j] = x + y;

}

I was thinking about making a routine where the user could register a
callback in each part so they could "hook" into the boundaries and put
whatever computation they wanted so they wouldn't have to actually write out
all the loops... but all those function calls would be slow unnecessarily I
suppose that since each side is mirrored inwards I could optmize it so that
there is only "one side" in some sense but I think that might be impossible
to actually do because it would depend on figuring how the indicies used in
the callback or end up with a lot of checks.

Anyways, chances are I'll have to recode this loop every time I want a tight
loop but I feel for my problem that if I can abstract things efficiently
that its going to get messy real quick.

Maybe someone has some ideas though?

Thanks,

Jon



public unsafe static ScalarField2d Laplacian(ScalarField2d SF, ref
ScalarField2d SF2)

{

double x, y, dx2 = SF.dx * SF.dx, dy2 = SF.dy * SF.dy;

SF2.dx = SF.dx; SF2.dy = SF.dy;

// Evaluate the Laplacian at corners using single sided differences

x = (-SF[3, 0] + 4 * SF[2, 0] - 5 * SF[1, 0] + 2 * SF[0, 0]) / dx2;

y = (-SF[0, 3] + 4 * SF[2, 0] - 5 * SF[1, 0] + 2 * SF[0, 0]) / dy2;

SF2[0, 0] = x + y;

x = (2 * SF[SF.Nx - 1, 0] - 5 * SF[SF.Nx - 1 - 1, 0] + 4 * SF[SF.Nx - 1 - 2,
0] - SF[SF.Nx - 1 - 3, 0]) / dx2;

y = (-SF[SF.Nx - 1, 3] + 4 * SF[SF.Nx - 1, 2] - 5 * SF[SF.Nx - 1, 1] + 2 *
SF[SF.Nx - 1, 0]) / dy2;

SF2[SF.Nx - 1, 0] = x + y;

x = (2 * SF[SF.Nx - 1, SF.Ny - 1] - 5 * SF[SF.Nx - 1 - 1, SF.Ny - 1] + 4 *
SF[SF.Nx - 1 - 2, SF.Ny - 1] - SF[SF.Nx - 1 - 3, SF.Ny - 1]) / dx2;

y = (2 * SF[SF.Nx - 1, SF.Ny - 1] - 5 * SF[SF.Nx - 1, SF.Ny - 1 - 1] + 4 *
SF[SF.Nx - 1, SF.Ny - 1 - 2] - SF[SF.Nx - 1, SF.Ny - 1 - 3]) / dy2;

SF2[SF.Nx - 1, SF.Ny - 1] = x + y;

x = (2 * SF[0, SF.Ny - 1] - 5 * SF[1, SF.Ny - 1] + 4 * SF[2, SF.Ny - 1] -
SF[3, SF.Ny - 1]) / dx2;

y = (2 * SF[0, SF.Ny - 1] - 5 * SF[0, SF.Ny - 1 - 1] + 4 * SF[0, SF.Ny - 1 -
2] - SF[0, SF.Ny - 1 - 3]) / dy2;

SF2[0, SF.Ny - 1] = x + y;

// left and right column

for (int j = 1; j < SF.Ny - 1; j++)

{

x = (-SF[3, j] + 4 * SF[2, j] - 5 * SF[1, j] + 2 * SF[0, j]) / dx2;

y = (SF[0, j + 1] + SF[0, j - 1] - 2 * SF[0, j]) / dy2;

SF2[0, j] = x + y;

x = (2 * SF[SF.Nx - 1, j] - 5 * SF[SF.Nx - 1 - 1, j] + 4 * SF[SF.Nx - 1 - 2,
j] - SF[SF.Nx - 1 - 3, j]) / dx2;

y = (SF[SF.Nx - 1, j + 1] + SF[SF.Nx - 1, j - 1] - 2 * SF[SF.Nx - 1, j]) /
dy2;

SF2[SF.Nx - 1, j] = x + y;

}

// top and bottom row

for (int i = 1; i < SF.Nx - 1; i++)

{

x = (SF[i + 1, 0] + SF[i - 1, 0] - 2 * SF[i, 0]) / dx2;

y = (-SF[i, 3] + 4 * SF[i, 2] - 5 * SF[i, 1] + 2 * SF[i, 0]) / dy2;

SF2[i, 0] = x + y;

x = (SF[i + 1, SF.Ny - 1] + SF[i - 1, SF.Ny - 1] - 2 * SF[i, SF.Ny - 1]) /
dx2;

y = (-SF[i, SF.Ny - 1 - 3] + 4 * SF[i, SF.Ny - 1 - 2] - 5 * SF[i, SF.Ny -
1 - 1] + 2 * SF[i, SF.Ny - 1]) / dy2;

SF2[i, SF.Ny - 1] = x + y;

}

// Inner Matrix

for (int j = 1; j < SF.Ny - 1; j++)

for (int i = 1; i < SF.Nx - 1; i++)

{

x = (SF[i + 1, j] + SF[i - 1, j] - 2 * SF[i, j]) / dx2;

y = (SF[i, j + 1] + SF[i, j - 1] - 2 * SF[i, j]) / dy2;

SF2[i, j] = x + y;

}

return SF2;

} // Laplacian


Sep 12 '07 #6

This discussion thread is closed

Replies have been disabled for this discussion.