440,665 Members | 1,900 Online 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
5 Replies

 P: n/a On Sep 12, 5:00 am, "Jon Slaughter" 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" I wrote a routine to replace Math's Exp method but it turns out to bealmost 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" I wrote a routine to replace Math's Exp method but it turns out to bealmost 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]" "Jon Slaughter" >I wrote a routine to replace Math's Exp method but it turns out to bealmost 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. 