text:
/* -*- linux-c -*- */
/* sigma_binsingle.c

   Copyright (C) 2002-2004 John M. Fregeau

   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 2 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

#include <stdio.h>
#include <stddef.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <getopt.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_rng.h>
#include "fewbody.h"
#include "sigma_binsingle.h"

/* print the usage */
void print_usage(FILE *stream)
{
	fprintf(stream, "USAGE:\n");
	fprintf(stream, "  sigma_binsingle [options...]\n");
	fprintf(stream, "\n");
	fprintf(stream, "OPTIONS:\n");
	fprintf(stream, "  -m --m0 <m0/MSUN>             : set mass of single star [%.6g]\n", FB_M0/FB_CONST_MSUN);
	fprintf(stream, "  -n --m10 <m10/MSUN>           : set mass of star 0 of binary [%.6g]\n", FB_M10/FB_CONST_MSUN);
	fprintf(stream, "  -o --m11 <m11/MSUN>           : set mass of star 1 of binary [%.6g]\n", FB_M11/FB_CONST_MSUN);
	fprintf(stream, "  -r --r0 <r0/RSUN>             : set radius of single star [%.6g]\n", FB_R0/FB_CONST_RSUN);
	fprintf(stream, "  -g --r10 <r10/RSUN>           : set radius of star 0 of binary [%.6g]\n", FB_R10/FB_CONST_RSUN);
	fprintf(stream, "  -i --r11 <r11/RSUN>           : set radius of star 1 of binary [%.6g]\n", FB_R11/FB_CONST_RSUN);
	fprintf(stream, "  -a --a1 <a1/AU>               : set semimajor axis of binary [%.6g]\n", FB_A1/FB_CONST_AU);
	fprintf(stream, "  -e --e1 <e1>                  : set eccentricity of binary 0 [%.6g]\n", FB_E1);
	fprintf(stream, "  -v --vinf <vinf/v_crit>       : set velocity at infinity [%.6g]\n", FB_VINF);
	fprintf(stream, "  -P --precision <dsigma/sigma> : set fractional precision on cross\n");
	fprintf(stream, "                                  section [%.6g]\n", FB_PRECISION);
	fprintf(stream, "  -t --tstop <tstop/t_dyn>      : set stopping time [%.6g]\n", FB_TSTOP);
	fprintf(stream, "  -D --dt <dt/t_dyn>            : set approximate output dt [%.6g]\n", FB_DT);
	fprintf(stream, "  -c --tcpustop <tcpustop/sec>  : set cpu stopping time [%.6g]\n", FB_TCPUSTOP);
	fprintf(stream, "  -A --absacc <absacc>          : set integrator's absolute accuracy [%.6g]\n", FB_ABSACC);
	fprintf(stream, "  -R --relacc <relacc>          : set integrator's relative accuracy [%.6g]\n", FB_RELACC);
	fprintf(stream, "  -N --ncount <ncount>          : set number of integration steps between calls\n");
	fprintf(stream, "                                  to fb_classify() [%d]\n", FB_NCOUNT);
	fprintf(stream, "  -z --tidaltol <tidaltol>      : set tidal tolerance [%.6g]\n", FB_TIDALTOL);
	fprintf(stream, "  -x --fexp <f_exp>             : set expansion factor of merger product [%.6g]\n", FB_FEXP);
	fprintf(stream, "  -k --ks                       : turn K-S regularization on or off [%d]\n", FB_KS);
	fprintf(stream, "  -s --seed                     : set random seed [%ld]\n", FB_SEED);
	fprintf(stream, "  -d --debug                    : turn on debugging\n");
	fprintf(stream, "  -V --version                  : print version info\n");
	fprintf(stream, "  -h --help                     : display this help text\n");
}

/* calculate the units used */
int calc_units(fb_obj_t *obj[2], fb_units_t *units)
{
	units->v = sqrt(FB_CONST_G*(obj[0]->m + obj[1]->m)/(obj[0]->m * obj[1]->m) * \
			(obj[1]->obj[0]->m * obj[1]->obj[1]->m / obj[1]->a));
	units->l = obj[1]->a;
	units->t = units->l / units->v;
	units->m = units->l * fb_sqr(units->v) / FB_CONST_G;
	units->E = units->m * fb_sqr(units->v);

	return(0);
}

/* the main attraction */
int main(int argc, char *argv[])
{
	int i, j, k, res, err, sid, bid, precisionmet=0;
	long l, nb;
	unsigned long int seed;
	double m0, m10, m11, r0, r10, r11, a1, e1, precision;
	double rtid, vinf, rperi, b, b0, bmin, bmax, db, bxlast, t, vc;
	fb_sigma_t sigmacurr[11], sigmaprev[11];
	char outcome[11][FB_MAX_STRING_LENGTH];
	/* char string1[FB_MAX_STRING_LENGTH], string2[FB_MAX_STRING_LENGTH]; */
	fb_hier_t hier;
	fb_input_t input;
	fb_ret_t retval;
	fb_units_t units;
	gsl_rng *rng;
	const gsl_rng_type *rng_type=gsl_rng_mt19937;
	const char *short_opts = "m:n:o:r:g:i:a:e:v:P:t:D:c:A:R:N:z:x:k:s:dVh";
	const struct option long_opts[] = {
		{"m0", required_argument, NULL, 'm'},
		{"m10", required_argument, NULL, 'n'},
		{"m11", required_argument, NULL, 'o'},
		{"r0", required_argument, NULL, 'r'},
		{"r10", required_argument, NULL, 'g'},
		{"r11", required_argument, NULL, 'i'},
		{"a1", required_argument, NULL, 'a'},
		{"e1", required_argument, NULL, 'e'},
		{"vinf", required_argument, NULL, 'v'},
		{"precision", required_argument, NULL, 'P'},
		{"tstop", required_argument, NULL, 't'},
		{"dt", required_argument, NULL, 'D'},
		{"tcpustop", required_argument, NULL, 'c'},
		{"absacc", required_argument, NULL, 'A'},
		{"relacc", required_argument, NULL, 'R'},
		{"ncount", required_argument, NULL, 'N'},
		{"tidaltol", required_argument, NULL, 'z'},
		{"fexp", required_argument, NULL, 'x'},
		{"ks", required_argument, NULL, 'k'},
		{"seed", required_argument, NULL, 's'},
		{"debug", no_argument, NULL, 'd'},
		{"version", no_argument, NULL, 'V'},
		{"help", no_argument, NULL, 'h'},
		{NULL, 0, NULL, 0}
	};

	/* set parameters to default values */
	m0 = FB_M0;
	m10 = FB_M10;
	m11 = FB_M11;
	r0 = FB_R0;
	r10 = FB_R10;
	r11 = FB_R11;
	a1 = FB_A1;
	e1 = FB_E1;
	vinf = FB_VINF;
	precision = 1.0e-2;
	input.ks = FB_KS;
	input.tstop = FB_TSTOP;
	input.Dflag = 0;
	input.dt = FB_DT;
	input.tcpustop = FB_TCPUSTOP;
	input.absacc = FB_ABSACC;
	input.relacc = FB_RELACC;
	input.ncount = FB_NCOUNT;
	input.tidaltol = FB_TIDALTOL;
	input.fexp = FB_FEXP;
	seed = FB_SEED;
	fb_debug = FB_DEBUG;

	while ((i = getopt_long(argc, argv, short_opts, long_opts, NULL)) != -1) {
		switch (i) {
		case 'm':
			m0 = atof(optarg) * FB_CONST_MSUN;
			break;
		case 'n':
			m10 = atof(optarg) * FB_CONST_MSUN;
			break;
		case 'o':
			m11 = atof(optarg) * FB_CONST_MSUN;
			break;
		case 'r':
			r0 = atof(optarg) * FB_CONST_RSUN;
			break;
		case 'g':
			r10 = atof(optarg) * FB_CONST_RSUN;
			break;
		case 'i':
			r11 = atof(optarg) * FB_CONST_RSUN;
			break;
		case 'a':
			a1 = atof(optarg) * FB_CONST_AU;
			break;
		case 'e':
			e1 = atof(optarg);
			if (e1 >= 1.0) {
				fprintf(stderr, "e0 must be less than 1\n");
				return(1);
			}
			break;
		case 'v':
			vinf = atof(optarg);
			if (vinf < 0.0) {
				fprintf(stderr, "vinf must be non-negative\n");
				return(1);
			}
			break;
		case 'P':
			precision = atof(optarg);
			if (precision <= 0.0) {
				fprintf(stderr, "precision must be > 0\n");
				return(1);
			} else if (precision > 1.0) {
				fprintf(stderr, "precision must be < 1\n");
				return(1);
			}
			break;
		case 't':
			input.tstop = atof(optarg);
			break;
		case 'D':
			input.Dflag = 1;
			input.dt = atof(optarg);
			break;
		case 'c':
			input.tcpustop = atof(optarg);
			break;
		case 'A':
			input.absacc = atof(optarg);
			break;
		case 'R':
			input.relacc = atof(optarg);
			break;
		case 'N':
			input.ncount = atoi(optarg);
			break;
		case 'z':
			input.tidaltol = atof(optarg);
			break;
		case 'x':
			input.fexp = atof(optarg);
			break;
		case 'k':
			input.ks = atoi(optarg);
			break;
		case 's':
			seed = atol(optarg);
			break;
		case 'd':
			fb_debug = 1;
			break;
		case 'V':
			fb_print_version(stdout);
			return(0);
		case 'h':
			fb_print_version(stdout);
			fprintf(stdout, "\n");
			print_usage(stdout);
			return(0);
		default:
			break;
		}
	}

	/* check to make sure there was nothing crazy on the command line */
	if (optind < argc) {
		print_usage(stdout);
		return(1);
	}

	/* initialize a few things for integrator */
	hier.nstarinit = 3;
	fb_malloc_hier(&hier);

	/* initialize GSL rng */
	gsl_rng_env_setup();
	rng = gsl_rng_alloc(rng_type);
	gsl_rng_set(rng, seed);

	/* zero everything out before summing */
	for (i=0; i<11; i++) {
		sigmacurr[i].sigma_r = 0.0;
		sigmacurr[i].dsigma2_r_minus = 0.0;
		sigmacurr[i].dsigma2_r_plus = 0.0;
		sigmacurr[i].sigma_nr = 0.0;
		sigmacurr[i].dsigma2_nr_minus = 0.0;
		sigmacurr[i].dsigma2_nr_plus = 0.0;
		sigmaprev[i].sigma_r = 0.0;
		sigmaprev[i].dsigma2_r_minus = 0.0;
		sigmaprev[i].dsigma2_r_plus = 0.0;
		sigmaprev[i].sigma_nr = 0.0;
		sigmaprev[i].dsigma2_nr_minus = 0.0;
		sigmaprev[i].dsigma2_nr_plus = 0.0;
	}
	/* loop through and calcluate cross sections */
	nb = (long) (1.0/precision);
	while (!precisionmet) {
		/* select b0 that gives r_p=2a */
		vc = sqrt(FB_CONST_G * (m0+m10+m11) / (m0*(m10+m11)) * (m10*m11/a1));
		rperi = 2.0 * a1;
		/* b0 is in code units here (vinf is in units of vc) */
		b0 = sqrt(fb_sqr(rperi) + 2.0 * FB_CONST_G * (m0+m10+m11) * rperi / fb_sqr(vinf*vc)) / a1;
		db = b0/((double) nb);
		bxlast = GSL_POSINF;

		l = 0;
		bmin = ((double) l) * db;
		bmax = bmin + db;
		while (bmax <= b0 || bmax <= 2.0*bxlast) {
			bmin = ((double) l) * db;
			bmax = bmin + db;
			/* choose b uniformly in area from bmin to bmax */
			b = sqrt(gsl_rng_uniform(rng) * (fb_sqr(bmax)-fb_sqr(bmin)) + fb_sqr(bmin));

			/* do scattering experiment */
			/* initialize a few things for integrator */
			t = 0.0;
			hier.nstar = 3;
			fb_init_hier(&hier);

			/* create binary */
			hier.hier[hier.hi[2]+0].obj[0] = &(hier.hier[hier.hi[1]+1]);
			hier.hier[hier.hi[2]+0].obj[1] = &(hier.hier[hier.hi[1]+2]);
			hier.hier[hier.hi[2]+0].t = t;

			/* give the objects some properties */
			for (j=0; j<hier.nstar; j++) {
				hier.hier[hier.hi[1]+j].ncoll = 1;
				hier.hier[hier.hi[1]+j].id[0] = j;
				snprintf(hier.hier[hier.hi[1]+j].idstring, FB_MAX_STRING_LENGTH, "%d", j);
				hier.hier[hier.hi[1]+j].n = 1;
				hier.hier[hier.hi[1]+j].obj[0] = NULL;
				hier.hier[hier.hi[1]+j].obj[1] = NULL;
				hier.hier[hier.hi[1]+j].Eint = 0.0;
				hier.hier[hier.hi[1]+j].Lint[0] = 0.0;
				hier.hier[hier.hi[1]+j].Lint[1] = 0.0;
				hier.hier[hier.hi[1]+j].Lint[2] = 0.0;
			}

			hier.hier[hier.hi[1]+0].R = r0;
			hier.hier[hier.hi[1]+1].R = r10;
			hier.hier[hier.hi[1]+2].R = r11;

			hier.hier[hier.hi[1]+0].m = m0;
			hier.hier[hier.hi[1]+1].m = m10;
			hier.hier[hier.hi[1]+2].m = m11;

			hier.hier[hier.hi[2]+0].m = m10 + m11;

			hier.hier[hier.hi[2]+0].a = a1;
			hier.hier[hier.hi[2]+0].e = e1;

			hier.obj[0] = &(hier.hier[hier.hi[1]+0]);
			hier.obj[1] = &(hier.hier[hier.hi[2]+0]);
			hier.obj[2] = NULL;

			/* get the units and normalize */
			calc_units(hier.obj, &units);
			fb_normalize(&hier, units);

			/* move hierarchies analytically in from infinity along hyperbolic orbit */
			rtid = pow(2.0*(hier.obj[0]->m + hier.obj[1]->m) / (hier.obj[1]->m * input.tidaltol), 1.0/3.0) *
				hier.obj[1]->a * (1.0 + hier.obj[1]->e);
			fb_init_scattering(hier.obj[0], hier.obj[1], vinf, b, rtid);

			/* trickle down the binary properties, then back up */
			fb_randorient(&(hier.hier[hier.hi[2]+0]), rng);
			fb_downsync(&(hier.hier[hier.hi[2]+0]), t);
			/* fb_upsync(&(hier.hier[hier.hi[2]+0]), t); */

			/* call fewbody! */
			retval = fewbody(input, &hier, &t);

			k = -1;
			res = 0;
			err = 0;

			/* analyze outcome */
			if (retval.retval == 1 &&
			    FB_MIN(fabs(retval.DeltaEfrac), fabs(retval.DeltaE)) <= 1.0e-4 &&
			    FB_MIN(fabs(retval.DeltaLfrac), fabs(retval.DeltaL)) <= 1.0e-4) {
				/* check for resonance */
				if (retval.Nosc >= 1) {
					res = 1;
				} else {
					res = 0;
				}
				/* classify outcome */
				if (hier.nstar == 3) {
					if (hier.nobj == 1) {
						/* triples physically forbidden */
						err = 1;
					} else if (hier.nobj == 2) {
						if (fb_n_hier(hier.obj[0])==2) {
							bid = 0;
							sid = 1;
						} else {
							bid = 1;
							sid = 0;
						}
						if (hier.obj[sid]->id[0] == 0) {
							k = 0;
						} else if (hier.obj[sid]->id[0] == 1) {
							k = 1;
						} else {
							k = 2;
						}
					} else { /* hier.nobj == 3 */
						k = 3;
					}
                                } else if (hier.nstar == 2) {
					if (hier.nobj == 1) {
						if ((hier.obj[0]->obj[0]->ncoll==1 && hier.obj[0]->obj[0]->id[0]==0) ||
						    (hier.obj[0]->obj[1]->ncoll==1 && hier.obj[0]->obj[1]->id[0]==0)) {
							k = 4;
						} else if ((hier.obj[0]->obj[0]->ncoll==1 && hier.obj[0]->obj[0]->id[0]==1) ||
							   (hier.obj[0]->obj[1]->ncoll==1 && hier.obj[0]->obj[1]->id[0]==1)) {
							k = 5;
						} else {
							k = 6;
						}
					} else { /* hier.nobj == 2 */
						if ((hier.obj[0]->ncoll==1 && hier.obj[0]->id[0]==0) ||
						    (hier.obj[1]->ncoll==1 && hier.obj[1]->id[0]==0)) {
							k = 7;
						} else if ((hier.obj[0]->ncoll==1 && hier.obj[0]->id[0]==1) ||
							   (hier.obj[1]->ncoll==1 && hier.obj[1]->id[0]==1)) {
							k = 8;
						} else {
							k = 9;
						}
					}
				} else if (hier.nstar == 1) {
					k = 10;
				} else {
					fprintf(stderr, "ERROR: hier.nstar!=1,2,3!\n");
					exit(1);
				}
                        } else { /* bad outcome */
				err = 1;
			}

			/* tally up cross sections */
			if (!err) {
				for (i=0; i<11; i++) {
					if (i==k) {
						if (res) {
							sigmacurr[i].sigma_r += FB_CONST_PI * (fb_sqr(bmax)-fb_sqr(bmin));
							sigmacurr[i].dsigma2_r_minus += fb_sqr(0.1587 * FB_CONST_PI * (fb_sqr(bmax)-fb_sqr(bmin)));
							sigmacurr[i].dsigma2_nr_plus += fb_sqr(0.8413/21.0 * FB_CONST_PI * (fb_sqr(bmax)-fb_sqr(bmin)));
						} else {
							sigmacurr[i].sigma_nr += FB_CONST_PI * (fb_sqr(bmax)-fb_sqr(bmin));
							sigmacurr[i].dsigma2_nr_minus += fb_sqr(0.1587 * FB_CONST_PI * (fb_sqr(bmax)-fb_sqr(bmin)));
							sigmacurr[i].dsigma2_r_plus += fb_sqr(0.8413/21.0 * FB_CONST_PI * (fb_sqr(bmax)-fb_sqr(bmin)));
						}
					} else {
						sigmacurr[i].dsigma2_r_plus += fb_sqr(0.8413/21.0 * FB_CONST_PI * (fb_sqr(bmax)-fb_sqr(bmin)));
						sigmacurr[i].dsigma2_nr_plus += fb_sqr(0.8413/21.0 * FB_CONST_PI * (fb_sqr(bmax)-fb_sqr(bmin)));
					}
				}
			} else {
				for (i=0; i<11; i++) {
					sigmacurr[i].dsigma2_r_plus += fb_sqr(0.8413/22.0 * FB_CONST_PI * (fb_sqr(bmax)-fb_sqr(bmin)));
					sigmacurr[i].dsigma2_nr_plus += fb_sqr(0.8413/22.0 * FB_CONST_PI * (fb_sqr(bmax)-fb_sqr(bmin)));
				}
			}

			/* if the outcome was not a weak flyby, then record last impact parameter */
			if (!err && (k!=0 || (k==0 && res))) {
				bxlast = bmax;
			}

			/*
			fprintf(stderr, "outcome:  retval=%d  %s  (%s)\n", retval.retval,
				fb_sprint_hier(hier, string1), fb_sprint_hier_hr(hier, string2));
			fprintf(stderr, "  DeltaL/L0=%.6g  DeltaL=%.6g\n", retval.DeltaLfrac, retval.DeltaL);
			fprintf(stderr, "  DeltaE/E0=%.6g  DeltaE=%.6g\n", retval.DeltaEfrac, retval.DeltaE);
			fprintf(stderr, "  l=%ld b0=%g bmin=%g bmax=%g bxlast=%g k=%d res=%d err=%d\n", l, b0, bmin, bmax, bxlast, k, res, err);
			*/

			l++;
		}

		/* add current and previous sigmas, and transfer current to previous */

		/* increase resolution geometrically */
		nb *= 2;

		precisionmet = 1;
	}

	snprintf(outcome[0], FB_MAX_STRING_LENGTH, "[1 2] 0 (preservation)");
	snprintf(outcome[1], FB_MAX_STRING_LENGTH, "[0 2] 1 (exchange_1)");
	snprintf(outcome[2], FB_MAX_STRING_LENGTH, "[0 1] 2 (exchange_2)");
	snprintf(outcome[3], FB_MAX_STRING_LENGTH, "0 1 2 (ionization)");
	snprintf(outcome[4], FB_MAX_STRING_LENGTH, "[1:2 0] (merger_binary_12)");
	snprintf(outcome[5], FB_MAX_STRING_LENGTH, "[0:2 1] (merger_binary_02)");
	snprintf(outcome[6], FB_MAX_STRING_LENGTH, "[0:1 2] (merger_binary_01)");
	snprintf(outcome[7], FB_MAX_STRING_LENGTH, "1:2 0 (merger_12)");
	snprintf(outcome[8], FB_MAX_STRING_LENGTH, "0:2 1 (merger_02)");
	snprintf(outcome[9], FB_MAX_STRING_LENGTH, "0:1 2 (merger_01)");
	snprintf(outcome[10], FB_MAX_STRING_LENGTH, "0:1:2 (triple_merger)");

	/* precision has now been met, so print out results */
	fprintf(stdout, "outcome  sigma_r  dsigma_r_minus  dsigma_r_plus  sigma_nr  dsigma_nr_minus  dsigma_nr_plus\n");
	for (i=0; i<11; i++) {
		fprintf(stdout, "%27s  %6.6g  %6.6g  %6.6g  %6.6g  %6.6g  %6.6g\n", outcome[i],
			sigmacurr[i].sigma_r/FB_CONST_PI,
			sqrt(sigmacurr[i].dsigma2_r_minus)/FB_CONST_PI,
			sqrt(sigmacurr[i].dsigma2_r_plus)/FB_CONST_PI,
			sigmacurr[i].sigma_nr/FB_CONST_PI,
			sqrt(sigmacurr[i].dsigma2_nr_minus)/FB_CONST_PI,
			sqrt(sigmacurr[i].dsigma2_nr_plus)/FB_CONST_PI);
	}

	/* free GSL stuff */
	gsl_rng_free(rng);

	/* free our own stuff */
	fb_free_hier(hier);

	/* done! */
	return(0);
}

max_stars_repo_name: gnodvi/cosmos

max_issues_repo_name: gnodvi/cosmos

max_forks_repo_name: gnodvi/cosmos


text:
/* vector/gsl_vector_uchar.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman, Brian Gough
 *
 * 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.
 */

#ifndef __GSL_VECTOR_UCHAR_H__
#define __GSL_VECTOR_UCHAR_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_inline.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_block_uchar.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

typedef struct
{
  size_t size;
  size_t stride;
  unsigned char *data;
  gsl_block_uchar *block;
  int owner;
}
gsl_vector_uchar;

typedef struct
{
  gsl_vector_uchar vector;
} _gsl_vector_uchar_view;

typedef _gsl_vector_uchar_view gsl_vector_uchar_view;

typedef struct
{
  gsl_vector_uchar vector;
} _gsl_vector_uchar_const_view;

typedef const _gsl_vector_uchar_const_view gsl_vector_uchar_const_view;

/* Allocation */

GSL_FUN gsl_vector_uchar *gsl_vector_uchar_alloc (const size_t n);
GSL_FUN gsl_vector_uchar *gsl_vector_uchar_calloc (const size_t n);

GSL_FUN gsl_vector_uchar *gsl_vector_uchar_alloc_from_block (gsl_block_uchar * b,
                                                     const size_t offset,
                                                     const size_t n,
                                                     const size_t stride);

GSL_FUN gsl_vector_uchar *gsl_vector_uchar_alloc_from_vector (gsl_vector_uchar * v,
                                                      const size_t offset,
                                                      const size_t n,
                                                      const size_t stride);

GSL_FUN void gsl_vector_uchar_free (gsl_vector_uchar * v);

/* Views */

GSL_FUN _gsl_vector_uchar_view
gsl_vector_uchar_view_array (unsigned char *v, size_t n);

GSL_FUN _gsl_vector_uchar_view
gsl_vector_uchar_view_array_with_stride (unsigned char *base,
                                         size_t stride,
                                         size_t n);

GSL_FUN _gsl_vector_uchar_const_view
gsl_vector_uchar_const_view_array (const unsigned char *v, size_t n);

GSL_FUN _gsl_vector_uchar_const_view
gsl_vector_uchar_const_view_array_with_stride (const unsigned char *base,
                                               size_t stride,
                                               size_t n);

GSL_FUN _gsl_vector_uchar_view
gsl_vector_uchar_subvector (gsl_vector_uchar *v,
                            size_t i,
                            size_t n);

GSL_FUN _gsl_vector_uchar_view
gsl_vector_uchar_subvector_with_stride (gsl_vector_uchar *v,
                                        size_t i,
                                        size_t stride,
                                        size_t n);

GSL_FUN _gsl_vector_uchar_const_view
gsl_vector_uchar_const_subvector (const gsl_vector_uchar *v,
                                  size_t i,
                                  size_t n);

GSL_FUN _gsl_vector_uchar_const_view
gsl_vector_uchar_const_subvector_with_stride (const gsl_vector_uchar *v,
                                              size_t i,
                                              size_t stride,
                                              size_t n);

/* Operations */

GSL_FUN void gsl_vector_uchar_set_zero (gsl_vector_uchar * v);
GSL_FUN void gsl_vector_uchar_set_all (gsl_vector_uchar * v, unsigned char x);
GSL_FUN int gsl_vector_uchar_set_basis (gsl_vector_uchar * v, size_t i);

GSL_FUN int gsl_vector_uchar_fread (FILE * stream, gsl_vector_uchar * v);
GSL_FUN int gsl_vector_uchar_fwrite (FILE * stream, const gsl_vector_uchar * v);
GSL_FUN int gsl_vector_uchar_fscanf (FILE * stream, gsl_vector_uchar * v);
GSL_FUN int gsl_vector_uchar_fprintf (FILE * stream, const gsl_vector_uchar * v,
                              const char *format);

GSL_FUN int gsl_vector_uchar_memcpy (gsl_vector_uchar * dest, const gsl_vector_uchar * src);

GSL_FUN int gsl_vector_uchar_reverse (gsl_vector_uchar * v);

GSL_FUN int gsl_vector_uchar_swap (gsl_vector_uchar * v, gsl_vector_uchar * w);
GSL_FUN int gsl_vector_uchar_swap_elements (gsl_vector_uchar * v, const size_t i, const size_t j);

GSL_FUN unsigned char gsl_vector_uchar_max (const gsl_vector_uchar * v);
GSL_FUN unsigned char gsl_vector_uchar_min (const gsl_vector_uchar * v);
GSL_FUN void gsl_vector_uchar_minmax (const gsl_vector_uchar * v, unsigned char * min_out, unsigned char * max_out);

GSL_FUN size_t gsl_vector_uchar_max_index (const gsl_vector_uchar * v);
GSL_FUN size_t gsl_vector_uchar_min_index (const gsl_vector_uchar * v);
GSL_FUN void gsl_vector_uchar_minmax_index (const gsl_vector_uchar * v, size_t * imin, size_t * imax);

GSL_FUN int gsl_vector_uchar_add (gsl_vector_uchar * a, const gsl_vector_uchar * b);
GSL_FUN int gsl_vector_uchar_sub (gsl_vector_uchar * a, const gsl_vector_uchar * b);
GSL_FUN int gsl_vector_uchar_mul (gsl_vector_uchar * a, const gsl_vector_uchar * b);
GSL_FUN int gsl_vector_uchar_div (gsl_vector_uchar * a, const gsl_vector_uchar * b);
GSL_FUN int gsl_vector_uchar_scale (gsl_vector_uchar * a, const unsigned char x);
GSL_FUN int gsl_vector_uchar_add_constant (gsl_vector_uchar * a, const double x);
GSL_FUN int gsl_vector_uchar_axpby (const unsigned char alpha, const gsl_vector_uchar * x, const unsigned char beta, gsl_vector_uchar * y);
GSL_FUN unsigned char gsl_vector_uchar_sum (const gsl_vector_uchar * a);

GSL_FUN int gsl_vector_uchar_equal (const gsl_vector_uchar * u,
                            const gsl_vector_uchar * v);

GSL_FUN int gsl_vector_uchar_isnull (const gsl_vector_uchar * v);
GSL_FUN int gsl_vector_uchar_ispos (const gsl_vector_uchar * v);
GSL_FUN int gsl_vector_uchar_isneg (const gsl_vector_uchar * v);
GSL_FUN int gsl_vector_uchar_isnonneg (const gsl_vector_uchar * v);

GSL_FUN INLINE_DECL unsigned char gsl_vector_uchar_get (const gsl_vector_uchar * v, const size_t i);
GSL_FUN INLINE_DECL void gsl_vector_uchar_set (gsl_vector_uchar * v, const size_t i, unsigned char x);
GSL_FUN INLINE_DECL unsigned char * gsl_vector_uchar_ptr (gsl_vector_uchar * v, const size_t i);
GSL_FUN INLINE_DECL const unsigned char * gsl_vector_uchar_const_ptr (const gsl_vector_uchar * v, const size_t i);

#ifdef HAVE_INLINE

INLINE_FUN
unsigned char
gsl_vector_uchar_get (const gsl_vector_uchar * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_VAL ("index out of range", GSL_EINVAL, 0);
    }
#endif
  return v->data[i * v->stride];
}

INLINE_FUN
void
gsl_vector_uchar_set (gsl_vector_uchar * v, const size_t i, unsigned char x)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_VOID ("index out of range", GSL_EINVAL);
    }
#endif
  v->data[i * v->stride] = x;
}

INLINE_FUN
unsigned char *
gsl_vector_uchar_ptr (gsl_vector_uchar * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (unsigned char *) (v->data + i * v->stride);
}

INLINE_FUN
const unsigned char *
gsl_vector_uchar_const_ptr (const gsl_vector_uchar * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (const unsigned char *) (v->data + i * v->stride);
}
#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_VECTOR_UCHAR_H__ */

ext: h

max_stars_repo_name: mgreter/astrometrylib

max_issues_repo_name: zzpwahaha/VimbaCamJILA

max_forks_repo_name: zzpwahaha/VimbaCamJILA


text:
/* ieee-utils/fp-aix.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Tim Mooney
 *
 * 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 2 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 <math.h>
#include <fptrap.h>
#include <float.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_errno.h>

int
gsl_ieee_set_mode (int precision, int rounding, int exception_mask)
{
  fptrap_t   mode = 0 ;
  fprnd_t    rnd  = 0 ;

  switch (precision)
    {

    /* I'm not positive about AIX only supporting default precision rounding,
     * but this is the best assumption until it's proven otherwise. */

    case GSL_IEEE_SINGLE_PRECISION:
      GSL_ERROR ("AIX only supports default precision rounding",
                 GSL_EUNSUP) ;
      break ;
    case GSL_IEEE_DOUBLE_PRECISION:
      GSL_ERROR ("AIX only supports default precision rounding",
                 GSL_EUNSUP) ;
      break ;
    case GSL_IEEE_EXTENDED_PRECISION:
      GSL_ERROR ("AIX only supports default precision rounding",
                 GSL_EUNSUP) ;
      break ;
    }

  switch (rounding)
    {
    case GSL_IEEE_ROUND_TO_NEAREST:
      rnd = FP_RND_RN ;
      fp_swap_rnd (rnd) ;
      break ;
    case GSL_IEEE_ROUND_DOWN:
      rnd = FP_RND_RM ;
      fp_swap_rnd (rnd) ;
      break ;
    case GSL_IEEE_ROUND_UP:
      rnd = FP_RND_RP ;
      fp_swap_rnd (rnd) ;
      break ;
    case GSL_IEEE_ROUND_TO_ZERO:
      rnd = FP_RND_RZ ;
      fp_swap_rnd (rnd) ;
      break ;
    default:
      rnd = FP_RND_RN ;
      fp_swap_rnd (rnd) ;
    }

  /* Turn on all the exceptions apart from 'inexact' */

  mode = TRP_INVALID | TRP_DIV_BY_ZERO | TRP_OVERFLOW | TRP_UNDERFLOW ;

  if (exception_mask & GSL_IEEE_MASK_INVALID)
    mode &= ~ TRP_INVALID ;

  if (exception_mask & GSL_IEEE_MASK_DENORMALIZED)
    {
      /* do nothing */
    }
  else
    {
      GSL_ERROR ("AIX does not support the denormalized operand exception. "
                 "Use 'mask-denormalized' to work around this.",
                 GSL_EUNSUP) ;
    }

  if (exception_mask & GSL_IEEE_MASK_DIVISION_BY_ZERO)
    mode &= ~ TRP_DIV_BY_ZERO ;

  if (exception_mask & GSL_IEEE_MASK_OVERFLOW)
    mode &= ~ TRP_OVERFLOW ;

  if (exception_mask & GSL_IEEE_MASK_UNDERFLOW)
    mode &=  ~ TRP_UNDERFLOW ;

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

  /* AIX appears to require two steps -- first enable floating point traps
   * in general... */
  fp_trap(FP_TRAP_SYNC);

  /* next, enable the traps we're interested in */
  fp_enable(mode);

  return GSL_SUCCESS ;

}

ext: c

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
/* Original simple matrix class taken from http://www.algarcia.org/nummeth/Cpp/Matrix.h */

#ifndef MATRIX_H
#define MATRIX_H

#include <assert.h>  // Defines the assert function.
//#define assert(x) {} //to avoid bound checking
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_errno.h>

#include <string.h>
#include <stdio.h>

#define BUFF_SIZE 4096

class Matrix {

public:

// Default Constructor. Creates a 1 by 1 matrix; sets value to zero.
Matrix () {
  ma = gsl_matrix_alloc( 1, 1);// Allocate memory
  ve.owner = 0;
  ve.block = ma->block;
  set(0.0);                // Set value of data_[0] to 0.0

  header = NULL;
}

// Regular Constructor. Creates an nR by nC matrix; sets values to zero.
// If number of columns is not specified, it is set to 1.
Matrix(int nR, int nC = 1, char *head=NULL) {
  assert(nR > 0 && nC > 0);    // Check that nC and nR both > 0.
  ma = gsl_matrix_alloc( nR, nC);// Allocate memory
  ve.owner = 0;
  ve.block = ma->block;
  set(0.0);                    // Set values of data_[] to 0.0

  if (head != NULL)
	header = strdup(head);
  else
	header = NULL;
}

//copy constructor
Matrix(Matrix &mat)  {
printf("Matrix::Matrix 2 ...\n");
  ma = gsl_matrix_alloc( mat.nRow(), mat.nCol());// Allocate memory
  ve.owner = 0;
  ve.block = ma->block;
  copy(mat);

  if (mat.Header() != NULL)
	header = strdup(mat.Header());
  else
	header = NULL;

}

//Identity (preferably square) matrix
void Iden() {
  //assert(nRow() == nCol());    // Check that square matrix.
  set(0.0);                    // Set values of data_[] to 0.0
  for (int i=0; i<nRow(); i++)
	ma->data[ (ma->tda)*i + i] = 1.0;
}

// Destructor. Called when a Matrix object goes out of scope or deleted.
~Matrix() {
	//printf("Matrix::~Matrix\n");
	if (ma != NULL)
		  gsl_matrix_free(ma);   // Release allocated memory
	if (header != NULL) {
		free(header);
	}
}

// Assignment operator function.
// Overloads the equal sign operator to work with
// Matrix objects.
Matrix& operator=(const Matrix& mat) {
  if( this == &mat ) return *this;  // If two sides equal, do nothing.
  this->copy(mat);                  // Copy right hand side to l.h.s.
  return *this;
}

// += operator function.
// Overloads the += sign operator to work with
// Matrix objects.
Matrix& operator+=(const Matrix& mat) {
  gsl_matrix_add( this->ma, mat.ma);
  return *this;
}

// -= operator function.
// Overloads the -= sign operator to work with
// Matrix objects.
Matrix& operator-=(const Matrix& mat) {
  gsl_matrix_sub( this->ma, mat.ma);
  return *this;
}

// *= operator function.
// Overloads the *= sign operator to work with
// Matrix objects.
Matrix& operator*=(const double a) {
  gsl_matrix_scale( this->ma, a);
  return *this;
}

// Set function. Sets all elements of a matrix to a given value.
void set(double value) {
   gsl_matrix_set_all( ma, value);
}

/*// Set function. Sets all elements of a matrix to a given value.
void eval(double value, double (*fct)(double x)) {
   to be done
}*/

//info
// Simple "get" functions. Return number of rows or columns.
int nRow() const { return ma->size1; }
int nCol() const { return ma->size2; }

//min and max
double Max() { gsl_matrix_max(ma); }
double Min() { gsl_matrix_min(ma); }

// Parenthesis operator function.
// Allows access to values of Matrix via (i,j) pair.
// Example: a(1,1) = 2*b(2,3);
// If column is unspecified, take as 1.
double& operator() (int i, int j = 0) {
  assert(i >= 0 && i < ma->size1);          // Bounds checking for rows
  assert(j >= 0 && j < ma->size2);          // Bounds checking for columns
  return ma->data[ (ma->tda)*i + j];  // Access appropriate value
}

// Parenthesis operator function (const version).
const double& operator() (int i, int j = 0) const{
  assert(i >= 0 && i < ma->size1);          // Bounds checking for rows
  assert(j >= 0 && j < ma->size2);          // Bounds checking for columns
  return ma->data[(ma->tda)*i + j];  // Access appropriate value
}

// Allows access to values of Matrix via ele(i,j) pair.
// If column is unspecified, take as 1.
double ele(int i, int j = 0) {
  assert(i >= 0 && i < ma->size1);          // Bounds checking for rows
  assert(j >= 0 && j < ma->size2);          // Bounds checking for columns
  return ma->data[ (ma->tda)*i + j];  // Access appropriate value
}

//i/o
//Print function
void print( FILE *F, char *name=NULL, int fw=11, int pres=4, char *sep=" ", int rownums=0, int colnums=0)
{
	if (name == NULL)
		if (header == NULL)
			fprintf( F, "\n");
		else
			fprintf( F, "%s\n", header);
	else
		fprintf( F, "%s", name);
	for (int i=0; i<ma->size1; i++)
	{
		if ((i == 0) && (colnums)) {
			fprintf( F, "%*s%s", fw, "", sep);
			for (int j=0; j<ma->size2; j++)
				fprintf( F, "%*d%s", fw, j, sep);
			fprintf(F, "\n");
		}
		for (int j=0; j<ma->size2; j++) {
			if ((j == 0) && (rownums))
				fprintf( F, "%3d%s", i, sep);
			fprintf( F, "%*.*g%s", fw, pres, ele(i,j), sep);
		}
		fprintf( F, "\n");
	}
}
void print(char *name=NULL) { print( stdout, name, 11, 4, " ", 0, 0); }
void printnums(char *name=NULL) { print( stdout, name, 11, 4, " ", 1, 1); }

void fileprint(char *fnam) {
	FILE *F;

	if ((F = fopen( fnam, "w+")) == NULL)
	{
		printf( "Could not open file %s for writing\n", fnam);

		assert(F == NULL);
	}
	else
	{
		print( F, NULL, 11, 4, " ");

		fclose(F);
	}
}

void printrow(int i) {
	assert(i < nRow());

	for (int j=0; j<ma->size2; j++) {
		printf( "%*.*g%s", 11, 4, ele(i,j), " ");
	}
	printf( "\n");
}

int filescan(char *fnam, int file_header=0) {
	FILE *F;

	if ((F = fopen( fnam, "r")) == NULL)
	{
		printf( "File %s not found\n", fnam);

		return 0;
	}
	else
	{
		if (file_header == 1) {
			header = (char *) malloc((size_t) BUFF_SIZE);

			header = fgets( header, (size_t) BUFF_SIZE, F);
			header[strlen(header)-1] = '\0'; //remove the \n
		}

		gsl_matrix_fscanf( F, ma);

		fclose(F);
		return 1;
	}
}

/*********************** Views *********************/
//Get the gsl matrix
gsl_matrix *Ma() { return ma; }

//Column gsl Vector
gsl_vector *AsColVec(int Col)   {
	assert(Col >= 0 && Col < ma->size2);          // Bounds checking for columns
	ve.size = ma->size1;
	ve.stride = ma->tda;
	ve.data = ma->data + Col;

	return &ve;
}
//First column is the defalut
gsl_vector *Ve()   {
	ve.size = ma->size1;
	ve.stride = ma->tda;
	ve.data = ma->data;

	return &ve;
}

//Row Vector
gsl_vector *AsRowVec(int Row, int ChopOff=0)   {
	assert(Row >= 0 && Row < ma->size1);          // Bounds checking for rows
	ve.size = ma->size2-ChopOff;
	ve.stride = 1;
	ve.data = ma->data + (ma->tda)*Row;

	return &ve;
}

// Copy function.
// Copies values from one Matrix object to another.
void Copy(const Matrix *mat) {

	assert((nRow() == mat->nRow()) && (nCol() == mat->nCol()));
    gsl_matrix_memcpy( ma, mat->ma);
}

// Copy function.
// Copies values from one Matrix object to another.
void copy(const Matrix& mat) {
    gsl_matrix_memcpy( ma, mat.ma);
}

const char *Header() { return header; }

//*********************************************************************
protected:

// Matrix data.
gsl_matrix *ma;
gsl_vector ve;

char *header;

}; // Class Matrix

class SubMatrix : public Matrix {

private:
	gsl_matrix_view view;
	Matrix *Parent;

public:
	SubMatrix(Matrix &mat, size_t k1, size_t k2, size_t n1, size_t n2) {
		view = gsl_matrix_submatrix ( mat.Ma(), k1, k2, n1, n2);
		ma = &view.matrix;
		Parent = &mat;
		if (mat.Header() != NULL)
			header = strdup(mat.Header());
		else
			header = NULL;
	}

	SubMatrix( Matrix& mat, size_t n1, size_t n2) {
		view = gsl_matrix_submatrix ( mat.Ma(), (size_t) 0, (size_t) 0, n1, n2);
		ma = &view.matrix;
		Parent = &mat;
		if (mat.Header() != NULL)
			header = strdup(mat.Header());
		else
			header = NULL;
	}

	SubMatrix() {
	//delay the opening of this submatrix
		ma = NULL;
		Parent = NULL;
		header = NULL;
	}

	~SubMatrix() {
		ma = NULL; //we avoid the base class call to free with this
		if (header != NULL) {
			free(header);
			header=NULL;
		}

	}

	void Set(Matrix *mat, size_t k1, size_t k2, size_t n1, size_t n2) {
		view = gsl_matrix_submatrix ( mat->Ma(), k1, k2, n1, n2);
		ma = &view.matrix;
		Parent = mat;
		if (mat->Header() == NULL)
			header = NULL;
		else
			header = strdup(mat->Header());
	}

	const char *SetHeader(char *head) {
		if (header != NULL)
			free(header);
		header == NULL;
		if (head != NULL)
			header = strdup(head);
		return header;
	}

// Assignment operator function.
// Overloads the equal sign operator to work with
// Matrix objects.
SubMatrix& operator=(const SubMatrix& mat) {
  if( this == &mat ) return *this;  // If two sides equal, do nothing.
  this->copy(mat);                  // Copy right hand side to l.h.s.
  return *this;
}

// += operator function.
// Overloads the += sign operator to work with
// Matrix objects.
SubMatrix& operator+=(const SubMatrix& mat) {
  gsl_matrix_add( this->ma, mat.ma);
  return *this;
}

// -= operator function.
// Overloads the -= sign operator to work with
// Matrix objects.
SubMatrix& operator-=(const SubMatrix& mat) {
  gsl_matrix_sub( this->ma, mat.ma);
  return *this;
}

// *= operator function.
// Overloads the *= sign operator to work with
// Matrix objects.
SubMatrix& operator*=(const double a) {
  gsl_matrix_scale( this->ma, a);
  return *this;
}

	void Set( Matrix *mat, size_t n1, size_t n2) {
		view = gsl_matrix_submatrix ( mat->Ma(), (size_t) 0, (size_t) 0, n1, n2);
		ma = &view.matrix;
		Parent = mat;
	}

	void ReSize(size_t k1, size_t k2, size_t n1, size_t n2) {
		view = gsl_matrix_submatrix( Parent->Ma(), k1, k2, n1, n2);
		ma = &view.matrix;
	}

	void ReSize(size_t n1, size_t n2) {
		view = gsl_matrix_submatrix ( Parent->Ma(), (size_t) 0, (size_t) 0, n1, n2);
		ma = &view.matrix;
	}

};

/*Auxiliary functions*/
//C = A^trA B^trB
int MatMul(Matrix &C, int trA, Matrix &A, int trB, Matrix &B);
int MatMul(Matrix &C, Matrix &A, Matrix &B);
double CuadForm(Matrix &A, Matrix &b);
int MatMulSym(Matrix &C, Matrix &A, Matrix &B);
//inverse of real (symetric) positive defined matrix
//as a by product we may obtain the square root of the determinat of A
int InvRPD(Matrix &R, Matrix &A, double *LogDetSqrt=NULL, Matrix *ExternChol=NULL);
//I have seen problems with this when dest and src are rows or cls of the same matrix
int VecCopy(gsl_vector * dest, const gsl_vector * src);

#endif

max_stars_repo_name: SFlantua/Workshops

max_issues_repo_name: SFlantua/Workshops

max_forks_repo_name: SFlantua/Workshops


text:
/* specfunc/bessel_zero.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * 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 2 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.
 */

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_airy.h>
#include <gsl/gsl_sf_pow_int.h>
#include <gsl/gsl_sf_bessel.h>

#include "error.h"

#include "bessel_olver.h"

/* For Chebyshev expansions of the roots as functions of nu,
 * see [G. Nemeth, Mathematical Approximation of Special Functions].
 * This gives the fits for all nu and s <= 10.
 * I made the fits for other values of s myself [GJ].
 */

/* Chebyshev expansion: j_{nu,1} = c_k T_k*(nu/2), nu <= 2 */
static const double coef_jnu1_a[] = {
  3.801775243633476,
  1.360704737511120,
 -0.030707710261106,
  0.004526823746202,
 -0.000808682832134,
  0.000159218792489,
 -0.000033225189761,
  0.000007205599763,
 -0.000001606110397,
  0.000000365439424,
 -0.000000084498039,
  0.000000019793815,
 -0.000000004687054,
  0.000000001120052,
 -0.000000000269767,
  0.000000000065420,
 -0.000000000015961,
  0.000000000003914,
 -0.000000000000965,
  0.000000000000239,
 -0.000000000000059,
  0.000000000000015,
 -0.000000000000004,
  0.000000000000001
};

/* Chebyshev expansion: j_{nu,1} = nu c_k T_k*((2/nu)^(2/3)), nu >= 2 */
static const double coef_jnu1_b[] = {
  1.735063412537096,
  0.784478100951978,
  0.048881473180370,
 -0.000578279783021,
 -0.000038984957864,
  0.000005758297879,
 -0.000000327583229,
 -0.000000003853878,
  0.000000002284653,
 -0.000000000153079,
 -0.000000000000895,
  0.000000000000283,
  0.000000000000043,
  0.000000000000010,
 -0.000000000000003
};

/* Chebyshev expansion: j_{nu,2} = c_k T_k*(nu/2), nu <= 2 */
static const double coef_jnu2_a[] = {
  6.992370244046161,
  1.446379282056534,
 -0.023458616207293,
  0.002172149448700,
 -0.000246262775620,
  0.000030990180959,
 -0.000004154183047,
  0.000000580766328,
 -0.000000083648175,
  0.000000012317355,
 -0.000000001844887,
  0.000000000280076,
 -0.000000000042986,
  0.000000000006658,
 -0.000000000001039,
  0.000000000000163,
 -0.000000000000026,
  0.000000000000004,
 -0.000000000000001
};

/* Chebyshev expansion: j_{nu,2} = nu c_k T_k*((2/nu)^(2/3)), nu >= 2 */
static const double coef_jnu2_b[] = {
  2.465611864263400,
  1.607952988471069,
  0.138758034431497,
 -0.003687791182054,
 -0.000051276007868,
  0.000045113570749,
 -0.000007579172152,
  0.000000736469208,
 -0.000000011118527,
 -0.000000011919884,
  0.000000002696788,
 -0.000000000314488,
  0.000000000008124,
  0.000000000005211,
 -0.000000000001292,
  0.000000000000158,
 -0.000000000000004,
 -0.000000000000003,
  0.000000000000001
};

/* Chebyshev expansion: j_{nu,3} = c_k T_k*(nu/3), nu <= 3 */
static const double coef_jnu3_a[] = {
  10.869647065239236,
   2.177524286141710,
  -0.034822817125293,
   0.003167249102413,
  -0.000353960349344,
   0.000044039086085,
  -0.000005851380981,
   0.000000812575483,
  -0.000000116463617,
   0.000000017091246,
  -0.000000002554376,
   0.000000000387335,
  -0.000000000059428,
   0.000000000009207,
  -0.000000000001438,
   0.000000000000226,
  -0.000000000000036,
   0.000000000000006,
  -0.000000000000001
};

/* Chebyshev expansion: j_{nu,3} = nu c_k T_k*((3/nu)^(2/3)), nu >= 3 */
static const double coef_jnu3_b[] = {
  2.522816775173244,
  1.673199424973720,
  0.146431617506314,
 -0.004049001763912,
 -0.000039517767244,
  0.000048781729288,
 -0.000008729705695,
  0.000000928737310,
 -0.000000028388244,
 -0.000000012927432,
  0.000000003441008,
 -0.000000000471695,
  0.000000000025590,
  0.000000000005502,
 -0.000000000001881,
  0.000000000000295,
 -0.000000000000020,
 -0.000000000000003,
  0.000000000000001
};

/* Chebyshev expansion: j_{nu,4} = c_k T_k*(nu/4), nu <= 4 */
static const double coef_jnu4_a[] = {
  14.750310252773009,
   2.908010932941708,
  -0.046093293420315,
   0.004147172321412,
  -0.000459092310473,
   0.000056646951906,
  -0.000007472351546,
   0.000001031210065,
  -0.000000147008137,
   0.000000021475218,
  -0.000000003197208,
   0.000000000483249,
  -0.000000000073946,
   0.000000000011431,
  -0.000000000001782,
   0.000000000000280,
  -0.000000000000044,
   0.000000000000007,
  -0.000000000000001
};

/* Chebyshev expansion: j_{nu,4} = nu c_k T_k*((4/nu)^(2/3)), nu >= 4 */
static const double coef_jnu4_b[] = {
  2.551681323117914,
  1.706177978336572,
  0.150357658406131,
 -0.004234001378590,
 -0.000033854229898,
  0.000050763551485,
 -0.000009337464057,
  0.000001029717834,
 -0.000000037474196,
 -0.000000013450153,
  0.000000003836180,
 -0.000000000557404,
  0.000000000035748,
  0.000000000005487,
 -0.000000000002187,
  0.000000000000374,
 -0.000000000000031,
 -0.000000000000003,
  0.000000000000001
};

/* Chebyshev expansion: j_{nu,5} = c_k T_k*(nu/5), nu <= 5 */
static const double coef_jnu5_a[] = {
  18.632261081028211,
   3.638249012596966,
  -0.057329705998828,
   0.005121709126820,
  -0.000563325259487,
   0.000069100826174,
  -0.000009066603030,
   0.000001245181383,
  -0.000000176737282,
   0.000000025716695,
  -0.000000003815184,
   0.000000000574839,
  -0.000000000087715,
   0.000000000013526,
  -0.000000000002104,
   0.000000000000330,
  -0.000000000000052,
   0.000000000000008,
  -0.000000000000001
};

/* Chebyshev expansion: j_{nu,5} = nu c_k T_k*((5/nu)^(2/3)), nu >= 5 */
/* FIXME: There is something wrong with this fit, in about the
 * 9th or 10th decimal place.
 */
static const double coef_jnu5_b[] = {
  2.569079487591442,
  1.726073360882134,
  0.152740776809531,
 -0.004346449660148,
 -0.000030512461856,
  0.000052000821080,
 -0.000009713343981,
  0.000001091997863,
 -0.000000043061707,
 -0.000000013779413,
  0.000000004082870,
 -0.000000000611259,
  0.000000000042242,
  0.000000000005448,
 -0.000000000002377,
  0.000000000000424,
 -0.000000000000038,
 -0.000000000000002,
  0.000000000000002
};

/* Chebyshev expansion: j_{nu,6} = c_k T_k*(nu/6), nu <= 6 */
static const double coef_jnu6_a[] = {
  22.514836143374042,
   4.368367257557198,
  -0.068550155285562,
   0.006093776505822,
  -0.000667152784957,
   0.000081486022398,
  -0.000010649011647,
   0.000001457089679,
  -0.000000206105082,
   0.000000029894724,
  -0.000000004422012,
   0.000000000664471,
  -0.000000000101140,
   0.000000000015561,
  -0.000000000002416,
   0.000000000000378,
  -0.000000000000060,
   0.000000000000009,
  -0.000000000000002
};

/* Chebyshev expansion: j_{nu,6} = nu c_k T_k*((6/nu)^(2/3)), nu >= 6 */
static const double coef_jnu6_b[] = {
  2.580710285494837,
  1.739380728566154,
  0.154340696401691,
 -0.004422028860168,
 -0.000028305272624,
  0.000052845975269,
 -0.000009968794373,
  0.000001134252926,
 -0.000000046841241,
 -0.000000014007555,
  0.000000004251816,
 -0.000000000648213,
  0.000000000046728,
  0.000000000005414,
 -0.000000000002508,
  0.000000000000459,
 -0.000000000000043,
 -0.000000000000002,
  0.000000000000002
};

/* Chebyshev expansion: j_{nu,7} = c_k T_k*(nu/7), nu <= 7 */
static const double coef_jnu7_a[] = {
  26.397760539730869,
   5.098418721711790,
  -0.079761896398948,
   0.007064521280487,
  -0.000770766522482,
   0.000093835449636,
  -0.000012225308542,
   0.000001667939800,
  -0.000000235288157,
   0.000000034040347,
  -0.000000005023142,
   0.000000000753101,
  -0.000000000114389,
   0.000000000017564,
  -0.000000000002722,
   0.000000000000425,
  -0.000000000000067,
   0.000000000000011,
  -0.000000000000002
};

/* Chebyshev expansion: j_{nu,7} = nu c_k T_k*((7/nu)^(2/3)), nu >= 7 */
static const double coef_jnu7_b[] = {
  2.589033335856773,
  1.748907007612678,
  0.155488900387653,
 -0.004476317805688,
 -0.000026737952924,
  0.000053459680946,
 -0.000010153699240,
  0.000001164804272,
 -0.000000049566917,
 -0.000000014175403,
  0.000000004374840,
 -0.000000000675135,
  0.000000000050004,
  0.000000000005387,
 -0.000000000002603,
  0.000000000000485,
 -0.000000000000047,
 -0.000000000000002,
  0.000000000000002
};

/* Chebyshev expansion: j_{nu,8} = c_k T_k*(nu/8), nu <= 8 */
static const double coef_jnu8_a[] = {
  30.280900001606662,
   5.828429205461221,
  -0.090968381181069,
   0.008034479731033,
  -0.000874254899080,
   0.000106164151611,
  -0.000013798098749,
   0.000001878187386,
  -0.000000264366627,
   0.000000038167685,
  -0.000000005621060,
   0.000000000841165,
  -0.000000000127538,
   0.000000000019550,
  -0.000000000003025,
   0.000000000000472,
  -0.000000000000074,
   0.000000000000012,
  -0.000000000000002
};

/* Chebyshev expansion: j_{nu,8} = nu c_k T_k*((8/nu)^(2/3)), nu >= 8 */
static const double coef_jnu8_b[] = {
  2.595283877150078,
  1.756063044986928,
  0.156352972371030,
 -0.004517201896761,
 -0.000025567187878,
  0.000053925472558,
 -0.000010293734486,
  0.000001187923085,
 -0.000000051625122,
 -0.000000014304212,
  0.000000004468450,
 -0.000000000695620,
  0.000000000052500,
  0.000000000005367,
 -0.000000000002676,
  0.000000000000505,
 -0.000000000000050,
 -0.000000000000002,
  0.000000000000002
};

/* Chebyshev expansion: j_{nu,9} = c_k T_k*(nu/9), nu <= 9 */
static const double coef_jnu9_a[] = {
  34.164181213238386,
   6.558412747925228,
  -0.102171455365016,
   0.009003934361201,
  -0.000977663914535,
   0.000118479876579,
  -0.000015368714220,
   0.000002088064285,
  -0.000000293381154,
   0.000000042283900,
  -0.000000006217033,
   0.000000000928887,
  -0.000000000140627,
   0.000000000021526,
  -0.000000000003326,
   0.000000000000518,
  -0.000000000000081,
   0.000000000000013,
  -0.000000000000002
};

/* Chebyshev expansion: j_{nu,9} = nu c_k T_k*((9/nu)^(2/3)), nu >= 9 */
static const double coef_jnu9_b[] = {
  2.600150240905079,
  1.761635491694032,
  0.157026743724010,
 -0.004549100368716,
 -0.000024659248617,
  0.000054291035068,
 -0.000010403464334,
  0.000001206027524,
 -0.000000053234089,
 -0.000000014406241,
  0.000000004542078,
 -0.000000000711728,
  0.000000000054464,
  0.000000000005350,
 -0.000000000002733,
  0.000000000000521,
 -0.000000000000052,
 -0.000000000000002,
  0.000000000000002
};

/* Chebyshev expansion: j_{nu,10} = c_k T_k*(nu/10), nu <= 10 */
static const double coef_jnu10_a[] = {
  38.047560766184647,
   7.288377637926008,
  -0.113372193277897,
   0.009973047509098,
  -0.001081019701335,
   0.000130786983847,
  -0.000016937898538,
   0.000002297699179,
  -0.000000322354218,
   0.000000046392941,
  -0.000000006811759,
   0.000000001016395,
  -0.000000000153677,
   0.000000000023486,
  -0.000000000003616,
   0.000000000000561,
  -0.000000000000095,
   0.000000000000027,
  -0.000000000000013,
   0.000000000000005
};

/* Chebyshev expansion: j_{nu,10} = nu c_k T_k*((10/nu)^(2/3)), nu >= 10 */
static const double coef_jnu10_b[] = {
  2.604046346867949,
  1.766097596481182,
  0.157566834446511,
 -0.004574682244089,
 -0.000023934500688,
  0.000054585558231,
 -0.000010491765415,
  0.000001220589364,
 -0.000000054526331,
 -0.000000014489078,
  0.000000004601510,
 -0.000000000724727,
  0.000000000056049,
  0.000000000005337,
 -0.000000000002779,
  0.000000000000533,
 -0.000000000000054,
 -0.000000000000002,
  0.000000000000002
};

/* Chebyshev expansion: j_{nu,11} = c_k T_k*(nu/22), nu <= 22 */
static const double coef_jnu11_a[] = {
  49.5054081076848637,
  15.33692279367165101,
 -0.33677234163517130,
  0.04623235772920729,
 -0.00781084960665093,
  0.00147217395434708,
 -0.00029695043846867,
  0.00006273356860235,
 -0.00001370575125628,
  3.07171282012e-6,
 -7.0235041249e-7,
  1.6320559339e-7,
 -3.843117306e-8,
  9.15083800e-9,
 -2.19957642e-9,
  5.3301703e-10,
 -1.3007541e-10,
  3.193827e-11,
 -7.88605e-12,
  1.95918e-12,
 -4.9020e-13,
  1.2207e-13,
 -2.820e-14,
  5.25e-15,
 -1.88e-15,
  2.80e-15,
 -2.45e-15
};

/* Chebyshev expansion: j_{nu,12} = c_k T_k*(nu/24), nu <= 24 */
static const double coef_jnu12_a[] = {
  54.0787833216641519,
  16.7336367772863598,
 -0.36718411124537953,
  0.05035523375053820,
 -0.00849884978867533,
  0.00160027692813434,
 -0.00032248114889921,
  0.00006806354127199,
 -0.00001485665901339,
  3.32668783672e-6,
 -7.5998952729e-7,
  1.7644939709e-7,
 -4.151538210e-8,
  9.87722772e-9,
 -2.37230133e-9,
  5.7442875e-10,
 -1.4007767e-10,
  3.437166e-11,
 -8.48215e-12,
  2.10554e-12,
 -5.2623e-13,
  1.3189e-13,
 -3.175e-14,
  5.73e-15,
  5.6e-16,
 -8.7e-16,
 -6.5e-16
};

/* Chebyshev expansion: j_{nu,13} = c_k T_k*(nu/26), nu <= 26 */
static const double coef_jnu13_a[] = {
  58.6521941921708890,
  18.1303398137970284,
 -0.39759381380126650,
  0.05447765240465494,
 -0.00918674227679980,
  0.00172835361420579,
 -0.00034800528297612,
  0.00007339183835188,
 -0.00001600713368099,
  3.58154960392e-6,
 -8.1759873497e-7,
  1.8968523220e-7,
 -4.459745253e-8,
  1.060304419e-8,
 -2.54487624e-9,
  6.1580214e-10,
 -1.5006751e-10,
  3.679707e-11,
 -9.07159e-12,
  2.24713e-12,
 -5.5943e-13,
  1.4069e-13,
 -3.679e-14,
  1.119e-14,
 -4.99e-15,
  3.43e-15,
 -2.85e-15,
  2.3e-15,
 -1.7e-15,
  8.7e-16
};

/* Chebyshev expansion: j_{nu,14} = c_k T_k*(nu/28), nu <= 28 */
static const double coef_jnu14_a[] = {
  63.2256329577315566,
  19.5270342832914901,
 -0.42800190567884337,
  0.05859971627729398,
 -0.00987455163523582,
  0.00185641011402081,
 -0.00037352439419968,
  0.00007871886257265,
 -0.00001715728110045,
  3.83632624437e-6,
 -8.7518558668e-7,
  2.0291515353e-7,
 -4.767795233e-8,
  1.132844415e-8,
 -2.71734219e-9,
  6.5714886e-10,
 -1.6005342e-10,
  3.922557e-11,
 -9.66637e-12,
  2.39379e-12,
 -5.9541e-13,
  1.4868e-13,
 -3.726e-14,
  9.37e-15,
 -2.36e-15,
  6.0e-16
};

/* Chebyshev expansion: j_{nu,15} = c_k T_k*(nu/30), nu <= 30 */
static const double coef_jnu15_a[] = {
  67.7990939565631635,
  20.9237219226859859,
 -0.45840871823085836,
  0.06272149946755639,
 -0.01056229551143042,
  0.00198445078693100,
 -0.00039903958650729,
  0.00008404489865469,
 -0.00001830717574922,
  4.09103745566e-6,
 -9.3275533309e-7,
  2.1614056403e-7,
 -5.075725222e-8,
  1.205352081e-8,
 -2.88971837e-9,
  6.9846848e-10,
 -1.7002946e-10,
  4.164941e-11,
 -1.025859e-11,
  2.53921e-12,
 -6.3128e-13,
  1.5757e-13,
 -3.947e-14,
  9.92e-15,
 -2.50e-15,
  6.3e-16
};

/* Chebyshev expansion: j_{nu,16} = c_k T_k*(nu/32), nu <= 32 */
static const double coef_jnu16_a[] = {
  72.3725729616724770,
  22.32040402918608585,
 -0.48881449782358690,
  0.06684305681828766,
 -0.01124998690363398,
  0.00211247882775445,
 -0.00042455166484632,
  0.00008937015316346,
 -0.00001945687139551,
  4.34569739281e-6,
 -9.9031173548e-7,
  2.2936247195e-7,
 -5.383562595e-8,
  1.277835103e-8,
 -3.06202860e-9,
  7.3977037e-10,
 -1.8000071e-10,
  4.407196e-11,
 -1.085046e-11,
  2.68453e-12,
 -6.6712e-13,
  1.6644e-13,
 -4.168e-14,
  1.047e-14,
 -2.64e-15,
  6.7e-16
};

/* Chebyshev expansion: j_{nu,17} = c_k T_k*(nu/34), nu <= 34 */
static const double coef_jnu17_a[] = {
  76.9460667535209549,
  23.71708159112252670,
 -0.51921943142405352,
  0.07096442978067622,
 -0.01193763559341369,
  0.00224049662974902,
 -0.00045006122941781,
  0.00009469477941684,
 -0.00002060640777107,
  4.60031647195e-6,
 -1.04785755046e-6,
  2.4258161247e-7,
 -5.691327087e-8,
  1.350298805e-8,
 -3.23428733e-9,
  7.8105847e-10,
 -1.8996825e-10,
  4.649350e-11,
 -1.144205e-11,
  2.82979e-12,
 -7.0294e-13,
  1.7531e-13,
 -4.388e-14,
  1.102e-14,
 -2.78e-15,
  7.0e-16
};

/* Chebyshev expansion: j_{nu,18} = c_k T_k*(nu/36), nu <= 36 */
static const double coef_jnu18_a[] = {
  81.5195728368096659,
  25.11375537470259305,
 -0.54962366347317668,
  0.07508565026117689,
 -0.01262524908033818,
  0.00236850602019778,
 -0.00047556873651929,
  0.00010001889347161,
 -0.00002175581482429,
  4.85490251239e-6,
 -1.10539483940e-6,
  2.5579853343e-7,
 -5.999033352e-8,
  1.422747129e-8,
 -3.40650521e-9,
  8.2233565e-10,
 -1.9993286e-10,
  4.891426e-11,
 -1.203343e-11,
  2.97498e-12,
 -7.3875e-13,
  1.8418e-13,
 -4.608e-14,
  1.157e-14,
 -2.91e-15,
  7.4e-16
};

/* Chebyshev expansion: j_{nu,19} = c_k T_k*(nu/38), nu <= 38 */
static const double coef_jnu19_a[] = {
  86.0930892477047512,
  26.51042598308271729,
 -0.58002730731948358,
  0.07920674321589394,
 -0.01331283320930301,
  0.00249650841778073,
 -0.00050107453900793,
  0.00010534258471335,
 -0.00002290511552874,
  5.10946148897e-6,
 -1.16292517157e-6,
  2.6901365037e-7,
 -6.306692473e-8,
  1.495183048e-8,
 -3.57869025e-9,
  8.6360410e-10,
 -2.0989514e-10,
  5.133439e-11,
 -1.262465e-11,
  3.12013e-12,
 -7.7455e-13,
  1.9304e-13,
 -4.829e-14,
  1.212e-14,
 -3.05e-15,
  7.7e-16
};

/* Chebyshev expansion: j_{nu,20} = c_k T_k*(nu/40), nu <= 40 */
static const double coef_jnu20_a[] = {
  90.6666144195163770,
  27.9070938975436823,
 -0.61043045315390591,
  0.08332772844325554,
 -0.01400039260208282,
  0.00262450494035660,
 -0.00052657891389470,
  0.00011066592304919,
 -0.00002405432778364,
  5.36399803946e-6,
 -1.22044976064e-6,
  2.8222728362e-7,
 -6.614312964e-8,
  1.567608839e-8,
 -3.75084856e-9,
  9.0486546e-10,
 -2.1985553e-10,
  5.375401e-11,
 -1.321572e-11,
  3.26524e-12,
 -8.1033e-13,
  2.0190e-13,
 -5.049e-14,
  1.267e-14,
 -3.19e-15,
  8.0e-16,
 -2.0e-16
};

static const double * coef_jnu_a[] = {
  0,
  coef_jnu1_a,
  coef_jnu2_a,
  coef_jnu3_a,
  coef_jnu4_a,
  coef_jnu5_a,
  coef_jnu6_a,
  coef_jnu7_a,
  coef_jnu8_a,
  coef_jnu9_a,
  coef_jnu10_a,
  coef_jnu11_a,
  coef_jnu12_a,
  coef_jnu13_a,
  coef_jnu14_a,
  coef_jnu15_a,
  coef_jnu16_a,
  coef_jnu17_a,
  coef_jnu18_a,
  coef_jnu19_a,
  coef_jnu20_a
};

static const size_t size_jnu_a[] = {
  0,
  sizeof(coef_jnu1_a)/sizeof(double),
  sizeof(coef_jnu2_a)/sizeof(double),
  sizeof(coef_jnu3_a)/sizeof(double),
  sizeof(coef_jnu4_a)/sizeof(double),
  sizeof(coef_jnu5_a)/sizeof(double),
  sizeof(coef_jnu6_a)/sizeof(double),
  sizeof(coef_jnu7_a)/sizeof(double),
  sizeof(coef_jnu8_a)/sizeof(double),
  sizeof(coef_jnu9_a)/sizeof(double),
  sizeof(coef_jnu10_a)/sizeof(double),
  sizeof(coef_jnu11_a)/sizeof(double),
  sizeof(coef_jnu12_a)/sizeof(double),
  sizeof(coef_jnu13_a)/sizeof(double),
  sizeof(coef_jnu14_a)/sizeof(double),
  sizeof(coef_jnu15_a)/sizeof(double),
  sizeof(coef_jnu16_a)/sizeof(double),
  sizeof(coef_jnu17_a)/sizeof(double),
  sizeof(coef_jnu18_a)/sizeof(double),
  sizeof(coef_jnu19_a)/sizeof(double),
  sizeof(coef_jnu20_a)/sizeof(double)
};

static const double * coef_jnu_b[] = {
  0,
  coef_jnu1_b,
  coef_jnu2_b,
  coef_jnu3_b,
  coef_jnu4_b,
  coef_jnu5_b,
  coef_jnu6_b,
  coef_jnu7_b,
  coef_jnu8_b,
  coef_jnu9_b,
  coef_jnu10_b
};

static const size_t size_jnu_b[] = {
  0,
  sizeof(coef_jnu1_b)/sizeof(double),
  sizeof(coef_jnu2_b)/sizeof(double),
  sizeof(coef_jnu3_b)/sizeof(double),
  sizeof(coef_jnu4_b)/sizeof(double),
  sizeof(coef_jnu5_b)/sizeof(double),
  sizeof(coef_jnu6_b)/sizeof(double),
  sizeof(coef_jnu7_b)/sizeof(double),
  sizeof(coef_jnu8_b)/sizeof(double),
  sizeof(coef_jnu9_b)/sizeof(double),
  sizeof(coef_jnu10_b)/sizeof(double)
};

/* Evaluate Clenshaw recurrence for
 * a T* Chebyshev series.
 * sizeof(c) = N+1
 */
static double
clenshaw(const double * c, int N, double u)
{
  double B_np1 = 0.0;
  double B_n   = c[N];
  double B_nm1;
  int n;
  for(n=N; n>0; n--) {
    B_nm1 = 2.0*(2.0*u-1.0) * B_n - B_np1 + c[n-1];
    B_np1 = B_n;
    B_n   = B_nm1;
  }
  return B_n - (2.0*u-1.0)*B_np1;
}

/* correction terms to leading McMahon expansion
 * [Abramowitz+Stegun 9.5.12]
 * [Olver, Royal Society Math. Tables, v. 7]
 * We factor out a beta, so that this is a multiplicative
 * correction:
 *   j_{nu,s} = beta(s,nu) * mcmahon_correction(nu, beta(s,nu))
 *   macmahon_correction --> 1 as s --> Inf
 */
static double
mcmahon_correction(const double mu, const double beta)
{
  const double eb   = 8.0*beta;
  const double ebsq = eb*eb;

  if(mu < GSL_DBL_EPSILON) {
    /* Prevent division by zero below. */
    const double term1 =  1.0/ebsq;
    const double term2 = -4.0*31.0/(3*ebsq*ebsq);
    const double term3 =  32.0*3779.0/(15.0*ebsq*ebsq*ebsq);
    const double term4 = -64.0*6277237.0/(105.0*ebsq*ebsq*ebsq*ebsq);
    const double term5 =  512.0*2092163573.0/(315.0*ebsq*ebsq*ebsq*ebsq*ebsq);
    return 1.0 + 8.0*(term1 + term2 + term3 + term4 + term5);
  }
  else {
    /* Here we do things in terms of 1/mu, which
     * is purely to prevent overflow in the very
     * unlikely case that mu is really big.
     */
    const double mi   = 1.0/mu;
    const double r  = mu/ebsq;
    const double n2 = 4.0/3.0    * (7.0 - 31.0*mi);
    const double n3 = 32.0/15.0  * (83.0 + (-982.0 + 3779.0*mi)*mi);
    const double n4 = 64.0/105.0 * (6949.0 + (-153855.0 + (1585743.0 - 6277237.0*mi)*mi)*mi);
    const double n5 = 512.0/315.0 * (70197.0 + (-2479316.0 + (48010494.0 + (-512062548.0 + 2092163573.0*mi)*mi)*mi)*mi);
    const double n6 = 2048.0/3465.0 * (5592657.0 + (-287149133.0 + (8903961290.0 + (-179289628602.0 + (1982611456181.0 - 8249725736393.0*mi)*mi)*mi)*mi)*mi);
    const double term1 = (1.0 - mi) * r;
    const double term2 = term1 * n2 * r;
    const double term3 = term1 * n3 * r*r;
    const double term4 = term1 * n4 * r*r*r;
    const double term5 = term1 * n5 * r*r*r*r;
    const double term6 = term1 * n6 * r*r*r*r*r;
    return 1.0 - 8.0*(term1 + term2 + term3 + term4 + term5 + term6);
  }
}

/* Assumes z >= 1.0 */
static double
olver_b0(double z, double minus_zeta)
{
  if(z < 1.02) {
    const double a = 1.0-z;
    const double c0 =  0.0179988721413553309252458658183;
    const double c1 =  0.0111992982212877614645974276203;
    const double c2 =  0.0059404069786014304317781160605;
    const double c3 =  0.0028676724516390040844556450173;
    const double c4 =  0.0012339189052567271708525111185;
    const double c5 =  0.0004169250674535178764734660248;
    const double c6 =  0.0000330173385085949806952777365;
    const double c7 = -0.0001318076238578203009990106425;
    const double c8 = -0.0001906870370050847239813945647;
    return c0 + a*(c1 + a*(c2 + a*(c3 + a*(c4 + a*(c5 + a*(c6 + a*(c7 + a*c8)))))));
  }
  else {
    const double abs_zeta = minus_zeta;
    const double t = 1.0/(z*sqrt(1.0 - 1.0/(z*z)));
    return -5.0/(48.0*abs_zeta*abs_zeta) + t*(3.0 + 5.0*t*t)/(24.0*sqrt(abs_zeta));
  }
}

inline
static double
olver_f1(double z, double minus_zeta)
{
  const double b0 = olver_b0(z, minus_zeta);
  const double h2 = sqrt(4.0*minus_zeta/(z*z-1.0)); /* FIXME */
  return 0.5 * z * h2 * b0;
}

int
gsl_sf_bessel_zero_J0_e(unsigned int s, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(s == 0){
    result->val = 0.0;
    result->err = 0.0;
    GSL_ERROR ("error", GSL_EINVAL);
  }
  else {
    /* See [F. Lether, J. Comp. Appl .Math. 67, 167 (1996)]. */

    static const double P[] = { 1567450796.0/12539606369.0,
                                8903660.0/2365861.0,
                                10747040.0/536751.0,
                                17590991.0/1696654.0
                              };
    static const double Q[] = { 1.0,
                                29354255.0/954518.0,
                                76900001.0/431847.0,
                                67237052.0/442411.0
                              };

    const double beta = (s - 0.25) * M_PI;
    const double bi2  = 1.0/(beta*beta);
    const double R33num = P[0] + bi2 * (P[1] + bi2 * (P[2] + P[3] * bi2));
    const double R33den = Q[0] + bi2 * (Q[1] + bi2 * (Q[2] + Q[3] * bi2));
    const double R33 = R33num/R33den;
    result->val = beta + R33/beta;
    result->err = fabs(3.0e-15 * result->val);
    return GSL_SUCCESS;
  }
}

int
gsl_sf_bessel_zero_J1_e(unsigned int s, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(s == 0) {
    result->val = 0.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else {
    /* See [M. Branders et al., J. Comp. Phys. 42, 403 (1981)]. */

    static const double a[] = { -0.362804405737084,
                                 0.120341279038597,
                                 0.439454547101171e-01,
                                 0.159340088474713e-02
                              };
    static const double b[] = {  1.0,
                                -0.325641790801361,
                                -0.117453445968927,
                                -0.424906902601794e-02
                              };

    const double beta = (s + 0.25) * M_PI;
    const double bi2  = 1.0/(beta*beta);
    const double Rnum = a[3] + bi2 * (a[2] + bi2 * (a[1] + bi2 * a[0]));
    const double Rden = b[3] + bi2 * (b[2] + bi2 * (b[1] + bi2 * b[0]));
    const double R = Rnum/Rden;
    result->val = beta * (1.0 + R*bi2);
    result->err = fabs(2.0e-14 * result->val);
    return GSL_SUCCESS;
  }
}

int
gsl_sf_bessel_zero_Jnu_e(double nu, unsigned int s, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(nu <= -1.0) {
    DOMAIN_ERROR(result);
  }
  else if(s == 0) {
    result->val = 0.0;
    result->err = 0.0;
    if (nu == 0.0) {
      GSL_ERROR ("no zero-th root for nu = 0.0", GSL_EINVAL);
    }
    return GSL_SUCCESS;
  }
  else if(nu < 0.0) {
    /* This can be done, I'm just lazy now. */
    result->val = 0.0;
    result->err = 0.0;
    GSL_ERROR("unimplemented", GSL_EUNIMPL);
  }
  else if(s == 1) {
    /* Chebyshev fits for the first positive zero.
     * For some reason Nemeth made this different from the others.
     */
    if(nu < 2.0) {
      const double * c = coef_jnu_a[s];
      const size_t   L = size_jnu_a[s];
      const double arg = nu/2.0;
      const double chb = clenshaw(c, L-1, arg);
      result->val = chb;
      result->err = 2.0e-15 * result->val;
    }
    else {
      const double * c = coef_jnu_b[s];
      const size_t   L = size_jnu_b[s];
      const double arg = pow(2.0/nu, 2.0/3.0);
      const double chb = clenshaw(c, L-1, arg);
      result->val = nu * chb;
      result->err = 2.0e-15 * result->val;
    }
    return GSL_SUCCESS;
  }
  else if(s <= 10) {
    /* Chebyshev fits for the first 10 positive zeros. */
    if(nu < s) {
      const double * c = coef_jnu_a[s];
      const size_t   L = size_jnu_a[s];
      const double arg = nu/s;
      const double chb = clenshaw(c, L-1, arg);
      result->val = chb;
      result->err = 2.0e-15 * result->val;
    }
    else {
      const double * c = coef_jnu_b[s];
      const size_t   L = size_jnu_b[s];
      const double arg = pow(s/nu, 2.0/3.0);
      const double chb = clenshaw(c, L-1, arg);
      result->val = nu * chb;
      result->err = 2.0e-15 * result->val;

      /* FIXME: truth in advertising for the screwed up
       * s = 5 fit. Need to fix that.
       */
      if(s == 5) {
        result->err *= 5.0e+06;
      }
    }
    return GSL_SUCCESS;
  }
  else if(s > 0.5*nu && s <= 20) {
    /* Chebyshev fits for 10 < s <= 20. */
    const double * c = coef_jnu_a[s];
    const size_t   L = size_jnu_a[s];
    const double arg = nu/(2.0*s);
    const double chb = clenshaw(c, L-1, arg);
    result->val = chb;
    result->err = 4.0e-15 * chb;
    return GSL_SUCCESS;
  }
  else if(s > 2.0 * nu) {
    /* McMahon expansion if s is large compared to nu. */
    const double beta = (s + 0.5*nu - 0.25) * M_PI;
    const double mc   = mcmahon_correction(4.0*nu*nu, beta);
    gsl_sf_result rat12;
    gsl_sf_pow_int_e(nu/beta, 14, &rat12);
    result->val  = beta * mc;
    result->err  = 4.0 * fabs(beta) * rat12.val;
    result->err += 4.0 * fabs(GSL_DBL_EPSILON * result->val);
    return GSL_SUCCESS;
  }
  else {
    /* Olver uniform asymptotic. */
    gsl_sf_result as;
    const int stat_as = gsl_sf_airy_zero_Ai_e(s, &as);
    const double minus_zeta = -pow(nu,-2.0/3.0) * as.val;
    const double z  = gsl_sf_bessel_Olver_zofmzeta(minus_zeta);
    const double f1 = olver_f1(z, minus_zeta);
    result->val  = nu * (z + f1/(nu*nu));
    result->err  = 0.001/(nu*nu*nu);
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return stat_as;
  }
}

/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"

double gsl_sf_bessel_zero_J0(unsigned int s)
{
  EVAL_RESULT(gsl_sf_bessel_zero_J0_e(s, &result));
}

double gsl_sf_bessel_zero_J1(unsigned int s)
{
  EVAL_RESULT(gsl_sf_bessel_zero_J1_e(s, &result));
}

double gsl_sf_bessel_zero_Jnu(double nu, unsigned int s)
{
  EVAL_RESULT(gsl_sf_bessel_zero_Jnu_e(nu, s, &result));
}


text:
#pragma once

#include <gsl\gsl>
#include <winrt\base.h>
#include <d3d11.h>
#include "DrawableGameComponent.h"
#include "BasicTessellationMaterial.h"
#include "RenderStateHelper.h"

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

		bool UseUniformTessellation() const;
		void SetUseUniformTessellation(bool useUniformTessellation);
		void ToggleUseUniformTessellation();
		bool ShowQuadTopology() const;
		void SetShowQuadTopology(bool showQuadTopology);
		void ToggleTopology();

		gsl::span<const float> EdgeFactors() const;
		void SetUniformEdgeFactors(float factor);
		void SetEdgeFactor(float factor, std::uint32_t index);
		void SetInsideFactor(float factor, std::uint32_t index);

		gsl::span<const float> InsideFactors() const;

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

	private:
		/*void UpdateEdgeFactors();
		void UpdateInsideEdgeFactors();*/

		Library::RenderStateHelper mRenderStateHelper;
		BasicTessellationMaterial mMaterial;
		winrt::com_ptr<ID3D11Buffer> mTriVertexBuffer;
		winrt::com_ptr<ID3D11Buffer> mQuadVertexBuffer;
		bool mUpdateMaterial{ true };
		bool mUseUniformTessellation{ true };
	};
}


text:
/* -----------------------------------------------------------------------------
 * Copyright 2021 Jonathan Haigh
 * SPDX-License-Identifier: MIT
 * ---------------------------------------------------------------------------*/

#ifndef SQ_INCLUDE_GUARD_core_narrow_h_
#define SQ_INCLUDE_GUARD_core_narrow_h_

#include "core/typeutil.h"

#include <concepts>
#include <cstddef>
#include <gsl/gsl>
#include <string_view>

namespace sq {

/**
 * Wrapper around gsl::narrow that provides better diags on failure.
 *
 * @param value value to convert
 * @param format_args optional arguments, passed to fmt::format that describe
 *        what kind of thing is being converted.
 *
 * E.g.
 *    narrow<int>(0.5, "number of words in file {}", path)
 * might throw an exception with a message like:
 *    number of words in file "/some/file" (0.5) does not fit in type int
 *
 * If the optional format_args args aren't given then the value's type is used
 * as a description instead. E.g.
 *    narrow<int>(0.5)
 * might throw an exception with a message like:
 *    float (0.5) does not fit in type int
 */
template <typename T, typename U, typename... FormatArgs>
SQ_ND constexpr T narrow(U value, FormatArgs &&...format_args);

/**
 * Convert an integral value to a gsl::index.
 *
 * Throws if the value can't be represented by a gsl::index.
 */
SQ_ND constexpr gsl::index to_index(auto value, auto &&...format_args);

/**
 * Convert an integral value to a std::size_t.
 *
 * Throws if the value can't be represented by a std::size_t.
 */
SQ_ND constexpr std::size_t to_size(auto value, auto &&...format_args);

} // namespace sq

#include "core/narrow.inl.h"

#endif // SQ_INCLUDE_GUARD_core_narrow_h_

max_stars_repo_name: jonathanhaigh/sq

max_issues_repo_name: jonathanhaigh/sq

max_forks_repo_name: jonathanhaigh/sq


text:
// Copyright (c) Microsoft Corporation.
// Licensed under the MIT License.

#pragma once

#include "arcana/bob/bob_helpers.h"

#include <cstdint>
#include <string>
#include <vector>

#include <cereal/cereal.hpp>
#include <cereal/types/vector.hpp>

#include <gsl/gsl>

namespace mira
{
    namespace bob
    {
        namespace v1
        {
            namespace
            {
                constexpr std::uint32_t VERSION = 1;
            }

            //
            // manifest and stream_file_header are both types found in the json
            // manifest file. This file describes in a human readable format what the bob
            // stream contains and which data it depends on.
            //
            // entry_instance is the header data of each instance of an entry stream in the
            // binary stream file (.bob).
            //

            //
            //  Describes what the binary stream file contains
            //
            class stream_file_header
            {
            public:
                //
                //  An entry describes the properties of one type of packet found in the stream
                //
                class entry
                {
                public:
                    entry()
                        : m_name{}
                        , m_type{}
                        , m_version{}
                    {}

                    entry(std::string name, ::uint32_t version)
                        : m_name{ std::move(name) }
                        , m_type{}
                        , m_version{ version }
                    {}

                    entry(std::string name, std::string type, ::uint32_t version)
                        : m_name{ std::move(name) }
                        , m_type{ std::move(type) }
                        , m_version{ version }
                    {}

                    const std::string& name() const
                    {
                        return m_name;
                    }

                    const std::string& type() const
                    {
                        return m_type;
                    }

                    std::uint32_t version() const
                    {
                        return m_version;
                    }

                    template<typename Archive>
                    void serialize(Archive& ar)
                    {
                        ar(cereal::make_nvp("Name", m_name), cereal::make_nvp("Type", m_type), cereal::make_nvp("Version", m_version));
                    }

                private:
                    std::string m_name;
                    std::string m_type;
                    std::uint32_t m_version;
                };

                stream_file_header() = default;

                explicit stream_file_header(std::vector<entry> entries)
                    : m_entries(std::move(entries))
                {}

                std::uint32_t version() const
                {
                    return VERSION;
                }

                const std::vector<entry>& entries() const
                {
                    return m_entries;
                }

                template<typename Archive>
                void serialize(Archive& ar, const std::uint32_t version)
                {
                    if (version != VERSION)
                        throw_invalid_version();

                    ar(cereal::make_nvp("Entries", m_entries));
                }

            private:
                std::vector<entry> m_entries;
            };

            class manifest
            {
            public:
                manifest() = default;

                manifest(stream_file_header header)
                    : m_stream{ std::move(header) }
                {}

                std::uint32_t version() const
                {
                    return VERSION;
                }

                const stream_file_header& stream() const
                {
                    return m_stream;
                }

                template<typename Archive>
                void serialize(Archive& ar, const std::uint32_t version)
                {
                    if (version != VERSION)
                        throw_invalid_version();

                    ar(cereal::make_nvp("Stream", m_stream));
                }

            private:
                stream_file_header m_stream;
            };

            //
            // describes an entry in the binary stream, the data for the entry
            // is contained between this block and the next.
            //
            struct entry_instance
            {
                std::int64_t Timestamp;
                std::int64_t NextBlockOffset;
                std::int32_t EntryIdx;

                template<typename Archive>
                void serialize(Archive& ar, const std::uint32_t version)
                {
                    if (version != VERSION)
                        throw_invalid_version();

                    ar(CEREAL_NVP(Timestamp), CEREAL_NVP(NextBlockOffset), CEREAL_NVP(EntryIdx));
                }
            };
        }

        using namespace v1;
    }
}

CEREAL_CLASS_VERSION(::mira::bob::v1::manifest, ::mira::bob::v1::VERSION);
CEREAL_CLASS_VERSION(::mira::bob::v1::stream_file_header, ::mira::bob::v1::VERSION);

max_stars_repo_name: syntheticmagus/mageslam

max_issues_repo_name: syntheticmagus/mageslam

max_forks_repo_name: syntheticmagus/mageslam


text:
#include <cblas.h>
#include <stdlib.h>
#include <string.h>
#include "loss.h"

void expit(const double *x, double *out, int x_len) {
    for (int i = 0; i < x_len; i++) {
        if (x[i] > 0) {
            out[i] = 1. / (1. + exp(-x[i]));
        } else {
            out[i] = 1. - 1. / (1. + exp(x[i]));
        }
    }
}

void log_logistic(const double *x, double *out, int x_len) {
    for (int i = 0; i < x_len; i++) {
        if (x[i] > 0.0) {
            out[i] = -log(1.0 + exp(-x[i]));
        } else {
            out[i] = x[i] - log(1.0 + exp(x[i]));
        }
    }
}

void logistic(const double *x, double *out, int x_len) {
    for (int i = 0; i < x_len; i++) {
        if (x[i] > 0) {
            out[i] = 1. / (1. + exp(-x[i]));
        } else {
            out[i] = 1. - 1. / (1. + exp(x[i]));
        }
    }
}

double log_sum_exp(const double *x, int x_len) {
    double max_x = x[0], out = 0.0;
    for (int i = 1; i < x_len; i++) {
        if (x[i] > max_x) {
            max_x = x[i];
        }
    }
    for (int i = 0; i < x_len; i++) {
        out += exp(x[i] - max_x);
    }
    return max_x + log(out);
}

void logistic_loss_grad(const double *w,
                        const double *x_tr,
                        const double *y_tr,
                        double *loss_grad,
                        double eta,
                        int n_samples,
                        int n_features) {
    for (int i = 0; i < n_features + 1; i++) {
        if (isnan(w[i])) {
            printf("%f\n", w[i]);
            printf("warning: loss grad error!\n");
            break;
        }
    }
    int i, n = n_samples, p = n_features;
    double intercept = w[p], sum_z0 = 0.0;
    loss_grad[0] = 0.0;
    double *yz = malloc(sizeof(double) * n);
    double *z0 = malloc(sizeof(double) * n);
    double *logistic = malloc(sizeof(double) * n);
    for (i = 0; i < n; i++) { /** calculate yz */
        yz[i] = intercept;
    }
    //x_tr^T*w+
    cblas_dgemv(CblasRowMajor, CblasNoTrans, n, p, 1., x_tr, p, w, 1, 1., yz, 1);
    for (i = 0; i < n; i++) {
        yz[i] *= y_tr[i];
    }
    expit(yz, z0, n); // calculate z0 and final intercept
    /** calculate logistic logistic[i] = 1/(1+exp(-y[i]*(xi^T*w+c)))*/
    log_logistic(yz, logistic, n);
    /** calculate loss of data fitting part*/
    for (i = 0; i < n; i++) {
        z0[i] = (z0[i] - 1.) * y_tr[i];
        sum_z0 += z0[i];
        loss_grad[0] -= logistic[i];
    }
    /**calculate loss of regularization part (it does not have intercept)*/
    loss_grad[0] += 0.5 * eta * cblas_ddot(p, w, 1, w, 1);
    /** calculate gradient of coefficients*/
    memcpy(loss_grad + 1, w, sizeof(double) * p);
    /** x^T*z0 + eta*w, where z0[i]=(logistic[i] - 1.)*yi*/
    cblas_dgemv(CblasRowMajor, CblasTrans,
                n, p, 1., x_tr, p, z0, 1, eta, loss_grad + 1, 1);
    /** calculate gradient of intercept part*/
    loss_grad[p + 1] = sum_z0; // intercept part
    free(logistic);
    free(z0);
    free(yz);
}

void logistic_loss_grad_sparse(const double *w,
                               const double *x_tr,
                               const double *y_tr,
                               double *loss_grad,
                               double eta,
                               int n_samples,
                               int n_features) {
    //sub_p is number of nonzeros in w.
    int i, n = n_samples, p = n_features, sub_p = 0;
    double intercept = w[p], sum_z0 = 0.0;
    loss_grad[0] = 0.0;
    double *yz = malloc(sizeof(double) * n);
    double *z0 = malloc(sizeof(double) * n);
    double *logistic = malloc(sizeof(double) * n);
    for (i = 0; i < n; i++) { /** calculate yz */
        yz[i] = intercept;
    }
    double nonzero_w[p];
    int nonzero_w_ind[p];
    for (i = 0; i < p; i++) {
        if (w[i] != 0.0) {
            nonzero_w_ind[sub_p] = i;
            nonzero_w[sub_p] = w[i];
            sub_p++;
        }
    }
    //x_tr^T*w+ to use the sparsity of w
    for (i = 0; i < n; i++) {
        double tmp_val = 0.0;
        for (int j = 0; j < sub_p; j++) {
            tmp_val += x_tr[i * p + nonzero_w_ind[j]] * nonzero_w[j];
        }
        yz[i] = yz[i] + tmp_val;
    }
    for (i = 0; i < n; i++) {
        yz[i] *= y_tr[i];
    }
    expit(yz, z0, n); // calculate z0 and final intercept
    /** calculate logistic logistic[i] = 1/(1+exp(-y[i]*(xi^T*w+c)))*/
    log_logistic(yz, logistic, n);
    /** calculate loss of data fitting part*/
    for (i = 0; i < n; i++) {
        z0[i] = (z0[i] - 1.) * y_tr[i];
        sum_z0 += z0[i];
        loss_grad[0] -= logistic[i];
    }
    /**calculate loss of regularization part (it does not have intercept)*/
    loss_grad[0] += 0.5 * eta * cblas_ddot(p, w, 1, w, 1);
    /** calculate gradient of coefficients*/
    memcpy(loss_grad + 1, w, sizeof(double) * p);
    /** x^T*z0 + eta*w, where z0[i]=(logistic[i] - 1.)*yi*/
    cblas_dgemv(CblasRowMajor, CblasTrans,
                n, p, 1., x_tr, p, z0, 1, eta, loss_grad + 1, 1);
    /** calculate gradient of intercept part*/
    loss_grad[p + 1] = sum_z0; // intercept part
    free(logistic);
    free(z0);
    free(yz);
}

void logistic_predict(const double *x_te,
                      const double *wt,
                      double *pred_prob,
                      double *pred_label,
                      double threshold,
                      int n,
                      int p) {
    openblas_set_num_threads(1);
    int i;
    cblas_dgemv(CblasRowMajor, CblasNoTrans,
                n, p, 1., x_te, p, wt, 1, 0., pred_prob, 1);
    for (i = 0; i < n; i++) {
        pred_prob[i] += wt[p]; // intercept
    }
    expit(pred_prob, pred_prob, n);
    for (i = 0; i < n; i++) {
        if (pred_prob[i] >= threshold) {
            pred_label[i] = 1.;
        } else {
            pred_label[i] = -1.;
        }
    }
}

void least_square_loss_grad(const double *w,
                            const double *x_tr,
                            const double *y_tr,
                            double *loss_grad,
                            double eta,
                            int n_samples,
                            int n_features) {
    for (int i = 0; i < n_features + 1; i++) {
        if (isnan(w[i])) {
            printf("%f\n", w[i]);
            printf("warning: loss grad error!\n");
            break;
        }
    }
    int i, n = n_samples, p = n_features;
    double *yz = malloc(sizeof(double) * n);
    double *w0 = malloc(sizeof(double) * n);
    cblas_dcopy(n, y_tr, 1, yz, 1); // y_tr --> yz
    cblas_dscal(p + 2, 0.0, loss_grad, 1);
    cblas_daxpy(p, eta, w, 1, loss_grad + 1, 1); // l2_lambda*w --> loss_grad
    //Order,TransA, M, N, alpha, A, lda, x, incX, beta, Y, incY
    // Y<-alpha*AX + beta*Y, where A=MxN, Y=N
    cblas_dgemv(CblasRowMajor, CblasNoTrans,
                n, p, 1., x_tr, p, w, 1, -1., yz, 1); //Xw - y
    for (i = 0; i < n; i++) { w0[i] = w[p]; }
    cblas_daxpy(n, 1., w0, 1, yz, 1);
    loss_grad[0] = cblas_ddot(n, yz, 1, yz, 1) / (2. * n);
    loss_grad[0] += (.5 * eta) * cblas_ddot(p, w, 1, w, 1); //loss
    cblas_dgemv(CblasRowMajor, CblasTrans,
                n, p, 1. / n, x_tr, p, yz, 1, 1., loss_grad + 1, 1);
    for (i = 0; i < n; i++) {
        w0[i] = w[p];
        loss_grad[p + 1] = yz[i];
    }
    loss_grad[p + 1] /= n;
    free(w0), free(yz);
}

void least_square_predict(const double *x_te,
                          const double *wt,
                          double *pred_prob,
                          double *pred_label,
                          double threshold,
                          int n_samples,
                          int p_features) {
    openblas_set_num_threads(1);
    int i, n = n_samples, p = p_features;
    cblas_dgemv(CblasRowMajor, CblasNoTrans,
                n, p, 1., x_te, p, wt, 1, 0., pred_prob, 1);
    for (i = 0; i < n; i++) {
        pred_prob[i] += wt[p]; // intercept
    }
    // this is not useful.
    for (i = 0; i < n; i++) {
        if (pred_prob[i] >= threshold) {
            pred_label[i] = 1.;
        } else {
            pred_label[i] = -1.;
        }
    }
}

max_stars_repo_name: baojianzhou/sparse-auc

max_issues_repo_name: baojianzhou/sparse-auc

max_forks_repo_name: baojianzhou/sparse-auc


text:
/* matrix/gsl_matrix_float.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman, Brian Gough
 *
 * 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.
 */

#ifndef __GSL_MATRIX_FLOAT_H__
#define __GSL_MATRIX_FLOAT_H__

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_inline.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_float.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

typedef struct
{
  size_t size1;
  size_t size2;
  size_t tda;
  float * data;
  gsl_block_float * block;
  int owner;
} gsl_matrix_float;

typedef struct
{
  gsl_matrix_float matrix;
} _gsl_matrix_float_view;

typedef _gsl_matrix_float_view gsl_matrix_float_view;

typedef struct
{
  gsl_matrix_float matrix;
} _gsl_matrix_float_const_view;

typedef const _gsl_matrix_float_const_view gsl_matrix_float_const_view;

/* Allocation */

gsl_matrix_float *
gsl_matrix_float_alloc (const size_t n1, const size_t n2);

gsl_matrix_float *
gsl_matrix_float_calloc (const size_t n1, const size_t n2);

gsl_matrix_float *
gsl_matrix_float_alloc_from_block (gsl_block_float * b,
                                   const size_t offset,
                                   const size_t n1,
                                   const size_t n2,
                                   const size_t d2);

gsl_matrix_float *
gsl_matrix_float_alloc_from_matrix (gsl_matrix_float * m,
                                    const size_t k1,
                                    const size_t k2,
                                    const size_t n1,
                                    const size_t n2);

gsl_vector_float *
gsl_vector_float_alloc_row_from_matrix (gsl_matrix_float * m,
                                        const size_t i);

gsl_vector_float *
gsl_vector_float_alloc_col_from_matrix (gsl_matrix_float * m,
                                        const size_t j);

void gsl_matrix_float_free (gsl_matrix_float * m);

/* Views */

_gsl_matrix_float_view
gsl_matrix_float_submatrix (gsl_matrix_float * m,
                            const size_t i, const size_t j,
                            const size_t n1, const size_t n2);

_gsl_vector_float_view
gsl_matrix_float_row (gsl_matrix_float * m, const size_t i);

_gsl_vector_float_view
gsl_matrix_float_column (gsl_matrix_float * m, const size_t j);

_gsl_vector_float_view
gsl_matrix_float_diagonal (gsl_matrix_float * m);

_gsl_vector_float_view
gsl_matrix_float_subdiagonal (gsl_matrix_float * m, const size_t k);

_gsl_vector_float_view
gsl_matrix_float_superdiagonal (gsl_matrix_float * m, const size_t k);

_gsl_vector_float_view
gsl_matrix_float_subrow (gsl_matrix_float * m, const size_t i,
                         const size_t offset, const size_t n);

_gsl_vector_float_view
gsl_matrix_float_subcolumn (gsl_matrix_float * m, const size_t j,
                            const size_t offset, const size_t n);

_gsl_matrix_float_view
gsl_matrix_float_view_array (float * base,
                             const size_t n1,
                             const size_t n2);

_gsl_matrix_float_view
gsl_matrix_float_view_array_with_tda (float * base,
                                      const size_t n1,
                                      const size_t n2,
                                      const size_t tda);

_gsl_matrix_float_view
gsl_matrix_float_view_vector (gsl_vector_float * v,
                              const size_t n1,
                              const size_t n2);

_gsl_matrix_float_view
gsl_matrix_float_view_vector_with_tda (gsl_vector_float * v,
                                       const size_t n1,
                                       const size_t n2,
                                       const size_t tda);

_gsl_matrix_float_const_view
gsl_matrix_float_const_submatrix (const gsl_matrix_float * m,
                                  const size_t i, const size_t j,
                                  const size_t n1, const size_t n2);

_gsl_vector_float_const_view
gsl_matrix_float_const_row (const gsl_matrix_float * m,
                            const size_t i);

_gsl_vector_float_const_view
gsl_matrix_float_const_column (const gsl_matrix_float * m,
                               const size_t j);

_gsl_vector_float_const_view
gsl_matrix_float_const_diagonal (const gsl_matrix_float * m);

_gsl_vector_float_const_view
gsl_matrix_float_const_subdiagonal (const gsl_matrix_float * m,
                                    const size_t k);

_gsl_vector_float_const_view
gsl_matrix_float_const_superdiagonal (const gsl_matrix_float * m,
                                      const size_t k);

_gsl_vector_float_const_view
gsl_matrix_float_const_subrow (const gsl_matrix_float * m, const size_t i,
                               const size_t offset, const size_t n);

_gsl_vector_float_const_view
gsl_matrix_float_const_subcolumn (const gsl_matrix_float * m, const size_t j,
                                  const size_t offset, const size_t n);

_gsl_matrix_float_const_view
gsl_matrix_float_const_view_array (const float * base,
                                   const size_t n1,
                                   const size_t n2);

_gsl_matrix_float_const_view
gsl_matrix_float_const_view_array_with_tda (const float * base,
                                            const size_t n1,
                                            const size_t n2,
                                            const size_t tda);

_gsl_matrix_float_const_view
gsl_matrix_float_const_view_vector (const gsl_vector_float * v,
                                    const size_t n1,
                                    const size_t n2);

_gsl_matrix_float_const_view
gsl_matrix_float_const_view_vector_with_tda (const gsl_vector_float * v,
                                             const size_t n1,
                                             const size_t n2,
                                             const size_t tda);

/* Operations */

void gsl_matrix_float_set_zero (gsl_matrix_float * m);
void gsl_matrix_float_set_identity (gsl_matrix_float * m);
void gsl_matrix_float_set_all (gsl_matrix_float * m, float x);

int gsl_matrix_float_fread (FILE * stream, gsl_matrix_float * m) ;
int gsl_matrix_float_fwrite (FILE * stream, const gsl_matrix_float * m) ;
int gsl_matrix_float_fscanf (FILE * stream, gsl_matrix_float * m);
int gsl_matrix_float_fprintf (FILE * stream, const gsl_matrix_float * m, const char * format);

int gsl_matrix_float_memcpy(gsl_matrix_float * dest, const gsl_matrix_float * src);
int gsl_matrix_float_swap(gsl_matrix_float * m1, gsl_matrix_float * m2);

int gsl_matrix_float_swap_rows(gsl_matrix_float * m, const size_t i, const size_t j);
int gsl_matrix_float_swap_columns(gsl_matrix_float * m, const size_t i, const size_t j);
int gsl_matrix_float_swap_rowcol(gsl_matrix_float * m, const size_t i, const size_t j);
int gsl_matrix_float_transpose (gsl_matrix_float * m);
int gsl_matrix_float_transpose_memcpy (gsl_matrix_float * dest, const gsl_matrix_float * src);

float gsl_matrix_float_max (const gsl_matrix_float * m);
float gsl_matrix_float_min (const gsl_matrix_float * m);
void gsl_matrix_float_minmax (const gsl_matrix_float * m, float * min_out, float * max_out);

void gsl_matrix_float_max_index (const gsl_matrix_float * m, size_t * imax, size_t *jmax);
void gsl_matrix_float_min_index (const gsl_matrix_float * m, size_t * imin, size_t *jmin);
void gsl_matrix_float_minmax_index (const gsl_matrix_float * m, size_t * imin, size_t * jmin, size_t * imax, size_t * jmax);

int gsl_matrix_float_equal (const gsl_matrix_float * a, const gsl_matrix_float * b);

int gsl_matrix_float_isnull (const gsl_matrix_float * m);
int gsl_matrix_float_ispos (const gsl_matrix_float * m);
int gsl_matrix_float_isneg (const gsl_matrix_float * m);
int gsl_matrix_float_isnonneg (const gsl_matrix_float * m);

int gsl_matrix_float_add (gsl_matrix_float * a, const gsl_matrix_float * b);
int gsl_matrix_float_sub (gsl_matrix_float * a, const gsl_matrix_float * b);
int gsl_matrix_float_mul_elements (gsl_matrix_float * a, const gsl_matrix_float * b);
int gsl_matrix_float_div_elements (gsl_matrix_float * a, const gsl_matrix_float * b);
int gsl_matrix_float_scale (gsl_matrix_float * a, const double x);
int gsl_matrix_float_add_constant (gsl_matrix_float * a, const double x);
int gsl_matrix_float_add_diagonal (gsl_matrix_float * a, const double x);

/***********************************************************************/
/* The functions below are obsolete                                    */
/***********************************************************************/
int gsl_matrix_float_get_row(gsl_vector_float * v, const gsl_matrix_float * m, const size_t i);
int gsl_matrix_float_get_col(gsl_vector_float * v, const gsl_matrix_float * m, const size_t j);
int gsl_matrix_float_set_row(gsl_matrix_float * m, const size_t i, const gsl_vector_float * v);
int gsl_matrix_float_set_col(gsl_matrix_float * m, const size_t j, const gsl_vector_float * v);
/***********************************************************************/

/* inline functions if you are using GCC */

INLINE_DECL float   gsl_matrix_float_get(const gsl_matrix_float * m, const size_t i, const size_t j);
INLINE_DECL void    gsl_matrix_float_set(gsl_matrix_float * m, const size_t i, const size_t j, const float x);
INLINE_DECL float * gsl_matrix_float_ptr(gsl_matrix_float * m, const size_t i, const size_t j);
INLINE_DECL const float * gsl_matrix_float_const_ptr(const gsl_matrix_float * m, const size_t i, const size_t j);

#ifdef HAVE_INLINE
INLINE_FUN
float
gsl_matrix_float_get(const gsl_matrix_float * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_VAL("first index out of range", GSL_EINVAL, 0) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_VAL("second index out of range", GSL_EINVAL, 0) ;
        }
    }
#endif
  return m->data[i * m->tda + j] ;
}

INLINE_FUN
void
gsl_matrix_float_set(gsl_matrix_float * m, const size_t i, const size_t j, const float x)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_VOID("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_VOID("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  m->data[i * m->tda + j] = x ;
}

INLINE_FUN
float *
gsl_matrix_float_ptr(gsl_matrix_float * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  return (float *) (m->data + (i * m->tda + j)) ;
}

INLINE_FUN
const float *
gsl_matrix_float_const_ptr(const gsl_matrix_float * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  return (const float *) (m->data + (i * m->tda + j)) ;
}

#endif

__END_DECLS

#endif /* __GSL_MATRIX_FLOAT_H__ */

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: ielomariala/Hex-Game