Why not read Frye's article? Here is a digest. You do need modular arithmetic to sieve candidates. Mod 5 works like charm (200 times reduction of the search space). You need to split it in A^4+B^4 and D^4-C^4. The program itself is written in Lisp.
The implementation is reasonably straightforward, while the underlying math is neat. First He sieves out some numbers using modular arithmetics. Then he studies A^4 + B^4 and D^4 - C^4 turning the problem into O(n^2)
First, he notes that mod 5 is very restrictive.
A A^4 A^4 mod 5
0 0 0
1 1 1
2 16 1
3 81 1
4 256 1
Hence, a counterexample must satisfy these two requirements:
A mod 5 = B mod 5 = 0 //A and B are divisible by 5
C^4 mod 5 = D^4 mod 5 = 1 // the reminder of division is 1
Note that D is the largest number and is not divisible by 5, C is not divisible by 5. A < B, A and B are divisible by 5. This reduce the number of possible pairs by a factor of 200. D must be odd using similar considerations for mod 4.
Mod 9 produces only 4 biquadratic yielding 6 congruences of interest. You gain 1.7 improvement from mod 9; 1.4 improvement from mod 13; 1.4 improvement from 29.
If D and C have a common prime factor F then ((AF)^4+(BF)^4 can be eliminated (it was treated earlier). The impact is small, but it helps, nevertheless.
Sieving D^4-C^4 is less efficient, but we still can do certain things. First, we can divide it by 625 before decomposition.
Let N be the number to be decomposed. The second constraint is to limit A, B
//I copy-pasted the rest of the article. Don't read it if you want to reinvent the implementation.
(defun constrained-search (L U)
"Find solution to A' + B' + C ' = D'
with L <= D < U."
(decompose-each (constraint-sieve L U)))
This first pass function constraint-sieve iterates on the D variable over the values between L and
U which are prime with respect to 10. For each D, it iterates on the C variable over the values less than
D which are in proper modulo 625 correspondence with D. When a (D, C) pair passes both the set of
modularity tests and the common kctor and prime-factor tests, then the value N = (D' - C')/625 is
put into a list. The lists for each choice of D are appended together and returned as the vdue of the
function. Here is the code for the first pass function:
(defun constraint-sieve (L U)
"Apply modular and factor constraints to D and C.
Return list of candidates which pass tests."
(loop for D in (prime-10 L U)
append
(loop for C in (good-e-for-d D)
when (and (mod-ok-p D C)
(factor-ok-p D C))
collect
(/ (- (biquad D) (biquad C))
625))))
The second pass function decompose-each takes the list of candidates generated by the first pass
as its argument. It iterates on the variable N over the candidates in the list. For each candidate, it
computes the proper limits on A and iterates between them. When it finds an A such that N - A' is
a biquadrate, it escapes to the outer level and returns some information about the winner. If a winner
is encountered, another function would be needed to fully identify it. Here is the code for the second
pass function:
(defun decompose-each (candidates)
"Attempt to decompose each of the candidates.
Return successful decomposition."
(loop with hd-root = (biquad-rt 0.5)
for N in candidates
for a-limit = (i-biquad-rt N)
for a-min = (ceiling (* halflroot a-lim't))
do
(loop for A from a-min below a-limit
when (biquadrate-p (- N (biquad A)))
do
(return-from decompose-each (list N A))))))
How many candidates can we expect the first pass to produce?
How many A values does the second pass have to scan for the average candidate? Assume that the
candidates are evenly distributed below U', then the average upper limit on A is about 415. (U')'/' =
4U/5 and only 16% of that range is scanned.