Precision-trimming expressions for cdo ยป trimprec.patch
| src/expr.cc | ||
|---|---|---|
|
static double f_deg(double x) { return x * 180.0 / M_PI; }
|
||
|
static double f_isMissval(double x) { (void)x; return 0.0; }
|
||
|
static double trimrel(double x, double keepbits) { // https://gmd.copernicus.org/articles/14/377/2021/
|
||
|
int maskbits;
|
||
|
int64_t *b;
|
||
|
int64_t mask, half_quantum1;
|
||
|
//
|
||
|
// We might wish to add something here that will check that 0 <= keepbits <= 52
|
||
|
b = (int64_t *) &x;
|
||
|
maskbits = 52 - (int) keepbits; // explicit bits in double precision
|
||
|
mask = (0xFFFFFFFFFFFFFFFF>>maskbits)<<maskbits;
|
||
|
half_quantum1 = ( 1<<(maskbits-1) ) - 1;
|
||
|
*b += ((*b>>maskbits) & 1) +half_quantum1;
|
||
|
*b &= mask;
|
||
|
return x;
|
||
|
}
|
||
|
static double trimabs(double x, double max_abs_err) {//Sec. A6 https://gmd.copernicus.org/articles/14/377/2021/
|
||
|
int maskbits;
|
||
|
double quantum;
|
||
|
int64_t *b;
|
||
|
const int64_t zerofrac = ~0ULL << 52; // zeros in place of mantissa in IEEE754 double
|
||
|
const double maxval = double(1ULL << 52); // max quantizable value in quanta
|
||
|
//
|
||
|
// We might wish to add something here that will check that max_abs_err > 0
|
||
|
//
|
||
|
if (abs(x) < maxval*max_abs_err) { //quantisation is meaningful, no integer overflow
|
||
|
b = (int64_t *) (& max_abs_err); // access it as int64
|
||
|
*b &= zerofrac; // zero fraction bits of xmax_abs_err: reduce it to benefit from compression
|
||
|
quantum = max_abs_err * 2; //quantum -- twice of that
|
||
|
x = nearbyint (x / quantum) * quantum;
|
||
|
}
|
||
|
return x;
|
||
|
}
|
||
|
static double pt_definegrid(const ParamEntry *const p) { printf("pt_definegrid: >%s<\n", p->name.c_str()); return -1;/*return cdo_define_grid(xxx);*/ }
|
||
|
static double pt_grid(const ParamEntry *const p) { return p->gridID; }
|
||
|
static double pt_ngp(const ParamEntry *const p) { return p->ngp; }
|
||
| ... | ... | |
|
{ FT_STD, 0, "pow", (void (*)(void))((double (*)(double, double)) std::pow) },
|
||
|
{ FT_STD, 0, "hypot", (void (*)(void))((double (*)(double, double)) std::hypot) },
|
||
|
{ FT_STD, 0, "atan2", (void (*)(void))((double (*)(double, double)) std::atan2) },
|
||
|
{ FT_STD, 0, "trimrel", (void (*)(void))((double (*)(double, double)) (&trimrel)) },
|
||
|
{ FT_STD, 0, "trimabs", (void (*)(void))((double (*)(double, double)) (&trimabs)) },
|
||
|
// scalar functions y=func(x)
|
||
|
{ FT_STD, 0, "float", reinterpret_cast<void (*)(void)>(&f_float) },
|
||
|
{ FT_STD, 0, "int", reinterpret_cast<void (*)(void)>(&f_int) },
|
||
| src/operator_help.cc | ||
|---|---|---|
|
" pow(x,y) " " Power function",
|
||
|
" hypot(x,y) " " Euclidean distance function, sqrt(x*x + y*y)",
|
||
|
" atan2(x,y) " " Arc tangent function of y/x, using signs to determine quadrants",
|
||
|
" trimrel(x,kb)" " Trim relative precision to kb keep-bits. Max relative error to a value",
|
||
|
" " " | result - x | / x < 2^(-1-kb) for any finite x. ",
|
||
|
" " " trimrel is a non-decreasing funciorn of x ",
|
||
|
" " " Loosely follows A5 of https://gmd.copernicus.org/articles/14/377/2021/",
|
||
|
" trimabs(x,err)" " Trim absolute precision introducing given max. absolute error",
|
||
|
" " " | result - x | < err, trimabs is a non-decreasing funciorn of x",
|
||
|
" " " Loosely follows A6 of https://gmd.copernicus.org/articles/14/377/2021/",
|
||
|
" ",
|
||
|
" Coordinates:",
|
||
|
" ",
|
||