You may use Kahan summation algorithm [1]
The idea is to reschedule the sum operations in such a way precision loss is limited. The code is very simple (reproduced from [1] below).
If this does not suffice, you may use multiprecision representations, such as quad doubles [2]. They are supported by several languages / compilers (including GNU c). Finally, if this does not suffice, you may use multi-precision arithmetics based on expansions [3]. I have written an implementation with a C++ class that can be used like a number [4]. Now if the underflow occurs within a single term (i.e., if the exp() underflows), then you may need to use floating point representations with larger number of bits. For instance, you can use MPFR [5]. It has implementations of all classical operators and functions with arbitrary precision (very good library, implemented by friends of mine :-).
function KahanSum(input)
var sum = 0.0
var y,t // Temporary values.
var c = 0.0 // A running compensation for lost low-order bits.
for i = 1 to input.length do
y = input[i] - c // So far, so good: c is zero.
t = sum + y // Alas, sum is big, y small, so low-order digits of y are lost.
c = (t - sum) - y // (t - sum) recovers the high-order part of y; subtracting y recovers -(low part of y)
sum = t // Algebraically, c should always be zero. Beware overly-aggressive optimizing compilers!
next i // Next time around, the lost low part will be added to y in a fresh attempt.
return sum
[1] https://en.wikipedia.org/wiki/Kahan_summation_algorithm
[2] https://en.wikipedia.org/wiki/Quadruple-precision_floating-point_format
[3] J. Shewchuk, Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates, Discrete and Computational Geometry (Impact Factor: 0.69). 07/1996; 18(3). DOI: 10.1007/PL00009321
[4] https://gforge.inria.fr/frs/?group_id=5833, MultiPrecision_psm
[5] MPFR: http://www.mpfr.org/