Doc/Articles/Play103
At Play With J ... Table of Contents ... Previous Chapter ... Next Chapter
2. Factoring a Number
. By Eugene McDonnell. First published in Vector, 10, 3, (January 1994), 100-105.
You may be among the readers of Vector whose stomachs start churning or whose eyes glaze over when they read the words "tacit definition" in articles about J. This column is meant for you. It is going to take the approach that you know how to write an APL function, and might not be averse to learning how to write a J function in a similar fashion.
[Footnote [1] in the 2nd Edition of APWJ: Even if you don't know APL, Gene's treatment is still readable as an introduction to developing a substantial iterative algorithm in J, notwithstanding the fact that since this article was written computers have become faster and J has acquired a prime factorization primitive (q:). IanClark]
The example I shall use is fairly short -- just 12 lines -- and was shown to me by Joey Tuttle many years ago. It's called factors and it factors a positive integer n. The result of factors n is an ordered list of primes such that n is the product over the list, that is, n <-> */ factors n. It found the factors of the 18-digit integer (2^59)-1 (576460752303423487), which are 179951 and 3203431780337, in 3 min 37.1 s on a Macintosh Quadra 700.
[Footnote [2]: Nowadays factors 576460752303423487x takes 5 s to replicate this result on my Mac Mini, bought in 2005. IanClark]
Various methods of finding the prime factors of a number are given in Knuth, Seminumerical Algorithms, Section 4.5.4 "Factoring into primes", pp. 338-360. The method used here is the simplest one, Algorithm A, which divides the number by increasingly larger primes until all factors have been found, but the method has been vectorized so that more than one factor may be found at once.
The argument must be a positive integer not greater than max, where max is a number derived from the floating-point characteristics of your computer. This discussion assumes that your computer uses IEEE floating point arithmetic, which is the case for PCs, Macintoshes, and most Unix machines. On PCs and Unix machines the maximum is the 16-digit integer (2^53)-1 (9007199254740991). For Macintoshes it is the 19-digit (2^63)-1 (9223372036854775807), or 1024 times as large. These are the largest integers which can be represented to full accuracy on the hardware of those machines.
First I'll show the function definition, then give explanations of its lines, paying special attention to those things that differ from APL. Let's assume for the moment that you created the function below.
factors=: 3 : 0
f=. i. 0 NB. Line 0
t=. 2 3 5 7 11 13 17 19 23 29 31 37 NB. Line 1
o=. +/\ 432 $ 4 2 4 2 4 6 2 6 NB. Line 2
whilst. y >: *: {. t do. NB. Line 3
whilst. #m do. NB. Line 4
m=. (-. * t |!.0 y) # t NB. Line 5
f=. f , m NB. Line 6
y=. y % */ m NB. Line 7
end. NB. Line 8
t=. o + {: t NB. Line 9
end. NB. Line 10
/:~ f , y -. 1 NB. Line 11
)
The line numbers are shown in comments (starting with NB.). These are ignored when the function is executed.
The first thing to notice is that the header line does not indicate whether there is an explicit result or not, nor what its name is if there is one, nor what its valence is, nor what the argument is named, nor which variables are local. Here is how each of these questions is resolved:
Explicit result. Every J program has an explicit result, and it is the value of the last expression executed.
Function name. In J a function is named the same way that a variable is named, by assignment, denoted by (=:). We use the "definition" operator, symbolized by the colon (:), whose result is a function, not a variable. The left argument of the definition operator is a number indicating the result type -- 3 means a function. The rest of the text up to the line with a single ")" is the function definition. This creates the function named factors.
Valence. For an ambivalent function, the monadic case is given by text up to a line containing : and the dyadic case is given by the remaining text. Either block of text may be empty. There is no such thing as a niladic function in J.
Argument name. By convention, the argument to a monadic function is (y). In a dyadic function the left argument is named (x) and the right argument is named (y).
Local variables. J makes a distinction between local assignment (=.) and global assignment (=:). This removes the need for a list of local names. The rule is that on the first assignment of a name using (=.) the name is made local.
Here beginneth the detailed description of the function's lines:
Line 0 sets the initial factor list to empty. J's "iota" is denoted by (i.). By the way, if the argument is 1, this will be the result as well, since 1 is not a prime and has no prime factors. However, it will still be true that n <-> */ factors n, since the product over an empty vector is 1.
Line 1 creates as the initial list of trial divisors the first 12 primes. This enables the function to use just one iteration for all arguments less than 1370 (37^2 is 1369).
Line 2 forms the list of offsets to be used in creating a new list of trial divisors, after there are no more values left in the current list that divide the current value of the number being factored.
In J, the scan operator is defined differently from the APL scan operator. In APL the scan operator implies the reduction operator. In J the function argument to scan is monadic, not dyadic, so that one has to use "sum" (+/) not "add" (+) if we want the sumscan.
There are 432 items in this vector: 54 repeats of the eight items 4 2 4 2 4 6 2 6. The purpose of this vector is to remove multiples of 2 3 5 from consideration as trial divisiors, since these can't be primes. For example,
37 + +/\ 4 2 4 2 4 6 2 6 ... 41 43 47 49 53 59 61 67 ...
Some of the items in these new trial divisors will not be primes. In the list above 49 is composite. Eliminating multiples of 2 3 5 reduces the number of trial divisors by over 73%.
Line 3 starts with the whilst. control word. The control block is the lines up to the matching end. control word. Following the whilst. word is a test that is executed at the end of the control block, which is described with Line 10 below.
Line 4 starts with another whilst. control word and test, which is described with Line 8 below.
Line 5 begins the iteration, and forms the vector m of newly-found factors of the argument. The factors are found by using the vector of trial divisors of the argument as the left argument to the residue function. However, the residue function in J is fuzzed, just as it is in some APLs, in order to permit proper-looking results when used with near-integers, and also to permit the use of decimal rational numbers as arguments to residue. We are dealing with exact integer arguments, so in order to extend the domain of the residue function we require that it not be fuzzed. To accomplish this in APL systems, the comparison tolerance system variable, ⎕ct, is set to 0. In J this is done by explicitly modifying the residue function with the "fit" or "customize" operator (!.) using 0 as right argument. Thus, instead of writing
t | y
We write
t | !. 0 y
and this allows us to work with a residue that has zero fuzz.
Curiously, the "Dictionary of J", which says that the "fit" operator "modifies certain verbs in ways prescribed in their definitions", doesn't describe this use of it with respect to residue. Tsk tsk. The advantage of having the modification of the verb occur in direct connection with its use is obvious: one doesn't have to remember whether or not ⎕ct has been modified, nor run into the hazard of not restoring it when it should be. The use is direct and immediate, and applies only to the function in question. Furthermore, anyone reading the function knows immediately that it is unfuzzed.
After finding the residues, their signum (*) is taken, yielding a vector of 0s and 1s, and these are negated (-.) complementing them to 1s and 0s. This boolean vector is used to select (#) the corresponding items from t, giving in m the new factors to be appended to the result. In APL there has always been a confusion about slash (/): is it an operator or a function? If we write +/1 2 3 it acts like an operator, but when we write 1 0 1/1 2 3 it acts like a function. In J the "#" function is used for the functional case, as in 1 0 1#1 2 3.
Line 6 appends the new factors to the result vector.
Line 7 factors the argument, by dividing (%) it by the product (*/) of the new factors, thus diminishing it.
Line 8 ends the innermost control block. The test #m is false if 0, and true otherwise. Thus the control block is repeated until m is empty.
Line 9 forms the new vector of trial divisors, by adding the offset vector to the last ({:) item of the current vector.
Line 10 ends the outermost control block. Execution is continued if the reduced argument is greater than or equal to (>:) the square (*:) of the first item ({.) of the new vector of trial divisors, since in this case there may be more factors. If it is less, this means that there can't be any new factors in the reduced argument, and thus it is either 1 or a prime.
Line 11 removes 1 from the argument (the result of a -. b is a with items equal to b removed; 17 -. 1 is 17; 1 -. 1 is empty). This leaves it unaltered if it is a prime, or makes it empty otherwise. It then appends the argument to the list of factors (essentially doing nothing if it had been one), sorts (/:~) the list into ascending order, and terminates.
In J the semantics of dyadic upgrade and downgrade have been changed. They no longer have the significance of using the left argument as a collating sequence to produce a grade with reference to it. Instead, they are used to sort the left argument into an order specified by the right argument. The definition of dyadic upgrade (a /: b) is
(/: b) { a
that is, the permutation that puts b in ascending order is used to permute the items of a.
The most frequent use of these sort verbs is with the left and right arguments identical, in which case the result is the sorted argument, either ascending or descending. The reflexive operator (~) applies to a monadic verb to produce a dyadic verb with left argument the same as the right argument. For example, +~1.2 is 2.4, and
/:~ 2 7 1 8 2 8 1 2 2 7 8 8
Below are some examples of the use of the factors function to factor some ten-digit numbers to give you some idea of how it behaves. Notice that the time to factor a number is longest when the argument is a prime, and fairly long also when there are two numbers roughly equal to the square root of the argument.
[Footnote: The reader is invited to compare these results with those obtained on an up-to-date computer. For the final example 6307059911 (a prime number), my Mac Mini (2005) timed the result at 0.011236 s (and 0.055676 s with q: in place of factors). IanClark]
time =: 6!:2 NB. yields seconds to execute its string argument
fmt =: ":!.20 NB. formats (":) numbers to 20 places (!.20)
time 'k=: factors 6307059899'
1.11667
fmt k
7 19 47421503
time 'k=: factors 6307059901'
0.85
fmt k
379 16641319
time 'k=: factors 6307059903'
0.983333
fmt k
3 127 3461 4783
time 'k=: factors 6307059907'
0.65
fmt k
1201 5251507
time 'k=: factors 6307059909'
3.51667
fmt k
3 24749 84947
time 'k=: factors 6307059911'
10.0167
fmt k
6307059911
Vector Vol.10 No. 3
. This article has been updated to reflect the current version of J.