I have checked & tested, such vectorization won't be faster than loop, unless input arrays are really long (5M-10M). We can do the long arrays faster by parallelism. But short arrays won't benefit at all.

Please use loop or arrayfun:

B=arrayfun(@(n,k)nchoosek(n,k),mp(2:4),mp(1:3))

Could you please elaborate a bit more how this should work and what output should be generated? 

This case is ambiguous and can be treated differently: one to one, one to all. Also this will make the function incompatible with MTLAB's built-in.

This is not an option either. Putting the slowness aside (applying the seed/state is a non-trivial process), there is another reason:

Random state must be reset to default on MATLAB startup (when user starts up the MATLAB/starts work with toolbox).

MATLAB doesn't allow differentiating startup vs. 'clear all' events. These events must be treated differently, but we can choose only one. We choose the most flexible: "startup" behavior and allow the user to take care of random state as one wishes.  User just need to call mp.RandState when needed (as it was designed).

This doesn't lead to any slowdowns, undocumented hacks/hooks to MATLAB kernel, etc. etc. etc.

Update: I edited the comment after posting it (to remove logically incorrect statement).

@"Matlab’s built-in random number generator never treats “clear all” as a command to reset or reseed its generator."

MATLAB doesn't notify the toolboxes (nor MEX modules) when user calls 'clear' function on it. So that toolbox has no way to react to this event in some special way (e.g. to store/pass the seed to the next call). The 'clear' just kills the toolbox from executing, unloads and clears its memory. Next call to toolbox starts it fresh from default state.

Please send MT request to allow flexible handling of 'clear' function in 3rd party toolboxes/MEX.

Now toolbox provides two ways of using random number generators:

(1) Using the standard MATLAB functions, e.g. mp(rand(..)). The seed functionality is controlled by MATLAB. But precision is limited to double precision numbers.

(2) The generation of a full-precision random numbers requires special syntax (rand(..,'mp')) and a separate seed control  function (mp.RandState). The separate seed control function is needed because of the MATLAB's restrictions on overloading of 'rng' and 'clear' functions (it is not possible to do or there is no notification mechanism provided to handle the events).

User can choose the (1) or the (2). In case of (2) manual control of a seed is required using the special function mp.RandState. 

The 'clear all' wields full power over the MEX - as it just kills all the MEX modules (including mp-toolbox).
The mp.RandState is provided exactly for such cases when precise and explicit control over generator is needed. Please use it.

Jon, I have just released new version of toolbox - 5.3.6.15927.

It doesn't use the MATLAB's implementation of nchoosek anymore.

Instead it is now implemented directly in toolbox core.

>> nchoosek(mp(76), mp(32))
ans = 
    2695592391875730827550

Please update your environment.

I have just checked the source code of the nchoosek in the newest versions of MATLAB.
MathWorks changed it to do all the computations in double. Here is some code from the function: 

% Do the computation in doubles.  
nd = full(double(n)); 
kd = full(double(k));           

I am not sure, but the code might start work fine after replacing all the "double" to "mp".

Looks like MW has no plans to support OOP/custom classes properly at all.

Meanwhile I will try to implement the function in toolbox's core.

Jon, Happy New Year!

We have just released new version 5.3.6.15902: https://bit.ly/3mvxXqg
Now, its least-squares solver generates bitwise deterministic results every time.

Would appreciate it if you would update me on your test results.

Repeated runs of the same problem in the same environment will yield the same bitwise result unless something in the environment changes; e.g, CPU type, number of threads,...

I think with more advanced algorithms and many-core CPUs this will change in the future.

For example. Even though we use the same number of threads in a solver - it might split the overall work in different number  of jobs and distribute them differently among cores on each call. This might happen because some CPU cores might be busy with other processes on the computer. More flexible/adaptive algorithms take this into account (as we did in least-square solver before 5.3.6).

Thank you for all your support and contributions!

Best wishes for the New Year!

Pavel

Jon, thank you very much for the report!

Generally, some degree of variability of the solution is to be expected (due to matrix conditioning, multi-threading, flow of divide-and-conquer/bulge chasing computations, etc.).

The condition number of a matrix is log10(cond(Amat)) = ~21, roughly meaning that last ~21 digits are expected to be lost during computation of the solution.

However, even if we use extended precision, the last ~21 digits of the computed solution are not accurate (=can vary every time we solve the system). What counts is the relative accuracy of the system delivered by the solution. 

I added few changes to your script:

clear all

ndigits=34;
mp.Digits(ndigits);

load vL
vL = mp(vL,ndigits);
Amat=mp(Amat,ndigits);

%vL = double(vL);
%Amat = double(Amat);

n0 = norm(Amat,1)*norm(vL, 1); % PH: norms to compute the relative accuracy

xi = Amat\vL;  %backslash error is nonzero when differenced with same backslash later on

for i=1:10
    notabug = Amat\vL;
    
    acc = norm(Amat*notabug-vL,1) / n0;     % PH: system rel. accuracy w.r.t solution vector.
    xacc = norm(xi - notabug,1)/norm(xi,1); % PH: rel. difference from the first solution.
    [m,idx] = max(abs(xi - notabug));
    
    fprintf("rel.acc(sys) = %.5e, rel.acc(x) = %.5e, [xi(%d) = %.8f, notabug(%d) = %.8f, rel.acc = %.5e]\n",acc,xacc,idx,xi(idx),idx,notabug(idx),abs(xi(idx)-notabug(idx))/abs(xi(idx)));
    
    %sum(abs(xi - notabug))
endv

Here is its results for precisions: 34, 50, 100 and 500 decimal digits.

% ndigits = 34
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 2.54932e-19, [xi(40) = 33835.88114407, notabug(40) = 33835.88114407, rel.acc = 3.10975e-19]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 1.24177e-10, [xi(41) = 5382.66065303, notabug(41) = 5382.66064798, rel.acc = 9.37861e-10]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 8.72727e-11, [xi(41) = 5382.66065303, notabug(41) = 5382.66065016, rel.acc = 5.33988e-10]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 1.04572e-10, [xi(41) = 5382.66065303, notabug(41) = 5382.66064912, rel.acc = 7.26616e-10]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 8.19852e-11, [xi(41) = 5382.66065303, notabug(41) = 5382.66064979, rel.acc = 6.01897e-10]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 1.23477e-10, [xi(41) = 5382.66065303, notabug(41) = 5382.66064801, rel.acc = 9.32217e-10]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 1.23477e-10, [xi(41) = 5382.66065303, notabug(41) = 5382.66064801, rel.acc = 9.32217e-10]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 1.24177e-10, [xi(41) = 5382.66065303, notabug(41) = 5382.66064798, rel.acc = 9.37861e-10]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 7.31347e-11, [xi(41) = 5382.66065303, notabug(41) = 5382.66065016, rel.acc = 5.33995e-10]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 1.24177e-10, [xi(41) = 5382.66065303, notabug(41) = 5382.66064798, rel.acc = 9.37861e-10]

% ndigits = 50
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 3.92138e-30, [xi(40) = 33835.88114392, notabug(40) = 33835.88114392, rel.acc = 3.26593e-30]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 9.09178e-27, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 6.17458e-26]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 5.67413e-27, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 3.50551e-26]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 5.65223e-27, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 3.38773e-26]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 5.21849e-27, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 3.40105e-26]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 3.92138e-30, [xi(40) = 33835.88114392, notabug(40) = 33835.88114392, rel.acc = 3.26593e-30]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 6.09511e-27, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 4.11961e-26]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 5.21855e-27, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 3.40105e-26]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 6.12574e-27, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 4.11961e-26]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 4.69512e-27, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 2.96222e-26]

% ndigits = 100
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 5.74485e-81, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 4.02796e-80]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 1.42945e-79, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 9.44583e-79]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 5.76549e-81, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 3.97130e-80]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 1.42444e-77, [xi(40) = 33835.88114392, notabug(40) = 33835.88114392, rel.acc = 1.04581e-77]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 1.42410e-77, [xi(40) = 33835.88114392, notabug(40) = 33835.88114392, rel.acc = 1.14507e-77]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 1.57315e-79, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 1.02690e-78]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 1.42394e-77, [xi(40) = 33835.88114392, notabug(40) = 33835.88114392, rel.acc = 1.14487e-77]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 1.99963e-77, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 1.18875e-76]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 1.99972e-77, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 1.18875e-76]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 1.47611e-77, [xi(40) = 33835.88114392, notabug(40) = 33835.88114392, rel.acc = 1.15290e-77]

% ndigits = 500
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 2.88681e-477, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 1.15224e-476]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 3.48807e-477, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 2.16831e-476]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 2.88681e-477, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 1.15224e-476]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 0.00000e+00, [xi(1) = 0.00000000, notabug(1) = 0.00000000, rel.acc = 0.00000e+00]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 0.00000e+00, [xi(1) = 0.00000000, notabug(1) = 0.00000000, rel.acc = 0.00000e+00]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 2.88681e-477, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 1.15224e-476]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 2.88681e-477, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 1.15224e-476]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 0.00000e+00, [xi(1) = 0.00000000, notabug(1) = 0.00000000, rel.acc = 0.00000e+00]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 0.00000e+00, [xi(1) = 0.00000000, notabug(1) = 0.00000000, rel.acc = 0.00000e+00]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 2.88681e-477, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 1.15224e-476]

First column is the system accuracy w.r.t current solution vector, second - current solution vector difference from the first one. Lastly in square brackets we have the current solution vector element with max difference from the first vector and its relative difference.

The relative accuracy of all the solutions are the same (1.55559e-15), this is clearly the maximum achievable accuracy for the system. 

The solution vectors differ from the first one by last ~21 decimal digits. This coincides with the matrix conditioning (as we expected). The solution vector elements with max. difference also differs by the same (very small) margin from the first vector.

Overall, this perfectly matches the expected behavior of matrix solution.

***

Could you please send me references to MATLAB's requirements for solutions of overdetermined system to be bit-by-bit deterministic? In this particular case the MATLAB doesn't use any advanced QR algorithms/multithreading and probably this leads to exactly the same solution each time. But, mathematically it is not guaranteed or expected/justified. 

In toolbox we use advanced algorithms for the case (multi-threadfinfg + D&C + bulge chasing with different execution paths each time). Nonetheless the solutions provide the same relative accuracy and differ from each other by expected margin (which is log10(cond(A)) = ~21 digits). This is how it is expected to behave. Let me know your comments.

I have just released new toolbox version (5.3.5.15874) with new function mp.RandState. It allows precise control of the random generator state. Its usage syntax and semantic is the same as for built-in 'rng', e.g.: 

mp.RandState(seed)
mp.RandState('default')
mp.RandState('shuffle')
mp.RandState(seed,'twister')
mp.RandState(state)

The mp.RandState returns full state of the generator, therefore it is exactly reproducible.

Simple example can be found here: https://www.advanpix.com/documentation/version-history/