This paper deals with efficient numerical representation and manipulation of differential and integral operators as symbols in phase-space, i.e., functions of space @math and frequency @math . The symbol smoothness conditions obeyed by many operators in connection to smooth linear partial differential equations allow to write fast-converging, non-asymptotic expansions in adequate systems of rational Chebyshev functions or hierarchical splines. The classical results of closedness of such symbol classes under multiplication, inversion and taking the square root translate into practical iterative algorithms for realizing these operations directly in the proposed expansions. Because symbol-based numerical methods handle operators and not functions, their complexity depends on the desired resolution @math very weakly, typically only through @math factors. We present three applications to computational problems related to wave propagation: 1) preconditioning the Helmholtz equation, 2) decomposing wavefields into one-way components and 3) depth-stepping in reflection seismology.
 A class of matrices ( ( H )-matrices) is introduced which have the following properties. (i) They are sparse in the sense that only few data are needed for their representation. (ii) The matrix-vector multiplication is of almost linear complexity. (iii) In general, sums and products of these matrices are no longer in the same set, but their truncations to the ( H )-matrix format are again of almost linear complexity. (iv) The same statement holds for the inverse of an ( H )-matrix. An algorithm is presented for the rapid direct solution of the Laplace equation on regions with fractal boundaries. In a typical application, the numerical simulation has to be on a very large scale involving at least tens of thousands of equations with as many unknowns, in order to obtain any meaningful results. Attempts to use conventional techniques have encountered insurmountable difficulties, due to excessive CPU time requirements of the computations involved. Indeed, conventional direct algorithms for the solution of linear systems require order O(N3) operations for the solution of an N � N-problem, while classical iterative methods require order O(N2) operations, with the constant strongly dependent on the problem in question. In either case, the computational expense is prohibitive for large-scale problems. The direct algorithm of the present paper requires O(N) operations with a constant dependent only on the geometry of the boundary, making it considerably more practical for large-scale problems encountered in the computation of harmonic measure of fractals, complex iteration theory, potential theory, and growth phenomena such as crystallization, electrodeposition, viscous fingering, and diffusion-limited aggregation.