mp.Digits(N) returns current precision (before changing it to N) in order to be used in a constructs like this:

prevN=mp.Digits(newN)
% Compute with precision "newN"
...
mp.Digits(prevN) % restore the precision 

Our documentation explains it as following:

"Setup the default precision. All subsequent multiprecision calculations will be conducted using the newly specified number of decimal digits. The function returns the previously used default precision. If called without argument, the function simply returns the current default precision in use."

Let me know if it is too vague or not accurate - will try to re-word it in some way. 

Thank you very much for the bug report and for the details on your work, very interesting!

Let me know if you need my help with anything else. I am marking the bug as fixed.

I have fixed the issue (see email). It was related to the same mix of indices/values in spdiags code (in different place).

The new spdiags will be included in the next toolbox update. Please use the file I sent you until then.

Overall, it is pretty common for M-language to mix the integers and floating point values in the same array.

Please let me know if you will find similar situations. 

The low precision of 1 decimal digit is quite an extreme case, you are the first user we know who uses it.

Could you share what problem are you solving?

Thank you for the report!

It is an interesting and tricky bug. Usual spdiags uses double arrays to store both - the indices and the numeric values.

If we convert the code to support mp-type, then mp array is used to store the indices and numeric values too.

Everything is Ok as long as precision is capable of storing the indices accurately. 

In the case of mp.Digits(1) this is no longer valid - and indices becomes inaccurate after hitting two digits.

I have sent you the fixed mpspdiags code via email.

Simple switching to using hardware "double" instructions is impossible. MATLAB, Intel, etc. have been working on implementing "double" precision math functions for decades. To make it efficient, fast, with small memory footprint, to utilize CPU cores, etc. etc. This is enormous work, very difficult to replicate.

We are focusing on doing the same for arbitrary precision computations - to make them as fast as possible on modern CPU. This requires a lot of efforts from deriving new algorithms capable of doing extended precision to low-level software optimizations. The algorithms and code are inherently created for arbitrary precision. No simple switching to "double" is possible. We have been working on this for more than 10 years with long todo list for the next 10 years :)

No special hardware instructions exist for quadruple precision. It is all emulated in software using 64-bit integer CPU instructions. We created it exactly because quadruple is not implemented in hardware.

Double precision is already implemented in hardware, and its functionality has been polished very well. So that, for adequate comparison it is better to use existing "double" precision libraries, provided in MATLAB, MKL, etc.

It is possible to write precision-independent code, so that it can run with "double" or "mp" scalars. Please check the page for more details: https://www.advanpix.com/2016/07/21/how-to-write-precision-independent-code-in-matlab/

The mp.Digits(16) is not equivalent to "double" since we are using wider exponent range and our arithmetic/basic math. functions are done with guaranteed accuracy. In other words, mp.Digits(16) is delivering more accurate results compared to native "double".

I agree, the extended precision must be applied carefully, to the parts where ill-conditioning is observed (cancellation, accumulating rounding errors, etc.). Or to verify the results and study asymptotic properties of the algorithms. 

Hello Michal,

Good timing! We have just finished all the features we planned for major release, please download 4.9.2.

It includes all the optimization work for (full) matrix computations we have been working on lately.

In this regard the 4.9.1 was transitional with partial results integrated (basic matrix operations only).

Please use the 4.9.2.

Jon,

Thank you very much for the excellent question.

The main overhead in such loops are indexed access to the elements of array.
MATLAB makes this operation slow for third-party toolboxes. In short, it requires the toolbox to make a copy of the whole array, even if only one element is requested (MEX API restrictions).

That is why we always suggest vectorizing code as much as possible (avoid loops and process data as a whole).
However application of IIR filter cannot be easily vectorized. We have been planning to add it to the toolbox core with all possible optimizations for a long time. Your request motivated us to prioritize this. Thanks!

Please download a new version of the toolbox (4.9.0.14753) from our website. Among other changes it includes the function 'filter' which implements IIR filter.
Thus your example becomes:

tic
r1 = filter([0], [1 -cosom], y1, [y1(1)]);
r2 = filter([0], [1 -cosom], y2, [y2(1)]);
toc

However, if you need to apply same filter on several arrays simultaneously, we would suggest to process them altogether:

y = [y1 y2];
tic
r = filter([0], [1 -cosom], y, [y1(1) y2(1)]);
toc

This will allow the toolbox process each array (column) in parallel, using different CPU cores.

@"The fft() respects input precision of the arguments."
All functions provided in toolbox respect precision of the input arguments.

However, user-written functions might behave differently - this is responsibility of the user who writes the code. See my previous reply for the cases when precision might be different inside user-written function (e.g. when mix of different precisions is used).

Could you please provide minimal, complete, and verifiable example of the improper behavior of the toolbox?

This would help a lot.

Edit:

@"This leads to an expectation of consistency; i.e, a subroutine (user-written or otherwise) should always be executed in precision of the arguments."

All functions in toolbox provide this consistency. But toolbox has no control of user-written code. All user's M-code is executed by MATLAB interpreter by following the object programming rules (it calls functions from toolbox when it encounters mp-variables). We cannot change this and force the MATLAB interpreter to behave differently.

In order to have this consistency in user-written functions, please use the same working precision everywhere. This is the only simple rule.