#includes <vector>
#include "gammaCDFderiv.hpp"
real lgammaJacLinf(double x[N], double alpha, double lambda){
		vector[num_elements(x)] prob;
		real my_target = 0;
		real num_elem = num_elements(x);
		for (i in 2:num_elements(x)){
			for(j in 1:i){	
				real gam_1 = gammaCDFderiv(alpha, lambda*x[j]);
				real gam_2 = gammaCDFderiv(alpha, lambda*x[i]);
				real gamma_term = tgamma(alpha + 1);
				real denom_1 = pow(lambda*x[j], alpha) * exp(-lambda*x[j]);
				real denom_2 = pow(lambda*x[i], alpha) * exp(-lambda*x[i]);
				real term_1 = gamma_term * gam_1 / denom_1;
				real term_2 = gamma_term * gam_2 / denom_2;
				
				my_target += x[i] * x[j] * fabs(term_1 -term_2);
			}
		}
		return log(pow(lambda*alpha, -1) * my_target);
}