text:
#include <math.h>
#include <stdio.h>
#include <mpfr.h>
#include <gmp.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <gsl/gsl_sf_legendre.h>
#include <petsc.h>

#undef __FUNCT__
#define __FUNCT__ "DEwk"
PetscErrorCode DEwk(PetscInt k, mpfr_t h, mpfr_t *wk)
{
  PetscErrorCode ierr;
  mpfr_t kh, msinh, mcosh, piOver2, tmp;
  PetscFunctionBegin;

  mpfr_inits(kh, msinh, mcosh, piOver2, tmp, NULL);

  //init piOver2
  mpfr_const_pi(piOver2, MPFR_RNDN);
  mpfr_mul_d(piOver2, piOver2, 0.5, MPFR_RNDN);

  mpfr_set_d(kh, (double) k, MPFR_RNDN);
  mpfr_mul(kh, kh, h, MPFR_RNDN);
  mpfr_set(*wk, h, MPFR_RNDN);
  mpfr_sinh_cosh(msinh, mcosh, kh, MPFR_RNDN);
  mpfr_mul(msinh, msinh, piOver2, MPFR_RNDN);
  mpfr_mul(mcosh, mcosh, piOver2, MPFR_RNDN);
  mpfr_cosh(tmp, msinh, MPFR_RNDN);
  mpfr_sqr(tmp, tmp, MPFR_RNDN);
  mpfr_mul(*wk, *wk, mcosh, MPFR_RNDN);
  mpfr_div(*wk, *wk, tmp, MPFR_RNDN);

  PetscFunctionReturn(0);
}

/* have to make single exponential */
#undef __FUNCT__
#define __FUNCT__ "SEwk"
PetscErrorCode SEwk(PetscInt k, mpfr_t h, mpfr_t *wk)
{
  PetscErrorCode ierr;
  mpfr_t kh, msinh, mcosh, piOver2, tmp;
  PetscFunctionBegin;

  mpfr_inits(kh, msinh, mcosh, piOver2, tmp, NULL);

  //init piOver2
  mpfr_const_pi(piOver2, MPFR_RNDN);
  mpfr_mul_d(piOver2, piOver2, 0.5, MPFR_RNDN);

  mpfr_set_d(kh, (double) k, MPFR_RNDN);
  mpfr_mul(kh, kh, h, MPFR_RNDN);
  //printf("wow kh: %15.15f\n", mpfr_get_d(kh, MPFR_RNDN));
  mpfr_tanh(*wk, kh, MPFR_RNDN);
  //printf("wow wk: %15.15f\n", mpfr_get_d(*wk, MPFR_RNDN));
  mpfr_mul(*wk, *wk, *wk, MPFR_RNDN);
  mpfr_d_sub(*wk, 1.0, *wk, MPFR_RNDN);
  mpfr_mul(*wk, *wk, h, MPFR_RNDN);
  //mpfr_set(*wk, h, MPFR_RNDN);
  //mpfr_sinh_cosh(msinh, mcosh, kh, MPFR_RNDN);
  //mpfr_mul(msinh, msinh, piOver2, MPFR_RNDN);
  //mpfr_mul(mcosh, mcosh, piOver2, MPFR_RNDN);
  //mpfr_cosh(tmp, msinh, MPFR_RNDN);
  //mpfr_sqr(tmp, tmp, MPFR_RNDN);
  //mpfr_mul(*wk, *wk, mcosh, MPFR_RNDN);
  //mpfr_div(*wk, *wk, tmp, MPFR_RNDN);

  PetscFunctionReturn(0);
}

/* have to make single exponential */
#undef __FUNCT__
#define __FUNCT__ "ERFwk"
PetscErrorCode ERFwk(PetscInt k, mpfr_t h, mpfr_t *wk)
{
  PetscErrorCode ierr;
  mpfr_t kh, msinh, mcosh, piOver2, tmp;
  PetscFunctionBegin;

  mpfr_inits(kh, msinh, mcosh, piOver2, tmp, NULL);

  // init constant
  mpfr_const_pi(piOver2, MPFR_RNDN);
  mpfr_sqrt(piOver2, piOver2, MPFR_RNDN);
  mpfr_d_div(piOver2, 2.0, piOver2, MPFR_RNDN);

  mpfr_set_d(kh, (double) k, MPFR_RNDN);
  mpfr_mul(kh, kh, h, MPFR_RNDN);
  mpfr_mul(tmp, kh, kh, MPFR_RNDN);
  mpfr_mul_d(tmp, tmp, -1.0, MPFR_RNDN);
  mpfr_exp(*wk, tmp, MPFR_RNDN);
  mpfr_mul(*wk, *wk, h, MPFR_RNDN);
  mpfr_mul(*wk, *wk, piOver2, MPFR_RNDN);

  PetscFunctionReturn(0);
}

/* generates n weights spaced h apart in the positive direction */
/* includes "center" n=0 weight as first entry */
#undef __FUNCT__
#define __FUNCT__ "DExkwk"
PetscErrorCode DExkwk(PetscInt n, mpfr_t h, mpfr_t **xk, mpfr_t **wk)
{
  PetscErrorCode ierr;
  PetscInt k;
  mpfr_t kh, msinh, mcosh, piOver2, tmp;
  PetscReal alp, bet;
  PetscInt flopCount;
  PetscFunctionBegin;
  flopCount = 0;

  mpfr_inits(kh, msinh, mcosh, piOver2, tmp, NULL);

  // init piOver2
  mpfr_const_pi(piOver2, MPFR_RNDN);
  mpfr_mul_d(piOver2, piOver2, 0.5, MPFR_RNDN); flopCount ++;

  for(k=0; k<n; ++k) {
    /* calculate wk */
    mpfr_set_d(kh, (double) k, MPFR_RNDN);
    mpfr_mul(kh, kh, h, MPFR_RNDN);
    mpfr_set((*wk)[k], h, MPFR_RNDN);
    mpfr_sinh_cosh(msinh, mcosh, kh, MPFR_RNDN);
    mpfr_mul(msinh, msinh, piOver2, MPFR_RNDN);
    mpfr_mul(mcosh, mcosh, piOver2, MPFR_RNDN);
    mpfr_cosh(tmp, msinh, MPFR_RNDN);
    mpfr_sqr(tmp, tmp, MPFR_RNDN);
    mpfr_mul((*wk)[k], (*wk)[k], mcosh, MPFR_RNDN);
    mpfr_div((*wk)[k], (*wk)[k], tmp, MPFR_RNDN); flopCount += 9;

    /* calculate xk */
    mpfr_set_d((*xk)[k], 1.0, MPFR_RNDZ);
    mpfr_cosh(tmp, msinh, MPFR_RNDN);
    mpfr_div((*xk)[k], (*xk)[k], tmp, MPFR_RNDZ);
    mpfr_exp(tmp, msinh, MPFR_RNDN);
    mpfr_div((*xk)[k], (*xk)[k], tmp, MPFR_RNDZ); flopCount += 4;

  }
  mpfr_clears(kh, msinh, mcosh, piOver2, tmp, NULL);

  ierr = PetscLogFlops(flopCount);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

/* generates n weights spaced h apart in the positive direction */
/* includes "center" n=0 weight as first entry */
#undef __FUNCT__
#define __FUNCT__ "SExkwk"
PetscErrorCode SExkwk(PetscInt n, mpfr_t h, mpfr_t **xk, mpfr_t **wk)
{
  PetscErrorCode ierr;
  PetscInt k;
  mpfr_t kh, msinh, mcosh, piOver2, tmp;
  PetscReal alp, bet;
  PetscInt flopCount;
  PetscFunctionBegin;
  flopCount = 0;

  mpfr_inits(kh, msinh, mcosh, piOver2, tmp, NULL);

  // init piOver2
  mpfr_const_pi(piOver2, MPFR_RNDN);
  mpfr_mul_d(piOver2, piOver2, 0.5, MPFR_RNDN); flopCount ++;

  for(k=0; k<n; ++k) {
    /* calculate wk */
    mpfr_set_d(kh, (double) k, MPFR_RNDN);
    mpfr_mul(kh, kh, h, MPFR_RNDN);
    mpfr_tanh((*wk)[k], kh, MPFR_RNDN);
    mpfr_mul((*wk)[k], (*wk)[k], (*wk)[k], MPFR_RNDN);
    mpfr_d_sub((*wk)[k], 1.0, (*wk)[k], MPFR_RNDN);
    mpfr_mul((*wk)[k], (*wk)[k], h, MPFR_RNDN); flopCount += 5;

    /* calculate xk */
    mpfr_tanh((*xk)[k], kh, MPFR_RNDN); flopCount += 1;
    //mpfr_set_d((*xk)[k], 1.0, MPFR_RNDZ);
    //mpfr_cosh(tmp, msinh, MPFR_RNDN);
    //mpfr_div((*xk)[k], (*xk)[k], tmp, MPFR_RNDZ);
    //mpfr_exp(tmp, msinh, MPFR_RNDN);
    //mpfr_div((*xk)[k], (*xk)[k], tmp, MPFR_RNDZ); flopCount += 4;
    //mpfr_d_sub(tmp, 1.0, (*xk)[k], MPFR_RNDN);
    //printf("xk: %15.15e\n1-xk: %15.15e\nwk: %3.3e\n\n", mpfr_get_d((*xk)[k], MPFR_RNDN), mpfr_get_d(tmp, MPFR_RNDN), mpfr_get_d((*wk)[k], MPFR_RNDN));
  }
  mpfr_clears(kh, msinh, mcosh, piOver2, tmp, NULL);

  ierr = PetscLogFlops(flopCount);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

/* generates n weights spaced h apart in the positive direction */
/* includes "center" n=0 weight as first entry */
#undef __FUNCT__
#define __FUNCT__ "ERFxkwk"
PetscErrorCode ERFxkwk(PetscInt n, mpfr_t h, mpfr_t **xk, mpfr_t **wk)
{
  PetscErrorCode ierr;
  PetscInt k;
  mpfr_t kh, msinh, mcosh, piOver2, tmp;
  PetscReal alp, bet;
  PetscInt flopCount;
  PetscFunctionBegin;
  flopCount = 0;

  mpfr_inits(kh, msinh, mcosh, piOver2, tmp, NULL);

  // init piOver2
  mpfr_const_pi(piOver2, MPFR_RNDN);
  mpfr_sqrt(piOver2, piOver2, MPFR_RNDN);
  mpfr_d_div(piOver2, 2.0, piOver2, MPFR_RNDN);
  //mpfr_mul_d(piOver2, piOver2, 0.5, MPFR_RNDN); flopCount ++;

  for(k=0; k<n; ++k) {
    /* calculate wk */
    mpfr_set_d(kh, (double) k, MPFR_RNDN);
    mpfr_mul(kh, kh, h, MPFR_RNDN);
    mpfr_mul(tmp, kh, kh, MPFR_RNDN);
    mpfr_mul_d(tmp, tmp, -1.0, MPFR_RNDN);
    mpfr_exp((*wk)[k], tmp, MPFR_RNDN);
    mpfr_mul((*wk)[k], (*wk)[k], h, MPFR_RNDN);
    mpfr_mul((*wk)[k], (*wk)[k], piOver2, MPFR_RNDN);

    /* calculate xk */
    mpfr_erf((*xk)[k], kh, MPFR_RNDN);
  }
  mpfr_clears(kh, msinh, mcosh, piOver2, tmp, NULL);

  ierr = PetscLogFlops(flopCount);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "inte"
PetscErrorCode inte(mpfr_t *x, mpfr_t *val, void *ctx)
{
  PetscErrorCode ierr;
  PetscInt flopCount;
  PetscFunctionBegin;
  flopCount = 0;

  //mpfr_mul(*val, *x, *x, MPFR_RNDN);
  //mpfr_mul(*val, *val, *x, MPFR_RNDN); flopCount += 2;
  mpfr_const_pi(*val, MPFR_RNDN);
  mpfr_mul_d(*val, *val, 2.0, MPFR_RNDN);
  mpfr_mul(*val, *val, *x, MPFR_RNDN);
  mpfr_sin(*val, *val, MPFR_RNDN);
  mpfr_mul(*val, *val, *x, MPFR_RNDN);
  //printf("xk: %15.15f\n", mpfr_get_d(*x, MPFR_RNDN));
  //printf("val: %15.15f\n", mpfr_get_d(*val, MPFR_RNDN));
  //mpfr_set_d(*val, 1.0, MPFR_RNDN);
  //mpfr_set(*val, *x, MPFR_RNDN);

  ierr = PetscLogFlops(flopCount);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "FindStepSize"
PetscErrorCode FindStepSize(PetscInt prec, PetscInt nPts, mpfr_t *step)
{
  PetscErrorCode ierr;
  const PetscReal tol = 1e-10;
  mpfr_t wk, tmp, h, hl, hr, lwk, rwk, swk, slwk, srwk, error;
  PetscInt cont;
  PetscFunctionBegin;

  if(prec < 1) {
    printf("prec should be positive\n");
    PetscFunctionReturn(1);
  }

  mpfr_inits(wk, tmp, h, hl, hr, lwk, rwk, swk, slwk, srwk, error, NULL);

  mpfr_set_d(hl, 0.0, MPFR_RNDN);
  mpfr_set_d(hr, 2.0, MPFR_RNDN);

  /* make hr large enough */
  mpfr_set_d(hr, 1.0, MPFR_RNDN);
  cont = 1;
  while(cont == 1) {
    ierr = DEwk(nPts, hr, &rwk);CHKERRQ(ierr);
    mpfr_log10(srwk, rwk, MPFR_RNDN);
    if(mpfr_get_d(srwk, MPFR_RNDN) < -2*prec)
      cont = 0;
    else
      mpfr_mul_d(hr, hr, 2.0, MPFR_RNDN); // make hr larger if value not small
  }
  //printf("hr: %15.15f\n", mpfr_get_d(hr, MPFR_RNDN));
  /* make hl small enough */
  mpfr_mul_d(hl, hr, 0.5, MPFR_RNDN);
  cont = 1;
  while(cont == 1) {
    ierr = DEwk(nPts, hl, &lwk);CHKERRQ(ierr);
    mpfr_log10(slwk, lwk, MPFR_RNDN);
    if(mpfr_get_d(slwk, MPFR_RNDN) > -2*prec)
      cont = 0;
    else
      mpfr_mul_d(hl, hl, 0.5, MPFR_RNDN); // make hl larger if value not small
  }
  //printf("hl: %15.15f\n", mpfr_get_d(hl, MPFR_RNDN));
  //ierr = DEwk(nPts, hl, &lwk);CHKERRQ(ierr);
  //ierr = DEwk(nPts, hr, &rwk);CHKERRQ(ierr);

  cont = 1;
  while(cont == 1) {
    // average hl and hr
    mpfr_add(h, hl, hr, MPFR_RNDN);
    mpfr_mul_d(h, h, 0.5, MPFR_RNDN);

    // calculate weight and size
    ierr = DEwk(nPts, h, &wk);CHKERRQ(ierr);
    ierr = DEwk(nPts, hl, &lwk);CHKERRQ(ierr);
    ierr = DEwk(nPts, hr, &rwk);CHKERRQ(ierr);
    mpfr_log10(swk, wk, MPFR_RNDN);
    mpfr_log10(slwk, lwk, MPFR_RNDN);
    mpfr_log10(srwk, rwk, MPFR_RNDN);
    if( mpfr_get_d(swk, MPFR_RNDN) < -2*prec)
      mpfr_set(hr, h, MPFR_RNDN);
    else
      mpfr_set(hl, h, MPFR_RNDN);
    // calculate maximum possible error
    mpfr_sub(error, slwk, srwk, MPFR_RNDN);
    mpfr_abs(error, error, MPFR_RNDN);
    if(mpfr_get_d(error, MPFR_RNDN) < tol)
       cont = 0;
  }

  mpfr_set(*step, h, MPFR_RNDN);

  mpfr_clears(wk, tmp, h, hl, hr, lwk, rwk, swk, slwk, srwk, error, NULL);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "FindStepSizeSE"
PetscErrorCode FindStepSizeSE(PetscInt prec, PetscInt nPts, mpfr_t *step)
{
  PetscErrorCode ierr;
  const PetscReal tol = 1e-25;
  mpfr_t wk, tmp, h, hl, hr, lwk, rwk, swk, slwk, srwk, error;
  PetscInt cont;
  PetscFunctionBegin;

  if(prec < 1) {
    printf("prec should be positive\n");
    PetscFunctionReturn(1);
  }

  mpfr_inits(wk, tmp, h, hl, hr, lwk, rwk, swk, slwk, srwk, error, NULL);

  mpfr_set_d(hl, 0.0, MPFR_RNDN);
  mpfr_set_d(hr, 2.0, MPFR_RNDN);

  /* make hr large enough */
  mpfr_set_d(hr, 1.0, MPFR_RNDN);
  cont = 1;
  while(cont == 1) {
    ierr = SEwk(nPts, hr, &rwk);CHKERRQ(ierr);
    mpfr_log10(srwk, rwk, MPFR_RNDN);
    if(mpfr_get_d(srwk, MPFR_RNDN) < -2*prec)
      cont = 0;
    else
      mpfr_mul_d(hr, hr, 2.0, MPFR_RNDN); // make hr larger if value not small
  }
  //printf("hr: %15.15f\n", mpfr_get_d(hr, MPFR_RNDN));
  /* make hl small enough */
  mpfr_mul_d(hl, hr, 0.5, MPFR_RNDN);
  cont = 1;
  while(cont == 1) {
    ierr = SEwk(nPts, hl, &lwk);CHKERRQ(ierr);
    mpfr_log10(slwk, lwk, MPFR_RNDN);
    if(mpfr_get_d(slwk, MPFR_RNDN) > -2*prec)
      cont = 0;
    else
      mpfr_mul_d(hl, hl, 0.5, MPFR_RNDN); // make hl larger if value not small
  }
  //printf("hl: %15.15f\n", mpfr_get_d(hl, MPFR_RNDN));
  ierr = SEwk(nPts, hl, &lwk);CHKERRQ(ierr);
  ierr = SEwk(nPts, hr, &rwk);CHKERRQ(ierr);
  //printf("lwk: %3.3e\n", mpfr_get_d(lwk, MPFR_RNDN));
  //printf("rwk: %3.3e\n", mpfr_get_d(rwk, MPFR_RNDN));

  cont = 1;
  while(cont == 1) {
    // average hl and hr
    mpfr_add(h, hl, hr, MPFR_RNDN);
    mpfr_mul_d(h, h, 0.5, MPFR_RNDN);

    // calculate weight and size
    ierr = SEwk(nPts, h, &wk);CHKERRQ(ierr);
    ierr = SEwk(nPts, hl, &lwk);CHKERRQ(ierr);
    ierr = SEwk(nPts, hr, &rwk);CHKERRQ(ierr);
    mpfr_log10(swk, wk, MPFR_RNDN);
    mpfr_log10(slwk, lwk, MPFR_RNDN);
    mpfr_log10(srwk, rwk, MPFR_RNDN);
    if( mpfr_get_d(swk, MPFR_RNDN) < -2*prec)
      mpfr_set(hr, h, MPFR_RNDN);
    else
      mpfr_set(hl, h, MPFR_RNDN);
    // calculate maximum possible error
    mpfr_sub(error, slwk, srwk, MPFR_RNDN);
    mpfr_abs(error, error, MPFR_RNDN);
    if(mpfr_get_d(error, MPFR_RNDN) < tol)
       cont = 0;
    //printf("h: %15.15f\nlwk: %3.3e\nrwk: %3.3e\n", mpfr_get_d(h, MPFR_RNDN), mpfr_get_d(lwk, MPFR_RNDN), mpfr_get_d(rwk, MPFR_RNDN));
  }

  mpfr_set(*step, h, MPFR_RNDN);

  mpfr_clears(wk, tmp, h, hl, hr, lwk, rwk, swk, slwk, srwk, error, NULL);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "FindStepSizeERF"
PetscErrorCode FindStepSizeERF(PetscInt prec, PetscInt nPts, mpfr_t *step)
{
  PetscErrorCode ierr;
  const PetscReal tol = 1e-25;
  mpfr_t wk, tmp, h, hl, hr, lwk, rwk, swk, slwk, srwk, error;
  PetscInt cont;
  PetscFunctionBegin;

  if(prec < 1) {
    printf("prec should be positive\n");
    PetscFunctionReturn(1);
  }

  mpfr_inits(wk, tmp, h, hl, hr, lwk, rwk, swk, slwk, srwk, error, NULL);

  mpfr_set_d(hl, 0.0, MPFR_RNDN);
  mpfr_set_d(hr, 2.0, MPFR_RNDN);

  /* make hr large enough */
  mpfr_set_d(hr, 1.0, MPFR_RNDN);
  cont = 1;
  while(cont == 1) {
    ierr = ERFwk(nPts, hr, &rwk);CHKERRQ(ierr);
    mpfr_log10(srwk, rwk, MPFR_RNDN);
    if(mpfr_get_d(srwk, MPFR_RNDN) < -2*prec)
      cont = 0;
    else
      mpfr_mul_d(hr, hr, 2.0, MPFR_RNDN); // make hr larger if value not small
  }
  //printf("hr: %15.15f\n", mpfr_get_d(hr, MPFR_RNDN));
  /* make hl small enough */
  mpfr_mul_d(hl, hr, 0.5, MPFR_RNDN);
  cont = 1;
  while(cont == 1) {
    ierr = ERFwk(nPts, hl, &lwk);CHKERRQ(ierr);
    mpfr_log10(slwk, lwk, MPFR_RNDN);
    if(mpfr_get_d(slwk, MPFR_RNDN) > -2*prec)
      cont = 0;
    else
      mpfr_mul_d(hl, hl, 0.5, MPFR_RNDN); // make hl larger if value not small
  }
  //printf("hl: %15.15f\n", mpfr_get_d(hl, MPFR_RNDN));
  ierr = ERFwk(nPts, hl, &lwk);CHKERRQ(ierr);
  ierr = ERFwk(nPts, hr, &rwk);CHKERRQ(ierr);
  //printf("lwk: %3.3e\n", mpfr_get_d(lwk, MPFR_RNDN));
  //printf("rwk: %3.3e\n", mpfr_get_d(rwk, MPFR_RNDN));

  cont = 1;
  while(cont == 1) {
    // average hl and hr
    mpfr_add(h, hl, hr, MPFR_RNDN);
    mpfr_mul_d(h, h, 0.5, MPFR_RNDN);

    // calculate weight and size
    ierr = ERFwk(nPts, h, &wk);CHKERRQ(ierr);
    ierr = ERFwk(nPts, hl, &lwk);CHKERRQ(ierr);
    ierr = ERFwk(nPts, hr, &rwk);CHKERRQ(ierr);
    mpfr_log10(swk, wk, MPFR_RNDN);
    mpfr_log10(slwk, lwk, MPFR_RNDN);
    mpfr_log10(srwk, rwk, MPFR_RNDN);
    if( mpfr_get_d(swk, MPFR_RNDN) < -2*prec)
      mpfr_set(hr, h, MPFR_RNDN);
    else
      mpfr_set(hl, h, MPFR_RNDN);
    // calculate maximum possible error
    mpfr_sub(error, slwk, srwk, MPFR_RNDN);
    mpfr_abs(error, error, MPFR_RNDN);
    if(mpfr_get_d(error, MPFR_RNDN) < tol)
       cont = 0;
    //printf("h: %15.15f\n", mpfr_get_d(h, MPFR_RNDN));
  }
  //printf("h: %15.15f\n", mpfr_get_d(h, MPFR_RNDN));
  //printf("h: %15.15f\nlwk: %3.3e\nrwk: %3.3e\n", mpfr_get_d(h, MPFR_RNDN), mpfr_get_d(lwk, MPFR_RNDN), mpfr_get_d(rwk, MPFR_RNDN));
  mpfr_set(*step, h, MPFR_RNDN);

  mpfr_clears(wk, tmp, h, hl, hr, lwk, rwk, swk, slwk, srwk, error, NULL);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "DEQuad"
PetscErrorCode DEQuad(PetscErrorCode (*f)(mpfr_t*,mpfr_t*,void*), mpfr_t a, mpfr_t b, PetscInt prec, PetscInt nPts, PetscReal *integral, void *ctx)
{
  PetscErrorCode ierr;
  PetscInt k;
  PetscInt nSide;
  mpfr_t alpha, beta;
  mpfr_t leftXk, rightXk;
  mpfr_t wk;
  mpfr_t *wkVals;
  mpfr_t *xkVals;
  mpfr_t lx, rx;
  mpfr_t leftFx, rightFx;
  mpfr_t pi2;
  mpfr_t kh, h;
  mpfr_t sum;
  PetscInt flopCount;
  PetscFunctionBegin;
  flopCount = 0;

  if(nPts % 2 == 0) {
    nSide = nPts/2; flopCount++;
  }
  if(nPts % 2 == 1) {
    nSide = (nPts-1)/2; flopCount++;
  }

  /* initializations */
  mpfr_inits(alpha, beta, leftXk, rightXk, wk, leftFx, rightFx, pi2, kh, h, sum, lx, rx, NULL);
  wkVals = (mpfr_t*) malloc(sizeof(mpfr_t)*(nSide+1));
  xkVals = (mpfr_t*) malloc(sizeof(mpfr_t)*(nSide+1));
  for(k=0; k<nSide+1; ++k)
    mpfr_inits(xkVals[k], wkVals[k], NULL);

  // init alpha, beta
  mpfr_set(alpha, b, MPFR_RNDN);
  mpfr_sub(alpha, alpha, a, MPFR_RNDN);
  mpfr_mul_d(alpha, alpha, 0.5, MPFR_RNDN);
  mpfr_set(beta, b, MPFR_RNDN);
  mpfr_add(beta, beta, a, MPFR_RNDN);
  mpfr_mul_d(beta, beta, 0.5, MPFR_RNDN); flopCount += 4;

  /* determine right step size */
  ierr = FindStepSize(prec, nSide, &h); // not including in flop count
  //printf("DE step size: %15.15f\n", mpfr_get_d(h, MPFR_RNDN));
  //mpfr_set_d(h, 0.04, MPFR_RNDN);

  /* calculate abscissas, weights */
  ierr = DExkwk(nSide+1, h, &xkVals, &wkVals);CHKERRQ(ierr);

  /* center term */
  mpfr_set_d(sum, 0.0, MPFR_RNDN);
  mpfr_set(lx, b, MPFR_RNDN);
  mpfr_add(lx, lx, a, MPFR_RNDN);
  mpfr_mul_d(lx, lx, .5, MPFR_RNDN);
  f(&lx, &sum, ctx);
  mpfr_mul(sum, sum, wkVals[0], MPFR_RNDN);
  mpfr_mul(sum, sum, alpha,     MPFR_RNDN); flopCount += 4;

  for(k=1; k <= nSide; ++k) {
    // adjust xk to [a,b]
    mpfr_sub_d(lx, xkVals[k], 1.0, MPFR_RNDZ);
    mpfr_mul(lx, lx, alpha, MPFR_RNDU);
    mpfr_add(lx, lx, beta, MPFR_RNDU);
    mpfr_d_sub(rx, 1.0, xkVals[k], MPFR_RNDZ);
    mpfr_mul(rx, rx, alpha, MPFR_RNDD);
    mpfr_add(rx, rx, beta , MPFR_RNDD); flopCount += 6;

    double left = mpfr_get_d(lx, MPFR_RNDN);
    double right = mpfr_get_d(rx, MPFR_RNDN);
    double xk = mpfr_get_d(xkVals[k], MPFR_RNDN);

    /* calculate function values */
    f(&lx, &leftFx , ctx);
    f(&rx, &rightFx, ctx);

    /* update sum */
    mpfr_mul(leftFx, leftFx, wkVals[k], MPFR_RNDN);
    mpfr_mul(rightFx, rightFx, wkVals[k], MPFR_RNDN);
    mpfr_mul(leftFx, leftFx, alpha, MPFR_RNDN);
    mpfr_mul(rightFx, rightFx, alpha, MPFR_RNDN);
    //double l1 = mpfr_get_d(leftFx, MPFR_RNDN);
    //double l2 = mpfr_get_d(rightFx, MPFR_RNDN);
    //printf("leftval: %15.15f\n", l1);
    //printf("rightval: %15.15f\n", l2);
    mpfr_add(sum, sum, leftFx, MPFR_RNDN);
    mpfr_add(sum, sum, rightFx, MPFR_RNDN); flopCount += 6;
  }
  PetscReal prog = mpfr_get_d(sum, MPFR_RNDN);
  //printf("sum: %15.15f\n", prog);

  *integral = prog;

  mpfr_clears(alpha, beta, leftXk, rightXk, wk, leftFx, rightFx, pi2, kh, h, sum, lx, rx, NULL);
  for(k=0; k<nSide+1; ++k)
    mpfr_clears(wkVals[k], xkVals[k], NULL);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "SEQuad"
PetscErrorCode SEQuad(PetscErrorCode (*f)(mpfr_t*,mpfr_t*,void*), mpfr_t a, mpfr_t b, PetscInt prec, PetscInt nPts, PetscReal *integral, void *ctx)
{
  PetscErrorCode ierr;
  PetscInt k;
  PetscInt nSide;
  mpfr_t alpha, beta;
  mpfr_t leftXk, rightXk;
  mpfr_t wk;
  mpfr_t *wkVals;
  mpfr_t *xkVals;
  mpfr_t lx, rx;
  mpfr_t leftFx, rightFx;
  mpfr_t pi2;
  mpfr_t kh, h;
  mpfr_t sum;
  PetscInt flopCount;
  PetscFunctionBegin;
  flopCount = 0;

  if(nPts % 2 == 0) {
    nSide = nPts/2; flopCount++;
  }
  if(nPts % 2 == 1) {
    nSide = (nPts-1)/2; flopCount++;
  }

  /* initializations */
  mpfr_inits(alpha, beta, leftXk, rightXk, wk, leftFx, rightFx, pi2, kh, h, sum, lx, rx, NULL);
  wkVals = (mpfr_t*) malloc(sizeof(mpfr_t)*(nSide+1));
  xkVals = (mpfr_t*) malloc(sizeof(mpfr_t)*(nSide+1));
  for(k=0; k<nSide+1; ++k)
    mpfr_inits(xkVals[k], wkVals[k], NULL);

  // init alpha, beta
  mpfr_set(alpha, b, MPFR_RNDN);
  mpfr_sub(alpha, alpha, a, MPFR_RNDN);
  mpfr_mul_d(alpha, alpha, 0.5, MPFR_RNDN);
  mpfr_set(beta, b, MPFR_RNDN);
  mpfr_add(beta, beta, a, MPFR_RNDN);
  mpfr_mul_d(beta, beta, 0.5, MPFR_RNDN); flopCount += 4;

  /* determine right step size */
  ierr = FindStepSizeSE(prec, nSide, &h); // not including in flop count
  //printf("step size: %15.15f\n", mpfr_get_d(h, MPFR_RNDN));
  //printf("step size: %15.15f\n", mpfr_get_d(h, MPFR_RNDN));
  //mpfr_set_d(h, 0.04, MPFR_RNDN);

  /* calculate abscissas, weights */
  ierr = SExkwk(nSide+1, h, &xkVals, &wkVals);CHKERRQ(ierr);

  /* center term */
  mpfr_set_d(sum, 0.0, MPFR_RNDN);
  mpfr_set(lx, b, MPFR_RNDN);
  mpfr_add(lx, lx, a, MPFR_RNDN);
  mpfr_mul_d(lx, lx, .5, MPFR_RNDN);
  f(&lx, &sum, ctx);
  mpfr_mul(sum, sum, wkVals[0], MPFR_RNDN);
  mpfr_mul(sum, sum, alpha,     MPFR_RNDN); flopCount += 4;

  for(k=1; k <= nSide; ++k) {
    // adjust xk to [a,b]
    //mpfr_sub_d(lx, xkVals[k], 1.0, MPFR_RNDZ);
    mpfr_mul_d(lx, xkVals[k], -1.0, MPFR_RNDN);
    mpfr_mul(lx, lx, alpha, MPFR_RNDU);
    mpfr_add(lx, lx, beta, MPFR_RNDU);
    //mpfr_d_sub(rx, 1.0, xkVals[k], MPFR_RNDZ);
    mpfr_set(rx, xkVals[k], MPFR_RNDN);
    mpfr_mul(rx, rx, alpha, MPFR_RNDD);
    mpfr_add(rx, rx, beta , MPFR_RNDD); flopCount += 6;

    double left = mpfr_get_d(lx, MPFR_RNDN);
    double right = mpfr_get_d(rx, MPFR_RNDN);
    double xk = mpfr_get_d(xkVals[k], MPFR_RNDN);

    /* calculate function values */
    f(&lx, &leftFx , ctx);
    f(&rx, &rightFx, ctx);

    /* update sum */
    mpfr_mul(leftFx, leftFx, wkVals[k], MPFR_RNDN);
    mpfr_mul(rightFx, rightFx, wkVals[k], MPFR_RNDN);
    mpfr_mul(leftFx, leftFx, alpha, MPFR_RNDN);
    mpfr_mul(rightFx, rightFx, alpha, MPFR_RNDN);
    //double l1 = mpfr_get_d(leftFx, MPFR_RNDN);
    //double l2 = mpfr_get_d(rightFx, MPFR_RNDN);
    //printf("leftval: %15.15f\n", l1);
    //printf("rightval: %15.15f\n", l2);
    mpfr_add(sum, sum, leftFx, MPFR_RNDN);
    mpfr_add(sum, sum, rightFx, MPFR_RNDN); flopCount += 6;
    //printf("sum[%d] = %15.15f\n", k, mpfr_get_d(sum, MPFR_RNDN));
    //printf("wk[%d] = %3.3e\n\n", k, mpfr_get_d(wkVals[k], MPFR_RNDN));
  }
  PetscReal prog = mpfr_get_d(sum, MPFR_RNDN);
  //printf("sum: %15.15f\n", prog);

  *integral = prog;

  mpfr_clears(alpha, beta, leftXk, rightXk, wk, leftFx, rightFx, pi2, kh, h, sum, lx, rx, NULL);
  for(k=0; k<nSide+1; ++k)
    mpfr_clears(wkVals[k], xkVals[k], NULL);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "ERFQuad"
PetscErrorCode ERFQuad(PetscErrorCode (*f)(mpfr_t*,mpfr_t*,void*), mpfr_t a, mpfr_t b, PetscInt prec, PetscInt nPts, PetscReal *integral, void *ctx)
{
  PetscErrorCode ierr;
  PetscInt k;
  PetscInt nSide;
  mpfr_t alpha, beta;
  mpfr_t leftXk, rightXk;
  mpfr_t wk;
  mpfr_t *wkVals;
  mpfr_t *xkVals;
  mpfr_t lx, rx;
  mpfr_t leftFx, rightFx;
  mpfr_t pi2;
  mpfr_t kh, h;
  mpfr_t sum;
  PetscInt flopCount;
  PetscFunctionBegin;
  flopCount = 0;

  if(nPts % 2 == 0) {
    nSide = nPts/2; flopCount++;
  }
  if(nPts % 2 == 1) {
    nSide = (nPts-1)/2; flopCount++;
  }

  /* initializations */
  mpfr_inits(alpha, beta, leftXk, rightXk, wk, leftFx, rightFx, pi2, kh, h, sum, lx, rx, NULL);
  wkVals = (mpfr_t*) malloc(sizeof(mpfr_t)*(nSide+1));
  xkVals = (mpfr_t*) malloc(sizeof(mpfr_t)*(nSide+1));
  for(k=0; k<nSide+1; ++k)
    mpfr_inits(xkVals[k], wkVals[k], NULL);

  // init alpha, beta
  mpfr_set(alpha, b, MPFR_RNDN);
  mpfr_sub(alpha, alpha, a, MPFR_RNDN);
  mpfr_mul_d(alpha, alpha, 0.5, MPFR_RNDN);
  mpfr_set(beta, b, MPFR_RNDN);
  mpfr_add(beta, beta, a, MPFR_RNDN);
  mpfr_mul_d(beta, beta, 0.5, MPFR_RNDN); flopCount += 4;

  /* determine right step size */
  ierr = FindStepSizeERF(prec, nSide, &h); // not including in flop count
  //mpfr_mul_d(h, h, 0.3, MPFR_RNDN);
  //printf("step size: %15.15f\n", mpfr_get_d(h, MPFR_RNDN));
  //printf("step size: %15.15f\n", mpfr_get_d(h, MPFR_RNDN));
  //mpfr_set_d(h, 0.04, MPFR_RNDN);

  /* calculate abscissas, weights */
  ierr = ERFxkwk(nSide+1, h, &xkVals, &wkVals);CHKERRQ(ierr);

  /* center term */
  mpfr_set_d(sum, 0.0, MPFR_RNDN);
  mpfr_set(lx, b, MPFR_RNDN);
  mpfr_add(lx, lx, a, MPFR_RNDN);
  mpfr_mul_d(lx, lx, .5, MPFR_RNDN);
  f(&lx, &sum, ctx);
  mpfr_mul(sum, sum, wkVals[0], MPFR_RNDN);
  mpfr_mul(sum, sum, alpha,     MPFR_RNDN); flopCount += 4;

  for(k=1; k <= nSide; ++k) {
    // adjust xk to [a,b]
    //mpfr_sub_d(lx, xkVals[k], 1.0, MPFR_RNDZ);
    mpfr_mul_d(lx, xkVals[k], -1.0, MPFR_RNDN);
    mpfr_mul(lx, lx, alpha, MPFR_RNDU);
    mpfr_add(lx, lx, beta, MPFR_RNDU);
    //mpfr_d_sub(rx, 1.0, xkVals[k], MPFR_RNDZ);
    mpfr_mul(rx, xkVals[k], alpha, MPFR_RNDD);
    mpfr_add(rx, rx, beta , MPFR_RNDD); flopCount += 6;

    double left = mpfr_get_d(lx, MPFR_RNDN);
    double right = mpfr_get_d(rx, MPFR_RNDN);
    double xk = mpfr_get_d(xkVals[k], MPFR_RNDN);
    //printf("left: %15.15f\n", left);
    //printf("right: %15.15f\n", right);
    //printf("\n");
    /* calculate function values */
    f(&lx, &leftFx , ctx);
    f(&rx, &rightFx, ctx);

    /* update sum */
    mpfr_mul(leftFx, leftFx, wkVals[k], MPFR_RNDN);
    mpfr_mul(rightFx, rightFx, wkVals[k], MPFR_RNDN);
    mpfr_mul(leftFx, leftFx, alpha, MPFR_RNDN);
    mpfr_mul(rightFx, rightFx, alpha, MPFR_RNDN);
    //double l1 = mpfr_get_d(leftFx, MPFR_RNDN);
    //double l2 = mpfr_get_d(rightFx, MPFR_RNDN);
    //printf("leftval: %15.15f\n", l1);
    //printf("rightval: %15.15f\n", l2);
    mpfr_add(sum, sum, leftFx, MPFR_RNDN);
    mpfr_add(sum, sum, rightFx, MPFR_RNDN); flopCount += 6;
    //printf("sum[%d] = %15.15f\n", k, mpfr_get_d(sum, MPFR_RNDN));
    //printf("wk[%d] = %3.3e\n\n", k, mpfr_get_d(wkVals[k], MPFR_RNDN));
  }
  PetscReal prog = mpfr_get_d(sum, MPFR_RNDN);
  //printf("sum: %15.15f\n", prog);

  *integral = prog;

  mpfr_clears(alpha, beta, leftXk, rightXk, wk, leftFx, rightFx, pi2, kh, h, sum, lx, rx, NULL);
  for(k=0; k<nSide+1; ++k)
    mpfr_clears(wkVals[k], xkVals[k], NULL);
  PetscFunctionReturn(0);
}

max_stars_repo_name: tom-klotz/ellipsoid-solvation

max_issues_repo_name: tom-klotz/ellipsoid-solvation

max_forks_repo_name: tom-klotz/ellipsoid-solvation


text:
#ifndef LIC_RUNTIME_TYPE
#define LIC_RUNTIME_TYPE

#include <string>
#include <gsl/gsl>

#include <format/Opcodes.h>

namespace lic
{

class RuntimeType
{
public:
    virtual ~RuntimeType();

    virtual std::string Name() const = 0;
    virtual std::string ToString(gsl::byte* data) const = 0;

    virtual void PerformUnaryOp(Opcode op, gsl::byte* data) const;
    virtual void PerformBinaryOp(Opcode op, gsl::byte* leftData, RuntimeType* rightType, gsl::byte* rightData) const;

protected:
    RuntimeType();
};

}

#endif // !LIC_RUNTIME_TYPE

ext: h

max_stars_repo_name: roberthusak/lic

max_issues_repo_name: roberthusak/lic

max_forks_repo_name: roberthusak/lic


text:
/* interpolation/interp_poly.c
 *
 * Copyright (C) 2001 DAN, HO-JIN
 * Copyright (C) 2013 Patrick Alken
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/* Modified for standalone use in polynomial directory, B.Gough 2001 */

#include <config.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_poly.h>

int
gsl_poly_dd_init (double dd[], const double xa[], const double ya[],
                  size_t size)
{
  size_t i, j;

  /* Newton's divided differences */

  dd[0] = ya[0];

  for (j = size - 1; j >= 1; j--)
    {
      dd[j] = (ya[j] - ya[j - 1]) / (xa[j] - xa[j - 1]);
    }

  for (i = 2; i < size; i++)
    {
      for (j = size - 1; j >= i; j--)
        {
          dd[j] = (dd[j] - dd[j - 1]) / (xa[j] - xa[j - i]);
        }
    }

  return GSL_SUCCESS;
}

int
gsl_poly_dd_taylor (double c[], double xp,
                    const double dd[], const double xa[], size_t size,
                    double w[])
{
  size_t i, j;

  for (i = 0; i < size; i++)
    {
      c[i] = 0.0;
      w[i] = 0.0;
    }

  w[size - 1] = 1.0;

  c[0] = dd[0];

  for (i = size - 1; i-- > 0;)
    {
      w[i] = -w[i + 1] * (xa[size - 2 - i] - xp);

      for (j = i + 1; j < size - 1; j++)
        {
          w[j] = w[j] - w[j + 1] * (xa[size - 2 - i] - xp);
        }

      for (j = i; j < size; j++)
        {
          c[j - i] += w[j] * dd[size - i - 1];
        }
    }

  return GSL_SUCCESS;
}

/*
gsl_poly_dd_hermite_init()
  Compute divided difference representation of data
for Hermite polynomial interpolation

Inputs: dd   - (output) array of size 2*size containing
               divided differences, dd[k] = f[z_0,z_1,...,z_k]
        za   - (output) array of size 2*size containing
               z values
        xa   - x data
        ya   - y data
        dya  - dy/dx data
        size - size of xa,ya,dya arrays

Return: success
*/

int
gsl_poly_dd_hermite_init (double dd[], double za[], const double xa[], const double ya[],
                          const double dya[], const size_t size)
{
  const size_t N = 2 * size;
  size_t i, j;

  /* Hermite divided differences */

  dd[0] = ya[0];

  /* compute: dd[j] = f[z_{j-1},z_j] for j \in [1,N-1] */
  for (j = 0; j < size; ++j)
    {
      za[2*j] = xa[j];
      za[2*j + 1] = xa[j];

      if (j != 0)
        {
          dd[2*j] = (ya[j] - ya[j - 1]) / (xa[j] - xa[j - 1]);
          dd[2*j - 1] = dya[j - 1];
        }
    }

  dd[N - 1] = dya[size - 1];

  for (i = 2; i < N; i++)
    {
      for (j = N - 1; j >= i; j--)
        {
          dd[j] = (dd[j] - dd[j - 1]) / (za[j] - za[j - i]);
        }
    }

  return GSL_SUCCESS;
} /* gsl_poly_dd_hermite_init() */

ext: c

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#pragma once

#include <gsl\gsl>
#include <winrt\Windows.Foundation.h>
#include <d3d11.h>
#include "DrawableGameComponent.h"
#include "MatrixHelper.h"
#include "DirectionalLight.h"

namespace Library
{
	class ProxyModel;
	class AnimationPlayer;
	class Model;
	class AnimationClip;
}

namespace Rendering
{
	class SkinnedModelMaterial;

	class AnimationDemo final : public Library::DrawableGameComponent
	{
	public:
		AnimationDemo(Library::Game& game, const std::shared_ptr<Library::Camera>& camera);
		AnimationDemo(const AnimationDemo&) = delete;
		AnimationDemo(AnimationDemo&&) = default;
		AnimationDemo& operator=(const AnimationDemo&) = default;
		AnimationDemo& operator=(AnimationDemo&&) = default;
		~AnimationDemo();

		bool ManualAdvanceEnabled() const;
		void SetManualAdvanceEnabled(bool enabled);
		void ToggleManualAdvance();

		const std::unique_ptr<Library::AnimationPlayer>& AnimationPlayer() const;
		void TogglePause();
		void RestartClip();

		void IncrementCurrentKeyframe();
		void DecrementCurrentKeyframe();

		bool InterpolationEnabled() const;
		void SetInterpolationEnabled(bool enabled);
		void ToggleInterpolation();

		float AmbientLightIntensity() const;
		void SetAmbientLightIntensity(float intensity);

		float DirectionalLightIntensity() const;
		void SetDirectionalLightIntensity(float intensity);

		const DirectX::XMFLOAT3& LightDirection() const;
		void RotateDirectionalLight(DirectX::XMFLOAT2 amount);

		virtual void Initialize() override;
		virtual void Update(const Library::GameTime& gameTime) override;
		virtual void Draw(const Library::GameTime& gameTime) override;

	private:
		inline static const float RotationRate{ DirectX::XM_PI };

		std::shared_ptr<SkinnedModelMaterial> mMaterial;
		DirectX::XMFLOAT4X4 mWorldMatrix{ Library::MatrixHelper::Identity };
		winrt::com_ptr<ID3D11Buffer> mVertexBuffer;
		winrt::com_ptr<ID3D11Buffer> mIndexBuffer;
		std::uint32_t mIndexCount{ 0 };
		Library::DirectionalLight mDirectionalLight;
		std::unique_ptr<Library::ProxyModel> mProxyModel;
		bool mUpdateMaterial{ true };
		bool mUpdateMaterialBoneTransforms{ true };
		std::unique_ptr<Library::AnimationPlayer> mAnimationPlayer;
		bool mManualAdvanceEnabled{ false };
		std::shared_ptr<Library::Model> mSkinnedModel;
	};
}


text:
//
//  View.hpp
//  PretendToWork
//
//  Created by tqtifnypmb on 14/12/2017.
//  Copyright © 2017 tqtifnypmb. All rights reserved.
//

#pragma once

#include "../rope/Range.h"
#include "../types.h"
#include "../editor/Editor.h"

#include <memory>
#include <utility>
#include <gsl/gsl>
#include <map>
#include <string>

namespace brick
{

class View {
public:

    using UpdateCb = std::function<void(View* view, const Editor::DeltaList&)>;

    View(size_t viewId, UpdateCb cb);
    View(size_t viewId, View* parent, UpdateCb cb);
    View(size_t viewId, const std::string& filePath, UpdateCb cb);

    ~View();

    void scroll(size_t begRow, size_t endRow);

    template<class Converter>
    void insert(gsl::span<const char> bytes);
    void insert(const detail::CodePointList& cplist);

    void erase();
    void undo();
    void select(Range sel);
    Range selection() const {
        return sel_;
    }

    void save();
    void save(const std::string& filePath);

    template <class Converter>
    std::map<size_t, std::string> region() {
        return region<Converter>(visibleRange_.first, visibleRange_.second);
    }

    template <class Converter>
    std::map<size_t, std::string> region(size_t begRow, size_t endRow);

    size_t viewId() const {
        return viewId_;
    }

    size_t viewSize() const {
        return viewSize_;
    }

    bool hasChildren() const {
        return !children_.empty();
    }

    void removeChild(const View* child) {
        auto found = std::find_if(children_.begin(), children_.end(), [child](const auto& v) { return v->viewId_ == child->viewId_; });
        if (found != children_.end()) {
            children_.erase(found);
        }
    }

    const std::vector<View*>& children() const {
        return children_;
    }

    const std::string filePath() const {
        return filePath_;
    }

    bool hasParent() const {
        return parent() != nullptr;
    }

    const View* parent() const {
        return parent_;
    }

private:
    void editor_sync_cb(const Editor::DeltaList& deltaList);
    void update(std::vector<View*>& src);

    std::map<size_t, detail::CodePointList> regionImpl(size_t begRow, size_t endRow);
    std::map<size_t, detail::CodePointList> regionImpl() {
        return regionImpl(visibleRange_.first, visibleRange_.second);
    }

    UpdateCb update_cb_;

    std::pair<size_t, size_t> visibleRange_;
    size_t viewId_;
    size_t viewSize_;
    View* parent_;
    Range sel_;
    std::vector<View*> children_;
    std::unique_ptr<Editor> editor_;
    std::string filePath_;
};

template<class Converter>
void View::insert(gsl::span<const char> bytes) {
    insert(Converter::encode(bytes));
}

template <class Converter>
std::map<size_t, std::string> View::region(size_t begRow, size_t endRow) {
    auto lines = regionImpl(begRow, endRow);
    std::map<size_t, std::string> ret;
    for (const auto& line : lines) {
        ret[line.first] = Converter::decode(line.second);
    }
    return ret;
}

}

max_stars_repo_name: tqtifnypmb/brick

max_issues_repo_name: tqtifnypmb/brick

max_forks_repo_name: tqtifnypmb/brick


text:
/*
 * Copyright 2020 Makani Technologies LLC
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#ifndef SIM_MATH_INTERP_H_
#define SIM_MATH_INTERP_H_

#include <gsl/gsl_interp.h>
#include <math.h>
#include <vector>

#include "common/c_math/vec3.h"
#include "common/macros.h"

namespace sim {

// A TimeseriesInterpolator performs linear interpolation at timepoints within
// the bounds of t, and 0th order extrapolation past the endpoints.
class TimeseriesInterpolator {
 public:
  explicit TimeseriesInterpolator(const std::vector<double> &t);
  ~TimeseriesInterpolator();

  double Interp(double t, const std::vector<double> &data);
  Vec3 InterpVec3(double t, const std::vector<Vec3> &data);

  double WrapInterp(double t, const std::vector<double> &data, double lower,
                    double upper);
  Vec3 WrapInterpVec3(double t, const std::vector<Vec3> &data, double lower,
                      double upper);

 private:
  double WrapInterpInternal(const std::vector<double> &values,
                            const std::vector<double> &weights, double lower,
                            double upper);

  struct InterpParams {
    InterpParams(const std::vector<int32_t> &indices__,
                 const std::vector<double> &weights__)
        : indices(indices__), weights(weights__) {}
    ~InterpParams() {}

    std::vector<int32_t> indices;
    std::vector<double> weights;
  };

  const std::vector<double> &t_;
  gsl_interp_accel *interp_accel_;

  InterpParams GetInterpParams(double t);

  DISALLOW_COPY_AND_ASSIGN(TimeseriesInterpolator);
};

}  // namespace sim

#endif  // SIM_MATH_INTERP_H_


text:
/*****************************************************************************
*
* Rokko: Integrated Interface for libraries of eigenvalue decomposition
*
* Copyright (C) 2012-2017 by Synge Todo <wistaria@phys.s.u-tokyo.ac.jp>
*
* Distributed under the Boost Software License, Version 1.0. (See accompanying
* file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
*
*****************************************************************************/

#include <cblas.h>
#include <lapacke.h>
#include <rokko/cmatrix.h>

int imin(int x, int y) { return (x < y) ? x : y; }
int imax(int x, int y) { return (x > y) ? x : y; }

int main(int argc, char** argv) {
  int m = 3;
  int n = 5;
  int r, lwork, i, j, info;
  double norm;
  double *s, *work;
  double **a, **u, **vt, **t;

  if (argc > 2) {
    m = atoi(argv[1]);
    n = atoi(argv[2]);
  }
  r = imin(m, n);

  /* generate matrix */
  a = alloc_dmatrix(m, n);
  for (j = 0; j < n; ++j) {
    for (i = 0; i < m; ++i) {
      mat_elem(a, i, j) = n - imax(i, j);
    }
  }
  printf("Matrix A: ");
  fprint_dmatrix(stdout, m, n, a);

  /* singular value decomposition */
  u = alloc_dmatrix(m, r);
  vt = alloc_dmatrix(r, n);
  s = alloc_dvector(r);
  t = alloc_dmatrix(m, n);
  cblas_dcopy(m * n, mat_ptr(a), 1, mat_ptr(t), 1);
  lwork = imax(3 * r + imax(m, n), 5 * r);
  work = alloc_dvector(lwork);
  info = LAPACKE_dgesvd_work(LAPACK_COL_MAJOR, 'S', 'S', m, n, mat_ptr(t), m,
                             vec_ptr(s), mat_ptr(u), m, mat_ptr(vt), r,
                             vec_ptr(work), lwork);
  if (info != 0) {
    fprintf(stderr, "Error: dgesvd fails\n");
    exit(255);
  }
  free_dmatrix(t);
  printf("Matrix U: ");
  fprint_dmatrix(stdout, m, r, u);
  printf("Singular values: ");
  fprint_dvector(stdout, r, s);
  printf("Matrix Vt: ");
  fprint_dmatrix(stdout, r, n, vt);

  /* orthogonality check */
  t = alloc_dmatrix(r, r);
  cblas_dgemm(CblasColMajor, CblasTrans, CblasNoTrans, r, r, m, 1, mat_ptr(u), m,
              mat_ptr(u), m, 0, mat_ptr(t), r);
  for (i = 0; i < r; ++i) mat_elem(t, i, i) -= 1;
  norm = LAPACKE_dlange(LAPACK_COL_MAJOR, 'F', r, r, mat_ptr(t), r);
  printf("|| U^t U - I || = %e\n", norm);
  if (norm > 1e-10) {
    fprintf(stderr, "Error: orthogonality check\n");
    exit(255);
  }

  cblas_dgemm(CblasColMajor, CblasNoTrans, CblasTrans, r, r, n, 1, mat_ptr(vt), r,
              mat_ptr(vt), r, 0, mat_ptr(t), r);
  for (i = 0; i < r; ++i) mat_elem(t, i, i) -= 1;
  norm = LAPACKE_dlange(LAPACK_COL_MAJOR, 'F', r, r, mat_ptr(t), r);
  printf("|| V^t V - I || = %e\n", norm);
  if (norm > 1e-10) {
    fprintf(stderr, "Error: orthogonality check\n");
    exit(255);
  }
  free_dmatrix(t);

  /* solution check */
  t = alloc_dmatrix(m, r);
  cblas_dgemm(CblasColMajor, CblasNoTrans, CblasTrans, m, r, n, 1, mat_ptr(a), m,
              mat_ptr(vt), r, 0, mat_ptr(t), m);
  free_dmatrix(vt);
  vt = alloc_dmatrix(r, r);
  cblas_dgemm(CblasColMajor, CblasTrans, CblasNoTrans, r, r, m, 1, mat_ptr(u), m,
              mat_ptr(t), m, 0, mat_ptr(vt), r);
  for (i = 0; i < r; ++i) mat_elem(vt, i, i) -= s[i];
  norm = LAPACKE_dlange(LAPACK_COL_MAJOR, 'F', r, r, mat_ptr(vt), r);
  printf("|| U^t A V - diag(S) || = %e\n", norm);
  if (norm > 1e-10) {
    fprintf(stderr, "Error: solution check\n");
    exit(255);
  }
  free_dmatrix(t);

  free_dmatrix(a);
  free_dmatrix(u);
  free_dmatrix(vt);
  free_dvector(s);
  return 0;
}

max_stars_repo_name: t-sakashita/rokko

max_issues_repo_name: t-sakashita/rokko

max_forks_repo_name: t-sakashita/rokko


text:
/* ============================================================ *
 * cosebi_demo.c						*
 * Martin Kilbinger, Liping Fu 2012.				*
 *  2013. Tomography.						*
 * Reference:							*
 * Schneider, Eifler, Krause 2010                               *
 * ============================================================ */

#include <stdio.h>
#include <stdlib.h>
#include "errorlist.h"
#include "io.h"
#include "maths.h"
#include "cosmo.h"
#include "lensing.h"
#include "nofz.h"
#include "decomp_eb.h"

#include <getopt.h>
#include <gsl/gsl_sf_gamma.h>

void write_tpmlog(int n, double thmin, double thmax, const double *c, error **err)
{
   FILE *F;
   double z, theta, tp, tm, dth;
   char name[128];

   sprintf(name, "T_%d_%.3f_%.3f.dat", n, thmin/arcmin, thmax/arcmin);

   F = fopen_err(name, "w", err); forwardError(*err, __LINE__,);

   fprintf(F, "# z=log(thmax/thmin) theta[arcmin] T_+^log T_-^log, (n, thmin, thmax = %d, %g, %g)\n",
	   n, thmin/arcmin, thmax/arcmin);

   dth = 1.01;
   for (theta=thmin; theta<=thmax; theta*=dth) {

      z = log(theta / thmin);

      tp = Tplog_c(z, c, n, err);
      forwardError(*err, __LINE__,);

      tm = Tmlog(z, c, n, err);
      forwardError(*err, __LINE__,);

      fprintf(F, "%f %g % g % g\n", z, theta/arcmin, tp, tm);
   }

   fclose(F);
}

void write_xi_pm_EB(int N, double thmin, double thmax, const double *c_cosebi,
		    const double *E_cos, const double *B_cos, error **err)
{
   double xi_pm_EB[4], theta, fth;
   const char name[128] = "xi_pm_EB.dat";
   int n;
   FILE *F;

   F = fopen_err(name, "w", err);       forwardError(*err, __LINE__,);
   fprintf(F, "# theta[arcmin] E+ E- B+ B-\n");

   fth = 1.1;
   for (theta=thmin; theta<=thmax; theta *= fth) {
      xipmEB(theta, thmin, thmax, c_cosebi, E_cos, B_cos, N, xi_pm_EB, err);
      forwardError(*err, __LINE__,);

      fprintf(F, "%12.3g  ", theta/arcmin);
      for (n=0; n<4; n++) fprintf(F, " % g", xi_pm_EB[n]);
      fprintf(F, "\n");
   }

   fclose(F);
}

double E_cosebi_data(const double *xip, const double *xim, const double *theta, int Nxi,
		     int n, double Psimin, double Psimax, const double *c_cosebi,
		     double *B_cosebi, error **err)
{
   double z, delta_theta, rp, rm, summand;
   int i;

   testErrorRetVA(Psimin < theta[0], mr_range, "COSEBI theta_min = %g' is smaller than xi theta_min = %g'",
		  *err, __LINE__, -1.0, Psimin / arcmin, theta[0] / arcmin);
   testErrorRetVA(Psimax > theta[Nxi-1], mr_range, "COSEBI theta_max = %g' is larger than xi theta_max = %g'",
		  *err, __LINE__, -1.0, Psimax / arcmin, theta[Nxi-1] / arcmin);

   rp = rm = 0.0;
   for (i=0; i<Nxi; i++) {

      if (theta[i] < Psimin) continue;
      if (theta[i] > Psimax) break;

      if (i > 0) {
	 delta_theta = theta[i] - theta[i-1];
      } else {
	 delta_theta = theta[1] - theta[0];
      }

      z = log(theta[i] / Psimin);

      summand  = xip[i];
      summand *= theta[i] * delta_theta;
      summand *= Tplog_c(z, c_cosebi, n, err);
      forwardError(*err, __LINE__, -1.0);
      rp      += summand;

      summand  = xim[i];
      summand *= theta[i] * delta_theta;
      summand *= Tmlog(z, c_cosebi, n, err);
      forwardError(*err, __LINE__, -1.0);
      rm      += summand;

      //printf("%g %g   %g %g\n", theta[i]/arcmin, z, xip[i], xim[i]);

   }

   if (B_cosebi != NULL) {
      *B_cosebi = 0.5 * (rp - rm);
   }

   return 0.5 * (rp + rm);
}

void usage(int ex)
{
   fprintf(stderr, "Calculates the COSEBIs from Schneider, Eifler, Krause 2010\n");
   fprintf(stderr, "Usage: coseb_demo [OPTIONS]\n");
   fprintf(stderr, "OPTIONS\n");
   fprintf(stderr, "  -N N       Maximum COSEBIs mode N\n");
   fprintf(stderr, "  -m THMIN   Minimum COSEBIs scale THMIN (in arcmin; default 1)\n");
   fprintf(stderr, "  -M THMAX   Maximum COSEBIs scale THMAX (in arcmin; default 400)\n");
   fprintf(stderr, "  -d XI      Data file XI (cosmo_pmc input ['xipm'] format; default none)\n");
   fprintf(stderr, "  -T         Write filter functions to files 'T_n_thmin_thmax.dat'\n");
   fprintf(stderr, "  -X         Write EB correlation functions to files 'xEB_n_thmin_thmax.dat'\n");
   fprintf(stderr, "  -P PATH    Path for coefficient files ('cosebi_tplog_rN_N_thmin_thmax';\n");
   fprintf(stderr, "              default $COSMOPMC/par_files/COSEBIs)\n");
   fprintf(stderr, "  -h         This message\n");

   if (ex >= 0) exit(ex);
}

int main(int argc, char* argv[])
{
   error *myerr = NULL, **err;
   double *E_cos, *B_cos;
   double Psimin, Psimax;
   int c, N, n, Nxi, write_T, write_X, i_bin, j_bin, Nzbin;
   cosmo_lens *model;
   char *xi_name, *path, *strp;
   FILE *F;
   datcov *dc;
   double *xip, *xim, *theta, *theta2, *c_cosebi;

   fprintf(stderr, "# cosebi_demo (M. Kilbinger, L. Fu, 2012)\n");

   err = &myerr;

   /* Command line arguments */
   Psimin  = 0.174;
   Psimax  = 200.0;
   N       = 10;
   xi_name = NULL;
   write_T = write_X = 0;
   path    = NULL;
   while (1) {

      static struct option long_options[] = {
	  {"", required_argument, 0, 'N'},
	  {"", required_argument, 0, 'd'},
	  {"", required_argument, 0, 'P'},
	  {0, 0, 0, 0}
	};

      int option_index = 0;

      c = getopt_long(argc, argv, "N:d:m:M:TXP:h", long_options, &option_index);
      switch (c) {
	 case 'N' :
	    N = atoi(optarg);
	    break;
	 case 'd' :
	    xi_name = optarg;
	    break;
	 case 'P' :
	    path = optarg;
	    break;
	 case 'm' :
       Psimin = atof(optarg);
	    break;
	 case 'M' :
	    Psimax = atof(optarg);
	    break;
	 case 'T' :
	    write_T = 1;
	    break;
	 case 'X' :
	    write_X = 1;
	    break;
	 case 'h' :
	    usage(0);
	 case -1 :
	    goto end_arg;
	 default :
	    usage(1);
      }

   }
 end_arg:

   if (N<0 || N>NMAX_COSEBI) {
      fprintf(stderr, "N=%d out of range [0; %d]", N, NMAX_COSEBI);
      exit(4);
   }

   if (path == NULL) {
      strp = getenv("COSMOPMC");
      if (strp == NULL) {
         strp = ".";
      }
      path = malloc_err(sizeof(char) * 1024, err);
      quitOnError(*err, __LINE__, stderr);
      sprintf(path, "%s/par_files/COSEBIs", strp);
   }

   E_cos = malloc_err(sizeof(double) * (N+1), err);      quitOnError(*err, __LINE__, stderr);
   B_cos = malloc_err(sizeof(double) * (N+1), err);      quitOnError(*err, __LINE__, stderr);

   if (xi_name == NULL) {

      /* COSEBIs from theoretical model */

      F = fopen_err("cosmo_lens.par", "r", err);
      quitOnError(*err, __LINE__, stderr);
      read_cosmological_parameters_lens(&model, F, err);
      quitOnError(*err, __LINE__, stderr);
      fclose(F);
      Nzbin = model->redshift->Nzbin;

      dump_param(model->cosmo, stdout);
      fprintf(stdout, "# n E_n^{ij} B_n^{ij}\n# n");
      for (i_bin=0; i_bin<Nzbin; i_bin++) {
         for (j_bin=i_bin; j_bin<Nzbin; j_bin++) {
            if (model->tomo==tomo_auto_only && i_bin!=j_bin) continue;
            if (model->tomo==tomo_cross_only && i_bin==j_bin) continue;
            fprintf(stdout, "        E_n^{%d%d}        B_n^{%d%d}", i_bin, j_bin, i_bin, j_bin);
         }
      }
      fprintf(stdout, "\n");

      for (n=1; n<=N; n++) {
         fprintf(stdout, "%3d", n);

         for (i_bin=0; i_bin<Nzbin; i_bin++) {
            for (j_bin=i_bin; j_bin<Nzbin; j_bin++) {
               if (model->tomo==tomo_auto_only && i_bin!=j_bin) continue;
               if (model->tomo==tomo_cross_only && i_bin==j_bin) continue;
               E_cos[n] = E_cosebi(model, n, Psimin*arcmin, Psimax*arcmin, i_bin, j_bin, path, B_cos+n, err);
               quitOnError(*err, __LINE__, stderr);
               fprintf(stdout, " % 15.6g % 15.6g", E_cos[n], B_cos[n]);
            }
         }

         fprintf(stdout, "\n");
      }

      c_cosebi = model->c_cosebi;

   } else {

      /* COSEBIs from data */
      c_cosebi = read_zeros_norm_cosebi_auto_check(Psimin*arcmin, Psimax*arcmin, path, err);
      quitOnError(*err, __LINE__, stderr);

      dc = malloc_err(sizeof(datcov), err);
      quitOnError(*err, __LINE__, stderr);
      dc->format = angle_center;
      /* Reading xi file */
      read_data_tomo(dc, xi_name, 0, second_order, err);
      quitOnError(*err, __LINE__, stderr);

      Nzbin    = dc->Nzbin;
      dc->type = xipm;

      /* Header */
      fprintf(stdout, "# n E_n^{ij} B_n^{ij}\n# n   ");
      for (i_bin=0; i_bin<Nzbin; i_bin++) {
         for (j_bin=i_bin; j_bin<Nzbin; j_bin++) {
            fprintf(stdout, "E_n^{%d%d}        B_n^{%d%d}          ", i_bin, j_bin, i_bin, j_bin);
         }
      }
      fprintf(stdout, "\n");

      for (n=1; n<=N; n++) {
	 fprintf(stdout, "%3d", n);

	 for (i_bin=0; i_bin<Nzbin; i_bin++) {
	    for (j_bin=i_bin; j_bin<Nzbin; j_bin++) {
	       datcov2xipm(dc, i_bin, j_bin, &xip, &xim, &theta, &theta2, &Nxi, err);
	       quitOnError(*err, __LINE__, stderr);

	       E_cos[n] = E_cosebi_data(xip, xim, theta, Nxi, n, Psimin*arcmin, Psimax*arcmin, c_cosebi, B_cos+n, err);
	       quitOnError(*err, __LINE__, stderr);
	       fprintf(stdout, "   % .9g % .9g", E_cos[n], B_cos[n]);
	    }
	 }
	 fprintf(stdout, "\n");

	 free(xip);
	 free(xim);
	 free(theta);
      }

      del_data_cov(&dc);

   }

   if (write_T) {

      /* Write filter functions to file */

      for (n=1; n<=N; n++) {
         write_tpmlog(n, Psimin*arcmin, Psimax*arcmin, c_cosebi, err);
         quitOnError(*err, __LINE__, stderr);
      }

   }

   if (write_X) {

      /* Calculate xipmEB according to SEK 10 eq (40) */

      write_xi_pm_EB(N, Psimin*arcmin, Psimax*arcmin, c_cosebi, E_cos, B_cos, err);

   }

   free(E_cos);
   free(B_cos);

   return 0;
}

max_stars_repo_name: danielgruen/ccv

max_issues_repo_name: danielgruen/ccv

max_forks_repo_name: danielgruen/ccv


text:
/* Copyright 2014-2018 The PySCF Developers. All Rights Reserved.

   Licensed under the Apache License, Version 2.0 (the "License");
    you may not use this file except in compliance with the License.
    You may obtain a copy of the License at

        http://www.apache.org/licenses/LICENSE-2.0

    Unless required by applicable law or agreed to in writing, software
    distributed under the License is distributed on an "AS IS" BASIS,
    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    See the License for the specific language governing permissions and
    limitations under the License.
 */

#include<stdio.h>
#include<stdlib.h>
#include <omp.h>
//#include <cblas.h>

/*
!
!
!Compute Y += A*X for CSR matrix A and dense vectors X,Y
! From scipy/sparse/sparsetools/csr.h
!
!
! Input Arguments:
!   I  n_row         - number of rows in A
!   I  n_col         - number of columns in A
!   I  Ap[n_row+1]   - row pointer
!   I  Aj[nnz(A)]    - column indices
!   T  Ax[nnz(A)]    - nonzeros
!   T  Xx[n_col]     - input vector
!
! Output Arguments:
!  T  Yx[n_row]     - output vector
!
! Note:
!   Output array Yx must be preallocated
!
!   Complexity: Linear.  Specifically O(nnz(A) + n_row)
*/
void scsr_matvec(int nrow, int ncol, int nnz, int *Ap, int *Aj,
    float *Ax, float *Xx, float *Yx)
{

  int i, jj;
  float sum = 0.0;

  # pragma omp parallel \
  shared (nrow, Yx, Ap, Ax, Xx, Aj) \
  private (i, jj, sum)
  {
    #pragma omp for
    for(i = 0; i < nrow; i++){
      sum = Yx[i];
      for(jj = Ap[i]; jj < Ap[i+1]; jj++){
        sum += Ax[jj] * Xx[Aj[jj]];
      }
      Yx[i] = sum;
    }
  }
}

void dcsr_matvec(int nrow, int ncol, int nnz, int *Ap, int *Aj,
    double *Ax, double *Xx, double *Yx)
{

  int i, jj;
  double sum = 0.0;

  # pragma omp parallel \
  shared (nrow, Yx, Ap, Ax, Xx, Aj) \
  private (i, jj, sum)
  {
    #pragma omp for
    for(i = 0; i < nrow; i++){
      sum = Yx[i];
      for(jj = Ap[i]; jj < Ap[i+1]; jj++){
        sum += Ax[jj] * Xx[Aj[jj]];
      }
      Yx[i] = sum;
    }
  }
}

/*
 * Compute Y += A*X for CSC matrix A and dense vectors X,Y
 * From scipy/sparse/sparsetools/csc.h
 *
 *
 * Input Arguments:
 *   I  n_row         - number of rows in A
 *   I  n_col         - number of columns in A
 *   I  Ap[n_row+1]   - column pointer
 *   I  Ai[nnz(A)]    - row indices
 *   T  Ax[n_col]     - nonzeros
 *   T  Xx[n_col]     - input vector
 *
 * Output Arguments:
 *   T  Yx[n_row]     - output vector
 *
 * Note:
 *   Output array Yx must be preallocated
 *
 *   Complexity: Linear.  Specifically O(nnz(A) + n_col)
 *
 */
void scsc_matvec(int n_row, int n_col, int nnz,
            int *Ap, int *Ai, float *Ax, float *Xx, float *Yx)
{
    int col_start, col_end, j, ii, i;

    for( j = 0; j < n_col; j++){
        col_start = Ap[j];
        col_end   = Ap[j+1];

        for( ii = col_start; ii < col_end; ii++){
            i    = Ai[ii];
            Yx[i] += Ax[ii] * Xx[j];
        }
    }
}

void dcsc_matvec(int n_row, int n_col, int nnz,
            int *Ap, int *Ai, double *Ax, double *Xx, double *Yx)
{
    int col_start, col_end, j, ii, i;

    for( j = 0; j < n_col; j++){
        col_start = Ap[j];
        col_end   = Ap[j+1];

        for( ii = col_start; ii < col_end; ii++){
            i    = Ai[ii];
            Yx[i] += Ax[ii] * Xx[j];
        }
    }
}

/*
 * Compute Y += A*X for CSC matrix A and dense block vectors X,Y
 * From scipy/sparse/sparsetools/csc.h
 *
 *
 * Input Arguments:
 *   I  n_row            - number of rows in A
 *   I  n_col            - number of columns in A
 *   I  n_vecs           - number of column vectors in X and Y
 *   I  Ap[n_row+1]      - row pointer
 *   I  Aj[nnz(A)]       - column indices
 *   T  Ax[nnz(A)]       - nonzeros
 *   T  Xx[n_col,n_vecs] - input vector
 *
 * Output Arguments:
 *   T  Yx[n_row,n_vecs] - output vector
 *
 * Note:
 *   Output array Yx must be preallocated
 *
*/
void scsc_matvecs(int n_row, int n_col, int n_vecs,
      int *Ap, int *Ai, float *Ax, float *Xx, float *Yx)
{
  int i, j, ii;
  printf("dcsc_matvecs not working in m_sparsetools.c\n");
  exit(1);
  /*
  # pragma omp parallel \
  shared (n_row, n_col, n_vecs, Ap, Ai, Ax, Xx, Yx) \
  private (i, ii, j)
  {
    #pragma omp for
    for( j = 0; j < n_col; j++){
      for( ii = Ap[j]; ii < Ap[j+1]; ii++){
        i = Ai[ii];
        //axpy(n_vecs, Ax[ii], Xx + (int)n_vecs * j, Yx + (int)n_vecs * i);
        cblas_saxpy (n_vecs, Ax[ii], &Xx[n_vecs*j], 1, &Yx[n_vecs*i], 1);
      }
    }
  //}
    */
}

void dcsc_matvecs(int n_row, int n_col, int n_vecs,
      int *Ap, int *Ai, double *Ax, double *Xx, double *Yx)
{
  int i, j, ii;
  printf("dcsc_matvecs not working in m_sparsetools.c\n");
  exit(1);
  /*
  # pragma omp parallel \
  shared (n_row, n_col, n_vecs, Ap, Ai, Ax, Xx, Yx) \
  private (i, ii, j)
  {
    #pragma omp for
    for( j = 0; j < n_col; j++){
      for( ii = Ap[j]; ii < Ap[j+1]; ii++){
        i = Ai[ii];
        //axpy(n_vecs, Ax[ii], Xx + (int)n_vecs * j, Yx + (int)n_vecs * i);
        cblas_daxpy (n_vecs, Ax[ii], &Xx[n_vecs*j], 1, &Yx[n_vecs*i], 1);
      }
    }
  //}
  */
}

max_stars_repo_name: robert-anderson/pyscf

max_issues_repo_name: robert-anderson/pyscf

max_forks_repo_name: robert-anderson/pyscf


text:
#ifndef __e23673_LIBRARIES_H
#define __e23673_LIBRARIES_H

#include <stdio.h>
#include <stdlib.h>
#include <fcntl.h>
#include <unistd.h>
#include <limits.h>

/* GNU Scientific library: */
#include <gsl/gsl_rng.h>

#include "parsing.h"
#include "dev_random.h"
#include "mersenne_twister.h"
#include "random_numbers.h"

#endif

max_stars_repo_name: gmagannaDevelop/TousAntiCovid

max_issues_repo_name: gmagannaDevelop/TousAntiCovid

max_forks_repo_name: gmagannaDevelop/TousAntiCovid


text:
/* gsl_histogram_calloc_range.c
 * Copyright (C) 2000  Simone Piccardi
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License as
 * published by the Free Software Foundation; either version 3 of the
 * License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License along
 * with this library; if not, write to the Free Software Foundation, Inc.,
 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */
/***************************************************************
 *
 * File gsl_histogram_calloc_range.c:
 * Routine to create a variable binning histogram providing
 * an input range vector. Need GSL library and header.
 * Do range check and allocate the histogram data.
 *
 * Author: S. Piccardi
 * Jan. 2000
 *
 ***************************************************************/
#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_histogram.h>

gsl_histogram *
gsl_histogram_calloc_range (size_t n, double *range)
{
  size_t i;
  gsl_histogram *h;

  /* check arguments */

  if (n == 0)
    {
      GSL_ERROR_VAL ("histogram length n must be positive integer",
                        GSL_EDOM, 0);
    }

  /* check ranges */

  for (i = 0; i < n; i++)
    {
      if (range[i] >= range[i + 1])
        {
          GSL_ERROR_VAL ("histogram bin extremes  must be "
                            "in increasing order", GSL_EDOM, 0);
        }
    }

  /* Allocate histogram  */

  h = (gsl_histogram *) malloc (sizeof (gsl_histogram));

  if (h == 0)
    {
      GSL_ERROR_VAL ("failed to allocate space for histogram struct",
                        GSL_ENOMEM, 0);
    }

  h->range = (double *) malloc ((n + 1) * sizeof (double));

  if (h->range == 0)
    {
      /* exception in constructor, avoid memory leak */
      free (h);
      GSL_ERROR_VAL ("failed to allocate space for histogram ranges",
                        GSL_ENOMEM, 0);
    }

  h->bin = (double *) malloc (n * sizeof (double));

  if (h->bin == 0)
    {
      /* exception in constructor, avoid memory leak */
      free (h->range);
      free (h);
      GSL_ERROR_VAL ("failed to allocate space for histogram bins",
                        GSL_ENOMEM, 0);
    }

  /* initialize ranges */

  for (i = 0; i <= n; i++)
    {
      h->range[i] = range[i];
    }

  /* clear contents */

  for (i = 0; i < n; i++)
    {
      h->bin[i] = 0;
    }

  h->n = n;

  return h;
}


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/time.h>
#include <math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

void read_matrix(int** index, int** matrix, double scaling, int N_kw, char* input_fileName)
{
    FILE *fp = fopen(input_fileName, "r");
    fscanf(fp, "%*[^\n]\n");

    for (int ii = 0; ii < N_kw; ii++)
        fscanf(fp, "%d,", &((*index)[ii]));
    fscanf(fp, "%*[^\n]\n");

    int tmp;
    for (int ii = 0; ii < N_kw; ii++)
    {
        for (int jj = 0; jj < N_kw; jj++)
        {
            fscanf(fp, "%d,", &tmp);
            (*matrix)[ii*N_kw + jj] = (int) (tmp * scaling);
        }
        fscanf(fp, "%*[^\n]\n");
    }
    fclose(fp);
}

void pad_matrix(int** matrix_padded, int** matrix, double lambda, int N_kw, int N_doc)
{
    // Initialising RNG
    const gsl_rng_type * T;
    gsl_rng *r1, *r2;
    gsl_rng_env_setup();
    r1 = gsl_rng_alloc(gsl_rng_default);
    r2 = gsl_rng_alloc (gsl_rng_taus);

    // perform padding
    int ii, jj;

    int max_resp_len = 0;
    for (ii = 0; ii < N_kw; ii++)
        if ((*matrix)[ii*N_kw + ii] > max_resp_len)
            max_resp_len = (*matrix)[ii*N_kw + ii];

    #pragma omp parallel for private(ii)
    for (ii = 0; ii < N_kw; ii++)
        (*matrix_padded)[ii*N_kw + ii] = (int) (lambda * (max_resp_len)) + gsl_rng_uniform_int(r2, max_resp_len+1);

    #pragma omp parallel for private(ii, jj)
    for (ii = 0; ii < N_kw; ii++)
    {
        for (jj = 0; jj < N_kw; jj++)
        {
            if (ii > jj)
            {
                int n1, n2, n3;
                n3 = 0;
                if ((*matrix_padded)[ii*N_kw + ii] > (*matrix)[ii*N_kw + ii])
                {
                    n1 = (*matrix)[ii*N_kw + jj];
                    n3 += gsl_ran_hypergeometric(r1, (*matrix_padded)[ii*N_kw + ii] - (*matrix)[ii*N_kw + ii], N_doc, (*matrix_padded)[jj*N_kw + jj]);
                }
                else
                    n1 = gsl_ran_hypergeometric(r1, (*matrix)[ii*N_kw + jj], (*matrix)[ii*N_kw + ii] - (*matrix)[ii*N_kw + jj], (*matrix_padded)[ii*N_kw + ii] - (*matrix)[ii*N_kw + jj]);

                if ((*matrix_padded)[jj*N_kw + jj] > (*matrix)[jj*N_kw + jj])
                {
                    n2 = n1;
                    n3 += gsl_ran_hypergeometric(r1, (*matrix_padded)[jj*N_kw + jj] - (*matrix)[jj*N_kw + jj], N_doc, (*matrix_padded)[ii*N_kw + ii] - (*matrix)[ii*N_kw + jj]);
                }
                else
                    n2 = gsl_ran_hypergeometric(r1, n1, (*matrix)[jj*N_kw + jj] - n1, (*matrix_padded)[jj*N_kw + jj]);

                (*matrix_padded)[ii*N_kw + jj] = n2 + n3;
                (*matrix_padded)[jj*N_kw + ii] = n2 + n3;
            }
        }
    }
    gsl_rng_free(r1);
    gsl_rng_free(r2);
}

void observe_matrix(gsl_matrix* matrix_obs, int** matrix_padded, int N_kw)
{
    // perform observed count generation
    for (int ii = 0; ii < N_kw; ii++)
        for (int jj = 0; jj < N_kw; jj++)
            gsl_matrix_set(matrix_obs, ii, jj, (double) ((*matrix_padded)[ii*N_kw + jj]));
}

void permutation_generation(int* idx1, int* idx2, int** permutation_tmp, int** permutation, int** permutation_inv, int N_kw, int N_obs)
{
    int count = 0;
    *idx1 = rand() % N_obs;
    *idx2 = -1;
    int idx_old = (*permutation)[*idx1];
    int idx_new = rand() % N_kw;

    (*permutation_tmp)[*idx1] = idx_new;
    if ((*permutation_inv)[idx_new] >= 0)
    {
        *idx2 = (*permutation_inv)[idx_new];
        (*permutation_tmp)[*idx2] = idx_old;
    }
}

max_stars_repo_name: RethinkingSSE/Attacks-on-SSE

max_issues_repo_name: RethinkingSSE/Attacks-on-SSE

max_forks_repo_name: RethinkingSSE/Attacks-on-SSE


text:
// @file ker2col_conv.c
//
//  \date Created on: Sep 24, 2017
//  \author Gopalakrishna Hegde
//
//   Description:
//
//
//

#include <stdbool.h>
#include <cblas.h>
#include <assert.h>
#include <stdlib.h>
#include <string.h>
#include "common_types.h"
#include "data_reshape.h"
#include "utils.h"
#include "conv_layers.h"

bool Kn2ColConvLayer(const float *in_data, const float *filters,
                         const float *bias, TensorDim in_dim,
                         TensorDim filt_dim, int stride, int pad, int group,
                         float *output) {
  // Currently we have limited support.
  assert(group == 1);
  assert((pad == 0) || (pad == filt_dim.w / 2));
  assert(in_dim.n == 1);
  assert(filt_dim.h == filt_dim.w);
  assert(stride == 1);

  // Output dimensions.
  TensorDim out_dim;
  out_dim.w = (in_dim.w + (pad + pad) - filt_dim.w) / stride + 1;
  out_dim.h = (in_dim.h + (pad + pad) - filt_dim.h) / stride + 1;
  out_dim.c = filt_dim.n;
  out_dim.n = in_dim.n;

  // Reshape filters in CHWN (ker_size x ker_size x no_in_ch x no_out_ch) format
  float *kkcm_filters = malloc(filt_dim.n * filt_dim.c * filt_dim.h *
                               filt_dim.w * sizeof(float));
  NCHW2HWCN(filters, filt_dim.n, filt_dim.c, filt_dim.h, filt_dim.w,
            kkcm_filters);

  // Just for convenience
  int H = in_dim.h;
  int W = in_dim.w;
  float alpha = 1.0;
  float beta = 0.0;

  // We need separate buffer because GEMM output will have width = H*W even
  // if there is no padding (pad = 0).
  float *gemm_output = malloc(out_dim.c * H * W * sizeof(float));
  float *nchw_gemm_output = malloc(out_dim.c * H * W * sizeof(float));
  // Prefill output buffer with bias if present else set to zero.
  if (bias) {
    for (int m = 0; m < out_dim.c; ++m) {
      for (int a = 0; a < out_dim.h * out_dim.w; ++a) {
        output[m * out_dim.h * out_dim.w + a] = bias[m];
      }
      // For batch size > 1
      for (int b = 1; b < out_dim.n; ++b) {
        memcpy(output + b * out_dim.c * out_dim.h * out_dim.w,
               output, out_dim.c * out_dim.h * out_dim.w * sizeof(float));
      }
    }
  } else {
    memset(output, 0, out_dim.n * out_dim.c * out_dim.h * out_dim.w *
           sizeof(float));
  }

  for (int kr = 0; kr < filt_dim.h; kr++) {
    int row_shift = pad - kr;
    for (int kc = 0; kc < filt_dim.w; kc++) {
      int group_no = kr * filt_dim.w + kc;
      int col_shift = pad - kc;
      // Matrix dimensions - A -> mxk B -> kxn  C --> mxn
      int m = in_dim.h * in_dim.w;
      int k = filt_dim.c;
      int n = filt_dim.n;

      // output is in H x W x C format.
      cblas_sgemm(CblasRowMajor, CblasTrans, CblasNoTrans,
                  m, n, k, alpha, in_data, m,
                  kkcm_filters + group_no * filt_dim.c * filt_dim.n, n, beta,
                  gemm_output, n);

      // convert to CxHxW format.
      // FIXME: this will be slow. Need to find other ways :(
      NHWC2NCHW(gemm_output, 1, filt_dim.n, H, W, nchw_gemm_output);

      for (int omap = 0; omap < filt_dim.n; omap++) {
        MatrixShiftAdd(output + omap * out_dim.h * out_dim.w,
                        out_dim.h, out_dim.w,
                        nchw_gemm_output + omap * H * W,
                        H, W, row_shift, col_shift);
      }
    }
  }
  free(kkcm_filters);
  free(gemm_output);
  free(nchw_gemm_output);
  return true;
}

max_stars_repo_name: gplhegde/convolution-flavors

max_issues_repo_name: chayitw/convolution-flavors

max_forks_repo_name: chayitw/convolution-flavors


text:
#ifndef LIC_OPCODES
#define LIC_OPCODES

#include <gsl/gsl>

namespace lic
{

constexpr uint8_t MethodHeaderFormatMask = 0x3;

enum class HeaderFormat : uint8_t
{
    Tiny = 0x2,
    Fat = 0x3
};

enum class Opcode : uint8_t
{
    Nop = 0x00,
    Break = 0x01,
    Ldarg_0 = 0x02,
    Ldarg_1 = 0x03,
    Ldarg_2 = 0x04,
    Ldarg_3 = 0x05,
    Ldloc_0 = 0x06,
    Ldloc_1 = 0x07,
    Ldloc_2 = 0x08,
    Ldloc_3 = 0x09,
    Stloc_0 = 0x0a,
    Stloc_1 = 0x0b,
    Stloc_2 = 0x0c,
    Stloc_3 = 0x0d,
    Ldarg_S = 0x0e,
    Ldarga_S = 0x0f,
    Starg_S = 0x10,
    Ldloc_S = 0x11,
    Ldloca_S = 0x12,
    Stloc_S = 0x13,
    Ldnull = 0x14,
    Ldc_I4_M1 = 0x15,
    Ldc_I4_0 = 0x16,
    Ldc_I4_1 = 0x17,
    Ldc_I4_2 = 0x18,
    Ldc_I4_3 = 0x19,
    Ldc_I4_4 = 0x1a,
    Ldc_I4_5 = 0x1b,
    Ldc_I4_6 = 0x1c,
    Ldc_I4_7 = 0x1d,
    Ldc_I4_8 = 0x1e,
    Ldc_I4_S = 0x1f,
    Ldc_I4 = 0x20,
    Ldc_I8 = 0x21,
    Ldc_R4 = 0x22,
    Ldc_R8 = 0x23,
    Dup = 0x25,
    Pop = 0x26,
    Jmp = 0x27,
    Call = 0x28,
    Calli = 0x29,
    Ret = 0x2a,
    Br_S = 0x2b,
    Brfalse_S = 0x2c,
    Brtrue_S = 0x2d,
    Beq_S = 0x2e,
    Bge_S = 0x2f,
    Bgt_S = 0x30,
    Ble_S = 0x31,
    Blt_S = 0x32,
    Bne_Un_S = 0x33,
    Bge_Un_S = 0x34,
    Bgt_Un_S = 0x35,
    Ble_Un_S = 0x36,
    Blt_Un_S = 0x37,
    Br = 0x38,
    Brfalse = 0x39,
    Brtrue = 0x3a,
    Beq = 0x3b,
    Bge = 0x3c,
    Bgt = 0x3d,
    Ble = 0x3e,
    Blt = 0x3f,
    Bne_Un = 0x40,
    Bge_Un = 0x41,
    Bgt_Un = 0x42,
    Ble_Un = 0x43,
    Blt_Un = 0x44,
    Switch = 0x45,
    Ldind_I1 = 0x46,
    Ldind_U1 = 0x47,
    Ldind_I2 = 0x48,
    Ldind_U2 = 0x49,
    Ldind_I4 = 0x4a,
    Ldind_U4 = 0x4b,
    Ldind_I8 = 0x4c,
    Ldind_I = 0x4d,
    Ldind_R4 = 0x4e,
    Ldind_R8 = 0x4f,
    Ldind_Ref = 0x50,
    Stind_Ref = 0x51,
    Stind_I1 = 0x52,
    Stind_I2 = 0x53,
    Stind_I4 = 0x54,
    Stind_I8 = 0x55,
    Stind_R4 = 0x56,
    Stind_R8 = 0x57,
    Add = 0x58,
    Sub = 0x59,
    Mul = 0x5a,
    Div = 0x5b,
    Div_Un = 0x5c,
    Rem = 0x5d,
    Rem_Un = 0x5e,
    And = 0x5f,
    Or = 0x60,
    Xor = 0x61,
    Shl = 0x62,
    Shr = 0x63,
    Shr_Un = 0x64,
    Neg = 0x65,
    Not = 0x66,
    Conv_I1 = 0x67,
    Conv_I2 = 0x68,
    Conv_I4 = 0x69,
    Conv_I8 = 0x6a,
    Conv_R4 = 0x6b,
    Conv_R8 = 0x6c,
    Conv_U4 = 0x6d,
    Conv_U8 = 0x6e,
    Callvirt = 0x6f,
    Cpobj = 0x70,
    Ldobj = 0x71,
    Ldstr = 0x72,
    Newobj = 0x73,
    Castclass = 0x74,
    Isinst = 0x75,
    Conv_R_Un = 0x76,
    Unbox = 0x79,
    Throw = 0x7a,
    Ldfld = 0x7b,
    Ldflda = 0x7c,
    Stfld = 0x7d,
    Ldsfld = 0x7e,
    Ldsflda = 0x7f,
    Stsfld = 0x80,
    Stobj = 0x81,
    Conv_Ovf_I1_Un = 0x82,
    Conv_Ovf_I2_Un = 0x83,
    Conv_Ovf_I4_Un = 0x84,
    Conv_Ovf_I8_Un = 0x85,
    Conv_Ovf_U1_Un = 0x86,
    Conv_Ovf_U2_Un = 0x87,
    Conv_Ovf_U4_Un = 0x88,
    Conv_Ovf_U8_Un = 0x89,
    Conv_Ovf_I_Un = 0x8a,
    Conv_Ovf_U_Un = 0x8b,
    Box = 0x8c,
    Newarr = 0x8d,
    Ldlen = 0x8e,
    Ldelema = 0x8f,
    Ldelem_I1 = 0x90,
    Ldelem_U1 = 0x91,
    Ldelem_I2 = 0x92,
    Ldelem_U2 = 0x93,
    Ldelem_I4 = 0x94,
    Ldelem_U4 = 0x95,
    Ldelem_I8 = 0x96,
    Ldelem_I = 0x97,
    Ldelem_R4 = 0x98,
    Ldelem_R8 = 0x99,
    Ldelem_Ref = 0x9a,
    Stelem_I = 0x9b,
    Stelem_I1 = 0x9c,
    Stelem_I2 = 0x9d,
    Stelem_I4 = 0x9e,
    Stelem_I8 = 0x9f,
    Stelem_R4 = 0xa0,
    Stelem_R8 = 0xa1,
    Stelem_Ref = 0xa2,
    Conv_Ovf_I1 = 0xb3,
    Conv_Ovf_U1 = 0xb4,
    Conv_Ovf_I2 = 0xb5,
    Conv_Ovf_U2 = 0xb6,
    Conv_Ovf_I4 = 0xb7,
    Conv_Ovf_U4 = 0xb8,
    Conv_Ovf_I8 = 0xb9,
    Conv_Ovf_U8 = 0xba,
    Refanyval = 0xc2,
    Ckfinite = 0xc3,
    Mkrefany = 0xc6,
    Ldtoken = 0xd0,
    Conv_U2 = 0xd1,
    Conv_U1 = 0xd2,
    Conv_I = 0xd3,
    Conv_Ovf_I = 0xd4,
    Conv_Ovf_U = 0xd5,
    Add_Ovf = 0xd6,
    Add_Ovf_Un = 0xd7,
    Mul_Ovf = 0xd8,
    Mul_Ovf_Un = 0xd9,
    Sub_Ovf = 0xda,
    Sub_Ovf_Un = 0xdb,
    Endfinally = 0xdc,
    Leave = 0xdd,
    Leave_S = 0xde,
    Stind_I = 0xdf,
    Conv_U = 0xe0,

    LongOpcodePrefix = 0xfe
};

enum class LongOpcodeSuffix: uint8_t
{
    Arglist = 0x00,
    Ceq = 0x01,
    Cgt = 0x02,
    Cgt_Un = 0x03,
    Clt = 0x04,
    Clt_Un = 0x05,
    Ldftn = 0x06,
    Ldvirtftn = 0x07,
    Ldarg = 0x09,
    Ldarga = 0x0a,
    Starg = 0x0b,
    Ldloc = 0x0c,
    Ldloca = 0x0d,
    Stloc = 0x0e,
    Localloc = 0x0f,
    Endfilter = 0x11,
    Unaligned_ = 0x12,
    Volatile_ = 0x13,
    Tail_ = 0x14,
    Initobj = 0x15,
    Cpblk = 0x17,
    Initblk = 0x18,
    Rethrow = 0x1a,
    Sizeof = 0x1c,
    Refanytype = 0x1d,

    InvalidOpcode = 0xff
};

const char* OpcodeName(Opcode first, LongOpcodeSuffix second = LongOpcodeSuffix::InvalidOpcode);

}

#endif // !LIC_OPCODES


text:
#ifndef VISCOUSVORTEXDOMAINSOLVER_H
#define VISCOUSVORTEXDOMAINSOLVER_H

#include <stdio.h>
#include <stdlib.h>
#include <common_functions.h>
#include <string.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>
#include <vector>
#include <algorithm>
//#include <parallel/algorithm>
#include "pstl/execution"
#include "pstl/algorithm"
#include "tbb/parallel_for_each.h"
#include <mutex>

#include <omp.h>

#include <streamline.h>
#include <vortex.h>
#include <vortexmotion.h>

//ToDo: add resumption of calculations mechanism
class ViscousVortexDomainSolver
{
    int NumberOfPanels;
    int MaxVortexes;
    int Niterations;

    double nu;
    double tau;
    double rho;
    double vx_inf;
    double vy_inf;

    double  *PanelLength;
    double  **PanelNodes,
            **PanelMids,
            **PanelNorms,
            **PanelTaus,
            **DeployPoints;

    gsl_matrix *OriginalVortexGenerationMatrix;
    gsl_matrix *VortexGenerationMatrix;
    gsl_matrix *VortexGenerationMatrix_Inversed;
    gsl_vector *VelocityProjections;
    gsl_vector *NewVorticities;
    gsl_permutation *Permutation;

    std::vector <Vortex> Flow;
    std::vector <VortexMotion> FlowEvolution;
    std::mutex SyncFlowEvolution;
    std::vector <Vortex> InBodyVortexes;

    std::vector <ControlPointOfFlow> StreamLine1;
    std::vector <ControlPointOfFlow> StreamLine2;
    void StepInTimeForStreamLines();
    void UpdateControlPointsOfStreamLine(std::vector<ControlPointOfFlow> &StreamLine);
    void MoveSingleControlPoint(ControlPointOfFlow &UpdatingPoint);

    int isVortexInBody(double x, double y); //ToDo: update to be virtual

    void DivideProfileToPanels();
    void CompletingGeneratingMatrix();
    void ComputeAssociatedVortexes();
    void AppendAssociatedVortexesToFlow();
    void UpdateVotexPositions();
    void ComputeVortexMotion(Vortex &Considered);
    double UpdateEpsilon(Vortex Considered, double InitialEpsilon);
    int LeakageControl(Vortex Checking);

public:

    ViscousVortexDomainSolver();
    ~ViscousVortexDomainSolver();

    void Solve();

    double xVelocityAt(double x, double y);
    double yVelocityAt(double x, double y);

    void Output_ParaView_Field(const char *FileName);
    void Output_ParaView_Field_with_exclusions(const char *FileName);
    void Output_ParaView_AllVortexes(const char *FileName);
    void Output_StreamLines(const char *FileName);
};

#endif // VISCOUSVORTEXDOMAINSOLVER_H


text:
/* ieee-utils/fp-darwin86.c
 *
 * Copyright (C) 2006 Erik Schnetter
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_errno.h>

/* Here is the dirty part. Set up your 387 through the control word
 * (cw) register.
 *
 *     15-13    12  11-10  9-8     7-6     5    4    3    2    1    0
 * | reserved | IC | RC  | PC | reserved | PM | UM | OM | ZM | DM | IM
 *
 * IM: Invalid operation mask
 * DM: Denormalized operand mask
 * ZM: Zero-divide mask
 * OM: Overflow mask
 * UM: Underflow mask
 * PM: Precision (inexact result) mask
 *
 * Mask bit is 1 means no interrupt.
 *
 * PC: Precision control
 * 11 - round to extended precision
 * 10 - round to double precision
 * 00 - round to single precision
 *
 * RC: Rounding control
 * 00 - rounding to nearest
 * 01 - rounding down (toward - infinity)
 * 10 - rounding up (toward + infinity)
 * 11 - rounding toward zero
 *
 * IC: Infinity control
 * That is for 8087 and 80287 only.
 *
 * The hardware default is 0x037f which we use.
 */

/* masking of interrupts */
#define _FPU_MASK_IM  0x01
#define _FPU_MASK_DM  0x02
#define _FPU_MASK_ZM  0x04
#define _FPU_MASK_OM  0x08
#define _FPU_MASK_UM  0x10
#define _FPU_MASK_PM  0x20

/* precision control */
#define _FPU_EXTENDED 0x300	/* libm requires double extended precision.  */
#define _FPU_DOUBLE   0x200
#define _FPU_SINGLE   0x0

/* rounding control */
#define _FPU_RC_NEAREST 0x0    /* RECOMMENDED */
#define _FPU_RC_DOWN    0x400
#define _FPU_RC_UP      0x800
#define _FPU_RC_ZERO    0xC00

#define _FPU_RESERVED 0xF0C0  /* Reserved bits in cw */

/* The fdlibm code requires strict IEEE double precision arithmetic,
   and no interrupts for exceptions, rounding to nearest.  */

#define _FPU_DEFAULT  0x037f

/* IEEE:  same as above.  */
#define _FPU_IEEE     0x037f

/* Type of the control word.  */
typedef unsigned int fpu_control_t __attribute__ ((__mode__ (__HI__)));

/* Macros for accessing the hardware control word.

   Note that the use of these macros is no sufficient anymore with
   recent hardware.  Some floating point operations are executed in
   the SSE/SSE2 engines which have their own control and status register.  */
#define _FPU_GETCW(cw) __asm__ __volatile__ ("fnstcw %0" : "=m" (*&cw))
#define _FPU_SETCW(cw) __asm__ __volatile__ ("fldcw %0" : : "m" (*&cw))

/* Default control word set at startup.  */
extern fpu_control_t __fpu_control;

#define _FPU_GETMXCSR(cw_sse) asm volatile ("stmxcsr %0" : "=m" (cw_sse))
#define _FPU_SETMXCSR(cw_sse) asm volatile ("ldmxcsr %0" : : "m" (cw_sse))

int
gsl_ieee_set_mode (int precision, int rounding, int exception_mask)
{
  fpu_control_t mode, mode_sse;

  _FPU_GETCW (mode) ;
  mode &= _FPU_RESERVED ;

  switch (precision)
    {
    case GSL_IEEE_SINGLE_PRECISION:
      mode |= _FPU_SINGLE ;
      break ;
    case GSL_IEEE_DOUBLE_PRECISION:
      mode |= _FPU_DOUBLE ;
      break ;
    case GSL_IEEE_EXTENDED_PRECISION:
      mode |= _FPU_EXTENDED ;
      break ;
    default:
      mode |= _FPU_EXTENDED ;
    }

  switch (rounding)
    {
    case GSL_IEEE_ROUND_TO_NEAREST:
      mode |= _FPU_RC_NEAREST ;
      break ;
    case GSL_IEEE_ROUND_DOWN:
      mode |= _FPU_RC_DOWN ;
      break ;
    case GSL_IEEE_ROUND_UP:
      mode |= _FPU_RC_UP ;
      break ;
    case GSL_IEEE_ROUND_TO_ZERO:
      mode |= _FPU_RC_ZERO ;
      break ;
    default:
      mode |= _FPU_RC_NEAREST ;
    }

  if (exception_mask & GSL_IEEE_MASK_INVALID)
    mode |= _FPU_MASK_IM ;

  if (exception_mask & GSL_IEEE_MASK_DENORMALIZED)
    mode |= _FPU_MASK_DM ;

  if (exception_mask & GSL_IEEE_MASK_DIVISION_BY_ZERO)
    mode |= _FPU_MASK_ZM ;

  if (exception_mask & GSL_IEEE_MASK_OVERFLOW)
    mode |= _FPU_MASK_OM ;

  if (exception_mask & GSL_IEEE_MASK_UNDERFLOW)
    mode |= _FPU_MASK_UM ;

  if (exception_mask & GSL_IEEE_TRAP_INEXACT)
    {
      mode &= ~ _FPU_MASK_PM ;
    }
  else
    {
      mode |= _FPU_MASK_PM ;
    }

  _FPU_SETCW (mode) ;

  _FPU_GETMXCSR (mode_sse) ;
  mode_sse &= 0xFFFF0000 ;

  if (exception_mask & GSL_IEEE_MASK_INVALID)
    mode_sse |= _FPU_MASK_IM << 7 ;

  if (exception_mask & GSL_IEEE_MASK_DENORMALIZED)
    mode_sse |= _FPU_MASK_DM << 7 ;

  if (exception_mask & GSL_IEEE_MASK_DIVISION_BY_ZERO)
    mode_sse |= _FPU_MASK_ZM << 7 ;

  if (exception_mask & GSL_IEEE_MASK_OVERFLOW)
    mode_sse |= _FPU_MASK_OM << 7 ;

  if (exception_mask & GSL_IEEE_MASK_UNDERFLOW)
    mode_sse |= _FPU_MASK_UM << 7 ;

  if (exception_mask & GSL_IEEE_TRAP_INEXACT)
    {
      mode_sse &= ~ _FPU_MASK_PM << 7 ;
    }
  else
    {
      mode_sse |= _FPU_MASK_PM << 7 ;
    }

  _FPU_SETMXCSR (mode_sse) ;

  return GSL_SUCCESS ;
}


text:
//=========================================
// Main driver for the proper time integral
//=========================================

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include <builtin_types.h>
#include "intT.h"

double* transformT(double t, void* p, int* order)
//Transforms the integral [0,infinity) to [0,1]
//x=(1-t)/t
{
	struct Wparams *params = (struct Wparams *) p;
	double* tfunc;
	double x = (1.0-t)/t;
	int i;

	tfunc = Tfunc(x,p, order);
	//loop over groups of worldlines
	for(i = 0; i < params->ng; i++)
	{
		tfunc[i] = (tfunc[i]/t)/t;

	}
	return tfunc;

}

void testval(double *val, int n, double t)
//used for outputting debugging information
{
  int i;
  double avg = 0.0;
  for(i = 0; i < n; i++)
  {
	if(isnan(val[i])) printf("for t=%f val[%d]=nan\n",t, i);
	avg += val[i];
  }
  avg /= n;

  printf("mean for Tt=%f, %f\n",t,avg);
}

void Tintout(double *val, int n, double rho, int fermion)
//Produces output about the T-integral
{
	int i;
	double avg = 0.0, SEM = 0.0;
  	for(i = 0; i < n; i++)
  	{
		if(isnan(val[i])) printf("for rho=%f val[%d]=nan\n",rho, i);
		avg += val[i];
  	}
  	avg /= n;
	for(i = 0; i < n; i++)
 	{
		SEM += (val[i]-avg)*(val[i]-avg);
  	}
  	SEM = sqrt(SEM/((double)n*((double)n-1.0)));
  	//rho Leff SEM-Leff

	if(fermion == 1)
  		printf("Leffvsrho: %e %e %e\n", rho, avg/(4.0*pi), SEM/(4.0*pi));
	else
  		printf("Leffvsrho: %e %e %e\n", rho, avg/(-2.0*pi), SEM/(2.0*pi));
}

int* shuffleWL(int Nl)
//Returns a vector representing a shuffled order of worldlines
{
  gsl_rng * r;
  int i;

  // select random number generator
  r = gsl_rng_alloc (gsl_rng_mt19937);

  int* order=(int*)malloc(Nl*sizeof(*order));

  for (i = 0; i < Nl; i++)
  {
  	order[i] = i;
  }

  //gsl_ran_shuffle (r, order, Nl, sizeof (int));
  free(r);

  return order;
}

double* integrateT(int n, void * p)
//Performs the proper time integral using Simpson's Method
//n = number of points to use for Simpson's method
//p = Wparams parameters for the integrand
{
        double* val;    //value of the integral for each group
	double* trt;    //Integrands evaluated at three points
	double* trtphalf;
	double* trtp1;
	int* order;     //Stores the order of worldlines
	struct Wparams *params = (struct Wparams *) p;
	int inx, ig, i;
	val = (double *)calloc(params->ng,sizeof(*val));
	double h = 1.0/((double) n);   //
	double t;
	if(!val) printf("integrate(): malloc failed\n");
	//initialize val to zero
	for(ig = 0; ig < params->ng; ig++)
	{
		val[ig] = 0.0;
	}
	//main integration loop
	for(inx = 0; inx < n; inx++)
	{
		printf("T integrate inx=%d\n",inx);
		//shuffle worldlines
		order=shuffleWL(params->Nl);
		t=h*((double)inx+1.0);
		trt=transformT(t,p,order);
		trtphalf=transformT(h*((double)inx + 0.5), p, order);
		trtp1=transformT(h*((double) inx+1.0), p, order);

		//A kludge to ignore possible infinities or nans
		for(ig = 0; ig < params->ng; ig++)
		{
			if(isinf(trt[ig]))
			{
				printf("Warning: trt[%d]=%f.  Setting to zero.\n",ig,trt[ig]);
				trt[ig]=0.0;
			}
			if(isinf(trtphalf[ig]))
			{
				printf("Warning: trtphalf[%d]=%f.  Setting to zero.\n",ig,trtphalf[ig]);
				trtphalf[ig]=0.0;
			}
			if(isinf(trtp1[ig]))
			{
				printf("Warning: trtp1[%d]=%f.  Setting to zero.\n",ig,trtp1[ig]);
				trtp1[ig]=0.0;
			}

        		val[ig]+=h/6.0*(
				trt[ig]
				+ 4.0*trtphalf[ig]
				+trtp1[ig]);
			if(isinf(val[ig])){
				printf("trt[%d]=%e\n",ig, trt[ig]);
				printf("trtphalf[%d]=%e\n",ig,trtphalf[ig]);
				printf("trtp1[%d]=%e\n",ig,trtp1[ig]);
				printf("val[%d]=%e\n",ig,val[ig]);

			}
			//printf("val[%d]=%e\n",ig,val[ig]);
		}
		testval(val, params->ng, t);
		free(order);

	}
	Tintout(val, params->ng, params->xcm.x, params->fermion);
	free(trt);
	free(trtphalf);
	free(trtp1);

        return val;

}

max_stars_repo_name: QEDan/Worldline-Flux-Tube-Lattice

max_issues_repo_name: QEDan/Worldline-Flux-Tube-Lattice

max_forks_repo_name: QEDan/Worldline-Flux-Tube-Lattice


text:
#ifndef ARRUS_CORE_DEVICES_US4R_US4ROUTPUTBUFFER_H
#define ARRUS_CORE_DEVICES_US4R_US4ROUTPUTBUFFER_H

#include <mutex>
#include <condition_variable>
#include <gsl/span>
#include <chrono>
#include <iostream>

#include "arrus/core/api/common/types.h"
#include "arrus/core/api/common/exceptions.h"
#include "arrus/common/asserts.h"
#include "arrus/common/format.h"
#include "arrus/core/common/logging.h"
#include "arrus/core/api/framework/DataBuffer.h"

namespace arrus::devices {

using ::arrus::framework::Buffer;
using ::arrus::framework::BufferElement;

class Us4ROutputBuffer;

/**
 * Buffer element owns the data arrrays, which then are returned to user.
 */
class Us4ROutputBufferElement : public BufferElement {
public:
    using AccumulatorType = uint16;

    using SharedHandle = std::shared_ptr<Us4ROutputBufferElement>;

    Us4ROutputBufferElement(int16 *address, size_t size,
                            const framework::NdArray::Shape &elementShape,
                            const framework::NdArray::DataType elementDataType,
                            AccumulatorType filledAccumulator,
                            size_t position)
        : data(address, elementShape, elementDataType,
               DeviceId(DeviceType::Us4R, 0)), size(size),
          filledAccumulator(filledAccumulator),
          position(position)
          {}

    void release() override {
        std::unique_lock<std::mutex> guard(mutex);
        this->accumulator = 0;
        releaseFunction();
    }

    int16 *getAddress() {
        validateState();
        return data.get<int16>();
    }

    framework::NdArray &getData() override {
        validateState();
        return data;
    }

    size_t getSize() override {
        return size;
    }

    size_t getPosition() override {
        return position;
    }

    void registerReleaseFunction(std::function<void()> &func) {
        releaseFunction = func;
    }

    [[nodiscard]] bool isElementReady() {
        std::unique_lock<std::mutex> guard(mutex);
        return accumulator == filledAccumulator;
    }

    void signal(Ordinal n) {
        std::unique_lock<std::mutex> guard(mutex);
        AccumulatorType us4oemPattern = 1ul << n;
        if((accumulator & us4oemPattern) != 0) {
            throw IllegalStateException("Detected data overflow, buffer is in invalid state.");
        }
        accumulator |= us4oemPattern;
    }

    void resetState() {
        accumulator = 0;
    }

    void markAsInvalid() {
        this->isInvalid = true;
    }

    void validateState() {
        if(this->isInvalid) {
            throw ::arrus::IllegalStateException(
                "The buffer is in invalid state "
                "(probably some data transfer overflow happened).");
        }
    }

private:
    std::mutex mutex;
    framework::NdArray data;
    size_t size;
    AccumulatorType accumulator;
    /** A pattern of the filled accumulator, which indicates that the
     * whole element is ready. */
    AccumulatorType filledAccumulator;
    std::function<void()> releaseFunction;
    bool isInvalid{false};
    size_t position;
};

/**
 * Us4R system's output circular FIFO buffer.
 *
 * The buffer has the following relationships:
 * - buffer contains **elements**
 * - the **element** is filled by many us4oems (with given ordinal)
 *
 * A single element is the output of a single data transfer (the result of running a complete sequence once).
 *
 * The state of each buffer element is determined by the field accumulators:
 * - accumulators[element] == 0 means that the buffer element was processed and is ready for new data from the producer.
 * - accumulators[element] > 0 && accumulators[element] != filledAccumulator means that the buffer element is partially
 *   confirmed by some of us4oems
 * - accumulators[element] == filledAccumulator means that the buffer element is ready to be processed by a consumer.
 *
 * The assumption is here that each element of the buffer has the same size (and the same us4oem offsets).
 */
class Us4ROutputBuffer : public framework::DataBuffer {
public:
    static constexpr size_t DATA_ALIGNMENT = 4096;
    using DataType = int16;

    /**
     * Buffer's constructor.
     *
     * @param us4oemOutputSizes number of bytes to allocate for each of the
     *  us4oem output. That is, the i-th value describes how many bytes will
     *  be written by i-th us4oem to generate a single buffer element.
     */
    Us4ROutputBuffer(const std::vector<size_t> &us4oemOutputSizes,
                     const framework::NdArray::Shape &elementShape,
                     const framework::NdArray::DataType elementDataType,
                     const unsigned nElements)
        : elementSize(0) {
        ARRUS_REQUIRES_TRUE(us4oemOutputSizes.size() <= 16,
                            "Currently Us4R data buffer supports up to 16 us4oem modules.");

        size_t nus4oems = us4oemOutputSizes.size();
        Us4ROutputBufferElement::AccumulatorType filledAccumulator((1ul << nus4oems) - 1);
        // Calculate us4oem write offsets for each buffer element.
        size_t us4oemOffset = 0;
        Ordinal us4oemOrdinal = 0;
        for(auto s : us4oemOutputSizes) {
            us4oemOffsets.emplace_back(us4oemOffset);
            us4oemOffset += s;
            if(s == 0) {
                // We should not expect any response from modules, which do not acquire any data.
                filledAccumulator &= ~(1ul << us4oemOrdinal);
            }
            ++us4oemOrdinal;
        }
        elementSize = us4oemOffset;
        // Allocate buffer with an appropriate size.
        dataBuffer = reinterpret_cast<DataType *>(operator new[](elementSize * nElements,
                                                                 std::align_val_t(DATA_ALIGNMENT)));
        getDefaultLogger()->log(LogSeverity::DEBUG,
                                ::arrus::format(
                                    "Allocated {} ({}, {}) bytes of memory, address: {}",
                                    elementSize * nElements, elementSize, nElements, (size_t) dataBuffer));

        for(unsigned i = 0; i < nElements; ++i) {
            auto elementAddress = reinterpret_cast<DataType *>(reinterpret_cast<int8 *>(dataBuffer) + i * elementSize);
            elements.push_back(std::make_shared<Us4ROutputBufferElement>(elementAddress, elementSize,
                                                                         elementShape, elementDataType,
                                                                         filledAccumulator,
                                                                         i));
        }
        this->initialize();
    }

    ~Us4ROutputBuffer() override {
        ::operator delete(dataBuffer, std::align_val_t(DATA_ALIGNMENT));
        getDefaultLogger()->log(LogSeverity::DEBUG, "Released the output buffer.");
    }

    void registerOnNewDataCallback(framework::OnNewDataCallback &callback) override {
        this->onNewDataCallback = callback;
    }

    [[nodiscard]] const framework::OnNewDataCallback &getOnNewDataCallback() const {
        return this->onNewDataCallback;
    }

    void registerOnOverflowCallback(framework::OnOverflowCallback &callback) override {
        this->onOverflowCallback = callback;
    }

    void registerShutdownCallback(framework::OnShutdownCallback &callback) override {
        this->onShutdownCallback = callback;
    }

    [[nodiscard]] size_t getNumberOfElements() const override {
        return elements.size();
    }

    BufferElement::SharedHandle getElement(size_t i) override {
        return std::static_pointer_cast<BufferElement>(elements[i]);
    }

    uint8 *getAddress(uint16 elementNumber, Ordinal us4oem) {
        return reinterpret_cast<uint8 *>(this->elements[elementNumber]->getAddress()) + us4oemOffsets[us4oem];
    }

    /**
     * Returns a total size of the buffer, the number of bytes.
     */
    [[nodiscard]] size_t getElementSize() const override {
        return elementSize;
    }

    /**
     * Signals the readiness of new data acquired by the n-th Us4OEM module.
     *
     * This function should be called by us4oem interrupt callbacks.
     *
     * @param n us4oem ordinal number
     *
     *  @return true if the buffer signal was successful, false otherwise (e.g. the queue was shut down).
     */
    bool signal(Ordinal n, uint16 elementNr) {
        std::unique_lock<std::mutex> guard(mutex);
        if(this->state != State::RUNNING) {
            getDefaultLogger()->log(LogSeverity::DEBUG, "Signal queue shutdown.");
            return false;
        }
        this->validateState();
        auto &element = this->elements[elementNr];
        try {
            element->signal(n);
        } catch(const IllegalArgumentException &e) {
            this->markAsInvalid();
            throw e;
        }
        if(element->isElementReady()) {
            guard.unlock();
            onNewDataCallback(elements[elementNr]);
        } else {
            guard.unlock();
        }
        return true;
    }

    void markAsInvalid() {
        std::unique_lock<std::mutex> guard(mutex);
        if(this->state != State::INVALID) {
            this->state = State::INVALID;
            for(auto &element: elements) {
                element->markAsInvalid();
            }
            this->onOverflowCallback();
        }
    }

    void shutdown() {
        std::unique_lock<std::mutex> guard(mutex);
        this->onShutdownCallback();
        this->state = State::SHUTDOWN;
        guard.unlock();
    }

    void resetState() {
        this->state = State::INVALID;
        this->initialize();
        this->state = State::RUNNING;
    }

    void initialize() {
        for(auto &element: elements) {
            element->resetState();
        }
    }

    void registerReleaseFunction(size_t element, std::function<void()> &releaseFunction) {
        this->elements[element]->registerReleaseFunction(releaseFunction);
    }

private:
    std::mutex mutex;
    /** A size of a single element IN number of BYTES. */
    size_t elementSize;
    /**  Total size in the number of elements. */
    int16 *dataBuffer;
    /** Host buffer elements */
    std::vector<Us4ROutputBufferElement::SharedHandle> elements;
    /** Relative addresses where us4oem modules will write. IN NUMBER OF BYTES. */
    std::vector<size_t> us4oemOffsets;
    // Callback that should be called once new data arrive.
    framework::OnNewDataCallback onNewDataCallback;
    framework::OnOverflowCallback onOverflowCallback{[]() {}};
    framework::OnShutdownCallback onShutdownCallback{[]() {}};

    // State management
    enum class State {
        RUNNING, SHUTDOWN, INVALID
    };
    State state{State::RUNNING};

    /**
     * Throws IllegalStateException when the buffer is in invalid state.
     *
     * @return true if the queue execution should continue, false otherwise.
     */
    void validateState() {
        if(this->state == State::INVALID) {
            throw ::arrus::IllegalStateException(
                "The buffer is in invalid state "
                "(probably some data transfer overflow happened).");
        } else if(this->state == State::SHUTDOWN) {
            throw ::arrus::IllegalStateException(
                "The data buffer has been turned off.");
        }
    }
};

}

#endif //ARRUS_CORE_DEVICES_US4R_US4ROUTPUTBUFFER_H


text:
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_poly.h>

void mgsl_poly_complex_eval(const double c[], const int len, const gsl_complex *z, gsl_complex *ret)
{
  *ret = gsl_poly_complex_eval(c, len, *z);
}

void mgsl_complex_poly_complex_eval(gsl_complex **c, const int len, const gsl_complex *z, gsl_complex *ret)
{
  gsl_complex *zc = calloc(len, sizeof(gsl_complex));
  int i;
  for(i=0; i<len; i++) zc[i] = *c[i];
  *ret = gsl_complex_poly_complex_eval(zc, len, *z);
}

max_stars_repo_name: frithnanth/raku-Math-Libgsl-Polynomial

max_issues_repo_name: frithnanth/raku-Math-Libgsl-Polynomial

max_forks_repo_name: frithnanth/raku-Math-Libgsl-Polynomial