[epsilon-devel] Tagged floating-point corrections, and SPARC [Was: Re: Unboxed tagged floating-point data in C: a (new?) efficient and portable solution]

[Luca Saiu]

On 2019-04-13 at 13:05 +0200, Luca Saiu wrote:

> On 2019-04-13 at 12:52 +0200, Luca Saiu wrote:
>
>> It compiles and works with every architecture I have played with
>> (x86_64, i386, mips, powerpc, riscv, sh, alpha, where applicable on
>> different word size and endianness).
>
> I had always forgotten to test on arm and aarch64, for some reason.
> [...]

Another architecture I had forgotten to test was SPARC.  Old versions
had the same limitation as i386 but I've found about the new
movstosw/movdtox instructions in recent specifications, and GCC
generates these at least with -mcpu=m8.  Even if I cannot easily test M8
code with QEmu the code looks correct.

About the first message, here are a few obvious corrections: [FIXME: did
the anonymous friend point out these as well?  Give thanks in public if
he did]

#define f2u(a)                   \
  (((union conv) {.f = (a)}).u)

#define u2f(a)                   \
  (((union conv) {.u = (a)}).f)

I'm attaching an updated copy of the testing file, including the fixes
above plus others.

-- 
Luca Saiu
* GNU epsilon:           http://www.gnu.org/software/epsilon
* My personal web site:  http://ageinghacker.net

I support everyone's freedom of mocking any opinion or belief, no
matter how deeply held, with open disrespect and the same unrelented
enthusiasm of a toddler who has just learned the word "poo".

/* Written by Luca Saiu  http://ageinghacker.net
   I, the author, hereby place this code into the public domain
   up to the extent of the applicable law. */

#define _GNU_SOURCE
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>
#include <fenv.h>
#include <limits.h>
#include <assert.h>
#include <signal.h>

#define FTYPE_IS_FLOAT

#ifdef FTYPE_IS_FLOAT
# define FTYPE float
# define UTYPE uint32_t
#else
# define FTYPE double
# define UTYPE uint64_t
#endif

#define TAG_BIT_NO 4
#define TAGMASK    ((UTYPE) (1 << TAG_BIT_NO) - 1)
#define TAG        TAGMASK // All tag bits set to 1.

union conv
{
  FTYPE f;
  UTYPE u;
};

void
print_bits_u (UTYPE a)
{
  int i;
  for (i = sizeof (UTYPE) * CHAR_BIT - 1; i >= 0; i --)
    {
      UTYPE mask = 1LLU << i;
      if (a & mask)
        printf ("1");
      else
        printf ("0");
    }
}

void
print_bits_f (FTYPE a)
{
  union conv c0;
  c0.f = a;
  print_bits_u (c0.u);
}

__attribute__ ((const))
FTYPE
untag__ (UTYPE a)
{
  union conv c;
  c.u = a;
  /* This is correct. */
  //c.u &= ~ ((1 << TAG_BIT_NO) - 1);

  /* This alternative is equally correct, assuming that the argument is actually
     tagged, as it should be, and faster on MIPS and SH.

     On MIPS the immediate variant of the bitwise-anding instruction takes a
     zero-extended, rather than sign-extended, argument; therefore it cannot be
     used here where the mask would need sign-extension to fit in a 16-bit
     immediate.  GCC solves this the obvious way, by loading a literal into a
     temporary with one separate instruction before and-ing two registers into
     another register.  On the other hand, any sensible architecture, MIPS
     included, has a sum or subtract instruction taking a sign-extended literal.

     SH is even more restricted: it allows and-ing an immediate (again,
     zero-extended) only when the operand (source and destination) is r0.
     Replacing the and with an add (sign-extended immediate, no register
     restriction) solves the problem. */
  c.u -= TAG;
  return c.f;
}

/* This is correct in Standard C99, relying on compound literals, designated
   initializers and non-constant initializers, but not on GNU C expression
   statements. */

#define f2u(a)  \
 (((union conv) {.f = (a)}).u)
#define u2f(a)  \
 (((union conv) {.u = (a)}).f)

#define untag_alt(a)      \
  (((union conv)          \
    {.u = (a) - TAG}).f)

#define untag_stmtexp(a)  \
  ({                      \
    union conv _c;        \
    _c.f = (a);           \
    _c.u &= ~ TAGMASK;    \
    _c.f;                 \
  })

#define untag(a) \
  (u2f (a - TAG))
#define tag(a) \
  ((f2u (a) & ~ TAGMASK) | TAG)

__attribute__ ((const))
UTYPE
tag_ (FTYPE a)
{
  return tag (a);
}

__attribute__ ((const))
FTYPE
untag_ (UTYPE a)
{
  return untag (a);
}

void print_float (FTYPE a)
{
  printf ("%-26.15f ", a);
  print_bits_f (a);
  //feclearexcept (FE_ALL_EXCEPT);
  if (! iscanonical (a))
    printf (" NC");
  if (! isfinite (a) && ! isnan (a))
    printf (" I");
  //if (isnormal (a))
  //  printf (" n");
  if (issignaling (a))
    printf (" S");
  if (issubnormal (a))
    printf (" s");
  if (isnan (a))
    printf (" N");
  if (iszero (a))
    printf (" Z");
fflush (stdout);
}

double
worst_error = 0.0;

__attribute__ ((noinline, noclone))
void
test (FTYPE a)
{
  // assert (untag (tag (a)) == untag (a));
  FTYPE t = untag (tag (a));
  printf ("*   "); print_float (a); printf ("\n");
  printf ("  ->"); print_float (t); printf ("\n");
  if (isfinite (a)
      && isfinite (t)
      && a != 0)
    {
      double error = fabs ((a - t) / a);
      if (error > worst_error)
        worst_error = error;
      if ((a - t) != 0)
        printf ("      (error: %.20f%%)\n", error * 100);
      else
        printf ("      (error: zero)\n");
    }
}

int
main (void)
{
#ifdef FE_SNANS_ALWAYS_SIGNAL
  printf ("FE_SNANS_ALWAYS_SIGNAL is defined.\n");
#endif // #ifdef FE_SNANS_ALWAYS_SIGNAL
  assert (sizeof (FTYPE) == sizeof (UTYPE));

  /*
  // FIXME: with POSIX threads I should use pthread_sigmask instead of 
sigprocmask.
  sigset_t signal_set;
  sigprocmask (0, NULL, & signal_set);
  printf ("SIGNAL MASK BEFORE: %lu\n", * (unsigned long *) &signal_set);
  sigaddset (& signal_set, SIGFPE);
  printf ("SIGNAL MASK AFTER:  %lu\n", * (unsigned long *) &signal_set);
  sigprocmask (SIG_BLOCK, & signal_set, NULL);
  */

  test (0.0);
  test (- 0.0);
  test (1.0);
  test (- 1.0);
  test (M_PI);
  test (- M_PI);
  test (INFINITY);
  test (- INFINITY);
  test (NAN);
  test (SNAN);
  test (1234567890.);
  test (123456789.0);
  test (12345678.90);
  test (1234567.890);
  test (123456.7890);
  test (12345.67890);
  test (1234.567890);
  test (123.4567890);
  test (12.34567890);
  test (1.234567890);
  test (.1234567890);
  test (.01234567890);
  test (.001234567890);
  test (.0001234567890);
  test (.00001234567890);
  test (.000001234567890);
  test (.0000001234567890);
  test (.00000001234567890);
  test (.000000001234567890);
  test ((double) (1LL << 23) + (1 << TAG_BIT_NO) - 1); // worst possible 
relative error on 32-bit float
  test ((double) (1LL << 52) + (1 << TAG_BIT_NO) - 1); // worst possible 
relative error on 64-bit double
  printf ("Worst error: %.20f%%\n", worst_error * 100);
  return 0;
}

signature.asc
Description: PGP signature