0

This question has been edited to clarify it.

I have the following matrix, defined in page 1 of [reteam.org/papers/e59.pdf], written in R notation:

m1 = matrix(c(207560540,956631177,1,956631177,2037688522,1,2037688522,1509348670,1),ncol=3,byrow=T)

The determinant of m1 should be a integer multiple of 2^31 -1.

As indicated in accepted answer, det(m1) should be -1564448852668574749.

However, in R, I got

> det(m1)
[1] -1564448852668573184

and, using a simple equation by hand:

> m1[1,1]*(m1[2,2]-m1[3,2]) - m1[2,1]*(m1[1,2] - m1[3,2]) + m1[3,1]*(m1[1,2]- m1[2,2])
[1] -1564448852668574720

As indicated in accepted answer, the correct determinant is obtained and checked by:

#include <inttypes.h>
#include <stdio.h>

int main() {

int64_t m1[3][3] = {{INT64_C(207560540) , INT64_C(956631177)   , INT64_C(1)},{ INT64_C(956631177), INT64_C(2037688522),       INT64_C(1)},{INT64_C(2037688522), INT64_C(1509348670) ,   INT64_C(1)}};                                                         

int64_t dm1 = m1[0][0]*(m1[1][1]-m1[2][1]) - m1[1][0]*(m1[0][1] - m1[2][1]) + m1[2][0]*(m1[0][1]- m1[1][1]);
int64_t divisor = (INT64_C(1)<<31) -1;
int64_t tmp = dm1/divisor;
int64_t check = tmp * divisor;

printf("dm1 == %" PRIu64"\n",dm1);
printf("(dm1/(2^31 -1))* %" PRIu64 " == %" PRIu64 "\n", divisor, check);

}

The following text is the old question. The main error was using a unsigned type.

My old minimum non-working code example was:

  #include <inttypes.h>
  #include <stdio.h>

  int main() {

  uint64_t dm1 = 1564448852668573184;
  uint64_t divisor = (UINT64_C(1)<<31) -1; //powl(2,31)-1;
  uint64_t tmp = dm1/divisor;
  uint64_t check = tmp*divisor;                                                                                      

  printf("dm1 == %" PRIu64"\n",dm1);
  printf("(dm1/(2^31 -1))* %" PRIu64 " == %" PRIu64 "\n", divisor, check);
}

Its output is 

dm1 == 1564448852668573184

(dm1/(2^31 -1))* 2147483647 == 1564448850521091102

The problem is that the value of the second line should be equal to the one in the first line.

What is my mistake? How can I make this work?

Marc Kees
  • 206
  • 2
  • 15

2 Answers2

0

For the same reason that (using integers) you get 10 / 3 * 3 = 9. There is a remainder of the division. 10 / 3 = 3, with remainder 1. When you multiply by 3, the remainder is lost, so you get 9 instead of 10.

In this case, your remainder is 2147482082, which, when added to 1564448850521091102, gives 1564448852668573184. Try this:

uint64_t dm1 = 1564448852668573184;
uint64_t divisor = (UINT64_C(1)<<31) -1; //powl(2,31)-1;
uint64_t tmp = dm1/divisor;
uint64_t remainder = dm1%divisor;
uint64_t check = tmp*divisor+remainder;

And you should get the correct result.

BlueMoon93
  • 2,910
  • 22
  • 39
0

The numerator is not an exact multiple of the divisor, so there is a remainder, and the quotient is truncated.

1564448852668573184 / 2147483647 = 728503266 remainder 2147482082

Multiplying back,

2147483647 * 728503266 + 2147482082 = 1564448852668573184

EDIT:

The determinant of the 3x3 matrix shown in your linked reference is -1564448852668574749. This is exactly divisible by 2147483647 to give -728503267.

So you have an arithmetic overflow somewhere.

ANSWER:

The value of the matrix determinant in your linked example is negative. Please use int64_t instead of uint64_t.

Weather Vane
  • 33,872
  • 7
  • 36
  • 56
  • I got this number 1564448852668573184 from section 2 [reteam.org/papers/e59.pdf]. I computed the determinant dm1 of the first matrix in there. The author says this is to be an integer multiple of 2^31 -1. – Marc Kees Apr 04 '16 at 14:18
  • The number `1564448852668573184` does not appear in the paper you quoted. My calculator says it does not divide exactly. – Weather Vane Apr 04 '16 at 14:20
  • No it doesn't appear. I computed it by the determinant from the matrix appearing in the first page. – Marc Kees Apr 04 '16 at 14:22
  • Before I was computing determinant outside C. Here is now how I compute: `uint64_t m1[3][3] = {{UINT64_C(207560540) , UINT64_C(956631177) , UINT64_C(1)},{ UINT64_C(956631177), UINT64_C(2037688522), UINT64_C(1)},{UINT64_C(2037688522), UINT64_C(1509348670) , UINT64_C(1)}};` This results in determinant dm1=16882295221040976867 which still is not multiple of 2^31 -1. I suspect I have overflow problems while computing my determinant. – Marc Kees Apr 04 '16 at 14:54
  • Use `int64_t` not unsigned. The values just fit the range, – Weather Vane Apr 04 '16 at 15:00
  • @WeatherVane DO NOT use `int64_t`, signed integer overflow is UB in C – Severin Pappadeux Apr 04 '16 at 16:48
  • @SeverinPappadeux What can I use instead of int64_t to avoid UB? – Marc Kees Apr 04 '16 at 16:52
  • @SeverinPappadeux did you downvote? The signed arithmetic did not overflow. That is how I found the determinant and I checked it with a bigint library too. The sum of the positive terms is `3816150093823856864` and the (positive) sum of the negative terms is `5380598946492431613`. Neither overflows 64 bit signed int. – Weather Vane Apr 04 '16 at 16:52
  • @MarcKees just ignore him. He is wrong in this case, but with larger numbers you'd need a big int library. I was proving the point that your determinant calculation was in error. – Weather Vane Apr 04 '16 at 16:59
  • @WeatherVane Before I thought this would be a good question. However thanks to you, now I understand what happened. It seems my question does not make sense anymore. Should I delete it? – Marc Kees Apr 04 '16 at 17:03
  • Please don't delete it. I answered your problem by saying you used unsigned instead of signed, when the result is negative. Regardless of the size of the operands, this is true. – Weather Vane Apr 04 '16 at 17:08
  • @WeatherVane it might be ok for this particular data nput, but when there is a general expression which could overflow NEVER use signed integers, code would be undefuned – Severin Pappadeux Apr 04 '16 at 17:15
  • @SeverinPappadeux it's not a general expression, it was a particular expression for this particular problem, which solved it. So what type would you suggest? Should I use your all-embracing type on every calculation, since you say I must never use signed arithmetic? Do you imply that unsigned arithmetic is ok? I know its wrapping is legal, but that does not prevent it giving a wrong amswer. – Weather Vane Apr 04 '16 at 17:16
  • @SeverinPappadeux good grief. I am struggling to take you seriously. We all know that arithmetic can overflow. This doesn't. – Weather Vane Apr 04 '16 at 17:24
  • @WeatherVane the stated problem is integer matrix DETERMINANT calculation. This is wrong kind of advice for determinant calculation to switch to integer type which could cause UB in case of overflow - pure and simple. Dixi! – Severin Pappadeux Apr 04 '16 at 17:38
  • @SeverinPappadeux you haven't even said what data type you would use. Please explain how, if `double`, when the significance exceed 63 bits, the modulus could be a whole number? Come on, post a better answer or undo the downvote! – Weather Vane Apr 04 '16 at 17:40
  • @MarcKees I briefly looked into paper - you need integer math, so use of doubles is pretty much no-no. Either continue with with `uint64_t` with explicit overflow control, or something like GNU MP. – Severin Pappadeux Apr 04 '16 at 18:58
  • @SeverinPappadeux OP did not use `double`. And `uint64_t` was the cause of his troubles, because he needed a signed type. Please get up to speed. My answer uses integer math. I ask again, what would you have done? You are talking utter nonsense. – Weather Vane Apr 04 '16 at 19:04