/**
 * Initialize Dense Inverse Log Polar Table
 * Fill each DILP[x][y] element with N <ecc,ang,dist> values so a
 * weighted sum can be computed for reconstrucion
 */
void LPolar::initDataDILP(int ecc, int ang, int s)
{
	float base,mod,lmod,alfa,lalfa,x,y;
	float radius= (float)s/2.;
	  		
	base = exp(log(radius)/((float)ecc));  //Such that pow(base,ecc) = radius

	for(int i=0; i<s; i++)
		for(int j=0; j<s; j++)
		{
			x = i - radius;
			y = j - radius;
			//Compute mod
			mod = sqrt(x*x +y*y);
			if(mod  >= radius) mod = radius-1;
			//Compute log-mod from base using log-base conversion. Result in 0..ecc range as computed in base
			lmod = log(mod)/log(base);
			//compute angle from atan2. Result in -Pi and Pi range
			alfa = atan2(y,x);
			//take alfa to 0..ang range
			lalfa = alfa*(float)ang/(2.* M_PI) + (ang/2);
			//closest ecc to lmod is (int)lmod.
			DILP[i][j].ecc = (int)rintf(lmod);
			DILP[i][j].ang = (int)rintf(lalfa);
			
		}
		
}

static void
sur_run(LV2_Handle instance, uint32_t n_samples)
{
	LV2meter* self = (LV2meter*)instance;
	uint32_t cors = self->chn > 3 ? 4 : 3;

	for (uint32_t c = 0; c < cors; ++c) {
		uint32_t in_a = rintf (*self->surc_a[c]);
		uint32_t in_b = rintf (*self->surc_b[c]);
		if (in_a >= self->chn) in_a = self->chn - 1;
		if (in_b >= self->chn) in_b = self->chn - 1;
		self->cor4[c]->process (self->input[in_a], self->input[in_b], n_samples);
		*self->surc_c[c] = self->cor4[c]->read();
	}

	for (uint32_t c = 0; c < self->chn; ++c) {
		float m, p;

		float* const input  = self->input[c];
		float* const output = self->output[c];

		self->mtr[c]->process(input, n_samples);

		static_cast<Kmeterdsp*>(self->mtr[c])->read(m, p);

		*self->level[c] = m;
		*self->peak[c]  = p;

		if (input != output) {
			memcpy(output, input, sizeof(float) * n_samples);
		}
	}
}

void addpoint(float x, float y, float z)
{
	if (points == NULL)
	{
		points = (vector3f*) malloc(sizeof(vector3f));
		count_points = 1;
	}
	else
	{
		count_points++;
		points = (vector3f*)realloc(points, count_points * sizeof(vector3f));
	}

	if( pieVersion == 2)
	{
		points[count_points-1].x = rintf(x);
		points[count_points-1].y = rintf(y);
		points[count_points-1].z = rintf(z);
	}
	else // pieVersion == 3
	{
		points[count_points-1].x = x;
		points[count_points-1].y = y;
		points[count_points-1].z = z;
	}
}

BRect
WorkspacesView::_WindowFrame(const BRect& workspaceFrame,
	const BRect& screenFrame, const BRect& windowFrame,
	BPoint windowPosition)
{
	BRect frame = windowFrame;
	frame.OffsetTo(windowPosition);

	// scale down the rect
	float factor = workspaceFrame.Width() / screenFrame.Width();
	frame.left = frame.left * factor;
	frame.right = frame.right * factor;

	factor = workspaceFrame.Height() / screenFrame.Height();
	frame.top = frame.top * factor;
	frame.bottom = frame.bottom * factor;

	// offset by the workspace fame position
	// and snap to integer coordinates without distorting the size too much
	frame.OffsetTo(rintf(frame.left + workspaceFrame.left),
		rintf(frame.top + workspaceFrame.top));
	frame.right = rintf(frame.right);
	frame.bottom = rintf(frame.bottom);

	return frame;
}

intvec * mri_bport_make_indexes( float dt, float fbot, float ftop, int ntime )
{
   int nfbot , nftop , nev=(ntime%2==0) , nth=ntime/2 , ff ;
   float df ;
   intvec *ivv ;

ENTRY("mri_bport_make_indexes") ;

   if( dt   <= 0.0f ) dt   = 1.0f ;
   if( fbot <  0.0f ) fbot = 0.0f ;
   if( ntime < 9 || ftop < fbot ) RETURN(NULL) ;

   df = 1.0f / (ntime *dt) ;  /* frequency spacing */

   nfbot = (int)rintf(fbot/df+0.1666666f) ;
   if( nfbot > nth      ) nfbot = nth ;
   if( nfbot < BP_ffbot ) nfbot = BP_ffbot ;

   nftop = (int)rintf(ftop/df-0.1666666f) ;
   if( nftop < nfbot   ) nftop = nfbot ;
   if( nftop > ntime/2 ) nftop = nth ;

   MAKE_intvec(ivv,nftop-nfbot+1) ;

   for( ff=nfbot ; ff <= nftop ; ff++ ) ivv->ar[ff-nfbot] = ff ;

#if 0
fprintf(stderr,"indexes: fbot=%g ftop=%g df=%g ntime=%d\n",fbot,ftop,df,ntime) ;
for( ff=nfbot ; ff <= nftop ; ff++ ) fprintf(stderr," %d",ivv->ar[ff-nfbot]) ;
fprintf(stderr,"\n") ;
#endif

   RETURN(ivv) ;
}

MRI_IMAGE * mri_bport_contig( float dt , float fbot , float ftop ,
                              int ntime , int nbefore , int nafter )
{
   MRI_IMAGE *fim ; float *far , *fii , *fip ;
   int nfbot , nftop , ff , ii,jj , nev=(ntime%2==0) , ncol,nrow ;
   float df , freq ;

ENTRY("mri_bport_contig") ;

   if( dt   <= 0.0f ) dt   = 1.0f ;
   if( fbot <  0.0f ) fbot = 0.0f ;
   if( ntime < 9 || ftop < fbot ) RETURN(NULL) ;
   if( nbefore < 0 ) nbefore = 0 ;
   if( nafter  < 0 ) nafter  = 0 ;

   df  = 1.0f / (ntime * dt) ;  /* frequency spacing */

   nfbot = (int)rintf(fbot/df+0.1666666f) ;
   if( nfbot > ntime/2    ) nfbot = ntime/2 ;
   if( nfbot < BP_ffbot   ) nfbot = BP_ffbot ;

   nftop = (int)rintf(ftop/df-0.1666666f) ;
   if( nftop < nfbot   ) nftop = nfbot   ;
   if( nftop > ntime/2 ) nftop = ntime/2 ;
#if 1
   ININFO_message("Frequency indexes: blocklen=%d nfbot=%d nftop=%d",ntime,nfbot,nftop) ;
#endif

   ncol = 2*(nftop-nfbot-1) ; if( ncol < 0 ) ncol = 0 ;
   ncol += (nfbot == 0 || (nfbot == ntime/2 && nev==1) ) ? 1 : 2 ;
   if( nftop > nfbot )
     ncol += (nftop == ntime/2 && nev) ? 1 : 2 ;

   if( ncol <= 0 ){
     ININFO_message("Failure :-(") ; RETURN(NULL) ;  /* should never happen */
   }

   nrow = ntime + nbefore + nafter ;
   fim  = mri_new( nrow , ncol , MRI_float ) ;
   far  = MRI_FLOAT_PTR(fim) ;

   for( ii=0,ff=nfbot ; ff <= nftop ; ff++ ){
     fii = far + (ii*nrow + nbefore) ;
     if( ff == 0 ){
       for( jj=0 ; jj < ntime ; jj++ ) fii[jj] = 1.0f ;
       ii++ ;
     } else if( ff == ntime/2 && nev ){
       for( jj=0 ; jj < ntime ; jj++ ) fii[jj] = 2*(jj%2)-1 ;
       ii++ ;
     } else {
       fip = fii + nrow ; freq = ff*df * (2.0f*3.141593f) * dt ;
       for( jj=0 ; jj < ntime ; jj++ ){
         fii[jj] = cos(freq*jj) ; fip[jj] = sin(freq*jj) ;
       }
       ii += 2 ;
     }
   }

   RETURN(fim) ;
}

TEST(math, rint) {
  auto guard = make_scope_guard([]() {
    fesetenv(FE_DFL_ENV);
  });

  fesetround(FE_UPWARD); // rint/rintf/rintl obey the rounding mode.
  feclearexcept(FE_ALL_EXCEPT); // rint/rintf/rintl do set the FE_INEXACT flag.
  ASSERT_EQ(1234.0, rint(1234.0));
  ASSERT_TRUE((fetestexcept(FE_ALL_EXCEPT) & FE_INEXACT) == 0);
  ASSERT_EQ(1235.0, rint(1234.01));
  ASSERT_TRUE((fetestexcept(FE_ALL_EXCEPT) & FE_INEXACT) != 0);

  feclearexcept(FE_ALL_EXCEPT); // rint/rintf/rintl do set the FE_INEXACT flag.
  ASSERT_EQ(1234.0f, rintf(1234.0f));
  ASSERT_TRUE((fetestexcept(FE_ALL_EXCEPT) & FE_INEXACT) == 0);
  ASSERT_EQ(1235.0f, rintf(1234.01f));
  ASSERT_TRUE((fetestexcept(FE_ALL_EXCEPT) & FE_INEXACT) != 0);

  feclearexcept(FE_ALL_EXCEPT); // rint/rintf/rintl do set the FE_INEXACT flag.
  ASSERT_EQ(1234.0, rintl(1234.0L));
  ASSERT_TRUE((fetestexcept(FE_ALL_EXCEPT) & FE_INEXACT) == 0);
  ASSERT_EQ(1235.0, rintl(1234.01L));
  ASSERT_TRUE((fetestexcept(FE_ALL_EXCEPT) & FE_INEXACT) != 0);

  fesetround(FE_TOWARDZERO); // rint/rintf obey the rounding mode.
  ASSERT_EQ(1234.0, rint(1234.01));
  ASSERT_EQ(1234.0f, rintf(1234.01f));
  ASSERT_EQ(1234.0, rintl(1234.01L));
}

static void
focusblur_fft_buffer_update_depth_table (FblurFftBuffer *fft,
                                         gint            division,
                                         gint            slide)
{
  FblurFftDepthTable    *table;
  gfloat                 dfac, fval;
  gfloat                 r, f, c;
  gint                   ri, fi, ci;
  gint                   i;

  g_assert (division > 0);

  if (division == fft->depth.division &&
      slide    == fft->depth.slide)
    return;

  dfac = (gfloat) FBLUR_DEPTH_MAX / division;

  for (i = 0; i <= FBLUR_DEPTH_MAX; i ++)
    {
      table = &(fft->depth.table[i]);
      fval = (gfloat) (i - slide) / dfac;

      r = rintf (fval);

      ri = rintf (r * dfac) + slide;
      ri = CLAMP (ri, 0, FBLUR_DEPTH_MAX);

      table->round = ri;

      if (fabsf (r - fval) < 0.001f)
        {
          table->floor = table->ceil = ri;
          table->diff = 0.0f;
        }
      else
        {
          f = floorf (fval);
          c = ceilf (fval);

          fi = rintf (f * dfac) + slide;
          ci = rintf (c * dfac) + slide;

          fi = CLAMP (fi, 0, FBLUR_DEPTH_MAX);
          ci = CLAMP (ci, 0, FBLUR_DEPTH_MAX);

          table->floor = fi;
          table->ceil  = ci;

          table->diff  = (ci > fi) ? (gfloat) (i - fi) / (ci - fi) : 0.0f;
        }
    }

  fft->depth.division = division;
  fft->depth.slide = slide;
  fft->depth.count = 0;
}

static float nearest_neighbor_interpolation_at(float *x,
        int w, int h, float p, float q)
{
	int ip = rintf(p);
	int iq = rintf(q);
	float r = getpixel_1(x, w, h, ip, iq);
    if (r < -1000) return NAN;
	return r;
}

void deth(int nr, TYPE *h) {
    const int hmax = 20;
    const int kmax = 30;
    const int lmax = 15;
    int i;
    for (i = 0; i < nr; i++) {
        h[DIM2_H * i + 0] = rintf(2 * hmax * (TYPE) random() / (TYPE) RAND_MAX - hmax);
        h[DIM2_H * i + 1] = rintf(2 * kmax * (TYPE) random() / (TYPE) RAND_MAX - kmax);
        h[DIM2_H * i + 2] = rintf(2 * lmax * (TYPE) random() / (TYPE) RAND_MAX - lmax);
    }
}

ISRF isPPostive2(SACHEAD *hdr, float *trace, float interval)
{
	int pre, suf;
	int index ;
	index = rintf( (hdr->o - hdr->b) / (hdr->delta) );
	pre = rintf( (hdr->o - hdr->b - interval ) / (hdr->delta) );
	suf = rintf( (hdr->o - hdr->b + interval ) / (hdr->delta) );
	if( integf(trace, pre, suf, hdr->delta) > 0.0f &&  trace[index] < 1.0f )
		return GOOD;
	else
		return BAD;
}

static int_pair zpadax_pm( int nx_super , float xorg_super ,
                           int nx_input , float xorg_input , float dx )
{
   int_pair pm ; float ts , ti ;

   ts   = xorg_super + nx_super * dx ;
   ti   = xorg_input + nx_input * dx ;
   pm.i = (int)rintf((xorg_input-xorg_super)/dx) ;
   pm.j = (int)rintf((ts-ti)/dx) ;

   return pm ;
}

/* check if we are still running and fill out the position 
	return == 0 we're still walking ... else we have reached a point */
inline int
ai_checkpos (_player * pl, _point * pos)
{
    _pointf _p;

    _p.x = CUTINT (pl->pos.x);
    _p.y = CUTINT (pl->pos.y);

    pos->x = rintf (pl->pos.x);
    pos->y = rintf (pl->pos.y);

    return ((_p.x < 0.15f || _p.x > 0.85f) && (_p.y < 0.15f || _p.y > 0.85f));
};

static int approximately_equal(float a, float b)
{
#ifdef STARPU_HAVE_NEARBYINTF
	int ai = (int) nearbyintf(a * 1000.0);
	int bi = (int) nearbyintf(b * 1000.0);
#elif defined(STARPU_HAVE_RINTF)
	int ai = (int) rintf(a * 1000.0);
	int bi = (int) rintf(b * 1000.0);
#else
#error "Please define either nearbyintf or rintf."
#endif
	return ai == bi;
}

static void plot_data_fft(SFSUI* ui) {
	cairo_t* cr;
	cr = cairo_create (ui->sf_dat);

	rounded_rectangle (cr, SS_BORDER, SS_BORDER, SS_SIZE, SS_SIZE, SS_BORDER);
	cairo_clip_preserve (cr);

	const float persistence = robtk_dial_get_value(ui->screen);
	float transp;
	cairo_set_operator (cr, CAIRO_OPERATOR_OVER);
	if (persistence > 0) {
		cairo_set_source_rgba(cr, 0, 0, 0, .25 - .0025 * persistence);
		transp = 0.05;
	} else {
		cairo_set_source_rgba(cr, 0, 0, 0, 1.0);
		transp = .5;
	}
	cairo_fill(cr);

	cairo_set_line_cap(cr, CAIRO_LINE_CAP_ROUND);
	cairo_set_line_width (cr, 1.0);

	const float xmid = rintf(SS_BORDER + SS_SIZE *.5) + .5;
	const float dnum = SS_SIZE / ui->log_base;
	const float denom = ui->log_rate / (float)ui->fft_bins;

	cairo_set_operator (cr, CAIRO_OPERATOR_OVER);
	for (uint32_t i = 1; i < ui->fft_bins-1 ; ++i) {
		if (ui->level[i] < 0) continue;

		const float level = MAKEUP_GAIN + fftx_power_to_dB(ui->level[i]);
		if (level < -80) continue;

		const float y  = rintf(SS_BORDER + SS_SIZE - dnum * fast_log10(1.0 + i * denom)) + .5;
		const float y1 = rintf(SS_BORDER + SS_SIZE - dnum * fast_log10(1.0 + (i+1) * denom)) + .5;
		const float pk = level > 0.0 ? 1.0 : (80 + level) / 80.0;
		const float a_lr = ui->lr[i];

		float clr[3];
		hsl2rgb(clr, .70 - .72 * pk, .9, .3 + pk * .4);
		cairo_set_source_rgba(cr, clr[0], clr[1], clr[2], transp  + pk * .2);
		cairo_set_line_width (cr, MAX(1.0, (y - y1)));

		cairo_move_to(cr, xmid, y);
		cairo_line_to(cr, SS_BORDER + SS_SIZE * a_lr, y);
		cairo_stroke(cr);

	}
	cairo_destroy (cr);
}

static void print_hz (char *t, float hz) {
	//printf("%f ", hz);
	hz = 5 * rintf(hz / 5.f);
	if (hz >= 990) {
		int dec = ((int)rintf (hz / 100.f)) % 10;
		if (dec != 0) {
			snprintf(t, 8, "%.0fK%d", floor(hz / 1000.f), dec);
		} else {
			snprintf(t, 8, "%.0fK", hz / 1000.f);
		}
	} else {
		snprintf(t, 8, "%.0f", hz);
	}
	//printf("-> %f -> %s\n", hz, t);
}

TEST (void)
{
  float a[NUM];
  float r[NUM];
  int i;

  init_src (a);

  for (i = 0; i < NUM; i++)
    r[i] = rintf (a[i]);

  /* check results:  */
  for (i = 0; i < NUM; i++)
    if (r[i] != rintf (a[i]))
      abort();
}

float EDIT_scale_misfit( int nxyz , float fac , short *sar , float *far )
{
    float sf , ff , sum=0.0f , df ;
    int ii , nf=0 ;

    ENTRY("EDIT_scale_misfit") ;

    if( nxyz <= 0 || sar == NULL || far == NULL ) RETURN(0.0f) ;

    if( fac == 0.0f ) fac = 1.0f ;
    df = 1.0f / fac ;

    for( ii=0 ; ii < nxyz ; ii++ ) {
        if( BAD(ii) ) continue ;
        ff = far[ii] ;
        if( ff == 0.0f ) continue ;
        sf = (short)rintf(fac*sar[ii]) ;
        if( sf == 0.0f ) {
            if( fabsf(ff) < df ) sum += fabsf(ff)*fac ;
            else                 sum += 1.0f ;
        } else {
            sf = fabsf((sf-ff)/ff) ;
            if( sf > 1.0f ) sf = 1.0f ;
            sum += sf ;
        }
        nf++ ;
    }

    if( nf > 0 ) sum /= nf ;
    RETURN(sum) ;
}

void
fun_C(int x, int y, int z)
{
	fesetround(FE_DOWNWARD);
	b=rintf(a);
	printf("%f\n", b);
}

void
enc_slinear_8(int fd, audio_encoding_t *enc, int chans, int rate)
{
	audio_info_t inf;
	struct ausrate rt;
	int8_t *samples = NULL, *p;
	int i, j;

	AUDIO_INITINFO(&inf);
	inf.play.precision = enc->precision;
	inf.play.encoding = enc->encoding;
	inf.play.channels = chans;
	inf.play.sample_rate = rate;; 

	if (ioctl(fd, AUDIO_SETINFO, &inf) == -1) {
		printf("[%s]", strerror(errno));
		goto out;
	}

	if (ioctl(fd, AUDIO_GETINFO, &inf) == -1) {
		printf("[getinfo: %s]", strerror(errno));
		goto out;
	}
	rt.r_rate = inf.play.sample_rate;
	rt.s_rate = inf.play.sample_rate;
	rt.bps = 1 * chans;
	rt.bytes = inf.play.sample_rate * chans * PLAYSECS;

	samples = (int8_t *)malloc(inf.play.sample_rate * chans);
	if (samples == NULL) {
		warn("malloc");
		goto out;
	}

	for (i = 0, p = samples; i < inf.play.sample_rate; i++) {
		float d;
		int8_t v;

		d = 127.0 * sinf(((float)i / (float)inf.play.sample_rate) *
		    (2 * M_PI * playfreq));
		d = rintf(d);
		v = d;

		for (j = 0; j < chans; j++) {
			*p = v;
			p++;
		}
	}

	mark_time(&rt.tv_begin);
	for (i = 0; i < PLAYSECS; i++)
		write(fd, samples, inf.play.sample_rate * chans);
	audio_wait(fd);
	mark_time(&rt.tv_end);
	check_srate(&rt);

out:
	if (samples != NULL)
		free(samples);
}