In this article we report an exact numerical method to calculate the partition function of Benham model which describes the stress-induced denaturation of cloublestrand DNA (dsDNA). Generating function approach is employed in the algorithm. It successes to avoid the catastrophic sign cancellation which makes Benham's algorithm impractical to analyze DNA sequences longer than 10 kb. By our algorithm, the transition profile of base-unpairing probability can be efficiently computed for very long sequences on PC (32-bit, CPU frequency 1.4 GHz). The algorithm requires O(N-3) operations and O(Nlog(2)N) memory, and by contrast Benham's algorithm needs O(N(3)log(2)N) operations and O(N-2) memory. Examples are given on the performance of the implementation on single CPU and multi-CPU cluster. (c) 2005 Elsevier B.V All rights reserved.