I have a paper explaining the algorithms used here.
Given integers p[0], p[1], ..., p[n-1] and q[0], q[1], ..., q[n-1], sortedpq enumerates pairs (i,j) in increasing order of p[i]+q[j]. It uses extra space for just n integers and 2n indices. sortedpp is a faster version of sortedpq when p = q; it enumerates pairs (i,j) with i >= j, in increasing order of p[i]+p[j].
solvepqrs enumerates (i,j,k,l) such that p[i]+q[j] = r[k]+s[l], in increasing order of p[i]+q[j]. solvepppp, solvepprs, solvepqpq, and solvepprr are faster versions of solvepqrs.
Several sample applications are included in the sortedsums package:
This work was supported by the National Science Foundation under grant DMS-9600083.
The fourth power of 8707481 is a sum of three positive fourth powers. Same for 12197457, 16003017, and 16430513. Other small known solutions are 422481 (Frye), 2813001 (MacLeod), 20615673 (Elkies), and 638523249 (MacLeod). There are no other primitive solutions below 21000000. This computation used some observations by Morgan Ward.
48988659276962496 (discovered independently by Wilson), 391909274215699968, and 490593422681271000 can be written in 5 ways as the sum of two positive cubes.