b + i j + i (t 1)

f ( ) (2eb / )g The integral in g is calculated by using the trapezoidal rule approximation with a division of the integration interval (0, 2) into 2m subintervals of length 1/m for an appropriately chosen value of m It is recommended to take m = 4M This gives 2m 1 f + f2m 1 p + 0 , (F3) fp cos g 2m m 2

b+i j +i (t 1)

j =1

b + i j + i (p/m 1)

The approximation of (2eb / )g to f ( ) is extraordinarily accurate Rather than calculating from (F3) the constants g for = 0, 1, , M 1 by direct summation, it is much better to use the discrete Fast Fourier Transform method to calculate the constants g for = 0, 1, , 2m 1 More important than speeding up the calculations, the discrete FFT method has the advantage of its numerical stability To see how to apply the discrete FFT method to (F3), de ne gk by gk =

1 2 (f0 fk , + f2m ),

(F4)

for = 0, 1, , 2m 1 The discrete FFT method can be applied to this representation Applying the inverse discrete FFT method to the vector (g0 , , g2m 1 ) yields the sought vector (g0 , , g2m 1 ) Here is a summary of the algorithm: Input: M, , b, n and m Output: f (k ) for k = 0, 1, , M 1 Step 1: Calculate for p = 0, 1, , 2m and 1 j n,

b + i j + i (p/m 1)

1 2 (f0 + f2m )

n Next calculate fp = j =1 j fjp for p = 0, 1, , 2m Let g0 = for k = 1, , 2m 1 and gk = fk

Step 2: Apply the inverse discrete FFT method to the vector (g0 , , g2m 1 ) in order to obtain the desired vector (g0 , , g2m 1 ) Step 3: Let f ( ) = (2e b / )g for 0 M 1 In step 3 of the algorithm g is multiplied by eb In order to avoid numerical instability, it is important to choose b not too large Assuming that the ratio m/M is large enough, say 4, numerical experiments indicate that b = 22/m gives results that are almost of machine accuracy 2E 16 (in general, it is best to choose b somewhat larger than ln( )/(2m) where is the machine precision) If f is suf ciently smooth, it usually suf ces to take n = 8, otherwise n = 16 is

recommended The parameter M is taken as a power of 2 (say, M = 32 or M = 64) while the parameter m is chosen equal to 4M The choices of M and are not particularly relevant when f is smooth enough (theoretically, the accuracy increases when gets smaller) In practice it is advisable to apply the algorithm for and 1 to see whether or not the results are affected by the choice of 2 Non-smooth functions The Den Iseger algorithm may also perform unsatisfactorily when f or its derivative has discontinuities In such cases the numerical dif culties may be circumvented by using a simple modi cation of the algorithm To do this, assume that f (s) can be represented as (F5) f (s) = v(s, e x0 s ) for some real scalar x0 and some function v(s, u) with the property that for any xed u the function v(s, u) is the Laplace transform of a smooth function As an example, consider the complementary waiting-time distribution f (t) = P {Wq > t} in the M/D/1 queue with deterministic service times D and service in order of arrival; see 9 This function f (t) is continuous but is not differentiable at the points t = D, 2D, The Laplace transform f (s) of f (t) is given by f (s) = s + e sD , s[s + e sD ] (F6)

where is the average arrival rate and = D < 1 Then (F5) applies with x0 = D and v(s, u) = s + u s(s + u)

In this example we have indeed that for any xed u the function v(s, u) is the Laplace transform of an analytic (and hence smooth) function In the modi ed Den Iseger algorithm the basic relation (F2) should be modi ed as n 1 eb j v j (t) cos( (t + 1)) dt (F7) f( )

j =1 1