core
/
eccodes
mirror of https://github.com/ecmwf/eccodes.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity
eccodes/tests/bits.c

327 lines
12 KiB
C
Raw Blame History

/*
* (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);
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;
}
Reference in New Issue View Git Blame Copy Permalink