text:
/* specfunc/airy.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 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.
 */

/* Author:  G. Jungman */

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

#include "error.h"
#include "check.h"

#include "chebyshev.h"
#include "cheb_eval_mode.c"

/*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/

/* chebyshev expansions for Airy modulus and phase
   based on SLATEC r9aimp()

 Series for AM21       on the interval -1.25000D-01 to  0.
                                        with weighted error   2.89E-17
                                         log weighted error  16.54
                               significant figures required  14.15
                                    decimal places required  17.34

 Series for ATH1       on the interval -1.25000D-01 to  0.
                                        with weighted error   2.53E-17
                                         log weighted error  16.60
                               significant figures required  15.15
                                    decimal places required  17.38

 Series for AM22       on the interval -1.00000D+00 to -1.25000D-01
                                        with weighted error   2.99E-17
                                         log weighted error  16.52
                               significant figures required  14.57
                                    decimal places required  17.28

 Series for ATH2       on the interval -1.00000D+00 to -1.25000D-01
                                        with weighted error   2.57E-17
                                         log weighted error  16.59
                               significant figures required  15.07
                                    decimal places required  17.34
*/

static double am21_data[37] = {
  0.0065809191761485,
  0.0023675984685722,
  0.0001324741670371,
  0.0000157600904043,
  0.0000027529702663,
  0.0000006102679017,
  0.0000001595088468,
  0.0000000471033947,
  0.0000000152933871,
  0.0000000053590722,
  0.0000000020000910,
  0.0000000007872292,
  0.0000000003243103,
  0.0000000001390106,
  0.0000000000617011,
  0.0000000000282491,
  0.0000000000132979,
  0.0000000000064188,
  0.0000000000031697,
  0.0000000000015981,
  0.0000000000008213,
  0.0000000000004296,
  0.0000000000002284,
  0.0000000000001232,
  0.0000000000000675,
  0.0000000000000374,
  0.0000000000000210,
  0.0000000000000119,
  0.0000000000000068,
  0.0000000000000039,
  0.0000000000000023,
  0.0000000000000013,
  0.0000000000000008,
  0.0000000000000005,
  0.0000000000000003,
  0.0000000000000001,
  0.0000000000000001
};
static cheb_series am21_cs = {
  am21_data,
  36,
  -1, 1,
  20
};

static double ath1_data[36] = {
  -0.07125837815669365,
  -0.00590471979831451,
  -0.00012114544069499,
  -0.00000988608542270,
  -0.00000138084097352,
  -0.00000026142640172,
  -0.00000006050432589,
  -0.00000001618436223,
  -0.00000000483464911,
  -0.00000000157655272,
  -0.00000000055231518,
  -0.00000000020545441,
  -0.00000000008043412,
  -0.00000000003291252,
  -0.00000000001399875,
  -0.00000000000616151,
  -0.00000000000279614,
  -0.00000000000130428,
  -0.00000000000062373,
  -0.00000000000030512,
  -0.00000000000015239,
  -0.00000000000007758,
  -0.00000000000004020,
  -0.00000000000002117,
  -0.00000000000001132,
  -0.00000000000000614,
  -0.00000000000000337,
  -0.00000000000000188,
  -0.00000000000000105,
  -0.00000000000000060,
  -0.00000000000000034,
  -0.00000000000000020,
  -0.00000000000000011,
  -0.00000000000000007,
  -0.00000000000000004,
  -0.00000000000000002
};
static cheb_series ath1_cs = {
  ath1_data,
  35,
  -1, 1,
  15
};

static double am22_data[33] = {
 -0.01562844480625341,
  0.00778336445239681,
  0.00086705777047718,
  0.00015696627315611,
  0.00003563962571432,
  0.00000924598335425,
  0.00000262110161850,
  0.00000079188221651,
  0.00000025104152792,
  0.00000008265223206,
  0.00000002805711662,
  0.00000000976821090,
  0.00000000347407923,
  0.00000000125828132,
  0.00000000046298826,
  0.00000000017272825,
  0.00000000006523192,
  0.00000000002490471,
  0.00000000000960156,
  0.00000000000373448,
  0.00000000000146417,
  0.00000000000057826,
  0.00000000000022991,
  0.00000000000009197,
  0.00000000000003700,
  0.00000000000001496,
  0.00000000000000608,
  0.00000000000000248,
  0.00000000000000101,
  0.00000000000000041,
  0.00000000000000017,
  0.00000000000000007,
  0.00000000000000002
};
static cheb_series am22_cs = {
  am22_data,
  32,
  -1, 1,
  15
};

static double ath2_data[32] = {
   0.00440527345871877,
  -0.03042919452318455,
  -0.00138565328377179,
  -0.00018044439089549,
  -0.00003380847108327,
  -0.00000767818353522,
  -0.00000196783944371,
  -0.00000054837271158,
  -0.00000016254615505,
  -0.00000005053049981,
  -0.00000001631580701,
  -0.00000000543420411,
  -0.00000000185739855,
  -0.00000000064895120,
  -0.00000000023105948,
  -0.00000000008363282,
  -0.00000000003071196,
  -0.00000000001142367,
  -0.00000000000429811,
  -0.00000000000163389,
  -0.00000000000062693,
  -0.00000000000024260,
  -0.00000000000009461,
  -0.00000000000003716,
  -0.00000000000001469,
  -0.00000000000000584,
  -0.00000000000000233,
  -0.00000000000000093,
  -0.00000000000000037,
  -0.00000000000000015,
  -0.00000000000000006,
  -0.00000000000000002
};
static cheb_series ath2_cs = {
  ath2_data,
  31,
  -1, 1,
  16
};

/* Airy modulus and phase for x < -1 */
static
int
airy_mod_phase(const double x, gsl_mode_t mode, gsl_sf_result * mod, gsl_sf_result * phase)
{
  gsl_sf_result result_m;
  gsl_sf_result result_p;
  double m, p;
  double sqx;

  if(x < -2.0) {
    double z = 16.0/(x*x*x) + 1.0;
    cheb_eval_mode_e(&am21_cs, z, mode, &result_m);
    cheb_eval_mode_e(&ath1_cs, z, mode, &result_p);
  }
  else if(x <= -1.0) {
    double z = (16.0/(x*x*x) + 9.0)/7.0;
    cheb_eval_mode_e(&am22_cs, z, mode, &result_m);
    cheb_eval_mode_e(&ath2_cs, z, mode, &result_p);
  }
  else {
    mod->val = 0.0;
    mod->err = 0.0;
    phase->val = 0.0;
    phase->err = 0.0;
    GSL_ERROR ("x is greater than 1.0", GSL_EDOM);
  }

  m =  0.3125 + result_m.val;
  p = -0.625  + result_p.val;

  sqx = sqrt(-x);

  mod->val   = sqrt(m/sqx);
  mod->err  = fabs(mod->val) * (GSL_DBL_EPSILON + fabs(result_m.err/result_m.val));
  phase->val = M_PI_4 - x*sqx * p;
  phase->err = fabs(phase->val) * (GSL_DBL_EPSILON + fabs(result_p.err/result_p.val));

  return GSL_SUCCESS;
}

/* Chebyshev series for Ai(x) with x in [-1,1]
   based on SLATEC ai(x)

 series for aif        on the interval -1.00000d+00 to  1.00000d+00
                                   with weighted error   1.09e-19
                                    log weighted error  18.96
                          significant figures required  17.76
                               decimal places required  19.44

 series for aig        on the interval -1.00000d+00 to  1.00000d+00
                                   with weighted error   1.51e-17
                                    log weighted error  16.82
                          significant figures required  15.19
                               decimal places required  17.27
 */
static double ai_data_f[9] = {
  -0.03797135849666999750,
   0.05919188853726363857,
   0.00098629280577279975,
   0.00000684884381907656,
   0.00000002594202596219,
   0.00000000006176612774,
   0.00000000000010092454,
   0.00000000000000012014,
   0.00000000000000000010
};
static cheb_series aif_cs = {
  ai_data_f,
  8,
  -1, 1,
  8
};

static double ai_data_g[8] = {
   0.01815236558116127,
   0.02157256316601076,
   0.00025678356987483,
   0.00000142652141197,
   0.00000000457211492,
   0.00000000000952517,
   0.00000000000001392,
   0.00000000000000001
};
static cheb_series aig_cs = {
  ai_data_g,
  7,
  -1, 1,
  7
};

/* Chebvyshev series for Bi(x) with x in [-1,1]
   based on SLATEC bi(x)

 series for bif        on the interval -1.00000d+00 to  1.00000d+00
                                        with weighted error   1.88e-19
                                         log weighted error  18.72
                               significant figures required  17.74
                                    decimal places required  19.20

 series for big        on the interval -1.00000d+00 to  1.00000d+00
                                        with weighted error   2.61e-17
                                         log weighted error  16.58
                               significant figures required  15.17
                                    decimal places required  17.03
 */
static double data_bif[9] = {
  -0.01673021647198664948,
   0.10252335834249445610,
   0.00170830925073815165,
   0.00001186254546774468,
   0.00000004493290701779,
   0.00000000010698207143,
   0.00000000000017480643,
   0.00000000000000020810,
   0.00000000000000000018
};
static cheb_series bif_cs = {
  data_bif,
  8,
  -1, 1,
  8
};

static double data_big[8] = {
   0.02246622324857452,
   0.03736477545301955,
   0.00044476218957212,
   0.00000247080756363,
   0.00000000791913533,
   0.00000000001649807,
   0.00000000000002411,
   0.00000000000000002
};
static cheb_series big_cs = {
  data_big,
  7,
  -1, 1,
  7
};

/* Chebyshev series for Bi(x) with x in [1,8]
   based on SLATEC bi(x)
 */
static double data_bif2[10] = {
  0.0998457269381604100,
  0.4786249778630055380,
  0.0251552119604330118,
  0.0005820693885232645,
  0.0000074997659644377,
  0.0000000613460287034,
  0.0000000003462753885,
  0.0000000000014288910,
  0.0000000000000044962,
  0.0000000000000000111
};
static cheb_series bif2_cs = {
  data_bif2,
  9,
  -1, 1,
  9
};

static double data_big2[10] = {
  0.033305662145514340,
  0.161309215123197068,
  0.0063190073096134286,
  0.0001187904568162517,
  0.0000013045345886200,
  0.0000000093741259955,
  0.0000000000474580188,
  0.0000000000001783107,
  0.0000000000000005167,
  0.0000000000000000011
};
static cheb_series big2_cs = {
  data_big2,
  9,
  -1, 1,
  9
};

/* chebyshev for Ai(x) asymptotic factor
   based on SLATEC aie()

 Series for AIP        on the interval  0.          to  1.00000D+00
                   with weighted error   5.10E-17
                    log weighted error  16.29
          significant figures required  14.41
               decimal places required  17.06

 [GJ] Sun Apr 19 18:14:31 EDT 1998
 There was something wrong with these coefficients. I was getting
 errors after 3 or 4 digits. So I recomputed this table. Now I get
 double precision agreement with Mathematica. But it does not seem
 possible that the small differences here would account for the
 original discrepancy. There must have been something wrong with my
 original usage...
*/
static double data_aip[36] = {
 -0.0187519297793867540198,
 -0.0091443848250055004725,
  0.0009010457337825074652,
 -0.0001394184127221491507,
  0.0000273815815785209370,
 -0.0000062750421119959424,
  0.0000016064844184831521,
 -0.0000004476392158510354,
  0.0000001334635874651668,
 -0.0000000420735334263215,
  0.0000000139021990246364,
 -0.0000000047831848068048,
  0.0000000017047897907465,
 -0.0000000006268389576018,
  0.0000000002369824276612,
 -0.0000000000918641139267,
  0.0000000000364278543037,
 -0.0000000000147475551725,
  0.0000000000060851006556,
 -0.0000000000025552772234,
  0.0000000000010906187250,
 -0.0000000000004725870319,
  0.0000000000002076969064,
 -0.0000000000000924976214,
  0.0000000000000417096723,
 -0.0000000000000190299093,
  0.0000000000000087790676,
 -0.0000000000000040927557,
  0.0000000000000019271068,
 -0.0000000000000009160199,
  0.0000000000000004393567,
 -0.0000000000000002125503,
  0.0000000000000001036735,
 -0.0000000000000000509642,
  0.0000000000000000252377,
 -0.0000000000000000125793
/*
  -.0187519297793868
  -.0091443848250055,
   .0009010457337825,
  -.0001394184127221,
   .0000273815815785,
  -.0000062750421119,
   .0000016064844184,
  -.0000004476392158,
   .0000001334635874,
  -.0000000420735334,
   .0000000139021990,
  -.0000000047831848,
   .0000000017047897,
  -.0000000006268389,
   .0000000002369824,
  -.0000000000918641,
   .0000000000364278,
  -.0000000000147475,
   .0000000000060851,
  -.0000000000025552,
   .0000000000010906,
  -.0000000000004725,
   .0000000000002076,
  -.0000000000000924,
   .0000000000000417,
  -.0000000000000190,
   .0000000000000087,
  -.0000000000000040,
   .0000000000000019,
  -.0000000000000009,
   .0000000000000004,
  -.0000000000000002,
   .0000000000000001,
  -.0000000000000000
*/
};
static cheb_series aip_cs = {
  data_aip,
  35,
  -1, 1,
  17
};

/* chebyshev for Bi(x) asymptotic factor
   based on SLATEC bie()

 Series for BIP        on the interval  1.25000D-01 to  3.53553D-01
                   with weighted error   1.91E-17
                    log weighted error  16.72
          significant figures required  15.35
               decimal places required  17.41

 Series for BIP2       on the interval  0.          to  1.25000D-01
                   with weighted error   1.05E-18
                    log weighted error  17.98
          significant figures required  16.74
               decimal places required  18.71
*/
static double data_bip[24] = {
  -0.08322047477943447,
   0.01146118927371174,
   0.00042896440718911,
  -0.00014906639379950,
  -0.00001307659726787,
   0.00000632759839610,
  -0.00000042226696982,
  -0.00000019147186298,
   0.00000006453106284,
  -0.00000000784485467,
  -0.00000000096077216,
   0.00000000070004713,
  -0.00000000017731789,
   0.00000000002272089,
   0.00000000000165404,
  -0.00000000000185171,
   0.00000000000059576,
  -0.00000000000012194,
   0.00000000000001334,
   0.00000000000000172,
  -0.00000000000000145,
   0.00000000000000049,
  -0.00000000000000011,
   0.00000000000000001
};
static cheb_series bip_cs = {
  data_bip,
  23,
  -1, 1,
  14
};

static double data_bip2[29] = {
  -0.113596737585988679,
   0.0041381473947881595,
   0.0001353470622119332,
   0.0000104273166530153,
   0.0000013474954767849,
   0.0000001696537405438,
  -0.0000000100965008656,
  -0.0000000167291194937,
  -0.0000000045815364485,
   0.0000000003736681366,
   0.0000000005766930320,
   0.0000000000621812650,
  -0.0000000000632941202,
  -0.0000000000149150479,
   0.0000000000078896213,
   0.0000000000024960513,
  -0.0000000000012130075,
  -0.0000000000003740493,
   0.0000000000002237727,
   0.0000000000000474902,
  -0.0000000000000452616,
  -0.0000000000000030172,
   0.0000000000000091058,
  -0.0000000000000009814,
  -0.0000000000000016429,
   0.0000000000000005533,
   0.0000000000000002175,
  -0.0000000000000001737,
  -0.0000000000000000010
};
static cheb_series bip2_cs = {
  data_bip2,
  28,
  -1, 1,
  10
};

/* assumes x >= 1.0 */
inline static int
airy_aie(const double x, gsl_mode_t mode, gsl_sf_result * result)
{
  double sqx = sqrt(x);
  double z   = 2.0/(x*sqx) - 1.0;
  double y   = sqrt(sqx);
  gsl_sf_result result_c;
  cheb_eval_mode_e(&aip_cs, z, mode, &result_c);
  result->val = (0.28125 + result_c.val)/y;
  result->err = result_c.err/y + GSL_DBL_EPSILON * fabs(result->val);
  return GSL_SUCCESS;
}

/* assumes x >= 2.0 */
static int airy_bie(const double x, gsl_mode_t mode, gsl_sf_result * result)
{
  const double ATR =  8.7506905708484345;
  const double BTR = -2.0938363213560543;

  if(x < 4.0) {
    double sqx = sqrt(x);
    double z   = ATR/(x*sqx) + BTR;
    double y   = sqrt(sqx);
    gsl_sf_result result_c;
    cheb_eval_mode_e(&bip_cs, z, mode, &result_c);
    result->val = (0.625 + result_c.val)/y;
    result->err = result_c.err/y + GSL_DBL_EPSILON * fabs(result->val);
  }
  else {
    double sqx = sqrt(x);
    double z   = 16.0/(x*sqx) - 1.0;
    double y   = sqrt(sqx);
    gsl_sf_result result_c;
    cheb_eval_mode_e(&bip2_cs, z, mode, &result_c);
    result->val = (0.625 + result_c.val)/y;
    result->err = result_c.err/y + GSL_DBL_EPSILON * fabs(result->val);
  }

  return GSL_SUCCESS;
}

/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/

int
gsl_sf_airy_Ai_e(const double x, const gsl_mode_t mode, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x < -1.0) {
    gsl_sf_result mod;
    gsl_sf_result theta;
    gsl_sf_result cos_result;
    int stat_mp  = airy_mod_phase(x, mode, &mod, &theta);
    int stat_cos = gsl_sf_cos_err_e(theta.val, theta.err, &cos_result);
    result->val  = mod.val * cos_result.val;
    result->err  = fabs(mod.val * cos_result.err) + fabs(cos_result.val * mod.err);
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    return GSL_ERROR_SELECT_2(stat_mp, stat_cos);
  }
  else if(x <= 1.0) {
    const double z = x*x*x;
    gsl_sf_result result_c0;
    gsl_sf_result result_c1;
    cheb_eval_mode_e(&aif_cs, z, mode, &result_c0);
    cheb_eval_mode_e(&aig_cs, z, mode, &result_c1);
    result->val  = 0.375 + (result_c0.val - x*(0.25 + result_c1.val));
    result->err  = result_c0.err + fabs(x*result_c1.err);
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    double x32 = x * sqrt(x);
    double s   = exp(-2.0*x32/3.0);
    gsl_sf_result result_aie;
    int stat_aie = airy_aie(x, mode, &result_aie);
    result->val  = result_aie.val * s;
    result->err  = result_aie.err * s + result->val * x32 * GSL_DBL_EPSILON;
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    CHECK_UNDERFLOW(result);
    return stat_aie;
  }
}

int
gsl_sf_airy_Ai_scaled_e(const double x, gsl_mode_t mode, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x < -1.0) {
    gsl_sf_result mod;
    gsl_sf_result theta;
    gsl_sf_result cos_result;
    int stat_mp  = airy_mod_phase(x, mode, &mod, &theta);
    int stat_cos = gsl_sf_cos_err_e(theta.val, theta.err, &cos_result);
    result->val  = mod.val * cos_result.val;
    result->err  = fabs(mod.val * cos_result.err) + fabs(cos_result.val * mod.err);
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    return GSL_ERROR_SELECT_2(stat_mp, stat_cos);
  }
  else if(x <= 1.0) {
    const double z = x*x*x;
    gsl_sf_result result_c0;
    gsl_sf_result result_c1;
    cheb_eval_mode_e(&aif_cs, z, mode, &result_c0);
    cheb_eval_mode_e(&aig_cs, z, mode, &result_c1);
    result->val  = 0.375 + (result_c0.val - x*(0.25 + result_c1.val));
    result->err  = result_c0.err + fabs(x*result_c1.err);
    result->err += GSL_DBL_EPSILON * fabs(result->val);

    if(x > 0.0) {
      const double scale = exp(2.0/3.0 * sqrt(z));
      result->val *= scale;
      result->err *= scale;
    }

    return GSL_SUCCESS;
  }
  else {
    return airy_aie(x, mode, result);
  }
}

int gsl_sf_airy_Bi_e(const double x, gsl_mode_t mode, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */
  if(x < -1.0) {
    gsl_sf_result mod;
    gsl_sf_result theta;
    gsl_sf_result sin_result;
    int stat_mp  = airy_mod_phase(x, mode, &mod, &theta);
    int stat_sin = gsl_sf_sin_err_e(theta.val, theta.err, &sin_result);
    result->val  = mod.val * sin_result.val;
    result->err  = fabs(mod.val * sin_result.err) + fabs(sin_result.val * mod.err);
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    return GSL_ERROR_SELECT_2(stat_mp, stat_sin);
  }
  else if(x < 1.0) {
    const double z = x*x*x;
    gsl_sf_result result_c0;
    gsl_sf_result result_c1;
    cheb_eval_mode_e(&bif_cs, z, mode, &result_c0);
    cheb_eval_mode_e(&big_cs, z, mode, &result_c1);
    result->val  = 0.625 + result_c0.val + x*(0.4375 + result_c1.val);
    result->err  = result_c0.err + fabs(x * result_c1.err);
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x <= 2.0) {
    const double z = (2.0*x*x*x - 9.0)/7.0;
    gsl_sf_result result_c0;
    gsl_sf_result result_c1;
    cheb_eval_mode_e(&bif2_cs, z, mode, &result_c0);
    cheb_eval_mode_e(&big2_cs, z, mode, &result_c1);
    result->val  = 1.125 + result_c0.val + x*(0.625 + result_c1.val);
    result->err  = result_c0.err + fabs(x * result_c1.err);
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    const double y = 2.0*x*sqrt(x)/3.0;
    const double s = exp(y);

    if(y > GSL_LOG_DBL_MAX - 1.0) {
      OVERFLOW_ERROR(result);
    }
    else {
      gsl_sf_result result_bie;
      int stat_bie = airy_bie(x, mode, &result_bie);
      result->val  = result_bie.val * s;
      result->err  = result_bie.err * s + fabs(1.5*y * (GSL_DBL_EPSILON * result->val));
      result->err += GSL_DBL_EPSILON * fabs(result->val);
      return stat_bie;
    }
  }
}

int
gsl_sf_airy_Bi_scaled_e(const double x, gsl_mode_t mode, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x < -1.0) {
    gsl_sf_result mod;
    gsl_sf_result theta;
    gsl_sf_result sin_result;
    int stat_mp  = airy_mod_phase(x, mode, &mod, &theta);
    int stat_sin = gsl_sf_sin_err_e(theta.val, theta.err, &sin_result);
    result->val  = mod.val * sin_result.val;
    result->err  = fabs(mod.val * sin_result.err) + fabs(sin_result.val * mod.err);
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    return GSL_ERROR_SELECT_2(stat_mp, stat_sin);
  }
  else if(x < 1.0) {
    const double z = x*x*x;
    gsl_sf_result result_c0;
    gsl_sf_result result_c1;
    cheb_eval_mode_e(&bif_cs, z, mode, &result_c0);
    cheb_eval_mode_e(&big_cs, z, mode, &result_c1);
    result->val  = 0.625 + result_c0.val + x*(0.4375 + result_c1.val);
    result->err  = result_c0.err + fabs(x * result_c1.err);
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    if(x > 0.0) {
      const double scale = exp(-2.0/3.0 * sqrt(z));
      result->val *= scale;
      result->err *= scale;
    }
    return GSL_SUCCESS;
  }
  else if(x <= 2.0) {
    const double x3 = x*x*x;
    const double z  = (2.0*x3 - 9.0)/7.0;
    const double s  = exp(-2.0/3.0 * sqrt(x3));
    gsl_sf_result result_c0;
    gsl_sf_result result_c1;
    cheb_eval_mode_e(&bif2_cs, z, mode, &result_c0);
    cheb_eval_mode_e(&big2_cs, z, mode, &result_c1);
    result->val  = s * (1.125 + result_c0.val + x*(0.625 + result_c1.val));
    result->err  = s * (result_c0.err + fabs(x * result_c1.err));
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    return airy_bie(x, mode, result);
  }
}

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

#include "eval.h"

double gsl_sf_airy_Ai(const double x, gsl_mode_t mode)
{
  EVAL_RESULT(gsl_sf_airy_Ai_e(x, mode, &result));
}

double gsl_sf_airy_Ai_scaled(const double x, gsl_mode_t mode)
{
  EVAL_RESULT(gsl_sf_airy_Ai_scaled_e(x, mode, &result));
}

double gsl_sf_airy_Bi(const double x, gsl_mode_t mode)
{
  EVAL_RESULT(gsl_sf_airy_Bi_e(x, mode, &result));
}

double gsl_sf_airy_Bi_scaled(const double x, gsl_mode_t mode)
{
  EVAL_RESULT(gsl_sf_airy_Bi_scaled_e(x, mode, &result));
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/* vectorviews.c */

// 8.4.13 Example programs for matrices of gsl-ref.pdf
// GNU GSL GNU Scientific Library reference book
// pp. 99,100

/*
  Compiling advice:
  gcc vectorviews.c -lgsl -lgslcblas -lm

 */

#include <math.h>
#include <stdio.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>

int main (void)
{
  size_t i,j;

  gsl_matrix *m = gsl_matrix_alloc(10,10);

  for (i = 0; i < 10; i++)
    for (j = 0; j < 10; j++)
      gsl_matrix_set (m, i, j, sin(i) + cos(j));

  for (j = 0; j < 10; j++)
    {
      gsl_vector_view column = gsl_matrix_column(m, j);
      double d;

      d = gsl_blas_dnrm2 (&column.vector);

      printf("matrix column %zu, norm = %g\n", j, d);
    }

  gsl_matrix_free(m);

  return 0;
}

ext: c

max_stars_repo_name: ernestyalumni/CompPhys

max_issues_repo_name: ernestyalumni/CompPhys

max_forks_repo_name: ernestyalumni/CompPhys


text:
#include <math.h>
#include <gsl/gsl_math.h>
#include "egsl_macros.h"

int main() {
	egsl_push();

		/* Creates a vector from a double array */
		double p[2] = {1,2};
		val vp = egsl_vFa(2,p);

		/* Creates a matrix from a double array */
		double md[4] = {
			1, 2,
			0, -1
		};
		val m = egsl_vFda(2,2,md);

		/* Creates a rotation matrix */
		val R = rot(M_PI/2);

		/* Multiplies the three together */
		val vrot = m3(R, m, vp);

		/* Displays the results */
		egsl_print("R", R);
		egsl_print("vp", vp);
		egsl_print("vrot", vrot);

		/* Create a semidifinite matrix (symmetric part of R) */
		val A = sc(0.5, sum(m, tr(m)) );

 		/* Displays spectrum */
		egsl_print("A",A);
		egsl_print_spectrum("A",A);

	egsl_pop();

	return 0;
}


text:
#ifndef OPENMC_MATERIAL_H
#define OPENMC_MATERIAL_H

#include <memory> // for unique_ptr
#include <string>
#include <unordered_map>
#include <vector>

#include <gsl/gsl>
#include <hdf5.h>
#include "pugixml.hpp"
#include "xtensor/xtensor.hpp"

#include "openmc/constants.h"
#include "openmc/bremsstrahlung.h"
#include "openmc/particle.h"

namespace openmc {

//==============================================================================
// Global variables
//==============================================================================

class Material;

namespace model {

extern std::unordered_map<int32_t, int32_t> material_map;
extern std::vector<std::unique_ptr<Material>> materials;

} // namespace model

//==============================================================================
//! A substance with constituent nuclides and thermal scattering data
//==============================================================================

class Material
{
public:
  //----------------------------------------------------------------------------
  // Types
  struct ThermalTable {
    int index_table; //!< Index of table in data::thermal_scatt
    int index_nuclide; //!< Index in nuclide_
    double fraction; //!< How often to use table
  };

  //----------------------------------------------------------------------------
  // Constructors, destructors, factory functions
  Material() {};
  explicit Material(pugi::xml_node material_node);
  ~Material();

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

  void calculate_xs(Particle& p) const;

  //! Assign thermal scattering tables to specific nuclides within the material
  //! so the code knows when to apply bound thermal scattering data
  void init_thermal();

  //! Set up mapping between global nuclides vector and indices in nuclide_
  void init_nuclide_index();

  //! Finalize the material, assigning tables, normalize density, etc.
  void finalize();

  //! Write material data to HDF5
  void to_hdf5(hid_t group) const;

  //! Add nuclide to the material
  //
  //! \param[in] nuclide Name of the nuclide
  //! \param[in] density Density of the nuclide in [atom/b-cm]
  void add_nuclide(const std::string& nuclide, double density);

  //! Set atom densities for the material
  //
  //! \param[in] name Name of each nuclide
  //! \param[in] density Density of each nuclide in [atom/b-cm]
  void set_densities(const std::vector<std::string>& name,
    const std::vector<double>& density);

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

  //! Get density in [atom/b-cm]
  //! \return Density in [atom/b-cm]
  double density() const { return density_; }

  //! Get density in [g/cm^3]
  //! \return Density in [g/cm^3]
  double density_gpcc() const { return density_gpcc_; }

  //! Get name
  //! \return Material name
  const std::string& name() const { return name_; }

  //! Set name
  void set_name(const std::string& name) { name_ = name; }

  //! Set total density of the material
  //
  //! \param[in] density Density value
  //! \param[in] units Units of density
  void set_density(double density, gsl::cstring_span units);

  //! Get nuclides in material
  //! \return Indices into the global nuclides vector
  gsl::span<const int> nuclides() const { return {nuclide_.data(), nuclide_.size()}; }

  //! Get densities of each nuclide in material
  //! \return Densities in [atom/b-cm]
  gsl::span<const double> densities() const { return {atom_density_.data(), atom_density_.size()}; }

  //! Get ID of material
  //! \return ID of material
  int32_t id() const { return id_; }

  //! Assign a unique ID to the material
  //! \param[in] Unique ID to assign. A value of -1 indicates that an ID
  //!   should be automatically assigned.
  void set_id(int32_t id);

  //! Get whether material is fissionable
  //! \return Whether material is fissionable
  bool fissionable() const { return fissionable_; }

  //! Get volume of material
  //! \return Volume in [cm^3]
  double volume() const;

  //! Get temperature of material
  //! \return Temperature in [K]
  double temperature() const;

  //----------------------------------------------------------------------------
  // Data
  int32_t id_ {C_NONE}; //!< Unique ID
  std::string name_; //!< Name of material
  std::vector<int> nuclide_; //!< Indices in nuclides vector
  std::vector<int> element_; //!< Indices in elements vector
  xt::xtensor<double, 1> atom_density_; //!< Nuclide atom density in [atom/b-cm]
  double density_; //!< Total atom density in [atom/b-cm]
  double density_gpcc_; //!< Total atom density in [g/cm^3]
  double volume_ {-1.0}; //!< Volume in [cm^3]
  bool fissionable_ {false}; //!< Does this material contain fissionable nuclides
  bool depletable_ {false}; //!< Is the material depletable?
  std::vector<bool> p0_; //!< Indicate which nuclides are to be treated with iso-in-lab scattering

  // To improve performance of tallying, we store an array (direct address
  // table) that indicates for each nuclide in data::nuclides the index of the
  // corresponding nuclide in the nuclide_ vector. If it is not present in the
  // material, the entry is set to -1.
  std::vector<int> mat_nuclide_index_;

  // Thermal scattering tables
  std::vector<ThermalTable> thermal_tables_;

  std::unique_ptr<Bremsstrahlung> ttb_;

private:
  //----------------------------------------------------------------------------
  // Private methods

  //! Calculate the collision stopping power
  void collision_stopping_power(double* s_col, bool positron);

  //! Initialize bremsstrahlung data
  void init_bremsstrahlung();

  //! Normalize density
  void normalize_density();

  void calculate_neutron_xs(Particle& p) const;
  void calculate_photon_xs(Particle& p) const;

  //----------------------------------------------------------------------------
  // Private data members
  gsl::index index_;

  //! \brief Default temperature for cells containing this material.
  //!
  //! A negative value indicates no default temperature was specified.
  double temperature_ {-1};
};

//==============================================================================
// Non-member functions
//==============================================================================

//! Calculate Sternheimer adjustment factor
double sternheimer_adjustment(const std::vector<double>& f, const
  std::vector<double>& e_b_sq, double e_p_sq, double n_conduction, double
  log_I, double tol, int max_iter);

//! Calculate density effect correction
double density_effect(const std::vector<double>& f, const std::vector<double>&
  e_b_sq, double e_p_sq, double n_conduction, double rho, double E, double tol,
  int max_iter);

//! Read material data from materials.xml
void read_materials_xml();

void free_memory_material();

} // namespace openmc
#endif // OPENMC_MATERIAL_H

ext: h

max_stars_repo_name: Hit-Weixg/openmc

max_issues_repo_name: Hit-Weixg/openmc

max_forks_repo_name: Hit-Weixg/openmc


text:
// Copyright (c) Microsoft Corporation. All rights reserved.
// Licensed under the MIT license.

#pragma once

#include <type_traits>
#include <algorithm>
#include <limits>
#include <iostream>
#include <cstring>
#include "seal/memorymanager.h"
#include "seal/serialization.h"
#include "seal/util/pointer.h"
#include "seal/util/defines.h"
#include "seal/util/common.h"
#ifdef SEAL_USE_MSGSL_SPAN
#include <gsl/span>
#endif

namespace seal
{
    class Ciphertext;

    // Forward-declaring the deflate_size_bound function
    namespace util::ztools
    {
        SEAL_NODISCARD std::size_t deflate_size_bound(std::size_t in_size) noexcept;
    }

    /**
    A resizable container for storing an array of arithmetic data types or
    SEAL_BYTE types (e.g., std::byte). The allocations are done from a memory
    pool. The IntArray class is mainly intended for internal use and provides
    the underlying data structure for Plaintext and Ciphertext classes.

    @par Size and Capacity
    IntArray allows the user to pre-allocate memory (capacity) for the array
    in cases where the array is known to be resized in the future and memory
    moves are to be avoided at the time of resizing. The size of the IntArray
    can never exceed its capacity. The capacity and size can be changed using
    the reserve and resize functions, respectively.

    @par Thread Safety
    In general, reading from IntArray is thread-safe as long as no other thread
    is concurrently mutating it.
    */
    template<typename T_,
        typename = std::enable_if_t<std::is_arithmetic<T_>::value ||
            std::is_same<typename std::decay<T_>::type, SEAL_BYTE>::value>>
    class IntArray
    {
        friend class Ciphertext;

    public:
        using T = typename std::decay<T_>::type;

        /**
        Creates a new IntArray. No memory is allocated by this constructor.

        @param[in] pool The MemoryPoolHandle pointing to a valid memory pool
        @throws std::invalid_argument if pool is uninitialized
        */
        IntArray(MemoryPoolHandle pool = MemoryManager::GetPool()) :
            pool_(std::move(pool))
        {
            if (!pool_)
            {
                throw std::invalid_argument("pool is uninitialized");
            }
        }

        /**
        Creates a new IntArray with given size.

        @param[in] size The size of the array
        @param[in] pool The MemoryPoolHandle pointing to a valid memory pool
        @throws std::invalid_argument if pool is uninitialized
        */
        explicit IntArray(std::size_t size,
            MemoryPoolHandle pool = MemoryManager::GetPool()) :
            pool_(std::move(pool))
        {
            if (!pool_)
            {
                throw std::invalid_argument("pool is uninitialized");
            }

            // Reserve memory, resize, and set to zero
            resize(size);
        }

        /**
        Creates a new IntArray with given capacity and size.

        @param[in] capacity The capacity of the array
        @param[in] size The size of the array
        @param[in] pool The MemoryPoolHandle pointing to a valid memory pool
        @throws std::invalid_argument if capacity is less than size
        @throws std::invalid_argument if pool is uninitialized
        */
        explicit IntArray(std::size_t capacity, std::size_t size,
            MemoryPoolHandle pool = MemoryManager::GetPool()) :
            pool_(std::move(pool))
        {
            if (!pool_)
            {
                throw std::invalid_argument("pool is uninitialized");
            }
            if (capacity < size)
            {
                throw std::invalid_argument("capacity cannot be smaller than size");
            }

            // Reserve memory, resize, and set to zero
            reserve(capacity);
            resize(size);
        }

        /**
        Creates a new IntArray with given size wrapping a given pointer. This
        constructor allocates no memory. If the IntArray goes out of scope, the
        Pointer object given here is destroyed. On resizing the IntArray to larger
        size, the data will be copied over to a new allocation from the memory pool
        pointer to by the given MemoryPoolHandle and the Pointer object given here
        will subsequently be destroyed. Unlike the other constructors, this one
        exposes the option of not automatically zero-filling the allocated memory.

        @param[in] ptr An initial Pointer object to wrap
        @param[in] capacity The capacity of the array
        @param[in] size The size of the array
        @param[in] fill_zero If true, fills ptr with zeros
        @param[in] pool The MemoryPoolHandle pointing to a valid memory pool
        @throws std::invalid_argument if ptr is null and capacity is positive
        @throws std::invalid_argument if capacity is less than size
        @throws std::invalid_argument if pool is uninitialized
        */
        explicit IntArray(util::Pointer<T> &&ptr,
            std::size_t capacity, std::size_t size, bool fill_zero,
            MemoryPoolHandle pool = MemoryManager::GetPool()) :
            pool_(std::move(pool)),
            capacity_(capacity)
        {
            if (!ptr && capacity)
            {
                throw std::invalid_argument("ptr cannot be null");
            }
            if (!pool_)
            {
                throw std::invalid_argument("pool is uninitialized");
            }
            if (capacity < size)
            {
                throw std::invalid_argument("capacity cannot be smaller than size");
            }

            // Grab the given Pointer
            data_ = std::move(ptr);

            // Resize, and optionally set to zero
            resize(size, fill_zero);
        }

        /**
        Creates a new IntArray with given size wrapping a given pointer. This
        constructor allocates no memory. If the IntArray goes out of scope, the
        Pointer object given here is destroyed. On resizing the IntArray to larger
        size, the data will be copied over to a new allocation from the memory pool
        pointer to by the given MemoryPoolHandle and the Pointer object given here
        will subsequently be destroyed. Unlike the other constructors, this one
        exposes the option of not automatically zero-filling the allocated memory.

        @param[in] ptr An initial Pointer object to wrap
        @param[in] size The size of the array
        @param[in] fill_zero If true, fills ptr with zeros
        @param[in] pool The MemoryPoolHandle pointing to a valid memory pool
        @throws std::invalid_argument if ptr is null and size is positive
        @throws std::invalid_argument if pool is uninitialized
        */
        explicit IntArray(util::Pointer<T> &&ptr, std::size_t size, bool fill_zero,
            MemoryPoolHandle pool = MemoryManager::GetPool()) :
            IntArray(std::move(ptr), size, size, fill_zero, std::move(pool))
        {
        }

        /**
        Constructs a new IntArray by copying a given one.

        @param[in] copy The IntArray to copy from
        */
        IntArray(const IntArray<T> &copy) :
            pool_(MemoryManager::GetPool()),
            capacity_(copy.size_),
            size_(copy.size_),
            data_(util::allocate<T>(copy.size_, pool_))
        {
            // Copy over value
            std::copy_n(copy.cbegin(), copy.size_, begin());
        }

        /**
        Constructs a new IntArray by moving a given one.

        @param[in] source The IntArray to move from
        */
        IntArray(IntArray<T> &&source) noexcept :
            pool_(std::move(source.pool_)),
            capacity_(source.capacity_),
            size_(source.size_),
            data_(std::move(source.data_))
        {
        }

        /**
        Destroys the IntArray.
        */
        ~IntArray()
        {
            release();
        }

        /**
        Returns a pointer to the beginning of the array data.
        */
        SEAL_NODISCARD inline T *begin() noexcept
        {
            return data_.get();
        }

        /**
        Returns a constant pointer to the beginning of the array data.
        */
        SEAL_NODISCARD inline const T *cbegin() const noexcept
        {
            return data_.get();
        }

        /**
        Returns a pointer to the end of the array data.
        */
        SEAL_NODISCARD inline T *end() noexcept
        {
            return size_ ? begin() + size_ : begin();
        }

        /**
        Returns a constant pointer to the end of the array data.
        */
        SEAL_NODISCARD inline const T *cend() const noexcept
        {
            return size_ ? cbegin() + size_ : cbegin();
        }
#ifdef SEAL_USE_MSGSL_SPAN
        /**
        Returns a span pointing to the beginning of the IntArray.
        */
        SEAL_NODISCARD inline gsl::span<T> span()
        {
            return gsl::span<T>(
                begin(), static_cast<std::ptrdiff_t>(size_));
        }

        /**
        Returns a span pointing to the beginning of the IntArray.
        */
        SEAL_NODISCARD inline gsl::span<const T> span() const
        {
            return gsl::span<const T>(
                cbegin(), static_cast<std::ptrdiff_t>(size_));
        }
#endif
        /**
        Returns a constant reference to the array element at a given index.
        This function performs bounds checking and will throw an error if
        the index is out of range.

        @param[in] index The index of the array element
        @throws std::out_of_range if index is out of range
        */
        SEAL_NODISCARD inline const T &at(std::size_t index) const
        {
            if (index >= size_)
            {
                throw std::out_of_range("index must be within [0, size)");
            }
            return data_[index];
        }

        /**
        Returns a reference to the array element at a given index. This
        function performs bounds checking and will throw an error if the
        index is out of range.

        @param[in] index The index of the array element
        @throws std::out_of_range if index is out of range
        */
        SEAL_NODISCARD inline T &at(std::size_t index)
        {
            if (index >= size_)
            {
                throw std::out_of_range("index must be within [0, size)");
            }
            return data_[index];
        }

        /**
        Returns a constant reference to the array element at a given index.
        This function does not perform bounds checking.

        @param[in] index The index of the array element
        */
        SEAL_NODISCARD inline const T &operator [](std::size_t index) const
        {
            return data_[index];
        }

        /**
        Returns a reference to the array element at a given index. This
        function does not perform bounds checking.

        @param[in] index The index of the array element
        */
        SEAL_NODISCARD inline T &operator [](std::size_t index)
        {
            return data_[index];
        }

        /**
        Returns whether the array has size zero.
        */
        SEAL_NODISCARD inline bool empty() const noexcept
        {
            return (size_ == 0);
        }

        /**
        Returns the largest possible array size.
        */
        SEAL_NODISCARD inline std::size_t max_size() const noexcept
        {
            return std::numeric_limits<std::size_t>::max();
        }

        /**
        Returns the size of the array.
        */
        SEAL_NODISCARD inline std::size_t size() const noexcept
        {
            return size_;
        }

        /**
        Returns the capacity of the array.
        */
        SEAL_NODISCARD inline std::size_t capacity() const noexcept
        {
            return capacity_;
        }

        /**
        Returns the currently used MemoryPoolHandle.
        */
        SEAL_NODISCARD inline MemoryPoolHandle pool() const noexcept
        {
            return pool_;
        }

        /**
        Releases any allocated memory to the memory pool and sets the size
        and capacity of the array to zero.
        */
        inline void release() noexcept
        {
            capacity_ = 0;
            size_ = 0;
            data_.release();
        }

        /**
        Sets the size of the array to zero. The capacity is not changed.
        */
        inline void clear() noexcept
        {
            size_ = 0;
        }

        /**
        Allocates enough memory for storing a given number of elements without
        changing the size of the array. If the given capacity is smaller than
        the current size, the size is automatically set to equal the new capacity.

        @param[in] capacity The capacity of the array
        */
        inline void reserve(std::size_t capacity)
        {
            std::size_t copy_size = std::min(capacity, size_);

            // Create new allocation and copy over value
            auto new_data(util::allocate<T>(capacity, pool_));
            std::copy_n(cbegin(), copy_size, new_data.get());
            std::swap(data_, new_data);

            // Set the coeff_count and capacity
            capacity_ = capacity;
            size_ = copy_size;
        }

        /**
        Reallocates the array so that its capacity exactly matches its size.
        */
        inline void shrink_to_fit()
        {
            reserve(size_);
        }

        /**
        Resizes the array to given size. When resizing to larger size the data
        in the array remains unchanged and any new space is initialized to zero
        if fill_zero is set to true; when resizing to smaller size the last
        elements of the array are dropped. If the capacity is not already large
        enough to hold the new size, the array is also reallocated.

        @param[in] size The size of the array
        @param[in] fill_zero If true, fills expanded space with zeros
        */
        inline void resize(std::size_t size, bool fill_zero = true)
        {
            if (size <= capacity_)
            {
                // Are we changing size to bigger within current capacity?
                // If so, need to set top terms to zero
                if (size > size_ && fill_zero)
                {
                    std::fill(end(), begin() + size, T(0));
                }

                // Set the size
                size_ = size;

                return;
            }

            // At this point we know for sure that size_ <= capacity_ < size so need
            // to reallocate to bigger
            auto new_data(util::allocate<T>(size, pool_));
            std::copy_n(cbegin(), size_, new_data.get());
            if (fill_zero)
            {
                std::fill(new_data.get() + size_, new_data.get() + size, T(0));
            }
            std::swap(data_, new_data);

            // Set the coeff_count and capacity
            capacity_ = size;
            size_ = size;
        }

        /**
        Copies a given IntArray to the current one.

        @param[in] assign The IntArray to copy from
        */
        inline IntArray<T> &operator =(const IntArray<T> &assign)
        {
            // Check for self-assignment
            if (this == &assign)
            {
                return *this;
            }

            // First resize to correct size
            resize(assign.size_);

            // Size is guaranteed to be OK now so copy over
            std::copy_n(assign.cbegin(), assign.size_, begin());

            return *this;
        }

        /**
        Moves a given IntArray to the current one.

        @param[in] assign The IntArray to move from
        */
        IntArray<T> &operator =(IntArray<T> &&assign) noexcept
        {
            capacity_ = assign.capacity_;
            size_ = assign.size_;
            data_ = std::move(assign.data_);
            pool_ = std::move(assign.pool_);

            return *this;
        }

        /**
        Returns an upper bound on the size of the IntArray, as if it was written
        to an output stream.

        @param[in] compr_mode The compression mode
        @throws std::invalid_argument if the compression mode is not supported
        @throws std::logic_error if the size does not fit in the return type
        */
        SEAL_NODISCARD inline std::streamoff save_size(
            compr_mode_type compr_mode) const
        {
            std::size_t members_size = Serialization::ComprSizeEstimate(
                util::add_safe(
                    sizeof(std::uint64_t), // size_
                    util::mul_safe(size_, sizeof(T))), // data_
                compr_mode);

            return util::safe_cast<std::streamoff>(util::add_safe(
                sizeof(Serialization::SEALHeader),
                members_size
            ));
        }

        /**
        Saves the IntArray to an output stream. The output is in binary format
        and not human-readable. The output stream must have the "binary" flag set.

        @param[out] stream The stream to save the IntArray to
        @param[in] compr_mode The desired compression mode
        @throws std::logic_error if the data to be saved is invalid, if compression
        mode is not supported, or if compression failed
        @throws std::runtime_error if I/O operations failed
        */
        inline std::streamoff save(
            std::ostream &stream,
            compr_mode_type compr_mode = Serialization::compr_mode_default) const
        {
            using namespace std::placeholders;
            return Serialization::Save(
                std::bind(&IntArray<T_>::save_members, this, _1),
                save_size(compr_mode_type::none),
                stream, compr_mode);
        }

        /**
        Loads a IntArray from an input stream overwriting the current IntArray.
        This function takes optionally a bound on the size for the loaded IntArray
        and throws an exception if the size indicated by the loaded metadata exceeds
        the provided value. The check is omitted if in_size_bound is zero.

        @param[in] stream The stream to load the IntArray from
        @param[in] in_size_bound A bound on the size of the loaded IntArray
        @throws std::logic_error if the loaded data is invalid, if the loaded size
        exceeds in_size_bound, or if decompression failed
        @throws std::runtime_error if I/O operations failed
        */
        inline std::streamoff load(
            std::istream &stream, std::size_t in_size_bound = 0)
        {
            using namespace std::placeholders;
            return Serialization::Load(
                std::bind(&IntArray<T_>::load_members, this, _1, in_size_bound),
                stream);
        }

        /**
        Saves the IntArray to a given memory location. The output is in binary
        format and not human-readable.

        @param[out] out The memory location to write the SmallModulus to
        @param[in] size The number of bytes available in the given memory location
        @param[in] compr_mode The desired compression mode
        @throws std::invalid_argument if out is null or if size is too small to
        contain a SEALHeader
        @throws std::logic_error if the data to be saved is invalid, if compression
        mode is not supported, or if compression failed
        @throws std::runtime_error if I/O operations failed
        */
        inline std::streamoff save(
            SEAL_BYTE *out,
            std::size_t size,
            compr_mode_type compr_mode = Serialization::compr_mode_default) const
        {
            using namespace std::placeholders;
            return Serialization::Save(
                std::bind(&IntArray<T_>::save_members, this, _1),
                save_size(compr_mode_type::none),
                out, size, compr_mode);
        }

        /**
        Loads a IntArray from a given memory location overwriting the current
        IntArray. This function takes optionally a bound on the size for the loaded
        IntArray and throws an exception if the size indicated by the loaded
        metadata exceeds the provided value. The check is omitted if in_size_bound
        is zero.

        @param[in] in The memory location to load the SmallModulus from
        @param[in] size The number of bytes available in the given memory location
        @param[in] in_size_bound A bound on the size of the loaded IntArray
        @throws std::invalid_argument if in is null or if size is too small to
        contain a SEALHeader
        @throws std::logic_error if the loaded data is invalid, if the loaded size
        exceeds in_size_bound, or if decompression failed
        @throws std::runtime_error if I/O operations failed
        */
        inline std::streamoff load(
            const SEAL_BYTE *in, std::size_t size, std::size_t in_size_bound = 0)
        {
            using namespace std::placeholders;
            return Serialization::Load(
                std::bind(&IntArray<T_>::load_members, this, _1, in_size_bound),
                in, size);
        }

    private:
        void save_members(std::ostream &stream) const
        {
            auto old_except_mask = stream.exceptions();
            try
            {
                // Throw exceptions on std::ios_base::badbit and std::ios_base::failbit
                stream.exceptions(std::ios_base::badbit | std::ios_base::failbit);

                std::uint64_t size64 = size_;
                stream.write(reinterpret_cast<const char*>(&size64), sizeof(std::uint64_t));
                if (size_)
                {
                    stream.write(reinterpret_cast<const char*>(cbegin()),
                        util::safe_cast<std::streamsize>(util::mul_safe(size_, sizeof(T))));
                }
            }
            catch (const std::ios_base::failure &)
            {
                stream.exceptions(old_except_mask);
                throw std::runtime_error("I/O error");
            }
            catch (...)
            {
                stream.exceptions(old_except_mask);
                throw;
            }
            stream.exceptions(old_except_mask);
        }

        void load_members(std::istream &stream, std::size_t in_size_bound)
        {
            auto old_except_mask = stream.exceptions();
            try
            {
                // Throw exceptions on std::ios_base::badbit and std::ios_base::failbit
                stream.exceptions(std::ios_base::badbit | std::ios_base::failbit);

                std::uint64_t size64 = 0;
                stream.read(reinterpret_cast<char*>(&size64), sizeof(std::uint64_t));

                // Check (optionally) that the size in the metadata does not exceed
                // in_size_bound
                if (in_size_bound && util::unsigned_gt(size64, in_size_bound))
                {
                    throw std::logic_error("unexpected size");
                }

                // Set new size; this is potentially unsafe if size64 was not checked
                // against expected_size
                resize(util::safe_cast<std::size_t>(size64));

                // Read data
                if (size_)
                {
                    stream.read(reinterpret_cast<char*>(begin()),
                        util::safe_cast<std::streamsize>(util::mul_safe(size_, sizeof(T))));
                }
            }
            catch (const std::ios_base::failure &)
            {
                stream.exceptions(old_except_mask);
                throw std::runtime_error("I/O error");
            }
            catch (...)
            {
                stream.exceptions(old_except_mask);
                throw;
            }
            stream.exceptions(old_except_mask);
        }

        MemoryPoolHandle pool_;

        std::size_t capacity_ = 0;

        std::size_t size_ = 0;

        util::Pointer<T> data_;
    };
}

max_stars_repo_name: nitrieu/SEAL

max_issues_repo_name: nitrieu/SEAL

max_forks_repo_name: nitrieu/SEAL


text:
/*
 * test_suite.h
 *
 *  Created on: Oct 16, 2013
 *      Author: vbonnici
 */

#ifndef TEST_SUITE_H_
#define TEST_SUITE_H_

/*
 * a set of functions for binomial tests, or other distribution
 */

//#define DEBUG_TEST_SUITE_H_

#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <math.h>

#include <vector>
#include <limits>

#include "data_ts.h"

#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>

namespace dolierlib{

/*
 * Pietro's binomial test
 */
inline
double dolier_binomial(const usize_t f_dcount, const usize_t f_total, const usize_t b_dcount, const usize_t b_total){
	double  p = (double)b_dcount / (double)b_total;
	double result = 0;
	for(usize_t i = f_dcount; i<=f_total; i++){
		result += gsl_ran_binomial_pdf(static_cast<unsigned int>(i), p, static_cast<unsigned int>(f_total));
	}
	return result;
}

/*
 * a more efficient binomial test, thanks to gsl
 */
inline
double dolier_binomial_Q(const usize_t f_dcount, const usize_t f_total, const usize_t b_dcount, const usize_t b_total){
	double  p = (double)b_dcount / (double)b_total;
	//return gsl_cdf_binomial_Q(static_cast<unsigned int>(f_dcount), p, static_cast<unsigned int>(f_total));
	return gsl_ran_binomial_pdf(static_cast<unsigned int>(f_dcount), p, static_cast<unsigned int>(f_total))
			+ gsl_cdf_binomial_Q(static_cast<unsigned int>(f_dcount), p, static_cast<unsigned int>(f_total));
}

/*
 * Pietro's binomial test with Bon Ferroni correction
 */
inline
double dolier_binomial_bonf(const usize_t f_dcount, const usize_t f_total, const usize_t b_dcount, const usize_t b_total, const usize_t nof_kmers){
	double p = (double)b_dcount / (double)b_total;
	double result = 0;
	for(usize_t i = f_dcount; i<=f_total; i++){
		result += gsl_ran_binomial_pdf(static_cast<unsigned int>(i), p, static_cast<unsigned int>(f_total));
	}
	return double(nof_kmers) * result;
}

/*
 * a more efficient binomial test with Bon Ferroni correction, thanks to gsl
 */
inline
double dolier_binomial_bonf_Q(const usize_t f_dcount, const usize_t f_total, const usize_t b_dcount, const usize_t b_total, const usize_t nof_kmers){
	double  p = (double)b_dcount / (double)b_total;
	//return gsl_cdf_binomial_Q(static_cast<unsigned int>(f_dcount), p, static_cast<unsigned int>(f_total)) * nof_kmers;
	return (
			gsl_ran_binomial_pdf(static_cast<unsigned int>(f_dcount), p, static_cast<unsigned int>(f_total))
			+
			gsl_cdf_binomial_Q(static_cast<unsigned int>(f_dcount), p, static_cast<unsigned int>(f_total))
			)
			* nof_kmers;
}

/*
 * other distributions tried during tests
 */

inline
double dolier_poisson_Q(const usize_t f_dcount, const usize_t f_total, const usize_t b_dcount, const usize_t b_total){
	double  p = (double)b_dcount / (double)b_total;
	double mu = p *  static_cast<double>(f_total);
	return gsl_ran_poisson_pdf(static_cast<unsigned int>(f_dcount),mu)
			+ gsl_cdf_poisson_Q(static_cast<unsigned int>(f_dcount), mu);
//	return gsl_ran_binomial_pdf(static_cast<unsigned int>(f_dcount), p, static_cast<unsigned int>(f_total))
//			+ gsl_cdf_binomial_Q(static_cast<unsigned int>(f_dcount), p, static_cast<unsigned int>(f_total));
}

inline
double dolier_poisson_bonf_Q(const usize_t f_dcount, const usize_t f_total, const usize_t b_dcount, const usize_t b_total, const usize_t nof_kmers){
	double  p = (double)b_dcount / (double)b_total;
	double mu = p *  static_cast<double>(f_total);
	if(mu==0) std::cout<<"#mu=0#\t"<<f_dcount<<"\t"<<f_total<<"\t"<<b_dcount<<"\t"<<b_total<<"\t"<<p<<"\n";
	return 		(gsl_ran_poisson_pdf(static_cast<unsigned int>(f_dcount),mu)
				+ gsl_cdf_poisson_Q(static_cast<unsigned int>(f_dcount), mu))
			* nof_kmers;
//	return (
//			gsl_ran_binomial_pdf(static_cast<unsigned int>(f_dcount), p, static_cast<unsigned int>(f_total))
//			+
//			gsl_cdf_binomial_Q(static_cast<unsigned int>(f_dcount), p, static_cast<unsigned int>(f_total))
//			)
//			* nof_kmers;
}

inline
double dolier_poisson_mu1_Q(const usize_t f_dcount, const usize_t f_total, const usize_t b_dcount, const usize_t b_total){
	return gsl_ran_poisson_pdf(static_cast<unsigned int>(f_dcount),1)
			+ gsl_cdf_poisson_Q(static_cast<unsigned int>(f_dcount), 1);
}
inline
double dolier_poisson_mu1_bonf_Q(const usize_t f_dcount, const usize_t f_total, const usize_t b_dcount, const usize_t b_total, const usize_t nof_kmers){
	return 		(gsl_ran_poisson_pdf(static_cast<unsigned int>(f_dcount),1)
				+ gsl_cdf_poisson_Q(static_cast<unsigned int>(f_dcount), 1))
			* nof_kmers;
}

const double SQRT_PI  = sqrt(3.141592653589793238462);

inline
double dolier_normal(double x, double sigma, double mu){
	return (1/(sigma * SQRT_PI))*(exp(-(((x-mu)*(x-mu)) / (2*sigma*sigma))));
}

inline
double dolier_approx_binomial(const usize_t f_dcount, const usize_t f_total, const usize_t b_dcount, const usize_t b_total){
	double  p = (double)b_dcount / (double)b_total;
	double result = 0;

	double np = static_cast<double>(f_total) * p;
	double npp = np * (1-p);

	for(usize_t i = f_dcount; i<=f_total; i++){
		result += dolier_normal(static_cast<double>(i), np, npp);
	}
	return result;
}

inline
double dolier_binomial(unsigned int x, double p, unsigned int total){
	return gsl_ran_binomial_pdf(x, p, total);
}

inline
double dolier_normall(unsigned int x, double p, unsigned int total){
	double np = static_cast<double>(total) * p;
	double npp = np * (1-p);
	return dolier_normal(static_cast<double>(x), np, npp);
}

}

#endif /* TEST_SUITE_H_ */

max_stars_repo_name: vbonnici/DoLiER

max_issues_repo_name: vbonnici/DoLiER

max_forks_repo_name: vbonnici/DoLiER


text:
/* cdf/negbinom.c
 *
 * Copyright (C) 2004 Jason H. Stover.
 *
 * 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 <config.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_cdf.h>

#include "error.h"

/*
 * Pr(X <= k) for a negative binomial random variable X, i.e.,
 * the probability of k or fewer failuers before success n.
 */

double
gsl_cdf_negative_binomial_P (const unsigned int k, const double p, const double n)
{
  double P;
  double a;
  double b;

  if (p > 1.0 || p < 0.0)
    {
      CDF_ERROR ("p < 0 or p > 1", GSL_EDOM);
    }

  if (n < 0)
    {
      CDF_ERROR ("n < 0", GSL_EDOM);
    }

  a = (double) n;
  b = (double) k + 1.0;
  P = gsl_cdf_beta_P (p, a, b);

  return P;
}

/*
 * Pr ( X > k ).
 */

double
gsl_cdf_negative_binomial_Q (const unsigned int k, const double p, const double n)
{
  double Q;
  double a;
  double b;

  if (p > 1.0 || p < 0.0)
    {
      CDF_ERROR ("p < 0 or p > 1", GSL_EDOM);
    }

  if (n < 0)
    {
      CDF_ERROR ("n < 0", GSL_EDOM);
    }

  a = (double) n;
  b = (double) k + 1.0;
  Q = gsl_cdf_beta_Q (p, a, b);

  return Q;
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


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

#include <config.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_vector.h>

/* Return the location of the i-th Greville abscissa */
double
gsl_bspline_greville_abscissa (size_t i, gsl_bspline_workspace *w)
{
  const size_t stride = w->knots->stride;
  size_t km1 = w->km1;
  double * data = w->knots->data + (i+1)*stride;
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= gsl_bspline_ncoeffs(w)))
    {
      GSL_ERROR_VAL ("Greville abscissa index out of range", GSL_EINVAL, 0);
    }
#endif

  if (km1 == 0)
    {
      /* Return interval midpoints in degenerate k = 1 case*/
      km1   = 2;
      data -= stride;
    }

  return gsl_stats_mean(data, stride, km1);
}

int
gsl_bspline_knots_greville (const gsl_vector *abscissae,
                            gsl_bspline_workspace *w,
                            double *abserr)
{
  /* Limited function: see https://savannah.gnu.org/bugs/index.php?34361 */

  int s;

  /* Check incoming arguments satisfy mandatory algorithmic assumptions */
  if (w->k < 2)
    GSL_ERROR ("w->k must be at least 2", GSL_EINVAL);
  else if (abscissae->size < 2)
    GSL_ERROR ("abscissae->size must be at least 2", GSL_EINVAL);
  else if (w->nbreak != abscissae->size - w->k + 2)
    GSL_ERROR ("w->nbreak must equal abscissae->size - w->k + 2", GSL_EINVAL);

  if (w->nbreak == 2)
    {
      /* No flexibility in abscissae values possible in this degenerate case */
      s = gsl_bspline_knots_uniform (
              gsl_vector_get (abscissae, 0),
              gsl_vector_get (abscissae, abscissae->size - 1), w);
    }
  else
    {
      double * storage;
      gsl_matrix_view A;
      gsl_vector_view tau, b, x, r;
      size_t i, j;

      /* Constants derived from the B-spline workspace and abscissae details */
      const size_t km2    = w->k - 2;
      const size_t M      = abscissae->size - 2;
      const size_t N      = w->nbreak - 2;
      const double invkm1 = 1.0 / w->km1;

      /* Allocate working storage and prepare multiple, zero-filled views */
      storage = (double *) calloc (M*N + 2*N + 2*M, sizeof (double));
      if (storage == 0)
        GSL_ERROR ("failed to allocate working storage", GSL_ENOMEM);
      A   = gsl_matrix_view_array (storage, M, N);
      tau = gsl_vector_view_array (storage + M*N,             N);
      b   = gsl_vector_view_array (storage + M*N + N,         M);
      x   = gsl_vector_view_array (storage + M*N + N + M,     N);
      r   = gsl_vector_view_array (storage + M*N + N + M + N, M);

      /* Build matrix from interior breakpoints to interior Greville abscissae.
       * For example, when w->k = 4 and w->nbreak = 7 the matrix is
       *   [   1,      0,      0,      0,      0;
       *     2/3,    1/3,      0,      0,      0;
       *     1/3,    1/3,    1/3,      0,      0;
       *       0,    1/3,    1/3,    1/3,      0;
       *       0,      0,    1/3,    1/3,    1/3;
       *       0,      0,      0,    1/3,    2/3;
       *       0,      0,      0,      0,      1  ]
       * but only center formed as first/last breakpoint is known.
       */
      for (j = 0; j < N; ++j)
        for (i = 0; i <= km2; ++i)
          gsl_matrix_set (&A.matrix, i+j, j, invkm1);

      /* Copy interior collocation points from abscissae into b */
      for (i = 0; i < M; ++i)
        gsl_vector_set (&b.vector, i, gsl_vector_get (abscissae, i+1));

      /* Adjust b to account for constraint columns not stored in A */
      for (i = 0; i < km2; ++i)
        {
          double * const v = gsl_vector_ptr (&b.vector, i);
          *v -= (1 - (i+1)*invkm1) * gsl_vector_get (abscissae, 0);
        }
      for (i = 0; i < km2; ++i)
        {
          double * const v = gsl_vector_ptr (&b.vector, M - km2 + i);
          *v -= (i+1)*invkm1 * gsl_vector_get (abscissae, abscissae->size - 1);
        }

      /* Perform linear least squares to determine interior breakpoints */
      s =  gsl_linalg_QR_decomp (&A.matrix, &tau.vector)
        || gsl_linalg_QR_lssolve (&A.matrix, &tau.vector,
                                  &b.vector, &x.vector, &r.vector);
      if (s)
        {
          free (storage);
          return s;
        }

      /* "Expand" solution x by adding known first and last breakpoints. */
      x = gsl_vector_view_array_with_stride (
          gsl_vector_ptr (&x.vector, 0) - x.vector.stride,
          x.vector.stride, x.vector.size + 2);
      gsl_vector_set (&x.vector, 0, gsl_vector_get (abscissae, 0));
      gsl_vector_set (&x.vector, x.vector.size - 1,
                      gsl_vector_get (abscissae, abscissae->size - 1));

      /* Finally, initialize workspace knots using the now-known breakpoints */
      s = gsl_bspline_knots (&x.vector, w);
      free (storage);
    }

  /* Sum absolute errors in the resulting vs requested interior abscissae */
  /* Provided as a fit quality metric which may be monitored by callers */
  if (!s && abserr)
    {
      size_t i;
      *abserr = 0;
      for (i = 1; i < abscissae->size - 1; ++i)
        *abserr += fabs (   gsl_bspline_greville_abscissa (i, w)
                          - gsl_vector_get (abscissae, i) );
    }

  return s;
}


text:
//Include guard
#ifndef SEQUENCE_H
#define SEQUENCE_H

//Forward declared dependencies

//Included dependencies
#include <vector>
#include <string>
#include <gsl/gsl_rng.h>

class Sequence{
private:
    //Internal variables
    std::vector<char> seq;

public:
    Sequence();
    Sequence(int length);
    Sequence(std::vector<char>& s);
    Sequence(std::string s);
    ~Sequence();

    void setSeq(std::vector<char>& s);
	std::string getSeqAsString();
    void setBase(int index, char b);
    void addBase(char b);
	void removeBase(int index);
    char getBase(int index);
    int getLength();
    void print();
    std::string printToString();
    std::string printToStringNoHyphens();
    Sequence subset(int start, int finish);
    void addBaseFront(char b);
};

#endif


text:
static char help[] =
"ODE system solver example using TS.  Solves N-dimensional system\n"
"    dy/dt = G(t,y)\n"
"with y(t0) = y0 to compute y(tf).  Serial only.\n"
"Sets TS type to Runge-Kutta.  The implemented example has\n"
"G_0 = y_1, G_1 = - y_0 + t, y_0(0) = 0, y_1(0) = 0.  The exact solution is\n"
"y_0(t) = t - sin(t), y_1(t) = 1 - cos(t).\n\n";

#include <petsc.h>

extern PetscErrorCode ExactSolution(double, Vec);
extern PetscErrorCode FormRHSFunction(TS, double, Vec, Vec, void*);

//STARTMAIN
int main(int argc,char **argv) {
  PetscErrorCode ierr;
  int       steps;
  double    t0 = 0.0, tf = 20.0, dt = 0.1, err;
  Vec       y, yexact;
  TS        ts;

  PetscInitialize(&argc,&argv,(char*)0,help);

  ierr = VecCreate(PETSC_COMM_WORLD,&y); CHKERRQ(ierr);
  ierr = VecSetSizes(y,PETSC_DECIDE,2); CHKERRQ(ierr);
  ierr = VecSetFromOptions(y); CHKERRQ(ierr);
  ierr = VecDuplicate(y,&yexact); CHKERRQ(ierr);

  ierr = TSCreate(PETSC_COMM_WORLD,&ts); CHKERRQ(ierr);
  ierr = TSSetProblemType(ts,TS_NONLINEAR); CHKERRQ(ierr);
  ierr = TSSetRHSFunction(ts,NULL,FormRHSFunction,NULL); CHKERRQ(ierr);
  ierr = TSSetType(ts,TSRK); CHKERRQ(ierr);

  // set time axis
  ierr = TSSetTime(ts,t0); CHKERRQ(ierr);
  ierr = TSSetMaxTime(ts,tf); CHKERRQ(ierr);
  ierr = TSSetTimeStep(ts,dt); CHKERRQ(ierr);
  ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP); CHKERRQ(ierr);
  ierr = TSSetFromOptions(ts); CHKERRQ(ierr);

  // set initial values and solve
  ierr = TSGetTime(ts,&t0); CHKERRQ(ierr);
  ierr = ExactSolution(t0,y); CHKERRQ(ierr);
  ierr = TSSolve(ts,y); CHKERRQ(ierr);

  // compute error and report
  ierr = TSGetStepNumber(ts,&steps); CHKERRQ(ierr);
  ierr = TSGetTime(ts,&tf); CHKERRQ(ierr);
  ierr = ExactSolution(tf,yexact); CHKERRQ(ierr);
  ierr = VecAXPY(y,-1.0,yexact); CHKERRQ(ierr);    // y <- y - yexact
  ierr = VecNorm(y,NORM_INFINITY,&err); CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,
              "error at tf = %.3f with %d steps:  |y-y_exact|_inf = %g\n",
              tf,steps,err); CHKERRQ(ierr);

  VecDestroy(&y);  VecDestroy(&yexact);  TSDestroy(&ts);
  return PetscFinalize();
}
//ENDMAIN

//STARTCALLBACKS
PetscErrorCode ExactSolution(double t, Vec y) {
    double *ay;
    VecGetArray(y,&ay);
    ay[0] = t - sin(t);
    ay[1] = 1.0 - cos(t);
    VecRestoreArray(y,&ay);
    return 0;
}

PetscErrorCode FormRHSFunction(TS ts, double t, Vec y, Vec g, void *ptr) {
    const double *ay;
    double       *ag;
    VecGetArrayRead(y,&ay);
    VecGetArray(g,&ag);
    ag[0] = ay[1];            // = G_1(t,y)
    ag[1] = - ay[0] + t;      // = G_2(t,y)
    VecRestoreArrayRead(y,&ay);
    VecRestoreArray(g,&ag);
    return 0;
}
//ENDCALLBACKS