  This paper describes a practical methodology for computing the Hardy function
Z(t), using just O(((t/epsilon)^(1/3))*(log(t))^(2+o(1)))) standard
computational operations, to a tolerance of epsilon in the relative error. The
methodology is analogous to work published relatively recently by G. H. Hiary,
although the details are very different. Initially the Hardy function is
formulated into sub-sequences of generalised quadratic Gaussian sums of ever
increasing length N. Adapting a theoretical framework formulated by R. B.
Paris, an algorithm for the computation of quadratic Gaussian sums using just
O(log(N)) operations is developed and tested. This algorithmic methodology is
itself incorporated (as a sub-program) into a somewhat larger algorithm
(algorithm ZT13) designed for the computation of Z(t), which is considerably
more efficient than classical O(sqrt(t)) methods, as exemplified by variants of
the Riemann-Siegel formula. Sample computations using algorithm ZT13 are
presented, the results of which conform to technical theoretical predictions
made regarding the limits of its performance (both in terms of its operational
count and degree of precision). Speculative ideas briefly floated in the
conclusions suggest that adaptations to this methodology have the potential to
deliver significant further computational savings in the future.
