Ilya Zakharevich on Fri, 12 Jul 2024 05:05:34 +0200


[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]

[BUG?] cannot vecsort without duplicates


BACKGROUND: I need to sort an array without duplicates — so that
	    consequent entries are not equal (in the sense that their
	    difference is not 0).  Mostly, I fail.  (I have a version
	    which passes all my tests, but I have no idea why it works¹⁾.)

  ¹⁾ I just do sorting as below, but twice.  (Do not know hwy this
     avoids the problems.)

BACKGROUND': sorting via lex() does not work by obvious reasons.  In
	     Math::Pari, I use gcmp() for <=> — but it is not
	     available in GP.  One can import() gcmp() — but
	     unfortunately, sorting with gcmp() STILL can produce
	     values with difference 0.  So I need to use a (silly?!)
	     sorting function below.

BACKGROUND'': in mid-70, William Kahan developed the math of IEEE
	      floating point.  The aim was to preserve as many laws of
	      arithmetic as possible — so that formulas of math which
	      work in ℝ have as wide range of applicability to FP as
	      possible.  I suspect that the rules of FP of PARI follow
	      pre-Kahan approach of “it looks cute”, leading to failures
	      much more often than what is known to be possible.

	        So: QUESTION (aside): were Kahan’s papers taken into account?

The failures are found by a quite complicated fuzzing process.  To
report the failing entry, I need a way to reproduce the failing
array. — However, PARI does not allow serialization except via
writebin() — and the result of writebin() segfaults on reading
(reported earlier today).

So today I decided to replace serialization by just reporting \x done
on the data POS1.  (See below.)  Then doing vecsort on POS1 as in:

  (19:27) gp > my(xRefl=1, yRefl=1, asX=1, asY=2, numcmp(x,y)=if(x==oo||y==oo,if(x<y,-1,x>y),x-y)); \
            POSR = vecsort(POS1, (x,y)->numcmp(x[asX],y[asX]), 8+4*(xRefl==-1))
  %5538 = List([[0.E-38, -0.75000000000000000000000000000000000000],
  	  	[0.E-38, 0.50000000000000000000000000000000000000],
  		[1.0000000000000000000000000000000000000, -0.25000000000000000000000000000000000000],
  		[2.0000000000000000000000000000000000000, -0.50000000000000000000000000000000000000]])

One can see that indeed 

  (19:29) gp > (POSR[1][1]-POSR[2][1])==0
  %5539 = 1
  
 — contrary to the sorting criterion.

Thanks,
Ilya

⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜⁜ The dump of POS1

(19:26) gp > POS1
%5537 = List([[0.E-38, -0.75000000000000000000000000000000000000], [2.0000000000000000000000000000000000000, -0.5000000000000000000000
0000000000000000], [1.0000000000000000000000000000000000000, -0.25000000000000000000000000000000000000], [0.E-38, 0.500000000000000000
00000000000000000000], [-3.673419846319648462 E-39, 0.99999999999999999999999999999999999999]])
(19:27) gp > \x
[&=00000000038b3090] LIST(lg=3,CLONE):2900000000000003 (subtyp=0,lmax=38):0000000004b80ce0 00000000047fdf20 00000000047fdfc0 000000000
4b80d70 00000000047fe060
  1st component = [&=0000000004b80ce0] VEC(lg=3,CLONE):2300000000000003 0000000004b80d18 0000000004b80cf8
    1st component = [&=0000000004b80d18] REAL(lg=2):0400000000000002 (0,expo=-128):1fffffffffffff80
    2nd component = [&=0000000004b80cf8] REAL(lg=4):0400000000000004 (-,expo=-1):dfffffffffffffff c000000000000000 0000000000000000
  2nd component = [&=00000000047fdf20] VEC(lg=3,CLONE):2300000000000003 00000000047fdf58 00000000047fdf38
    1st component = [&=00000000047fdf58] REAL(lg=4):0400000000000004 (+,expo=1):6000000000000001 8000000000000000 0000000000000000
    2nd component = [&=00000000047fdf38] REAL(lg=4):0400000000000004 (-,expo=-1):dfffffffffffffff 8000000000000000 0000000000000000
  3rd component = [&=00000000047fdfc0] VEC(lg=3,CLONE):2300000000000003 00000000047fdff8 00000000047fdfd8
    1st component = [&=00000000047fdff8] REAL(lg=4):0400000000000004 (+,expo=0):6000000000000000 8000000000000000 0000000000000000
    2nd component = [&=00000000047fdfd8] REAL(lg=4):0400000000000004 (-,expo=-2):dffffffffffffffe 8000000000000000 0000000000000000
  4th component = [&=0000000004b80d70] VEC(lg=3,CLONE):2300000000000003 0000000004b80da8 0000000004b80d88
    1st component = [&=0000000004b80da8] REAL(lg=2):0400000000000002 (0,expo=-129):1fffffffffffff7f
    2nd component = [&=0000000004b80d88] REAL(lg=4):0400000000000004 (+,expo=-2):5ffffffffffffffe ffffffffffffffff ffffffffffffffff
  5th component = [&=00000000047fe060] VEC(lg=3,CLONE):2300000000000003 00000000047fe098 00000000047fe078
    1st component = [&=00000000047fe098] REAL(lg=3):0400000000000003 (-,expo=-128):dfffffffffffff80 a000000000000000
    2nd component = [&=00000000047fe078] REAL(lg=4):0400000000000004 (+,expo=-1):5fffffffffffffff ffffffffffffffff fffffffffffffffb