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.