eccodes/tests/bits.cc

/*
 * (C) Copyright 2005- ECMWF.
 *
 * This software is licensed under the terms of the Apache Licence Version 2.0
 * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
 *
 * In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
 * virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
 */

#include <stdio.h>
#include <stdlib.h>

#include "grib_api_internal.h"

#define NO_PRINT 1
#define PRINT_ALL (1 << 1)
#define PRINT_ROUND_ERRORS (1 << 2)
#define PRINT_TRUNCATION (1 << 3)
#define NO_NEAREST_SMALLER_IBM_FLOAT (1 << 4)

void test_bits(double x, unsigned long mode)
{
    int err = 0;
    double y = 0, z = 0, w = 0, t = 0;
    double dy = 0, dz = 0, dw = 0, dt = 0;
    unsigned long i = 0, j = 0, k = 0, l = 0;

    i  = grib_ibm_to_long(x);
    y  = grib_long_to_ibm(i);
    dy = y - x;

    z  = grib_long_to_ibm(j);
    dz = z - x;

    k  = grib_ibm_nearest_smaller_to_long(x);
    w  = grib_long_to_ibm(k);
    dw = w - x;

    if (!(mode & NO_NEAREST_SMALLER_IBM_FLOAT)) {
        err  = grib_nearest_smaller_ibm_float(x, &t);
        ECCODES_ASSERT(!err);
        l  = grib_ibm_to_long(y);
        dt = t - x;
    }

    if (mode & PRINT_ALL) {
        printf("x=%.20e \n", x);
        printf("grib_ibm_to_long                  %lu (0x%lx) -> %.20e diff=%.20e\n", i, i, y, dy);
        printf("grib_nearest_smaller_ibm_float    %lu (0x%lx) -> %.20e diff=%.20e\n", l, l, t, dt);
        printf("grib_ibm_nearest_smaller_to_long  %lu (0x%lx) -> %.20e diff=%.20e\n", k, k, w, dw);
    }

    if (fabs(dz) < fabs(dy) && ((mode & PRINT_ROUND_ERRORS) | (mode & PRINT_ALL))) {
        printf("!!!!!!!!!!--------------  Rounding error \n");
        printf("x=%.20e \n", x);
        printf("grib_ibm_to_long                  %lu (0x%lx) -> %.20e diff=%.20e\n", i, i, y, dy);
        printf("grib_ibm_to_long_trunc            %lu (0x%lx) -> %.20e diff=%.20e\n", j, j, z, dz);
        printf("grib_nearest_smaller_ibm_float    %lu (0x%lx) -> %.20e diff=%.20e\n", l, l, t, dt);
        printf("grib_ibm_nearest_smaller_to_long  %lu (0x%lx) -> %.20e diff=%.20e\n", k, k, w, dw);
    }

    if ((w != t) && ((mode & PRINT_TRUNCATION) | (mode & PRINT_ALL))) {
        printf("!!!!!!!!!!---------------  Different nearest smaller \n");
        printf("x=%.20e \n", x);
        printf("grib_ibm_to_long                  %lu (0x%lx) -> %.20e diff=%.20e\n", i, i, y, dy);
        printf("grib_ibm_to_long_trunc            %lu (0x%lx) -> %.20e diff=%.20e\n", j, j, z, dz);
        printf("grib_nearest_smaller_ibm_float    %lu (0x%lx) -> %.20e diff=%.20e\n", l, l, t, dt);
        printf("grib_ibm_nearest_smaller_to_long  %lu (0x%lx) -> %.20e diff=%.20e\n", k, k, w, dw);
    }
}

void print_machine_parameters()
{
    int ibeta = 0, it = 0, irnd = 0, ngrd = 0, machep = 0, negep = 0, iexp = 0, minexp = 0, maxexp = 0;
    float eps = 0, epsneg = 0, xmin = 0, xmax = 0;

    printf("-- Machine parameters ---\n");
    printf("ibeta=%d\n", ibeta);
    printf("it=%d\n", it);
    printf("irnd=%d\n", irnd);
    printf("ngrd=%d\n", ngrd);
    printf("machep=%d\n", machep);
    printf("negep=%d\n", negep);
    printf("iexp=%d\n", iexp);
    printf("minexp=%d\n", minexp);
    printf("maxexp=%d\n", maxexp);
    printf("eps=%g\n", eps);
    printf("epsneg=%g\n", epsneg);
    printf("xmin=%.20e\n", xmin);
    printf("xmax=%.20e\n", xmax);
    printf("-------------------------\n\n");
}

int main(int argc, char* argv[])
{
    double x, y, z, d, eps, epst, epsh;
    unsigned long i, j, nc, c;
    /*unsigned char ibm[4];*/
    unsigned long iminp = 0x00100000, imaxp = 0x7fffffff;
    unsigned long iminn = 0x80100000, imaxn = 0xffffffff;
    double dminp = 0, dmaxp = 0, dminn = 0, dmaxn = 0;
    double dminpn = 0, dmaxpp = 0, dminnn = 0, dmaxnp = 0;
    unsigned long A;
    int k, cc;

    dminp  = grib_long_to_ibm(iminp);
    dminpn = grib_long_to_ibm(iminp + 1);
    dmaxp  = grib_long_to_ibm(imaxp);
    dmaxpp = grib_long_to_ibm(imaxp - 1);
    dminn  = grib_long_to_ibm(iminn);
    dminnn = grib_long_to_ibm(iminn + 1);
    dmaxn  = grib_long_to_ibm(imaxn);
    dmaxnp = grib_long_to_ibm(imaxn - 1);

    x = -1.33642013258103330000e-83;
    i = grib_ibm_to_long(x);
    j = grib_ibm_nearest_smaller_to_long(x);
    printf("grib_ibm_to_long(%.20e)=0x%lX grib_long_to_ibm(0x%lX)=%.20e grib_ibm_nearest_smaller_to_long(%.20e)=0x%lX\n",
           x, i, i, grib_long_to_ibm(i), x, j);
    printf("grib_long_to_ibm(grib_ibm_nearest_smaller_to_long(%.20e))=%.20e\n", x, grib_long_to_ibm(j));
    /* exit(0); */

    /*
  if (argc > 1 ) {
    x=atof(argv[1]);
    test_bits(x,PRINT_ALL);
    printf("-----------------\n");
    exit(0);
  }
  */

    /*
  i=0x49800000;
  x=grib_long_to_ibm(i);
  j=grib_ibm_to_long(x);
  printf(" grib_long_to_ibm(0x%lX)=%.8e  grib_ibm_to_long(grib_long_to_ibm(0x%lX))=0x%lX\n",i,x,i,j);
  exit(0);
*/

    printf("-------------------------------- ibm float -------------------------------\n");
    printf("Min positive 0x%lX   ->  %.8e ", iminp, dminp);
    printf("-- next  0x%lX   ->  %.8e  diff=%.3e\n", iminp + 1, dminpn, dminp - dminpn);
    printf("Max positive 0x%lX ->  %.8e ", imaxp, dmaxp);
    printf("-- prev  0x%lX ->  %.8e  diff= %.3e\n", imaxp - 1, dmaxpp, dmaxp - dmaxpp);
    printf("Min negative 0x%lX -> %.8e ", iminn, dminn);
    printf("-- next  0x%lX -> %.8e  diff= %.3e\n", iminn + 1, dminnn, dminn - dminnn);
    printf("Max negative 0x%lX -> %.8e ", imaxn, dmaxn);
    printf("-- prev  0x%lX -> %.8e  diff=%.3e\n", imaxn - 1, dmaxnp, dmaxn - dmaxnp);
    printf("%lu negative values, %lu positive values\n", (imaxn - iminn), (imaxp - iminp));
    printf("--------------------------------------------------------------------------\n");
    fflush(stdout);
#if 0
  if (0) {
  printf("---- Test --------- ibm_to_long(long_to_ibm()) ---------------------------\n");
  nc=((imaxn-iminp))/100;
  c=0;
  for (i=iminp;i<imaxn;i++) {
   if (c==nc) { printf("%lu \n",i-iminp);fflush(stdout); c=0; }
     if ( (i & 0x00800000) != 0 &&  i != grib_ibm_to_long(grib_long_to_ibm(i))) {
     printf("test failed: 0x%lX != ibm_to_long(long_to_ibm(0x%lX)) \n",grib_ibm_to_long(grib_long_to_ibm(i)),i);
     exit(1);
   }
   c++;
  }
  printf("%lu \n",i-iminp);
  printf("---- Test passed ---------------------------------------------------------\n");
  }

  printf("\n--- Negative Underflow test  --------------\n");

  x=1;
  y=fabs(x)/2;

  while ( grib_long_to_ibm(grib_ibm_to_long(x)) != 0 ) {
    /* test_bits(x,PRINT_ALL); */
    test_bits(x,PRINT_TRUNCATION | PRINT_ROUND_ERRORS);
    x*=y;
  }
  test_bits(x,PRINT_ALL);
  printf("----------------------------------\n");


  printf("\n--- Positive Underflow test  --------------\n");

  x=1;
  y=fabs(x)/2;
  w=0;
  while ( grib_long_to_ibm(grib_ibm_to_long(x)) != 0 ) {
    test_bits(x,PRINT_TRUNCATION | PRINT_ROUND_ERRORS);
    /* test_bits(x,PRINT_ALL); */
  w=x;
    x*=y;
  }
  test_bits(x,PRINT_ALL);
  printf("----------------------------------\n");

  printf("\n--- Positive Overflow test  --------------\n");
  x=grib_long_to_ibm(0x7fffffff);
  test_bits(x,PRINT_ALL);
  x=7.237005145973e75;
  test_bits(x,PRINT_ALL);
  x=7.237005145974e75;
  test_bits(x,PRINT_ALL);
  printf("----------------------------------\n");

  printf("\n--- Negative Overflow test  --------------\n");
  x=grib_long_to_ibm(0xffffffff);
  test_bits(x,PRINT_ALL );
  x=-7.237005145974e75;
  test_bits(x,PRINT_ALL | NO_NEAREST_SMALLER_IBM_FLOAT);
  printf("----------------------------------\n");
#endif

    printf("\n--- Rounding & nearest smaller test  --------------\n");
    fflush(stdout);
    nc = ((imaxn - iminp)) / 100;
    printf("imaxn=%lu iminp=%lu nc=%lu\n", imaxn, iminp, nc);
    fflush(stdout);
    c  = 0;
    cc = 1;
    i  = 0xCC100000;
    /* for (i=0x7F100000 ;i<imaxn;i++)  */
    for (i = iminp; i < imaxn; i++) {
        if ((i & 0x00100000) == 0 &&
            (i & 0x00200000) == 0 &&
            (i & 0x00400000) == 0 &&
            (i & 0x00800000) == 0) {
            c++;
            continue;
        }
        if (c >= nc) {
            printf("%d - %lu 0x%lX\n", cc, i - iminp, i);
            fflush(stdout);
            c = 0;
            cc++;
        }

        A   = (i & 0x7f000000) >> 24;
        /*e   = A - 70;*/
        eps = 1;
        /*m   = (i & 0xffffff);*/
        /* printf("---m=0x%lX\n",m); */

        eps = grib_ibm_table_e(A);

        epsh = eps / 2;
        epst = eps / 3;

        y = grib_long_to_ibm(i);
        if (grib_ibm_to_long(y) != i) {
            printf("------ grib_long_to_ibm(grib_ibm_to_long()) Failure -----------------\n");
            printf("i=0x%lX\n", i);
            printf("grib_long_to_ibm(0x%lX)=%.20e \n", i, grib_long_to_ibm(i));
            printf("grib_long_to_ibm(0x%lX)=%.20e \n", i - 1, grib_long_to_ibm(i - 1));
            printf("grib_long_to_ibm(0x%lX)=%.20e \n", i + 1, grib_long_to_ibm(i + 1));
            printf("eps=%.20e\n", eps);
            printf("y=%.20e grib_ibm_to_long(y)=0x%lX grib_long_to_ibm(grib_ibm_to_long(y))=%.20e\n", y, j, z);
            printf("fabs(grib_long_to_ibm(grib_ibm_to_long(y))-y)=%.20e \n", d);
            printf("--------------------------------------------------------\n");
            fflush(stdout);
            exit(1);
        }

        for (k = 1; k < 3; k++) {
            if (i == 0x7FFFFFFF) continue;
            if (y > 0)
                y += epst;
            else
                y -= epst;

            /* rounding test */
            j = grib_ibm_to_long(y);
            z = grib_long_to_ibm(j);
            d = fabs(z - y);
            if (d > epsh) {
                printf("---------------- Rounding Failure ----------------------\n");
                printf("i=0x%lX\n", i);
                printf("grib_long_to_ibm(0x%lX)=%.20e \n", i, grib_long_to_ibm(i));
                printf("grib_long_to_ibm(0x%lX)=%.20e \n", i - 1, grib_long_to_ibm(i - 1));
                printf("grib_long_to_ibm(0x%lX)=%.20e \n", i + 1, grib_long_to_ibm(i + 1));
                printf("eps=%.20e\n", eps);
                printf("y=%.20e grib_ibm_to_long(y)=0x%lX grib_long_to_ibm(grib_ibm_to_long(y))=%.20e\n", y, j, z);
                printf("fabs(grib_long_to_ibm(grib_ibm_to_long(y))-y)=%.20e \n", d);
                printf("--------------------------------------------------------\n");
                fflush(stdout);
                exit(1);
            }

            /* nearest smaller test */
            j = grib_ibm_nearest_smaller_to_long(y);
            z = grib_long_to_ibm(j);
            d = z - y;
            if (d > 0) {
                printf("---------------- Nearest smaller Failure ---------------\n");
                printf("i=0x%lX\n", i);
                printf("grib_long_to_ibm(0x%lX)=%.20e \n", i - 1, grib_long_to_ibm(i - 1));
                printf("grib_long_to_ibm(0x%lX)=%.20e \n", i, grib_long_to_ibm(i));
                printf("grib_long_to_ibm(0x%lX)=%.20e \n", i + 1, grib_long_to_ibm(i + 1));
                printf("eps=%.20e d=%.20e\n", eps, d);
                printf("                      y=%.20e grib_ibm_nearest_smaller_to_long(y)=0x%lX\n", y, j);
                printf("grib_long_to_ibm(0x%lX)=%.20e\n", j, z);
                printf("--------------------------------------------------------\n");
                fflush(stdout);
                exit(1);
            }
            if (abs(d) > eps) {
                printf("---------------- Nearest smaller warning ---------------\n");
                printf("i=0x%lX\n", i);
                printf("grib_long_to_ibm(0x%lX)=%.20e \n", i - 1, grib_long_to_ibm(i - 1));
                printf("grib_long_to_ibm(0x%lX)=%.20e \n", i, grib_long_to_ibm(i));
                printf("grib_long_to_ibm(0x%lX)=%.20e \n", i + 1, grib_long_to_ibm(i + 1));
                printf("eps=%.20e d=%.20e\n", eps, d);
                printf("                      y=%.20e grib_ibm_nearest_smaller_to_long(y)=0x%lX\n", y, j);
                printf("grib_long_to_ibm(0x%lX)=%.20e\n", j, z);
                printf("--------------------------------------------------------\n");
                fflush(stdout);
            }
        }

        c++;
    }
    printf("\n--- test finished ---------------------------------\n");

    return 0;
}