回调函数有一个奇怪的行为，我想用二分法计算一个函数的根，要解决的系统涉及到一个函数的积分，其中上限是另一个函数的函数，这不是什么大不了的事情，当我把数值放进去的时候，程序结束得很好，但是如果我做函数的回调，终端杀死了进程，程序运行得更慢。顺便说一下，这些是我用c++写的第一段代码。
代码块是下一个;

#include <gsl/gsl_integration.h>
#include <functional>
#include <iostream>
#include <cmath>
struct Ninj_Params{
   double p;
   double Q0;
   double gamma_m;
};
double c = 1;
double eV = 1.0;
double MeV = pow(10,6) * eV;
double me = 0.510998902 * MeV;
double e = 8.542 * pow(10,-2);
double fm = 1/197.3269602 * 1/MeV;
double cm = pow(10,13) * fm;
double Gauss = 5.788381749 *  pow(10,-15) *  MeV * 2 * me/e;
double meter = pow(10,2) * cm;
double sec = 299792458*meter;
// Function to find the root using the bisection method
double bisection(const std::function<double(double)> &func, double a, double b, double tolerance){
   if (func(a) * func(b) >= 0){
       printf("Bisection method cannot guarantee convergence in this interval.\n");
       return NAN; // Not a Number to indicate failure
   }
   double c = a;
   while ((b - a) >= tolerance) {
       // Find middle point
       c = (a + b) / 2;
       // Check if the root is found
       if (func(c) == 0.0)
           break;
       // Update the interval based on the sign of the function at the middle point
       if (func(c) * func(a) < 0)
           b = c;
       else
           a = c;
   }
   return c;
}
double beta(double Gamma){
   
   return sqrt(1. - pow(Gamma, -2.));
}
double R(double t, double Rs, double Gamma){
   
   return (Rs + beta(Gamma)*c*Gamma*t)*1./cm;
}
double B(double t, double Rs, double R0, double Gamma, double B0, double a){
   
   return B0*pow((R(t, Rs, Gamma)*cm)/R0, -a)*1./Gauss;
   
}
double Q(double gammae, void * args){
   
   Ninj_Params& params = *(Ninj_Params*)args;
   
   if(params.gamma_m <= gammae){
       return params.Q0*pow(gammae/params.gamma_m, -params.p)*sec;
   }else{
       return 0.0;
   }
}
double Ninj(double Q0, double p, double gamma_min, double gamma_max){
   
   struct Ninj_Params params = {
       .p = 2.5,
       .Q0 = Q0,
       .gamma_m = gamma_min,
   };
   
   gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000);
   double result, error;
   gsl_function F;
   
   F.function = &Q;
   F.params   = &params;
   
   gsl_integration_qags (&F, gamma_min, gamma_max, 0, 1e-7, 1000, w, &result, &error);
   return result;
}
double Q0rate(double Ninj_Gest, double t, double Rs, double R0, double Gamma, double B0, double a){
   
   //double gamma_max = 1e8*pow(B(t, Rs, R0, Gamma, B0, a), -0.5);
   double gamma_max = 7.957e+06;
   printf("**** %.3e\n", gamma_max);
   double p = 2.5;
   
   auto root2solve = [Ninj_Gest, p, &gamma_max](double Q0) -> double{
       return Ninj_Gest - Ninj(Q0*1/sec, p, 1e5, gamma_max);
   };
   
   double root = bisection(root2solve, 1e40*1./sec, 1e53*1./sec, 1e-7);
   return root*sec;
}
void TEST_ROOT(double t, double Rs, double R0, double Gamma, double gamma_m, double B0, double a){
   
   double Nj = 1e47*1/sec;
   printf("%.5e\n",Q0rate(Nj, t, Rs, R0, Gamma, B0, a));
}
int main(){
   
   double t = 10*sec;
   double Rs = 1e14*cm;
   double R0 = 1e15*cm;
   double B0 = 0.3*1e2*Gauss;
   double Gamma = 300;
   double gamma_m = 1e5;
   double a = 1.0;
   TEST_ROOT(t, Rs, R0, Gamma, gamma_m, B0, a);
   
}

字符串
有问题的行是105（和注解的106）。When设置为：

double gamma_max = 1e8*pow(B(t, Rs, R0, Gamma, B0, a), -0.5);

型
输出是（具有缓慢的行为）：

• 7.957e+06 zsh：killed ./Main

如果更改为：

double gamma_max = 7.957e+06;

型
程序以期望的根值结束：
% ./主要 * 7.957e+06 1.50212e+42