```
Mtest[a_,b_,z_,trials_]:=
Timing[Do[M[a,b,z],{trials-1}];
M[a,b,z]]
```

Thanks again for your experiments. I just tested with the above Mathematica code. The result gives the timing (in secs) to do the computation 'trials' times, and then the M value.

Taking 10,000 trials:

Mtest[1.,4.,50 I,10000]

{1.185608,0.00241259 +0.0599983 I}

Mtest[3.,10.,30+100. I,10000]

{1.747211,-13811.6-953.439 I}

Mtest[15.,20.,200. I,10000]

{1.326008,-2.82866*10^-6-3.27839*10^-6 I}

Mtest[4.,1.,50. I,10000]

{0.046800,-9044.44-18975.1 I}

Mtest[10.,3.,30+100. I,10000]

{0.078001,-1.8468*10^20+9.77876*10^21 I}

Mtest[4000.,4200.,50. I,10000]

{0.249602,-0.868055-0.468886 I}

So, roughly, the results are [0.05 - 1.75] secs to do 10,000 evals. To make a port to C/C++ (or FORTRAN or whatever) worth the effort, I would need a factor of 10x improvement on these times or better, as well as a similar speedup with the Bessel [$]I_{\nu}(z)[$] function. (Also, [$]M(a,b,z)[$] needs to admit complex 'a' and 'z', and [$]I_{\nu}(z)[$] needs complex z). I know you haven't paid attention to speed, but that's where I am on the ballpark requirements.