Vos commentaires

Dear Mohsen,


You are absolutely right - tests show that interp1 is unacceptably slow!

The slowest code is index lookup of interval where x belongs while interpolation.

Now it is simple linear search and of course should be improved.


Basically we have to implement histcounts, which is not trivial.


Will try to improve this in near days and let you know about the result.

Updated version has been released, download link: http://goo.gl/pMXV3


Now I am working on speed optimization of Bessel functions for special cases.

For example new version already has much faster bessely(n,x) when n - integer, x - non-negative real in quadruple precision.


Please try it (if you use it in your work). I am finishing besselj - will be available very soon as well.

Bug has been fixed, updated version will be released in a few hours (after compilation finishes).

All relational operations were affected, and also element-wise power x.^y.

Dear Michael,


Thank you for your report, will fix it asap.

(Recently I have added support for sparse matrices in relational operations - probably something went wrong).

Just want to let you know, that I have extended toolbox with support of "shallow" complex numbers:


>> x = mp(complex(12))
x = 
    12 +     0i

>> isreal(x)
ans =
     0

But, in comparison to MATLAB, toolbox supports complex numbers much more correctly, including branch cuts, signed zeros as imaginary part, etc.


If you are interested in details, I wrote short post with examples for demonstration:

http://www.advanpix.com/2016/04/28/branch-cuts-and-signed-zeros-in-matlab/

Toolbox allows working with floating-point numbers from 9.5303e-323228497 to 5.2464e+323228495 with arbitrary precision (arbitrary length of mantissa):


>> format longG
>> mp.Digits(34);
>> mp('pi')
ans = 
    3.141592653589793238462643383279503

>> mp.Digits(100);
>> mp('pi')
ans = 
    3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862
8034825342117068

>> mp.Digits(1000);
>> mp('pi')
ans = 
    3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862
80348253421170679821480865132823066470938446095505822317253594081284811174502841027019385
21105559644622948954930381964428810975665933446128475648233786783165271201909145648566923
46034861045432664821339360726024914127372458700660631558817488152092096282925409171536436
78925903600113305305488204665213841469519415116094330572703657595919530921861173819326117
93105118548074462379962749567351885752724891227938183011949129833673362440656643086021394
94639522473719070217986094370277053921717629317675238467481846766940513200056812714526356
08277857713427577896091736371787214684409012249534301465495853710507922796892589235420199
56112129021960864034418159813629774771309960518707211349999998372978049951059731732816096
31859502445945534690830264252230825334468503526193118817101000313783875288658753320838142
06171776691473035982534904287554687311595628638823537875937519577818577805321712268066130
01927876611195909216420199

The range 9.5303e-323228497 to 5.2464e+323228495 is fixed but mantissa(precision) can be of any length.

The range covers pretty much all atoms in universe.


Could you show example where you need wider exponent range?

Yes, max/min numbers are the same for all precisions except IEEE quadruple.


Range of numbers depends on floating point exponent.

It is small for IEEE quadruple. For all other precisions we use the same large exponent.

IEEE 2008-754 quadruple 34-digits precision is not enough to hold the 1e+10000 numbers. Please use higher precision.


You can always check the max/min number for current precision using functions: mp.realmax, mp.realmin:


>> mp.Digits(34);
>> mp.realmax, mp.realmin
ans = 
    1.1897e+4932
ans = 
    3.3621e-4932

>> mp.Digits(100);
>> mp.realmax, mp.realmin
ans = 
    5.2464e+323228495
ans = 
    9.5303e-323228497

The code mp(1e+500) is incorrect. The 1e+500 is converted to double precision first (truncated to 1e+308).

Use mp('1e+500') instead.


mp.Digits(1000000) is too high. The usual mp.Digits(34) is enough to handle such small number as 1e+500 :).


Also please check the following chapter in our User Guide, it has examples on the issue:

http://www.advanpix.com/documentation/users-manual/#More_on_Existing_Code_Porting

Of course, toolbox supports large/small numbers:


>> mp.Digits(34);

>> besselk(0,1000)  % double
ans =
     0

>> besselk(0,mp(1000)) % multiprecision
ans = 
    2.0115e-436

>> besselk(0,mp(10000))
ans = 
    1.4231e-4345

Do you have a simple example demonstrating your issue?


We have User's Manual with many details on how to use the toolbox: http://www.advanpix.com/documentation/users-manual/