# How to write precision independent code in MATLAB

by on July 21, 2016

One of the main design goals of toolbox is ability to run existing scripts in extended precision with minimum changes to code itself. Thanks to object-oriented programming the goal is accomplished to great extent. Thus, for example, MATLAB decides which function to call (its own or from toolbox) based on type of input parameter, automatically and absolutely transparently to user:

```>> A = rand(100); % create double-precision random matrix >> B = mp(rand(100)); % create multi-precision random matrix   >> [U,S,V] = svd(A); % built-in functions are called for double-precision matrices >> norm(A-U*S*V',1) ans = 2.35377689561389e-13   >> [U,S,V] = svd(B); % Same syntax, but now functions from toolbox are used >> norm(B-U*S*V',1) ans = 3.01016831776648753720608552494953562e-31```

Syntax stays the same, allowing researcher to port code to multi-precision almost without modifications.

However there are several situations which are not handled automatically and it is not so obvious how to avoid manual changes:

• Conversion of constants, e.g. `1/3 -> mp('1/3')`, `pi -> mp('pi')`, `eps -> mp('eps')`.
• Creation of basic arrays, e.g. `zeros(...) -> mp(zeros(...))`, `ones(...) -> mp(ones(...))`.

In this post we want to show technique on how to handle these situations and write pure precision-independent code in MATLAB, so that no modifications are required at all. Code will be able to run with standard numeric types `'double'/'single'` as well as with multi-precision numeric type `'mp'` from toolbox.

## Unified handling of numeric types

For this purpose we introduce flexible wrapper function which unifies work with different numeric types:

```function r = numeric_t(expression) global class_t; if (nargin > 0) if(strcmpi(class_t,'mp')), r = mp(expression); else if isnumeric(expression) r = expression; else r = eval(expression); end; end; else r = class_t; end; end```

The function produces output of numeric type specified by global variable `'class_t'`.
Now we can evaluate literal expressions, create constants and basic arrays using the same syntax for `'double'` and `'mp'` types:

```>> global class_t;   >> class_t = 'double';   >> x = numeric_t('pi/2') x = 1.5707963267949   >> y = numeric_t('eps') y = 2.22044604925031e-16   >> B = zeros(10,numeric_t); >> A = magic(2,numeric_t);   >> whos Name Size Bytes Class Attributes A 2x2 32 double B 10x10 800 double x 1x1 8 double y 1x1 8 double```
```>> class_t = 'mp';   >> x = numeric_t('pi/2') x = 1.5707963267948966192313216916397514   >> y = numeric_t('eps') y = 1.92592994438723585305597794258492732e-34   >> B = zeros(10,numeric_t); >> A = magic(2,numeric_t);   >> whos Name Size Bytes Class Attributes A 2x2 296 mp B 10x10 1832 mp x 1x1 248 mp y 1x1 248 mp```

This technique allows us to write numeric-independent code, which automatically works in `'double'` or extended precision depending on value of `'class_t'` variable.

## Usage in functions

Same approach can be used inside functions without the need for global variable. Examples on how `class_t/numeric_t` can be used to create numeric-independent functions:

• Type of input parameter defines what numeric type is used inside the function:
```function r = foo(x)   class_t = class(x); function r = numeric_t(expression) if (nargin > 0) if(strcmpi(class_t,'mp')), r = mp(expression); else if isnumeric(expression) r = expression; else r = eval(expression); end; end; else r = class_t; end; end % numeric_t   % Function code relying on 'numeric_t' goes below % ... end```
• Numerical type is explicitly passed as parameter:
```function r = foo(n, class_t)   function r = numeric_t(expression) if (nargin > 0) if(strcmpi(class_t,'mp')), r = mp(expression); else if isnumeric(expression) r = expression; else r = eval(expression); end; end; else r = class_t; end; end % numeric_t   % Function code relying on 'numeric_t' goes below % ... end```

(The code of `'numeric_t'` is the same in last two examples – we repeat it on purpose so that user can easily copy-paste the corresponding function template in his/her code.)

In some simple cases, only passing `'class_t'` as a parameter is enough:

```function [c,r] = inverthilb(n, class_t) % INVERTHILB Numeric-independent inversion of Hilbert matrix   A = hilb(n,class_t); c = cond(A); r = norm(eye(n,class_t)-A*invhilb(n,class_t),1)/norm(A,1); % relative accuracy end```

Now we can use the function with `'double'` or extended precision without changing its code:

```% MATLAB/double precision: >> [c,r] = inverthilb(20,'double') c = 1.5316e+18 r = 2.1436e+11```
```% Advanpix/multiple-precision: >> mp.Digits(50); >> [c,r] = inverthilb(20,'mp') c = 2.4522e+28 r = 9.4122e-25   >> mp.Digits(100); >> [c,r] = inverthilb(20,'mp') c = 2.4522e+28 r = 2.2689e-74```

As expected, high level of precision is required to invert Hilbert matrix accurately.