Inconsitent Euclidean operations in Z

[fieker]

On Singular-Spielwiese-Shell:

(-7) / 5 -> (warning), -1 (-7) div 5 -> -2

I propose to enforce and document for numbers that nIntDiv, nDiv, nMod (and the still missing nDivMod) are all completely compatible, ie. for a fixed choice of Euclidean division in Z, be it positive remainders, or symmetric (I don't care but it should be documented and enforced) I'd like

a = qb + r

and nDiv(a,b) = nIntDiv(a,b) = q = nDiv(a-r,b)

r = a - nDiv(a, b)* b

and nDivMod(a,b, &q, &r) to always, regardless of signs, return the same (compatible) results. Looking at the code, it might be that this actually stems from mpz in which case the singular interface should correct for it.

I also propose that the intended behaviour is documented in some header file.

I am guessing that nExactDiv will throw an error if there is a remainder.

I suggest to rename n_Neg to nInpNeg (and to write a proper nNeg) to be "in-line" with other Inp* functions.

I'd also be happy to actually (try to) do it myself.

[Oleksandr]

Thanks for your feedback and proposal! You are welcome to propose fixes and submit patches (for separate issues).

Besides, we discussed this issue yesterday on Singular-Spielwiese meeting and decided that I am to correct nIntDiv so that it is compatible with nIntMod for bigints. But I am still playing around with it.

[Oleksandr]

I'm testing Singular using the following code:

list L7 = list( int(7), bigint(7), bigint(7000000000000000), int(2^28), bigint(2^28) );
list L5 = list( int(5), bigint(5), bigint(5000000000000000), int(-1), bigint(-1) );
proc Test(a, b)
{
  def q = a div b;
  def r = a mod b;

  if( r < 0 ) // never fails!
  {
    typeof(a), a, "/", typeof(b), b, "::", "/:", a/b, ", q: ", q, ", r: ", r;
    "r0: ", r >= 0;
  }

  if( (a - b*q) != r ) // sometimes fails!
  {
    typeof(a), a, "/", typeof(b), b, "::", "/:", a/b, ", q: ", q, ", r: ", r;
    "==: ", (a - b*q) == r, ", diff: ", (a - q * b), " <> r: ", r; 
    "q must be: ", (a - r) div b;
  }
}
int i,j;
for( i = 1; i <= size(L7); i++ )
{
  for( j = 1; j <= size(L5); j++ )
  {
    Test( L7[i], L5[j] );
    Test( -L7[i], L5[j] );
    Test( L7[i], -L5[j] );
    Test( -L7[i], -L5[j] );
  }
}

the result is as follows:

int -7 / bigint 5 :: /: -1 , q:  -1 , r:  3
==:  0 , diff:  -2  <> r:  3
q must be:  -2

int -7 / bigint -5 :: /: 1 , q:  1 , r:  3
==:  0 , diff:  -2  <> r:  3
q must be:  2

int -7 / bigint 5000000000000000 :: /: 0 , q:  0 , r:  4999999999999993
==:  0 , diff:  -7  <> r:  4999999999999993
q must be:  -1

int -7 / bigint -5000000000000000 :: /: 0 , q:  0 , r:  4999999999999993
==:  0 , diff:  -7  <> r:  4999999999999993
q must be:  1

bigint -7 / int 5 :: /: -1 , q:  -1 , r:  3
==:  0 , diff:  -2  <> r:  3
q must be:  -2

bigint -7 / int -5 :: /: 1 , q:  1 , r:  3
==:  0 , diff:  -2  <> r:  3
q must be:  2

bigint -7 / bigint 5 :: /: -1 , q:  -1 , r:  3
==:  0 , diff:  -2  <> r:  3
q must be:  -2

bigint -7 / bigint -5 :: /: 1 , q:  1 , r:  3
==:  0 , diff:  -2  <> r:  3
q must be:  2

bigint -7 / bigint 5000000000000000 :: /: 0 , q:  0 , r:  4999999999999993
==:  0 , diff:  -7  <> r:  4999999999999993
q must be:  -1

bigint -7 / bigint -5000000000000000 :: /: 0 , q:  0 , r:  4999999999999993
==:  0 , diff:  -7  <> r:  4999999999999993
q must be:  1

bigint 7000000000000000 / bigint -5000000000000000 :: /: -2 , q:  -2 , r:  2000000000000000
==:  0 , diff:  -3000000000000000  <> r:  2000000000000000
q must be:  -1

bigint -7000000000000000 / bigint -5000000000000000 :: /: 1 , q:  1 , r:  3000000000000000
==:  0 , diff:  -2000000000000000  <> r:  3000000000000000
q must be:  2

int -268435456 / bigint 5 :: /: -53687091 , q:  -53687091 , r:  4
==:  0 , diff:  -1  <> r:  4
q must be:  -53687092

int 268435456 / bigint -5 :: /: -53687092 , q:  -53687092 , r:  1
==:  0 , diff:  -4  <> r:  1
q must be:  -53687091

int -268435456 / bigint -5 :: /: 53687091 , q:  53687091 , r:  4
==:  0 , diff:  -1  <> r:  4
q must be:  53687092

int -268435456 / bigint 5000000000000000 :: /: 0 , q:  0 , r:  4999999731564544
==:  0 , diff:  -268435456  <> r:  4999999731564544
q must be:  -1

int 268435456 / bigint -5000000000000000 :: /: -1 , q:  -1 , r:  268435456
==:  0 , diff:  -4999999731564544  <> r:  268435456
q must be:  0

int -268435456 / bigint -5000000000000000 :: /: 0 , q:  0 , r:  4999999731564544
==:  0 , diff:  -268435456  <> r:  4999999731564544
q must be:  1

bigint -268435456 / int 5 :: /: -53687091 , q:  -53687091 , r:  4
==:  0 , diff:  -1  <> r:  4
q must be:  -53687092

bigint 268435456 / int -5 :: /: -53687092 , q:  -53687092 , r:  1
==:  0 , diff:  -4  <> r:  1
q must be:  -53687091

bigint -268435456 / int -5 :: /: 53687091 , q:  53687091 , r:  4
==:  0 , diff:  -1  <> r:  4
q must be:  53687092

bigint -268435456 / bigint 5 :: /: -53687091 , q:  -53687091 , r:  4
==:  0 , diff:  -1  <> r:  4
q must be:  -53687092

bigint 268435456 / bigint -5 :: /: -53687092 , q:  -53687092 , r:  1
==:  0 , diff:  -4  <> r:  1
q must be:  -53687091

bigint -268435456 / bigint -5 :: /: 53687091 , q:  53687091 , r:  4
==:  0 , diff:  -1  <> r:  4
q must be:  53687092

bigint -268435456 / bigint 5000000000000000 :: /: 0 , q:  0 , r:  4999999731564544
==:  0 , diff:  -268435456  <> r:  4999999731564544
q must be:  -1

bigint 268435456 / bigint -5000000000000000 :: /: -1 , q:  -1 , r:  268435456
==:  0 , diff:  -4999999731564544  <> r:  268435456
q must be:  0

bigint -268435456 / bigint -5000000000000000 :: /: 0 , q:  0 , r:  4999999731564544
==:  0 , diff:  -268435456  <> r:  4999999731564544
q must be:  1

[Oleksandr]

the problem is in nlIntDiv which is not compatible with nlIntMod.

moreover it is to be fixed at 2 places: for small "ints" and for GMP numbers, where we used to use mpz_tdiv_q.

ps: in the above output switch to mpz_fdiv_q.

see also ​http://gmplib.org/manual/Integer-Division.html for GMP division functions

[Oleksandr]

Looking at the code of nlIntMod I can see an additional timebomb: as far as I know, maybe it's different in C++, the behaviour of % is not completely defined. Ie. you cannot rely on % to give positive results on all architectures....

For future reference let me copypaste ​http://en.wikipedia.org/wiki/Modulo_operation:

C++ ISO 1998 Result of % is Implementation defined:

ISO/IEC 14882:2003 : Programming languages -- C++. 5.6.4: ISO, IEC. 2003. "the binary % operator yields the remainder from the division of the first expression by the second. .... If both operands are nonnegative then the remainder is nonnegative; if not, the sign of the remainder is implementation-defined".

C++ ISO 2011: Result of % has the same sign as the Dividend

[Oleksandr]

see also ​http://www.singular.uni-kl.de:8002/trac/ticket/442