Chris Jones wrote: gcc -ffast-math ? -- Rex

Chris Jones wrote: [snip] If the math package is written at all well, then you aren't going to speed it up by lowering the precision, as you are already getting hardware boost. The Intel architecture (I presume that's what you are using) includes hardware computation of logs etc. For example, there is a single instruction FPATAN which does an arctangent. Writing some code to do a truncated series with lower precision is going to be SLOWER. Mike -- p="p=%c%s%c;main(){printf(p,34,p,34);}";main(){printf(p,34,p,34);} Oppose globalization and One World Governments like the UN. This message made from 100% recycled bits. You have found the bank of Larn. I can explain it for you, but I can't understand it for you. I speak only for myself, and I am unanimous in that!

Hi, Thanks for your feedback. Yes, I have profiled the code, quite extensively, using the valgrind/calltree application. From this I know this that I'm know I've tidying up this to the point where its hard to find big improvements, the cpu time is fairly well spread around, not isolated in a few places. So am now looking a a few places where math calls are taking more time than I would hope. I'm not going to get factors in speed in the overall application, but I hope in a few places things can be improvemed a lot locally. Also, the project is not small, massive in fact, and I'm only writting one small part. If you are interested you can find it here http://lhcb-release-area.web.cern.ch/LHCb-release-area/DOC/brunel/release... It also has to be supported on a *lot* of hardware. Basically gcc 3.2.3 based Scientific Linux 3 machines, gcc 3.4.6 SL4 machines (32 and 64 bit) and (not my decision), windows VC 7.1. I cannot rely on for instance SSE math calls etc. Taking an example from another thread, one place I'm trying to understand is where I use atan2 see http://www.hep.phy.cam.ac.uk/~jonesc/atan2.png for the profiler output. atan2 is taking 50% of the time of this method. Not here I don't need that much precision on the result - say +- O(2*pi/100). Anything you can suggest here - The code is here http://www.hep.phy.cam.ac.uk/~jonesc/RichPhotonRecoUsingCKEstiFromRadius.cpp ( note though its full of internal classes etc...) cheers Chris

Out of interest, whats needs to be done to use these, some gcc flag ? My profiling suggests that atan2(x,y) is taking ~ 50% of some method speed. Yes, the methods is already fast, but still is this reasonable if hardware instructions are being used. Is there some way I can check to see exactly what is used ? Yes, I considered coding up my own methods along this line, but figured getting it correct was not a small amount of effort, and was hoping to find some project that had done it for me. cheers Chris

On Sat, 2007-09-01 at 09:41 -0500, Michael Hennebry wrote: That would be PJ Plauger's The Standard C Library, Prentice Hall, 1992 0-13-131509-9. Most of his math lib implementation is based on Cody and Waite, Software Manual for the Elementary Functions, Prentice Hall, 1980 (sorry, he doesn't give the ISBN). -- Matthew Saltzman Clemson University Math Sciences mjs AT clemson DOT edu http://www.math.clemson.edu/~mjs

On Sat, 2007-09-01 at 20:27 -0700, Les wrote: A wonderful book in it's day, but far out of date now. I'm curious (but we should take it off-list) what you do that you use G&N for. -- Matthew Saltzman Clemson University Math Sciences mjs AT clemson DOT edu http://www.math.clemson.edu/~mjs

On Tue, 2007-09-04 at 18:17 -0500, Mike McCarty wrote: Early NR, especially NR in C, had some real laughers. (BTW, does Fortran 77 have an equivalent to fmod()?. If not, that might explain where they picked up the repeated subtraction idea. They really didn't know much C back in the early days...) I haven't seen any reviews yet, but the third edition is just now coming out, after 15 years. More information at http://www.nr.com. (I'm not endorsing anything here. I have no affiliation with anyone involved with NR.) Still may depend on the processor and exactly what instructions are executed doing the approximation. Just because it's hardware, that doesn't mean it will necessarily take only a few clock cycles. -- Matthew Saltzman Clemson University Math Sciences mjs AT clemson DOT edu http://www.math.clemson.edu/~mjs

On Tue, 2007-09-04 at 18:17 -0500, Mike McCarty wrote: I do have to agree with your assessment of their algorithms. But having a working algorithm means I only have to find the optimizations. And sometimes what seems archaic may be able to take advantage of compiler and processor optimizations to achieve faster results. The business of subtraction is one cycle per subtraction, where as fmod is multiple cycles to begin with plus call and return overhead. If the iterative is known to be some small number repeated subtraction may be faster. Only some practical work with the algorithm will tell you the real results. The same is true of Floating point operations vs integers. When floats had to be calculated by loops with an integer processor, they were expensive and integers were faster. Now with high speed floating point units, simple float operations are quite fast if done in line. Ditto for doubles. It costs no more to calculate doubles than singles except when you store and retrieve them on a 32 bit machine (if your blocking is set right.) On a 64 bit machine, doubles may actually be faster since you don't have to truncate or do the store offsets (note that this depends on the hardware implementations inside the processor and the microcode used to achieve the doubles and storage calculations). With some operations, the operation of the ram may be important, impacting due to cyclic overhead. Processors with high i/o bandwidth work well as single cpu, but suffer a hit when in dual cpu due to inability to overlap cycles as effectively. They may well be faster, but it depends a lot on the algorithm, and relative timing of the calculation and results, some of which can be controlled by the programmer directly, and some of which may be limited in the processor design. Ditto for memory access. DDRR ram can do i/o overlap if the processor and mb electronics can handle it. So dual processors which have certain addressing setups can both be full speed and overlapped if the other i/o functions can support it. When discussing algorithm timing, only the algorithm being used and its variants can be discussed for realtime applications. This is one reason that benchmarks are basically useless in choosing a processor for real time applications, unless you are using the bench mark algorithm in your specific application. I have spent many many hours optimizing code. I have gradually come to the conclusion that I can optimize my own code by 10 to 20 percent each cycle I put into it, along with some additional benefit from hardware advances in each cycle. This is especially true now that the hardware cycle is into the 13-14 month timeframe. Sometimes just taking a different perspective on the problem will help. I have chased various trig functions for times then figured out a relative way to achieve the same effective results with less overhead or sometimes in hardware that cost less than $40.00. Using a phase comparator to find angular offsets is one example. Sometimes a good step back, and a look at the ultimate goal will help as well. You may find that the results are more about the value of the peak than the angular offset for example. Or an FIR filter may be faster than an fft in some cases to find a specific freq value. There are more ways to tackle problems than most of us can imagine. That is why we call upon the local wizards to help us. Regards, Les H

Bruno Wolff III wrote: Some models of CRAYS had no fixed point instructions at all. Today, the 486 and later machines in the Intel architecture are FASTER at floating point than with integers. The best way to get the speed this fellow is looking for is to ensure that the fp hardware is actually being used. Mike -- p="p=%c%s%c;main(){printf(p,34,p,34);}";main(){printf(p,34,p,34);} Oppose globalization and One World Governments like the UN. This message made from 100% recycled bits. You have found the bank of Larn. I can explain it for you, but I can't understand it for you. I speak only for myself, and I am unanimous in that!

While we are at this, can someone please point me to some C code for integrating a function numerically? Something that is good? Many thanks and best wishes, Trotter ____________________________________________________________________________________ Building a website is a piece of cake. Yahoo! Small Business gives you all the tools to get online. http://smallbusiness.yahoo.com/webhosting

On Fri, 2007-09-07 at 08:40 -0700, Globe Trotter wrote: One thing I have found is that specifying the range over which I want integration to occur is important. In some cases the function will generate a near infinite list of terms, and your system may be tied up for a long time. (doesn't complete). Regards, Les H

On Sat, 1 Sep 2007 11:28:28 +0100 Chris Jones <jonesc(a)hep.phy.cam.ac.uk&gt; wrote: You can use the binutils tools to disassemble the code if you want, or gdb to trace the instructions and step them (stepi) one at a time. The Intel hardware FP is pretty fast especially on the later processors and there's probably not a lot you can do to beat it. The AMD processors have some limited low precision trig functions in the 3DNow extensions designed for fast 3D gaming (before the cards started doing all the work ;)) Failing that small lookup tables are about the only option although if your precision is high the cache miss cost rapidly exceeds the gain from avoiding the FP math. Alan

On Saturday 1 September 2007 9:36:02 am Mogens Kjaer wrote: Talyor expansions are valid when your argument is 'small' i.e. sin(x) ~ x tan(x) ~ x cos(x) ~ (1-(x)*(x)/2) etc. only work when x is small, and the error increases as x does (since the size of the truncated terms become more important). Yes you can include more terms but that only allows you to go to large x before the errors explode. I'm well aware of these series, and already use them when appropriate. In this case, I need tan(x) with fixed errors, between 0 and 2pi, the whole range. for instance, a really course estimate for atan2(x,y) can be made, by just comparing the signs of x and y, i.e. if x>0 and y>0, atan2(x,y) is between 0 and pi/2 ... Chris

On Saturday 1 September 2007 12:53:36 pm Mogens Kjaer wrote: yes, anthing is possible. However, before I go to the effort on implementing something myself, I was hoping to find some OSS solution. It struck me as such a common issue there must be some solution out there already... Chris

Mogens Kjaer wrote: ... You could try this piece of code: static double atan2x(double y, double x) { double z; if(x==0) { if(y<0) return -M_PI_2; else if(y>0) return M_PI_2; else return 0; } else { z = y/x; if(z>=-1 && z<=1) { z=z*(15.0+4.0*z*z)/(15.0+9.0*z*z); if(x>0) return z; else { if(y>=0) return M_PI+z; else return z-M_PI; } } else { z=x/y; z=M_PI_2-z*(15.0+4.0*z*z)/(15.0+9.0*z*z); if(y<=0) return z-M_PI; else return z; } } } Testing it with: maxe=0; for(i=0;i<10001;i++) for(j=0;j<10001;j++) { x=0.0001*i-0.5; y=0.0001*j-0.5; a1=atan2(y,x); a2=atan2x(y,x); e=fabs(a1-a2); if(e>maxe) maxe=e; } printf("maxe: %g\n", maxe); gives maxe: 0.0062685 Timing (Xeon 5160 3 GHz) for the loop with only one function call: atan2: 7.028u atan2x: 1.586u I'm not an expert on this, so use it at your own risk :-) Mogens -- Mogens Kjaer, Carlsberg A/S, Computer Department Gamle Carlsberg Vej 10, DK-2500 Valby, Denmark Phone: +45 33 27 53 25, Fax: +45 33 27 47 08 Email: mk(a)crc.dk Homepage: http://www.crc.dk

On Sat, Sep 01, 2007 at 11:23:38 +0100, Chris Jones <jonesc(a)hep.phy.cam.ac.uk&gt; wrote: Just for the record (i.e. you might already know this) you can do taylor expansions around any point, not just the places where you get the nice looking series. This is still a significant restriction, but in cases where you know you will be doing all of your evaluations near a specific point, this a possibility.

On Sat, 2007-09-01 at 11:23 +0100, Chris Jones wrote: Hi, Chris, When I need limited precision for such a value (fixed, some limited resolution), I generally implement it with a lookup table. If implemented in assembly, it is about 6 to 10 instructions to retrieve the value, and add a multiply and truncate to get the resolution you need from a floating point result, and an add if you have negative results. This is often faster than even hardware processing for complex values. Regards, Les H

On Tue, 4 Sep 2007, Mike McCarty wrote: The size of the interval over which an approximation is valid usually depends on the desired accuracy whether or not the approximation is Taylor series. That said, given the desire for a maximum error, a polynomial approximation and a given interval, one should probably use Chebyshev polynomials. Better yet would be to use the hardware instructions. According to at least one speed test, cosines and arctangents each take about as long as a dozen multiplications. http://www.obliquity.com/computer/speedtest.html -- Mike hennebry(a)web.cs.ndsu.NoDak.edu "Horse guts never lie." -- Cherek Bear-Shoulders

Thanks for the tips. I'll take a look. Chris

On Fri, Aug 31, 2007 at 03:48:45PM +0100, Chris Jones wrote: The AMD Core Math library ( http://developer.amd.com/acml.jsp ) is free beer, not F/OSS, but you might want to compare your own results to the documented scalar and vector performance of their transcendental function implementations. [Intel has a more extensive library, but I've never looked through their documentation.] Below are some performance figures excerpted from the ACML v3.6 documentation. In addition to the scalar and general vector versions, there are also specialized versions for 2-vectors and 4-vectors. The documentation for each of these functions says: "the routine is accurate to better than 1 ulp over the valid input range." [ulp = unit in the last place, i.e., the trailing bit.] fastexp: fast double precision exponential function : double fastexp (double x) * Performance: + 75 cycles for most valid inputs. vrda_exp: Array double precision exponential function : void vrda_exp (int n, double *x, double *y) * Performance: + 33 cycles per value for valid inputs, n = 24. fastlog: fast double precision natural logarithm function : double fastlog (double x) * Performance: + 97 cycles for most valid inputs. 86 cycles for .97 < x < 1.03 vrda_log: Array double precision natural logarithm : void vrda_log (int n, double *x, double *y) * Performance: + 51 cycles per value for valid inputs, n = 24. Regards, Bill Rugolsky