arf_get_d() is rounding incorrectly

[saraedum]

I might be doing something wrong here but it seems to me that arf_get_d() is getting some roundings wrong with ARF_RND_NEAR.

In the following example I take an arf element a that came up randomly in https://github.com/videlec/e-antic/pull/114 and round it with arf_get_d(ARF_RND_NEAR) which gives x. However, the double y is closer to a so I think arf_get_d should have returned that double instead.

#include <arf.h>

const char* dump = "39d570dfa5946f0ee839305ed4a4767b8e3fd7640c0d3b3aaed2eb693a69ec9d8be507545d11486e31a814a2750b7af22c942c855aaa8963adfe1c90563ed1d57218c1d9da25d41df26d17a23962eeef7d268d35b607cc9bd7bb4d5f457d3 -2e7";

int main() {
  arf_t a;

  arf_init(a);
  arf_load_str(a, dump);

  printf("a = ");
  arf_debug(a);
  printf(" = ");
  arf_print(a);
  printf("\n");

  double d = arf_get_d(a, ARF_RND_NEAR);

  printf("get_d(a) = %a\n", d);

  {
      arf_t ax, ay;
      double x = 8139377218865719./4398046511104.;
      double y = 1017422152358215./549755813888.;

      arf_init(ax);
      arf_init(ay);

      arf_set_d(ax, x);
      arf_set_d(ay, y);

      printf("x = %a = ", x);
      arf_print(ax);
      printf(", y = %a = ", y);
      arf_print(ay);
      printf(", i.e., get_d(a) == x\n");

      arf_sub(ax, a, ax, 2048, ARF_RND_NEAR);
      arf_sub(ay, ay, a, 2048, ARF_RND_NEAR);

      printf("a - ax = ");
      arf_print(ax);
      printf("\n");
      printf("ay - a = ");
      arf_print(ay);
      printf("\n");

      printf("i.e., a is closer to ay so get_d(a) should have been y.");
  }
}

The output of this program is:

a = {exp=11; size=12; sgnbit=0; digits=[ 6840182220401197056 17625458189021426287 14534345839491136445 14439392970997715063 13256471128156358485 12848965541941618062 14312530367707212744 3428397452016034232 13496071632152146550 4107104259050630378 11593604076168927726 16669444544236993595]} = (85628460650401049251127855502029414724934253826929216416186643198844579534836066518350052226030633770627919100392797002027482544691237435538601045912922419435288490067010919856733222580924798494640070427510850119619385654269907 * 2^-743)
get_d(a) = 0x1.ceab86fd2ca38p+10
x = 0x1.ceab86fd2ca37p+10 = (8139377218865719 * 2^-42), y = 0x1.ceab86fd2ca38p+10 = (1017422152358215 * 2^-39), i.e., get_d(a) == x
a - ax = (5566438600596433074713521695163990960379759211712698385474929106159744847177898243282594905260413291056794690186073761436849961449387020494086468491054171774665449146745841172172380715633413360899934303884171219 * 2^-743)
ay - a = (4953833202500313939768458070596266370720920393933649333521632699977719461416263400944738167294763611397171247526282673989188903050980587232169161050249589925244998195511048024210946800135232073642652199587391533 * 2^-743)

[saraedum]

Hm…on the last run of the program it suddenly works. Sorry for the noise. Let me investigate what I was doing wrong.

[saraedum]

Sorry, my bad. I was comparing pointers instead of values in an earlier version of the above reproducer. Everything seems to work fine.