Initial variables:

A : Angle of elevation, in radians
V : Velocity at release, in m/s
M : Mass of projectile, in kg
C : Coefficient of drag, pure number
D : Diameter of head of projectile, in m
F : Drag factor, in m^(-1)
T : Time increment, in s
H : Height at release, in m
G : Acceleration due to gravity, in m/s^2

F = -0.473595 * C * D^2 / M
T = 0.004 : for most computations
T = 0.001 : to be painfully accurate
G = -9.80665

Loop variables:

v : Velocity, in m/s
x : Horizontal velocity, in m/s
y : Vertical velocity, in m/s
s : Range, in m
h : Height, in m
t : Flight time, in s

Prior to the first iteration:

v = V
x = V * cos ( A )
y = V * sin ( A )
s = 0
h = H
t = 0

Iterate generating new values for the loop variables:

sA : Sine of angle of travel
cA : Cosine of angle of travel
d : Acceleration due to drag
dx : Horizontal acceleration
dy : Vertical acceleration

sA = y / v : avoid doing trig functions
cA = x / v : avoid doing trig functions
d = F * v^2 : drag depends on square of velocity at these speeds

dx = d * cA : drag slows horizontally
dy = G + d * sA : drag slows vertically, plus gravity slows on ascent and accelerates on descent

s = s + ( x * T ) + ( 0.5 * dx * T^2 ) : s = ut + 1/2 at^2
h = h + ( y * T ) + ( 0.5 * dy * T^2 ) : s = ut + 1/2 at^2

x = x + ( dx * T ) : v = u + at
y = y + ( dy * T ) : v = u + at
v = sqrt ( x^2 + y^2 ) : avoid doing trig functions

t = t + T : tick

After the final iteration, correct s, h, and t so that they match the exact moment of impact with ground:

dh : Correction to h
ds : Correction to s
dt : Correction to t

dh = h
h = h - dh
dt = dh / y
t = t - dt
ds = x * dt
s = s - ds

And that's it.

(From: James Prescot)