EPFLSTILOBMarcel LeuteneggerAbout MATLAB toolbox Double class |
---|

MATLAB toolbox

**Double class**

Extended class

Single class

Flex99-12C correlator

Flex02-0xD correlator

Vector functions

MEX-Builder

- Introduction
- Requirements
- Installation
- Comments
- Copyright
- Warranty
- History
- Known bugs
- Download
- Trademarks
- x86 Linux port

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.

- An Intel-based computer architecture with an Intel Pentium II or better processor.
- MATLAB 6.0 or newer.

Unpack the folder '@double' from the archive.

- 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,`

to call the original version of`arguments`)`func`

. Please note that this works only for built-in functions - this means that`which func`

reports "func is a built-in function". - 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))`

.

**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.

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(`

with real `x`)

always produces the expected value `x``sin(rem(`

at the same accuracy. The reminder is computed explicitly and prescaled to a 64bit mantissa.`x`,2*pi))

**Example 1:**

```
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.

Initial release

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:

- The output is well defined. In cases with more than one possible solution, the function limit towards that value has been used.
- 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.

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.

Dimensions of `function(empty matrix,scalar)`

do not match with MATLAB's behaviour.

Downloading these files, you accept the copyright terms.

- Double class and source code for MATLAB 6.0 or newer

Digital signature

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