16 unsigned int num_spec) {
17 if (num_spec <= 0)
return 0;
19 time_deriv->production_rates =
20 (
long double *)malloc(num_spec *
sizeof(
long double));
21 if (time_deriv->production_rates == NULL)
return 0;
23 time_deriv->loss_rates =
24 (
long double *)malloc(num_spec *
sizeof(
long double));
25 if (time_deriv->loss_rates == NULL) {
26 free(time_deriv->production_rates);
30 time_deriv->num_spec = num_spec;
33 time_deriv->last_max_loss_precision = 0.0;
47 double *deriv_est,
unsigned int output_precision) {
48 long double *r_p = time_deriv.production_rates;
49 long double *r_l = time_deriv.loss_rates;
52 time_deriv.last_max_loss_precision = 1.0;
55 for (
unsigned int i_spec = 0; i_spec < time_deriv.num_spec; ++i_spec) {
56 double prec_loss = 1.0;
57 if (*r_p + *r_l != 0.0) {
59 long double scale_fact;
62 (1.0 / (*r_p + *r_l) + MAX_PRECISION_LOSS / fabsl(*r_p - *r_l));
64 scale_fact * (*r_p - *r_l) + (1.0 - scale_fact) * (*deriv_est);
66 *dest_array = *r_p - *r_l;
69 if (*r_p != 0.0 && *r_l != 0.0) {
70 prec_loss = *r_p > *r_l ? 1.0 - *r_l / *r_p : 1.0 - *r_p / *r_l;
71 if (prec_loss < time_deriv.last_max_loss_precision)
72 time_deriv.last_max_loss_precision = prec_loss;
81 if (deriv_est) ++deriv_est;
83 if (output_precision == 1) {
84 printf(
"\nspec %d prec_loss %le", i_spec, -log(prec_loss) / log(2.0));