I regularily deal with a huge number of elementary functions as rounding, transcendentals, exponential and logarithm and was unsatisfied with the average performance. So, I started to provide my own assembler subroutines executing with full floating-point performance on any Intel Pentium II+ or compatible system. On recent Pentium 4 computers, I was able to improve the average performance by two to five times.
This archive is dedicated to any number cruncher using MATLAB 6.0 or newer.
Unpack the folder '@double' from the archive.
builtin(func,arguments)to call the original version of
func. Please note that this works only for built-in functions - this means that
which funcreports "func is a built-in function".
ffunc. Existing MATLAB code will work with the built-in functions unless you call
If you work with MATLAB R13 or newer, you may replace 'xor.dll' by 'xur.dll' for matching the representation of logical values (
The libraries should always reside in a subfolder called '@double' to make sure they are not called for any data except of type
double. They do not check the data type of passed arguments.
Due to the overhead for calling external functions, the built-in functions work faster if called for matrices with less than about 4 to 8 elements. Therefore, if you know the size of your matrix, you can choose the appropriate function.
xor are not built-in but MATLAB scripts. Therefore, the functions within this package work much faster in any case as shown in the following figure.
xur is an even faster version of
xor. It can be used in any situation where its output, a logical (boolean)
uint8 matrix, will not be transformed to
double for computation. It serves in particular as a value selector in a statement like
Figure 1: Benchmark on a 2GHz Pentium 4 Mobile system
The results of the external functions slightly differ from the built-in functions. The accuracy of the external functions benefits from the floating-point registers with a 64bit mantissa on Intel processors. Intermediate values are kept in floating-point registers such that rounding takes place mostly once - when writing the result into the output matrix with a 53bit mantissa.
In general, a statement of
value with a relative error of less than 1.0E-13. The roundoff error of about 1.1E-16 leads to relatively important deviations in exponentiation. Note also that in particular addition/subtraction in the argument of a logarithm are critical operations due to a relative amplification of the rounding error. The logarithm itself works accurate over the full complex plane R x iR. For the sake of performance, the inverse transcendental functions are currently not implemented in that way. See also the summary about complex functions.
An exception is made for every periodic transcendental function, where the constant 2pi is truncated to a 53bit mantissa. This guarantees that a statement of the form
sin(x) with real
x always produces the expected value
sin(rem(x,2*pi)) at the same accuracy. The reminder is computed explicitly and prescaled to a 64bit mantissa.
x=pow2(pi,0:256); plot(x,builtin('sin',x),'r',x,sin(x),'b'); set(gca,'XScale','log');
Copyright © Marcel Leutenegger, 2003-2007, École Polytechnique Fédérale de Lausanne (EPFL), Laboratoire d'Optique Biomédicale (LOB), BM - Station 17, 1015 Lausanne, Switzerland. This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation; version 2.1 of the License. This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with this library; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. See '/FILES' in the source archive for a list of the original package contents.
Any warranty is strictly refused. Don't rely on any financial or technical support in case of malfunction or damage.
Comments are welcome. I will try to track reported problems and fix bugs.
Bug fixed in
rem(matrix,value): the result was stored back to the
value causing an assertion failure.
Service release thanks to a bug reported by Tom Minka in
exp(matrix): the output was NaN for infinite input.
This bug fix made me think about affine inputs. They are now all handled as particular values for two reasons:
Service release. Bux fixed in [c,s]=cis(matrix): the second output was set to an invalid dimension. Version information included.
Source code released under GNU Lesser General Public License (LGPL) version 2.1.
Optimized routine for calculating the exponential according to Agner Fog, "Optimizing subroutines in assembly language," at Copenhagen University College of Engineering.
Dimensions of (empty matrix*scalar) matched on MATLAB's behaviour. Thanks to Paolo Bardella for reporting the issue.
Bug report by Tom Minka: with MATLAB 7, the
double functions are also called for
sparse input returning a broken output. As a workaround, you may rename the functions or include/exclude them in/from the MATLAB path.
function(empty matrix,scalar) do not match with MATLAB's behaviour.
Downloading these files, you accept the copyright terms.
MATLAB is a registered trademark of The MathWorks, Inc. Pentium II is a registered trademark of Intel Corporation. Other product or brand names are trademarks or registered trademarks of their respective holders.
© 2011 École Polytechnique Fédérale de Lausanne