EPFLSTILOBMarcel LeuteneggerAbout

MATLAB toolbox

Double class

MATLAB toolbox

Data types

Double class
Extended class
Single class


Flex99-12C correlator
Flex02-0xD correlator
Vector functions


C-MEX wrapper
FPU instructions


Error function
Hankel transform


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.

  1. The subfolder '@double' contains libraries that directly overload the built-in functions. This has the advantage, that existing MATLAB code automatically benefits from the increased performance. If '@double' is in a folder in your MATLAB path, all MATLAB functions benefit. If it is in a 'private' project folder, the project functions benefit. It is save to use only a subset of the class functions. Just say builtin(func,arguments) to call the original version of func. Please note that this works only for built-in functions - this means that which func reports "func is a built-in function".
  2. You may rename the functions func into ffunc. Existing MATLAB code will work with the built-in functions unless you call ffunc.

If you work with MATLAB R13 or newer, you may replace 'xor.dll' by 'xur.dll' for matching the representation of logical values (uint8).


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.

The functions angle, mod and 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.

The function 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 values=matrix(xur(a,b)).

Benchmark of elementary MATLAB functions on a 2GHz Pentium 4 Mobile

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 inverseFunction(function(value)) produces 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.

Transcendental functions: pi

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.

Example 1:



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.


January 18, 2004

Initial release

February 28, 2004

Bug fixed in rem(matrix,value): the result was stored back to the value causing an assertion failure.

June 20, 2004

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:

  1. The output is well defined. In cases with more than one possible solution, the function limit towards that value has been used.
  2. The performance does not degreade but increases considerably (table look-up instead of calculation). Any floating-point operation producing an affine result tries to throw an exception. Even if the exception is masked as within MATLAB, the processor calls up an internal assist slowing down the computation to about 10%-20% of normal performance.
May 2, 2005

Service release. Bux fixed in [c,s]=cis(matrix): the second output was set to an invalid dimension. Version information included.

June 28, 2007

Source code released under GNU Lesser General Public License (LGPL) version 2.1.

May 12, 2008

Optimized routine for calculating the exponential according to Agner Fog, "Optimizing subroutines in assembly language," at Copenhagen University College of Engineering.

September 17, 2008

Dimensions of (empty matrix*scalar) matched on MATLAB's behaviour. Thanks to Paolo Bardella for reporting the issue.

Known bugs

May 23, 2005

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.

September 17, 2008

Dimensions of 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