void do_boundaries_fcs(int nloc, fcs_float *pos)
{
  int l, i;

  /* PBC in x direction */
  if (1==pbc_dirs.x)
    for (l=0; l<3*nloc; l+=3) {
      i = -FLOOR( SPRODFCS(pos+l,tbox_x) );
      pos[l  ] += i * box_x.x;
      pos[l+1] += i * box_x.y;
      pos[l+2] += i * box_x.z;
    }

  /* PBC in y direction */
  if (1==pbc_dirs.y)
    for (l=0; l<3*nloc; l+=3) {
      i = -FLOOR( SPRODFCS(pos+l,tbox_y) );
      pos[l  ] += i * box_y.x;
      pos[l+1] += i * box_y.y;
      pos[l+2] += i * box_y.z;
    }

  /* PBC in z direction */
  if (1==pbc_dirs.z)
    for (l=0; l<3*nloc; l+=3) {
      i = -FLOOR( SPRODFCS(pos+l,tbox_z) );
      pos[l  ] += i * box_z.x;
      pos[l+1] += i * box_z.y;
      pos[l+2] += i * box_z.z;
    }

}

void
aubio_pitchmcomb_do (aubio_pitchmcomb_t * p, cvec_t * fftgrain, fvec_t * output)
{
  uint_t j;
  smpl_t instfreq;
  fvec_t *newmag = (fvec_t *) p->newmag;
  //smpl_t hfc; //fe=instfreq(theta1,theta,ops); //theta1=theta;
  /* copy incoming grain to newmag */
  for (j = 0; j < newmag->length; j++)
    newmag->data[j] = fftgrain->norm[j];
  /* detect only if local energy > 10. */
  //if (aubio_level_lin (newmag) * newmag->length > 10.) {
  //hfc = fvec_local_hfc(newmag); //not used
  aubio_pitchmcomb_spectral_pp (p, newmag);
  aubio_pitchmcomb_combdet (p, newmag);
  //aubio_pitchmcomb_sort_cand_freq(p->candidates,p->ncand);
  //return p->candidates[p->goodcandidate]->ebin;
  j = (uint_t) FLOOR (p->candidates[p->goodcandidate]->ebin + .5);
  instfreq = aubio_unwrap2pi (fftgrain->phas[j]
      - p->theta->data[j] - j * p->phasediff);
  instfreq *= p->phasefreq;
  /* store phase for next run */
  for (j = 0; j < p->theta->length; j++) {
    p->theta->data[j] = fftgrain->phas[j];
  }
  //return p->candidates[p->goodcandidate]->ebin;
  output->data[0] =
      FLOOR (p->candidates[p->goodcandidate]->ebin + .5) + instfreq;
  /*} else {
     return -1.;
     } */
}

/**
 *@brief fold space to extend the domain over which we can take
 * noise values.
 *
 *@n x, y, z are set to the noise space location for the source point.
 *@n ix, iy, iz are the integer lattice point (integer portion of x, y, z)
 *@n fx, fy, fz are the fractional lattice distance above ix, iy, iz
 *
 * The noise function has a finite domain, which can be exceeded when
 * using fractal textures with very high frequencies.  This routine is
 * designed to extend the effective domain of the function, albeit by
 * introducing periodicity.	-FKM 4/93
 */
static void
filter_args(fastf_t *src, fastf_t *p, fastf_t *f, int *ip)
{
    register int i;
    point_t dst;
    static unsigned long max2x = ~((unsigned long)0);
    static unsigned long max = (~((unsigned long)0)) >> 1;

    for (i=0; i < 3; i++) {
	/* assure values are positive */
	if (src[i] < 0) dst[i] = -src[i];
	else dst[i] = src[i];


	/* fold space */
	while (dst[i] > max || dst[i]<0) {
	    if (dst[i] > max) {
		dst[i] = max2x - dst[i];
	    } else {
		dst[i] = -dst[i];
	    }
	}

    }

    p[X] = dst[0];	ip[X] = FLOOR(p[X]);	f[X] = p[X] - ip[X];
    p[Y] = dst[1];	ip[Y] = FLOOR(p[Y]);	f[Y] = p[Y] - ip[Y];
    p[Z] = dst[2];	ip[Z] = FLOOR(p[Z]);	f[Z] = p[Z] - ip[Z];
}

static int
_fill_polygon_gray8(uint8_t *img,
        int npoints, Edge *e, uint8_t intensity, 
        int w, int h, int rowstride)
{
    int i, j;
    float *xx;
    int ymin, ymax;
    float y;

    if (npoints <= 0) return 0;

    /* Find upper and lower polygon boundary (within image) */
    ymin = e[0].ymin;
    ymax = e[0].ymax;
    for (i = 1; i < npoints; i++) {
        if (e[i].ymin < ymin) ymin = e[i].ymin;
        if (e[i].ymax > ymax) ymax = e[i].ymax;
    }

    if (ymin < 0) ymin = 0;
    if (ymax >= h) ymax = h-1;

    /* Process polygon edges */
    xx = malloc(npoints * sizeof(float));
    if (!xx) return -1;

    for (;ymin <= ymax; ymin++) {
        y = ymin+0.5F;
        for (i = j = 0; i < npoints; i++) {
            if (y >= e[i].ymin && y <= e[i].ymax) {
                if (e[i].d == 0)
                    draw_hline_gray8(img, e[i].xmin, ymin, e[i].xmax, 
                            intensity, w, h, rowstride);
                else
                    xx[j++] = (y-e[i].y0) * e[i].dx + e[i].x0;
            }
        }
        if (j == 2) {
            if (xx[0] < xx[1])
                draw_hline_gray8(img, CEIL(xx[0]-0.5), ymin, FLOOR(xx[1]+0.5),
                        intensity, w, h, rowstride);
            else
                draw_hline_gray8(img, CEIL(xx[1]-0.5), ymin, FLOOR(xx[0]+0.5),
                        intensity, w, h, rowstride);
        } else {
            qsort(xx, j, sizeof(float), x_cmp);
            for (i = 0; i < j-1 ; i += 2)
                draw_hline_gray8(img, CEIL(xx[i]-0.5), ymin, FLOOR(xx[i+1]+0.5),
                        intensity, w, h, rowstride);
        }
    }

    free(xx);

    return 0;
}

void JelloMesh::GetCell(int idx, int& i, int &j, int& k) const
{
    float rows = m_rows+1;
    float cols = m_cols+1;
    float stacks = m_stacks+1;

    // derived from idx = cols*(rows*k + i) + j
    float tmp = FLOOR(idx/cols);
    j = (int) ROUND(cols*(FRACT(idx/cols)));
    i = (int) ROUND(rows*(FRACT(tmp/rows)));
    k = (int) FLOOR(tmp/rows);
}

double fprec(double x, double digits)
{
    double l10, pow10, sgn, p10, P10;
    int e10, e2, do_round, dig;
    /* Max.expon. of 10 (=308.2547) */
    const double max10e = numeric_limits<double>::max_exponent * M_LOG10_2;

#ifdef IEEE_754
    if (ISNAN(x) || ISNAN(digits))
	return x + digits;
    if (!R_FINITE(x)) return x;
    if (!R_FINITE(digits)) {
	if(digits > 0) return x;
	else return 0;
    }
#endif
    if(x == 0) return x;
    dig = (int)FLOOR(digits+0.5);
    if (dig > MAX_DIGITS) {
	return x;
    } else if (dig < 1)
	dig = 1;

    sgn = 1.0;
    if(x < 0.0) {
	sgn = -sgn;
	x = -x;
    }
    l10 = log10(x);
    e10 = (int)(dig-1-FLOOR(l10));
    if(fabs(l10) < max10e - 2) {
	p10 = 1.0;
	if(e10 > max10e) {
	    p10 =  std::pow(10., e10-max10e);
	    e10 = static_cast<int>(max10e);
	} else if(e10 < - max10e) {
	    p10 =  std::pow(10., e10+max10e);
	    e10 = static_cast<int>(-max10e);	
	}
	pow10 = std::pow(10., e10);
	return(sgn*(FLOOR((x*pow10)*p10+0.5)/pow10)/p10);
    } else { /* -- LARGE or small -- */
	do_round = max10e - l10	 >= std::pow(10., -dig);
	e2 = dig + ((e10>0)? 1 : -1) * MAX_DIGITS;
	p10 = std::pow(10., e2);	x *= p10;
	P10 = std::pow(10., e10-e2);	x *= P10;
	/*-- p10 * P10 = 10 ^ e10 */
	if(do_round) x += 0.5;
	x = FLOOR(x) / p10;
	return(sgn*x/P10);
    }
}

/* So, all of the encodings point to the ->screen->frameBuffer,
 * We need to change this!
 */
void rfbScaledCorrection(rfbScreenInfoPtr from, rfbScreenInfoPtr to, int *x, int *y, int *w, int *h, char *function)
{
    double x1,y1,w1,h1, x2, y2, w2, h2;
    double scaleW = ((double) to->width) / ((double) from->width);
    double scaleH = ((double) to->height) / ((double) from->height);


    /*
     * rfbLog("rfbScaledCorrection(%p -> %p, %dx%d->%dx%d (%dXx%dY-%dWx%dH)\n",
     * from, to, from->width, from->height, to->width, to->height, *x, *y, *w, *h);
     */

    /* If it's the original framebuffer... */
    if (from==to) return;

    x1 = ((double) *x) * scaleW;
    y1 = ((double) *y) * scaleH;
    w1 = ((double) *w) * scaleW;
    h1 = ((double) *h) * scaleH;


    /*cast from double to int is same as "*x = floor(x1);" */
    x2 = FLOOR(x1);
    y2 = FLOOR(y1);

    /* include into W and H the jitter of scaling X and Y */
    w2 = CEIL(w1 + ( x1 - x2 ));
    h2 = CEIL(h1 + ( y1 - y2 ));

    /*
     * rfbLog("%s (%dXx%dY-%dWx%dH  ->  %fXx%fY-%fWx%fH) {%dWx%dH -> %dWx%dH}\n",
     *    function, *x, *y, *w, *h, x2, y2, w2, h2,
     *    from->width, from->height, to->width, to->height);
     */

    /* simulate ceil() without math library */
    *x = (int)x2;
    *y = (int)y2;
    *w = (int)w2;
    *h = (int)h2;

    /* Small changes for a thumbnail may be scaled to zero */
    if (*w==0) (*w)++;
    if (*h==0) (*h)++;
    /* scaling from small to big may overstep the size a bit */
    if (*x+*w > to->width)  *w=to->width - *x;
    if (*y+*h > to->height) *h=to->height - *y;
}

// Init random .............................................................
void FileName::initRandom(int length)
{
    randomize_random_generator();
    *this = "";
    for (int i = 0; i < length; i++)
        *this += 'a' + FLOOR(rnd_unif(0, 26));
}

int tabler_init(CSOUND *csound, TABL *p) {

  int ndx, len;
  int mask;

  if (UNLIKELY((p->ftp = csound->FTnp2Find(csound, p->ftable)) == NULL))
      return csound->InitError(csound,
                               Str("table: could not find ftable %d"),
                               (int) *p->ftable);
  mask = p->ftp->lenmask;
  p->np2 = mask ? 0 : 1;
  len = p->ftp->flen;

  if (*p->mode)
      p->mul = len;
    else
      p->mul = 1;

  ndx = FLOOR((*p->ndx + *p->offset)*p->mul);
  if (*p->wrap) {
    if (p->np2) {
      while(ndx >= len) ndx -= len;
      while(ndx < 0)  ndx += len;
    }
    else ndx &= mask;
  } else {
    if (UNLIKELY(ndx >= len)) ndx = len - 1;
    else if (UNLIKELY(ndx < 0)) ndx = 0;
  }
  *p->sig = p->ftp->ftable[ndx];
  return OK;
}

glm::vec3 WrappedFbmTexture::getColor(const ShadeRec& sr)const
{
	float noiseV = expansionNumber * noise->valueFBM(sr.localHitPoint);
	float value = noiseV - FLOOR(noiseV);

	return lerp(value, colorMin, colorMax);
}

int bestPrecision(float F, int _width)
{
    // If it is 0
    if (F == 0)
        return 1;

    // Otherwise
    int exp = FLOOR(log10(ABS(F)));
    int advised_prec;

    if (exp >= 0)
        if (exp > _width - 3)
            advised_prec = -1;
        else
            advised_prec = _width - 2;
    else
    {
        advised_prec = _width + (exp - 1) - 3;
        if (advised_prec <= 0)
            advised_prec = -1;
    }

    if (advised_prec < 0)
        advised_prec = -1; // Choose exponential format

    return advised_prec;
}

void aubio_hist_dyn_notnull (aubio_hist_t *s, fvec_t *input) {
  uint_t i;
  sint_t tmp = 0;
  smpl_t ilow = fvec_min(input);
  smpl_t ihig = fvec_max(input);
  smpl_t step = (ihig-ilow)/(smpl_t)(s->nelems);

  /* readapt */
  aubio_scale_set_limits (s->scaler, ilow, ihig, 0, s->nelems);

  /* recalculate centers */
  s->cent->data[0] = ilow + 0.5f * step;
  for (i=1; i < s->nelems; i++)
    s->cent->data[i] = s->cent->data[0] + i * step;

  /* scale */
  aubio_scale_do(s->scaler, input);

  /* reset data */
  fvec_zeros(s->hist);
  /* run accum */
  for (i=0;  i < input->length; i++) {
    if (input->data[i] != 0) {
      tmp = (sint_t)FLOOR(input->data[i]);
      if ((tmp >= 0) && (tmp < (sint_t)s->nelems))
        s->hist->data[tmp] += 1;
    }
  }
}

static smpl_t interp_2(fvec_t *input, smpl_t pos) {
    uint_t idx = (uint_t)FLOOR(pos);
    smpl_t frac = pos - (smpl_t)idx;
    smpl_t a = input->data[idx];
    smpl_t b = input->data[idx + 1];
    return a + frac * ( b - a );
}

/* execute tempo detection function on iput buffer */
void aubio_tempo_do(aubio_tempo_t *o, fvec_t * input, fvec_t * tempo)
{
  uint_t i;
  uint_t winlen = o->winlen;
  uint_t step   = o->step;
  fvec_t * thresholded;
  aubio_pvoc_do (o->pv, input, o->fftgrain);
  aubio_specdesc_do (o->od, o->fftgrain, o->of);
  /*if (usedoubled) {
    aubio_specdesc_do(o2,fftgrain, onset2);
    onset->data[0] *= onset2->data[0];
  }*/
  /* execute every overlap_size*step */
  if (o->blockpos == (signed)step -1 ) {
    /* check dfframe */
    aubio_beattracking_do(o->bt,o->dfframe,o->out);
    /* rotate dfframe */
    for (i = 0 ; i < winlen - step; i++ ) 
      o->dfframe->data[i] = o->dfframe->data[i+step];
    for (i = winlen - step ; i < winlen; i++ ) 
      o->dfframe->data[i] = 0.;
    o->blockpos = -1;
  }
  o->blockpos++;
  aubio_peakpicker_do (o->pp, o->of, o->onset);
  tempo->data[1] = o->onset->data[0];
  thresholded = aubio_peakpicker_get_thresholded_input(o->pp);
  o->dfframe->data[winlen - step + o->blockpos] = thresholded->data[0];
  /* end of second level loop */
  tempo->data[0] = 0; /* reset tactus */
  i=0;
  for (i = 1; i < o->out->data[0]; i++ ) {
    /* if current frame is a predicted tactus */
    if (o->blockpos == FLOOR(o->out->data[i])) {
      tempo->data[0] = o->out->data[i] - FLOOR(o->out->data[i]); /* set tactus */
      /* test for silence */
      /*
      if (aubio_silence_detection(input, o->silence)==1) {
        tempo->data[0] = 0; // unset beat if silent
      }
      */
      o->last_beat = o->total_frames + (uint_t)ROUND(tempo->data[0] * o->hop_size);
    }
  }
  o->total_frames += o->hop_size;
  return;
}

void CL_ActorStateChange (const eventRegister_t *self, struct dbuffer *msg)
{
	le_t *le;
	int entnum, state;
	character_t *chr;

	NET_ReadFormat(msg, self->formatString, &entnum, &state);

	le = LE_Get(entnum);
	if (!le)
		LE_NotFoundError(entnum);

	if (!LE_IsActor(le)) {
		Com_Printf("StateChange message ignored... LE is no actor (number: %i, state: %i, type: %i)\n",
			entnum, state, le->type);
		return;
	}

	/* If standing up or crouching down remove the reserved-state for crouching. */
	if (((state & STATE_CROUCHED) && !LE_IsCrouched(le)) ||
		 (!(state & STATE_CROUCHED) && LE_IsCrouched(le))) {
		if (CL_ActorUsableTUs(le) < TU_CROUCH && CL_ActorReservedTUs(le, RES_CROUCH) >= TU_CROUCH) {
			/* We have not enough non-reserved TUs,
			 * but some reserved for crouching/standing up.
			 * i.e. we only reset the reservation for crouching if it's the very last attempt. */
			CL_ActorReserveTUs(le, RES_CROUCH, 0); /* Reset reserved TUs (0 TUs) */
		}
	}

	/* killed by the server: no animation is played, etc. */
	if ((state & STATE_DEAD) && LE_IsLivingActor(le)) {
		le->state = state;
		FLOOR(le) = NULL;
		LE_SetThink(le, NULL);
		VectorCopy(player_dead_maxs, le->maxs);
		CL_ActorRemoveFromTeamList(le);
		return;
	} else {
		le->state = state;
		LE_SetThink(le, LET_StartIdle);
	}

	/* save those states that the actor should also carry over to other missions */
	chr = CL_ActorGetChr(le);
	if (!chr)
		return;

	chr->state = (le->state & STATE_REACTION);

	/* change reaction button state */
	if (!(le->state & STATE_REACTION)) {
		UI_ExecuteConfunc("disable_reaction");
	} else {
		UI_ExecuteConfunc("startreaction");
	}

	/* state change may have affected move length */
	CL_ActorConditionalMoveCalc(le);
}

static void quilted (const Vector3d& EPoint, const TNORMAL *Tnormal, Vector3d& normal)
{
    Vector3d value;
    DBL t;

    value[X] = EPoint[X]-FLOOR(EPoint[X])-0.5;
    value[Y] = EPoint[Y]-FLOOR(EPoint[Y])-0.5;
    value[Z] = EPoint[Z]-FLOOR(EPoint[Z])-0.5;

    t = value.length();

    t = quilt_cubic(t, dynamic_cast<QuiltedPattern*>(Tnormal->pattern.get())->Control0, dynamic_cast<QuiltedPattern*>(Tnormal->pattern.get())->Control1);

    value *= t;

    normal += (DBL)Tnormal->Amount * value;
}

//#define DEBUG
blobtype best_blob(double alpha_0, double alpha_F, double inc_alpha,
                   double a_0, double a_F, double inc_a, double w, double *target_att,
                   int target_length)
{
    blobtype retval;
    retval.order = 2;
    int alpha_size = FLOOR((alpha_F - alpha_0) / inc_alpha + 1);
    int a_size = FLOOR((a_F - a_0) / inc_a + 1);
    Image<double> att, ops;
    att().resize(alpha_size, a_size);
    ops().resize(alpha_size, a_size);
    int i, j;
    double a, alpha, best_a = -1, best_alpha = -1;
    double best_ops = 1e10, best_att = 0;
    for (i = 0, alpha = alpha_0; i < alpha_size; alpha += inc_alpha, i++)
        for (j = 0, a = a_0; j < a_size; a += inc_a, j++)
        {
            retval.radius = a;
            retval.alpha = alpha;
            A2D_ELEM(att(), i, j) = blob_att(w, retval);
            A2D_ELEM(ops(), i, j) = blob_ops(w, retval);
            if (j > 0)
                for (int n = target_length - 1; n >= 0; n--)
                    if (A2D_ELEM(att(), i, j - 1) > target_att[n] &&
                        A2D_ELEM(att(), i, j) < target_att[n])
                    {
                        A2D_ELEM(att(), i, j) = 0;
                        if (A2D_ELEM(ops(), i, j - 1) < best_ops &&
                            A2D_ELEM(att(), i, j - 1) >= best_att)
                        {
                            best_alpha = alpha;
                            best_a = a - inc_a;
                            best_ops = A2D_ELEM(ops(), i, j - 1);
                            best_att = target_att[n];
                        }
                    }
        }
#ifdef DEBUG
    att.write((std::string)"att" + floatToString(w) + ".xmp");
    ops.write((std::string)"ops" + floatToString(w) + ".xmp");
#endif

    retval.radius = best_a;
    retval.alpha = best_alpha;
    return retval;
}

ssfloat ExpMgr::gnoise(ssfloat x, ssfloat y, ssfloat z)
{
    int ix, iy, iz;
    ssfloat fx0, fx1, fy0, fy1, fz0, fz1;
    ssfloat wx, wy, wz;
    ssfloat vx0, vx1, vy0, vy1, vz0, vz1;
    static int initialized = 0;

    if (!initialized) {
        gradientTabInit(665);
        initialized = 1;
    }

    ix = FLOOR(x);
    fx0 = x - ix;
    fx1 = fx0 - 1;
    wx = SMOOTHSTEP(fx0);

    iy = FLOOR(y);
    fy0 = y - iy;
    fy1 = fy0 - 1;
    wy = SMOOTHSTEP(fy0);

    iz = FLOOR(z);
    fz0 = z - iz;
    fz1 = fz0 - 1;
    wz = SMOOTHSTEP(fz0);

    vx0 = glattice(ix,iy,iz,fx0,fy0,fz0);
    vx1 = glattice(ix+1,iy,iz,fx1,fy0,fz0);
    vy0 = LERP(wx, vx0, vx1);
    vx0 = glattice(ix,iy+1,iz,fx0,fy1,fz0);
    vx1 = glattice(ix+1,iy+1,iz,fx1,fy1,fz0);
    vy1 = LERP(wx, vx0, vx1);
    vz0 = LERP(wy, vy0, vy1);

    vx0 = glattice(ix,iy,iz+1,fx0,fy0,fz1);
    vx1 = glattice(ix+1,iy,iz+1,fx1,fy0,fz1);
    vy0 = LERP(wx, vx0, vx1);
    vx0 = glattice(ix,iy+1,iz+1,fx0,fy1,fz1);
    vx1 = glattice(ix+1,iy+1,iz+1,fx1,fy1,fz1);
    vy1 = LERP(wx, vx0, vx1);
    vz1 = LERP(wy, vy0, vy1);

    return LERP(wz, vz0, vz1);
}

void fmat_rev(fmat_t *s) {
  uint_t i,j;
  for (i=0; i< s->height; i++) {
    for (j=0; j< FLOOR(s->length/2); j++) {
      ELEM_SWAP(s->data[i][j], s->data[i][s->length-1-j]);
    }
  }
}

static int64_t alignTimeDown(int64_t ts, int unit, int count)
{
	switch (unit) {
	case FT_SECOND:
		ts = FLOOR(ts, count * 1000);
		break;
	case FT_MINUTE:
		ts = FLOOR(ts, count * 60 * 1000);
		break;
	case FT_HOUR:
		ts = FLOOR(ts, count * 3600 * 1000);
	default:
		// TODO: more FT_XXX
		break;
	}

	return ts;
}