text:
#include <petsc.h>
#include <petsc/private/fortranimpl.h>

#define MEMORY_ALIGNMENT
//#define OFFLOAD_MODE
//#define HAVE_REDUNDANT_WORK
//#define HAVE_EDGE_COLORING

#ifdef _OPENMP
#include <omp.h>
#ifdef OFFLOAD_MODE
#include <offload.h>
#endif
#ifndef HAVE_EDGE_COLORING
#include <metis.h>
#endif
#endif

#define max_colors    200
#define MAX_FILE_SIZE 256
#define mem_alignment 64

#define REAL          double

#if defined(MEMORY_ALIGNMENT)
#define ICALLOC(size, y) ierr = \
        posix_memalign((void **) y, mem_alignment, PetscMax(size, 1) * sizeof(int)); \
        CHKERRQ(ierr);
#define FCALLOC(size, y) ierr = \
        posix_memalign((void **) y, mem_alignment, PetscMax(size, 1) * sizeof(double)); \
        CHKERRQ(ierr);
#else
#define ICALLOC(size, y) ierr = \
        PetscMalloc((PetscMax(size, 1)) * sizeof(int), y);\
        CHKERRQ(ierr);
#define FCALLOC(size, y) ierr = \
        PetscMalloc((PetscMax(size, 1)) * sizeof(double), y);\
        CHKERRQ(ierr);
#endif

#if defined(OFFLOAD_MODE)
#pragma offload_attribute(push, target(mic))
int max_threads;
#pragma offload_attribute(pop)
#else
int max_threads;
#endif

#ifdef OFFLOAD_MODE
#define  ALLOC  alloc_if(1)
#define  REUSE  alloc_if(0)
#define  FREE   free_if(1)
#define  KEEP   free_if(0)
#endif

int     rank;                 // MPI process rank
int     size;                 // Number of MPI processes
int     rstart;               // Used for identifying the ghost points in the mesh

typedef struct gxy{                                   /* GRID STRUCTURE             */
  int         nnodes;                                 /* Number of nodes            */
  int         ncell;                                  /* Number of cells            */
  int         nedge;                                  /* Number of edges            */
  int         nnbound;
  int         nvbound;
  int         nfbound;
  int         nnfacet;
  int         nvfacet;
  int         nffacet;
  int         ntte;
  int         nsface;                                 /* Total # of solid faces     */
  int         nvface;                                 /* Total # of viscous faces   */
  int         nfface;                                 /* Total # of far field faces */
  int         nsnode;                                 /* Total # of solid nodes     */
  int         nvnode;                                 /* Total # of viscous nodes   */
  int         nfnode;                                 /* Total # of far field nodes */
  double      *rxy;
# if defined(OFFLOAD_MODE)                            /* Offloading mode            */
# pragma offload_attribute(push, target(mic))         /* Define in the coprocessor  */
  int         nnodesLoc;                              /* Local nodes                */
  int         nedgeLoc;                               /* Local edges                */
  int         nvertices;                              /* Local nodes + ghost nodes  */
  double      *xyz;                                   /* Node Coordinates           */
  int         *eptr;                                  /* Edge pointers              */
# ifndef HAVE_EDGE_COLORING                           /* Using MeTiS                */
  int         nedgeAllThr;
  int         *part_thr;
  int         *nedge_thr;
  int         *edge_thr;
  double      *xyzn_thr;
# else                                                /* Using edge coloring        */
  int         ncolor;                                 /* Number of colors           */
  int         ncount[max_colors];                     /* No. of faces in color      */
# endif
# pragma offload_attribute(pop)
# else                                                /* Not Offloading mode        */
  int         nnodesLoc;                              /* Local nodes                */
  int         nedgeLoc;                               /* Local edges                */
  int         nvertices;                              /* Local nodes + ghost nodes  */
  double      *xyz;                                   /* Node Coordinates           */
  int         *eptr;                                  /* Edge pointers              */
# ifndef HAVE_EDGE_COLORING
  int         nedgeAllThr;
  int         *part_thr;
  int         *nedge_thr;
  int         *edge_thr;
  double      *xyzn_thr;
# else
  int         ncolor;
  int         ncount[max_colors];
# endif
# endif

# ifdef HAVE_REDUNDANT_WORK
  double      *resd;
# endif
  int         jvisc;                                  /* 0 = Euler                  */
                                                      /* 1 = laminar no visc LHS    */
                                                      /* 2 = laminar w/ visc LHS    */
                                                      /* 3 = turb B-B no visc LHS   */
                                                      /* 4 = turb B-B w/ visc LHS   */
                                                      /* 5 = turb Splrt no visc LHS */
                                                      /* 6 = turb Splrt w/ visc LHS */
  int         ileast;                                 /* 1 = Lst square gradient    */
  int         *isface;                                /* Face # of solid faces      */
  int         *ifface;                                /* Face # of far field faces  */
  int         *ivface;                                /* Face # of viscous faces    */
  int         *isford;                                /* Copies of isface, ifface,  */
  int         *ifford;                                /*  and ivface used for       */
  int         *ivford;                                /*  ordering                  */
  int         *isnode;                                /* Node # of solid nodes      */
  int         *ivnode;                                /* Node # of viscous nodes    */
  int         *ifnode;                                /* Node # of far field nodes  */
  int         *c2n;                                   /* Cell-to-node pointers      */
  int         *c2e;                                   /* Cell-to-edge pointers      */
  int         *ia, *ja;                               /* Stuff for ILU(0)           */
  int         *loc2pet;                               /* local to PETSc mapping     */
  int         *loc2glo;                               /* local to global mapping    */
  double      *area;                                  /* Area of control volumes    */
  int         *nntet,*nnpts;
  int         *f2ntn;
  double      *sxn, *syn, *szn, *sa;                  /* Normals at solid nodes     */
  int         *v2p;				                            /* Vertex to processor mapping*/
  int         *sface_bit, *vface_bit;
  int         *nvtet,*nvpts,*nftet,*nfpts;
  int         *f2ntv,*f2ntf;
  double      *vxn, *vyn, *vzn, *va;                  /* Normals at viscous nodes   */
  double      *fxn, *fyn, *fzn, *fa;                  /* Normals at far field nodes */
  double      *xyzn;                                  /* Normal to faces and length */
  double      *cdt;                                   /* Local time step            */
  double      *phi;                                   /* Flux limiter               */
  double      *dft1, *dft2;                           /* Turb mod linearization     */
  double      *slen;                                  /* Generalized distance       */
  double      *turbre;                                /* nu x turb Reynolds #       */
  double      *amut;                                  /* Turbulent mu (viscosity)   */
  double      *turbres;                               /* Turbulent residual         */
  double      *us, *vs, *as;                          /* For linearizing viscous    */
  int         nsnodeLoc, nvnodeLoc, nfnodeLoc;
  int         nnfacetLoc, nvfacetLoc, nffacetLoc;

  /* PETSc datastructures and other related info */
  Vec         qnode;                                  /* Global distributed solution vector        */
  Vec         qnodeLoc;                               /* Local sequential solution vector          */
  Vec         qold;                                   /* Global distributed solution vector        */
  Vec         res;                                    /* Residual                                  */
  Vec         grad;                                   /* Gradient Vector	                         */
  Vec         gradLoc;                  	            /* Local Gradient Vector	                   */
  Vec         B;                                      /* Right hand side                           */
  Mat         A;                                      /* Left hand side                            */
  VecScatter  scatter, gradScatter;                   /* Scatter between local and global vectors  */

} GRID;

typedef struct{                                 /* GENERAL INFORMATION        */
  double      title[20];                        /* Title line                 */
# ifdef OFFLOAD_MODE
# pragma offload_attribute(push, target(mic))
  double      beta;                             /* Artificial Compress. Param */
# pragma offload_attribute(pop)
# else
  double      beta;                             /* Artificial Compress. Param */
# endif
  double      alpha;                            /* Angle of attack            */
  double      Re;                               /* Reynolds number            */
  double      dt;                               /* Input cfl                  */
  double      tot;                              /* total computer time        */
  double      res0;                             /* Begining residual          */
  double      resc;                             /* Current residual           */
  int         ntt;                              /* A counter                  */
  int         mseq;                             /* Mesh sequencing            */
  int         ivisc;                            /* 0 = Euler                  */
                                                /* 1 = laminar no visc LHS    */
                                                /* 2 = laminar w/ visc LHS    */
                                                /* 3 = turb BB no visc LHS    */
                                                /* 4 = turb BB w/ visc LHS    */
                                                /* 5 = turb SA w/ visc LHS    */
                                                /* 6 = turb SA w/ visc LHS    */
  int         irest;                            /* for restarts irest = 1     */
  int         icyc;                             /* iterations completed       */
  int         ihane;                            /* ihane = 0 for van leer fds */
                                                /*       = 1 for hanel flux   */
                                                /*       = 2 for Roe's fds    */
  int         ntturb;                           /* Counter for turbulence     */
} CINFO;                                        /* COMMON INFO                */

typedef struct {                                /* FLOW SOLVER VARIABLES      */
  double      cfl1;                             /* Starting CFL number        */
  double      cfl2;                             /* Ending   CFL number        */
  int         nsmoth;                           /* How many its for Res smooth*/
  int         iflim;                            /* 1=use limiter 0=no limiter */
  int         itran;                            /* 1=transition (spalart only)*/
  int         nbtran;                           /* No. of transition points   */
  int         jupdate;                          /* For freezing Jacobians */
  int         nstage;                           /* Number of subiterations    */
  int         ncyct;                            /* Subiterations for turb mod */
  int         iramp;                            /* Ramp CFL over iramp iters  */
  int         nitfo;                            /* Iterations first order     */
} CRUNGE;                                       /* COMMON RUNGE               */

typedef struct{                                 /*============================*/
  double      gtol;                             /* linear system tolerence    */
  int         icycle;                           /* Number of GMRES iterations */
  int         nsrch;                            /* Dimension of Krylov        */
  int         ilu0;                             /* 1 for ILU(0)               */
  int         ifcn;                             /* 0=fcn2 1=fcneval(nwt Krlv) */
} CGMCOM;                                       /* COMMON GMCOM               */

// Time stepping
typedef struct {
  Vec         qnew, qold, func;
  double      fnorm_ini, dt_ini, cfl_ini;
  double      ptime;
  double      cfl_max, max_time;
  double      fnorm, dt, cfl;
  double      fnorm_ratio;
  int         ires, iramp, itstep;
  int         max_steps, print_freq;
  int         LocalTimeStepping;
} TstepCtx;

// Encapsulate the grid data structure information and time step
// data structure information
typedef struct {
  GRID        *grid;      // Pointer to grid info
  TstepCtx    *tsCtx;     // Pointer to time stepping context
  PetscBool   PreLoading;
} AppCtx;

// C routines
// =============================
int GetLocalOrdering(GRID *);
#ifdef _OPENMP
#ifdef HAVE_EDGE_COLORING
int EdgeColoring(int nnodes, int nedge, int *e2n, int *eperm, int *ncolor, int *ncount);
#endif
#endif
int SetPetscDS(GRID *, TstepCtx *);
int FormJacobian(SNES, Vec, Mat, Mat, void *);
int FormFunction(SNES, Vec, Vec, void *);
int FormInitialGuess(SNES, GRID *);
int Update(SNES, void *);
int ComputeTimeStep(SNES, int, void *);

/* =================  Fortran routines called from C ======================= */

#if defined(PETSC_HAVE_FORTRAN_UNDERSCORE)
#define f77name(ucase,lcase,lcbar) lcbar
#elif defined(PETSC_HAVE_FORTRAN_CAPS)
#define f77name(ucase,lcase,lcbar) ucase
#else
#define f77name(ucase,lcase,lcbar) lcase
#endif
#define f77FORLINK  f77name(FORLINK,forlink,forlink_)
#define f77CLINK    f77name(CLINK,clink,clink_)
#define f77GETIA    f77name(GETIA,getia,getia_)
#define f77GETJA    f77name(GETJA,getja,getja_)
#define f77SORTER   f77name(SORTER,sorter,sorter_)
#define f77SUMGS    f77name(SUMGS,sumgs,sumgs_)
#define f77INIT     f77name(INIT,init,init_)
#define f77LSTGS    f77name(LSTGS,lstgs,lstgs_)
#define f77GETRES   f77name(GETRES,getres,getres_)
#define f77FILLA    f77name(FILLA,filla,filla_)
#define f77FORCE    f77name(FORCE,force,force_)

EXTERN_C_BEGIN

extern void PETSC_STDCALL f77FORLINK(void);
extern void PETSC_STDCALL f77GETIA(int*,int*,int*,int*,int*,int*);
extern void PETSC_STDCALL f77GETJA(int*,int*,int*,int*,int*,int*,int*);
extern void PETSC_STDCALL f77SUMGS(int*,int*,int*,PetscScalar*,PetscScalar*,int*,int*);
extern void PETSC_STDCALL f77INIT(int*,PetscScalar*,PetscScalar*,PetscScalar*,int*,int*,
                                  int*);
extern void PETSC_STDCALL f77LSTGS(int*,int*,int*,PetscScalar*,PetscScalar*,PetscScalar*,
                                   PetscScalar*,int*,int*);
extern void PETSC_STDCALL f77GETRES(int*,int*,int*,int*,int*,int*,int*,int*,int*,int*,int*,
                                   int*,int*,int*,int*,int*,int*,int*,int*,int*,int*,int*,
                                   int*,int*,int*,int*,int*,PetscScalar*,PetscScalar*,
                                   PetscScalar*,PetscScalar*,PetscScalar*,PetscScalar*,
                                   PetscScalar*,PetscScalar*,PetscScalar*,PetscScalar*,
                                   PetscScalar*,PetscScalar*,PetscScalar*,PetscScalar*,
                                   PetscScalar*,PetscScalar*,PetscScalar*,PetscScalar*,int*,
                                   int*,PetscScalar*,PetscScalar*,PetscScalar*,PetscScalar*,
                                   PetscScalar*,int*,
#                                  if defined(_OPENMP)
                                   int*,
#                                  if defined(HAVE_EDGE_COLORING)
                                   int*, int*,
#                                  elif defined(HAVE_REDUNDANT_WORK)
                                   PetscScalar*,
#                                  else
                                   int*,
                                   int*,int*,int*,PetscScalar*,
#                                  endif
#                                  ifdef OFFLOAD_MODE
                                   PetscScalar*, int*,
#                                  endif
#                                  endif
                                   int*,int*,int*
                                   );
extern void PETSC_STDCALL f77FILLA(int*,int*,int*,int*,int*,PetscScalar*,PetscScalar*,PetscScalar*,
                                   PetscScalar*,PetscScalar*,PetscScalar*,int*,int*,int*,int*,int*,
                                   int*,PetscScalar*,Mat*,PetscScalar*,PetscScalar*,PetscScalar*,
                                   PetscScalar*,int*,int*);
extern void PETSC_STDCALL f77FORCE(int*,int*,int*,int*,int*,int*,int*,int*,int*,int*,int*,
                                   PetscScalar*,PetscScalar*,int*,int*,PetscScalar*,PetscScalar*,
                                   PetscScalar*,int*,int*);
EXTERN_C_END

ext: h


text:
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <string.h>
#include "complex.h"
#include "api.h"

static void ke_complex_alloc(sml_t* sml) {
	int y = sml_pop_int(sml);
	int x = sml_pop_int(sml);
	sml_push_new_complex(sml);
	GSL_SET_COMPLEX(sml_get_complex_adr, (double)x, (double)y);
}

static void ke_complex_set_real(sml_t* sml) {
	double x = sml_pop_real(sml);
	gsl_complex *z = sml_pop_complex_adr(sml);
	GSL_SET_REAL(z, x);
}

static void ke_complex_set_imag(sml_t* sml) {
	double y = sml_pop_real(sml);
	gsl_complex *z = sml_pop_complex_adr(sml);
	GSL_SET_IMAG(z, y);
}

static void ke_complex_arg(sml_t* sml) {
	gsl_complex z = sml_pop_complex(sml);
	double r = gsl_complex_arg(z);
	sml_push_real(sml,r);
}

static void ke_complex_abs(sml_t* sml) {
	gsl_complex z = sml_pop_complex(sml);
	double r = gsl_complex_abs(z);
	sml_push_real(sml, r);
}

static void ke_complex_abs2(sml_t* sml) {
	gsl_complex z = sml_pop_complex(sml);
	double r = gsl_complex_abs2(z);
	sml_push_real(sml, r);
}

static void ke_complex_logabs(sml_t* sml) {
	gsl_complex z = sml_pop_complex(sml);
	double r = gsl_complex_logabs(z);
	sml_push_real(sml, r);
}

static void ke_complex_add(sml_t* sml) {
	gsl_complex b = sml_pop_complex(sml);
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_add(a, b);
	sml_push_complex(sml, z);
}

static void ke_complex_sub(sml_t* sml) {
	gsl_complex b = sml_pop_complex(sml);
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_sub(a, b);
	sml_push_complex(sml, z);
}

static void ke_complex_mul(sml_t* sml) {
	gsl_complex b = sml_pop_complex(sml);
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_mul(a, b);
	sml_push_complex(sml, z);
}

static void ke_complex_div(sml_t* sml) {
	gsl_complex b = sml_pop_complex(sml);
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_div(a, b);
	sml_push_complex(sml, z);
}

static void ke_complex_add_real(sml_t* sml) {
	double x = sml_pop_real(sml);
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_add_real(a, x);
	sml_push_complex(sml, z);
}

static void ke_complex_sub_real(sml_t* sml) {
	double x = sml_pop_real(sml);
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_sub_real(a, x);
	sml_push_complex(sml, z);
}

static void ke_complex_mul_real(sml_t* sml) {
	double x = sml_pop_real(sml);
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_mul_real(a, x);
	sml_push_complex(sml, z);
}

static void ke_complex_div_real(sml_t* sml) {
	double x = sml_pop_real(sml);
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_div_real(a, x);
	sml_push_complex(sml, z);
}

static void ke_complex_add_imag(sml_t* sml) {
	double y = sml_pop_real(sml);
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_add_imag(a, y);
	sml_push_complex(sml, z);
}

static void ke_complex_sub_imag(sml_t* sml) {
	double y = sml_pop_real(sml);
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_sub_imag(a, y);
	sml_push_complex(sml, z);
}

static void ke_complex_mul_imag(sml_t* sml) {
	double y = sml_pop_real(sml);
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_mul_imag(a, y);
	sml_push_complex(sml, z);
}

static void ke_complex_div_imag(sml_t* sml) {
	double y = sml_pop_real(sml);
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_div_imag(a, y);
	sml_push_complex(sml, z);
}

static void ke_complex_conjugate(sml_t* sml) {
	gsl_complex z = sml_pop_complex(sml);
	z = gsl_complex_conjugate(z);
	sml_push_complex(sml, z);
}

static void ke_complex_inverse(sml_t* sml) {
	gsl_complex z = sml_pop_complex(sml);
	z = gsl_complex_inverse(z);
	sml_push_complex(sml, z);
}

static void ke_complex_negative(sml_t* sml) {
	gsl_complex z = sml_pop_complex(sml);
	z = gsl_complex_negative(z);
	sml_push_complex(sml, z);
}

static void ke_complex_sqrt(sml_t* sml) {
	gsl_complex z = sml_pop_complex(sml);
	z = gsl_complex_sqrt(z);
	sml_push_complex(sml, z);
}

static void ke_complex_sqrt_real(sml_t* sml) {
	double x = sml_pop_real(sml);
	gsl_complex z = gsl_complex_sqrt_real(x);
	sml_push_complex(sml, z);
}

static void ke_complex_pow(sml_t* sml) {
	gsl_complex b = sml_pop_complex(sml);
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_pow(a, b);
	sml_push_complex(sml, z);
}

static void ke_complex_pow_real(sml_t* sml) {
	double b = sml_pop_real(sml);
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_pow_real(a, b);
	sml_push_complex(sml, z);
}

static void ke_complex_exp(sml_t* sml) {
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_exp(a);
	sml_push_complex(sml, z);
}

static void ke_complex_log(sml_t* sml) {
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_log(a);
	sml_push_complex(sml, z);
}

static void ke_complex_log10(sml_t* sml) {
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_log10(a);
	sml_push_complex(sml, z);
}

static void ke_complex_log_b(sml_t* sml) {
	gsl_complex b = sml_pop_complex(sml);
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_log_b(a, b);
	sml_push_complex(sml, z);
}

static void ke_complex_sin(sml_t* sml) {
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_sin(a);
	sml_push_complex(sml, z);
}

static void ke_complex_cos(sml_t* sml) {
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_cos(a);
	sml_push_complex(sml, z);
}

static void ke_complex_tan(sml_t* sml) {
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_tan(a);
	sml_push_complex(sml, z);
}

static void ke_complex_sec(sml_t* sml) {
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_sec(a);
	sml_push_complex(sml, z);
}

static void ke_complex_csc(sml_t* sml) {
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_csc(a);
	sml_push_complex(sml, z);
}

static void ke_complex_cot(sml_t* sml) {
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_cot(a);
	sml_push_complex(sml, z);
}

static void ke_complex_arcsin(sml_t* sml) {
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_arcsin(a);
	sml_push_complex(sml, z);
}

static void ke_complex_arcsin_real(sml_t* sml) {
	double x = sml_pop_real(sml);
	gsl_complex z = gsl_complex_arcsin_real(x);
	sml_push_complex(sml, z);
}

static void ke_complex_arccos(sml_t* sml) {
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_arccos(a);
	sml_push_complex(sml, z);
}

static void ke_complex_arccos_real(sml_t* sml) {
	double x = sml_pop_real(sml);
	gsl_complex z = gsl_complex_arccos_real(x);
	sml_push_complex(sml, z);
}

static void ke_complex_arctan(sml_t* sml) {
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_arctan(a);
	sml_push_complex(sml, z);
}

static void ke_complex_arcsec(sml_t* sml) {
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_arcsec(a);
	sml_push_complex(sml, z);
}

static void ke_complex_arcsec_real(sml_t* sml) {
	double x = sml_pop_real(sml);
	gsl_complex z = gsl_complex_arcsec_real(x);
	sml_push_complex(sml, z);
}

static void ke_complex_arccsc(sml_t* sml) {
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_arccsc(a);
	sml_push_complex(sml, z);
}

static void ke_complex_arccsc_real(sml_t* sml) {
	double x = sml_pop_real(sml);
	gsl_complex z = gsl_complex_arccsc_real(x);
	sml_push_complex(sml, z);
}

static void ke_complex_arccot(sml_t* sml) {
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_arccot(a);
	sml_push_complex(sml, z);
}

static void ke_complex_sinh(sml_t* sml) {
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_sinh(a);
	sml_push_complex(sml, z);
}

static void ke_complex_cosh(sml_t* sml) {
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_cosh(a);
	sml_push_complex(sml, z);
}

static void ke_complex_tanh(sml_t* sml) {
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_tanh(a);
	sml_push_complex(sml, z);
}

static void ke_complex_sech(sml_t* sml) {
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_sech(a);
	sml_push_complex(sml, z);
}

static void ke_complex_csch(sml_t* sml) {
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_csch(a);
	sml_push_complex(sml, z);
}

static void ke_complex_coth(sml_t* sml) {
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_coth(a);
	sml_push_complex(sml, z);
}

static void ke_complex_arcsinh(sml_t* sml) {
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_arcsinh(a);
	sml_push_complex(sml, z);
}

static void ke_complex_arccosh(sml_t* sml) {
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_arccosh(a);
	sml_push_complex(sml, z);
}

static void ke_complex_arctanh(sml_t* sml) {
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_arctanh(a);
	sml_push_complex(sml, z);
}

static void ke_complex_arcsech(sml_t* sml) {
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_arcsech(a);
	sml_push_complex(sml, z);
}

static void ke_complex_arccsch(sml_t* sml) {
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_arccsch(a);
	sml_push_complex(sml, z);
}

static void ke_complex_arccoth(sml_t* sml) {
	gsl_complex a = sml_pop_complex(sml);
	gsl_complex z = gsl_complex_arccoth(a);
	sml_push_complex(sml, z);
}

void ke_complex_hash(sml_t* sml) {
    ke_hash_add(sml, (fncp)&ke_complex_alloc, COMPLEX);
    ke_hash_add(sml, (fncp)&ke_complex_set_real, COMPLEX_SET_REAL);
    ke_hash_add(sml, (fncp)&ke_complex_set_imag, COMPLEX_SET_IMAG);
    ke_hash_add(sml, (fncp)&ke_complex_arg, COMPLEX_ARG);
    ke_hash_add(sml, (fncp)&ke_complex_abs, COMPLEX_ABS);
    ke_hash_add(sml, (fncp)&ke_complex_abs2, COMPLEX_ABS2);
    ke_hash_add(sml, (fncp)&ke_complex_logabs, COMPLEX_LOGABS);
    ke_hash_add(sml, (fncp)&ke_complex_add, COMPLEX_ADD);
    ke_hash_add(sml, (fncp)&ke_complex_sub, COMPLEX_SUB);
    ke_hash_add(sml, (fncp)&ke_complex_mul, COMPLEX_MUL);
    ke_hash_add(sml, (fncp)&ke_complex_div, COMPLEX_DIV);
    ke_hash_add(sml, (fncp)&ke_complex_add_real, COMPLEX_ADD_REAL);
    ke_hash_add(sml, (fncp)&ke_complex_sub_real, COMPLEX_SUB_REAL);
    ke_hash_add(sml, (fncp)&ke_complex_mul_real, COMPLEX_MUL_REAL);
    ke_hash_add(sml, (fncp)&ke_complex_div_real, COMPLEX_DIV_REAL);
    ke_hash_add(sml, (fncp)&ke_complex_add_imag, COMPLEX_ADD_IMAG);
    ke_hash_add(sml, (fncp)&ke_complex_sub_imag, COMPLEX_SUB_IMAG);
    ke_hash_add(sml, (fncp)&ke_complex_mul_imag, COMPLEX_MUL_IMAG);
    ke_hash_add(sml, (fncp)&ke_complex_div_imag, COMPLEX_DIV_IMAG);
    ke_hash_add(sml, (fncp)&ke_complex_conjugate, COMPLEX_CONJUGATE);
    ke_hash_add(sml, (fncp)&ke_complex_inverse, COMPLEX_INVERSE);
    ke_hash_add(sml, (fncp)&ke_complex_negative, COMPLEX_NEGATIVE);
    ke_hash_add(sml, (fncp)&ke_complex_sqrt, COMPLEX_SQRT);
    ke_hash_add(sml, (fncp)&ke_complex_sqrt_real, COMPLEX_SQRT_REAL);
    ke_hash_add(sml, (fncp)&ke_complex_pow, COMPLEX_POW);
    ke_hash_add(sml, (fncp)&ke_complex_pow_real, COMPLEX_POW_REAL);
    ke_hash_add(sml, (fncp)&ke_complex_exp, COMPLEX_EXP);
    ke_hash_add(sml, (fncp)&ke_complex_log, COMPLEX_LOG);
    ke_hash_add(sml, (fncp)&ke_complex_log10, COMPLEX_LOG10);
    ke_hash_add(sml, (fncp)&ke_complex_log_b, COMPLEX_LOG_B);
    ke_hash_add(sml, (fncp)&ke_complex_sin, COMPLEX_SIN);
    ke_hash_add(sml, (fncp)&ke_complex_cos, COMPLEX_COS);
    ke_hash_add(sml, (fncp)&ke_complex_tan, COMPLEX_TAN);
    ke_hash_add(sml, (fncp)&ke_complex_sec, COMPLEX_SEC);
    ke_hash_add(sml, (fncp)&ke_complex_csc, COMPLEX_CSC);
    ke_hash_add(sml, (fncp)&ke_complex_cot, COMPLEX_COT);
    ke_hash_add(sml, (fncp)&ke_complex_arcsin, COMPLEX_ARCSIN);
    ke_hash_add(sml, (fncp)&ke_complex_arcsin_real, COMPLEX_ARCSIN_REAL);
    ke_hash_add(sml, (fncp)&ke_complex_arccos, COMPLEX_ARCCOS);
    ke_hash_add(sml, (fncp)&ke_complex_arccos_real, COMPLEX_ARCCOS_REAL);
    ke_hash_add(sml, (fncp)&ke_complex_arctan, COMPLEX_ARCTAN);
    ke_hash_add(sml, (fncp)&ke_complex_arcsec, COMPLEX_ARCSEC);
    ke_hash_add(sml, (fncp)&ke_complex_arcsec_real, COMPLEX_ARCSEC_REAL);
    ke_hash_add(sml, (fncp)&ke_complex_arccsc, COMPLEX_ARCCSC);
    ke_hash_add(sml, (fncp)&ke_complex_arccsc_real, COMPLEX_ARCCSC_REAL);
    ke_hash_add(sml, (fncp)&ke_complex_arccot, COMPLEX_ARCCOT);
    ke_hash_add(sml, (fncp)&ke_complex_sinh, COMPLEX_SINH);
    ke_hash_add(sml, (fncp)&ke_complex_cosh, COMPLEX_COSH);
    ke_hash_add(sml, (fncp)&ke_complex_tanh, COMPLEX_TANH);
    ke_hash_add(sml, (fncp)&ke_complex_sech, COMPLEX_SECH);
    ke_hash_add(sml, (fncp)&ke_complex_csch, COMPLEX_CSCH);
    ke_hash_add(sml, (fncp)&ke_complex_coth, COMPLEX_COTH);
    ke_hash_add(sml, (fncp)&ke_complex_arcsinh, COMPLEX_ARCSINH);
    ke_hash_add(sml, (fncp)&ke_complex_arccosh, COMPLEX_ARCCOSH);
    ke_hash_add(sml, (fncp)&ke_complex_arctanh, COMPLEX_ARCTANH);
    ke_hash_add(sml, (fncp)&ke_complex_arcsech, COMPLEX_ARCSECH);
    ke_hash_add(sml, (fncp)&ke_complex_arccsch, COMPLEX_ARCCSCH);
    ke_hash_add(sml, (fncp)&ke_complex_arccoth, COMPLEX_ARCCOTH);
}

void ke_complex_print(sml_t* sml,token_t *k) {
    printf("Complex: %s\n", k->name);
    printf("%g,%g \n", GSL_REAL(k->obj.tcomplex), GSL_IMAG(k->obj.tcomplex));
}

ext: c

max_stars_repo_name: vinej/sml

max_issues_repo_name: vinej/sml

max_forks_repo_name: vinej/sml


text:
#pragma once
#include "cru/common/Base.h"
#include "cru/common/Event.h"
#include "cru/platform/graphics/Base.h"
#include "cru/platform/graphics/Brush.h"
#include "cru/platform/graphics/Geometry.h"
#include "cru/platform/graphics/Painter.h"
#include "cru/platform/gui/Cursor.h"
#include "cru/platform/gui/TimerHelper.h"
#include "cru/platform/gui/UiApplication.h"
#include "cru/ui/Base.h"
#include "cru/ui/controls/Control.h"
#include "cru/ui/helper/ClickDetector.h"

#include <gsl/pointers>
#include <memory>
#include <optional>
#include <unordered_map>

namespace cru::ui::render {
class ScrollRenderObject;

enum class ScrollKind { Absolute, Relative, Page, Line };

struct Scroll {
  Direction direction;
  ScrollKind kind;
  // For absolute, the new scroll position. Otherwise, offset.
  float value;
};

enum class ScrollBarAreaKind {
  UpArrow,    // Line up
  DownArrow,  // Line down
  UpSlot,     // Page up
  DownSlot,   // Page down
  Thumb
};

enum class ScrollBarBrushUsageKind { Arrow, ArrowBackground, Slot, Thumb };
enum class ScrollBarBrushStateKind { Normal, Hover, Press, Disable };

String CRU_UI_API GenerateScrollBarThemeColorKey(ScrollBarBrushUsageKind usage,
                                                 ScrollBarBrushStateKind state);

class CRU_UI_API ScrollBar : public Object {
 public:
  ScrollBar(gsl::not_null<ScrollRenderObject*> render_object,
            Direction direction);

  CRU_DELETE_COPY(ScrollBar)
  CRU_DELETE_MOVE(ScrollBar)

  ~ScrollBar() override;

 public:
  Direction GetDirection() const { return direction_; }

  bool IsEnabled() const { return is_enabled_; }
  void SetEnabled(bool value);

  bool IsExpanded() const { return is_expanded_; }
  void SetExpanded(bool value);

  void Draw(platform::graphics::IPainter* painter);

  IEvent<Scroll>* ScrollAttemptEvent() { return &scroll_attempt_event_; }

  void InstallHandlers(controls::Control* control);
  void UninstallHandlers() { InstallHandlers(nullptr); }

  gsl::not_null<std::shared_ptr<platform::graphics::IBrush>>
  GetCollapsedThumbBrush();
  // Brush could be nullptr to use the theme brush.
  void SetCollapsedThumbBrush(
      std::shared_ptr<platform::graphics::IBrush> brush);
  gsl::not_null<std::shared_ptr<platform::graphics::IBrush>> GetBrush(
      ScrollBarBrushUsageKind usage, ScrollBarBrushStateKind state);
  // Brush could be nullptr to use the theme brush.
  void SetBrush(ScrollBarBrushUsageKind usage, ScrollBarBrushStateKind state,
                std::shared_ptr<platform::graphics::IBrush> brush);

 protected:
  void OnDraw(platform::graphics::IPainter* painter, bool expand);

  virtual void DrawUpArrow(
      platform::graphics::IPainter* painter, const Rect& area,
      gsl::not_null<platform::graphics::IBrush*> arrow_brush,
      gsl::not_null<platform::graphics::IBrush*> background_brush) = 0;
  virtual void DrawDownArrow(
      platform::graphics::IPainter* painter, const Rect& area,
      gsl::not_null<platform::graphics::IBrush*> arrow_brush,
      gsl::not_null<platform::graphics::IBrush*> background_brush) = 0;

  std::optional<ScrollBarAreaKind> ExpandedHitTest(const Point& point);

  virtual bool IsShowBar() = 0;

  virtual std::optional<Rect> GetExpandedAreaRect(
      ScrollBarAreaKind area_kind) = 0;
  virtual std::optional<Rect> GetCollapsedTriggerExpandAreaRect() = 0;
  virtual std::optional<Rect> GetCollapsedThumbRect() = 0;

  virtual float CalculateNewScrollPosition(const Rect& thumb_original_rect,
                                           const Point& mouse_offset) = 0;

  virtual bool CanScrollUp() = 0;
  virtual bool CanScrollDown() = 0;

 private:
  void SetCursor();
  void RestoreCursor();

  void BeginAutoCollapseTimer();
  void StopAutoCollapseTimer();

  void OnMouseLeave();

  ScrollBarBrushStateKind GetState(ScrollBarAreaKind area);

 protected:
  gsl::not_null<ScrollRenderObject*> render_object_;

  std::unique_ptr<platform::graphics::IGeometry> arrow_geometry_;

 private:
  Direction direction_;

  bool is_enabled_ = true;

  bool is_expanded_ = false;

  std::shared_ptr<platform::graphics::IBrush> collapsed_thumb_brush_;
  std::unordered_map<
      ScrollBarBrushUsageKind,
      std::unordered_map<ScrollBarBrushStateKind,
                         std::shared_ptr<platform::graphics::IBrush>>>
      brushes_;

  Rect move_thumb_thumb_original_rect_;
  std::optional<Point> move_thumb_start_;

  std::optional<ScrollBarAreaKind> mouse_hover_;
  std::optional<ScrollBarAreaKind> mouse_press_;

  EventRevokerListGuard event_guard_;

  Event<Scroll> scroll_attempt_event_;

  bool cursor_overrided_ = false;

  platform::gui::TimerAutoCanceler auto_collapse_timer_canceler_;
};

class CRU_UI_API HorizontalScrollBar : public ScrollBar {
 public:
  explicit HorizontalScrollBar(
      gsl::not_null<ScrollRenderObject*> render_object);

  CRU_DELETE_COPY(HorizontalScrollBar)
  CRU_DELETE_MOVE(HorizontalScrollBar)

  ~HorizontalScrollBar() override = default;

 protected:
  void DrawUpArrow(
      platform::graphics::IPainter* painter, const Rect& area,
      gsl::not_null<platform::graphics::IBrush*> arrow_brush,
      gsl::not_null<platform::graphics::IBrush*> background_brush) override;
  void DrawDownArrow(
      platform::graphics::IPainter* painter, const Rect& area,
      gsl::not_null<platform::graphics::IBrush*> arrow_brush,
      gsl::not_null<platform::graphics::IBrush*> background_brush) override;

  bool IsShowBar() override;

  std::optional<Rect> GetExpandedAreaRect(ScrollBarAreaKind area_kind) override;
  std::optional<Rect> GetCollapsedTriggerExpandAreaRect() override;
  std::optional<Rect> GetCollapsedThumbRect() override;

  float CalculateNewScrollPosition(const Rect& thumb_original_rect,
                                   const Point& mouse_offset) override;

  bool CanScrollUp() override;
  bool CanScrollDown() override;
};

class CRU_UI_API VerticalScrollBar : public ScrollBar {
 public:
  explicit VerticalScrollBar(gsl::not_null<ScrollRenderObject*> render_object);

  CRU_DELETE_COPY(VerticalScrollBar)
  CRU_DELETE_MOVE(VerticalScrollBar)

  ~VerticalScrollBar() override = default;

 protected:
  void DrawUpArrow(
      platform::graphics::IPainter* painter, const Rect& area,
      gsl::not_null<platform::graphics::IBrush*> arrow_brush,
      gsl::not_null<platform::graphics::IBrush*> background_brush) override;
  void DrawDownArrow(
      platform::graphics::IPainter* painter, const Rect& area,
      gsl::not_null<platform::graphics::IBrush*> arrow_brush,
      gsl::not_null<platform::graphics::IBrush*> background_brush) override;

  bool IsShowBar() override;

  std::optional<Rect> GetExpandedAreaRect(ScrollBarAreaKind area_kind) override;
  std::optional<Rect> GetCollapsedTriggerExpandAreaRect() override;
  std::optional<Rect> GetCollapsedThumbRect() override;

  float CalculateNewScrollPosition(const Rect& thumb_original_rect,
                                   const Point& mouse_offset) override;

  bool CanScrollUp() override;
  bool CanScrollDown() override;
};

// A delegate to draw scrollbar and register related events.
class CRU_UI_API ScrollBarDelegate : public Object {
 public:
  explicit ScrollBarDelegate(gsl::not_null<ScrollRenderObject*> render_object);

  CRU_DELETE_COPY(ScrollBarDelegate)
  CRU_DELETE_MOVE(ScrollBarDelegate)

  ~ScrollBarDelegate() override = default;

 public:
  bool IsHorizontalBarEnabled() const { return horizontal_bar_.IsEnabled(); }
  void SetHorizontalBarEnabled(bool value) {
    horizontal_bar_.SetEnabled(value);
  }

  bool IsVerticalBarEnabled() const { return horizontal_bar_.IsEnabled(); }
  void SetVerticalBarEnabled(bool value) { horizontal_bar_.SetEnabled(value); }

  IEvent<Scroll>* ScrollAttemptEvent() { return &scroll_attempt_event_; }

  void DrawScrollBar(platform::graphics::IPainter* painter);

  void InstallHandlers(controls::Control* control);
  void UninstallHandlers() { InstallHandlers(nullptr); }

 private:
  gsl::not_null<ScrollRenderObject*> render_object_;

  HorizontalScrollBar horizontal_bar_;
  VerticalScrollBar vertical_bar_;

  Event<Scroll> scroll_attempt_event_;
};
}  // namespace cru::ui::render

max_stars_repo_name: crupest/Cru

max_issues_repo_name: crupest/Cru

max_forks_repo_name: crupest/Cru


text:
#pragma once

#include <halley/utils/utils.h>
#include <vector>
#include <gsl/gsl>
#include "halley/data_structures/flat_map.h"
#include "halley/maths/vector2.h"
#include "halley/maths/rect.h"
#include "halley/file/path.h"
#include <map>
#include <unordered_map>
#include <cstdint>
#include <utility>
#include <set>
#include "halley/data_structures/maybe.h"
#include "halley/maths/vector4.h"

namespace Halley {
	class String;

	class Serializer {
	public:
		Serializer();
		explicit Serializer(gsl::span<gsl::byte> dst);

		template <typename T, typename std::enable_if<std::is_convertible<T, std::function<void(Serializer&)>>::value, int>::type = 0>
		static Bytes toBytes(const T& f)
		{
			Serializer dry;
			f(dry);
			Bytes result(dry.getSize());
			Serializer s(gsl::as_writeable_bytes(gsl::span<Halley::Byte>(result)));
			f(s);
			return result;
		}

		template <typename T, typename std::enable_if<!std::is_convertible<T, std::function<void(Serializer&)>>::value, int>::type = 0>
		static Bytes toBytes(const T& value)
		{
			return toBytes([&value](Serializer& s) { s << value; });
		}

		size_t getSize() const { return size; }

		Serializer& operator<<(bool val) { return serializePod(val); }
		Serializer& operator<<(int8_t val) { return serializePod(val); }
		Serializer& operator<<(uint8_t val) { return serializePod(val); }
		Serializer& operator<<(int16_t val) { return serializePod(val); }
		Serializer& operator<<(uint16_t val) { return serializePod(val); }
		Serializer& operator<<(int32_t val) { return serializePod(val); }
		Serializer& operator<<(uint32_t val) { return serializePod(val); }
		Serializer& operator<<(int64_t val) { return serializePod(val); }
		Serializer& operator<<(uint64_t val) { return serializePod(val); }
		Serializer& operator<<(float val) { return serializePod(val); }
		Serializer& operator<<(double val) { return serializePod(val); }

		Serializer& operator<<(const std::string& str);
		Serializer& operator<<(const String& str);
		Serializer& operator<<(const Path& path);
		Serializer& operator<<(gsl::span<const gsl::byte> span);
		Serializer& operator<<(const Bytes& bytes);

		template <typename T>
		Serializer& operator<<(const std::vector<T>& val)
		{
			unsigned int sz = static_cast<unsigned int>(val.size());
			*this << sz;
			for (unsigned int i = 0; i < sz; i++) {
				*this << val[i];
			}
			return *this;
		}

		template <typename T, typename U>
		Serializer& operator<<(const FlatMap<T, U>& val)
		{
			*this << static_cast<unsigned int>(val.size());
			for (auto& kv : val) {
				*this << kv.first << kv.second;
			}
			return *this;
		}

		template <typename T, typename U>
		Serializer& operator<<(const std::map<T, U>& val)
		{
			*this << static_cast<unsigned int>(val.size());
			for (auto& kv : val) {
				*this << kv.first << kv.second;
			}
			return *this;
		}

		template <typename T, typename U>
		Serializer& operator<<(const std::unordered_map<T, U>& val)
		{
			std::map<T, U> m;
			for (auto& kv: val) {
				m[kv.first] = kv.second;
			}
			return (*this << m);
		}

		template <typename T>
		Serializer& operator<<(const std::set<T>& val)
		{
			unsigned int sz = static_cast<unsigned int>(val.size());
			*this << sz;
			for (auto& v: val) {
				*this << v;
			}
			return *this;
		}

		template <typename T>
		Serializer& operator<<(const Vector2D<T>& val)
		{
			return *this << val.x << val.y;
		}

		template <typename T>
		Serializer& operator<<(const Vector4D<T>& val)
		{
			return *this << val.x << val.y << val.z << val.w;
		}

		template <typename T>
		Serializer& operator<<(const Rect2D<T>& val)
		{
			return *this << val.getTopLeft() << val.getBottomRight();
		}

		template <typename T, typename U>
		Serializer& operator<<(const std::pair<T, U>& p)
		{
			return *this << p.first << p.second;
		}

		template <typename T>
		Serializer& operator<<(const Maybe<T>& p)
		{
			if (p) {
				return *this << true << p.get();
			} else {
				return *this << false;
			}
		}

		template <typename T>
		Serializer& operator<<(const Range<T>& p)
		{
			return *this << p.start << p.end;
		}

		template <typename T, std::enable_if_t<std::is_enum<T>::value == true, int> = 0>
		Serializer& operator<<(const T& val)
		{
			using B = typename std::underlying_type<T>::type;
			return *this << B(val);
		}

		template <typename T, std::enable_if_t<std::is_enum<T>::value == false, int> = 0>
		Serializer& operator<<(const T& val)
		{
			val.serialize(*this);
			return *this;
		}

	private:
		bool dryRun;
		size_t size = 0;
		gsl::span<gsl::byte> dst;

		template <typename T>
		Serializer& serializePod(T val)
		{
			if (!dryRun) {
				memcpy(dst.data() + size, &val, sizeof(T));
			}
			size += sizeof(T);
			return *this;
		}
	};

	class Deserializer {
	public:
		Deserializer(gsl::span<const gsl::byte> src);
		explicit Deserializer(const Bytes& src);

		template <typename T>
		static T fromBytes(const Bytes& src)
		{
			T result;
			Deserializer s(src);
			s >> result;
			return result;
		}

		template <typename T>
		static T fromBytes(gsl::span<const gsl::byte> src)
		{
			T result;
			Deserializer s(src);
			s >> result;
			return result;
		}

		template <typename T>
		static void fromBytes(T& target, const Bytes& src)
		{
			Deserializer s(src);
			s >> target;
		}

		template <typename T>
		static void fromBytes(T& target, gsl::span<const gsl::byte> src)
		{
			Deserializer s(src);
			s >> target;
		}

		Deserializer& operator>>(bool& val) { return deserializePod(val); }
		Deserializer& operator>>(int8_t& val) { return deserializePod(val); }
		Deserializer& operator>>(uint8_t& val) { return deserializePod(val); }
		Deserializer& operator>>(int16_t& val) { return deserializePod(val); }
		Deserializer& operator>>(uint16_t& val) { return deserializePod(val); }
		Deserializer& operator>>(int32_t& val) { return deserializePod(val); }
		Deserializer& operator>>(uint32_t& val) { return deserializePod(val); }
		Deserializer& operator>>(int64_t& val) { return deserializePod(val); }
		Deserializer& operator>>(uint64_t& val) { return deserializePod(val); }
		Deserializer& operator>>(float& val) { return deserializePod(val); }
		Deserializer& operator>>(double& val) { return deserializePod(val); }

		Deserializer& operator>>(std::string& str);
		Deserializer& operator>>(String& str);
		Deserializer& operator>>(Path& p);
		Deserializer& operator>>(gsl::span<gsl::byte>& span);
		Deserializer& operator>>(Bytes& bytes);

		template <typename T>
		Deserializer& operator>>(std::vector<T>& val)
		{
			unsigned int sz;
			*this >> sz;
			ensureSufficientBytesRemaining(sz); // Expect at least one byte per vector entry

			val.clear();
			val.reserve(sz);
			for (unsigned int i = 0; i < sz; i++) {
				val.push_back(T());
				*this >> val[i];
			}
			return *this;
		}

		template <typename T, typename U>
		Deserializer& operator>>(FlatMap<T, U>& val)
		{
			unsigned int sz;
			*this >> sz;
			ensureSufficientBytesRemaining(sz * 2); // Expect at least two bytes per map entry

			std::vector<std::pair<T, U>> tmpData(sz);
			for (unsigned int i = 0; i < sz; i++) {
				*this >> tmpData[i].first >> tmpData[i].second;
			}
			val = FlatMap<T, U>(boost::container::ordered_unique_range_t(), tmpData.begin(), tmpData.end());
			return *this;
		}

		template <typename T, typename U>
		Deserializer& operator >> (std::map<T, U>& val)
		{
			unsigned int sz;
			*this >> sz;
			ensureSufficientBytesRemaining(sz * 2); // Expect at least two bytes per map entry

			for (unsigned int i = 0; i < sz; i++) {
				T key;
				U value;
				*this >> key >> value;
				val[key] = std::move(value);
			}
			return *this;
		}

		template <typename T, typename U>
		Deserializer& operator >> (std::unordered_map<T, U>& val)
		{
			unsigned int sz;
			*this >> sz;
			ensureSufficientBytesRemaining(sz * 2); // Expect at least two bytes per map entry

			for (unsigned int i = 0; i < sz; i++) {
				T key;
				U value;
				*this >> key >> value;
				val[key] = std::move(value);
			}
			return *this;
		}

		template <typename T>
		Deserializer& operator>>(std::set<T>& val)
		{
			unsigned int sz;
			*this >> sz;
			ensureSufficientBytesRemaining(sz); // Expect at least one byte per set entry

			val.clear();
			for (unsigned int i = 0; i < sz; i++) {
				T v;
				*this >> v;
				val.insert(std::move(v));
			}
			return *this;
		}

		template <typename T>
		Deserializer& operator>>(Vector2D<T>& val)
		{
			*this >> val.x;
			*this >> val.y;
			return *this;
		}

		template <typename T>
		Deserializer& operator>>(Vector4D<T>& val)
		{
			*this >> val.x;
			*this >> val.y;
			*this >> val.z;
			*this >> val.w;
			return *this;
		}

		template <typename T>
		Deserializer& operator>>(Rect2D<T>& val)
		{
			Vector2D<T> p1, p2;
			*this >> p1;
			*this >> p2;
			val = Rect2D<T>(p1, p2);
			return *this;
		}

		template <typename T, typename U>
		Deserializer& operator>>(std::pair<T, U>& p)
		{
			return *this >> p.first >> p.second;
		}

		template <typename T>
		Deserializer& operator>>(Range<T>& p)
		{
			return *this >> p.start >> p.end;
		}

		template <typename T>
		Deserializer& operator>>(Maybe<T>& p)
		{
			bool present;
			*this >> present;
			if (present) {
				T tmp;
				*this >> tmp;
				p = tmp;
			} else {
				p = Maybe<T>();
			}
			return *this;
		}

		template <typename T, std::enable_if_t<std::is_enum<T>::value == true, int> = 0>
		Deserializer& operator>>(T& val)
		{
			typename std::underlying_type<T>::type tmp;
			*this >> tmp;
			val = T(tmp);
			return *this;
		}

		template <typename T, std::enable_if_t<std::is_enum<T>::value == false, int> = 0>
		Deserializer& operator>>(T& val)
		{
			val.deserialize(*this);
			return *this;
		}

		void setVersion(int version);
		int getVersion() const;

	private:
		size_t pos = 0;
		gsl::span<const gsl::byte> src;
		int version = 0;

		template <typename T>
		Deserializer& deserializePod(T& val)
		{
			ensureSufficientBytesRemaining(sizeof(T));
			memcpy(&val, src.data() + pos, sizeof(T));
			pos += sizeof(T);
			return *this;
		}

		void ensureSufficientBytesRemaining(size_t bytes);
		size_t getBytesRemaining() const;
	};
}

max_stars_repo_name: sunhay/halley

max_issues_repo_name: sunhay/halley

max_forks_repo_name: sunhay/halley


text:
#ifndef OPENMC_TALLIES_FILTER_UNIVERSE_H
#define OPENMC_TALLIES_FILTER_UNIVERSE_H

#include <cstdint>
#include <unordered_map>

#include <gsl/gsl-lite.hpp>

#include "openmc/tallies/filter.h"
#include "openmc/vector.h"

namespace openmc {

//==============================================================================
//! Specifies which geometric universes tally events reside in.
//==============================================================================

class UniverseFilter : public Filter {
public:
  //----------------------------------------------------------------------------
  // Constructors, destructors

  ~UniverseFilter() = default;

  //----------------------------------------------------------------------------
  // Methods

  std::string type() const override { return "universe"; }

  void from_xml(pugi::xml_node node) override;

  void get_all_bins(const Particle& p, TallyEstimator estimator,
    FilterMatch& match) const override;

  void to_statepoint(hid_t filter_group) const override;

  std::string text_label(int bin) const override;

  //----------------------------------------------------------------------------
  // Accessors

  void set_universes(gsl::span<int32_t> universes);

private:
  //----------------------------------------------------------------------------
  // Data members

  //! The indices of the universes binned by this filter.
  vector<int32_t> universes_;

  //! A map from universe indices to filter bin indices.
  std::unordered_map<int32_t, int> map_;
};

} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_UNIVERSE_H

max_stars_repo_name: norberto-schmidt/openmc

max_issues_repo_name: norberto-schmidt/openmc

max_forks_repo_name: norberto-schmidt/openmc


text:
/*
 * author: Pierre Schnizer
 * created: June 2003
 */

#define PyGSL_SET_GSL_ERROR_HANDLER_OFF 1
#include <pygsl/error_helpers.h>
#include <pygsl/arrayobject.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_nan.h>
#undef _PyGSL_UFUNC_SUPPORTED
#ifdef PyGSL_NUMPY
#include <numpy/arrayobject.h>
#include <numpy/ufuncobject.h>
#define _PyGSL_UFUNC_SUPPORTED
#endif
#ifdef PyGSL_NUMERIC
#include <Numeric/arrayobject.h>
#include <Numeric/ufuncobject.h>
#define _PyGSL_UFUNC_SUPPORTED
#endif
#ifndef _PyGSL_UFUNC_SUPPORTED
#error "UFUNCs not supported by the selected/found array module!"
#endif
/*
   #include <pygsl/utils.h>
   #include <pygsl/error_helpers.h>
*/

static PyObject* module = NULL; /*used by the backtrace*/

#define IMPORTALL 1

/* I add the evaluate functions types by hand to have an extra check */
#ifdef IMPORTALL
#include "sf__arrays.c"

#undef SF_ARRAY
#endif

#ifdef IMPORTALL
#include "sf__evals.c"
#include "sf__data.c"
#endif

static void * polar_to_rect_data [] = { (void *) gsl_sf_polar_to_rect, (void *) gsl_sf_polar_to_rect};
static void * rect_to_polar_data [] = { (void *) gsl_sf_rect_to_polar, (void *) gsl_sf_rect_to_polar};

static char PyGSL_sf_ufunc_qi_dd_D_one_types [] = { PyArray_FLOAT,  PyArray_FLOAT,  PyArray_CFLOAT,
                                                    PyArray_DOUBLE, PyArray_DOUBLE, PyArray_CDOUBLE };

static char PyGSL_sf_ufunc_qi_D_dd_one_types [] = { PyArray_CFLOAT,  PyArray_FLOAT,  PyArray_FLOAT,
                                                    PyArray_CDOUBLE, PyArray_DOUBLE, PyArray_DOUBLE };

typedef int PyGSL_sf_ufunc_qi_dd_D_one(double, double, gsl_sf_result *, gsl_sf_result *);
void  PyGSL_sf_ufunc_qi_dd_D(char **args, int *dimensions, int *steps, void *func)
{
        int is0=steps[0], is1=steps[1], os0=steps[2];
        char   *ip0=args[0], *ip1=args[1], *op0=args[2];
	double *dop;

	gsl_sf_result x, y;
	int tmp;
	PyGSL_array_index_t i;
	FUNC_MESS_BEGIN();
	DEBUG_MESS(2, "dimensions = %d %d %d", dimensions[0], dimensions[1],dimensions[2]);
	DEBUG_MESS(2, "steps = %d %d %d", steps[0], steps[1], steps[2]);
	DEBUG_MESS(2, "args = %p %p %p", args[0], args[1], args[2]);
        for(i = 0; i<dimensions[0]; i++, ip0+=is0, ip1+=is1, op0+=os0){
	     DEBUG_MESS(2, "i = %ld", (long)i);
	    dop = ((double *) op0);
	    tmp = ((PyGSL_sf_ufunc_qi_dd_D_one *) func)(*((double *)ip0), *((double *)ip1), &x, &y);
	    dop[0] = x.val;
	    dop[1] = y.val;
        }
	FUNC_MESS_END();
}
void  PyGSL_sf_ufunc_qi_dd_D_as_ff_F(char **args, int *dimensions, int *steps, void *func)
{
        int i,  is0=steps[0], is1=steps[1], os0=steps[2], tmp;
        char   *ip0=args[0], *ip1=args[1], *op0=args[2];
	gsl_sf_result x, y;

	FUNC_MESS_BEGIN();
	DEBUG_MESS(2, "dimensions = %d %d %d", dimensions[0], dimensions[1],dimensions[2]);
	DEBUG_MESS(2, "steps = %d %d %d", steps[0], steps[1], steps[2]);
	DEBUG_MESS(2, "args = %p %p %p", args[0], args[1], args[2]);
        for(i = 0; i<dimensions[0]; i++, ip0+=is0, ip1+=is1, op0+=os0){
	     DEBUG_MESS(2, "i = %d", i);
	     tmp = ((PyGSL_sf_ufunc_qi_dd_D_one *) func)((double) *((float *)ip0), (double) *((float *)ip1), &x, &y);
	     *((float *)op0) = (float) x.val;
	     *(((float *)op0)+1) = (float) y.val;
        }
	FUNC_MESS_END();
}
static PyUFuncGenericFunction PyGSL_sf_ufunc_qi_dd_D_one_data[] = {PyGSL_sf_ufunc_qi_dd_D_as_ff_F, PyGSL_sf_ufunc_qi_dd_D};

typedef int PyGSL_sf_ufunc_qi_D_dd_one(double, double, gsl_sf_result *, gsl_sf_result *);
void  PyGSL_sf_ufunc_qi_D_dd(char **args, int *dimensions, int *steps, void *func)
{
	int tmp;
	PyGSL_array_index_t i,  is0=steps[0], os0=steps[1], os1=steps[2];
        char   *ip0=args[0], *op0=args[1], *op1=args[2];
	double *dip;
	gsl_sf_result radius, angle;

	FUNC_MESS_BEGIN();
	DEBUG_MESS(2, "dimensions = %d %d %d", dimensions[0], dimensions[1],dimensions[2]);
	DEBUG_MESS(2, "steps = %d %d %d", steps[0], steps[1], steps[2]);
	DEBUG_MESS(2, "args = %p %p %p", args[0], args[1], args[2]);
	DEBUG_MESS(1, "rect_to_polar %p", (void *)func);
        for(i = 0; i<dimensions[0]; i++, ip0+=is0, op0+=os0, op1+=os1){
	     DEBUG_MESS(2, "i = %ld", (long)i);
	    dip = ((double *) ip0);
	    tmp = ((PyGSL_sf_ufunc_qi_D_dd_one *) func)(dip[0], dip[1], &radius, &angle);
	    *((double *)op0) = radius.val;
	    *((double *)op1) = angle.val;
        }
	FUNC_MESS_END();
}
void  PyGSL_sf_ufunc_qi_D_dd_as_F_ff(char **args, int *dimensions, int *steps, void *func)
{
        PyGSL_array_index_t i,  is0=steps[0], os0=steps[1],  os1=steps[2], tmp;
        char   *ip0=args[0],  *op0=args[1], *op1=args[2];
	gsl_sf_result radius, angle;
	FUNC_MESS_BEGIN();
	DEBUG_MESS(2, "dimensions = %d %d %d", dimensions[0], dimensions[1],dimensions[2]);
	DEBUG_MESS(2, "steps = %d %d %d", steps[0], steps[1], steps[2]);
	DEBUG_MESS(2, "args = %p %p %p", args[0], args[1], args[2]);
	DEBUG_MESS(1, "polar_to_rect %p", (void *)func);
        for(i = 0; i<dimensions[0]; i++, ip0+=is0, op0+=os0, op1+=os1){
	     tmp = ((PyGSL_sf_ufunc_qi_D_dd_one *) func)((double) *((float *)ip0), (double) *(((float *)ip0)+1), &radius, &angle);
	     DEBUG_MESS(2, "i = %ld", (long)i);
	      *((float *)op0) = (float) radius.val;
	      *((float *)op1) = (float) angle.val;
	}
	FUNC_MESS_END();
}
static PyUFuncGenericFunction PyGSL_sf_ufunc_qi_D_dd_one_data[] = {PyGSL_sf_ufunc_qi_D_dd_as_F_ff, PyGSL_sf_ufunc_qi_D_dd};

static char *PyGSL_polar_to_rect_doc = "Convert r, theta into a complex representation";
static char *PyGSL_rect_to_polar_doc = "Convert a complex into r and theta. Theta will be in the range [-pi,pi]";
static char * sf_module_doc =
"GSL Special Functions wrapped mostly as Numeric ufuncs. \n\
Use print sf.<function>.__doc__ to get more information about the ufunc.\n\
Help does not display the information stored there!\n";

DL_EXPORT(void) init_ufuncs(void)
{
  PyObject* sf_module=NULL, *sf_dict=NULL;
  PyObject* f=NULL, *item=NULL;

  FUNC_MESS_BEGIN();
  DEBUG_MESS(2, "Module compiled at %s %s\n", __DATE__, __TIME__);
  sf_module=Py_InitModule("_ufuncs", sf_array_functions);
  module = sf_module;

  import_array();
  import_ufunc();
  init_pygsl();

  /* here is an error check needed! */
  /*
  gsl_errors_module=PyImport_ImportModule ("pygsl.errors");

  gsl_errors_dict=PyModule_GetDict(gsl_errors_module);
  gsl_module_error=PyDict_GetItemString(gsl_errors_dict,"gsl_Error");
  */

  sf_dict=PyModule_GetDict(sf_module);
  if (!(item = PyString_FromString(sf_module_doc))){
       PyErr_SetString(PyExc_ImportError,
		       "I could not generate module doc string!");
       goto fail;
  }
  if (PyDict_SetItemString(sf_dict, "__doc__", item) != 0){
       PyErr_SetString(PyExc_ImportError,
		       "I could not add doc string to module dict!");
       goto fail;
  }

  PyDict_SetItemString(sf_dict, "PREC_DOUBLE", PyInt_FromLong((long) GSL_PREC_DOUBLE));
  PyDict_SetItemString(sf_dict, "PREC_SINGLE", PyInt_FromLong((long) GSL_PREC_SINGLE));
  PyDict_SetItemString(sf_dict, "PREC_APPROX", PyInt_FromLong((long) GSL_PREC_APPROX));

  DEBUG_MESS(1, "polar_to_rect %p", (void *)gsl_sf_polar_to_rect);
  f = PyUFunc_FromFuncAndData(PyGSL_sf_ufunc_qi_dd_D_one_data,
			      polar_to_rect_data,
			      PyGSL_sf_ufunc_qi_dd_D_one_types,
			      2, /* number of supported types */
			      2, /* in args */
			      1, /* out args */
			      PyUFunc_None,
			      "polar_to_rect",
			      PyGSL_polar_to_rect_doc,
			      0 /*check return*/);
  PyDict_SetItemString(sf_dict, "polar_to_rect", f);
  /* Py_DECREF(f); */

  DEBUG_MESS(1, "rect_to_polar %p", (void *)gsl_sf_rect_to_polar);
  f = PyUFunc_FromFuncAndData(PyGSL_sf_ufunc_qi_D_dd_one_data,
			      rect_to_polar_data,
			      PyGSL_sf_ufunc_qi_D_dd_one_types,
			       2, /* number of supported types */
			       1, /* in args */
			       2, /* out args */
			       PyUFunc_None,
			       "rect_to_polar",
			       PyGSL_rect_to_polar_doc,
			       0 /*check return*/);
  PyDict_SetItemString(sf_dict, "rect_to_polar", f);
  /* Py_DECREF(f); */

#ifdef IMPORTALL
#include "sf__objects.c"
#endif

  FUNC_MESS_END();
  return;

 fail:
     if(!PyErr_Occurred()){
	  PyErr_SetString(PyExc_ImportError, "I could not init sf module!");
     }
}

max_stars_repo_name: juhnowski/FishingRod

max_issues_repo_name: juhnowski/FishingRod

max_forks_repo_name: juhnowski/FishingRod


text:
#include <stdio.h>
#include <gsl/gsl_sf_bessel.h>
int main(void)
{
  double x = 5.0;
  printf("hello\n");
  double y = gsl_sf_bessel_J0(x);
  printf("J0(%g) = %.18e\n", x, y);
#ifdef WIN32
  system("pause");
#endif
  return 0;
}


text:
// sf: special functions

// == Literature ==
// [1] Richard P. Brent: "Fast Algorithms for High-Precision
// Computation of Elementary Functions", 2006
// [2] Semjon Adlaj: "An eloquent formula for the perimeter
// of an ellipse", Notices of the AMS 59(8) (2012)
// [3] Semjon Adlaj: "An arithmetic-geometric mean of a
// third kind!", 2015

#include <math.h>
#include <complex.h>
#include <stdio.h>
#include <stdlib.h>
#include <moss.h>
// #include <gsl/gsl_sf.h>

#define PI 3.14159265358979323846
#define SQRT_PI 1.772453850905516027

#ifndef CMPLX
    #define CMPLX(x,y) ((x)+(y)*I)
#endif

double mf_float(mt_object* x, int* error);
int mf_complex(mt_object* z, mt_object* x);

static double hzeta(double s, double a);
static double complex chzeta(double complex s, double complex a);
static
double sgngamma(double x){
    if(x>0) return 1;
    return copysign(1,sin(PI*x));
}

static
double binomial_coefficient(double n, double k){
    double s = sgngamma(n+1)*sgngamma(k+1)*sgngamma(n-k+1);
    return s*exp(lgamma(n+1)-lgamma(k+1)-lgamma(n-k+1));
}

// Lanczos approximation
static
double complex cgammapz(double complex z){
    static const double p[] = {
        0.99999999999980993, 676.5203681218851, -1259.1392167224028,
        771.32342877765313, -176.61502916214059, 12.507343278686905,
        -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7
    };
    z-=1;
    double complex y = p[0];
    int i;
    for(i=1; i<9; i++){
        y+=p[i]/(z+i);
    }
    double complex t = z+7.5;
    return sqrt(2*PI)*cpow(t,(z+0.5))*cexp(-t)*y;
}

// Use Euler's reflection formula
double complex mf_cgamma(double complex z){
    if(creal(z)<0.5){
        return PI/(csin(PI*z)*cgammapz(1-z));
    }else{
        return cgammapz(z);
    }
}

static
double Beta(double x, double y){
    return(sgngamma(x)*sgngamma(y)*sgngamma(x+y)
        *exp(lgamma(x)+lgamma(y)-lgamma(x+y))
    );
}

// Use Euler-Maclaurin summation formula applied
// to a partial sum of the harmonic series,
// reflection formula psi(x)=psi(1-x)-pi*cot(pi*x)
// and functional equation psi(x)=psi(x+1)-1/x.
static
double digamma_positive(double x){
    double s = 0;
    while(x<10){
        s-=1/x; x+=1;
    }
    double x2,x4,x8,x10;
    x2 = x*x; x4 = x2*x2; x8 = x4*x4; x10 = x8*x2;
    return s+(log(x)-0.5/x-1/(12*x2)
        +1/(120*x4)-1/(252*x4*x2)+1/(240*x8)
        -5/(660*x10)+691/(32760*x10*x2)
    );
}

static
double digamma(double x){
    if(x<0){
        return digamma_positive(1-x)-PI/tan(PI*x);
    }else{
        return digamma_positive(x);
    }
}

static
double complex cdigamma_positive(double complex x){
    double complex s=0;
    while(creal(x)<10){
        s-=1/x; x+=1;
    }
    double complex x2,x4,x8,x10;
    x2 = x*x; x4 = x2*x2; x8 = x4*x4; x10 = x8*x2;
    return s+(clog(x)-0.5/x-1/(12*x2)
        +1/(120*x4)-1/(252*x4*x2)+1/(240*x8)
        -5/(660*x10)+691/(32760*x10*x2)
    );
}

static
double complex cdigamma(double complex x){
    if(creal(x)<0){
        return cdigamma_positive(1-x)-PI/ctan(PI*x);
    }else{
        return cdigamma_positive(x);
    }
}

static
double polygamma(long n, double x){
    if(n>0){
        return cos(PI*(n+1))*tgamma(n+1)*hzeta(n+1,x);
    }else if(n==0){
        return digamma(x);
    }else{
        return NAN;
    }
}

static
double complex cpolygamma(long n, double complex z){
    if(n>0){
        return cos(PI*(n+1))*tgamma(n+1)*chzeta(n+1,z);
    }else if(n==0){
        return cdigamma(z);
    }else{
        return NAN+NAN*I;
    }
}

double cfGamma(double s, double z, int n){
    double x = 0;
    int k;
    for(k=n; k>0; k--){
        x = k*(s-k)/(2*k+1+z-s+x);
    }
    return exp(-z)/(1+z-s+x);
}

double complex ccfGamma(double complex s, double complex z, int n){
    double complex x = 0;
    int k;
    for(k=n; k>0; k--){
        x = k*(s-k)/(2*k+1+z-s+x);
    }
    return cexp(-z)/(1+z-s+x);
}

double sgamma(double s, double z, int n){
    double y = 0;
    double p = 1/s;
    int k;
    for(k=1; k<=n; k++){
        y+=p;
        p = p*z/(s+k);
    }
    return y*exp(-z);
}

double complex csgamma(double complex s, double complex z, int n){
    double complex y = 0;
    double complex p = 1/s;
    int k;
    for(k=1; k<=n; k++){
        y+=p;
        p = p*z/(s+k);
    }
    return y*cexp(-z);
}

static
double incGamma(double s, double z){
    if(s+4<z){
        return pow(z,s)*cfGamma(s,z,40);
    }else{
        return tgamma(s)-pow(z,s)*sgamma(s,z,60);
    }
}

static
double incgamma(double s, double z){
    if(s+4<z){
        return tgamma(s)-pow(z,s)*cfGamma(s,z,40);
    }else{
        return pow(z,s)*sgamma(s,z,60);
    }
}

static
double complex cincGamma(double complex s, double complex z){
    if(creal(s)+4<creal(z)){
        return cpow(z,s)*ccfGamma(s,z,40);
    }else{
        return mf_cgamma(s)-cpow(z,s)*csgamma(s,z,60);
    }
}

static
double complex cincgamma(double complex s, double complex z){
    if(creal(s)+4<creal(z)){
        return mf_cgamma(s)-cpow(z,s)*ccfGamma(s,z,40);
    }else{
        return cpow(z,s)*csgamma(s,z,60);
    }
}

// Use Euler-Maclaurin summation formula
static
double hzeta(double s, double a){
    double y = 0;
    int N = 18;
    double Npa = N+a;
    int k;
    for(k=0; k<N; k++){
        y+=pow(k+a,-s);
    }
    double s2 = s*(s+1)*(s+2);
    double s4 = s2*(s+3)*(s+4);
    y+=pow(Npa,1-s)/(s-1)+0.5*pow(Npa,-s);
    y+=s*pow(Npa,-s-1)/12;
    y-=s2*pow(Npa,-s-3)/720;
    y+=s4*pow(Npa,-s-5)/30240;
    y-=s4*(s+5)*(s+6)*pow(Npa,-s-7)/1209600;
    return y;
}

// Use reflection formula (Riemann functional equation)
static
double zeta(double s){
    if(s>-1){
        return hzeta(s,1);
    }else{
        double a = 2*pow(2*PI,s-1)*sin(0.5*PI*s);
        return a*tgamma(1-s)*hzeta(1-s,1);
    }
}

static
double complex chzeta(double complex s, double complex a){
    double complex y = 0;
    int N = 18;
    double complex Npa = N+a;
    int k;
    for(k=0; k<N; k++){
        y+=cpow(k+a,-s);
    }
    double complex s2 = s*(s+1)*(s+2);
    double complex s4 = s2*(s+3)*(s+4);
    y+=cpow(Npa,1-s)/(s-1)+0.5*cpow(Npa,-s);
    y+=s*cpow(Npa,-s-1)/12;
    y-=s2*cpow(Npa,-s-3)/720;
    y+=s4*cpow(Npa,-s-5)/30240;
    y-=s4*(s+5)*(s+6)*cpow(Npa,-s-7)/1209600;
    return y;
}

static
double complex czeta(double complex s){
    if(creal(s)>-1){
        return chzeta(s,1);
    }else{
        double complex a = 2*cpow(2*PI,s-1)*csin(PI*s/2);
        return a*mf_cgamma(1-s)*chzeta(1-s,1);
    }
}

static
double ALP_iter(long n, long m, double x){
    if(n==m){
        return SQRT_PI/tgamma(0.5-n)*pow(2*sqrt(1-x*x),n);
    }else if(n-1==m){
        return x*(2.0*n-1)*ALP_iter(m,m,x);
    }else{
        double a = ALP_iter(m,m,x);
        double b = ALP_iter(m+1,m,x);
        double h;
        long k;
        for(k=m+2; k<=n; k++){
            h = ((2.0*k-1)*x*b-(k-1.0+m)*a)/(k-m);
            a = b;
            b = h;
        }
        return b;
    }
}

static
double ALP(long n, long m, double x){
    if(n<0) n = -n-1;
    if(labs(m)>n){
        return 0;
    }else if(m<0){
        return((m%2==0?1:-1)
            *exp(lgamma(n+m+1)-lgamma(n-m+1))
            *ALP_iter(n,-m,x)
        );
    }else{
        return ALP_iter(n,m,x);
    }
}

static
double Hermite(long n, double x){
    if(n==0){
        return 1;
    }else if(n==1){
        return 2*x;
    }else{
        double a,b,h;
        long k;
        a=1; b=2*x;
        for(k=1; k<n; k++){
            h = 2*x*b-2*k*a;
            a=b; b=h;
        }
        return b;
    }
}

static
double DHermite(long n, long m, double x){
    if(m==0){
        return Hermite(n,x);
    }else if(n<m){
        return 0;
    }else if(m==1){
        return 2*n*Hermite(n-1,x);
    }else{
        return 2*n*DHermite(n-1,m-1,x);
    }
}

static
double PT(long n, double x){
    if(n==0){
        return 1;
    }else if(n==1){
        return x;
    }else{
        double a,b,h;
        a=1; b=x;
        long k;
        for(k=2; k<=n; k++){
            h = 2*x*b-a;
            a=b; b=h;
        }
        return b;
    }
}

static
double PU(long n, double x){
    if(n==0){
        return 1;
    }else if(n==1){
        return 2*x;
    }else{
        double a,b,h;
        a=1; b=2*x;
        long k;
        for(k=2; k<=n; k++){
            h = 2*x*b-a;
            a=b; b=h;
        }
        return b;
    }
}

// Arithmetic-geometric mean
static
double agm(double a, double b){
    double ah,bh;
    int i;
    for(i=0; i<20; i++){
        ah = (a+b)/2;
        bh = sqrt(a*b);
        a=ah; b=bh;
        if(fabs(a-b)<1E-15) break;
    }
    return a;
}

static
double complex cagm(double complex a, double complex b){
    double complex ah,bh;
    int i;
    for(i=0; i<20; i++){
        ah = (a+b)/2;
        bh = csqrt(a*b);
        a=ah; b=bh;
        if(cabs(a-b)<1E-15) break;
    }
    return a;
}

// Modified arithmetic-geometric mean, see
// Semjon Adlaj: "An eloquent formula for the perimeter
// of an ellipse", Notices of the AMS 59(8) (2012), p. 1094-1099
double magm(double x, double y){
    double z=0;
    double xh,yh,zh,r;
    int i;
    for(i=0; i<20; i++){
        xh = 0.5*(x+y);
        r = sqrt((x-z)*(y-z));
        yh = z+r;
        zh = z-r;
        x=xh; y=yh; z=zh;
        if(fabs(x-y)<2E-15) break;
    }
    return x;
}

// Branch cut Re(z)>0, Im(z)=0
double complex csqrta(double complex z){
    return cimag(z)>=0? csqrt(z): -csqrt(z);
}

// Branch cut Re(z)>0, Im(z)=0
double complex csqrtb(double complex z){
    return cimag(z)<=0? csqrt(z): -csqrt(z);
}

double complex cmagm(double complex x, double complex y){
    double complex z=0;
    double complex xh,yh,zh,r;
    double complex (*root)(double complex);
    if(signbit(cimag(y))){
        root = csqrtb;
    }else{
        root = csqrta;
    }

    int i;
    for(i=0; i<20; i++){
        xh = 0.5*(x+y);
        r = root((x-z)*(y-z));
        yh = z+r;
        zh = z-r;
        x=xh; y=yh; z=zh;
        if(cabs(x-y)<2E-15) break;
    }
    return x;
}

static
double eiK(double m){
    return 0.5*PI/agm(1,sqrt(1-m));
}

static
double complex ceiK(double complex m){
    return 0.5*PI/cagm(1,csqrt(1-m));
}

double eiE(double m){
    return 0.5*PI*magm(1,1-m)/agm(1,sqrt(1-m));
}

double complex ceiE(double complex m){
    if(m==1) return 1;
    return 0.5*PI*cmagm(1,1-m)/cagm(1,csqrt(1-m));
}

static
double Fmn(long m, long n, double* a, double* b, double x){
    long i,k;
    double s=1,p=1;
    for(k=0; k<100; k++){
        p = p*x/(k+1);
        for(i=0; i<m; i++) p = p*(a[i]+k);
        for(i=0; i<n; i++) p = p/(b[i]+k);
        s+=p;
    }
    return s;
}

static
double complex cerf_series(double complex z, int N){
    double complex y=z;
    double p=1;
    double m;
    int n;
    for(n=1; n<=N; n++){
        p = -p*n;
        m = 2*n+1;
        y+=cpow(z,m)/(m*p);
    }
    return y*2/SQRT_PI;
}

static
double complex cerf_cf(double complex z, int N){
    double complex t=0;
    int k;
    for(k=N; k>=1; k--){
        t = k/((k%2==0?z:2*z)+t);
    }
    return 1-cexp(-z*z)/(SQRT_PI*(z+t));
}

static
double complex cerf_cf2(double complex z, int N){
    double complex t=1, z2=z*z;
    int k;
    for(k=N; k>=1; k--){
        t = (2*k-1)+(k%2==0?2:-2)*k*z2/t;
    }
    return 2*z/SQRT_PI*cexp(-z2)/t;
}

static
double complex mf_cerf(double complex z){
    double s;
    if(creal(z)<0){
        z=-z; s=-1;
    }else{
        s=1;
    }
    if(cabs(z)<2){
        return s*cerf_series(z,30);
    }else if(creal(z)<2 && fabs(cimag(z))<6){
        return s*cerf_cf2(z,60);
    }else{
        return s*cerf_cf(z,40);
    }
}

static
int sf_bc(mt_object* x, int argc, mt_object* v){
    if(argc!=2){
        mf_argc_error(argc,2,2,"bc");
        return 1;
    }
    int e=0;
    double n = mf_float(v+1,&e);
    if(e){
        mf_type_error1("in bc(n,k): cannot convert n (type: %s) to float.",v+1);
        return 1;
    }
    double k = mf_float(v+2,&e);
    if(e){
        mf_type_error1("in bc(n,k): cannot convert k (type: %s) to float.",v+2);
        return 1;
    }
    x->type = mv_float;
    x->value.f = binomial_coefficient(n,k);
    return 0;
}

static
int sf_gamma(mt_object* x, int argc, mt_object* v){
    if(argc==1){
        if(v[1].type==mv_complex){
            double complex w = mf_cgamma(CMPLX(v[1].value.c.re,v[1].value.c.im));
            x->type = mv_complex;
            x->value.c.re = creal(w);
            x->value.c.im = cimag(w);
            return 0;
        }else{
            int e=0;
            double t = mf_float(v+1,&e);
            if(e){
                mf_type_error1("in sf.gamma(z): cannot convert z (type: %s) to complex.",v+1);
                return 1;
            }
            x->type = mv_float;
            x->value.f = tgamma(t);
            return 0;
        }
    }else if(argc==2){
        mt_object s,z;
        if(v[1].type==mv_complex || v[2].type==mv_complex){
            mf_complex(&s,v+1);
            mf_complex(&z,v+2);
            double complex w = cincgamma(
                CMPLX(s.value.c.re,s.value.c.im),
                CMPLX(z.value.c.re,z.value.c.im)
            );
            x->type = mv_complex;
            x->value.c.re = creal(w);
            x->value.c.im = cimag(w);
            return 0;
        }else{
            int e=0;
            s.value.f = mf_float(v+1,&e);
            z.value.f = mf_float(v+2,&e);
            x->type = mv_float;
            x->value.f = incgamma(s.value.f,z.value.f);
            return 0;
        }
    }else{
        mf_argc_error(argc,1,2,"sf.gamma");
        return 1;
    }
}

static
int sf_Gamma(mt_object* x, int argc, mt_object* v){
    if(argc==2){
        mt_object s,z;
        if(v[1].type==mv_complex || v[2].type==mv_complex){
            mf_complex(&s,v+1);
            mf_complex(&z,v+2);
            double complex w = cincGamma(
                CMPLX(s.value.c.re,s.value.c.im),
                CMPLX(z.value.c.re,z.value.c.im)
            );
            x->type = mv_complex;
            x->value.c.re = creal(w);
            x->value.c.im = cimag(w);
            return 0;
        }else{
            int e=0;
            s.value.f = mf_float(v+1,&e);
            z.value.f = mf_float(v+2,&e);
            x->type = mv_float;
            x->value.f = incGamma(s.value.f,z.value.f);
            return 0;
        }
    }else{
        mf_argc_error(argc,2,2,"sf.Gamma");
        return 1;
    }
}

int sf_Beta(mt_object* x, int argc, mt_object* v){
    if(argc!=2){
        mf_argc_error(argc,2,2,"Beta");
        return 1;
    }
    int e=0;
    double a = mf_float(v+1,&e);
    if(e){
        mf_type_error1("in Beta(x,y): cannot convert x (type: %s) to flaot.",v+1);
        return 1;
    }
    double b = mf_float(v+2,&e);
    if(e){
        mf_type_error1("in Beta(x,y): cannot convert y (type: %s) to flaot.",v+2);
        return 1;
    }
    x->type = mv_float;
    x->value.f = Beta(a,b);
    return 0;
}

static
int sf_digamma(mt_object* x, int argc, mt_object* v){
    if(argc!=1){
        mf_argc_error(argc,1,1,"digamma");
        return 1;
    }
    if(v[1].type==mv_complex){
        double complex w = cdigamma(CMPLX(v[1].value.c.re,v[1].value.c.im));
        x->type = mv_complex;
        x->value.c.re = creal(w);
        x->value.c.im = cimag(w);
        return 0;
    }else{
        int e=0;
        double t = mf_float(v+1,&e);
        if(e){
            mf_type_error1("in digamma(x): cannot convert x (type: %s) to complex.",v+1);
            return 1;
        }
        x->type = mv_float;
        x->value.f = digamma(t);
        return 0;
    }
}

static
int sf_polygamma(mt_object* x, int argc, mt_object* v){
    if(argc!=2){
        mf_argc_error(argc,2,2,"polygamma");
        return 1;
    }
    if(v[1].type!=mv_int){
        mf_type_error1("in polygamma(n,z): n (type: %s) is not an integer.",v+1);
        return 1;
    }
    long n = v[1].value.i;
    if(v[2].type==mv_complex){
        double complex w = cpolygamma(n,CMPLX(v[2].value.c.re,v[2].value.c.im));
        x->type = mv_complex;
        x->value.c.re = creal(w);
        x->value.c.im = cimag(w);
        return 0;
    }else{
        int e=0;
        double t = mf_float(v+2,&e);
        if(e){
            mf_type_error1("in polygamma(n,z): z (type: %s) is not a complex number.",v+2);
            return 1;
        }
        x->type = mv_float;
        x->value.f = polygamma(n,t);
        return 0;
    }
}

static
int sf_zeta(mt_object* x, int argc, mt_object* v){
    int e;
    if(argc==1){
        if(v[1].type==mv_complex){
            double complex w = czeta(CMPLX(v[1].value.c.re,v[1].value.c.im));
            x->type = mv_complex;
            x->value.c.re = creal(w);
            x->value.c.im = cimag(w);
            return 0;
        }else{
            e=0;
            double t = mf_float(v+1,&e);
            if(e){
                mf_type_error1("in zeta(s): cannot convert s (type: %s) to complex.",v+1);
                return 1;
            }
            x->type = mv_float;
            x->value.f = zeta(t);
            return 0;
        }
    }else if(argc==2){
        mt_object z,a;
        if(v[1].type==mv_complex || v[2].type==mv_complex){
            mf_complex(&z,v+1);
            mf_complex(&a,v+2);
            double complex w = chzeta(
                CMPLX(z.value.c.re,z.value.c.im),
                CMPLX(a.value.c.re,a.value.c.im)
            );
            x->type = mv_complex;
            x->value.c.re = creal(w);
            x->value.c.im = cimag(w);
            return 0;
        }else{
            e=0;
            z.value.f = mf_float(v+1,&e);
            a.value.f = mf_float(v+2,&e);
            x->type = mv_float;
            x->value.f = hzeta(z.value.f,a.value.f);
            return 0;
        }
    }else{
        mf_argc_error(argc,1,2,"zeta");
        return 1;
    }
}

static
int sf_BernoulliB(mt_object* x, int argc, mt_object* v){
    if(argc!=1){
        mf_argc_error(argc,1,1,"sf.B");
        return 1;
    }
    int e=0;
    double t = mf_float(v+1,&e);
    if(e){
        mf_type_error1("in sf.B(n): cannot convert n (type: %s) to float.",v+1);
        return 1;
    }
    x->type=mv_float;
    if(t==0){
        x->value.f = 1;
    }else{
        x->value.f = -t*zeta(1-t);
    }
    return 0;
}

static
int sf_PP(mt_object* x, int argc, mt_object* v){
    if(argc==2){
        if(v[1].type!=mv_int){
            mf_type_error1("in sf.PP(n,x): n (type: %s) is not an integer.",v+1);
            return 1;
        }
        int n = (int)v[1].value.i;
        int e=0;
        double t = mf_float(v+2,&e);
        if(e){
            mf_type_error1("in sf.PP(n,x): cannot convert x (type: %s) to float.",v+2);
            return 1;
        }
        x->type = mv_float;
        x->value.f = ALP(n,0,t);
        return 0;
    }else if(argc==3){
        if(v[1].type!=mv_int){
            mf_type_error1("in sf.PP(n,m,x): n (type: %s) is not an integer.",v+1);
            return 1;
        }
        if(v[2].type!=mv_int){
            mf_type_error1("in sf.PP(n,m,x): m (type: %s) is not an integer.",v+2);
            return 1;
        }
        long n = v[1].value.i;
        long m = v[2].value.i;
        int e=0;
        double t = mf_float(v+3,&e);
        if(e){
            mf_type_error1("in sf.PP(n,m,x): cannot convert x (type: %s) to float.",v+3);
            return 1;
        }
        x->type = mv_float;
        x->value.f = ALP(n,m,t);
        return 0;
    }else{
        mf_argc_error(argc,2,3,"sf.PP");
        return 1;
    }
}

static
int sf_PH(mt_object* x, int argc, mt_object* v){
    if(argc==2){
        if(v[1].type!=mv_int){
            mf_type_error1("in sf.PH(n,x): n (type: %s) is not an integer.",v+1);
            return 1;
        }
        int n = (int)v[1].value.i;
        int e=0;
        double t = mf_float(v+2,&e);
        if(e){
            mf_type_error1("in sf.PH(n,x): cannot convert x (type: %s) to float.",v+2);
            return 1;
        }
        x->type = mv_float;
        x->value.f = Hermite(n,t);
        return 0;
    }else if(argc==3){
        if(v[1].type!=mv_int){
            mf_type_error1("in sf.PH(n,m,x): n (type: %s) is not an integer.",v+1);
            return 1;
        }
        if(v[2].type!=mv_int){
            mf_type_error1("in sf.PH(n,m,x): m (type: %s) is not an integer.",v+2);
            return 1;
        }
        long n = v[1].value.i;
        long m = v[2].value.i;
        int e=0;
        double t = mf_float(v+3,&e);
        if(e){
            mf_type_error1("in sf.PH(n,m,x): cannot convert x (type: %s) to float.",v+3);
            return 1;
        }
        x->type = mv_float;
        x->value.f = DHermite(n,m,t);
        return 0;
    }else{
        mf_argc_error(argc,2,3,"sf.PH");
        return 1;
    }
}

static
int sf_PT(mt_object* x, int argc, mt_object* v){
    if(argc!=2){
        mf_argc_error(argc,2,2,"sf.PT");
        return 1;
    }
    if(v[1].type!=mv_int){
        mf_type_error1("in sf.PT(n,x): n (type: %s) is not an integer.",v+1);
        return 1;
    }
    long n = v[1].value.i;
    int e=0;
    double t = mf_float(v+2,&e);
    if(e){
        mf_type_error1("in sf.PT(n,x): cannot convert x (type: %s) to float.",v+2);
        return 1;
    }
    x->type = mv_float;
    x->value.f = PT(n,t);
    return 0;
}

static
int sf_PU(mt_object* x, int argc, mt_object* v){
    if(argc!=2){
        mf_argc_error(argc,2,2,"sf.PU");
        return 1;
    }
    if(v[1].type!=mv_int){
        mf_type_error1("in sf.PU(n,x): n (type: %s) is not an integer.",v+1);
        return 1;
    }
    long n = v[1].value.i;
    int e=0;
    double t = mf_float(v+2,&e);
    if(e){
        mf_type_error1("in sf.PU(n,x): cannot convert x (type: %s) to float.",v+2);
        return 1;
    }
    x->type = mv_float;
    x->value.f = PU(n,t);
    return 0;
}

static
int sf_eiK(mt_object* x, int argc, mt_object* v){
    if(argc!=1){
        mf_argc_error(argc,1,1,"sf.eiK");
        return 1;
    }
    if(v[1].type==mv_complex){
        double complex w = ceiK(CMPLX(v[1].value.c.re,v[1].value.c.im));
        x->type = mv_complex;
        x->value.c.re = creal(w);
        x->value.c.im = cimag(w);
        return 0;
    }else{
        int e=0;
        double t = mf_float(v+1,&e);
        if(e){
            mf_type_error1("in sf.eiK(z): cannot convert z (type: %s) to complex.",v+1);
            return 1;
        }
        x->type = mv_float;
        x->value.f = eiK(t);
        return 0;
    }
}

static
int sf_eiE(mt_object* x, int argc, mt_object* v){
    if(argc!=1){
        mf_argc_error(argc,1,1,"sf.eiE");
        return 1;
    }
    if(v[1].type==mv_complex){
        double complex w = ceiE(CMPLX(v[1].value.c.re,v[1].value.c.im));
        x->type = mv_complex;
        x->value.c.re = creal(w);
        x->value.c.im = cimag(w);
        return 0;
    }else{
        int e=0;
        double t = mf_float(v+1,&e);
        if(e){
            mf_type_error1("in sf.eiE(z): cannot convert z (type: %s) to complex.",v+1);
            return 1;
        }
        x->type = mv_float;
        x->value.f = eiE(t);
        return 0;
    }
}

static
int sf_erf(mt_object* x, int argc, mt_object* v){
    if(argc!=1){
        mf_argc_error(argc,1,1,"sf.erf");
        return 1;
    }
    if(v[1].type==mv_complex){
        double complex w = mf_cerf(CMPLX(v[1].value.c.re,v[1].value.c.im));
        x->type = mv_complex;
        x->value.c.re = creal(w);
        x->value.c.im = cimag(w);
        return 0;
    }else{
        int e=0;
        double t = mf_float(v+1,&e);
        if(e){
            mf_type_error1("in sf.erf(z): cannot convert z (type: %s) to complex.",v+1);
            return 1;
        }
        x->type = mv_float;
        x->value.f = erf(t);
        return 0;
    }
}

static
int sf_Fmn(mt_object* x, int argc, mt_object* v){
    if(argc!=3){
        mf_argc_error(argc,3,3,"sf.F");
        return 1;
    }
    if(v[1].type!=mv_list){
        mf_type_error1("in sf.F(a,b,x): a (type: %s) is not a list.",v+1);
        return 1;
    }
    if(v[2].type!=mv_list){
        mf_type_error1("in sf.F(a,b,x): b (type: %s) is not a list.",v+2);
        return 1;
    }
    mt_list* La = (mt_list*)v[1].value.p;
    mt_list* Lb = (mt_list*)v[2].value.p;
    int e=0;
    double t = mf_float(v+3,&e);
    if(e){
        mf_type_error1("in sf.F(a,b,x): cannot convert x (type: %s) to float.",v+3);
        return 1;
    }
    double* a = mf_malloc(sizeof(double)*La->size);
    double* b = mf_malloc(sizeof(double)*Lb->size);
    long i;
    for(i=0; i<La->size; i++){
        a[i] = mf_float(La->a+i,&e);
        if(e){
            mf_type_error1("in sf.F(a,b,x): cannot convert a[k] (type: %s) to float.",La->a+i);
            mf_free(a);
            return 1;
        }
    }
    for(i=0; i<Lb->size; i++){
        b[i] = mf_float(Lb->a+i,&e);
        if(e){
            mf_type_error1("in sf.F(a,b,x): cannot convert b[k] (type: %s) to float.",Lb->a+i);
            mf_free(a);
            mf_free(b);
            return 1;
        }
    }
    x->type = mv_float;
    x->value.f = Fmn(La->size,Lb->size,a,b,t);
    mf_free(a);
    mf_free(b);
    return 0;
}

static
int sf_test(mt_object* x, int argc, mt_object* v){
    int e;
    if(argc==1){
        if(v[1].type==mv_complex){
            double complex w = czeta(CMPLX(v[1].value.c.re,v[1].value.c.im));
            x->type = mv_complex;
            x->value.c.re = creal(w);
            x->value.c.im = cimag(w);
            return 0;
        }else{
            e=0;
            double t = mf_float(v+1,&e);
            if(e){
                mf_type_error1("in test(s): cannot convert s (type: %s) to complex.",v+1);
                return 1;
            }
            x->type = mv_float;
            x->value.f = t;
            // x->value.f = gsl_sf_psi(t);
            return 0;
        }
    }else if(argc==2){
        if(v[1].type==mv_complex || v[2].type==mv_complex){
            mt_object z,a;
            if(mf_complex(&z,v+1)){
                mf_type_error1("in test(x,y): cannot convert x (type: %s) to complex.",v+1);
                return 1;
            }
            if(mf_complex(&a,v+2)){
                mf_type_error1("in test(x,y): cannot convert y (type: %s) to complex.",v+2);
                return 1;
            }
            double complex w = chzeta(
                CMPLX(z.value.c.re,z.value.c.im),
                CMPLX(a.value.c.re,a.value.c.im)
            );
            x->type = mv_complex;
            x->value.c.re = creal(w);
            x->value.c.im = cimag(w);
            return 0;
        }else{
            double rz,ra;
            e=0;
            rz = mf_float(v+1,&e);
            if(e){
                mf_type_error1("in test(x,y): cannot convert x (type: %s) to float.",v+1);
                return 1;
            }
            ra = mf_float(v+2,&e);
            if(e){
                mf_type_error1("in test(x,y): cannot convert y (type: %s) to float.",v+2);
                return 1;
            }
            x->type = mv_float;
            if(ra>0 && rz>0){
                // x->value.f=tgamma(rz)*gsl_sf_gamma_inc_P(rz,ra);
            }else if(ra>0){
                // x->value.f=tgamma(rz)-gsl_sf_gamma_inc(rz,ra);
            }else{
                x->value.f = NAN;
            }
            return 0;
        }
    }else{
        mf_argc_error(argc,1,2,"test");
        return 1;
    }
}

mt_table* mf_sf_load(){
    mt_table* sf = mf_table(NULL);
    sf->m = mf_empty_map();
    mt_map* m = sf->m;
    mf_insert_function(m,2,2,"bc",sf_bc);
    mf_insert_function(m,1,2,"gamma",sf_gamma);
    mf_insert_function(m,2,2,"Gamma",sf_Gamma);
    mf_insert_function(m,2,2,"Beta",sf_Beta);
    mf_insert_function(m,1,1,"digamma",sf_digamma);
    mf_insert_function(m,2,2,"polygamma",sf_polygamma);
    mf_insert_function(m,1,2,"zeta",sf_zeta);
    mf_insert_function(m,1,1,"B",sf_BernoulliB);
    mf_insert_function(m,2,3,"PP",sf_PP);
    mf_insert_function(m,2,3,"PH",sf_PH);
    mf_insert_function(m,2,2,"PT",sf_PT);
    mf_insert_function(m,2,2,"PU",sf_PU);
    mf_insert_function(m,1,1,"eiK",sf_eiK);
    mf_insert_function(m,1,1,"eiE",sf_eiE);
    mf_insert_function(m,1,1,"erf",sf_erf);
    mf_insert_function(m,3,3,"F",sf_Fmn);
    mf_insert_function(m,1,2,"test",sf_test);
    m->frozen = 1;
    return sf;
}

max_stars_repo_name: JohnBSmith/moss-c

max_issues_repo_name: JohnBSmith/moss-c

max_forks_repo_name: JohnBSmith/moss-c


text:
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>

int
main ()
{
  double data[] = { 1.0  , 1/2.0, 1/3.0, 1/4.0,
                    1/2.0, 1/3.0, 1/4.0, 1/5.0,
                    1/3.0, 1/4.0, 1/5.0, 1/6.0,
                    1/4.0, 1/5.0, 1/6.0, 1/7.0 };

  gsl_matrix m = gsl_matrix_view(data, 4, 4);

  gsl_vector *eval = gsl_vector_alloc (4);
  gsl_matrix *evec = gsl_matrix_alloc (4, 4);

  gsl_eigen_symmv_workspace * w =
    gsl_eigen_symmv_alloc (4);

  gsl_eigen_symmv (&m, eval, evec, w);

  gsl_eigen_symmv_free(w);

  gsl_eigen_sort (eval, evec, GSL_EIGEN_SORT_ABS_DESC);

  {
    int i;

    for (i = 0; i < 4; i++)
      {
        gsl_vector v = gsl_matrix_column(evec, i);
        printf("eigenvalue = %g\n", gsl_vector_get(eval, i));
        printf("eigenvector = \n");
        gsl_vector_fprintf(stdout, &v, "%g");
      }
  }
}

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim