/* Last edited on 2024-11-19 05:48:20 by stolfi */ /* ---------------------------------------------------------------------- */ /* wt_median.c */ double wt_median(int32_t nx, double x[], int32_t ix, int32_t nw, int32_t w[], int32_t wsum, bool_t interp, int32_t nk, int32_t kx[]) { bool_t debug = FALSE; demand(wsum > 0, "total weight must be positive"); demand(wsum < INT32_MAX/2, "total weight is too large"); demand((nw >= 1) && ((nw % 2) == 1), "invalid window width {nw}"); int32_t hw = nw/2; demand((ix-hw >= 0) && (ix+hw < nx), "window spills outside {0..nx-1}"); if (nw == 1) { /* Trivial median: */ if (kx != NULL) { kx[0] = ix; } return x[ix]; } /* Get the set {tx[0..nw-1]} of indices of elems in window, sorted by {x} value: */ int32_t *tx = NULL; { int32_t nt; if (kx == NULL) { nt = 0; tx = talloc(nw, int32_t); } else { nt = nk; tx = kx; } int32_t nt_add = wt_median_index_set_update(nt, tx, nx, ix, nw); if (nt_add < 5) { wt_median_index_set_insertion_sort(nx, x, nw, tx); } else { wt_median_index_set_full_sort(nx, x, nw, tx); } } int32_t ia, ib; int32_t Fa, Fb; wt_median_find_crossing(nx, x, ix, nw, w, wsum, tx, &ia, &Fa, &ib, &Fb); if (debug) { fprintf(stderr, " ia = %d Fa = %+12d ib = %d Fb = %+12d", ia, Fa, ib, Fb); } double xm = NAN; if (ia >= ib) { assert((Fa == 0) && (Fb == 0)); /* Because {F} is monotonic non-decreasing. */ /* Any sample in {x[tx[ib..ia]]} is a median: */ int32_t im = (ia + ib)/2; if (im < 0) { im = 0; } else if (im >= nw) { im = nw-1; } else { im = (ia + ib + (im % 2))/2; /* Round to even if {ia+ib} is odd. */ } xm = x[tx[im]]; } else { assert(ib == ia+1); /* Because {F} is monotonic non-decreasing. */ double xa = x[tx[ia]]; double xb = x[tx[ib]]; if (xa == xb) { /* Imperfect median, but the value is {xa==xb} anyway: */ xm = xa; } else { /* Median is between {xa} and {xb}: */ assert(xa < xb); if (interp) { double f = ((double)Fb)/((double)(Fb - Fa)); xm = f*xa + (1-f)*xb; } else { /* Median is either at {ib} or at {ia}: */ if (fabs(Fa) < fabs(Fb)) { xm = xa; } else if (fabs(Fb) < fabs(Fa)) { xm = xb; } else { /* Choose the one with even window index: */ xm = ((ib % 2) == 0 ? xb : xa); } } } } /* Cleanup: */ if (tx != kx) { free(tx); } return xm; } int32_t wt_median_index_set_update(int32_t nt, int32_t tx[], int32_t nx, int32_t ix, int32_t nw); /* Assumes that, on entry, {tx[0..nt-1]} is a set of {nt} consecutive indices in {0..nx-1}. On exit, {tx[0..nw-1]} will contain the {nw} consecutive indices centered at {ix} namely {ix-hw..ix+hw} where {hw=nw/2}. The input set size {nt} may be bigger or smaller than the window width {nw}, which must be odd. In any case, the indices of the output set which are in the input set are merely copied, preserving their order. The procedure returns the number of new indices, which will be appended at the end. The range {ix-hw..ix+hw} must be a sub-range of {0..nx-1}. */ void wt_median_find_crossing ( int32_t nx, double x[], int32_t ix, int32_t nw, int32_t w[], int32_t wsum, int32_t kx[], int32_t *ia_P, int32_t *Fa_P, int32_t *ib_P, int32_t *Fb_P ); /* Assumes that {x[0..nx-1]} is a list of samples and {ix} is the index of the central sample of a window of {nw} consecutive samples, which have weights {w[0..nw-1]] with total sum {wsum}. The window width {nw} must be odd. Also assumes that {kx[0..nw-1]} is a permutation of the indices {ix-hw..ix+hw}, where {hw=nw/2}, and that the samples {x[kx[k]]} are non-decreasing. Finds the largest index {ia} in {0..nw-1} such that {F(x[kx[ia]])} is zero or negative, and the smallest index {ib} in {0..nw-1} such that {F(x[kx[ib]])} is zero or positive. Returns these indices in {*ia_P,*ib_P} and those two {F} values in {*Fa_P,*Fb_P}. */ int32_t wa = w[hw + tx[ib]-ix]; int32_t wb = w[hw + tx[ib]-ix]; int32_t da = abs(Fa + wa); /* What {|F|} would be if {xa} is chosen. */ int32_t db = abs(Fb - wb); /* What {|F|} would be if {xb} is chosen. */ assert(db + da > 0); /* ---------------------------------------------------------------------- */ /* wt_median.c (binary search version) */ /* Find min and max values in window: */ double xlo = +INF; double xhi = -INF; double sumw = 0; for (int32_t k = -hw; k <= +hw; k++) { double xk = x[ix+k]; demand(isfinite(xk), "can't handle {±INF} or {NAN}"); if (xk < xlo) { xlo = xk; } if (xk > xhi) { xhi = xk; } sumw += w[hw+k]; } double halfw = 0.5000000001*sumw; /* Cooked to hopefully handle roundoff. */ double xm = (isfinite(guess) && (guess >= xlo) && (guess <= xhi) ? guess : (xlo + xhi)/2); int32_t ncand = nw; /* Number of window elems in {[xlo _ xhi]}. */ while (TRUE) { if (debug) { fprintf(stderr, " xlo = %+14.8f xhi = %+14.8f ncand = %d xm = %+14.8f\n", xlo, xhi, ncand, xm); } double sumlo = 0, sumhi = 0; /* Sum of weights below and above {xm}. */ int32_t ixlo = -1, ixhi = -1; /* Extreme elems below and above {xm}. */ for (int32_t k = -hw; k <= +hw; k++) { double xk = x[ix + k]; if (xk < xm) { if ((ixlo == -1) || (xk > x[ixlo])) { ixlo = ix+k; } sumlo += w[hw + k]; } else if (xk > xm) { if ((ixhi == -1) || (xk < x[ixhi])) { ixhi = ix+k; } sumhi += w[hw + k]; } } if (debug) { fprintf(stderr, " sumlo = %12.8f ixlo = %3d sumhi = %12.8f ixhi = %3d\n", sumlo, ixlo, sumhi, ixhi); } if ((sumlo <= halfw) && (sumhi <= halfw)) { return xm; } assert(ncand > 0); /* Shrink search interval {xlo,xhi}: */ if (sumlo > halfw) { assert(ixlo != -1); assert(x[ixlo] < xhi); if (sumlo - w[hw+ixlo-ix] <= halfw) { return x[ixlo]; } else { xhi = x[ixlo]; } } else if (sumhi > halfw) { assert(ixhi != -1); assert(x[ixhi] > xlo); if (sumhi - w[hw+ixhi-ix] <= halfw) { return x[ixhi]; } else { xlo = x[ixhi]; } } else { assert(FALSE); } /* Recompute {ncand} and ensure it decreased: */ int32_t ncand1 = 0; for (int32_t k = -hw; k <= +hw; k++) { double xk = x[ix + k]; if ((xk >= xlo) && (xk <= xhi)) { ncand1++; } } assert(ncand1 < ncand); ncand = ncand1; /* Recompute {xm}: */ xm = (xlo + xhi)/2.0; } /* ---------------------------------------------------------------------- */ char *txtcat (const char *a, const char *b) { char *r = talloc(strlen(a) + strlen(b) + 1, char); affirm (r != NULL, "memory exhausted"); strcpy(r, a); strcat(r, b); return(r); } char *txtcat3 (const char *a, const char *b, const char *c) { char *r = talloc(strlen(a) + strlen(b) + strlen(c) + 1, char); affirm (r != NULL, "memory exhausted"); strcpy(r, a); strcat(r, b); strcat(r, c); return(r); } char *fmt_int(int64_t x, uint32_t wid) { #define fmt_int_BUFSIZE 1000 demand(wid <= fmt_int_BUFSIZE-1, "{wid} too big"); char buf[fmt_int_BUFSIZE]; int32_t rcode = snprintf(buf, fmt_int_BUFSIZE, "%0*ld", wid, x); assert(rcode >= 0); { uint64_t n = strlen(buf); char *r = talloc(n + 1, char); strcpy(r, buf); return r; } #undef fmt_int_BUFSIZE } /* ---------------------------------------------------------------------- * /* fget.c */ char *fget_to_delim(FILE *f, char delim) { return fget_to_delims(f, delim, NULL); } bool_t fget_test_comment_or_eol(FILE *f, char cmtc); /* Similar to {fget_comment_or_eol}, but returns {FALSE} if it finds anything (including end-of-file) before a newline or {cmtc}. In this case, that character is not consumed. Otherwise the procedure returns {TRUE}, after consuming all characteds up to and including the newline. The character {cmtc} should not be a space or {EOF}. */ bool_t fget_test_comment_or_eol(FILE *f, char cmtc) { fget_skip_spaces(f); int32_t r = fgetc(f); if (r == EOF) { return FALSE; } else if ((char)r == '\n') { return TRUE; } else if ((char)r == cmtc) { do { r = fgetc(f); } while ((r != EOF) && ((char)r != '\n')); demand(r != EOF, "no newline at end of file"); return TRUE; } else { if (r != EOF) { ungetc(r, f); } return FALSE; } } char fget_skip_spaces_to_something(FILE *f); /* Skips spaces, then requires that the next character exists but is not a formatting char. Consumes that character and returns it. */ int32_t cmax_dec, cmax_alo, cmax_ahi; /* Max chars for decmimal and alpha ranges. */ if (base <= 10) { cmax_dec = '0' + (base-1); /* Not used: */ cmax_alo = cmax_ahi = '*'; } else { cmax_dec = '9'; cmax_alo = 'a' + (base-11); cmax_ahi = 'A' + (base-11); } int32_t fget_digit(FILE *f, uint32_t base) { int32_t r = fgetc(f); if (r == EOF) { return -1; } char c = (char)r; if ((c >= '0') && (c <= '9')) { return c - '0'; } else if (alpha) { if ((c >= 'a') && (c <= 'z')) { return 10 + (c - 'a'); } else if ((c >= 'A') && (c <= 'Z')) { return 10 + (c - 'A'); } } else { ungetc(r, f); } return -1; } } int fget_int(FILE *f); unsigned int fget_uint(FILE *f, int base); int fget_int(FILE *f) { int64_t x = fget_int64(f); demand((x >= INT_MIN) && (x < INT_MAX), "integer does not fit in {int} type"); return (int)x; } unsigned int fget_uint(FILE *f, int32_t base) { uint64_t x = fget_uint64(f, base); demand(x < UINT_MAX, "integer does not fit in {unsigned int32_t} type"); return (unsigned int)x; } int nget_int(FILE *f, char *name); unsigned int nget_uint(FILE *f, char *name, int base); int nget_int(FILE *f, char *name) { nget_name_eq(f, name); return fget_int(f); } unsigned int nget_uint(FILE *f, char *name, int base) { nget_name_eq(f, name); return fget_uint(f, base); } /* ---------------------------------------------------------------------- */ /* hufftree.{h,c} */ /* THE PARENT ENCODING TABLE The procedure also returns in the vector {parent} a table that lets one encode the elements of {V} into variable-length binary strings, that can be decoded with the tree {T}. For this table, each node is represented by a /node index/ in {0..(maxval+1)+ni-1} where {ni} is the number of internal nodes created (which is at most {(maxval+1)-1}). The index of a leaf node {p} is {p.value}. The index of an internal node {p} is {(maxval+1)} plus the order in which it was created, namely {(maxval+1) + (p.value - MIN_VALUE)}. The root node will thus have index {(maxval+1)+ni-1}. If a node {q} (leaf or internal) with index {iq} is the child {p.child[r]} of an internal node {p} with index {ip}, then {parent[iq]} will be {2*ip + r}. Note that {ip} must be strictly greater than {iq}, so {parent[iq]} must be strictly greater than {iq} if {q} has a parent. Note also that the maximum index of a node other than the root is {2*MAX_VALUE-1}, so the maximum value stored in {parent} is {4*MAX_VALUE-1} with by definition is less than {UINT32_MAX}. If {q} has no parent, namely it is the root node, then {parent[iq]} is set to zero. Thus the bit string that encodes any value {val} in {V} can be obtained by starting with and {iv=val} and repeatedly computing {r = parent[iv]%2, iv = parent[iv]/2} until reaching the root. The bits {r} thus computed, in reverse order, will be the code of {val}. The entry {parent[val]} will be zero also for any value {val} in {0..(maxval+1)-1} that was excluded from the tree for having zero frequency. These values can be detected because they are in {0..(maxval+1)-1} whereas the true root node has index at least {(maxval+1)}. */ /* ---------------------------------------------------------------------- */ /* mst.h */ int mst_find_root(int R[], int u); /* The "find" operation of Union-Find: Finds the root of the tree that contains the vertex {u}, by following the shortcut map {R} until reaching a vertex {v} such that {R[v]=v}. In the process, changes {R[w]} of every visited vertex {w} to point directy to {u}, thus making {u} the root of the tree. */ /* ---------------------------------------------------------------------- */ /* spmat.c */ spmat_t *spmat_cast_ref(void *M); /* Casts the argument as an address to a {spmat_t}. Can be used to apply spmat_expand or spmat_trim to sparse matrices of specific types. We define it as a function in order to trigger type warnings when {M} is not an address (a common mistake). */ /* ---------------------------------------------------------------------- */ /* ulist.c */ /* Picking a random non-NULL item in a sequential hash table. */ ulist_item_t ulist_pick(ulist_t *S) { ulist_count_t sz = S->it.nel; if (sz == 0) { return NULL; } ulist_index_t k = abrandom(0, sz-1); if (S->it.el[k] != NULL) { return S->it.el[k]; } if (sz == 1) { return NULL; } /* Choose a good stride {d} rel prime to {sz}: */ ulist_count_t d = (int)(0.618034 * (double)sz); while (gcd(d,sz) != 1) { d--; } ulist_count_t i = k + d % sz; do { ulist_item_t a = S->it.el[i]; if (a != NULL) { return a; } i = (i + d) % sz; } while (i != k); return NULL; /* gawk \ ' BEGIN { \ for (n = 2; n < 1000; n++) \ { d0 = int(0.618034 * n); \ d = d0; \ while (gcd(d,n) > 1) { d--; } \ printf "%5d %5d %5d %5d\n", n, d0, d, d0 - d; \ } \ } \ function gcd(a,b, r) \ { while (b != 0) \ { r = a % b; a = b; b = r; } \ return a; \ } \ ' */ } /* ---------------------------------------------------------------------- */ /* ulist.c */ void ulist_clear_all_slots(ulist_item_vec_t *it); /* Sets all slots of {it} to NULL. */ void ulist_clear_all_slots(ulist_item_vec_t *it) { unsigned int i; for (i = 0; i < it->nel; i++) { it->el[i] = NULL; } } /* ---------------------------------------------------------------------- */ /* indexing.h, indexing.c */ bool_t ix_next_C ( ix_dim_t d, ix_index_t ix[], ix_size_t sz[], ix_step_t st[], ix_pos_t *p ); /* Set {ix[0..d-1]} to the next valid index tuple in the C/Pascal standard order (with the last index in the innermost loop). More precisely, tries to increment {ix[d-1]} first. If an index {ix[i]} to be incremented has reached its limit {sz[i]-1}, resets all indices {ix[i..d-1]} to 0, and tries to increment index {ix[i-1]} instead. If this `carry over' falls off the index tuple -- that is, if the tuple {ix} changes from {sz - (1,..1)} back to {(0,..0)} -- the procedure returns TRUE; otherwise it returns FALSE. The effect is undefined if any {sz[i]} is 0 or {ix[i]} is not in {0..sz[i]-1}. If {p} is not null, assumes that {*p} is the position corresponding to the index tuple {ix}, and adjusts it to account for the change. This procedure can be used to vary the index tuple {ix[0..d-1]} over all the elements of an array, in a single loop, e.g. | if (ix_assign_min(d,ix,sz)) | { ix_pos_t p = ix_position(d,ix,bp,st); | do { ... } while (! ix_next_C(d,ix,sz,st,&p)); | } */ bool_t ix_next_C ( ix_dim_t d, ix_index_t ix[], ix_size_t sz[], ix_step_t st[], ix_pos_t *p ) { /* Look for an index {ix[i]} that can be incremented: */ int i = d-1; ix_index_t *ixi = ix + i; while ((*ixi) >= sz[i]-1) { /* Index {ix[i]} is at its limit, reset to 0: */ if (p != NULL) { (*p) -= (*ixi)*st[i]; } (*ixi) = 0; /* Wrap around? */ if (i <= 0) { return TRUE; } /* Go to next index: */ i--; ixi--; } /* Increment index {ix[i]} */ (*ixi)++; if (p != NULL) { (*p) += st[i]; } return FALSE; } bool_t ix_prev_C ( ix_dim_t d, ix_index_t ix[], ix_size_t sz[], ix_step_t st[], ix_pos_t *p ) { /* Look for an index {ix[i]} that can be decremented: */ int i = d-1; ix_index_t *ixi = ix + i; while ((*ixi) <= 0) { /* Index {ix[i]} is at its limit, reset to {sz[i]-1}: */ (*ixi) = sz[i]-1; if (p != NULL) { (*p) += (*ixi)*st[i]; } /* Wrap around? */ if (i <= 0) { return TRUE; } /* Go to next index: */ i--; ixi--; } /* Decrement index {ix[i]} */ (*ixi)--; if (p != NULL) { (*p) -= st[i]; } return FALSE; } sign_t ix_compare_C ( ix_dim_t d, ix_index_t ixa[], ix_index_t ixb[] ) { int i; ix_index_t *ixai = ixa, *ixbi = ixb; for (i = 0; i < d; i++) { if ((*ixai) < (*ixbi)) { return NEG; } else { return POS; } ixai++; ixbi++; } return ZER; } /* ---------------------------------------------------------------------- */ /* jstime.c */ #if (defined(SunOS4)) #warning "implementing {now} with {getrusage(...)}." double now(void) { struct rusage ru; getrusage(RUSAGE_SELF, &ru); return(((double)ru.ru_utime.tv_sec)*1000000.0 + ((double)ru.ru_utime.tv_usec)); } #endif #if (defined(SunOS5)) #warning "implementing {now} with {times(...)}." extern long int __sysconf (int); double now(void) { struct tms buf; clock_t t; times(&buf); t = buf.tms_utime; return(1000000.0L * ((double) t)/((double) _sysconf(3))); } #endif #if (defined(OSF1V4)) #warning "implementing {now} with {times(...)}." double now(void) { struct tms buf; clock_t t; times(&buf); t = buf.tms_utime; return(1000000.0L * ((double) t)/((double) CLK_TCK)); } #endif /* ---------------------------------------------------------------------- */ /* Locate {jMin,jMax} near {iTry} s.t. */ /* {iMin <= jMin <= jMax <= iMax} and {f(jMin) < 0 <= f(jMax)}: */ int jMin = iTry-1; int jMax = iTry; int d = 1; if (TB_DEBUG) { fprintf(stderr, "!<%d:%d>", jMin,jMax); } double gMin = f(jMin); while ((jMin > iMin) && (gMin >= 0)) { jMax = jMin; jMin -= d; if (jMin < iMin) { jMin = iMin; } gMin = (jMin == iMin ? fMin : f(jMin)); d += d; if (TB_DEBUG) fprintf(stderr, "-"); if (TB_DEBUG) fprintf(stderr, "!<%d:%d>", jMin,jMax); } double gMax = f(jMax); while ((jMax < iMax) && (gMax < 0)) { jMin = jMax; jMax += d; if (jMax > iMax) { jMax = iMax; } gMax = (jMax == iMax ? fMax : f(jMax)); d += d; if (TB_DEBUG) fprintf(stderr, "+"); if (TB_DEBUG) fprintf(stderr, "!<%d:%d>", jMin,jMax); } if (jMin != iMin) { iMin = jMin; fMin = f(iMin); } if (jMax != iMax) { iMax = jMax; fMax = f(iMax); } demand (fMin <= fMax, "tbl out of order"); /* ---------------------------------------------------------------------- */ void srt_bins_sort(int *h, int n, int cmp(int x, int y), int sgn) { int m; for (m = 1; m < n; m++) { /* Insert {h[m]} among {h[0..m-1]}: */ int v = h[m], i = m; while ((i > 0) && (sgn*cmp(v, h[i-1]) < 0)) { h[i] = h[i-1]; i--; } h[i] = v; } } /* ---------------------------------------------------------------------- */ void ix_indices ( ix_dim_t d, ix_disp_t pos, ix_disp_t b, ix_disp_t s[], ix_index_t e[] ); /* Computes the indices {e[0,..d-1]} of an element given its position {pos}. Assumes that the steps {s[0..i-1]} are graded (see {ix_step_valid} below). May fail if {pos} is not a valid position computable from {b} and {s}. */ /* ---------------------------------------------------------------------- */ /* INTERNAL PROTOS */ void ix_progerror ( const char *msg, const char *file, const unsigned int line, const char* proc ); /* Prints {file ":" line ": (" *proc ")" *msg} to {stderr} and stops. */ #define ix_error(msg) \ ix_progerror((msg), __FILE__, __LINE__, __FUNCTION__) /* Prints {*msg} and the source location. */ /* ---------------------------------------------------------------------- */ /* ERROR MESSAGES */ void ix_progerror ( const char *msg, const char *file, const unsigned int line, const char* proc ) { fprintf (stderr, "%s:%u: *** (%s) %s\n", file, line, proc, msg); exit(1); } /* ---------------------------------------------------------------------- */ /* INDEX ROLL OPERATION - CORRECT BUT TOO CONFUSING TO BE USEFUL */ void ix_roll ( ix_dim_t d, ix_size_t sz[], ix_pos_t *bp, ix_step_t st[], ix_axis_t i, ix_axis_t j, int n ); /* Rotates indices {i..j} by {n} places forward. The procedure moves {sz[k]} to {sz[k']} and {st[k]} to {st[k']} for all {k} in {i..j} simultaneously, where {k'} is {k + n} reduced cyclically to the range {i..j} --- that is, {k'=i+(k-i+n)mod(j-i+1)}. Thus, {ix[k] = ix'[k']} for all {k} in {i..j}, and {ix[k]=ix'[k]} for all other {k}. The displacement {n} may be negative. For example, to eliminate a trivial index {i} from a matrix, use {ix_roll(d,sz,&bp,st,i,d-1,-1); d--;}. The procedure has no effect if {n==0} or {i==j}. Undefined if {i > j}. */ void ix_roll ( ix_dim_t d, ix_size_t sz[], ix_pos_t *bp, ix_step_t st[], ix_axis_t i, ix_axis_t j, int n ) { /* Check axes: */ demand(i < d, "bad axis i"); demand(j < d, "bad axis j"); demand(i <= j, "bad axis interval"); /* Check for trivial roll: */ if ((i == j) || (n == 0)) { return; } /* Number of elements to be rolled: */ int r = j - i + 1; /* Reduce {n} modulo {r}: */ n = n % r; if (n < 0) { n += r; } /* Remaining candidate cycle heads are {i+kini..i+klim-1}: */ int kini = 0, klim = r; while (kini < klim) { /* Save cycle head {i+kini}: */ ix_size_t tm = sz[i + kini]; ix_pos_t ts = st[i + kini]; /* Shift elems in cycle starting at {i+kini}: */ int ka = kini; int kb; do { kb = (ka + r - n) % r; sz[i + ka] = sz[i + kb]; st[i + ka] = st[i + kb]; ka = kb; /* Update {klim} (should happen in first cycle only): */ if (kb < klim) { klim = kb; } } while (ka != kini); /* Store cycle head: */ sz[i + kb] = tm; st[i + kb] = ts; /* Get next cycle head: */ kini++; } /* Exchange parameters of axes {i,j}: */ } /* ---------------------------------------------------------------------- */ /* box.c */ bool_t ety_BLO = (x <= LO(*Ba)); bool_t ety_BHI = (x >= HI(*Ba)); bool_t ety_BMD; if (LO(*Ba) == HI(*Ba)) { ety_BMD = (x != LO(*Ba)); } else { ety_BMD = (x <= LO(*Ba)) || (x >= HI(*Ba)); } if (debug) { fprintf(stderr, " ety BLO = %d BMD = %d BHI = %d\n", ety_BLO, ety_BMD, ety_BHI); } /* Return the parts: */ if (BLO != NULL) { if (ety_BLO) { box_empty(d, BLO); } else { for (int32_t i = 0; i < d; i++) { if (i == a) { LO(BLO[i]) = LO(B[i]); HI(BLO[i]) = x; assert(LO(BLO[i]) < HI(BLO[i])); } else { BLO[i] = B[i]; } } } } if (BMD != NULL) { if (ety_BMD) { box_empty(d, BMD); } else { for (int32_t i = 0; i < d; i++) { if (i == a) { LO(BMD[i]) = x; HI(BMD[i]) = x; } else { BMD[i] = B[i]; } } } } if (BHI != NULL) { if (ety_BHI) { box_empty(d, BHI); } else { for (int32_t i = 0; i < d; i++) if (i == a) { LO(BHI[i]) = x; HI(BHI[i]) = HI(B[i]); assert(LO(BHI[i]) < HI(BHI[i])); } else { BHI[i] = B[i]; } } } } /* ---------------------------------------------------------------------- */ /* interp_spline.{h,c} */ interp_spline_kind_O, /* O-splines (interpolating, min degree, not convex). */ #include case interp_spline_kind_O: return interp_spline_O_compute_num_samples(ord); case interp_spline_kind_O: interp_spline_O_get_weights(z, ord, nw, wt); break; /* ---------------------------------------------------------------------- */ /* wt_table_convolution.c */ bool_t debug = TRUE; if (debug) { fprintf(stderr, " i = %u k2min = %u k2max = %u", i, k2min, k2max); } if (debug) { fprintf(stderr, " k1min = %d k1max = %d\n", ((int32_t)i)-(int32_t)(k2max*stride), ((int32_t)i)-(int32_t)(k2min*stride)); }