Residual Arithmetic

Carrying out integer arithmetic as described will quickly lead to an integer overflow, i.e., the numbers arising in the analysis will soon be larger than the largest numbers that can be expressed in Java. So instead the arithmetic is carried out not on the number themselves but on their residuals with respect to large prime numbers. The largest such prime number is

p₀=2³¹-1 = 2,147,483,647

and the top 32 are listed in the box nearby. (If you really want to know look here for the top 1000.)

For two integers A and N we denote by

A mod N

the remainder of A divided by N. We call the number A mod N the residual of A with respect to (or modulo) N. For example,

38 mod 6 =2 since 38 divided by 6 equals 6 with a remainder 2. Thus the residual of 38 with respect to 6 equals 2. We also say the residual of 38 modulo 6 equals 2.

Taking remainders distributes over addition, subtraction, and multiplication. i.e.,

(x + y) mod N = ((x mod N) + (y mod N)) mod N

(x - y) mod N = ((x mod N) - (y mod N)) mod N

(x * y) mod N = ((x mod N) * (y mod N)) mod N

(It does not make sense for division since the result of dividing two integers may not be an integer.)

When analyzing minimal determining sets and the dependence of some Bézier ordinates on others it only matters whether certain entries in a certain matrix are zero or non-zero. The utility of using remainders instead of the integers themselves rests in the following simple observation:

If A=0 then A mod N = 0 for all integers N.

Thus if we find that the residual of a number Nmodulo a large prime number is non-zero then N cannot be zero. It almost works the other way too: If the residual of N modulo a large prime number is zero then it is very likely (but not certain) that N is zero.

Why should we use prime numbers? Prime numbers become significant when we use several to increase our confidence in treating numbers as zero when all its residuals are zero.

If the numbers with respect to which which we compute residuals have common factors then some combinations of residuals are impossible and the range of numbers that we can distinguish is decreased.

This is best explained with an example.

Suppose we consider residuals with respect to N₁=4 and N₂=6. The following table lists the smallest number that gives rise to a combination of residuals. The columns correspond to the remainder modulo 6, the rows to the remainders modulo 4.

For example,

9 mod 6 = 3 and 9 mod 4 = 1

and there is no number at all that has residual 1 mod 4 and residual 2 mod 6.

Compare this with the choice N₁=3 and N₅=6

The relevant table is now given by

Thus all combinations of residuals are possible and the range of numbers contained in the table is larger, even though the table itself is smaller than that for the previous choice N₁=4 and N₂=6.

These simple observations are made more precise in the celebrated

Chinese Remainder Theorem. Suppose N₁, N₂, ..., N_n are relatively prime. Let

N=N₁ * N₂ * ... * N_n.

Then for any set of residuals n₁, n₂, ..., n_n there exists a unique number

0 <= x < N

such that

x mod N_i = n_i, i = 1,...,n.

Note that a corresponding result does not hold for our choice N₁=4 and N₂=6. Indeed, since 12 has a zero residual both modulo 4 and modulo 6, adding 12 to any number in the first table above gives another number between 0 and 23 that has the same residuals.

A set of numbers is relatively prime if their greatest common factor is 1. Instead of the largest primes available in Java we could thus use the largest available numbers that are relatively prime. However, using prime numbers (which of course are relatively prime) entails only a small loss in efficiency. The reason for using prime numbers, instead of relatively prime numbers, is that prime numbers are widely known and independent of each other. Also note that p₀ = 2³¹ -1 is not only the largest prime number expressible as an integer in Java, it is also the largest integer itself.

Java does allow for a data type long which allocates 64 bits (instead of 32) to every integer. However, that data type is used for computing the residuals of products, by first computing the product of two integers as a long integer. Using long as the basic data type would be possible, but it would require additional programming to compute residuals of products and sums that overflow the long data type and presumably would run more slowly. An effect similar to using the long data type can be obtained by doubling the number of primes used by the applet.

Residuals of Negative Numbers

The residual of A modulo a positive integer N is defined more precisely to be the unique non-negative integer R less than N such that

A = mN + R.

It is thus a number that ranges from 0 to N-1. This is the notion of a residual that is required in applications like ours.

Unfortunately, Java, and other programming languages, botch this notion when computing the residuals of negative numbers.

For example,

-10 = 2 mod 6 since -10 = -2*6 + 2.

On the other hand, Java returns -4 as the residual of -10 mod 6. In general, when computing the residual of a negative number with respect to N one has to add N to the residual to get the correct value. In this example, adding 6 to -4 gives 2, which is the correct residual.

Goliath

Dave Eyre has written a program Goliath that employs residual arithmetic to analyze large sparse rational linear systems. The program has been used extensively in my spline research. For details examine items 33, 35, 38 in my bibliography. The FORTRAN code of Goliath is available (as ACM Collected Algorithm 701) from Netlib, an online software depository. Specifically, it is here.

Dependability of Residual Arithmetic

The following Table shows results computed on the generic double Clough-Tocher split. The columns contain:

• The value of r.
• The value of d.
• The number of equations.
• The number of variables.
• The rank of the matrix.
• The dimension of the spline space.
• The number of residuals required to compute and display the initial minimal determining set.
• The likelihood that no residual of a non-zero matrix entry is zero, based on the assumption that all residuals are equally likely, and using a single prime number 2³¹-1=2147483647

Even for the largest realistic problems that can be handled by the code, the likelihood of an accidentally zero residual is small. For the largest possible problem the likelihood is still much better than even. Even then an accidentally zero residual would be unlikely to cause the computed dimension to be wrong (since it would effect the computation only if it affects a pivot in the Gaussian elimination). It would not follow that the computed dimension is wrong. Moreover, the likelihood of failure can be further reduced, even to the certainty of success, by using several prime numbers.

Residual Arithmetic for Small Prime Numbers

It is instructive to consider the performance of the code when using small prime numbers since then failure may occur more readily.

The following table illustrates the performance of residual arithmetic. The dimension of the spline space on the generic double Clough-Tocher split, for r=1 and d=4 is 39. It was computed using three consecutive prime numbers. The table gives the top one of those prime numbers, and the color indicates the result, as follows:

• Red:

  All entries in the linear system have at least one of their residuals equal to zero. This makes Gaussian Elimination impossible and the matrix is considered to have rank 0.

• Purple:

  The computed dimension is too high. Some non-zero numbers in the linear system are treated as being zero.

• Cyan:

  The dimension is computed correctly. However, some entries in the linear systems have a mixture of zero and non-zero residuals. The program recognizes this as a potential problem. Numbers with mixed residuals are considered non-zero, but they cannot serve as pivots.

• Green:

  The dimension is computed correctly and all non-zero entries in the linear system have three non-zero residuals. This is the expected and desired situation. The smallest triple of primes having these properties are 83, 89, 97.

Notes

• For more complicated linear systems, the primes have to be larger for all the non-zero entries to lead to all non-zero residuals. Here are some examples:
  • r=2, d=7
  • r=3, d=10
• The Table contains all primes through 3449, and selected ranges of larger primes.
• The linear system being analyzed comprises 48 equations in 76 variables. Its rank is 40.
• It's not surprising that for the small prime numbers illustrated in the table the computed dimension is not always correct.
• On the other hand, note that the dimension is often computed correctly even though the residuals suggest that there is a problem.
• For large prime numbers the residuals are all non-zero, as one would expect.
• It is remarkable that the primes for which the dimension is overestimated occur in blocks.