INT64 arithmetic in MATLAB

Enables int64 Addition, subtraction, multiplication, division and modulus.

Description

This submission enables the following operations for the int64 and uint64 data types:
* Addition
* Subtraction
* Multiplication (element-wise and matrix)
* Division (element-wise only)
* mod, abs, bitshift

See the published file for more details.

• MATLAB release: MATLAB 7.4 (R2007a)
• Other requirements: C compiler supporting "signed long long".

[Bruno Luong]

After a quick look I have few comments that I wish the author could correct:

As stated, the C Mex can't be compiled by the popular MS VISUAL C (type LONG LONG must be replaced by __int64)

It is dangerous when no checking is performed after mxDuplicateArray that could returns NULL and results in a crash if there is not enough memory.

[Bruno Luong]

It is also puzzling this can work

a=int64([0 1; 2 3])
b=int64([4 5 6 7]')
a+b

[Petter]

I used Microsoft Visual C++ Express 2008 to compile the library. "long long" _is_ supported in later versions of Visual C++. Also, "long long" is more supported by other compilers and will probably end up in the next C++ standard as well.

[Bruno Luong]

What a relief to have finally a basic int64 arithmetic operators supported.

[Derek O'Connor]

This seems to be an excellent package. I used it to test matrix multiplication
and obtained very surprising results : multiplying two 10^4x10^4 int64 matrices took about 2 secs while two doubles took about 35 secs. The double mat-mult used the Math Kernel and all 8 cores while the int64 did not use the Math Kernel.

Also, a 2-degree polyfit of the int64 mat-mult times was very good, i.e., t(n) is O(n^2)!

Can anyone enlighten me? Here is the test function:

function times = TimesInt64(v,degree);
% Test matrix multiplication using int64arithmetic
% int64arithmetic compiled with Microsoft Visual C++ 2008
% Dell Precision 690, 2xQuadcore Xeon 5345, 2.3GHz, 16GB
% Windows Vista 64
% Example : >> times = TimesInt64(1:15,2);

% Derek O'Connor Sept 2009

nvals = v'*10^3;
times = zeros(length(nvals),1);
c = intmax('int64');
for k = 1:length(nvals)
    n = nvals(k);
    A64=int64(floor(c*rand(n,n)));
    B64=int64(floor(c*rand(n,n)));

    tic; A64*B64; times(k) = toc;
end;

p=polyfit(nvals,times,degree);
ptimes = polyval(p,nvals);
plot(nvals,times,'.',nvals,ptimes,'-')

% End Function TimesInt64(v,degree)

% n times (secs)
% 1000 0.0176
% 2000 0.06972
% 3000 0.15354
% 4000 0.34136
% 5000 0.50145
% 6000 0.63655
% 7000 0.9838
% 8000 1.2829
% 9000 1.5717
% 10000 2.0059
% 11000 2.4728
% 12000 3.0439
% 13000 3.4914
% 14000 4.2588
% 15000 4.9525

% p(x) = 0.098649 - 5.6075e-005*x + 2.5027e-008*x^2

[John D'Errico]

Derek - it points out that even though the function mtimes is provided, this package only claims to support ELEMENTWISE operations.

The quadratic behavior tells me that this code did not do a true matrix multiply. Did you check that A64*B64 was correct in your test? Element wise multiplication will take O(N^2) flops, whereas a matrix multiply (as * is supposed to generate) is O(N^3).

I would have preferred that an error be generated if you use * on a pair of matrices, when that operation will not return the proper result. This would normally make me downrate this code. Since I cannot test it without compiling it, I won't give any rating at all.

[Derek O'Connor]

John -- Thanks for pointing out that Int64 matrix multiplication is element-wise. I should have been more careful, given that the description clearly states that "multiplication and division are element-wise only". This of course explains the O(n^2) behavior.

On a more general point, I wish there was a consistent standard for matrix multiplicaton. For example, in O-Matrix, C= A*B is the same as Matlab, but C=A^2 gives cij = (aij)^2. This is really confusing. In Fortran 90, C=A*B gives cij = aij*bij.

[Petter]

John, the reason the * operator is included at all is because it is desirable to support operations like
a = int64(2);
b = int64( [1 2; 3 4] );
c = a * b

That is the reason * was included. The next update will issue an error if both a and b are matrices. Thank you for your suggestions.

[John D'Errico]

I can't test this myself, so I cannot rate it. But if it now prevents A*B when A and B are both matrices, neither of which is a scalar, then this is good if it would otherwise have done the wrong thing. Better of course is to implement array multiplication too (perhaps in a future version.) Regardless, there have been a few times recently when I wanted to have int64 operations defined, so this is a good addition to MATLAB.

[Petter]

Array multiplication .* is already implemented.

[Petter]

Never mind, I misunderstood what you meant. Implementing a crude A*B is very easy, so I might do that.

[Jan Simon]

A really needed and useful submission!
The efficiency can be improved significantly:
1. checkinputs: "if length(size(a)) ~= length(size(b)) ... if ~all(size(a) == size(b))" can be combined to: "if isequal(size(a), size(b))". Performing the check in the Mex and omitting the intermediate call of the M-functions would reduce the overhead.
2. (u)int64abs: At first the input is duplicated to plhs[0]. Then it is not necessary to copy positive values from the input to the output. A faster approach:
  void mexFunction(...) {
  signed long long *c, *cf;
  plhs[0] = mxDuplicateArray(prhs[0]);
  if (!plhs[0]) ...
  c = (signed long long *) mxGetData(plhs[0]); // not mxGetPr
  cf = c + mxGetNumberOfElements(plhs[0]);
  for ( ; c < cf; c++) {
     if (*c < 0L) {*c = -*c; }
  }
  }
3. Similar improvements for plus, minus, times, rdived: If the 1st input is duplicated already, the operation does not need to access the values from the inputs again. So these lines in (u)int64operations.c:
  c[i] = a[i] OPERATOR b[i]
can be changed to the faster:
  c[i] OPERATOR b[i]
using e.g. "+=" as operator instead of "+".
4. TMW suggests to use mxGetData instead of mxGetPr to get a pointer to non-DOUBLE arrays. But at least both method works currently.
5. Matlab's MOD can handle 0 as second input, and a user might expect this for the (U)INT64 version also.
6. "long long" is supported by the LCC compiler shipped with Matlab. For (all) other compilers, you can use uint64_T and int64_T as defined in "tmwtypes.h".
7. The check for division by zero can be included in the same loop as the operation to save time.

[Petter]

1. The reason the check is made in an .m file is that otherwise the MEX files would have to be quite verbose with lots of conversions from various types to int64.
2. I'll implement this optimization
3. Same here
4. OK, I'll change that
5. OK, I guess I could follow Matlab here
6. long long works on GCC, MSVC++ and LCC and feels more standard
7. True, I'll change that

Thanks for the suggestions.

[Igor]

Thanks a lot for this. Could you also add support for _colonobj so that one doesn't get a "Undefined function or method '_colonobj' for input arguments of type 'uint64'" when writing 1:n

[Jan Simon]

Thank you for the fast update!

[Jan Simon]

Although INT64 values can be very large, you can still get an overflow. While Matlab uses saturated arithmetic for integer types, your tools do not. It is important to mention this in the documentation. Implementing a saturation (means: value is the max or min INT value if it exceeds the valid range) is not trivial, and omitting it is no reason to reduce the 5 star rating.

You mention precompiled in the ReadMe - where can I find them?
Following Loren's block about inplace operations, I suggest to replace e.g.:
  function c = abs(a) c = a; end
by: function a = abs(a), end

[Vincent]

Hi Petter,
Thanks for your great package.

It didn't compile out of the box, however.
I tried to compile this on a Debian Stable machine, having gcc-4.3.2 installed.
The compilation of the first mex function aborts with the following error:

int64operation.c:38:1: error: pasting "OPERATOR" and "=" does not give a valid preprocessing token

I had to apply the following modifications:
Current situation:
One macro in each operator source files (#define OPERATOR +) and the OPERATOR_EQ defined as OPERATOR ## = in the (u)int64operation.c source files.
New situation:
Use two separate macro functions.
For example, the contens of int64plus.c is now:
---
#define UNARY_OPERATOR(a,b) ((a) += (b))
#define BINARY_OPERATOR(a,b) ((a) + (b))
#include "int64operation.c"
---
Where the OPERATOR and OPERATOR_EQ were used in int64operation.c:
---
c[i] OPERATOR_EQ b[i];
becomes
UNARY_OPERATOR(c[i], b[i]);
---
c[i] = a[0] OPERATOR b[i];
becomes
c[i] = BINARY_OPERATOR(a[0], b[i]);
---

I don't know if this works for every compiler around there, but for me, this works like a charm...
If you're interested, I can send you my modified version to the address specified in the README.txt file, so that you can more easily update your package.

[Petter]

Vincent: Interesting. I thought the standard allowed the pasting of operators. Is this not true?

[Jon]

Seamless use. Well done Petter and thanks!