weird conversion from QQ to RDF #14416

zimmermann6 · 2013-04-05T14:31:37Z

the following is weird:

sage: RDF(1/10)-RDF(1)/RDF(10)
-1.38777878078e-17

One would expect that the conversion from 1/10 to RDF is done as follows:

first convert 1 to RDF, which is exact
then convert 10 to RDF, which is exact
then divide RDF(1) by RDF(10)

For RR we get as a comparison:

sage: RR(1/10)-RR(1)/RR(10) 
0.000000000000000

More examples:

sage: for p in [1..10]:
....:     for q in [1..10]:
....:         if RDF(p/q) <> RDF(p)/RDF(q):
....:             print p, q
....:             
1 5
1 10
2 5
2 10
4 5
4 10
5 3
5 6
5 7
5 9
7 3
7 6
7 9
8 5
8 10
9 5
9 7
9 10
10 3
10 6
10 7
10 9

and for RR:

sage: for p in [1..10]: 
....:     for q in [1..10]:
....:         if RR(p/q) <> RR(p)/RR(q):
....:             print p, q
....:             
sage:

Apply attachment: 14416_QQ_to_RDF_v2.patch

Depends on #14335
Depends on #14336

CC: @robertwb

Component: basic arithmetic

Author: Paul Zimmermann, Jeroen Demeyer

Reviewer: Paul Zimmermann

Merged: sage-5.11.beta0

Issue created by migration from https://trac.sagemath.org/ticket/14416

The text was updated successfully, but these errors were encountered:

mwhansen · 2013-04-05T14:52:27Z

comment:1

Tracing things back, the conversion eventually gets done by:

mpq_get_d(self.value)

since Rational is just a wrapper around an mpq. This is where the weirdness seems to be.

zimmermann6 · 2013-04-05T14:58:19Z

comment:2

ok, then the explanation is that mpq_get_d (according to the GMP manual) rounds towards zero.

One could call mpz_get_d(numerator)/mpz_get(denominator) but then for large numerator and denominator one would get three roundings, instead of only one when calling mpq_get_d.

I propose to close this ticket.

Paul

kcrisman · 2013-04-05T17:21:52Z

comment:3

So it's okay that

sage: RDF(1/10)*10 == RDF(1)
False
sage: RDF(1/10)*10 - RDF(1)
-1.11022302463e-16

sage: RR(1/10)*10 == RR(1)  
True
sage: sage: RR(1/10)*10 - RR(1) 
0.000000000000000

as Thierry pointed out here? Just asking. Or is it RR that is behaving sub-optimally in this case? I assume that three roundings is ordinarily worse than one - my apologies for asking dumb questions.

tscrim · 2013-04-05T18:34:34Z

comment:4

RR is not the same as RDF in regards to how rounding is done, so I would say yes, it is okay. I'd suspect if you change the rounding mode, you should get a similar result, but rounding and computations in RDF are tied to your hardware (up to a certain precision) since I believe they are done in your native machine double type. In particular, the ALU (arithmetic logic unit) of whatever CPU you're using. I should state for the record I'm not 100% certain of this.

The reason why RR works uniformly is as it was noted in the ask sage, it is emulating a CPU in a program (one can think of it as constant hardware).

However instead of closing this ticket, I propose someone should implement a tutorial and/or improve the documentation and/or sage basics. I can do it if needbe.

tscrim · 2013-04-05T18:35:08Z

comment:5

Whoops, didn't mean to change the milestone.

sagetrac-tmonteil · 2013-04-05T19:13:28Z

comment:6

As Paul explained, the rounding is the same for RR and RDF: RDF rounds its computations to the nearest.

Actually, the problem is somewhere else: the conversion from QQ to RDF rounds towards zero, which is not consistent with the general behaviour of RDF.

IMHO, this conversion from QQ to RDF should be fixed, if possible.

tscrim · 2013-04-05T19:31:38Z

comment:7

Ah I see. I misread/misinterpreted Paul's explanation. I also agree that this should be fixed for consistency:

sage: RDF(RDF(1)/10) - RDF(1)/RDF(10)      
0.0
sage: RDF(RR(1)/10) - RDF(1)/RDF(10) 
0.0
sage: RDF(1/RDF(10)) - RDF(1)/RDF(10)
0.0
sage: RDF(1/float(10)) - RDF(1)/RDF(10)
0.0

zimmermann6 · 2013-04-05T20:50:14Z

comment:8

Karl-Dieter,

Replying to @kcrisman:

So it's okay that
sage: RDF(1/10)*10 == RDF(1)
False
sage: RDF(1/10)*10 - RDF(1)
-1.11022302463e-16

sage: RR(1/10)*10 == RR(1)  
True
sage: sage: RR(1/10)*10 - RR(1) 
0.000000000000000
as Thierry pointed out here? Just asking. Or is it RR that is behaving sub-optimally in this case? I assume that three roundings is ordinarily worse than one - my apologies for asking dumb questions.

what is annoying is that RDF and RR give different results. The fact that RR gives 0.0 in that case
is not that important, for example:

sage: RR(1/49)*49-RR(1)
-1.11022302462516e-16

Paul

zimmermann6 · 2013-04-05T21:06:23Z

comment:9

IMHO, this conversion from QQ to RDF should be fixed, if possible.

it is possible: first call the MPFR function mpfr_set_q on the rational fraction
(which is what RR(p/q) does) then convert back to double with mpfr_get_d, but this might be less efficient than the current code, since it would convert from QQ to RR, then from RR to RDF.

Converting separately the numerator and the denominator to RDF, then dividing in RDF does not work,
due to double rounding.
Consider for example the fraction p/q where p=2403806706169061971 and q=983883817941434958. If you first round p and q to RDF (say with mpz_get_d), then you get
2403806706169061888/983883817941435008, which is rounded to nearest to
5501555564164225/2251799813685248, whereas the direct rounding of p/q to nearest gives
2750777782082113/1125899906842624:

sage: p=2403806706169061971; q=983883817941434958; pp=RR(p); qq=RR(q)
sage: (pp/qq).exact_rational()
5501555564164225/2251799813685248
sage: RR(p/q).exact_rational()
2750777782082113/1125899906842624

Paul

robertwb · 2013-04-06T00:42:13Z

comment:10

RDF is for when you want your computations to be fast, at the expense of a little bit of accuracy/less control of rounding/platform dependence. As such, I think different behavior than RR is completely acceptable. It would be nice if mpq_get_d rounded towards nearest, but until someone implements that in a manner that's at least comparable in speed to mpq_get_d I think the current behavior is more in line with the philosophy of RDF than making things slow to get the last bit correct.

zimmermann6 · 2013-04-06T07:48:42Z

comment:11

Robert, would the following be ok in what concerns speed?

Let p/q the fraction to be converted to RDF, I assume p, q > 0.

multiply p or q by a power of 2 into pp and qq so that both have the same number of bits
(using mpz_sizeinbase and mpz_mul_2exp from GMP)
if pp < qq, multiply pp by 2⁵⁴, otherwise multiply pp by 2⁵³, using mpz_mul_2exp
now 2⁵³ <= pp/qq < 2⁵⁴, compute (using mpz_tdiv_qr) the quotient r = trunc(pp/qq) and the
remainder s; we know that r has exactly 54 bits
if r is even, return r*2^k where k is the normalizing constant (exact since r is exact on 53 bits)
if s is not zero, return (r+1)*2^k [r+1 is even, and exact on 53 bits, even in the
case where r+1 = 2⁵⁴]

5a) if s=0, return (r-1)*2^k if the bit of weight 1 of r is 0, otherwise (r+1)*2^k

I guess mpq_get_d does basically steps 1-4, thus it should be comparable in speed.

Paul

zimmermann6 · 2013-04-09T13:55:40Z

comment:12

here is some tentative code which converts from QQ to RDF with rounding to nearest.
In the usual case where both the numerator and the denominator are less than 2⁵³ in absolute value, it is about twice as fast as mpq_get_d.

double
mpq_get_d_nearest (mpq_t q)
{
  mpz_ptr a = mpq_numref (q);
  mpz_ptr b = mpq_denref (q);
  size_t sa = mpz_sizeinbase (a, 2);
  size_t sb = mpz_sizeinbase (b, 2);
  size_t na, nb;
  mpz_t aa, bb;
  double d;

  /* easy case: |a|, |b| < 2^53, no overflow nor underflow can occur */
  if (sa <= 53 && sb <= 53)
    return mpz_get_d (a) / mpz_get_d (b);

  /* same if a = m*2^e with m representable on 53 bits, idem for b, but beware  
     that both a and b do not give an overflow */
  na = sa - mpz_scan1 (a, 0);
  nb = sb - mpz_scan1 (b, 0);
  if (sa <= 1024 && na <= 53 && sb <= 1024 && nb <= 53)
    return mpz_get_d (a) / mpz_get_d (b);

  /* hard case */
  mpz_init (aa);
  mpz_init (bb);

  if (sa >= sb)
    {
      mpz_set (aa, a);
      mpz_mul_2exp (bb, b, sa - sb);
    }
  else
    {
      mpz_mul_2exp (aa, a, sb - sa);
      mpz_set (bb, b);
    }

  /* q = aa/bb*2^(sa-sb) */

  if (mpz_cmpabs (aa, bb) >= 0)
    {
      mpz_mul_2exp (bb, bb, 1);
      sa ++;
    }

  mpz_mul_2exp (aa, aa, 54);
  sb += 54;

  mpz_tdiv_qr (aa, bb, aa, bb);

  /* the quotient aa should have exactly 54 bits */

  if (mpz_tstbit (aa, 0) == 0)
    {
    }
  else if (mpz_cmp_ui (bb, 0) != 0)
    {
      if (mpz_sgn (aa) > 0)
        mpz_add_ui (aa, aa, 1);
      else
        mpz_sub_ui (aa, aa, 1);
    }
  else /* mid case: round to even */
    {
      if (mpz_tstbit (aa, 1) == 0)
        {
          if (mpz_sgn (aa) > 0)
            mpz_sub_ui (aa, aa, 1);
          else
            mpz_add_ui (aa, aa, 1);
        }
      else
        {
          if (mpz_sgn (aa) > 0)
            mpz_add_ui (aa, aa, 1);
          else
            mpz_sub_ui (aa, aa, 1);
        }
    }

  mpz_clear (aa);
  mpz_clear (bb);
  d = mpz_get_d (aa); /* exact */
  return ldexp (d, (long) sa - (long) sb);
}

However I don't know where to incorporate that code in Sage. Could someone help?

Paul

kcrisman · 2013-04-09T14:27:35Z

comment:13

Volker points to here in sage/rings/rational.pyx. I imagine that one could just replace that return mpq_get_d(self.value) with your code, or maybe make an auxiliary function since I suspect yours is pure C or C++ and this is a Cython file.

tscrim · 2013-04-09T14:45:16Z

comment:14

That looks like pure C code. Isn't there an mpz/mpq package with the underlying C(++) code, and wouldn't this go in there...?

jdemeyer · 2013-04-15T08:35:58Z

comment:48

Some tests might be machine-specific, but at least on most machines, all tests pass indeed.

jdemeyer · 2013-04-30T11:11:53Z

Dependencies: #14335, #14336

jdemeyer · 2013-04-30T20:33:04Z

comment:50

Rebased for doctest failures with #14336.

Paul: are you sure you don't want to review the patch? If the both of us look at the patch, that should be sufficient, no?

zimmermann6 · 2013-05-03T15:18:24Z

comment:51

I will try to review it next week. But if someone beats me, no problem!

Paul

zimmermann6 · 2013-05-07T06:19:50Z

Reviewer: Paul Zimmermann

zimmermann6 · 2013-05-07T06:19:50Z

comment:52

Jeroen, a few tiny comments:

in the following test, you can replace 20 by 13. Also, you could add all([RDF(q) == RR(q) for q in Q]) which would exercise the code for large numerators and denominators.

sage: Q = continued_fraction(pi, bits=3000).convergents()[20:]
sage: RDFpi = RDF(pi) 
sage: all([RDF(q) == RDFpi for q  in Q])

is the except? value really needed in mpq_get_d_nearest?
please replace round-to-even with round-to-nearest-even. Round to even and round to odd also exist.
please replace occured by occurred (several places, already mentioned)

Apart from that (and if tests still pass) I'm fine with the patch.

Paul

jdemeyer · 2013-05-07T06:43:38Z

comment:53

Replying to @zimmermann6:

Jeroen, a few tiny comments:

in the following test, you can replace 20 by 13.

Sure, but I wanted some safety margin.

Also, you could add all([RDF(q) == RR(q) for q in Q])

OK, but then for all covergents (before throwing away the first 20).

is the except? value really needed in mpq_get_d_nearest?

Yes, because sig_on() can throw exceptions.

zimmermann6 · 2013-05-07T07:47:55Z

comment:54

I could not apply the patch to 5.9 (after applying successfully the two dependencies):

----------------------------------------------------------------------
| Sage Version 5.9, Release Date: 2013-04-30                         |
| Type "notebook()" for the browser-based notebook interface.        |
| Type "help()" for help.                                            |
----------------------------------------------------------------------
sage: hg_sage.import_patch("/tmp/14416_QQ_to_RDF_v2.patch")
cd "/localdisk/tmp/sage-5.9/devel/sage" && sage --hg import   "/tmp/14416_QQ_to_RDF_v2.patch"
applying /tmp/14416_QQ_to_RDF_v2.patch
patching file sage/matrix/matrix_double_dense.pyx
Hunk #1 FAILED at 1018
Hunk #2 FAILED at 1947
Hunk #3 FAILED at 2087
3 out of 3 hunks FAILED -- saving rejects to file sage/matrix/matrix_double_dense.pyx.rej
patching file sage/plot/colors.py
Hunk #1 FAILED at 1301
1 out of 1 hunks FAILED -- saving rejects to file sage/plot/colors.py.rej
patching file sage/rings/contfrac.py
Hunk #1 FAILED at 632
1 out of 1 hunks FAILED -- saving rejects to file sage/rings/contfrac.py.rej
patching file sage/rings/rational.pyx
Hunk #1 FAILED at 73
Hunk #2 FAILED at 1999
Hunk #3 succeeded at 3660 with fuzz 2 (offset 228 lines).
2 out of 3 hunks FAILED -- saving rejects to file sage/rings/rational.pyx.rej
abort: patch failed to apply

Thus I cannot check all tests still pass. We'll have to rely on the testbot.

Paul

jdemeyer · 2013-05-07T08:58:45Z

comment:55

Attachment: 14416_QQ_to_RDF_v2.patch.gz

Replying to @zimmermann6:

I could not apply the patch to 5.9 (after applying successfully the two dependencies):

Since essentially every hunk fails, it looks like you're applying the patch on top of itself.

Anyway, I updated the patch with your suggestions.

jdemeyer · 2013-05-24T11:59:04Z

comment:57

Paul, any chance for a review again?

zimmermann6 · 2013-05-24T12:24:49Z

comment:58

Paul, any chance for a review again?

yes, will do.

zimmermann6 · 2013-05-27T08:26:58Z

comment:59

on top of Sage 5.9, I get doctest failures:

sage -t --long __init__.pyc  # AttributeError in doctesting framework
sage -t --long env.pyc  # AttributeError in doctesting framework
sage -t --long misc/interpreter.py  # 1 doctest failed
sage -t --long misc/trace.py  # 2 doctests failed
sage -t --long tests/cmdline.py  # 11 doctests failed
sage -t --long version.pyc  # AttributeError in doctesting framework
sage -t --long tests/interrupt.pyx  # Time out

Paul

jdemeyer · 2013-05-27T08:31:02Z

comment:60

None of these failures look related to this ticket, what does a "clean" Sage 5.9 (without any applied patches) give? Please attach the actual doctest failures, not just the summary at the end, otherwise it is impossible to find out what went wrong.

zimmermann6 · 2013-05-27T08:38:14Z

comment:61

maybe the doctest failures are due to some interaction with a spkg I installed in another clone of Sage (I installed the patches needed for #9880, and I believe there is side effect from one clone to the other ones for installed spkgs). Is there any way to get a "clean" Sage 5.9 without recompiling the sources again?

Anyway my remarks from comment [comment:52] are taken into account, thus provided all tests pass with the testbot, I give a positive review.

Paul

jdemeyer · 2013-06-06T12:32:33Z

Merged: sage-5.11.beta0

zimmermann6 added this to the sage-5.10 milestone Apr 5, 2013

zimmermann6 added c: basic arithmetic labels Apr 5, 2013

zimmermann6 assigned aghitza Apr 5, 2013

zimmermann6 added the s: needs review label Apr 5, 2013

tscrim added s: needs info and removed s: needs review labels Apr 5, 2013

tscrim modified the milestones: sage-5.10, sage-5.9 Apr 5, 2013

tscrim modified the milestones: sage-5.9, sage-5.10 Apr 5, 2013

tscrim added c: basic arithmetic and removed c: documentation labels Apr 5, 2013

jdemeyer added s: needs work and removed s: needs review labels Apr 30, 2013

jdemeyer added s: needs review and removed s: needs work labels Apr 30, 2013

zimmermann6 added s: needs work and removed s: needs review labels May 7, 2013

jdemeyer added s: needs review and removed s: needs work labels May 7, 2013

jdemeyer modified the milestones: sage-5.10, sage-5.11 May 30, 2013

jdemeyer added s: positive review and removed s: needs review labels May 30, 2013

jdemeyer removed the s: positive review label Jun 6, 2013

jdemeyer closed this as completed Jun 6, 2013

jdemeyer mentioned this issue May 29, 2014

Fix rounding ZZ -> Python float #16385

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

weird conversion from QQ to RDF #14416

weird conversion from QQ to RDF #14416

zimmermann6 commented Apr 5, 2013

mwhansen commented Apr 5, 2013

zimmermann6 commented Apr 5, 2013

kcrisman commented Apr 5, 2013

tscrim commented Apr 5, 2013

tscrim commented Apr 5, 2013

sagetrac-tmonteil mannequin commented Apr 5, 2013

tscrim commented Apr 5, 2013

zimmermann6 commented Apr 5, 2013

zimmermann6 commented Apr 5, 2013

robertwb commented Apr 6, 2013

zimmermann6 commented Apr 6, 2013

zimmermann6 commented Apr 9, 2013

kcrisman commented Apr 9, 2013

tscrim commented Apr 9, 2013

jdemeyer commented Apr 15, 2013

jdemeyer commented Apr 30, 2013

jdemeyer commented Apr 30, 2013

zimmermann6 commented May 3, 2013

zimmermann6 commented May 7, 2013

zimmermann6 commented May 7, 2013

jdemeyer commented May 7, 2013

zimmermann6 commented May 7, 2013

jdemeyer commented May 7, 2013

jdemeyer commented May 24, 2013

zimmermann6 commented May 24, 2013

zimmermann6 commented May 27, 2013

jdemeyer commented May 27, 2013

zimmermann6 commented May 27, 2013

jdemeyer commented Jun 6, 2013

weird conversion from QQ to RDF #14416

weird conversion from QQ to RDF #14416

Comments

zimmermann6 commented Apr 5, 2013

mwhansen commented Apr 5, 2013

zimmermann6 commented Apr 5, 2013

kcrisman commented Apr 5, 2013

tscrim commented Apr 5, 2013

tscrim commented Apr 5, 2013

sagetrac-tmonteil mannequin commented Apr 5, 2013

tscrim commented Apr 5, 2013

zimmermann6 commented Apr 5, 2013

zimmermann6 commented Apr 5, 2013

robertwb commented Apr 6, 2013

zimmermann6 commented Apr 6, 2013

zimmermann6 commented Apr 9, 2013

kcrisman commented Apr 9, 2013

tscrim commented Apr 9, 2013

jdemeyer commented Apr 15, 2013

jdemeyer commented Apr 30, 2013

jdemeyer commented Apr 30, 2013

zimmermann6 commented May 3, 2013

zimmermann6 commented May 7, 2013

zimmermann6 commented May 7, 2013

jdemeyer commented May 7, 2013

zimmermann6 commented May 7, 2013

jdemeyer commented May 7, 2013

jdemeyer commented May 24, 2013

zimmermann6 commented May 24, 2013

zimmermann6 commented May 27, 2013

jdemeyer commented May 27, 2013

zimmermann6 commented May 27, 2013

jdemeyer commented Jun 6, 2013