Ternary Search

#include <bits/stdc++.h>
using namespace std ;

double M_a , M_b ,  V_a , V_eq , K_a , K_w , V_m , V_p ;

double pH(double V_b)
{
    double beta , b , c , x , E ;
    x = (V_eq+(K_a/M_b)*(V_a+V_eq)) ;

    beta = (K_w*K_w*(V_a+V_eq)*(V_a+V_eq))/(2*M_b*M_b*x*x) ;
    beta *= (1 + sqrt(1 + (4*K_a*M_b*x)/(K_w*(V_a+V_eq)))) ;
    beta += (K_a*K_w*(V_a+V_eq))/(M_b*x) ;
    beta = 1/sqrt(beta) ;

    x = V_b + (K_a*(V_a+V_b))/M_b ;
    c = -(K_a*K_w*(V_a+V_b))/(M_b*x) ;
    b = (K_a*(V_b - (M_a/M_b)*V_a)) - (K_w/M_b)*(V_a+V_b) ;
    b /= x ;

    x = -(b/2) + sqrt(((b*b)/4)-c) ;
    E = ((beta*beta*x*x)-1)/(2*beta*x) ;

    double ret = log10(beta) - (asinh(E)*log10(exp(1))) ;

    return ret ;
}

double err(double u){
    M_a = u ;
    V_eq = (M_a*V_a)/M_b ;
    double mm = (3.1+4.4)/2 ;
    double mp = (8.2+10)/2 ;
    double ret = abs(mm - pH(V_m))+abs(mp-pH(V_p)) ;
    return ret ;
}

int main(int argc, char const *argv[])
{
    V_a = 0.05 ;
    M_b = 0.1 ;
    K_a = 1.8 * 1e-5 ;
    K_w = 1e-14 ;
    V_m = 0.02 ;
    V_p = 0.07 ;

    //ternary search
    double l = 0 , r = 1 ;

    while(r-l > 1e-6){
        double l1 = (l*2+r)/3 ;
        double l2 = (l+2*r)/3 ;

       if(err(l1) < err(l2))r = l2 ;
       else l = l1 ;

    }

    printf("%lf\n" , l) ;

    return 0;
}