Version History

August 1, 2019 –

  • Functions ISHERMITIAN, ISBANDED, BANDWIDTH and similar have been improved. Now all are (much) faster especially for large matrices.
  • Speed of all functions related to sparse matrices have been increased. This is because we switched to our new engine for computations with sparse matrices.
    The most notable speed-up is in sparse LU – decomposition is 2-3 times faster now, solver is up to 10 times depending on matrix size & structure.
  • New functions MP.GETWORDS64 and MP.SETWORDS64 for low-level access to quadruple number/array have been added.
    Toolbox represents quadruple numbers in 128-bit IEEE-754 format which can be considered as 2 sequentially stored 64-bit unsigned integers (high and low part). Function MP.GETWORDS64 allows to extract these integers, whereas MP.SETWORDS64 generates quadruple number/array from integer parts:

    >> x = mp('pi')
    x = 
    >> format hex
    >> [h,l] = mp.getwords64(x)   % Get the high & low integer parts of quadruple.
    h =
    l =
    >> y = mp.setwords64(h,l)  % Generate quadruple number from integer parts.
    y = 

    These functions can be used for cross-platform data exchange or to interface with other quadruple-enabled codes.

  • Added function to compute KOORNWINDER orthogonal polynomials on triangle.

March 28, 2019 –

  • Functions MAX and MIN have been completely re-implemented. Now both are faster and support extended set of options – 'omitnan', 'includenan' and 'all'. Summary on currently supported variants:
    C = max(A,B)
    C = max(A,B,nanflag)
    M = max(A)
    M = max(A,[],dim)
    M = max(A,[],nanflag)
    M = max(A,[],dim,nanflag)
    [M,I] = max(___)
    M = max(A,[],'all')
    M = max(A,[],'all',nanflag)

    At this moment, option 'vecdim' is not yet supported.

  • Added CUMMAX and CUMMIN functions.
  • All cumulative functions CUMPROD, CUMSUM, CUMMAX and CUMMIN have been enabled with the support for 'omitnan', 'includenan', 'forward' and 'reverse' options. Thanks to Abe Ellison for request.
  • Fixed incompatibility with older syntax of RAND and RANDN. Although it is now considered obsolete, some old code is still using it. Thanks to  Nick Higham for requesting this feature.

February 25, 2019 –

  • Added implicit expansion/broadcasting of dimensions to following operations:
    +       -       .*      ./
    .\      .^      <       <=
    >       >=      ==      ~=
    |       &       xor
    min     max     mod     rem
    hypot   atan2

    This feature can be enabled/disabled by special command mp.EnableImplicitExpansion:

    >> mp.EnableImplicitExpansion(true);
    >> mp([1 2 3]).*mp([1;2;3])
    ans = 
            1        2        3    
            2        4        6    
            3        6        9    
    >> mp.EnableImplicitExpansion(false);
    >> mp([1 2 3]).*mp([1;2;3])
    Error using  .*  (line 1677)
    Matrix dimensions must agree
  • Greatly improved load balancing among computational threads. Expected speed-up is up to 25% in all “heavy” matrix operations (depending on number of CPU cores, CPU load and operation).
  • Numerous smaller fixes and improvements.

January 17, 2019 –

  • Version 4.5.3 is available for all platforms now. Update is strongly recommended.

December 7, 2018 –

  • Numerical stability of EXPM has been improved (=higher accuracy in working with highly ill-conditioned matrices).
  • All service routines has been reviewed and improved: mp.NumberOfThreads, mp.GuardDigits, mp.ExtendConstAccuracy, mp.FollowMatlabNumericFormat and mp.OverrideDoubleBasicArrays.
  • Added support for multidimensional slice assignments with same number of elements but different shapes:
    a = zeros(10,10,10,'mp');
    b = ones(10,10,10,'mp');
    a(:,:,1) = b(:,1,:);

    Thanks to Taras Plakhotnik for requesting this feature.

November 1, 2018 –

  • Version 4.5.2 is available for all platforms now. Update is strongly recommended.

October 23, 2018 –

  • Stability and convergence of divide & conquer algorithm in SVD has been improved. Thanks to Tao Meng from Peking University for reporting the issue.

September 21, 2018 –

  • Major update of toolbox including new optimized code for operations in small to medium precision range.
    Now speed difference between true quadruple(34) and other precision of comparable level is reduced.  

    Timings of EIG(A), with random unsymmetric matrix of 200×200 size:

    % (before)
    20 digits, time = 6.01 sec
    25 digits, time = 6.26 sec
    30 digits, time = 6.82 sec
    34 digits, time = 3.39 sec
    35 digits, time = 7.07 sec
    40 digits, time = 9.35 sec
    45 digits, time = 9.33 sec
    50 digits, time = 9.74 sec
    55 digits, time = 10.3 sec
    % (now)
    20 digits, time = 4.51 sec
    25 digits, time = 4.77 sec
    30 digits, time = 5.06 sec
    34 digits, time = 3.20 sec
    35 digits, time = 5.24 sec
    40 digits, time = 7.03 sec
    45 digits, time = 7.01 sec
    50 digits, time = 7.27 sec
    55 digits, time = 7.81 sec

    The code to reproduce the timings:

    for p=[20:5:30 34 35:5:55]
        A = randn(200,'mp');
        s = tic; [V,D] = eig(A); e = toc(s);
        fprintf('%d digits, time = %.3g sec\trel.error = %g\n', p, e, norm(A*V - V*D,1)/norm(A,1));
        clear A V D e;

    Similar performance gain can be seen in all operations in arbitrary precision mode.

July 25, 2018 –

  • The CIRCSHIFT has been updated to support 3rd argument. This syntax was introduced in recent MATLAB versions. Thanks to Patrick Dorey for requesting this update.

April 18, 2018 –

  • Now SORT treats NaN as the largest value in array. This is needed for one-to-one compatibility with MATLAB. Thanks to Massimiliano Fasi for reporting this issue.

March 27, 2018 –

  • Fixed incompatibility with R2018a on all platforms. Our deepest thanks to Andre Van Moer, Michal Kvasnicka, Enrico Onofri and Siegfried Rump for their help with tests on various systems and MATLAB versions.
  • Improved QZ decomposition to avoid premature exit in case of extremely high precision requested. Thanks to Bor Plestenjak.
  • Switched to new TLS model on GNU Linux. This might help for some GNU Linux installations where famous TLS-issue still shows up. Thanks to Denis Tkachenko for detailed tests.

March 15, 2018 –

  • Added ERFINV, ERFCINV and NORMINV routines in quadruple precision.
  • Improved ORDQZ to avoid false positive alarm on numerical instability. Thanks to Denis Tkachenko.

December 25, 2017 –

  • Fast Fourier Transform (FFT) routines have been upated with various speed optimizations, including multi-core parallelism, extended set of small FFT with minimum number of arithmetic operations, etc.Now toolbox uses split-radix for power-of-two, mixed-radix Stockham framework for composite and Bluestein algorithm for prime lengths FFT. All algorithms have been optimized for parallel execution and quadruple/multi-precision modes.

    >> mp.Digits(34);
    >> x = rand(2^23,1,'mp');  % n = 8M/quadruple
    >> tic; y = fft(x); toc;
    Elapsed time is 16.759009 seconds.
    >> tic; z = ifft(y); toc;
    Elapsed time is 29.263380 seconds.
    >> mp.Digits(100);
    >> x = rand(2^20,1,'mp');  % n = 1M/100-digits
    >> tic; y = fft(x); toc;
    Elapsed time is 14.243092 seconds.
    >> tic; z = ifft(y); toc;
    Elapsed time is 21.368476 seconds.


    >> mp.Digits(34);
    >> x = rand(2^23,1,'mp');  % n = 8M/quadruple
    >> tic; y = fft(x); toc;   % ~5 times faster
    Elapsed time is 3.068695 seconds.
    >> tic; z = ifft(y); toc;  % ~6 times faster
    Elapsed time is 4.953689 seconds.
    >> mp.Digits(100);
    >> x = rand(2^20,1,'mp');  % n = 1M/100-digits
    >> tic; y = fft(x); toc;   % ~8 times faster
    Elapsed time is 1.753420 seconds.
    >> tic; z = ifft(y); toc;  % ~9 times faster
    Elapsed time is 2.293286 seconds.
  • The COSINT and SININT have been extended for negative arguments. Thanks to Tomoaki Okayama and Ryota Hara.
  • Fixed AIRY for pure imaginary arguments. Thanks to Tom Wallace for detailed bug report.
  • Fixed bug in SUBSAGN reported by Jon Vegard.

October 24, 2017 –

  • Improved storage format for cross-platform compatibility. Multiprecision variables stored in file using standard SAVE/LOAD commands can now be transferred to any platform without loss in accuracy (e.g. from GNU Linux to Windows). Please note, quadruple precision data were always cross-platform.
  • Fixed numerical instability in SYLVESTER_TRI. Thanks to Massimiliano Fasi.
  • Fixed accuracy loss in ELLIPKE and ELLIPJ. Thanks to Denis Silantyev.
  • Fixed spurious complex output for extreme arguments in BESSELJ/Y. Thanks to Taras Plakhotnik.
  • Fixed premature overflow/underflow in BESSELI for corner case arguments.

October 3, 2017 –

  • Added extended set of routines to compute orthogonal polynomials:
    chebyshevV Chebyshev polynomials of the third kind
    chebyshevW Chebyshev polynomials of the fourth kind
    laguerreL Generalized Laguerre polynomials
    gegenbauerC Gegenbauer (ultraspherical) polynomials
    jacobiP Jacobi polynomials
    zernikeR Radial Zernike polynomials

    All routines are optimized for parallel execution, quadruple and multi-precision modes, e.g.:

    >> x = randn(500);
    >> tic; laguerreL(10,vpa(x)); toc;
    Elapsed time is 308.504995 seconds.
    >> tic; laguerreL(10,mp(x)); toc;  % ~ 2500 times faster
    Elapsed time is 0.120173 seconds.
  • Improved stability of divide & conquer algorithm for SVD. Thanks to Bartosz Protas for detailed bug report.
  • Various small changes for better compatibility with plotting engine of MATLAB.

August 22, 2017 –

    1. Parallelism. Modern processors usually run two logical threads per one physical CPU core (hyper-threading). This improves performance in situations when user runs several independent tasks (e.g. usual desktop applications). However this leads to sub-optimal results in case of heavy-load scientific computations, when multi-threaded algorithms are specifically designed to be well-balanced and exhaust all available resources of each core.To alleviate the issue, now toolbox binds its threads one-to-one to physical cores on CPU, ignoring hyper-threaded (logical) cores. You might see performance boost up to 15% depending on algorithm and your CPU:
      >> mp.Digits(100);
      >> A = randn(500,'mp');
      % 4.3.6
      >> tic; [S,V,D] = svd(A); toc;
      Elapsed time is 33.177931 seconds.
      % 4.4.1
      >> tic; [S,V,D] = svd(A); toc;
      Elapsed time is 28.520263 seconds.

      This is experimental feature, implemented only in Windows version of toolbox. Other platforms will follow.

    2. Bessel functions. Whole family of Bessel and Hankel functions have been revised. Majority of improvements are related to the case of non-integer orders, but other cases have been polished too:
      • Fixed accuracy loss for half-integer orders when |v| >> |z|. Thanks to Mert Kalfa for detailed tests and reports.
      • Fixed accuracy loss in modified Bessel functions in case of pure imaginary arguments and large orders. Now we apply several different asymptotic approximations depending on parameters (one form is specifically tuned for imaginary arguments).
      • Added fast implementation for Hankel functions for real arguments and integer orders.
      • Added fast quadruple implementation for modified Bessel functions (real arguments and orders).

      Besides, now all Bessel functions rely on parallel execution, for all kinds of arguments.

    3. Miscellaneous.
      • Now toolbox includes pre-computed coefficients of Gauss-type quadrature for up to 512 order and up to 250 digits of precision. Which means retrieving quadrature nodes & weights now costs O(1)! All usual quadrature are covered – Legendre, Hermite, Gegenbauer, Jacobi and Laguerre.
        % 4.3.6
        >> mp.Digits(34);
        >> tic; [x,w] = mp.GaussJacobi(256); toc;
        Elapsed time is 3.265601 seconds.
        >> mp.Digits(100);
        >> tic; [x,w] = mp.GaussJacobi(256); toc;
        Elapsed time is 14.518643 seconds.
        % 4.4.1
        >> mp.Digits(34);
        >> tic; [x,w] = mp.GaussJacobi(256); toc;
        Elapsed time is 0.000994 seconds.
        >> mp.Digits(100);
        >> tic; [x,w] = mp.GaussJacobi(256); toc;
        Elapsed time is 0.002819 seconds.

        This improvement was inspired by communication with Massimiliano Fasi (Padé approximation for LOGM) and Enrico Onofri. We hope this will be helpful for various spectral methods and approximations.

      • New platform-independent storage format for multiprecision arrays of non-quadruple precision. Multiprecision variables stored in file using standard SAVE/LOAD commands can now be transferred to any platform without loss in accuracy (e.g. from GNU Linux to Windows). Please note, quadruple precision data were always cross-platform.
      • Fixed critical bug in computing Bernoulli coefficients (widely used for evaluation of some special functions). The bug caused race conditions in multi-threaded computations with possible MATLAB crash.

The large number of changes allows us to bump the version directly to 4.4.1 bypassing smaller revisions.

June 6, 2017 –

  • Improved accuracy of matrix logarithm – LOGM. Now it has better handling of near-singular matrices and refined Pade approximation. Thanks to Massimiliano Fasi who reported the issues and provided test cases.
  • The SORT function in toolbox has been enabled with undocumented features of MATLAB’s built-in SORT. Thanks to Daniele Prada from “Istituto di Matematica Applicata e Tecnologie Informatiche”, in Pavia, Italy.

May 12, 2017 –

  • Speed-up of EIGS routine. The heaviest part (modified Gram-Schmidt) has been moved to C++. Performance improvement is 4-8 times depending on a problem.
  • Various bug fixes in EIGS, ORSCHUR and in mixed-precision computations (especially with sparse matrices).

March 28, 2017 –

  • Added mp.NumberOfThreads(N) function to control how many threads toolbox can use in computations enabled with multi-core parallelism. Be default, toolbox uses all available CPU cores, which might not be optimal for all cases. Now user can adjust this flexibly. Requested by Clay Thompson.
    >> A = rand(1000,'mp');
    >> mp.NumberOfThreads(1);
    >> tic; exp(A); toc;
    Elapsed time is 1.297683 seconds.
    >> mp.NumberOfThreads(4);
    >> tic; exp(A); toc;
    Elapsed time is 0.425246 seconds.
  • Added support for second input argument in ROUND function. Requested by Michael Klanner.
  • Added support for second output argument in LINSOLVE function. Requested by Zhaolong Hu.
  • Added optimization to BESSELK for half-integer orders.
  • Fixed “freeze” bug in BESSELJN for high-orders and real quadruple arguments, only Windows version was affected. Reported by Maxime Dana.
  • Improved integration with warning system in MATLAB. In some cases toolbox showed warnings ignoring the fact that they are disabled in MATLAB. This is ongoing work.

January 25, 2017 –

  • Added ACCUMARRAY, CELL2MAT and MP.CELLFUN routines.
  • Fixed bug related to empty arrays in BSXFUN, reported by Vladimír Šedenka.
  • Added special functions-overloads to retrieve machine epsilon, smallest/largest floating numbers for a given precision:
    % Before
    >> mp.eps
    >> mp.eps(mp(1))
    >> mp.realmin
    >> mp.realmax
    % Now /
    >> eps('mp')
    >> eps(mp(1))
    >> realmin('mp')
    >> realmax('mp')

    We have replaced MP.EPS, MP.REALMIN and MP.REALMAX with functions without 'MP.' prefix. This is important change needed to write precision independent code.

January 4, 2017 –

  • Added EIGS routine for computing subset of eigenvalues and eigenvectors. All features are supported: standard and generalized eigenproblems, matrix and function handle inputs and various parts of the spectrum – 'SM', 'LM', 'SA', 'LA', 'SI', 'LI', 'SR' and 'LR'.
  • Added FFTN/IFFTN – routines for multi-dimensional Fourier transformation.
  • All Fourier-related routines have been updated with specific optimizations for quadruple precision. Expected speed-up is 2-4 times.
  • Fixed several bugs related to empty arrays in assignment and relational operators, reported by Clay Thompson.

December 6, 2016 –

  • Load balancing in multi-threading pool has been improved (especially on Windows).
  • ODE15S has been updated to include support for JPattern option, requested by Nicola Courtier.
  • Fixed crash in case when inputs to ORDSCHUR and ORDQZ have different complexity.

November 8, 2016 –

  • All elementary and special mathematical functions (exp, sqrt, sin, cos, gamma, bessel, etc. ) have been updated with parallel execution on multi-core CPUs. Speed-up is proportional to number of CPU cores.
    For example, timing reduction for Bessel Y0 on Core i7 990x (6 hardware cores):

    % Before
    >> A = rand(2000,'mp');
    >> tic; bessely(0,A); toc;
    Elapsed time is 4.278553 seconds.
    % Now / 4.3.0
    >> A = rand(2000,'mp');
    >> tic; bessely(0,A); toc;
    Elapsed time is 0.678716 seconds.
  • Bug fixes: memory leak in complex arithmetic when used in multi-threading environments. Quadruple precision mode wasn’t affected by the bug.

November 2, 2016 –

  • The right matrix division (MRDIVIDE) has been updated to match the speed of the left matrix division (MLDIVIDE). Reported by Massimiliano Fasi.

October 26, 2016 –

  • The eigensolver EIG(A[,B]) has been updated and optimized further. One of the important changes is that now we have special ultra-fast solver for (k,p)-diagonal matrices in case of generalized eigenproblem. In some cases our quadruple precision solver is faster than built-in double precision solver in MATLAB! Please see the bottom of the article for example – Architecture of eigenproblem solver.
  • Least square linear solver has been speeded up by at least 2 times (2000x1500). Actual speed-up grows with matrix size and number of CPU cores. Please note that now we use Rank-Revealing QR (RRQR) instead of SVD.
  • Bug fixes: memory leaks in generalized eigensolver (complex inputs), SVD for scalar inputs (reported by Denis Tkachenko) and MIN/MAX now ignore NaN values for sure (reported by Gerard Ketefian).

October 5, 2016 –

  • The linear system solver MLDIVIDE has been updated and optimized further. Now every decomposition in the solver (LU, LDLT and Cholesky) run faster by 10-25%. The final results and algorithm are outlined in our recent article Architecture of linear systems solver.
  • Computation of eigenvalues of generalized symmetric eigenproblem has been speeded up by 2 times. This is the result of refined parallel algorithm for our D&C solver.
  • Also we have added MRRR solver for generalized symmetric eigenproblem [V,D] = eig(A,B). It is slightly faster than D&C for computing eigenvectors, but most importantly – it provides better accuracy in some cases.
  • We have added special algorithm to solve n-diagonal/banded standard eigenproblems – eig(A). We are preparing article with detailed outline of our poly-algorithm for standard eigenproblems (similar to MLDIVIDE).

September 11, 2016 –

  • MLDIVIDE (linear system solver) has been updated with specialized solver for banded matrices.
    Timings have been substantially improved for such matrices, e.g. 5Kx5K pentadiagonal matrix:

    % Before:
    >> n = 5000;
    >> A = mp(diag(1:n,0)+diag(1:n-1,1)+diag(1:n-2,2)+diag(-(1:n-1),-1)+diag(-(1:n-2),-2));
    >> b = mp(rand(n,1));
    >> tic; x = A\b; toc;               % LU+pivoting, O(n^3)
    Elapsed time is 257.355805 seconds.
    >> norm(A*x-b,1)/norm(A,1)
    ans = 
    % Now:
    >> n = 5000;
    >> A = mp(diag(1:n,0)+diag(1:n-1,1)+diag(1:n-2,2)+diag(-(1:n-1),-1)+diag(-(1:n-2),-2));
    >> b = mp(rand(n,1));
    >> tic; x = A\b; toc;               % Specialized LU+pivoting within bandwidth, O(n*p*q)
    Elapsed time is 0.233104 seconds.   % ~x1000 times faster
    >> norm(A*x-b,1)/norm(A,1)
    ans = 

    The speed-up comes from the fact that now algorithms works with only few non-zero diagonals, instead of crunching full matrix.

  • Positive definite banded matrices have received specialized solver as well (only superdiagonals are used = less computations).

MLDIVIDE is a poly-algorithm which selects the best solver depending on matrix properties, basically we have rewritten it from scratch lately. Now it has (much)faster analysis and full set of specialized solvers + rcond estimators for dense matrices.

August 27, 2016 –

  • MLDIVIDE (linear system solver) has been updated with specialized solver for tridiagonal matrices.
    Now computation complexity for the case is just O(n) instead of O(n^3):

    % Before:
    >> n = 2000;
    >> A = mp(diag(1:n,0)+diag(1:n-1,1)+diag(1:n-1,-1));
    >> b = mp(rand(n,1));
    >> tic; x = A\b; toc;             % LU+pivoting, O(n^3)
    Elapsed time is 27.767122 seconds.
    >> norm(A*x-b,1)/norm(A,1)
    ans = 
    % Now:
    >> n = 2000;
    >> A = mp(diag(1:n,0)+diag(1:n-1,1)+diag(1:n-1,-1));
    >> b = mp(rand(n,1));
    >> tic; x = A\b; toc;              % Specialized tridiagonal LU+pivoting, O(n)
    Elapsed time is 0.079523 seconds.  % ~350 times faster
    >> norm(A*x-b,1)/norm(A,1)
    ans = 
  • SUBSASGN and DIAG have been improved to be compatible with MATLAB in special cases:
    % Column to row assignment:
    >> A = magic(3,'mp');
    >> x = [1;2;3];
    >> A(1,:) = x
    A = 
            1        2        3    
            3        5        7    
            4        9        2
    % Building matrix from empty diagonal:        
    >> diag(rand(1,0,'mp'),0)
    ans =
    >> diag(rand(1,0,'mp'),1)
    ans = 
    >> diag(rand(1,0,'mp'),2)
    ans = 
            0        0    
            0        0
  • Added functions for orthogonal polynomials: legendreP, chebyshevT, chebyshevU and hermiteH.Variable Precision Arithmetic (VPA)/MATLAB 2016a:
    >> N = 1000;   % polynomial degree
    >> M = 5000;   % number of points
    >> digits(25);
    >> x = vpa(rand(M,1));
    >> tic; v = chebyshevT(N,x); toc;
    Elapsed time is 2.600266 seconds.
    >> tic; v = chebyshevU(N,x); toc;
    Elapsed time is 1.908803 seconds.
    >> tic; v = hermiteH(N,x);   toc;
    Elapsed time is 34.765995 seconds.
    >> tic; v = legendreP(N,x);  toc;
    Elapsed time is 44.810674 seconds.

    Advanpix Multiprecision Toolbox:

    >> N = 1000;   % polynomial degree
    >> M = 5000;   % number of points
    >> mp.Digits(34);
    >> x = mp(rand(M,1));
    >> tic; v = chebyshevT(N,x); toc;
    Elapsed time is 0.426981 seconds.    % x6 times faster
    >> tic; v = chebyshevU(N,x); toc;    % x4 times faster
    Elapsed time is 0.450537 seconds.
    >> tic; v = hermiteH(N,x);   toc;    % x53 times faster
    Elapsed time is 0.653591 seconds.
    >> tic; v = legendreP(N,x);  toc;    % x43 times faster
    Elapsed time is 1.035238 seconds.

    Please note, our routines are not multi-core optimized (yet). In due course timings will be divided by number of CPU cores.

August 24, 2016 –

  • Added following new functions (by categories):Linear Equations
    SYLVESTER, LYAP, DLYAP – solvers for linear matrix equations of Lyapunov and Sylvester.
    Continuous, discrete-time and generalized forms are covered.  

    Matrix Decomposition
    Also GSVD for generalized singular value decomposition.

    Matrix Analysis

    All new functions (with few exceptions) are implemented in toolbox core using C++ for better performance.

  • Following functions have been revised and improved:
    CROSS and DOT now support multi-dimensional arrays.
    QR, CHOL, SVD and KRON now have better handling of corner cases (e.g. empty inputs, error reports, etc.)
  • Added support for differential algebraic equations (DAEs) to ODE15S.

The number of new functions and changes allow us to bump the version directly to 4.1 bypassing smaller revisions.

August 11, 2016 –

  • Improvements to linear least-squares solver. Instead of SVD, now we use rank-revealing QR decomposition (RRQR), heavily optimized for parallel execution on multi-core CPUs:
    % Before
    >> mp.Digits(34);
    >> A = rand(1000,500,'mp')-0.5;
    >> tic; x = A\A; toc;
    Elapsed time is 45.334321 seconds.
    % Now -
    >> mp.Digits(34);
    >> A = rand(1000,500,'mp')-0.5;
    >> tic; x = A\A; toc;
    Elapsed time is 12.059790 seconds.

    All cases are covered – real, complex, quadruple and arbitrary precision. Speed-up ratio scales with matrix size and number of CPU cores.

July 6, 2016 –

  • New functions – NEXTABOVE and NEXTBELOW have been added. They generate the next representable floating-point value towards positive/negative infinity:
    >> mp.Digits(34);
    >> nextabove(mp(1))
    ans = 
    >> nextbelow(mp(1))
    ans = 
    >> nextabove(mp(-1))
    ans = 
    >> nextbelow(mp(-1))
    ans = 

    The routines are able to generate all/any floating-point numbers representable in given precision, thus they are indispensable for accuracy checks of various algorithms. In particular, to investigate quality of approximation in close vicinity of roots or singularities. As an example, here is quick accuracy check of MATLAB's built-in gamma function:

    >> mp.Digits(34);
    >> mp.FollowMatlabNumericFormat(true);
    >> format longE;
    >> x = nextabove(-1)  % closest value to the pole
    x =
    >> gamma(mp(x))       % correct function value in quadruple precision
    ans = 
    >> gamma(x)           % MATLAB gives no correct digits at all!
    ans =
  • Occasional crashes on Mac OSX have been fixed. As it turned out, sometimes MATLAB (especially older versions) forget to delete MEX module from memory on “clear” command, even though at-exit handlers were called. This leaves MEX in uninitialized and unusable state! Next attempt to use any command from such MEX results in crash.Now we do the unloading procedure manually on all platforms to avoid this from happening again.
  • Prediction on maximum number of iterations needed to reach target accuracy level in Schur & SVD algorithms have been revisited. Now we rely on more pessimistic assumption to make sure algorithms converge even if it requires much higher number of iterations.

June 21, 2016 –

  • Matrix exponential – EXPM has been switched to use classic scaling and squaring algorithm. Although Schur-Parlett[Higham2003] might give more accurate results in some cases [Higham2009], it is approx. 5-10 times slower. Thus we decided to use fast scaling & squaring.
  • Eigen-decomposition re-ordering functions: ORDSCHUR, ORDEIG and ORDQZ have been updated with proper error handling in corner cases.

June 3, 2016 –

  • Whole set of numerical integration routines has been refreshed: INTEGRAL, INTEGRAL2, INTEGRAL3, QUAD2D, TRAPZ and CUMTRAPZ. The INTEGRALx routines are based on nested Gauss-Kronrod quadrature rule or its product for 2-3D.For some reason, MATLAB's built-in INTEGRAL doesn’t return the error bound. We consider this unacceptable and our version fixes the flaw:
    >>f = @(x)sin((1:5)*x);
    >>[q, errbnd] = integral(f,mp(0),mp(1),'ArrayValued',true,'RelTol',mp('eps*100'));
    >>[q' errbnd'] 
    ans =
       0.4596976941318602825990633925570232        2.422586612556771991014735044455665e-34    
       0.7080734182735711934987841147503806        3.709251321227161864242960518419781e-34    
       0.6633308322001484857571909315770862        3.510290962167396105950983616082283e-34    
       0.4134109052159029786597920457744374        2.247449605195361383728053642397239e-34    
       0.1432675629073547471066721656972884        9.201635067043439859020812453968865e-35

    Integration of vector-valued function with error bound for every component.

May 19, 2016 –

  • Added X = SYLVESTER_TRI(A,B,C) function for solving Sylvester equation A*X + X*B = C.
    The most important case is when A and B are upper quasi-triangular (real Schur canonical form) – needed for computing matrix functions (SQRTM, LOGM, etc.).
  • Speed-up of HYPOT and ATAN2, approx. 10-50 times each:
    >> A = 100*mp(rand(1000,1000)-0.5);
    >> B = 100*mp(rand(1000,1000)-0.5);
    >> tic; atan2(A,B); toc;
    Elapsed time is 18.056319 seconds.
    >> A = 100*mp(rand(1000,1000)-0.5);
    >> B = 100*mp(rand(1000,1000)-0.5);
    >> tic; atan2(A,B); toc;
    Elapsed time is 0.262123 seconds.   % ~70 times faster

May 17, 2016 –

  • Speed-up of SQRTM_TRI and INTERP1, approx. 50-100 times each:
    >> A = triu(100*mp(rand(100,100)-0.5));
    >> tic; sqrtm_tri(A); toc;
    Elapsed time is 0.058902 seconds.
    >> tic; old_sqrtm_tri(A); toc;
    Elapsed time is 3.936160 seconds.

    New SQRTM_TRI is implemented in C++, old one was implemented using MATLAB:

    function R = old_sqrtm_tri(T)
    %SQRTM_TRI Square root of upper triangular matrix.
       n = length(T);
       R = mp(zeros(n));
       for j=1:n
           R(j,j) = sqrt(T(j,j));
           for i=j-1:-1:1
               R(i,j) = (T(i,j) - R(i,i+1:j-1)*R(i+1:j-1,j))/(R(i,i) + R(j,j));

    Square root of upper triangular matrix is important for SQRTM and LOGM functions.

May 3, 2016 –

  • Speed-up of some (basic) special functions: gamma, gammaln, erf and erfc:
    % Before/
    >> mp.Digits(34);
    >> A = mp(10*(rand(1000)-0.5));
    >> tic; B = gamma(A); toc;
    Elapsed time is 73.195803 seconds.
    >> tic; B = erf(A); toc;
    Elapsed time is 31.863047 seconds.
    % Now/
    >> mp.Digits(34);
    >> A = mp(10*(rand(1000)-0.5));
    >> tic; B = gamma(A); toc;
    Elapsed time is 1.695979 seconds.  % x43 times faster
    >> tic; B = erf(A); toc;
    Elapsed time is 1.769617 seconds.  % x18 times faster

    Please note, in contrast to MATLAB, these functions support all kinds of input arguments (real negative, complex, sparse, etc.):

    % MATLAB
    >> gammaln(-0.25)
    'Error using gammaln'
    'Input must be nonnegative.'
    % Advanpix
    >> gammaln(mp(-0.25))
    ans = 
        1.589575312551185990315897214778783 +     3.141592653589793238462643383279503i
    >> A = mp(sprand(5,5,0.25));
    >> [x,y,s] = find(A);
    >> s = s + (rand(size(s))-0.5)*1i; % Add imaginary part
    >> A = sparse(x,y,s,5,5)
    A =
        (2,1)      0.1066527701805843886262437081313692 +     0.3173032206534329713321085364441387i
        (4,1)    0.004634224134067443934270613681292161 +     0.3686947053635096782642222024151124i
        (3,3)      0.9618980808550536831802446613437496 -     0.4155641544890896765807042356755119i
        (1,5)      0.4426782697754463313799533352721483 -     0.1002173509011035079652174317743629i
        (4,5)       0.774910464711502378065688390051946 -     0.2401295971493457859224918138352223i
    % Advanpix
    >> erf(A)
    ans =
        (2,1)      0.1324883666365941115619861448158018 +      0.3659492942681403578764987541731807i
        (4,1)    0.005990516950461289561597280148562875 +      0.4356625058347425700172815442490943i
        (3,3)       0.903034046175742903946764521185423 -      0.1758911959897686907097927081692329i
        (1,5)       0.472854450257120712106056672396968 -     0.09314858177203069504581523930686257i
        (4,5)      0.7550126399539313384996685128575232 -      0.1480116547266936627042332666803371i
    % MATLAB
    >> erf(double(A))
    'Error using erf'
    'Input must be real and full.'

April 26, 2016 –

  • We have finished several months of optimization work for elementary functions.
    Please see the results and comparisons: Performance of Elementary Functions.
  • Added support for signed zeros as imaginary part of complex numbers to expression evaluator. As a note, we have been supporting signed zeros from the start, with proper handling of branch cuts, etc. Some additional information: Branch Cuts and Signed Zeros in MATLAB.

April 12, 2016 –

  • Added element-wise relational & logical operations for sparse matrices:
    >> A = mp(sprand(5,5,0.2))
    A =
        (1,1)    0.4387443596563982417535498825600371
        (2,2)    0.7951999011370631809114684074302204
        (1,4)    0.3815584570930083962991830048849806
        (1,5)    0.7655167881490023695789659541333094
        (4,5)    0.1868726045543785962976812697888818
    >> A(A<0.5) = 0
    A =
        (2,2)    0.7951999011370631809114684074302204
        (1,5)    0.7655167881490023695789659541333094

April 6, 2016 –

  • Performance of power and related functions has been improved:
    % Before:
    >> mp.Digits(34);
    >> A = mp(rand(2000)-0.5);
    >> B = mp(rand(2000)-0.5);
    >> tic; C = A.^B; toc;
    Elapsed time is 82.506463 seconds.
    >> tic; C = log(A); toc;
    Elapsed time is 31.506463 seconds.
    % Now:
    >> mp.Digits(34);
    >> A = mp(rand(2000)-0.5);
    >> B = mp(rand(2000)-0.5);
    >> tic; C = A.^B; toc;
    Elapsed time is 5.363670 seconds. % x15 times
    >> tic; C = log(A); toc;
    Elapsed time is 0.843535 seconds. % x37 times

    Similarly for all related functions (log2, log10, sqrt, etc.)

  • Added support for sparse logical indices.
  • Improved support for sparse matrices in isequal, isnequal, isnan, isinf, isfinite, sqrt, expm1, log1p and many other basic functions.
  • Fixed spfun in case when function handle has nested calls to another functions.

March 23, 2016 –

  • Added sortrows, few minor issues has been fixed.
  • Fixed small incompatibilities with R2016a.

March 15, 2016 –

  • Added solver for systems of nonlinear equations – fsolve:

    x = fsolve(fun,x0)
    x = fsolve(fun,x0,options)
    x = fsolve(problem)
    [x,fval] = fsolve(fun,x0)
    [x,fval,exitflag] = fsolve(...)
    [x,fval,exitflag,output] = fsolve(...)
    [x,fval,exitflag,output,jacobian] = fsolve(...)

    All features are supported, including sparse formulation, Jacobian and three algorithms of optimization: trust-region-dogleg, trust-region-reflective and levenberg-marquardt.

    function test_fsolve
    x0 = mp([-5; -5]);
    options = optimset('Display','iter', 'TolFun', mp.eps,...
                                           'TolX', mp.eps,...
    [x,fval,exitflag,output] = fsolve(@myfun,x0,options)
        function F = myfun(x)
            F = [2*x(1) - x(2) - exp(-x(1));
                  -x(1) + 2*x(2) - exp(-x(2))];


    &gt;&gt; mp.Digits(34);
    &gt;&gt; test_fsolve
                                            First-Order                    Norm of 
     Iteration  Func-count    Residual       optimality      Lambda           step
         0           3         47071.2        2.29e+04         0.01
         1           6         6527.47        3.09e+03        0.001        1.45207
         2           9         918.374             418       0.0001        1.49186
         3          12         127.741            57.3        1e-05        1.55326
         4          15         14.9153            8.26        1e-06        1.57591
         5          18        0.779056            1.14        1e-07        1.27662
         6          21      0.00372457          0.0683        1e-08       0.484659
         7          24      9.2164e-08        0.000336        1e-09      0.0385554
         8          27     5.66104e-17        8.34e-09        1e-10    0.000193709
         9          30     2.35895e-34         1.7e-17        1e-11    4.80109e-09
        10          33     1.70984e-51        4.58e-26        1e-12    9.80056e-18
        11          36      1.8546e-68        1.51e-34        1e-13    2.63857e-26
    Equation solved.
    fsolve completed because the vector of function values is near zero
    as measured by the selected value of the function tolerance, and
    the problem appears regular as measured by the gradient.
    x = 
    fval = 
    exitflag =
    output = 
           iterations: 11
            funcCount: 36
             stepsize: [1x1 mp]
         cgiterations: []
        firstorderopt: [1x1 mp]
            algorithm: 'levenberg-marquardt'
              message: 'Equation solved.
    Equation solved. The sum of squared function values, r = 1.854603e-68, is less
    than sqrt(options.TolFun) = 1.387779e-17. The relative norm of the gradient of
    r, 6.583665e-39, is less than 1e-4*options.TolFun = 1.925930e-38.
    Optimization Metric                                        Options
    relative norm(grad r) =   6.58e-39           1e-4*TolFun =   2e-38 (selected)
    r =   1.85e-68                             sqrt(TolFun) =  1.4e-17 (selected)
  • Fixes for various minor bugs and other small improvements.

March 8, 2016 –

  • Added light wrappers for text and histogram routines. Now both accept mp-parameters without errors.
  • We have upgraded to new Intel C++ compiler, MPIR, MPFR and MPC. Any compatibility tests and reports would be highly appreciated.

March 3, 2016 –

  • Basic test matrices are enabled with classname='mp' support: hadamard, hilb, invhilb, magic, pascal, rosser and wilkinson.
    function [c,r] = inverthilb(n, classname)
    % INVERTHILB Class-independent inversion of Hilbert matrix
        A = hilb(n,classname);
        c = cond(A);
        r = norm(eye(n,classname)-A*invhilb(n,classname),1)/norm(A,1);

    Invert Hilbert matrix using different levels of precision:

    >> [c,r] = inverthilb(20,'double')
    c =
    r =
    >> mp.Digits(50);
    >> [c,r] = inverthilb(20,'mp')
    c = 
    r = 
    >> mp.Digits(100);
    >> [c,r] = inverthilb(20,'mp')
    c = 
    r = 
  • Added cast function for conversion of mp-entities to/from different data types.
    >> cast(magic(3),'like',mp('pi'))
    ans = 
            8        1        6    
            3        5        7    
            4        9        2
    >> format longE
    >> cast(rand(2,'double'),'mp')
    ans = 
         6.9907672265668596711662985399016179e-01     9.5929142520544430361439935950329527e-01    
         8.9090325253579849551499592053005472e-01     5.4721552996380307121171426842920482e-01
    >> cast(rand(2,'mp'),'double')
    ans =
         1.386244428286791e-01     2.575082541237365e-01
         1.492940055590575e-01     8.407172559836625e-01
  • Added solvers for Riccati equations: care, gcare, dare and gdare.

February 29, 2016 –

  • Added cplxpair and unwrap functions.
  • Performance of multiplication has been improved for the precision levels <= 385 of decimal digits.
    >> mp.Digits(350);
    >> A = mp(rand(500)-0.5); 
    >> tic; A*A; toc;
    Elapsed time is 22.412486 seconds.
    >> tic; A*A; toc;    % x4 times faster
    Elapsed time is 5.212486 seconds.

    Parameters of our arithmetic engine were tuned for old CPUs. Now it is updated for modern architectures.
    Thanks to Massimiliano Fasi who noticed and reported to us performance drop for small precisions.

  • Conjugate pairs computed by generalized eigen-solver are now guaranteed to match exactly (bit-to-bit).
    Before they matched up to machine epsilon in any precision. Reported by Stefan Güttel.
  • Now toolbox overloads global built-in functions for array creation: zeros, ones, eye, nan, inf, rand, randn and randi.The main goal is to embed support for ‘mp’ arrays in MATLAB out-of-the-box, without the need for mp(zeros(...)) wrapper:
    >> A = zeros(3,3,'mp'); % didn't work since built-in zeros doesn't know what 'mp' is.
    >> A = mp(zeros(3,3));  % manual modification was needed. 
    >> A = zeros(3,3,'mp'); % now works as expected
    >> whos A
      Name      Size            Bytes  Class    Attributes
      A         3x3               376  mp

    This change is very important, as it allows much easier conversion of existing scripts to multiprecision.

    In fact, we have added option to override behavior of such basic functions to always produce mp-outputs for floating-point types: double and single. Be careful with such powerful option (it is turned off by default).

    Please run doc mp.OverrideDoubleBasicArrays for more information on its usage, pros and cons.

February 22, 2016 –

  • Added arbitrary-precision overloads for the following functions: orth, rref, airy, beta, betaln, ellipj, ellipke and legendre.
  • The svd, null and norm have been updated to support more special cases.
    >> A = [1 2 3; 1 2 3; 1 2 3];
    >> Z = null(mp(A))
    Z = 
            -0.8728715609439695250643899416624781         0.4082482904638630163662140124509814    
            -0.2182178902359923812660974854156192        -0.8164965809277260327324280249019639    
             0.4364357804719847625321949708312385         0.4082482904638630163662140124509821
    >> norm(A*Z,1)
    ans = 
    >> null(mp(A),'r')
    ans = 
            -2        -3    
             1         0    
             0         1

    While working on MATLAB’s internal scripts we found several cases of undocumented syntax:

    [Q,Z] = svd(...)                 % two-outputs only
    norm(A,'inf')                    % Inf is passed as a string
    [US,TS, Success] = ordschur(...) % three-outputs, with status as last one.

    Now svd and norm in toolbox support these special cases for compatibility.

  • Basic matrix manipulation & analysis functions have been optimized for sparse matrices: diag, triu, tril, norm, max and min.
    >> A = mp(sprand(10000,10000,0.01));
    >> nnz(A)
    ans =
    >> tic; norm(A,1); toc;
    Elapsed time is 0.0531895 seconds.
    >> tic; max(sum(abs(A))); toc;
    Elapsed time is 0.103413 seconds.

    Not too bad for handling 1M of nonzero quadruples in sparse format!

February 15, 2016 –

  • Added LDLT decomposition for Hermitian indefinite matrices (dense). Real and complex cases are supported and optimized for multi-core CPUs.
    >> A = mp(rand(1000)-0.5);
    >> A = A+A';
    >> tic; [L,D,p] = ldl(A,'vector'); toc;
    Elapsed time is 7.412486 seconds.
    >> norm(A(p,p) - L*D*L',1)/norm(A,1)
    ans = 

    Still we consider this as a draft version – since there is further potential for speed-up.

February 9, 2016 –

January 30, 2016 –

  • Added spline, pchip, ppval, mkpp, unmkpp and interp1 routines.
  • Fixed accuracy loss in expm and crash in ordschur when U is real but T is complex matrix. Thanks to Massimiliano Fasi for tests and reports!
  • Indexing engine subsref has been enabled with all the ad-hoc rules of MATLAB in case when the first (and only) index is semi-empty matrix. This is needed to match the MATLAB behavior in rare cases, e.g. when empty matrices are used as indices in operands of arithmetic operations.

January 19, 2016 –

  • Added norm computation for sparse matrices and vectors. All norms are supported except the 2-norm for sparse matrices (since it needs svds). Please use 1, Inf or recommended 'fro' norm for matrices.
  • Added mp.GaussKronrod function for computation of Gauss-Kronrod nodes and weights.
  • Improved accuracy of eig in computing eigenvectors of real symmetric tridiagonal matrices.The method we used previously (inverse iteration) suffers from numerical instability for very ill-conditioned eigenvectors. Now we have switched to implicit QR.

December 12, 2015 –

  • Added concatenation for sparse matrices: cat, vertcat and horzcat.
  • Fixed find in case when no named output parameters is provided.
  • Changed error messages to be shorter and more informative. Default function from MEX API – mexErrMsgTxt shows intimidating messages with little information:
    % mexErrMsgTxt (Before):
    >> A = mp(magic(3));
    >> A(-1,0) = 10
    'Error using mpimpl'
    'Subscript indices must either be real positive integers or logicals.'
    'Error in mp/subsasgn (line 871)'
    '             [varargout{1:nargout}] = mpimpl(171, varargin{:});'
    % Using our custom workaround:
    >> A = mp(magic(3));
    >> A(-1,0) = 10
    'Error using mp/subsasgn (line 871)'
    'Subscript indices must either be real positive integers or logicals.'

    Now error messages two times shorter with all the necessary information.

December 2, 2015 –

  • Fixed rcond for better handling of singular matrices.
  • Improved generalized eigen-solver (unsymetric case). Now it tries to recover converged eigenvalues even if QZ has failed. Thanks to Mohammad Rahmanian for reporting and helping to reproduce the issue!

December 1, 2015 –

  • Added linsolve function as a simple wrapper over mldivide. Toolbox already has mature solver implemented in mldivide which analyses input matrix properties and applies the best suitable algorithm automatically. No need to re-implement linsolve separately.
  • Improved performance and fixed bug in power functions: .^ and ^.
  • Added functions for saving/loading mp-objects to/from text file:
    A = mp(rand(25));
    mp.write(A,'matrix.txt');  % Write mp-matrix to the text file
    B ='matrix.txt'); % Read it back
    norm(A-B,1) % check accuracy - difference should be 0

    Function mp.write converts mp-matrix to text format and stores it to file. To load the saved matrix back to MATLAB use Data is saved with enough precision to be restored without loss of accuracy.

November 6, 2015 –

  • Special functions – hypergeometric, gamma and whole family of Bessel functions have been updated. In particular, Bessel functions have been rewritten from scratch to support arbitrary orders and arguments, avoid instability regions and accuracy loss in some cases. Now we are not only the fastest in MATLAB world:
    % Symbolic Math Toolbox - R2015a
    >> digits(34);
    >> z = 150*vpa((rand(50)-0.5)+(rand(50)-0.5)*1i);
    >> tic; hypergeom([],vpa(1000-1000*i),z); toc;
    Elapsed time is 6.208366 seconds.
    >> tic; besseli(0,z); toc;
    Elapsed time is 8.635691 seconds.
    >> tic; besselk(0,z); toc;
    Elapsed time is 9.938277 seconds.
    >> tic; bessely(0,z); toc;
    Elapsed time is 17.444061 seconds.
    >> tic; besselj(0,z); toc;
    Elapsed time is 10.570827 seconds.
    % Advanpix Multiprecision Toolbox -
    >> mp.Digits(34);
    >> Z = sym2mp(z);
    >> tic; hypergeom([],mp('1000-1000*i'),Z); toc;
    Elapsed time is 0.155974 seconds.
    >> tic; besseli(0,Z); toc;
    Elapsed time is 2.131102 seconds.
    >> tic; besselk(0,Z); toc;
    Elapsed time is 2.163481 seconds.
    >> tic; bessely(0,Z); toc;
    Elapsed time is 4.301771 seconds.
    >> tic; besselj(0,Z); toc;
    Elapsed time is 3.228625 seconds.

    but also the most accurate:

    % Symbolic Math Toolbox - R2015a
    >> digits(34);
    >> z = 8-43*1i;
    >> besseli(1,vpa(z))
    ans =
    -17.2937_918159785... + 178.7197_1803657...i  % only 6-7 digits are correct!! 
    % Advanpix Multiprecision Toolbox -
    >> mp.Digits(34);
    >> besseli(1,mp(z))
    ans = 
    -17.29378620503294288287975364153844 + 178.7197375141310577798644760813918i % full precision
    Maple:       -17.29378620503294288287975364153844 + 178.7197375141310577798644760813918*I
    Mathematica: -17.29378620503294288287975364153844 + 178.7197375141310577798644760813918*I
    % One more example:
    >> bessely(0,vpa(1000*i))
    ans =
    -2.54099907...376797872e381 + 2.485686096075864174562771484145676e432i  % huge real part?
    >> bessely(0,mp(1000*i))
    ans = 
    -1.28057169...387105221e-436 + 2.485686096075864174562771484145677e+432i % not really

    This is just simplest examples, we have logged many other cases.

    IMPORTANT. Bessel functions from MathWorks Symbolic Math Toolbox do suffer from accuracy loss. Please avoid using it if you need high accuracy.

October 19, 2015 –

  • Added/improved following functions: gradient, linspace, logspace, meshgrid and ndgrid.
  • Fixed bugs in svd(X,0) and in diff.
  • Fixed issue with complex division with infinite components:
    % Before
    >> 1/mp(Inf + 0.1*1i)
    ans = 
        NaN -     NaNi
    % Now (correct)
    >> 1/mp(Inf + 0.1*1i)
    ans = 

Thanks to Ciprian Necula for bug reports!

October 16, 2015 –

  • Thread manager has been improved with better load balancing. Now majority of dense matrix operations are faster by 15-20% (even the matrix multiplication).
  • Fixed two issues reported here and here.
  • Added support for MATLAB’s line spacing setting (format compact/loose):
    >> mp.Digits(34)
    >> mp.FollowMatlabNumericFormat(true);
    >> A = mp(rand(2));
    >> format loose
    >> A
    A = 
            0.4853756487228412241918817926489282        0.1418863386272153359612957501667552    
            0.8002804688888001116708892368478701        0.4217612826262749914363325842714403    
    >> format compact
    >> A
    A = 
            0.4853756487228412241918817926489282        0.1418863386272153359612957501667552    
            0.8002804688888001116708892368478701        0.4217612826262749914363325842714403
    >> format longE
    >> A
    A = 
            4.8537564872284122419188179264892824e-01      1.4188633862721533596129575016675517e-01    
            8.0028046888880011167088923684787005e-01      4.2176128262627499143633258427144028e-01
    >> format shortE
    >> A
    A = 
            4.8538e-01        1.4189e-01    
            8.0028e-01        4.2176e-01

October 8, 2015 –

  • Speed of symmetric eigen-decompositions has been boosted up. All cases are covered: standard and generalized problems in quadruple and arbitrary precision modes and both complexities.Actual speed-up factor depends on number of cores and matrix size. On our 1st-gen Core i7 we see ~3 times ratio for 1Kx1K matrix:
    A = mp(rand(1000)-0.5);
    A = A+A';
    tic; eig(A); toc;
    Elapsed time is 48.154 seconds.
    tic; eig(A); toc;           % x3.4 times faster
    Elapsed time is 14.212 seconds.

    Speed-up is even higher for other precision levels:

    A = mp(rand(1000)-0.5);
    A = A+A';
    tic; eig(A); toc;
    Elapsed time is 204.431 seconds.
    tic; eig(A); toc;           % x5.3 times faster
    Elapsed time is 38.534 seconds.
  • In fact, all operations with dense symmetric/Hermitian and SPD/HPD matrices have became faster: chol, inv, det, solvers, etc. Speed-up ratio for inv is up to 10 times depending on matrix structure and number of cores.

Full comparison table is in preparation.

September 22, 2015 –

  • Estimation of reciprocal of the condition number has been updated for all matrix types in solver (\) and inversion(inv). Now it is much faster and always produce accurate values (before we encountered occasional 0‘s).Just for fun:
    A = magic(5);
    A(3,:) = A(2,:)+10*eps;      % make problem ill-conditioned in double precision
    b = rand(5,1);
    x = A\b;
    Warning: Matrix is close to singular or badly scaled. RCOND =  2.899200e-18. 
    norm(A*x-b,1)                % computed solution is useless (as expected)    
    ans =
    x = mp(A)\mp(b);             % problem is not ill-conditioned when we use 34-digits
    ans = 

    Although half of digits is lost (see the RCOND magnitude), our solver still gives us ~17 digits of accuracy.

September 21, 2015 –

  • Architecture of multi-core parallelism has been revised and improved. You may notice better timings in various dense operations, depending on number of cores and matrix size.For example, this is very beneficent for computations of eigenvectors in unsymmetric case. Which is faster up to 45% on Core i7 (4 hardware cores):
    A = mp(rand(300)-0.5 + 1i*(rand(300)-0.5));
    tic; [V,D]=eig(A); toc;
    Elapsed time is 71.154 seconds.
    tic; [V,D]=eig(A); toc;
    Elapsed time is 38.5 seconds.

    Both modes (quadruple & multiprecision) and complexities (real & complex) are improved.

September 13, 2015 –

  • Following functions have been updated: cond, svd, rank, norm, det, inv, trace, eps and null. Performance optimization, stability, special cases.

September 6, 2015 –

  • Improved matrix analysis stage in dense matrix solver.
  • Improved m-wrapper for conv and colon.
  • Improved compatibility of inv with MATLAB. Now it returns full Inf matrix in case of singular input matrix, estimates RCOND and provides warning if problem is near-singular.
  • Added check for NaN/Inf elements in input matrix in eig.

August 6, 2015 –

  • Added option 'valid' for conv.

August 2, 2015 –

  • Fixed critical issue in subsasgn in multiple-precision mode (when digits <>34).
  • Added workaround for occasional Matlab crashes due to incorrect unload of toolbox (after clear all).

Update is strongly recommended!

July 24, 2015 –

  • Added two-output variants of Cholesky decomposition:

    [R,p] = chol(A)
    [L,p] = chol(A,'lower')
    [R,p] = chol(A,'upper')
  • Added support for 'vector'/'matrix' option in LU decomposition.

July 22, 2015 –

  • Fixed bugs in chol and in auto-detection of numeric constants.
    >> mp.Digits(50);
    >> mp.ExtendConstAccuracy(false);    % auto-detection is disabled by default  
    >> mp(1/3)
    ans = 
    >> mp.ExtendConstAccuracy(true);
    >> mp(1/3)
    ans = 

July 6, 2015 –

  • Added spdiags for sparse matrices.
  • Improved multi-core parallelism in operations with really large matrices (N > 5000). Different thread-scheduling is required (and has been implemented) for the case.

June 4, 2015 –

  • Improved performance of pinv and null with optimized SVD code.
  • Fixed norm to prevent crashes in some cases when input matrix is empty.

May 28, 2015 –

  • Fixes in multi-dimensional array slicing and indexing of complex matrices.
    Thanks to Thomas Rinder, Stefan Güttel and Maxime Pigou for reports and help.
  • Improved performance of arithmetic operations when one of the arguments is complex scalar.

This is maintenance release. Meanwhile we continue focusing on adding sparse matrices functionality in development branch.

April 20, 2015 –

  • Changes in architecture for faster work with sparse matrices.
  • Ordering functions for sparse matrices (to reduce fill-in prior direct solvers) have been added: amd, colamd, symamd and symrcm.Efficient direct solvers are needed for spectral transformation in generalized eigenvalue solver for large matrices.
  • Improved performance of find for sparse matrices.

April 1, 2015 –

  • Minor speed-up in dense linear algebra (especially in SVD) due to more efficient memory management.
  • New analysis functions have been added: issymmetric, ishermitian, isdiag, istriu and istril.
    Syntax and semantic are fully compatible with MATLAB’s built-in functions.
  • Fixed issues with Not-a-Number (NaN) handling in relational operators.

March 19, 2015 –

  • Restrictions on maximum iterations in SVD has been changed to be dependent on precision (needed for the high-precision settings, e.g. 1000 digits or more).
  • Minor fixes in sub-scripted assignment operator, sum and prod.
  • Stability and accuracy of quadgk has been improved. Now it can be used in arbitrary precision mode:
    f = (@(x)sin(x));
    [q,errbnd] = quadgk(f,mp(0),mp(1),'RelTol',100*mp.eps,'AbsTol',100*mp.eps,'MaxIntervalCount',2000)
    q = 
    errbnd = 

March 1, 2015 –

  • System solver (\, mldivide) has been updated with all the recent optimizations.Solver is a meta-algorithm which automatically selects decomposition to use (LU, CHOL or SVD) depending on matrix structure and problem type/size. Now it is optimized for multi-core architectures and shows better performance overall (especially for large problems).
  • Minor updates to arbitrary precision arithmetic engine.

February 23, 2015 –

  • Arbitrary precision Cholesky decomposition, LU and QR have been updated with more optimizations (including multi-core). Speed-up depends on matrix size, larger matrix = higher speed-up.For instance, 1Kx1K shows x3 times better speed whereas 100x100 only 30%.

February 17, 2015 –

  • Fixed memory leak in coefficient-wise operations reported by Ito Kazuho (Thanks!).

February 10, 2015 –

  • Eigenvalue decomposition routines are switched to use “Small Bulge Multi-shift QR Algorithm with Aggressive Early Deflation” (Braman, Byers and Mathias).Speed gain is x2-x3 times for large dense matrices (>500) and decreases with smaller matrix size. The highest speed-up is achieved in multiprecision mode.

January 23, 2015 –

  • Further optimization of large matrix computations by using parallel algorithms on multi-core architectures. Now speed scales better with number of cores:
    >> A = mp(rand(2000));
    % 1 - core CPU:
    >> tic; lu(A); toc;
    Elapsed time is 24.14 seconds.
    % 4 - core CPU:
    >> tic; lu(A); toc;
    Elapsed time is 8.9 seconds.

    Solvers, decompositions and eigen-routines benefit from the update (quadruple & multiprecision mode).

January 15, 2015 –

  • Memory manager for multiprecision objects has been completely re-written with the focus on fast allocation and minimizing memory fragmentation for large matrices. Speed gain depends on matrix size and operation. For example, addition of two 3Kx3K complex matrices shows x3 times improvement:
    % - new version
    >> mp.Digits(50);
    >> A = mp(rand(3000)-0.5 + 1i*(rand(3000)-0.5));
    >> tic; A+A; toc;
    Elapsed time is 3.88 seconds.
    % - previous
    >> mp.Digits(50);
    >> A = mp(rand(3000)-0.5 + 1i*(rand(3000)-0.5));
    >> tic; A+A; toc;
    Elapsed time is 11.9 seconds.

    Speed-up is higher for larger matrices.

  • Several changes to formatted input/output, including new feature of whole-matrix conversion to mp-object:
    >> A = mp('[1/2, sqrt(-1); (-1)^(1/2) pi]')
    A = 
        0.5 +     0i                                          0 +     1i    
          0 +     1i        3.141592653589793238462643383279503 +     0i

    This applies only for matrices with constant elements – useful for scripts which store constant parameters in a matrix.

January 7, 2015 –

  • Scripts with high volume of small scale computations (scalars, small matrices, etc.) are faster by up to 2 times.Usually in such computations the most time was spent in communicating with MATLAB through MEX API functions, which are very slow and impose heavy overhead (e.g. it is not possible to access variable’s memory directly, only after making the deep copy).Today we have found long-awaited workaround for one of such issues – slow creation of ‘mp’-objects in MATLAB. It was affecting every single routine and now it is faster by 50%.
  • Linux version of toolbox has been updated with all the recent changes and improvements.

January 5, 2015 –

  • Multiple precision generalized eigenproblem and SVD functions are 1.5-3 times faster (eig(A,B), qz, svd). Speed-up is higher for larger matrices, e.g. for 2000 x 2000 we can expect 5 times improvement or more.
  • Matrix multiplication is faster by 1.5-2 times – in quadruple and multiprecision mode. Recent versions of MATLAB restrict number of threads allowed to be used in parallel computations. Now we choose this independently from MATLAB. As a result, we can use more cores.
  • Overall, we have been working on multi-core optimizations and you might see noticeable speed-up in large matrix computations.

December 26, 2014 –

  • Multiple precision complex computations are faster by up to 2 times. This affects all the routines including matrix computations – eig, svd, qr, lu, etc.
  • Optimized Kronecker product function kron has been added.

December 24, 2014 –

  • The functions schur, ordschur, qz, ordqz, ordeig, balance and rcond have been extended and optimized. Now we fully support 'real'/'complex' standard and generalized Schur decomposition together with re-ordering in arbitrary and quadruple precision.In a last two weeks the singular/eigenvalues module of the toolbox has been completely changed, 80% of code refactored for improved functionality, stability, speed and feature support. This is part of the work needed for adding EIGS function in next versions.

December 17, 2014 –

  • Eigen-decomposition routines for full matrices have been completely revised. Improved processing of symmetric, Hermitian and symmetric tridiagonal matrices by using MRRR and Divide&Conquer algorithms. Higher speed, better functionality.
  • Routine for generalized eigen-problem, eig(A,B), has been enabled with full multiprecision support without restrictions (before we had it in quadruple precision).

December 12, 2014 –

  • Implemented arbitrary precision divide & conquer SVD algorithm (we had it only for quadruple precision).
    Overall speed-up factor is 7 times for real and complex matrices.

December 3, 2014 –

  • Added element-wise arithmetic operations with scalar for sparse matrices.

November 25, 2014 –

  • Improved compatibility with Parallel Computing Toolbox (removed intra-process blocks).
  • Fixed issue in indexed assignment when LHS is empty matrix.

November 19, 2014 –

  • Added workaround for malfunctioning mxDuplicateArray.
  • Matrix inversion is sped-up in quadruple mode.
  • Small fixes and improvements.

November 14, 2014 –

  • Speed-up of array manipulation operations. All platforms are covered – Windows, Linux and Mac OSX.
  • Fixed bug in matrix power function (for complex matrices).

November 11, 2014 –

  • Optimized memory usage for mp objects after on-the-fly precision change.
  • Fixed minor bugs in colon and matrix power functions.

September 16, 2014 –

  • Performance of matrix computations are boosted up by additional 25-30% (Windows only).

September 9, 2014 –

  • Speed of matrix computations (solvers, decompositions, etc.) is boosted up by 30% for real dense, by 40% for complex dense and by 50% for sparse cases. The improvement is for pure multiple precision computations (not quadruple).

August 28, 2014 –

  • Added indexed assignment capability for sparse matrices (subsasgn).
  • Added nextprime, prevprime, sym2mp, mp2sym, isequal, isequaln, logical and islogical functions.
  • Fixed issues with empty sparse matrices support.

August 21, 2014 –

  • Improved code for generation of nodes and weights for Gaussian-family quadrature.

August 15, 2014 –

  • Added primes, factor, gcd, lcm and isprime functions in arbitrary precision.

August 8, 2014 –

  • Added condeig function.
  • Refined eig to produce matrix of left eigenvectors: [V,D,W] = eig(...).

August 6, 2014 –

  • Added polynomial functions: poly, polyeig, residue, polyder, polyval, polyvalm, deconv, polyfit and polyint.
  • Refined eig to support special flags ('vector'/'matrix', etc.) .

July 31, 2014 –

  • Added fast permute, ipermute, shiftdim, blkdiag, squeeze, circshift and rot90 functions.

July 25, 2014 –

  • Added bsxfun function (and kron as a consequence).
  • Enabled mp-objects to be used as indices in referencing operations.

July 24, 2014 –

  • Further speed-up of data exchange with MATLAB – total speed-up in all operations is 15-20%.

July 21, 2014 –

  • Speed-up of data exchange with MATLAB – all operations are faster now.
  • Fixed bugs related to memory leaks and memory alignment.

May 2, 2014 –

  • Updated hankel and toeplitz to make it compatible with recent changes to toolbox architecture. Thanks to Ilya Tyuryukanov for reporting the bug.

April 14, 2014 –

  • Fixed bug when scalar is passed to diag function.

March 3, 2014 –

  • Speed of sparse matrices computations is boosted up by 3-5 times on Linux platform.

January 22, 2014 –

  • Fixed incompatible exception handling with MATLAB (on Linux).
    Toolbox was catching all exceptions coming from withing the code. However some of the MEX functions throw their own exceptions of hidden (and unknown) type to toolbox. Any attempts to catch them led to MATLAB crash. Now this is fixed – thanks to Takashi Nishikawa for sending the crash dumps.

November 22, 2013 –

  • Improved performance of dense and sparse solvers thanks to re-designed memory layout & access patterns.
  • Fixed triangular solvers, improved matrix analysis routines (positive-definite, etc.) in solvers.

November 8, 2013 –

  • Added precomputeLU function targeted for use in iterative schemes for sparse matrices (e.g. Arnoldi process for computing eigenvalues).New function computes and stores LU factorization of a given sparse matrix directly in toolbox’s core. Then pre-computed LU can be used to solve system of equations with different right-hand side by standard "\". Here is simple example:
    A = rand(1000);
    A = mp(sparse((A>0.5).*A));
    F = precomputeLU(A);
    for k=1:10
     b = mp.rand(1000,1);
     x = F\b;            % re-use of LU with different b

    This approach has advantage over usual x = U\(L\(P*b)) by avoiding overhead of transferring data between MATLAB and toolbox.

  • Most trigonometric & exponential functions are sped up in favor to quadruple precision.
  • Fixed minor bug in sort.

November 1, 2013 –

  • Greatly improved performance of the dense matrix solver (both real & complex) in quadruple precision mode.Now we tear up famous competitors by even greater margin: Advanpix vs. VPA vs. Maple – Dense Solvers and Factorization.This improvement gives us right to claim that now we have the fastest quadruple precision on the planet 🙂
  • Lots of small performance-oriented improvements in toolbox core – avoiding temporaries where possible, better caching, etc.

October 20, 2013 –

  • Added routine for generalized Schur decomposition – qz. Syntax & functionality are equivalent to MATLAB’s routine with the exception that we do not support 'real' flag. Results are computed in 'complex‘ mode (default in MATLAB).
      A = mp(rand(30));                     
      B = mp(rand(30));
      [AA,BB,Q,Z,V,W] = qz(A,B);
      a = diag(AA);
      b = diag(BB);
      l = a./b;
      % Check absolute error
      ans = 
      ans = 
      norm(A*V - B*V*diag(l),1)
      ans = 
      norm(W'*A - diag(l)*W'*B,1)
      ans = 
  • Added new function mp.FollowMatlabNumericFormat(true | false). It allows user to choose numeric formatting preferences for the toolbox. If true, toolbox will obey numeric format settings in MATLAB (show limited number of digits) or display all digits of precision otherwise (by default).Default settings can be adjusted in mpstartup.m script as usual.
    >> mp.Digits(30);
    >> mp.FollowMatlabNumericFormat(false);
    >> mp('pi')
    ans = 
    >> mp.FollowMatlabNumericFormat(true);
    % Fixed-point formats
    >> format short
    >> mp('pi')
    ans = 
    >> format long
    >> mp('pi')
    ans = 
    % Scientific floating-point formats
    >> format shortE
    >> mp('pi')
    ans = 
    >> format longE
    >> mp('pi')
    ans = 
    % Fixed or scientific formats
    >> format shortG
    >> mp('pi')
    ans = 
    >> format longG
    >> mp('pi')
    ans = 
    % C99 hex float format
    >> format hex
    >> mp('pi')
    ans = 

October 16, 2013 –

  • (Preliminary) Added adaptive Gauss-Kronrod numerical integration routine – quadgk. It is completely compatible with default MATLAB function including options and all arguments.
    >> mp.Digits(34);
    >> f = @(x) exp(-x.^2).*log(x).^2;
    >> Q = quadgk(f,mp(0),mp(inf),'RelTol',mp('1e-15'),'AbsTol',mp('1e-30'))
    Q = 
  • Fixed bug in mp constructor to handle the case with recursive precision downgrade, e.g.: mp(mp(pi,50),10). Both cases – dense and sparse are covered.

September 26, 2013 –

  • New sparse matrices serialization, now it is much more efficient and supports NZMAX parameter to reserve space / lower chances for costly re-allocations during computations.
  • Indexed referencing (subsref) is now supported for sparse matrices. For example: A(:), A(2,3), A(1,:), etc. Please note – indexed assignment for sparse matrices is not yet implemented.
  • Fixed bug in diff. MATLAB ignores indexation overloads and calls built-in functions in methods of the class. This is quite an unpleasant surprise (another violation of classic OOP design). Anyway ODE routines now work correctly.

September 17, 2013 –

  • All direct solvers for sparse matrices (LDLT, SuperLU and QR) are optimized for quadruple precision. Speed gain is approx. x10 times. Tables with new timings can be found here.
  • Sparse QR has been fixed and improved, now we can solve undetermined systems as well.
  • Added new functions: conv and poly.
  • Fixed bug in subsref, related to special case of vector indexing.

September 5, 2013 –

  • Added direct solver for sparse matrices, operator "\".Solver includes sparse LDLT, SuperLU and QR decomposition enabled with arbitrary precision support. The most appropriate method is chosen automatically depending on matrix properties.
  • Added arithmetic operations for sparse matrices (+,-,*) and new functions: sparse and transpose. Usage syntax & semantic is completely compatible with MATLAB’s rules.
  • Improved support of empty matrices when used as indices.
  • Speed-up of sparse matrices handling, mixed complexity matrix multiplication, conversion to double, formatted output, etc.
  • Fixed bugs in subsasgn, subsref, norm and roots.

August 13, 2013 –

  • Completely re-designed arbitrary precision arithmetic engine with emphasis on performance. Now real arithmetic operations are up to 20% faster. Complex operations are x2 times faster.

This improvement has major impact on overall performance of the toolbox. Thus, Fourier transform has became x2 times faster (as direct consequence of faster complex arithmetic). Even matrix multiplication receives benefit of 20%-50% increase in performance.

August 9, 2013 –

  • Added functions for sparse matrix manipulations: nonzeros, spfun, spones, full, nnz, nzmax.
  • Improved support of empty matrices and fixed minor bug in numel (case of sparse matrices).

August 6, 2013 –

  • We have re-fined common array operations with improved multi-dimensional support: sum, prod, cumsum, cumprod, max, min, sort, fft, ifft.
  • Fourier transform speed-up by 25%.

July 31, 2013 –

The first beta version with rudimentary sparse matrices support.  

  • Very basic support of sparse matrices in multiple precision. Will extend it in upcoming versions. See Rudimentary Support for Sparse Matrices for more details.
  • Custom implementation of indexing, referencing and similar operations (subsref, subsasgn, cat, etc.).
  • Matrix computations x2-3 times speed up on Windows 64-bit.
  • Display & formatted output improvement

July 12, 2013 –

  • We have added mp.randn function for generating normally distributed pseudorandom numbers with arbitrary accuracy. All special cases are supported for full compatibility with MATLAB:
    r = mp.randn(n)
    r = mp.randn(m,n)
    r = mp.randn([m,n])
    r = mp.randn(m,n,p,...)
    r = mp.randn([m,n,p,...])
    r = mp.randn
    r = mp.randn(size(A))
    r = mp.randn(..., 'double')   %  'double' is ignored
    r = mp.randn(..., 'single')   %  'single' is ignored

    Also we have refined mp.rand, now it is faster and able to generate multidimensional arrays as well.

July 3, 2013 –

  • To resolve the problem with mixed usage of limited accuracy double precision constants in expressions with mp entities, we have introduced new global setting in the toolbox: mp.ExtendConstAccuracy().

    It enables/disables auto-detection and recomputing of the constants with higher precision to match toolbox’s settings, e.g.:

    >> mp.Digits(34);
    >> mp.ExtendConstAccuracy(false);
    >> sin(mp(pi))
    >> mp.ExtendConstAccuracy(true);
    >> sin(mp(pi))

    In the first example toolbox uses pi as it is, with double precision accuracy (at most 16 correct digits). In the second example, toolbox recognizes the pi constant and re-computes it with required precision of 34 digits.

    Be default (can be changed in mpstartup.m) we use mp.ExtendConstAccuracy(true).

June 26, 2013 –

  • After few reports from confused users we have removed auto-detection and recomputing of commonly used constants (pi, eps, etc). Now all double precision constants are converted to mp as it is – with at most 16 correct digits. Before toolbox was trying to recompute them in higher precision.

May 20, 2013 –

  • Improved support for empty arrays/matrices.
  • Optimized and improved rcond (including quadruple precision).
  • Speed up of operations with multidimensional arrays.
  • Fixed minor bugs of quadruple precision mode.

April 21, 2013 –

  • Fixed bugs in incomplete gamma function computation.
  • Speed up of determinant computation.

April 15, 2013 –

  • Real & complex Schur decomposition has been improved with more intelligent restriction on allowable number of Francis QR iterations. Thanks to Gokhan Tekeli from Bogazici University!

March 21, 2013 –

  • Fixed bug in expm. Thanks to Dmitry Smelov from Stanford!
  • Fixed bug and improved performance of colon.

March 8, 2013 –

New feature:  

  • Added support for multidimensional arrays and operations with them. We support coefficient-wise arithmetic operators, mathematical functions, basic information, array manipulation and other functions. Example of array manipulation:
    >> x = mp(eye(2));
    >> a = cat(3,x,2*x,3*x)
    (:,:,1) = 
    1    0    
    0    1    
    (:,:,2) = 
    2    0    
    0    2    
    (:,:,3) = 
    3    0    
    0    3    
    >> B = permute(a,[3 2 1]);
    >> C = ipermute(B,[3 2 1]);
    >> isequal(a,C)
    ans =

February 13, 2013 –


  • 20% performance increase in multiprecision linear algebra.
  • Fixed bug in Cholesky decomposition in quadruple precision mode.

February 6, 2013 –

New features:

  • We have re-implemented QR decomposition with optimizations for quadruple precision. This resulted in significant speed-up by 10-20 times. We support following variants of qr function:X = qr(A)
    X = qr(A,0)
    [Q,R] = qr(A)
    [Q,R] = qr(A,0)
    [Q,R,E] = qr(A)
    [Q,R,E] = qr(A,0)

    Few tests with timing:

    >> mp.Digits(34);                 % Quadruple precision
    >> mp.GuardDigits(0);
    >> A = rand(100,50);              % 100 x 50 test matrix
    >> X = mp(A);
    >> tic; [Q,R] = qr(X); toc;       % Full QR
    Elapsed time is 0.152928 seconds.
    >> norm(X-Q*R,1)                  
    >> tic; [Q,R] = qr(X,0); toc;     % Economic QR
    Elapsed time is 0.110708 seconds.
    >> norm(X-Q*R,1)                  

January 20, 2013 –

New features:

  • This version introduces very important feature we have been working for a few months. In order to boost performance we have made thorough speed optimization of our core engine for computations in quadruple precision (34 decimal digits).As a result we have 20-50 times better performance for matrix computations done in quadruple precision. Here is example of computing singular values of a 100×100 matrix:
    >> mp.Digits(34);                 % Switch to quadruple mode by using 34 decimal digits
    >> mp.GuardDigits(0);
    >> A = mp(rand(100));             % Use 100 x 100 matrix
    >> tic; svd(A); toc;
    Elapsed time is 0.133788 seconds.
    >> tic; [U,S,V] = svd(A); toc;
    Elapsed time is 0.668525 seconds.
    >> norm(A-U*S*V',1)               % Accuracy check 

    Basic matrix operations, eigensolvers and some decomposition routines are already updated to be much faster in quadruple precision. Further extensions are under development.

  • Routine eig() is extended to support arbitrary matrices A,B when computing solution of generalized eigenvalue problem with quadruple precision.
    >> mp.Digits(34);                 
    >> mp.GuardDigits(0);
    >> A = mp(rand(100));
    >> B = mp(rand(100));
    >> [V,D] = eig(A,B);
    >> norm(A*V - B*V*D,1)

January 9, 2013 –

New features:

  • Added new solver for ordinary differential equations – ode113.
  • Both ODE solvers ode45 and ode113 are enabled with the support of event functions.

Bug fix:

  • Fixed minor bug in operations involving mp-objects with different precisions.

December 28, 2012 –


  • Fixed incompatibility with the latest OpenMP version included in MATLAB R2012b. Our sincere gratitude goes to Kazuho Ito from University of Yamanashi for his excellent help in finding and investigating the problem.Unfortunately it comes with the price of dropping support of old R2008b. Now toolbox supports all versions of MATLAB starting from R2009b.

December 19, 2012 –


  • Performance has been improved by 20% – 50% in most of mathematical operations.
    We have switched to Intel C++/Fortran Compilers and their new OpenMP runtime provides this gain.  

    Side effect is that Windows & MATLAB R2008b users need to define special environment variable KMP_DUPLICATE_LIB_OK=TRUE in order to enable compatibility with older Intel Compilers (used to build R2008b).

December 1, 2012 –


  • Due to growing number of customers using Windows 8 we have added support for the OS.

November 7, 2012 –


  • Few months ago we have started re-implementation of linear algebra algorithms to benefit from multi-core parallelism. Basically this means that performance of the toolbox will be better on systems with more cores/CPUs.Today’s release is the first version with enabled parallelism in basic matrix operations, like multiplication.

Bug fix:

  • Fixed (similar)bugs in matrix left and right divide. We have used in-proper speed optimization for a special case when one matrix is real and other has complex elements.

October 31, 2012 –


  • Performance is increased by 2.5-3 times in majority of linear algebra functions including: qr, svd, eig, schur, hess.
  • Algorithm for automatic detection and re-calculation of common constants is improved to avoid false-positive errors.

October 24, 2012 –

Bug fix:

  • Fixed bug in matrix properties analysis stage of eig() function. Imaginary part of elements were erroneously stripped off in case of complex diagonal matrices.


  • Speed up of multi-precision arithmetic engine by 5-10%.
  • Dynamic memory manager is tuned to gain more speed on Windows 7.

October 3, 2012 –

New features:  

  • Added new function, balance aimed to improve accuracy of computation of eigenvalues and/or eigenvectors. It applies similarity transformation (permutation & scaling) to a matrix so that row and column norms becomes approximately equal.Algorithm is based on xGEBAL, xGEBAK from LAPACK library. All MATLAB’s features are supported:[T,B] = balance(A)
    [S,P,B] = balance(A)
    B = balance(A)
    B = balance(A,'noperm')

October 2, 2012 –

New features:  

  • Added new function, rcond to compute the 1-norm estimate of the reciprocal condition number.
    Algorithm is based on xGECON routines from LAPACK, both complex and real matrices are supported. Usage syntax is compatible with MATLAB:  

    c = rcond(A)

September 25, 2012 –

New features:  

  • Added new functions, mod and idivide. They are needed for some built-in functions to work correctly with multi-precision entities (e.g. unwrap).


  • Mathematical expression parsing & evaluation have been completely re-written to be more error-robust and flexible.

August 8, 2012 –

New feature:  

  • Now display format of mp-entities are controlled by MATLAB’s formatted output settings. We support the following formats: short, long, shortE, longE, shortG, longG, shortEng, longEng and hex.Short formats show only 4 digits after the decimal point. Long formats show full precision.
    >> mp.Digits(30);
    % Fixed-point formats
    >> format short
    >> mp('pi')
    >> format long
    >> mp('pi')
    % Scientific floating-point formats
    >> format shortE
    >> mp('pi')
    >> format longE
    >> mp('pi')
    % Fixed or scientific formats
    >> format shortG
    >> mp('pi')
    >> format longG
    >> mp('pi')
    % C99 hex float format
    >> format hex
    >> mp('pi')

July 26, 2012 –

New features:  

  • Added functions for conversion to integers and vice versa: int8, uint8, int16, uint16, int32, uint32, int64 and uint64:
    >> x = mp(int64(magic(3)))
    8    1    6    
    3    5    7    
    4    9    2    
    >> int64(x)
    ans =
                        8                    1                    6
                        3                    5                    7
                        4                    9                    2
    >> x = mp(intmax('uint64'))
    >> uint64(x)
    ans =

    An interesting fact is that MATLAB rounds floating point numbers to nearest integer in conversion:

    >> int32(1.2)
    ans =
    >> int32(1.5)
    ans =

    This is quite unexpected and contradicts the majority of other programming languages, where decimal parts are just truncated.

July 25, 2012 –

New features:  

  • Added functions for specialized matrices computing: compan, hankel, vander, toeplitz, mp.hilb and mp.invhilb.Hilbert matrix inversion test:
    % Multiprecision Computing Toolbox:
    >> mp.Digits(100);
    >> norm(mp.invhilb(20) - inv(mp.hilb(20)),1)
    % MATLAB:
    >> norm(invhilb(20) - inv(hilb(20)),1)
    Warning: Matrix is close to singular or badly scaled.
             Results may be inaccurate. RCOND = 1.155429e-019. 
    ans =

    Please note, integer-valued special matrices can be converted to mp objects directly, e.g.:


Bug fixes:

  • Fixed bug in element-wise power function when NaN was returned for small negative arguments.

July 19, 2012 –


  • Code for data interchange between toolbox and MATLAB has been re-written to be more generic. This is the next step towards sparse matrices support. This part is greatly optimized thanks to reduced number of heap memory (de-)allocations.

July 10, 2012 –


  • We have implemented adaptive computational load balancing between Pade approximation and ISS (inverse scaling and squaring) in matrix logarithm. Now logm() is more stable and much faster.

July 6, 2012 –

New features:  

  • Added functions for 2-D fast Fourier transform, fft2 and ifft2. All features and special cases are supported to provide full compatibility with standard functions:

    Y = fft2(X)
    Y = fft2(X,m,n)

    Y = ifft2(X)
    Y = ifft2(X,m,n)
    y = ifft2(..., 'symmetric')
    y = ifft2(..., 'nonsymmetric')

  • Added subsindex function for smooth usage of mp objects as indexes to matrix coefficients.

Bug fixes:

  • Fixed bug in ifft related to special case when a 'symmetric' option is supplied.

June 22, 2012 –

New features:  

  • Added fprintf() function. Now mp numbers can be printed the same way as standard floating point numbers:
    >> fprintf('double pi = %f and \n50-digits pi = %.50f\n', pi, mp(pi,50))
    double pi = 3.141593 and 
    50-digits pi = 3.14159265358979323846264338327950288419716939937511

    This combined with auto-recomputing of commonly used constants makes existing scripts porting to arbitrary precision even easier. There are much less modifications needed in code than before.


  • We have implemented new heap memory management for entities of the mp type. Memory of “disposed” objects is not freed immediately but can be re-used for new objects without slow system calls (allocation/deallocation). This gives up to two times a performance boost in all linear algebra functions.

June 20, 2012 –

New features:  

  • Automatic detection and re-calculation of commonly used double constants if they are encountered in expressions with multi-precision numbers. Usage of limited precision double constants (like pi, exp(1), sqrt(2), etc.) in arbitrary precision computations has always been one of the main source of low accuracy final results.General rule is that floating-point constants should be re-computed in high precision from the beginning:

    See More on Existing Code Porting for details.

    The new version of toolbox automatically detects and re-calculates common constants encountered in computations with arbitrary precision, including pi, sqrt(2), exp(1), log(2), eps:

    >> pi
    ans =
    >> mp(pi,50) % pi is automatically re-computed to have 50 correct digits
    >> mp(eps,50)  % eps is automatically recognized and adjusted to supplied precision
    >> mp(10*sqrt(2)/571,50) % fractions with constants are also supported

    Additionally toolbox correctly recognizes constants if they are used with fractional coefficients, e.g: 7*pi/3, 10*eps, sqrt(2)/2.

    Hopefully this new feature will make porting of existing programs to arbitrary precision even easier.

Bug fixes:

  • Extended dot product to support vectors of different shapes, e.g. row – column.
  • Fixed incompatibility of funm with older versions of MATLAB.

June 14, 2012 –

New features:  

Matrix functions (logm in particular) have been completely revised thanks to feedback of Numerical Linear Algebra Group from The University of Manchester. Now we use the state-of-the-art Schur-Parlett algorithm for computing general matrix function described in the following references:

  • P. I. Davies and N. J. Higham, A Schur-Parlett algorithm for computing matrix functions. SIAM J. Matrix Anal. Appl., 25(2):464-485, 2003.
  • N. J. Higham, Functions of Matrices: Theory and Computation, Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2008.

June 12, 2012 –

New features:  

  • Added ordschur() – eigenvalues reordering in complex Schur factorization. All special cases are supported:

    [US,TS] = ordschur(U,T,select)
    [US,TS] = ordschur(U,T,keyword)
    [US,TS] = ordschur(U,T,clusters)

    In some situations, order of eigenvalues within clusters might not coincide with what is generated by MATLAB. This is because we do additional minimization of permutations to gain more speed. Otherwise functionality is identical to MATLAB’s built-in ordschur().
  • Added solver for triangular Sylvester equation, A*X + X*B = C. Syntax is trisylv(A,B,C).

June 5, 2012 –

New features:  

  • Added matrix power function, mpower(). We use binary squaring for integer powers and combination of expm, logm for other cases (preliminary).
  • Added support for logical variables, now can be used with mp numbers in comparisons, conversions, and expressions:

    mp('pi') + true

Bug fixes:

  • Fixed premature stop in Schur factorization due to restriction on maximum iteration count in Francis QR step algorithm.
  • Fixed program crash when a user without administrator rights tries to install toolbox in restricted folders (e.g. Program Files)


  • Implemented workaround for libstdc++ & dynamic std::string problem on Mac OS X. No more segfaults because of this!
  • Improved performance in mp output routines thanks to removed dependency on boost::format which was extremely sloooow and buggy on Mac OS X.

Above ‘Version History’ doesn’t reflect development prior June 2012.

Idea for Multiprecision Computing Toolbox was born on September 13, 2010. Development started in February 2011 (one of the first commits was done during the Great East Japan Earthquake on March 11, 2011 in shaken building midst falling shelves, computers, books and panicked people). First public version was released on September 16, 2011. Major new versions are released twice a month, with minor updates – almost every day.