function CDF = digital_option_price_gpu(S0, K, T, CP, sv1, sv2, sv3, rf, zeta, N, L, rho1, v1, k1, s1, rho2, v2, k2, s2, c0, c1, c2, muy, sy, mv, rhoj)
% This function allows one to price digital options and get resulting risk-neutral CDFs
% for a given known underlying distribution (see the function: moment_gen_function.m).
% Updated to run on a GPU to speed up computations.
%
% The below is lightly adapted code from Nicola Fusari's Matlab Toolbox for
% Option Pricing (including comments). This code edits (with some
% notational hacks) his code to obtain the price of a digital option
% (should work for digital calls and puts, but I've tested it only for
% calls, so set CP=1) following Fang and Oosterlee (2008). The option pays
% K if S_T>K and 0 otherwise. It then transforms the resulting price to a
% risk-neutral CDF evaluated at the same strikes.
%
% Nicola's code:
%
%
% The option price approximation is based on the paper by Fang & Oosterlee 
% (2008), SIAM J. Sci. Comput. Vol 31, No. 2, pp. 826-848.
%
% I've tryied to keep the notation as close as possible to the one in the
% paper.
%
% OUTPUT:   RN CDF
% INPUT:    Distr ---> select the distribution of the underlying (see
%                      function: moment_gen_function.m). For GPU,
%                      hard-coded to run 2F_Gaussian to avoid character
%                      function input.
%           S0    ---> underlying initial value
%           K     ---> column vector of strikes
%           x     ---> vector that contains the parameters for the selected
%                      distribution (see function: moment_gen_function.m).
%                      For GPU version of function, I've extracted the
%                      individual elements.
%           T     ---> maturity (in YEARS)
%           CP    ---> Call/Put Flag  (call=1, put=-1)
%           sv1   ---> value of the first variance component (if needed)
%           sv2   ---> value of the first variance component (if needed)
%           rf    ---> risk-free rate used to discound the option price
%           d     ---> dividend yield
% REMARK: S0 and K must have the same length ---> COLUMN VECTORS
% Nicola Fusari (2011):
% You can distribute and use the code. Please cite the author if you use
% the code in in your research.

% Compute numerical derivatives of the MGF to computed MEAN and VARIANCE of
% the underlying distribution to decide WHERE to TRUNCATE the DENSITY
% approximation:
delta = 0.01; % 0.01

cf1 = log(moment_gen_function_gpu(rho1, v1, k1, s1, rho2, v2, k2, s2, c0, c1, c2, muy, sy, mv, rhoj, -1.*delta./2  , 0, 0, 0, sv1, sv2, sv3, T, rf, zeta));
cf0 = log(moment_gen_function_gpu(rho1, v1, k1, s1, rho2, v2, k2, s2, c0, c1, c2, muy, sy, mv, rhoj,            0  , 0, 0, 0, sv1, sv2, sv3, T, rf, zeta));
CF1 = log(moment_gen_function_gpu(rho1, v1, k1, s1, rho2, v2, k2, s2, c0, c1, c2, muy, sy, mv, rhoj,  1.*delta./2  , 0, 0, 0, sv1, sv2, sv3, T, rf, zeta));

%   http://mathworld.wolfram.com/FiniteDifference.html

d1 = (CF1 - cf1)./delta;
d2 = (CF1 + cf1 - 2.*cf0)./((delta/2)^2);  

a = real(d1 - L*sqrt(abs(d2)));
b = real(d1 + L*sqrt(abs(d2)));

if CP==1
    c = 0; d = b;
else
    c = a; d = 0;
end

dt = T;

X = log(S0./K);

VC_sum=0;
for u_j=1:N
    uu_j = (u_j-1)*pi./(b-a);
    
    % Psi_k
    if u_j==1
        PSI = d-c;
    else
        PSI = (sin(uu_j.*(d-a)) - sin(uu_j.*(c-a))).*(1./uu_j);
    end

    % Characteristic function of log(S_T/S_t)
    cf = moment_gen_function_gpu(rho1, v1, k1, s1, rho2, v2, k2, s2, c0, c1, c2, muy, sy, mv, rhoj,1i.*uu_j,0,0,0,sv1,sv2,sv3,dt,rf,zeta);

    % Coefficients
    vk = CP.*K.*PSI.*2./(b-a);
    fk = real(cf.*exp(1i.*uu_j.*(X-a)));
    if u_j==1
        fk = 0.5.*fk;
    end
    VC_sum = VC_sum + fk.*vk;
end
VC = max(exp(-rf*dt).*VC_sum,0.0);

CDF = 1-VC.*exp(rf*T).*(1./K);
