Two modes (i.e. ground rings) are possible in QFermat, rational arithmetic and modular.

A session with QFermat starts in rational mode. All arithmetic is that of rational numbers: to add you get a lowest common denominator, etc.

Changing Modes

To leave rational mode and enter modular mode, the command is   &p.

The interpreter will display the message "Changing arithmetic mode" and will wait for the user to enter an integer n which will become the modulus. n must be at least 2 and no more than 2^24 - 1 = 16,777,215. Arithmetic will be slightly faster if n is no more than 2^16 - 1 = 65535.

The current values on the environment stack will be converted as follows: for the rational number x = r/s, compute r mod n, then s mod n, then divide the two by inverting s. If s is not invertible mod n, x will be set to 0 and a warning printed.

To change back to rational mode, just enter   &p  again. You can change back and forth as often as you wish.

If you enter a constant from the keyboard in modular mode, it will be reduced modulo the modulus.

It is possible to temporarily turn off modular arithmetic in modular mode. This is useful for the following reason. Suppose you want a loop in a function to execute m times and you write

The variant form  &m( < expression > ) turns modular mode off, computes the expression, then restores modular mode.

Never do modular arithmetic, especially multiplication, with outsized values that have been created with &m. (They may be used in arithmetic that is enclosed by &m.) At the very least garbage will be introduced. The machine may crash. When modular mode is in effect, use these values only for special purposes, such as to control loops.

However, it is permissible and sometimes convenient to set a variable equal to the modulus, such as  mod : = &m(...). Constructions of the form x | mod are legal in modular mode.

Don't make the mistake of inadvertently turning modular mode on when you thought you were turning it off. For example, if you are in modular mode with modulus 7 and wish to access entry 7 in array x, the well intentioned command x[&m(7)] will fail, because Fermat automatically turns modular mode off when it computes array coordinates. You inadvertently turned it back on.

A message of acknowledgement will be printed if you execute  &m  from the keyboard, but not if it is executed inside a function.

Disabling modular mode in this way creates a small problem for saving and loading. If a variable is saved with a value of, say, -5, and the file later read during modular mode, the -5 will not be loaded, rather it will be reduced modulo the modulus. However, modular arithmetic may be disabled during polynomial read-in, so a polynomial with negative or outsized coefficients will be read as is. See the later discussion of "polynomial read-in" and &n.
 

Mode Conversion Automatically Stored

If you are in modular mode when you save to output, the command

&(p = < modulus > )

with the correct number will be inserted in the file. Thus, in a later session, when you read that file, conversion to modular mode will be done automatically. Of course, if you wish, you may insert this command yourself in an input file. The syntax is   &(p = < modulus > )  . When the file is read, the expression   < modulus >  will be evaluated.

Fast, Selective Mode Conversion

This is one of the "Dangerous Commands." See that later chapter.
 
 

Polynomials
 
 

Having chosen an arithmetic mode, the user may change the ground field (or ring) by converting it into a polynomial ring. This is done by adjoining variables - as many as desired - that remain unevaluated. The names of these variables follow the usual rules of variable names in Fermat. To adjoin the variable "t", enter the adjoin polynomial command, & ¶. You will be prompted for the name of the variable. ¶  is option-7 on the keyboard. Variables added later have higher precedence than those earlier. Among other things, this means that if you adjoin t and then u, the polynomial u*t is displayed as (t)u, not (u)t.

To place this command in an input file, use the imperative form &(¶ = t). If you have attached this polynomial variable and later enter the save command, the command &(¶ = t) will be automatically included in the saved file. Repeat for other variables.

A variable name can be read from an array, just as the name of a file can be read from an array during read and save commands. This allows the easy creation of many names. For example, to create variables x1, x2, x3, ..., x9, set an array, say [x] equal to 'x ' (note the blank). Then write a loop for i=1, 9 do x[2] : = i+48; &(¶ = [x])  od. (The rather puzzling "48" converts integers to ASCII.)

The polynomial variable t can be dropped or cancelled by entering & ¶  , followed by -t.

A name once chosen for a polynomial variable cannot be used for any other purpose. If a variable previously existed called "t", it will be inaccessible.

Exponents must be at least 0, unless the Laurent option has been chosen - see that chapter for more information. In any case, they must have absolute value less than 2²⁸. (Technically, an exponent can be up to 2³¹ - 1, but Fermat will not let you directly assign an exponent larger than 2²⁸ - 1.)

Producing an exponent larger than 2³¹ by multiplication or exponentiation will probably not result in an error message. But later results will be garbage!

Letting F denote the ground ring, suppose variables have been adjoined to create the polynomial ring F[t, u, v, ...]. Arithmetic is done in the obvious way, with the ordinary signs +, -, *, ^  . The div operator \ and the mod operator | work about as one would expect, at least if there is only one variable. Division, using / as in x/y, reduces the quotient by dividing out the greatest common divisor of the two polynomials. The result may be a polynomial, if y divides evenly into x, or, more likely, will be an element of the quotient field F(t, u, v, ...). Such elements are called quolynomials, and are written as fractions or quotients of polynomials. This topic is discussed in greater detail in the following chapter. There is one mathematical complication. In modular mode, F might not be an integral domain. Then neither is F[t, u, v,...]. Therefore, there is no quotient field, and writing quotients of polynomials is a bit of a sham. Fermat basically ignores this. When division is indicated, the same g.c.d. algorithms are invoked as in the field F case. Fermat leaves it to the user to make sure that the results computed have any meaning. Even more trouble occurs when polymods are used with division. See the later chapter on "Polymods."

The operator "div", \ is essentially obvious when applied to single variable polynomials, but not to multivariable ones. For in that case, there is no division, there is only pseudo-division. (See any text on abstract algebra, or Knuth volume 2.) Fermat will yield the result of pseudo-division in this case. Similarly, "mod", |, returns the pseudo-remainder. Ideally, x \y = q and x |y = r mean c^k x = q y + r, where c = the leading coefficient of y and k >= 0 is an integer. (QFermat adopts the convention that k = deg(x) - deg(y) + 1.) However, if y is at a lower level than x, the equation c^k x = q y + r will probably fail because y is divided into each coefficient of x and a different k will probably arise in each case.

Polynomial Built-in Functions

Many of the built-in functions cannot be given polynomials of positive degree. Those that can, like greatest integer, work upon all the coefficients.

There is a special operator for polynomial evaluation, denoted á, which is option-8 on the keyboard. Suppose first that only one polynomial variable, say t, exists. Then xáy replaces every occurrence of t in x with y. For example, if x = t³ + 3t² + 3t + 1, then xá2 = 27. y can be any expression evaluating to a polynomial.

If several polynomial variables exist, then xáy replaces the top-most polynomial variable in x with y, i.e., the latest created. To replace one of the other variables, say u, use the alternate form xá(u = y). The name of any polynomial variable may be used instead of u, and any expression may be used for y.

A fast shortcut form of evaluation called total evaluation exists. To evaluate x at every variable, use the syntax xá(v₁, v₂, ...), where all the v_i are numbers. There must be a number corresponding to each polynomial variable, in the precedence order - highest precedence (last attached) listed first. In QFermat the following more general syntax is allowed. Suppose there are five poly vars, e, d, c, b, a in that order (e last and highest). Then qá(d = w,x,y) will replace each d in q with w, each c with x, each b with y. e and a are untouched. Further, w, x, and y can be arbitrary quolynomials.

Fermat has a built-in function for polynomial coefficients. Syntax of use is ç(x, n), or ç(x, n₁, n₂, ...) if several polynomial variables have been adjoined. ç is option-c on the keyboard. x can be any expression. The n_i must be numbers. In modular mode, modular arithmetic is ignored while the n_i are evaluated.

There is a built-in function to compute the degree of a polynomial, using the symbol ¶. Syntax: ¶x or ¶(x, i). In the second form, the i names the i^th polynomial variable. ¶x returns the largest exponent within x of the specified polynomial variable. If no specification is given, the highest precedence variable is assumed. In modular mode, it returns an actual integer, not reduced modulo the modulus. ¶ is option-d on the keyboard.

Similarly, there is a built-in function to compute the "codegree" of a polynomial, using the symbol ^o. It returns the smallest exponent of the specified polynomial variable. Syntax as with ¶. ^o is option-k on the keyboard.

There is a built-in function to compute the leading term of a polynomial, using the symbol £  (equivalent, Lterm).  £ is option-3 on the keyboard. Follow it with the expression whose leading term is to be computed. By definition, the leading term of a polynomial is coef * uⁿ where u is the variable of highest precedence, n is the highest degree, and coef is the coefficient of uⁿ.

Fermat has a built-in function to compute the derivative of a polynomial, using the symbol ^¢. This is option-shift-] on the keyboard. Don't confuse it with apostrophe. Precede it with the expression whose derivative is to be computed, such as x^¢. Differentiation is with respect to the variable of highest precedence.

There is a function Remquot to return the (pseudo)remainder (r, say) of dividing d into x and the (pseudo)quotient, q. c^k x = q d + r, where c is the leading coefficient of d and k = deg(x) - deg(d)+ 1 (unless d is a number or c is invertible; then k = 0).

There are functions Factor, Sqfree, and Irred to factor one variable polynomials and test for irreducibility, both over finite fields. These are described in detail below.

See the earlier chapter on built-in functions for descriptions of Level, Height, Raise, Lower, Lcoef, Content, and Numcon.

Polynomial Read-in

Fermat has a feature that facilitates the reading of large polynomial constants, the polynomial read-in, ¶. It is used by simply putting ¶  in front of a polynomial constant that is written out fully following Fermat's conventions.

Example:

where the two polynomial variables t and u have been attached in that order (u has higher precedence, but there could be polynomial variables between them).

The point of this feature is speed. If the ¶  is omitted in the above example, Fermat will square t, multiply by 5, multiply 2 times t, add to the previous result 5t², then add 1, then square u, then multiply by the previous result, etc. - in other words, it will evaluate the expression normally. But since the expression is written in the canonical order, y can be pieced together without any calls to addition or multiplication at all. The ¶  signals Fermat to expect the canonical order and simply " attach" the terms as they come to y. This provides a very large saving of time, which is very noticeable even for polynomials of moderate size.

The canonical order is that of decreasing exponents, no zero exponents, coefficients written first in each term, and nesting according to the precedence of the polynomial variables. There must be no spurious parentheses, no multiplication signs *, no division /, no extra plus signs, like +7t-3, and no 0's (as stand-alone constants). There must be no variables or any other built-in functions. Here is another example:

s : = ¶  v³ + (3u - 2t + 3)v² + (3u² + (-4t + 6)u +4t² - 4t + 6)v + u³ + (-2t + 3)u² + (4t²+ 4t + 6)u - 7t³ + 4t² - 4t - 2

If you are unsure of some of the fine points of canonical form, create some polynomials and watch how Fermat displays them.

During polynomial read-in in modular mode, one sometimes wants modular arithmetic disabled. One may want the read command to load correctly (i.e., as is) data that was created when modular mode was off. Note, however, that non-polynomial data of that type will not be read "correctly." On the other hand one may indeed want the coefficients to be modded out. The toggle switch &n allows the user to choose.

You can pretty much ignore polynomial read-in, especially if you use only one polynomial variable and have no polynomials of more than 15 or 20 terms. However, Fermat will often insert it in saved files and on the screen when it dumps a list of variables, so you must at least be aware of its existence.

Factoring Polynomials

Fermat allows the factoring of monic one-variable polynomials over any finite field. The finite field is created by simply being in modular mode over a prime modulus, or by additionally modding out by irreducible polynomials to form a more complex finite field, as described in the later section "Polymods." Factoring into irreducibles or square-free polynomials is possible, or polynomials can just be checked for irreducibility. Square-free factoring applies to more general ground rings, not just finite fields.

For factoring, the syntax is Factor( <poly> , [x]) or Factor( <poly> , [x], < level > ). The factors of <poly> will be deposited into an array [x] having two columns and as many rows as necessary. (The number of factors (rows) is returned as the value of Factor.) In each row, the first entry is an irreducible polynomial p(t) and the second is the largest exponent e such that p(t)^e divides <poly> .

In the second form, < level > specifies the subfield to factor over. For example, suppose over the prime p = 30203 two polynomial variables t and u have been attached, and the irreducible polynomial t⁴ + t³ + t² + t + 1 has been modded out to form the field of order p⁴. Then Factor(u¹²+1, [x], 0) returns six quadratic irreducible factors - irreducible over F_p. So setting level = 0 suppresses the existence of the field of order p⁴. Factor(u¹²+1, [x], 1) produces a factorization over the ground ring F_p plus the first modded-out polynomial, i.e. the field of order p⁴. The result is twelve irreducible factors involving t. Note that the same twelve factors could be attained much more rapidly by factoring each of the six quadratics obtained at level 0. Similarly, if now a further field extension is created by modding out on u, the most efficient way to factor is level-by-level, to the extent possible. It is therefore best for factoring to use as many variables t, u, ... as possible in creating the field.

Sqfree is exactly like Factor except it produces factors that are square-free only, and can be done over more general ground rings.

Irred tells if its argument is irreducible, and, if not, describes the factorization. Syntax is Irred( <poly> ) or Irred( <poly> , < level > ) ("level" is explained above). The value returned is as follows:

    -1 means can't decide (too many variables, for instance).

     0 means the argument is a number or is a field element.

     1 means irreducible.

     n > 1 means the argument is the product of n distinct irreducibles of the same degree (and is therefore not irreducible).

     x, a poly, means x is a factor of the argument (which is therefore not irreducible).

Fermat uses the algorithms of Cohen, " A Course in Computational Number Theory," Springer Verlag, 1993, pages 123-130. There is some randomness in these algorithms, so the time it takes to factor a given example can vary.
 

Quolynomials
 
 
 

Letting F denote the ground ring, suppose variables have been adjoined to create the polynomial ring F[t, u, v, ...]. Division, using / as in x/y, reduces the quotient by dividing out the greatest common divisor of the two polynomials. The result may be a polynomial, if y divides evenly into x, or, more likely, will be an element of the quotient field F(t,u,v, ...). Such elements are called quolynomials, and are written as fractions or quotients of polynomials. Fermat automatically creates quolynomials whenever it is necessary.

Basic arithmetic of quolynomials is essentially obvious to anyone who remembers high school algebra.

There are three mathematical complications. In modular mode, F might not be an integral domain. Then neither is F[t, u, v, ...]. Therefore, there is no quotient field, and writing quotients of polynomials is a bit of a sham. Fermat basically ignores this. Certainly these polynomials can be added, subtracted, and multiplied with no problem. When division is indicated, there may or may not be an immediate problem. For example, 1/(2t+2) will cause an error message if F = Z₄, because 2 is not invertible in F. If some computation f/g does not cause this problem, the same g.c.d. algorithms are invoked as in the field F case. Fermat leaves it to the user to make sure that the results computed have any meaning.

Secondly, a similar problem arises when polynomials p, q, r, ... are used to mod out and form the quotient ring F[t, u, v, ...]  / < p, q, r, ... > of polymods. If the polynomials p, q, r, ... are not all irreducible, then this ring is not a field, and dividing and forming quolynomials is a sham. Fermat will not stop you from trying to do this, and there is no guarantee of the results. In particular, there will no longer be canonical forms, i.e., different looking expressions may in fact be equal. If they are subtracted, the result may or may not be an honest zero. The following chapter describes polymods.

Div and mod applied to quolynomials act only on the numerators.

There are built-in functions Numer and Denom that return the numerator and denominator of a quoly. The built-in function degree ¶ returns the difference between the degrees of the numerator and denominator. Leading coefficient ç ignores any denominator. Leading term £ produces an error if called on a non-polynomial.
 
 

Polymods
 
 
 

Choosing an arithmetic mode establishes the ground ring F. On top of this may be attached any number of unevaluated variables t,u, ..., thereby creating the polynomial ring F[t, u, ...] and its quotient field, the field of rational functions or quolynomials. Further, certain polynomials p, q, ... can be chosen to mod out with, creating the quotient ring F[t, u, ...]  / < p, q, ... > , whose elements are called polymods, and, implicitly, quolymods, which are expressed as quotients of polynomials, superficially the same as ordinary quolynomials. If the polynomial p involves the variable t, we say that t has a relation imposed on it.

For example, in rational mode, F[t]  / < t²+1 > is the mathematical ring Z(i), with t playing the role of ÷(-1).

The basic command to mod out by a polynomial is &P. You will be prompted for the quotient polynomial, which must be monic, say tⁿ + c₁t^n-1 + .... The imperative form of this command is &(P = tⁿ + c₁t^n-1 + ...).

In principle, any monic polynomial may be modded out, say tⁿ + c₁t^n-1 + .... However, QFermat is best at the case where the quotient ring becomes a field (note well, Q[t, u, ...]  / < p, q,... > is a field, Z[t, u, ...]  / < p, q, ... > is not). Specifically, suppose the polynomial variables have been attached in the temporal sequence t, u, v, .... Begin by modding out a monic irreducible polynomial p(t) such that F₁ = Z[t]/ < p > is an integral domain, and its field of fractions is Q[t]/ < p > . Then, if desired, mod out by a monic polynomial q(u,t) such that F₂ = F₁[u]/ < q > is an integral domain, and continue in this manner always creating an integral domain, and, by the same stroke, its field of fractions.

You must tell QFermat that a field will result, and it is your responsibility to check this. Do this at each step by adding a comma and 1, as &(P = tⁿ + c₁t^n-1 + ..., 1). In QFermat, you may append a list of primes q such that the modder p(t) remains irreducible mod q. For example, &(P = t³ + t² + t + 2, 1: 151, 167, 191,839, 859, 863, 907, 911, 991). Add as many as you like, at least 30 is desirable. Or, probably better, omit the list and QFermat will compute it for you. (This takes a while, so if you save the session to a file, QFermat will include the list in the saved file and just reload it next time.) QFermat then computes a second list of auxilliary primes of a different sort: modulo these primes the modding polynomial has a root. Both types of primes are used to speed up g.c.d. computations. If Fermat cannot find enough of either type, it will tell you and instruct you how to get more, using the commands &A and &B.

When you have implemented a field in the above manner, division is done correctly. For example, with the above p(t), 1/t² is computed as ( t² - t - 1 ) / 4 , which is correct. If you have not told QFermat that a field results, 1/t² is left as 1/t². Whether this has any meaning is, of course, your problem.

To later rescind modding out on t, enter &P again and then -t. To change the modding polynomial to another, you must first rescind the current one.

Note that polynomial evaluation (substitution) of a variable involved in a modding relation is not a homomorphism, and may not make sense. This is the user's responsibility.

If a field has been formed in this way, one can drop into modularmode mod any prime such that the modding polynomials remain irreducible. It is the user's responsibility to check this.
 
 

Laurent Polynomials
 
 
 

In Fermat a Laurent polynomial is one with negative exponents. If you wish to allow this, activate the toggle switch &±. (± is option-shift-= on the keyboard.) All of the variables you have created up to that point will be converted to the new format. For example, 1/(t²+2t) will become t^-1/(t+2). No negative exponents are left in the denominator of a quolynomial and all positives are factored out and moved to the numerator.

Unlike polymods, Laurent polynomials present no inherent mathematical difficulty. Be aware that in Laurent mode, the function GCD treats variables t as invertibles, and always presents answers with no negative exponents.

Laurent mode should be used if you expect a fair amount of results to have denominators that a variable could be factored out of. If a matrix has several entries such as 1/t, use this mode to compute its determinant or inverse.

Laurent polynomials implement computations in the group ring Z(Z).
 
 

Character Strings
 
 
 

Fermat allows the assignment of characters to scalar variables and character strings or literals to arrays. In the first case the syntax is  y : =  ¢a¢ to assign the character  ¢a¢ to the variable y. In the second case the syntax is  [x] : =  ¢sample¢. (Use the single quote, or apostrophe.)

What actually happens is that y gets the ASCII code of ¢a¢. It is still therefore "really" a number. If one simply types y or ¢a¢, the ASCII code is displayed. (This works in all arithmetic modes. In modular mode the result will not be reduced modulo the modulus.) To see the character 'a' instead, use the display command  !(y: char). Note the reserved word char.

Strings are more useful than characters, for which one needs arrays. If the user first creates an array [x] of size, say, six, and then enters the command [x] : =  ¢sample¢, each character is stored in the six slots of [x]. Entering ![x will display the ASCII codes of the six letters. To see the letters, one needs the char attribute, the syntax of which is  ![x: char. (Note that this is array short form.) In the above example, this will produce the original string sample. A spacing attribute can be added, as in  ![x: char : n  , which places each character in the left of a field n units wide. Compare the display command form,  !([x: char). One can also use array long form  ![x]: char, which produces one character per line, or  !([x]: char). Try them both.

Array features discussed in earlier chapters enable one to work with substrings, concatenate, find the length of a string, etc. For example, to concatenate a string stored in [x] with a string stored in [y], enter   [z] : = [x]f [y]. (This could also be executed more cumbersomely with subarrays.)

Character strings can easily be read into an array. The command is ?[a]:char , which will display on the next line > [a]  : =  ¢ . Note the prompt > and single quote. The user then enters any character string followed by the single quote and hits return. Then array [a] gets the ASCII codes for the characters.

There are various fine points, which are here explained:

(1) If the array [x] has some undefined elements, they will appear as ¢*¢ upon execution of ![x: char. Note the space after the *.

(2) If [x] has not yet been created, [x] : =  ¢sample¢ creates a 1 X 6 array x[1,6], not an array x[6]. Predefining [x] as one-dimensional x[6] solves the problem (see " array assignments" above).

(3) The command to display as a character string can be given to any array. Values which do not correspond to legal ASCII codes will result in the printing of a blank or of a little square box.

(4) These two are equivalent: [x] : =  ¢abc¢ and [x] : = [( ¢a¢, ¢b¢, ¢c¢)].

(5) When typing from the terminal, the continuation character ,, should not be in a literal. For comparison, enter these two statements:

[x] : =  ¢abc,,
def¢

and

[y] : =  ¢abc
def¢

(Both will have to be selected with the mouse and then entered.) Then display each as a character string.

(6) Displaying a two-dimensional array x[n, m] as a character string produces n rows of length m character strings.

(7) Constants like ¢e¢ have a somewhat ambiguous existence in Fermat - are they scalars or arrays? If you simply enter ¢e¢, the ASCII code for the letter e will be displayed, so it seems to be a scalar. If you enter [x] : =  ¢e¢ you will get a 1 X 1 array. If you enter [w] : =  ¢e¢f  [x]. where [x] is any array, the ¢e¢ will be treated as a scalar, so [w] will be a one dimensional array with first character e
 
 

New Built-in Functions in QFermat
 
 

Divides is implemented in QFermat. Syntax is Divides(n,m) or Divides(n,m,q). Does n divide evenly into m? n and m may not have denominators. As always, 0 means false, 1 means true. If q is present, the quotient will be returned if Divides is true. Various time-saving ideas based on homomorphisms are used.

PDivides(n,m) also determines if n divides evenly into m. The strategy is to substitute each polynomial variable except the highest with a constant. PDivides says true iff these reduced polynomials divide evenly. The constants are chosen with some care. Nonetheless, this is a probabalistic algorithm. An answer of False is always correct, but an answer of True has an infinitesimal probability of being wrong. SDivide(n, m) is similar but faster; see Appendix 4.

Irred, Factor, and Sqfree are described in the Polynomial chapter and in Appendix 4.

Powermod(x,n,m) computes xⁿ mod  m. x must be a polynomial or integer, n must be a positive integer, and m must be a monic polynomial or positive integer. You may omit the third argument if you are in modular mode or polymodding. Note that n often needs to be very large. In modular mode, this is a problem. The solution is that n must be either a constant or must involve only variables that have been created in rational mode while under " Selective Mode Conversion" - see that topic under "The Dangerous Commands." Otherwise QFermat will crash.

Adjoint[x] is the adjoint of matrix [x].

PRoot(x) returns the p^th root of x, when x  is in the ground ring, a field of characteristic p.

Totdeg(x, [a], n) returns the subpolynomial of x of total degree n in the variables listed by [a]. See Appendix Four.

In QFermat GCD(x,y) for polynomials is sophisticated. Hensel's Lemma and the Chinese Remainder Theorem are used. These techniques provide orders of magnitude improvements in speed when x and y are multivariate.

If x is a rational function ("quolynomial") and y a polynomial, then x | y will mod out the numerator and denominator of x by y, then form the quotient. This may not be reasonable mathematically - that's the user's problem.

Two more built-in functions are unique to QFermat. First is Rat(x), which returns 1 if the expression x is rational but not integral, otherwise it returns 0. The second is LCM([x]), which computes the least common multiple of all the denominators in matrix [x]. "Denominators" means those of rational numbers or of expressions like (t² + 3t + 1)/17 or 3/(2t). The denominator of 2/(3+2t) is ignored, since you can't clear it by multiplying [x] by any number (except 0!).

QFermat has an interpreter command &H, with imperative form &(H = ...) to affect the display of constants in modular mode. When turned on, modular numbers will be displayed in the form -m, -m+1,..., 0, 1, ..., n where m = (modulus- 1) div 2 and n = modulus div 2.

Deriv is implemented in QFermat. Deriv(x, t, n) returns the n^th derivative of x with respect to t. x is any (scalar) expression, t is a polynomial variable, and n >= 0.

In QFermat, with the degree-of-polynomial function ¶ there is a difference between ¶x and ¶(x). The second form (where x can be any expression) is slower because it will spend time duplicating x before it computes the degree. All other versions of Fermat will do the "needless" duplication in any event. In the first form, x must be a variable name.

The Characteristic Polynomial

Chpoly computes the characteristic polynomial of a square matrix, in terms of the polynomial variable of highest precedence. The syntax of the command and the method used depend on whether the matrix is sparse or "ordinary."

With the ordinary matrix storage scheme, LaGrange interpolation is used when the matrix consists of all numbers. It is to your advantage to clear the matrix of all numerical denominators before invoking Chpoly. Use LCM to do this. When using the LaGrange interpolation method on integer matrices, Fermat computes the many necessary determinants using the Chinese Remainder Theorem. To do so, it must make an initial estimate of the absolute value of the determinant. The estimate is often rather liberal, resulting in longer than necessary computation time. Now, the determinants in question are simply f(c_i), where f is the characteristic polynomial and { c_i } is a set of certain " sample points." You, the user, may be able to supply a far better bound on |f(c_i)|. For example, you may have some estimate of the location of the roots of f. For this reason, there is a second optional argument to Chpoly in QFermat, a polynomial g such that |f(t)|   <   |g(t)|   for all t. The syntax is Chpoly([x], g).

With sparse matrices, a clever way to compute the characteristic polynomial is the Leverrier-Faddeev method. (See example 4 in Appendix 3 for a short Fermat program that implements the algorithm.) This method often loses to the standard one, det([x] - t I), but it can be faster for matrices that contain quolynomials. The user selects the method with the second argument. If n = 0, Chpoly([x], n) will simply subtract the poly var from the diagonal and invoke determinant. The determinant method will depend on choices made for &K, &L, etc. as described elsewhere. If n <> 0 the Leverrier-Faddeev method is used.

The "modifed Mills method" is a stunningly fast probabalistic "black box" or Wiedeman algorithm that computes the minimal polynomial M(t) of a sparse matrix of integers, or, more precisely, a factor of the minimal polynomial. If one of the roots of M(t) is 0, the associated factor t of M(t) will not show up, but other factors may not show up either. This algorithm is built into Fermat via the command Minpoly. Syntax is Minpoly([a], level, bound). [a] is the matrix, which must be Sparse. level = 0, 1, 2,3, 4 is a switch to tell Minpoly how much effort to expend in its basic strategy. Larger levels take longer, but have a better chance of giving the entire minimal polynomial. bound is an integer at least as big as any coefficient in the minimal polynomial. This argument can be omitted, in which case Fermat will estimate it with the well-known Gershgorn's Theorem.

Repeated calls to Minpoly may return different answers. It may be worthwhile to run it several times and compute the l. c. m. of the answers.

Minpoly is available in rational or modular mode.
 
 

The Dangerous Commands
 
 
 

This chapter describes several commands that should not be used by novices nor the timid. Each one can produce impressive speedups in execution time. Each one, if misapplied, can crash the system, or worse, produce reasonable but erroneous output.

Nothing is free. The cost of speed is loss of flexibility, loss of error checking, and a certain amount of fragility.

Integer

Integer is a directive to Fermat that promises that every scalar quantity encountered during the execution of a function will be a number, not a polynomial. In QFermat, it promises that every quantity is not merely a number, but is an integer of absolute value less than 2²⁸. Beware therefore of overflow, especially when multiplying! Beware of computing something like Ç(50,20), which is larger than 2²⁸

Fermat does not check to see that quantities encountered really are (small) numbers; if one is not a number, the system may crash. Nor is overflow checked for.

Setting Integer in function F has no effect on a function G that F may invoke. If F calls itself recursively, the later invocations will not automatically inherit Integer from the earlier. For it to hold in any invocation, the command Integer must be executed in that invocation.

With Integer, the one character symbolic forms of built-in functions usually execute much faster than the spelled out English word forms (for example, " instead of Cos.)

To further enhance speed, while Integer is in effect some range checking is disabled. If you access an array via x[20] and [x] has less than 20 elements, the ordinary error message will not be generated. What actually will happen is completely unpredictable.

Integer will not catch the error of treating an array parameter as if it were a scalar parameter.

Implicit multiplication doesn't work under Integer. You must write 2*x, not 2x

When Integer is used in a function, the function must use Return to pass back its value.

Ideally, one should modularize a computation so that segments that can take advantage of this feature are moved into their own functions. Very impressive speedups are possible.

Compile

Like Integer, Compile is a directive to Fermat, and applies only to the function in which it is located. Compile implements one of the many functionalities of a genuine compiler. It is a promise to Fermat that the variables given to Compile will not change their identity throughout all the executions of the function. Note carefully: the identity must not change; it is perfectly fine to change the value of the variable. Essentially, " identity" here means " address" or " access path."

To use this feature, the command Compile must be the very first command in the function. (If Integer is to be used also, Integer should be the second command.)

For example, the command Compile(x, [y], p[2], z) promises that x, [y], p[2], and z will always refer to the same variables, everywhere within the function and in every call of the function. The first time the function is executed, when Compile is seen, Fermat looks up the names x, [y], p[2], and z. It then scans the entire function definition and augments every occurrence of these names with a pointer to the symbol table record of the name that exists at that instant. When a reference to one of the variables is encountered during execution, the need to look up the name in the symbol table is therefore obviated. This can provide a very impressive speedup. (Note: there is a speed advantage to compiling p[2], not just [p]. But make sure p[2] has been initialized before the Compile statement is executed. Unfortunately, two dimensional references like p[2,3] cannot be compiled. Such a reference can be converted to one dimensional form using column-major order.)

Only the programmer can tell for sure that it is safe to compile a certain name, say z. If z is a global variable one time the function is executed, but is some other function's local variable another time, it must not be listed. Furthermore, if z is another function's parameter or local variable, it is not the same z during different invocations of that function - recall that the symbol table entry for local variables is removed after invocation ends - and cannot be compiled. This is the most common mistake in using Compile. It is an insidious error.

The list of names following Compile must under no circumstances contain the parameters and local variables of the function. These are automatically compiled, when Compile is specified (but see next paragraph). Rather, the list should be of names created elsewhere, outside the function. If the list (and the parentheses) are omitted, only the parameters and local variables will be compiled.

It is possible not to compile the parameters, by using the syntax Compile*, i.e., placing a * after the word "Compile" and before the parenthesis. Then only the variables in the following list will be compiled. Make sure the list does not contain the name of a parameter! The system may crash if it does. A compiled function cannot be called recursively, unless you use Compile*. (Indeed, that is why Compile* exists as an option.) There will be no error message if you try it. Rather, garbage will be generated or the system will crash. Beware of this mistake when using Compile!

A subtlety arises involving compiled array references and the array initial indexing value, 0 or 1, which can be changed by the command &a. If you place x[5] in the list of variables to be compiled, the address of x[5] will be determined according to the array initial indexing value that exists at the instant of compilation. If you execute &a within the function before a reference to x[5], you have a problem - at the very least, the function will not execute the way it would without compilation.

This problem does not arise from compiled references of the type [x], as in Compile([x]). Later executions of &a are irrelevant. Therefore, a cure to the problem of the previous paragraph is to use Compile([x]) instead of Compile(x[5]) when &a occurs within the function being compiled. (Compile(x[5]) is, of course, slightly more efficient.)

Range checking of compiled array references is disabled. If you access an array via x[20] and [x] has less than 20 elements, the ordinary error message will not be generated. What actually will happen is completely unpredictable.

Compile and Sparse do not work together: Sparse arrays cannot be passed as arguments to a compiled array parameter. Nor can Sparse arrays appear in the list of variables after the Compile command.

Changes to a variable's identity before the compiled function is executed the first time are irrelevant. Compile is an executable statement, and is not noted at function definition time.

There is no point incurring the cost of compilation unless the function will be executed often, or contains a loop that will be executed often. If a function is dominated by polynomial or matrix operations, little relative improvement will result from compilation.

I have seen the use of Integer, Compile, and the increment command eliminate 90% of the execution time of some functions.

Selective Change of Arithmetic Mode

n, for various ns, and then constructing the answer in the field of rationals. For the user to do this well in Fermat it is necessary to drop in and out of modular mode quickly, converting only certain variables. This selective change of arithmetic mode is implemented, for example, with the syntax &(p = 0,19: (x, [y], z)) and &(p = ´: (x, [y], z)). The first command enters modular mode with modulus 19, converting only x, [y], and z. The second returns to rational mode. Other variables are never touched. Attempting to reference other pre-existing variables while in modular mode could lead to a crash. An assignment such as w : = 2 performed while in modular mode will be dangerous if w pre-existed. If not, and if w is added to the list when converting back to rational mode, all will be o.k.

In the example above the variables x, [y], and z will not have their initial values when rational mode is reinstated. For example, 21 would be converted to 2, which will be "converted" back to 2; the 21 will be lost unless the user has gone to the trouble of storing it somewhere else. To preserve the initial value of, say, x and z, use the syntax &(p = 0,19: (*x, [y], *z)) and &(p = ´: (*x, [y], *z)). No variable preserved in this fashion should be changed while in modular mode; think of such variables as read-only input to modular mode. To preserve all, use &(p = 0,19:(**)). Then in returning to rational mode, probably a new variable has been created (otherwise you could preserve results only by writing to a file), say x, so use &(p = ´:(x)). ¹

If you enter modular mode selectively, you must leave it selectively - i.e., with a command of the form &(p = ´:(...)) with something in the parenthses.

If you enter modular mode selectively, you must be cognizant of the fragility of the situation. While in modular mode, do not change Laurent. Do not adjoin or cancel a polynomial variable. Do not change the polymod situation.

Here is a subtlety with selective change of arithmetic mode. The built-in function Chpoly works quickly because certain polynomials are precomputed once and for all, and then used to compute the characteristic polynomial of any matrix given later. The computation of these hidden polynomials takes a fair amount of time, which is why the command to adjoin a polynomial variable, say t, takes longer when it is the first variable adjoined - the polynomials are being computed. When changing arithmetic mode, these polynomials have to be converted to the new format. That takes some time, and would defeat the purpose of quickly dropping in and out of modular mode.

The solution is an interpreter command &B (for Block) that sets a switch to block the conversion of these hidden polynomials. When switched on (the default), the polynomials will not be converted. So if Chpoly is called in modular mode, it will be done the straight-forward but relatively slow way of subtracting t from the diagonal and computing the determinant of the resulting matrix. &B cannot be executed while in modular mode.

Full Interrupt

Many of Fermat's more sophisticated procedures for g.c.d., factoring, etc., affect system-wide global status variables. For example, the g.c.d. routines in rational mode drop in and out of modular mode. Ordinarily, Fermat does not allow Command-interrupt at these times, for the very good reason that the user would be in modular mode, not know it, and may attempt to access some now undefined variable, with bizarre or devastating results. However, some of these procedures can be very time consuming, and the user often wants to interrupt them. The solution is to allow the user to set a special flag that allows interrupts even at these dangerous times. The command is &Z, a simple toggle switch. When switched on, and when &m has been set on, the user can get pretty fast interrupts just about all the time. It is a good idea to immediately do an &g to see if the system globals have changed. And be careful!
 
 

The Interface with C and Pascal
 
 
 

Complex algorithms often have significant parts that don't use any of Fermat's unlimited precision, polynomial or matrix features. Those subtasks could be written in simpler, compiled languages, such as Pascal and C. Assigning such a subtask its own function in Fermat, using Compile, and using Integer, will speed it up tremendously. Nonetheless, the subtask may still dominate the execution time of the entire program.

In these situations the MPW versions of Fermat allows the user to write the subtask in C or Pascal, where its execution will be much faster, and invoke the resulting exterior function from Fermat. For example, the arithmetic of complex numbers is easy to code up in Pascal, and becomes tedious only if there are many variables of complex type involved. Fermat programs investigating Julia sets and the Mandelbrot set could do most of the arithmetic in Pascal, computing and storing points in a file. Then the Fermat function reads the file and displays the points with Fermat's graphics commands. [See illustration below.]

The user must of course be a C or Pascal programmer and must have a fair familiarity with MPW's Linker. The user codes up the subtask as a, let us say, Pascal program, basically as she would write any other Pascal program. She creates a make file that will link together her Pascal file with the Fermat object code files. Building the make file produces a "customized" version of Fermat able to invoke the user's function.

There are several technicalities and protocols. The command in Fermat to invoke the exterior Pascal function is Pfunc; for the C function it is Cfunc. One longint parameter can be passed back and forth. Otherwise Fermat and the function communicate by reading and writing to files.

In the 68XXX architecture, the total size of all global variables must be less than 32K, a bound imposed by the Linker.² The several versions of Fermat use about 15 - 18K, leaving 14 - 17K for the user's variables. The Linker displays (or can be made to display) this number.

I have provided an example of all of this in Appendix 3, example 3. Again, this feature is available only in the MPW versions of Fermat.
 
 

Errors and Warnings
 
 
 

When an error occurs (dividing by 0, array index out of bounds, etc.) execution of the function containing the offending statement stops and the interpreter returns to the lowest interactive level and displays a message describing the error. If the error occurred during or because of a command entered from the terminal, the portion of the command successfully parsed is displayed. Then a function trace dump is displayed: the stack of function calls from most recent to earliest. However this feature is not perfect: if the stack of calls is very deep, only the first 500 (i.e, the 500 most recent) are to be believed.

If the error occurred during a function evaluation, the interpreter then displays two statements, the last one successfully executed followed by the expression it was working on when the interrupt occurred, up to the offending character. Usually the offending character is not displayed, but sometimes it is the last thing displayed. Sometimes the first of these statements (the last one successfully executed) is not displayed, usually because the interrupt occurred in the first line of a function. Keep in mind that the two statements may not be physically adjacent in the function definition.

If a great many function calls have been stacked up when an error occurs, then if the interpreter simply returned to the lowest level, all of those functions' parameters would be on the environment stack. The user would find this extremely inconvenient: his "real" variables would be buried. Thus, the interpreter also pops off all the now useless parameters, up to around 500 calls deep.

The command  &e  can be used to change the way errors are handled in some cases - see the next chapter.

Here is an example of an error that occurred during a function evaluation, taken from a Fermat session. Line numbers [n] added:

 >G2(10)

 *** Inappropriate symbol: '

[1]  error in function definition. ':' or ';' expected.

 *** Error occurred at this point:

[2]    G2(10);

 *** The most recent function evaluations were:

[3]  Undo

     Next

     G2

 *** Function was interrupted at:

[4]  n<0;

     varcount[n] '

The user invoked G2(10). In line [1], Fermat reports the basic problem, apparently a quote in place of a colon in some assignment statement. Line [2] reports the offending place in the input line (read from the terminal) that caused the error. Line [3] begins the trace dump of function evaluation calls.  G2 called Next which called Undo. The error is therefore inside Undo. Line [4] and the following are the last two commands that Fermat executed. n < 0 was executed, and execution terminated at the end of the next line, where we can see the offending quote.

Warnings
 

There are a few situations that produce warning messages, such as entering the command & > from the console. No error occurs, zero is returned, and normal execution proceeds. The more serious of these are called "Fermat Warnings," which occur when some internal condition that "must" be true is in fact false. For example, in the middle of certain G.C.D. computations, some polynomials "must" divide into others. If that fails, a Fermat Warning is issued. Usually the interpreter pops up to a new level.

Please report all "Fermat Warnings" to the author!

A commonly encountered but usually harmless warning message is " end of line found inside a comment." It occurs whenever Fermat reads a function definition and sees a multi-line comment. It is certainly allowable to have multi-line comments. Yet the warning message must be produced in order to catch the error of entirely omitting the end-of-comment symbol ( } ). If you do that, Fermat will be stuck in an infinite loop, and your only notification of that fact will be the above warning message.
 
 

Popping, Pushing, Debugging, and Panic Stops
 
 
 

Let's say that a computation seems to be going on too long, or for some other reason you want to stop it. Just press the command key - Fermat will halt. The interpreter will push on a new command level. You can now examine any of the variables, or do anything else - it's just as good as the original command level (almost - more later). This is a useful debugging feature. To resume the computation (and fall back to the original command level), enter the command &p  ("pop"). To stop the computation completely and fall back to the lowest level (in effect, a panic stop), enter &á

The ability to interrupt with the command key is convenient but can be a bit expensive in terms of execution time. Therefore, the user may disable this feature.

When off, the user's programs will execute most quickly, but the user cannot interrupt them.

When on, the user's program will almost always halt immediately after the command key is pressed. Just keep holding it down until Fermat halts. Excepted from this are some basic arithmetic operators. You will not be able to interrupt something like 6 ^ 99999.

To set the interrupt level, use the command &m, or its imperative form &(m= < level > ). With the first form, you will be prompted to enter a level. Enter 0 or 1, to stand for off or on. As usual, an expression may be entered that evaluates to 0 or 1.

Execution time savings vary. Typically you will save 20% by disabling the interrupt facility. Experiment!

An interrupt subtlety: if in the new command level you delete the function you were executing, when you fall back to that command level, the function will still be able to resume, i.e., it still exists internally until its execution terminates. However, since its name has been erased from the list of functions, no function can call it. Therefore, an attempt at recursion will result in an error.

If you interrupt a Fermat built-in function, do not invoke the SAME built-in at the higher command level.

The commands &p and &á can be inserted in a program, as can &P ("push"), which has the same effect as pressing the command key.

Command levels can be stacked up to level 9.

In command levels beyond the first errors are handled differently. When an error occurs, the interpreter tries to recover from it, prints a warning, and continues the computation, with nonguaranteed results. For example, division by 0 causes Fermat to print a warning and use 0 for the result. Obviously, this may spawn further errors. If no reasonable recovery is apparent, the Fermat interpreter is forced to bomb back to the lowest command level. There is then no way to restart the computation at the point that originally created the interrupt.

To somewhat overcome this annoyance, the command  &e  has been added to Fermat. This sets a flag that changes the way errors are handled during function executions and at higher command levels than the first. If you have previously executed &e, when an error occurs, Fermat will push on a new command level and halt. You can then examine the variables and try to rectify the error yourself, and then pop down to the previous level to allow the computation to continue. (You may find that new errors are being spawned. It may take several  &p  commands to get you down to the previous level, so keep trying.) Note however that this is still less then perfect. For instance, if the problem was a syntax error in a function, it will do no good to cancel the function and enter a patched up version - upon popping down to continue the previous level, Fermat will still be reading the previous version of the function. Furthermore, if your attempted patch is not good enough, Fermat may crash. It is safer, therefore, to do a panic stop instead of trying to resume the lower level.
 
 
 

Initializing with Ferstartup
 
 
 

Just as MPW has Startup and Userstartup files, Fermat has a "Ferstartup" file to provide for user-defined initializations. The last thing Fermat does while booting up is to read this file as if it were an input file and execute the statements it finds there. The file can contain any Fermat statements, including read commands, each followed by a semicolon. It should not have &x. It can be empty. Lines can be commented out by starting them with a semicolon.

Besides mere convenience, the initialization file provides the capability of invoking Fermat from MPW execution files, since it is now possible to direct Fermat without ever typing a command on the keyboard. For example, some program (another Fermat?) could create data for Fermat and store it in a file in Fermat-readable form, then insert the right read command in Ferstartup, then invoke Fermat.

Ferstartup must be at the top level inside the Fermat folder.
 
 

Hints and Observations
 
 
 

These comments resulted from questions by users of Fermat, mostly about the MPW version.

· All versions of Fermat needs a lot of stack space to run. With MPW, use Resedit to reset the MPW Shell resource HEXA 128 (or just plain HEXA in early versions of MPW) to a large value, say 00064000. In versions of MPW before 3.2, the largest possible is 00032000, which is sometimes not enough for Fermat³. The analagous setting for the stand-alone versions can be found as a resource and changed with Resedit. It's called " App Stack."

With system 7, 8, or 9, you will need at least 8 meg of RAM to run QFermat, and probably 12 is a more reasonable lower bound.

· When functions are defined, they are not parsed completely. Only certain key features are looked for, such as comments, loops, and if-statements. Therefore, many syntax errors are not detected until the function is executed.

· The order that the functions are defined in the input file is completely irrelevant, unless you put a statement in the input file that actually invokes one of the functions. Then, of course, that function must have been previously defined, and any function it calls must have been previously defined.

· The order in which most mode-changing commands are given is irrelevant. You may change precision, then set the Laurent flag, then change array indexing - in this order or another.

· Fermat is an interactive language, reading and writing to the terminal, i.e., to the same file. This is somewhat difficult to arrange. Other interactive languages require the user to type a special character, usually a semicolon, at the end of each line. I have used some of those languages and been continually annoyed by the necessity of adding semicolons. Therefore, I arranged it that Fermat does not require them - although they are required at the end of statements in an input file. But nothing is free, and hence the necessity for the continuation character ,, when entering more than one line from the terminal. The only place you may find this annoying is when you are editing a function and you forget to put a ,, at the end of every line. When you attempt to enter the new amended function, Fermat will detect this omission and print out the offending line.

· Functions can contain anything that can be entered from the terminal. In particular, they can contain the definitions of other functions, or even redefine "themselves."

· Very large numbers (those occupying ten or more lines) or polynomials take time to display on the screen. For example, if you enter 790!, more time will be spent computing the digits of this number in base 10 then in actually multiplying it out in the first place. For comparison, compute 790!/789! The timer command &t never counts the time it takes to display something.

· Every time a line is displayed on the terminal, roughly 10 bytes are used up from the memory available to MPW. This is beyond Fermat's control or ability to garbage collect. Therefore, if you leave a program running all night, don't produce hundreds of pages of output.

· The first command that you execute after starting Fermat takes about a half second longer that it otherwise would. This seems to be due to MPW initializing something or loading segments.

· Suppose you want to invert matrix [a]. Typing  [a]^-1  will work, but typing 1/[a] will not. In the first case, the initial "[" alerts Fermat to expect an array expression.

· If you execute a save command &s, and later the machine crashes for some reason, the material that was written to the file may not survive, because the file was never closed. To overcome this, enter &S and name a new file. The original file will be closed.

However, sometimes material is still lost after a bad crash. You can absolutely save something by interrupting Fermat with the command key, displaying the variables you want, grabbing the data with Copy (command-c), and then transferring to some word processor. Open a document there and save the material. Or use the Scrapbook.

· An annoying problem sometimes arises when highlighting several lines of text with the mouse and then hitting enter, as one would do in order to define a function or matrix on the MPW Worksheet during a Fermat session. If an error occurs, due to, for example, not supplying the continuation character, the next attempt to enter anything almost always yields an inappropriate error message. Apparently Fermat (or MPW) misreads the prompt character. This is beyond my control. Hitting enter a few times will straighten it all out.

· In reading and saving, the difference between &r and &R, and &s and &S, may be confusing. Use the small letter form either the first time you give a read or save command, or to continue reading or saving to the same file as you previously specified. Use the capital letter form to change to a different file, thereby closing the former one.

In either case, unless you give the imperative form of the command, Fermat will prompt you for the name of the file. This means that if you put &S in the middle of a function, you cannot leave the terminal until the function has gotten to that point in the program; someone must be there to respond to the prompt. The problem is solved by using the imperative form: &(S = < filename > ). This command can be used whether or not it is the first save. You can simply start up the function and leave it to run overnight.

You may also use a character string stored in an array to name the file, as in &(S = [x]). This allows file names to be computed within functions.

· Fermat has been designed to catch all sorts of errors, but no program can allow for all eventualities. Besides the features described in the chapter "The Dangerous Commands," one thing that will cause Fermat to crash is too much recursion. If a program calls itself several hundred times without terminating, all stack memory will be exhausted.

Another thing that will cause a crash is to execute a change of precision &p after you have pushed to a new command level by pressing the command key. Upon resuming the previous computation, Fermat will probably crash. There are also a few strange things you can do with object parameters to cause a crash, and with Powermod.

There's an old Henny Youngman joke: A man goes to see his doctor. He lifts his arm up and down in an odd fashion, saying "Doc, it hurts when I do this." The doctor says, "Don't do that."
 
 

Appendix 1. Installation Instructions for MPW
 
 
 

Of course it is not neccesary to use MPW, but for those doing so, your MPW folder must be at the top level on the desk top, and must contain a folder called "fermat," which in turn must contain pow2028, char_poly_data, char_poly_index, primeprods, primelist, all the errs files, the executable versions of Fermat, ferstartup, and the folder called " object files", which contains all the Fermat object files. (You need the object files only if you are going to invoke Pascal or C programs, as described in "The Interface with C and Pascal.")

Open your UserStartup file and add the commands:

Alias ferm "{Boot}MPW:fermat:ferm"

Export ferm

for every version of Fermat (here called "ferm") that you have in the fermat folder. This makes it possible to invoke Fermat from inside any folder; i.e., no matter what the default directory is.

On 68XXX machines, Fermat assumes 68020 or better processor with 68881 or better coprocessor.

Fermat needs a lot of stack space to run. Use Resedit to reset the MPW Shell resource HEXA 128 (or just plain HEXA) to 00064000 (00032000 in MPW 3.1 or earlier).

Fermat needs at least four meg of main memory. Eight meg is better, and you will need more if you do lots of graphics, especially movies (FFermat only).

Click on the MPW shell, then do command-I (information). Check the amount of space allocated to MPW. You want it to be as large as feasible, 8 meg anyway. This is especially important if you are using system 7. Otherwise, Fermat may crash when you attempt any serious computation. 12 meg is a more reasonable bare minimum.

Under system 6 (is anyone still using system 6?), make sure you have turned on the RAM cache with the control panel. 256K is enough.

If you have a MacIIci (what a great machine that was!), try turning off the color on your monitor and selecting only two gray scales with the control panel. This should provide a 30% speedup of all MPW tools, because there is no video card in the standard IIci. This is worth a try on other MacIIs as well.

Things will look best on the Worksheet if you choose the font Courier 12.
 
 
 

Please send constructive criticisms to Robert H. Lewis, Department of Mathematics, Fordham University, Bronx NY 10458, rlewis@fordham.edu.

January 20, 2002

(Some of the special symbols don't look quite right here, such as option-=, option-shift-7, and option-t.)
 
 

Symbol	Equivalent	Keystroke	Meaning	Use	Apply to arrays?
\	-	\	divide and truncate	x \y, [x]\y	yes
|	-	shift-\	modulo	x | y, [x] |y	yes
|...|	-	shift-\	absolute value	|(x+y)|	no
?	-	option-=	not equal to	x ? y, x <> y	no
>	-	option- >	greater than or equal	x > y, x >= y	no
<	-	option- <	less than or equal	x < y, x <= y	no
$	-	shift-4	greatest integer	$x	yes
^	-	shift-6	raise to power	x^n	yes
_	-	shift- -	suppress modular	_9087	no
:	-	shift- ;	suppress display; terminal only	10^9999 :	no
'	-	'	start, end quote	' ... '	no
¢	-	option-shift-]	derivative of poly. or quoly	x¢	no
p	-	option-p	pop	&p	no
P	-	option-shift-p	push	&P	no
P	-	option-shift-p	multiply out	P<i = 1,n > [... ]	no
÷	Sqrt	option-v	square root	÷x	no
!	-	shift-1	factorial	x!	no
!	-	shift-1	display	!(x,y,z), !!x, ![z	yes
?	-	shift-/	get from terminal	?s, ?[x]	yes
¿	-	option-shift-/	random number	¿, ¿¿	yes
t	-	option-t	previous value	x :=  t, [t], t[n]	yes
S	-	option-w	add up	S < i = 1,n > [... ]	no
·	-	option-8	polynomial evaluation	x·y , x·(u = y), x·[y]	yes
·	-	option-8	panic stop	&·	no
ç	-	option-c	coefficient of poly.	ç(x,n)	no
£	Lterm	option-3	leading term of poly.	£x	no
¶	Deg	option-d	degree of a polynomial or size of matrix	¶x, ¶[x]	yes
& oe	-	option-q	monomial list (QFermat)	& oe	no
D	Det	option-j	determinant	D[x], others	yes
¬	-	option-shift-v	diagonal of matrix	¬[x]	yes
TM	-	option-2	transpose	TM[x]	yes
¢	Cols	option-4	num. of columns	¢ [x]	yes
f	-	option-f	concatenate arrays	[x] f [y]	yes
±	-	option-shift-=	change Laurent	&±	no
@	-	shift-2	cancel	@[x], @F, @ < [x] >	yes
m	-	option-m	suppress modular	&m, &m(...)	no
¶	-	option-7	add (drop) poly. var.	&¶	no
{	-	shift-[	start comment	{ ... }	no
}	-	shift-]	end comment	{ ... }	no
}	-	shift-]	exit function	& }	no
>	-	shift-.	exit loop	& >	no
]	-	]	cycle (in loop)	&]	no
...	-	option-;	ellipsis	[x[2 ... , 4 ...]	yes
´	-	option-5	rational mode	&(p=´)	no
,,	-	option-shift-w	line continuation	only from terminal	no
§	-	option-6	system function	only from terminal	no
f	-	option-shift-7	array of arrays	f[1] := [x]	yes
o	Codeg	option-k	codegree; in poly.	o x, o (x, 2)	no
Ç	-	option-shift-c	binomial coefficient	Ç(n, r)	no
¥	Var	option-y	poly var, by ordinal	¥1   (= top level poly var)	no

Appendix Three

These examples are not necessarily "best possible" solutions.

Example 1.

Note: To appease the typesetting program that created this document, it was necessary to write &u in the following where &m was intended, # instead of á, Transpose instead of ^TM, f instead of f . Several other less noticeable compromises were made.

;  This set of functions was used in November 1990 by Robert H. Lewis to 
; answer an unsolved question in group theory. A certain property of some 
; groups called the Gkj groups was checked by counting the number of 
; homomorphisms from Gkj -> SL(2, Z/4). It turns out that for certain k and 
; j, this number can vary from the expected number. Gilbert Baumslag 
; presented the question in this form.  Gkj is a one relator group with two 
; generators; the relation involves the two generators raised to powers k 
; and j. To find a homomorphism from Gkj into a group it is necessary and 
; sufficient to find 2 elements of that group that satisfy the same relation. 
; SL(2, Z/4) is picked because it is finite, 2x2 matrices are easy to 
; compute with, and it has an appropriate lower central series.

; Fermat note: Fermat reads the first six statements and does them. Then 
; it reads the function definitions. Then it executes the last 11 statements. 
; There is no reason to have this particular order of statements in this file. 
; In particular, the function definitions can be in any order.

&(t = 1);
&(m = 0);      ;; Turn off the mouse interrupt capability. 

&(¶   = x);   &(¶   = y);   &(¶   = z);   &(¶   = w);         ;; Adjoin 4 polynomial variables.

 
;; Find [p] among the list of 2x2 matrices representing SL(2,Z/4).
 Function Find(i,size,j,ans) =  
    ans := 0; 
    for j = i, size do 
       if p[1] <> e[j] then &]
       fi; 
       if p[2] <> f[j] then &]
       fi; 
       if p[3] <> g[j] then &]
       fi; 
       if p[4] = h[j] then 
          ans := &u(j); 
          & > 
       fi
    od; 
    ans.;

;; Find where each element [p]'s inverse transpose is in SL(2,Z/4).
 Function SetInvTr(i) =  
    for i = 1, invs do 
       [p] := [(e[i],f[i],g[i],h[i])]; 
       { ith element of SL(2,Z/4)  }
       [p] := Transpose[p]; 
       [p] := [(p[4],-p[2],-p[3],p[1])]; 
       { fast inverse  }
       it[i] := Find(1, &u(invs))
    od.;

;; Insert is used by DelConj.
 Function Insert(x1,i,k) =  
    i := 1; 
    while i <= j do 
       if x1 > conjs[i] then 
          for k = j + 1, i + 1, - 1 do 
             conjs[k] := conjs[k-1]
          od; 
          conjs[i] := x1; 
          &u(j :+ ); 
          &}
       else
          if x1 = conjs[i] then 
             &}
          else
             &u(i :+ )
          fi
       fi
    od; 
    conjs[i] := x1; 
    &u(j :+ ).;

; Compute the conjugacy classes in SL(2,Z/4). Look at each element of 
; SL(2,Z/4). Delete from the list of SL(2,Z/4) all of its conjugates.  
; When finished, one representative is left in each class. 
  
 Function DelConj(i,j,k,n,spot) =  
    Array conjs[invs\2]; 
    Array a1[2,2]; 
    Array count[2]; 
    [count] := 0; 
    &u(size := invs); 
    i := 1; 
    while i <= size do 
       [a1] := [(e[i],f[i],g[i],h[i])]; 
       j := 0; 
       { Fill the array conjs to tell where the conjugates of [a1] are. }
       for k = 1, invs do 
          [p] := [(e1[k],f1[k],g1[k],h1[k])]; 
          [q] := [(p[4],-p[2],-p[3],p[1])]; 
          [p] := [q]*[a1]*[p]; 
          spot := Find(i, size); 
          if spot = 0 then 
             !!'spot is 0'
          fi; 
          if spot <> i then 
             Insert(spot)
          fi
       od; 
       !!('number of conjugates is ', &u(j + 1)); 
       [count] := [count] f & u(j + 1); 
       { Delete the conjugates of [a1].  Change [it] accordingly.   }
       for k = 1, j do 
          &u(size:-); 
          for n = 1, conjs[k] - 1 do 
             if it[n] = conjs[k] then 
                it[n] := i
             fi; 
             if it[n] > conjs[k] then 
                &u(it[n] := it[n] - 1)
             fi
          od; 
          for n = conjs[k], size do 
             e[n] := e[n+1]; 
             f[n] := f[n+1]; 
             g[n] := g[n+1]; 
             h[n] := h[n+1]; 
             it[n] := it[n+1]; 
             if it[n] = conjs[k] then 
                it[n] := i
             fi; 
             if it[n] > conjs[k] then 
                &u(it[n]:-)
             fi
          od
       od; 
       &u(i :+ )
    od; 
    [count] := [count[3∞]]; 
    @([conjs], [a1]).;
 
;; OneDet tests every possible 2x2 matrix over Z/4 and saves the 
;; invertible ones.  j1 is owned by CreateData.
 Function OneDet(i,j,k,l) =  
    for i = 0, dim - 1 do 
       for j = 0, dim - 1 do 
          for k = 0, dim - 1 do 
             for l = 0, dim - 1 do 
                if i*l - k*j = 1 then 
                   &u(j1 := j1 + 1); 
                   e[j1] := i; 
                   f[j1] := j; 
                   g[j1] := k; 
                   h[j1] := l
                fi
             od
          od
       od
    od.;

;; Count the number of all 2x2 matrices over Z/4 that have determinant 1.
 Function CountInv(i,j,k,l,j1) =  
    for i = 0, dim - 1 do 
       for j = 0, dim - 1 do 
          for k = 0, dim - 1 do 
             for l = 0, dim - 1 do 
                if i*l - k*j = 1 then 
                   &u(j1 :+ )
                fi
             od
          od
       od
    od; 
    j1.;

 Function CreateData(j1) =  
    Array e[invs]; 
    Array f[invs]; 
    Array g[invs]; 
    Array h[invs]; 
    OneDet; 
    [e1] := [e]; 
    [f1] := [f]; 
    [g1] := [g]; 
    [h1] := [h].;
            
;    Main tries one pair of exponents, exp1 and exp2.  It counts the number 
; of homomorphisms from the group they define into the 2x2 matrix group.
; The basic word that defines the group has been massaged to minimize
; the amount of work within the innermost loop.  That is also the reason 
; for using the polynomials -- part of the word can be evaluated in advance 
; " once and for all".
 Function Main(exp1,exp2,i,j,k,l,yes) =  
    time := &T; 
    hits := 0; 
    for i = 1, size do 
       [c] := [(e[i],f[i],g[i],h[i])]; 
       [ci] := [c]^exp1; 
       [cj] := [c]^exp2; 
       [prod1] := [cj]*[b]; 
       [prod2] := [b]*[cj]; 
       { If [c] is in the center, so that count[i] = 1,
         the answer is easy. Increment and cycle.   }
       if count[i] = 1 then 
          &u(hits := hits + invs); 
          &]
       fi; 
       for j = 1, invs do 
          &(m = 1); 
          { Allow mouse interrupts.   }
          [a] := [(e1[j],f1[j],g1[j],h1[j])]; 
          [p] := [ci]*[a]; 
          &(m = 0); 
          { Turn off mouse interrupts.   }
          [a1] := [(p[4],-p[2],-p[3],p[1])]; 
          [test] := [a1]*[a]*[p]; 
          [test] := [prod1] - [prod2]*[test]; 
          for k = 1, invs do 
             yes := 1; 
             for l = 1, 4 do 
                if test[l]#(h1[k], g1[k], f1[k], e1[k]) <> 0 then 
                   yes := 0; 
                   & > 
                fi
             od; 
             if yes = 1 then 
                &u(hits := hits + count[i])
             fi
          od
       od
    od; 
    !!('time is ', &u((&T - time)/60), 'seconds.'); 
    !!('i,  j,  hits are ', exp1:3, exp2:3, hits).;

; Driver is the 'main program'. It first counts the number of invertible 
; 2x2 matrices over Z/dim, using CountInv. (dim = 4 is the right choice, 
; for which invs is 48). CreateData then actually creates the 48 matrices. 
; Matrices are stored in 4 arrays [e], [f], [g], [h],one for the upper 
; left coordinates, one for the lower left, etc. (This is inelegant and 
; could be done better with Fermat's array of arrays.)  These are 
; duplicated into [e1], [f1], etc. Then SetInvTr finds where each matrix's 
; inverse transpose is. This data is stored in array [it]. DelConj breaks 
; the group into conjugacy classes. (This is done because the basic defining 
; word is invariant under conjugation, so we save some time this way.)  
; When it's finished, [e], [f], [g], [h] store 1 representative from each 
; class, and [it] says where the inverse transpose of each of these 
; representatives is.  But it turns out that it[j] = j for all j.  Had 
; this not been true, a further speed up would have been possible. Then a 
; loop indexed on i tests the exponents from the arrays [try1] and [try2] 
; thereby looking for maps from the group Gkj,where k = try1[i] and j = 
; try2[i]. Main counts the number of homomorphisms from the Gkj group into 
; SL(2,Z/4). 
 Function Driver(size,time,invs) =  
    time := &T; 
    Array a[2,2]; 
    Array b[2,2]; 
    Array c[2,2]; 
    Array ci[2,2]; 
    Array cj[2,2]; 
    Array p[2,2]; 
    Array a1[2,2]; 
    Array b1[2,2]; 
    Array q[2,2]; 
    Array test[2,2]; 
    [b] := [(x,z,y,w)]; 
    invs := CountInv; 
    !!('order of group is ', invs); 
    Array it[invs]; 
    CreateData; 
    SetInvTr; 
    !!('initial [it]'); 
    ![it;   !; 
    DelConj; 
    !!('time is ', &u((&T - time)/60), ' seconds.'); 
    !!('size is ', size); 
    { size = number of conjugacy classes.   }
    !!('Now [it] is '); 
    ![it[1...size];  !; 
    !!'[count] is '; 
    ![count;  !; 
    for i = 1, 100 do 
       Main(try1[i], try2[i])
    od.;
 
?dim;                      ; Ask user for dim.
&(p = dim);                ; Go to modular mode mod dim. In converting, dim
                           ;  becomes 0.  (of course)
dim := -1;                 
&u(dim := dim + 1);        ; dim now holds dim, not 0.
Array try1[100];           ; Create exponents to try for k and j.
Array try2[100];
[try1[1...7]] := 1;
[try1[8...14]] := 2;
&(a = 0);
[try2] := [( &u (i|7 + 2) )];   
&(a = 1);

&x;

Example 2.

These functions compute the determinant det of a matrix [x] with one-variable polynomial entries. This method is now part of the built-in determinant function. The strategy is to figure out a bound on the degree of the determinant det, then evaluate det at a set of points, then use interpolation to construct det.

The hard thing is to estimate in advance the degree of the determinant. Several methods are possible. This one has the advantages of being fairly fast and not often overestimating the degree. A bad overestimate will cost dearly in the final interpolation, here accomplished with the S command. This method also spotlights many of Fermat's features.

We assume that the matrix [x] has already been defined. 'Degree' estimates the degree of det by simulating what would happen to the terms of highest degree in the matrix if the determinant were actually computed with row and column manipulations.

Rational mode is assumed.

 Function Degree(n,i,j,k,diff,answer) =  
    { Create matrix of degrees of entries in [x]. n = dimension of [x].  }
    Array y[n,n]; 
    for i = 1, n*n do 
       if x[i] = 0 then  y[i] := -10^10   { " -infinity" }
          else y[i] := Deg(x[i]) fi
    od; 
    answer := 0; 
    { This loop simulates row and column manipulations in [x] by working
       with the degrees in [y].  }
    for i = 1, n do 
       { Try to factor out as many powers as you can.  }
       Rowscan(n, i); 
       Colscan(n, i); 
       { Find spot with lowest degree >= 0.  }
       j := Det_ + ([y[i...,i...]] + 1); 
       if j = -1 then  & > fi; 
       { Convert linear to (row, column) coordinates in [y].  }
       k := (j - 1)\(n - i + 1) + i; 
       j := (j - 1)|(n - i + 1) + i; 
       { If the remaining matrix is all constants, stop now.  }
       if Maxdegree(n, i) = 0 then & >  fi; 
       { Simulate row and column manipulations.  }
       Switchrow([y[,i...]], i, j); 
       Switchcol([y[i...,]], i, k); 
       for j = i + 1, n do 
          diff := y[j,i] - y[i,i]; 
          for k = i + 1, n do 
             [y[j,k]] := Max(y[j,k], diff + y[i,k])
          od
       od; 
       answer := answer + y[i,i]
    od; 
    answer.;
					
 Function Rowscan(n,i,r,j,deg) =  
    for r = i, n do 
       { Find position and value of smallest nonnegative entry in rth row
              in columns >= i.  }
       j := Det_ + ([y[r,i...]] + 1); 
       deg := y[r,j+i-1]; 
       { If deg is positive, factor it out, changing answer accordingly.  }
       if deg > 0 then 
          [y[r,i...]] := [y[r,i...]] - deg; 
          answer := answer + deg
       fi
    od.;
						
 Function Colscan(n,i,r,j,deg) =  
    for r = i, n do 
       { Find position and value of smallest nonnegative entry in rth column
              in rows >= i.  }
       j := Det_ + ([y[i...,r]] + 1); 
       deg := y[j+i-1,r]; 
       { If deg is positive, factor it out, changing answer accordingly.  }
       if deg > 0 then 
          [y[i...,r]] := [y[i...,r]] - deg; 
          answer := answer + deg
       fi
    od.;
					
 Function Maxdegree(n,i,r,s,deg) =  
    deg := 0; 
    for r = i, n do 
       for s = i, n do 
          if y[r,s] > deg then 
             deg := y[r,s]; 
             &> 
          fi
       od; 
       if deg > 0 then &>  fi
    od; 
    deg.;
		 
 Function Max(a,b) =  if a > b then a else b fi.;		 

{ Assume that t is the polynomial variable. }   

n := Degree(Cols[x]) \2;

{ Standard LaGrange interpolation at points -n, -n+1, ..., n-1, n: }

f := P [ t-i ];

det := S [Det([x]áj)*f/(t-j) / P (j-i) /P (j-i) ]

Example 3.

Included here are

• an external C program that could be invoked from within Fermat.

• a semantically equivalent external Pascal program that could be invoked from within Fermat.

• a make file to compile and link the two programs and create a customized version of Fermat, ferext, to invoke them.

I make no attempt here to explain the syntax of C, Pascal, or MPW.

First, the C program, in a file called ferc.c:

#include <string.h>
#include <stdio.h>
#include <stdlib.h>

##pragma segment [segext]

Cfunc(m)                                  { This name MUST be Cfunc. }
int *m;
{
 int n;
	
 char infname[12] = "instuff";
 char outfname[12] = "outstuff";	
 FILE *infile;
 FILE *outfile;	
 infile = fopen(infname, "r");
 outfile = fopen(outfname, "a");
 fscanf(infile, "%d", &n);
 fprintf(outfile,":x = %d %d \n", n, m[0]);
 m[0] = n+1;	
 fclose(outfile);
 fclose(infile);	
} 

Next, the Pascal program, in a file called ferp.p:

unit ferp_ext;  
															
interface

{$S[segext]}		
{$R-}

uses paslibintf, sane;
				
procedure Pfunc(var m: longint);          { This name MUST be Pfunc. }
	
implementation
	
procedure Pfunc;
   var n: integer;
       infname, outfname: text;
   begin
      open(infname, 'instuff');
      open(outfname, 'outstuff');
      rewrite(outfname);
      readln(infname, n);
      writeln(outfname, ':x = ', n, m);   { a Fermat readable statement }
      m := n + 1;                         { pass back n+1 to Fermat }
      close(infname);
      close(outfname)		
   end;
	
end.

Next, the make file, called ferext.make, that compiles the above programs and links them with the Fermat object files. This assumes that the default directory is ::mpw:fermat:, the mpw folder is at the top level on the desk top, the mpw folder contains all the right interface and library files, and that the object files for the above Pascal and C programs are in the top level inside the fermat folder:

# File: ferext.make

# Target: ferext

# Sources: ferc.c ferp.p

ferp.p.o ˘- ferext.make ferp.p

Pascal -MC68020 -MC68881 -D elems881=true -p -r ferp.p

ferc.c.o ˘- ferext.make ferc.c

C -MC68020 -MC68881 -D elems881=true -p -r ferc.c

SOURCES = ferc.c ferp.p

OBJECTS = ¶
::Libraries:Libraries:Interface.o œ
::Libraries:Libraries:Runtime.o œ
::Libraries:Libraries:ToolLibs.o œ
::PLibraries:PasLib.o ¶
::PLibraries:SANELib881.o ¶
":object files:ferm114.p.o" ¶
":object files:fer414.p.o" ¶
":object files:fer314.p.o" ¶
":object files:fer214.p.o" ¶
ferp.p.o ¶
ferc.c.o

ferext ˘-˘- ferext.make OBJECTS

Link -ss 200000 -p -w -t MPST -c 'MPS ' ¶
OBJECTS ¶
"{Libraries}" stubs.o ¶
"{CLibraries}"CRuntime.o ¶
"{Libraries}"Runtime.o ¶
"{CLibraries}" StdCLib.o ¶
"{CLibraries}"CSANELib.o ¶
"{CLibraries}"Math.o ¶
"{CLibraries}"CInterface.o ¶
"{CLibraries}"Complex.o ¶
"{Libraries}" Interface.o ¶
"{PLibraries}"PasLib.o ¶
"{PLibraries}" SANELib.o ¶
"{Libraries}"ToolLibs.o ¶
-o ferext

Example 4.

The Leverrier-Faddeev method for computing matrix inverses and characteristic polynomials (see "The College Mathematics Journal," May 1992, p. 196.) I present two slightly different programs for each case.

 Function CharPoly(a,n,i,pm) =  
    { Assume t is a polynomial variable.  }
    n := Cols[a]; 
    Array c[n]; 
    [q] := [a]; 
    c[1] := -Trace[a]; 
    for i = 2, n do 
       [q] := [q] + c[i-1]*[1]; 
       [q] := [a]*[q]; 
       c[i] := -Trace[q]/i
    od; 
    { Correct for odd n.  }
    if n|2 then 
       pm :=  - 1
    else
       pm := 1
    fi.;

 Function Invert(a,q,n,i) =  
    { Set [q] = the inverse of [a]  }
    n := Cols[a]; 
    Array c[n]; 
    [p] := [a]; 
    c[1] := -Trace[a]; 
    for i = 2, n do 
       [q] := [p] + c[i-1]*[1]; 
       [p] := [a]*[q]; 
       c[i] := -Trace[p]/i
    od; 
    [q] := [q]/(-c[n]); 
    @([p], [c]).;

If appropriate, the inclusion of Integer in either procedure will provide a very dramatic increase in speed. However, Integer cannot be placed in the first procedure as it stands, because of the presence of the polynomial variable t. Rather, all but the the last line could be placed in a seperate procedure that may contain Integer.

Appendix Four

Matrix and polynomial computations are the heart of Fermat. Here is a concise annotated list of relevant features and commands. Most are explained further earlier in the manual.

Polynomials

Recall that most built-in functions allow call-by-reference for parameters. For example, time and compare Terms(x) and Terms(^x) for a large x.

Adjoin polynomial variable: use the command &¶ (Fermat interrogates user for name), or &(¶  = t) (imperative to adjoin t).

Cancel (delete) polynomial variable: similar to above, &(¶   = - t).

Quotient rings and fields: One may choose to mod out by some of the polynomial variables to create quotient rings (or fields). The chapter on "Polymods" describes how. In principle, any monic polynomial may be modded out, say tⁿ + c₁t^n-1 + ....

However, QFermat is best at the case where the quotient ring becomes a field (note well, Q[t, u, ...]  / < p, q, ... > is a field, Z[t, u, ...]  / < p, q, ... > is not). Specifically, suppose the polynomial variables have been attached in the temporal sequence t, u, v, .... Begin by modding out a monic irreducible polynomial p(t) such that F₁ = Z[t]/ < p > is an integral domain, and its field of fractions is Q[t]/ < p > . Then, if desired, mod out by a monic polynomial q(u, t) such that F₂ = F₁[u]/ < q > is an integral domain, and continue in this manner always creating an integral domain, and, by the same stroke, its field of fractions.

You must tell QFermat that a field will result, and it is your responsibility to check this. Do this at each step by adding a comma and 1, as &(P = tⁿ + c₁t^n-1 + ..., 1). You may append a list of primes q such that the modder p(t) remains irreducible mod q. For example, &(P = t³ + t² + t + 2, 1:151, 167, 191, 839, 859, 863, 907, 911, 991). If you omit the list QFermat will compute it for you. (This takes a while, so if you save the session to a file, QFermat will include the list in the saved file and just reload it next time.) QFermat then computes a second list of auxiliary primes: modulo these primes the modding polynomial has a root. Both types of primes are used to speed up g.c.d. computations. If Fermat cannot find enough of either type, it will tell you and instruct you how to get more, using the commands &A and &B.

In all of the above the base field is Q. To make a finite field, first do the above then go into modular mode with &(p = some prime n), where the prime n and the polynomials p, q, ... have been chosen so that Z/n[t, u, ...]  / < p, q, ... > is a field.

Laurent Polynomials: a polynomial with negative exponents. To allow this, activate the toggle switch &±. All of the variables you have created up to that point will be converted so that no negative exponents are in the denominator of a quolynomial and all positives are factored out and moved to the numerator. For example, 1/(t²+2t) will become t^-1/(t+2).

Basic arithmetic: +,  -,  *,  /,  |,  \, ^. p/q creates a rational function ("quolynomial") unless q divides evenly into p. In any event, GCD(p,q) is computed and divided into p and q. p|q means p mod q and p\q means p div q, i.e. divide and truncate.

Mod and div may be "pseudo" results: if c is the leading coefficient of q and is not invertible, there exist polynomials r and y such that c^k p = y q + r, where 0 £ k £ deg(p) + deg(q). QFermat chooses the smallest possible k. p|q is r and p \q is y.

Remquot = remainder and quotient. Syntax is Remquot(x,y,q). q gets the quotient of dividing x by y and the function returns the remainder. Almost twice as fast as calling mod and div separately. Pseudo-remainder and quotient will be returned if necessary.

¶ = Deg. degree of a polynomial (or quolynomial). There are three variants: (1) Deg(x) computes the highest exponent in x (any expression) of the highest precedence polynomial variable. (2) ¶(x, i) computes the highest exponent in x of the i^th polynomial variable, where the highest level variable has the ordinal 1. (3) ¶(x, t) computes the highest exponent in x of the polynomial variable t. In modular mode, Deg returns an actual integer, not reduced modulo the modulus. For a quolynomial, it returns the degree of the numerator.

^o = Codeg; "codegree," just like Deg except it computes the lowest exponent.

á = polynomial evaluation. xáy replaces the highest precedence variable everywhere in x with y. xá(u = y) replaces the variable u with y. x and y could be quolynomials.

There is a fast shortcut form of evaluation called total evaluation. To evaluate x at every variable, use the syntax xá(v₁, v₂, ...), where all the v_i are numbers. There must be a number corresponding to each polynomial variable, in the precedence order - highest precedence (last attached) first.

Fermat also allows evaluation of polynomials at a square matrix. The syntax is xá[y]. The highest precedence polynomial variable in x is replaced with the matrix [y] and the resulting expression simplified. [y] can contain entries that are quolynomials.

In QFermat the following much more general syntax is allowed. Suppose there are five poly vars, e, d, c, b, a in that order (e last and highest). Then qá(d = w, x, y) will replace each d in q with w, each c with x, each b with y. e and a are untouched. Further, w, x, and y can be arbitrary quolynomials.

Numb = is the argument a number (as opposed to a polynomial or quolynomial)? If so the result is 1, else it's 0.

Numer = numerator of a quolynomial. In rational mode, also gives the numerator of a rational number.

Denom = denominator of a quolynomial. In rational mode, also gives the denominator of a rational number.

ç = coefficient in a polynomial (or quolynomial).

(1) Suppose first that only one polynomial variable t has been adjoined. Then the syntax of use is either ç(x) or ç(x, n). (ç is option-c on the keyboard.) x can be any expression. n, the desired exponent, must be a number. In the first form, without n, the leading coefficient is computed. ç(x, n) returns the coefficient of tⁿ in the polynomial x. If x is a quolynomial, the denominator is ignored.

To replace a coefficient, use Rcoef(x, n) : = y; the coefficient of tⁿ in x will be set to the expression y. y must be a number.

If there are several polynomial variables, the coefficient desired is specified by listing the exponents of the variables in precedence order, such as c(x, 1, 2).

In Rcoef(x, ...) = y, x must be a polynomial and ... must be suitable for y.

(2) If t is any polynomial variable, ç(x, t, n) computes the coefficient in x of tⁿ, as if t were the highest level variable. In other words, if x were written out as a sum of monomial terms, find all the terms containing exactly tⁿ and factor out the tⁿ. x must be a variable name, either a scalar name or an array reference x[i]. This form cannot be used on the left of an assignment. Especially useful for n = 0.

Killden(x) sets the denominator of x to 1 - it actually changes x.

Lterm(x) = leading term of polynomial x.

Lcoef(x) = leading numerical coefficient of polynomial x.

Flcoef(x) = leading field-element coefficient of polynomial x, when a quotient field has been created.

Lmon(z) = leading monomial of z. Lmon has an optional second argument. First, Lmon(z) returns the leading monomial of z. This is always an authenic monomial. If you think of a multivariate polynomial in nested (recursive) form, Lmon recursively finds the first term in each level and throws away all the other terms. The result is clearly a monomial. Lmon(z, x), where x is a poly var, stops the recursion at the level of x. For example, suppose u is the higher variable and t the lower variable. Let z = (t² + 3t + 5)u² + 5t*u + 7t - 2. Then Lmon(z) is t²*u² but Lmon(z, t) is (t² + 3t + 5)u². This allows one to have lower variables of a different " status," as if the lower variables are "parameters" and the upper ones are the "real" variables.

Mcoef(x, m1, m2) = monomially-oriented coefficient. m1 (and m2 if present) must evaluate to a monomial; their numerical coefficient is irrelevant. Factor out m1 from all the terms in x that contain it exactly, and return the factor. m2 (optional) specifies variables that are to occur with exponent 0 (as they cannot be included in m1!).

Mfact(x, m) = monomially-oriented factor. m must evaluate to a monomial; its numerical coefficient is irrelevant. Factor out m from all the terms in x that it divides into and return the factor. m may not contain any negative exponents.

Mons(x, [a]) dumps the monomials of x into a linear array [a]. If you want each monomial stripped of its numerical coefficient, use Mons(x, [a], 1).

PRoot(x) returns the p^th root of x, when x ë the ground ring, a field of characteristic p.

Terms(x) = if x were written out as a sum of monomial terms, the number of such terms. x must be a variable name, either a scalar name or an array reference x[i].

&oe: Display each polynomial as a list of monomials.

&c: Enable full Hensel checking. This one is quite technical. If this flag is on (the default) Fermat will double-check the results of certain Hensel Lemma GCD computations. Leaving it off will slightly speed up GCD but introduce an extremely minute probability of GCD giving the wrong answer. See &O next below.

&O: Toggle switch to disable the Hensel and Chinese Remainder Theorem (CRT) methods for polynomial gcd. This is a good idea only when you are working over Z_p for small primes, say p < 30, and the degrees in each variable are fairly small. For such small primes, the Hensel and CRT methods often fail, for "lack of room."

^¢ = derivative of a polynomial (or quolynomial) with respect to the highest precedence variable (the last attached), as in x^¢ (see also Deriv below).

GCD = greatest common divisor, as in GCD(x, y), x and y can be numbers or polynomials, but not quolynomials. If numbers, the result is always positive, except that GCD(0, 0) = 0. GCD(0, x) = |x| if x is a number not 0, and is 1 if x is a polynomial. If they are both polynomials, the result always has positive leading coefficient. If in rational mode, not complex, the result is normalized to have all coefficients integral and be of content 1. In cases where the ground ring is a field, the result has leading coefficient 1.

Content = content of a polynomial; i.e., the GCD of all its coefficients.

Numcon = numerical content, the GCD of all its numerical coefficients.

¯Y = Var. Followed by an expression that evaluates to a positive integer, as in Var(i) or ¯Y-i, returns the i^th polynomial variable, counting the highest (last created) as 1.

Height = the difference between the levels (ordinals) of the polynomial variables in an expression.

Level = the ordinal position of the highest precedence polynomial variable in an expression.

Raise = Two forms: Raise(x) and Raise(x, i). In the first, replace each polynomial variable with the variable one level higher. The second form allows the user to provide an expression i that evaluates to a positive integer, and raises x that many levels, if possible.

Lower = The inverse of Raise. See above.

Divides(n,m) = does n divide evenly into m?

PDivides(n, m) = does n divide evenly into m? When n and m are multivariable polynomials, this procedure attempts to answer quickly by substituting each polynomial variable except the highest with a constant. PDivides says true iff these reduced polynomials divide evenly. The constants are chosen with care. Nonetheless, this is a probabalistic algorithm. An answer of False is always correct, but an answer of True has an infinitesimal probability of being wrong.

SDivide(n,m) = does n divide evenly into m? " S" stands for " space-saving". To save space m is cannibalized. If n does divide m, m becomes the quotient; if not, m becomes 0. m must be a variable name, not an expression. Not probabalistic. Use when you are virtually certain that n divides m and you want the quotient in the fastest way.

Powermod(x, n, m) computes xⁿ  mod  m. x must be a polynomial or integer, n must be a positive integer, and m must be a monic polynomial or positive integer. You may omit the third argument if you are in modular mode or polymodding. Note that n often needs to be very large. In modular mode, this is a problem. The solution is that n must be either a constant or must involve only variables that have been created in rational mode while under " Selective Mode Conversion."

Deriv(x, t, n) returns the n^th derivative of x with respect to t, where t is one of the polynomial variables.

Totdeg(x, [a], n) returns the subpolynomial of x of total degree n in the variables listed in [a]. [a] is an existing array. Each entry should be a single polynomial variable, in no particular order. n is an integer. x is a polynomial; laurent is ok. &oe  is useful.

Factoring Polynomials: QFermat allows the factoring of monic one-variable polynomials over any finite field. The finite field is created by simply being in modular mode over a prime modulus, or by additionally modding out by irreducible polynomials to form a more complex finite field, as described in the section "Polymods." Factoring into irreducibles or square-free polynomials is possible, or polynomials can just be checked for irreducibility.

Factor( <poly> , [x]) or Factor( <poly> , [x], <level> ). The factors of <poly> will be deposited into an array [x] having two columns and as many rows as necessary. (The number of factors (rows) is returned as the value of Factor.) In each row, the first entry is an irreducible polynomial p(t) and the second is the largest exponent e such that p(t)^e divides <poly> . In the second form, <level> specifies the subfield to factor over. Examples are given earlier in this manual. It is best for factoring to use as many variables t, u, ... as possible in creating the field.

Sqfree is similar to Factor except it produces factors that are square-free only.

Sqfree works for any number of variables and over Q. Also, it works recursively by first extracting the content of its argument and factoring it. Over quotient fields, the product of all the factors in the answer may differ from the argument by an invertible factor.

Irred tells if its argument is irreducible, and, if not, describes the factorization. The syntax is Irred( <poly> ) or Irred( <poly> , <level> ) ("level" is explained above). The value returned is as follows:

     -1 means can't decide (too many variables, for instance).

     0 means the argument is a number or a field element.

     1 means irreducible.

     n > 1 means the argument is the product of n distinct irreducibles of the same degree.

     x, a poly, means x is a factor of the argument (which is therefore not irreducible).

Fermat uses the algorithms of H. Cohen, " A Course in Computational Number Theory," Springer Verlag, 1993, p. 123-130. There is some randomness in the algorithms, so the time it takes to factor can vary.

Matrices

Creation: An n X m matrix is created with the command   Array  x[n, m]. Access elements in such an array via x[i,j] or via x[k], which returns the k^th element in column-major order. To refer to an entire matrix, use the syntax [x].

Sparse Matrices: Sparse matrices are implemented in Fermat. This is an alternative mode of storing the data of the array. In an " ordinary" n X m matrix, nm adjacent spots in memory are allocated. If an array consists of mostly 0's, this is wasteful of space. In a Sparse implementation, only the non-zero entries are stored in a list structure.

A Sparse matrix is created by following the creation command with the keyword "Sparse," as in  Array  x[5, 5] Sparse. The only size limitation in QFermat is that the number of rows and the number of columns must each be less than 2²⁸. An array [x] already created can be converted to Sparse format with the command Sparse [x]. There is no requirement that [x] actually have a certain number of zeros.

Indexing: One has a choice of how to index the first element of an array. The default in Fermat is x[1]. This can be changed by entering the command &a, which switches the initial array index to 0. Entering &a again switches back to 1. Note that this is not a property of any particular array, but of how all arrays are indexed.

Dynamic Allocation of Arrays: Arrays that are no longer needed can be freed to provide space for new arrays. This is done with the cancel command, whose syntax is @[x], or, to free several, @([x], [y], [z]).

Arithmetic: Most of the ordinary arithmetic built-in functions can be applied to arrays. See Appendix 2, last column. For example, [x] + [y] is the sum. 2*[a], or [a]*2, multiplies every component of [a] by 2. [a] + 3 adds 3 to every component of [a], and so forth. [a] : = [1] sets an already existing square matrix [a] equal to the identity. [a] : = 1 sets every entry to 1. [z] : = [x] * [y] is the product. [z] : = 1 / [y] is the inverse. Matrix exponentiation (including inverse) is just like scalars (but see Altpower below), such as [z] : = [x]^n.

Arithmetical Expressions: Like numerical expressions, such as [z] : = [a]*([x] + [y] - [1]).

Parameters: Matrices may be parameters in functions.

Subarrays: Fermat allows subarray expressions. That is, part of an array [c] can be assigned part of an array [a]. For example,

In subarray assignments, a vector declared to be one-dimensional (like a[5]) is treated as a column vector, i.e., a[5,1].

Subarray cannot be used with Sparse matrices, but Minors can (see below).

Matrix Built-in Functions:

D, or Det, is used in several ways to compute a scalar from an array argument. If used by itself on a square matrix, D is determinant. D#([x] = a) returns the number of entries in [x] that equal a. Similarly D#([x] > a) and D#([x] < a) compute the number of entries of [x] larger or smaller than a. If any entry is a polynomial, an error results. D ^ô[x] returns the index of the largest element of [x] (in column major order if [x] is a matrix). D_[x] returns the index of the smallest element of [x]. D_ + [x] returns the index of the smallest nonzero element of [x], or -1 if there is no such element.

Determinant: Fermat uses three basic methods to compute determinant: expansion by minors, Gaussian elimination, and reducing modulo n for some n's. The last of these is used for matrices of integers or polynomials with integer coefficients. The actual determinant can be reconstructed from its values modulo n (for a "good" set of n's) by the Chinese Remainder Theorem (see Knuth volume 2). Alternatively, it is often possible to work modulo an easily computed "pseudo determinant" known to bound the actual determinant. Gaussian elimination is applicable in all situations - all one needs is the ability to invert any nonzero element in a matrix. If the matrix is small enough, expansion by minors is faster (see &D.) Gaussian elimination can be nontrivial and even problematical in modular arithmetic over a nonprime modulus, in polynomial rings, and in polynomial rings modulo a polynomial. Fermat has heuristics to guide its choice of method.

When the matrix has all polynomial entries, QFermat has two other methods. It may also compute the determinant with "Lagrangian interpolation:" constants are substituted for the highest polynomial variable everywhere in the matrix, and so on recursively. The algorithm is probabalistic, with very high probability of success. It is very fast, especially for two or more variable polynomials.

Nonetheless, if there are many polynomial variables and the matrix is sparse or has a regular pattern of zeros, expansion by minors can be by far the fastest method. Setting the determinant cutoff (with &D) at least as large as the number of rows will force Fermat to do this method. This is true of Sparse matrices as well as " ordinary."

Secondly, QFermat uses the well-known Gauss-Bareiss method (for a matrix of all polynomial entries).

QFermat has heuristics to choose among the methods, but the user may override them and force a particular method. Assuming an m X m matrix of all polynomial entries, if m is more than three and the user has left &D = -1, the default method is Lagrangian interpolation, unless the "mass" of the matrix is very small. The "mass" is estimated by a heuristic and compared to a cutoff. The user can change the cutoff with the command &L. The default is 5000. Therefore, to turn off Lagrangian interpolation, give a very large value (up to 2³¹-1). To then choose Gauss-Bareiss, set the command &K = 1. This is a bit confusing, so summary:

To force Gauss-Bareiss, set &K = 1 and &D = 2.

To force Gaussian elimination, set &K = 0 and &D = any d > 0. At the d X d stage, it will switch to expansion by minors.

If m Ò 4, to force Lagrangian interpolation, set &D = -1 and &L = 1.

LCM = the least common multiple of all the denominators in a matrix. "Denominators" means those of rational numbers or of expressions like (t² + 3t + 1)/17 or 3/(2t). Use this to clear a matrix of its integer denominators. The denominator of 2/(3+2t) is ignored, since you can't clear it by multiplying [x] by any number.

Adjoint = adjoint of a square matrix.

Chpoly = the characteristic polynomial of a square matrix. The syntax of the command and the method used depend on whether the matrix is sparse or "ordinary."

With the ordinary matrix storage scheme, LaGrange interpolation is used when the matrix consists of all numbers. It is to your advantage to clear the matrix of all numerical denominators before invoking Chpoly. When using the LaGrange interpolation method on integer matrices, Fermat computes the many necessary determinants using the Chinese Remainder Theorem. To do so, it must make an initial estimate of the absolute value of the determinant. The estimate is often rather liberal. The determinants in question are simply f(c_i), where f is the characteristic polynomial and { c_i } is a set of " sample points." You, the user, may be able to supply a far better bound on |f(c_i)|. For example, you may have some estimate of the location of the roots of f. For this reason, there is a second optional argument to Chpoly in QFermat, a polynomial g such that |f(t)|£ |g(t)| for all t. The syntax is Chpoly([x], g).

With sparse matrices, a clever way to compute the characteristic polynomial is the Leverrier-Faddeev method. (See example 4 in Appendix 3 for a short Fermat program that implements the algorithm.) This method often loses to the standard one, det([x] - t I), but it can be faster for matrices that contain quolynomials. The user selects the method with the second argument. If n = 0, Chpoly([x], n) will simply subtract the poly var from the diagonal and invoke determinant. The determinant method will depend on choices made for &K, &L, etc. as described elsewhere. If n <> 0 the Leverrier-Faddeev method is used.

Minpoly: The "modifed Mills method," a fast probabalistic "black box" or Wiedeman algorithm that computes the minimal polynomial M(t) of a sparse matrix of integers, or, more precisely, a factor of the minimal polynomial. If one of the roots of M(t) is 0, the associated factor t of M(t) will not show up, but other factors may not show up either. This algorithm is built into Fermat via the command Minpoly. Syntax of use is Minpoly([a], level, bound). [a] is the matrix, which must be Sparse. level = 0, 1, 2,3, 4 is a switch to tell Minpoly how much effort to expend in its basic strategy. Larger levels will take longer, but have a better chance of giving the entire minimal polynomial. bound is an integer at least as big as any coefficient in the minimal polynomial. This argument can be omitted, in which case Fermat will supply an estimate based on the well-known Gershgorn's Theorem.

Repeated calls to Minpoly may return different answers. It may be worthwhile to run it several times and compute the l. c. m. of the answers.

Sumup = add up the elements of an array.

Trace = trace of a matrix.

Altmult. Multiply two matrices using the algorithm of Knuth volume II, p. 481. A big time saver when multiplication in the ring is much slower than addition. Especially good for Polymods (see that chapter). Syntax is Altmult([x],[y]).

Altpower. Uses Altmult to take a matrix [x] to the power n. Syntax is Altpower([x],n).

MPowermod([x], n, m) computes [x]^n mod m, analagously to Powermod. [x] contains only polynomials or integers, n must be a positive integer, and m must be a monic polynomial or positive integer. You may omit the third argument if you are in modular mode or polymodding. Note that n often needs to be very large. In modular mode, this is a problem. The solution is that n must be either a constant or must involve only variables that have been created in rational mode while under "Selective Mode Conversion."
To see the effect of this command, try creating a 20 X 20 matrix of random integers, let n = 100 and m = 10^8. Compute [x]^n mod m both ways.

^TM = transpose matrix, as in [y] : = ^TM[x].

¬ refers to the diagonal of a matrix, as in ¬[y]: = [x]. [x] is considered a linear array. The diagonal of [y] becomes the entries of [x]. If the name [y] does not yet exist, a new square matrix will be created with off-diagonal entries 0. If square matrix [y] of the right size (i.e., rows equal to the number of entries of [x]) does exist then the off-diagonal elements are not changed.

Dually, ¬ can be used on the right side of an assignment, as in [y] : = ¬[x], which sets [y] equal to a linear array consisting of the diagonal elements of [x]. [x] does not have to be square.

To create a diagonal matrix with all entries equal to a constant, say 1, you can use the easier form  [x] : = [1], if [x] already exists as a square matrix.

¢[x] = Cols[x] = number of columns of array [x].

œ = Deg = number of elements in an array. œ[x] = total size of array [x] (rows X columns).

[t] = last computed array. Fermat has a hidden system array. If you type the command  [x] + [y], arrays [x] and [y] will be added and, since you didn't provide an assignment of the result, the result will go into the system array. You can later access it by typing, for example, [z]: = [z] + [t]. Subarrays cannot be used with [t]. The symbol is option-t on the keyboard.

f = concatenate arrays; glue two arrays together to form a larger one, as in [z] : = [x] f [y]. Neither array can be Sparse.

Iszero = is the argument (an array) entirely 0? If so, return 1, else return 0. Syntax: Iszero[x].

Switchrow = Interchange two rows in an array. Syntax: Switchrow([x], n, m).

Switchcol = Interchange two columns in an array. Syntax: Switchcol([x], n, m).

Normalize = convert to a diagonal matrix. The matrix must not be Sparse. If requested, Fermat will return the change of basis matrices used in normalizing. Possible invocations include Normalize([x]) and Normalize([x], [a], [b], [c], [d]). In the second case, matrices [a], [b], [c], and [d] will be returned that satisfy [a]*[x¢]*[b] = [x], where [x¢] is the original [x], and where [c] = [a]^-1 and [d] = [b]^-1. The value returned by Normalize is the rank of [x].

You can omit any of the change of basis matrices. For example, Normalize([x], ,[b], , [d]) and Normalize([x], [a], , [c]). Every comma promises that an argument will eventually follow.

Colreduce = Column reduce a matrix. The matrix may NOT be Sparse. By column manipulations, the argument is converted to a lower triangular matrix. If requested, Fermat will also return the change of basis (or conversion) matrices that it used in normalizing. Possible invocations include Colreduce([x]) and Colreduce([x], [a], [b], [c], [d]). In the second case, matrices [a], [b], [c], and [d] will be returned that satisfy [a]*[x¢]*[b] = [x], where [x¢] is the original [x], and where [c] = [a]^-1 and [d] = [b]^-1. The value returned by Colreduce is the rank of [x]. As with Normalize, you can omit any of the change of basis matrices.

Colreduce cannot be used on sparse arrays. In addition a function Pseudet is implemented. Pseudet([x]) computes (indirectly) a "pseudo-determinant," a nonzero determinant of a maximal rank submatrix. It returns the rank of the matrix and leaves the matrix in diagonal form (so [x] is changed). The product of the diagonal entries is (up to sign) the "pseudo-determinant." The optional form Pseudet([x], [rc]) returns a 2 X rank[x] matrix [rc] specifying the rows (first row of [rc]) and columns (second row of [rc]) that constitute the maximal rank submatrix. A second optional form is Pseudet([x], [rc], k) or Pseudet([x], , k). Integer k specifies the last row/col to pivot on. The entries beyond spot [k,k] are left.

Rowreduce = Row reduce a matrix. The matrix must be Sparse. Exactly like Colreduce but for sparse arrays and row reduction.

Smith = Put a matrix of integers into Smith normal form. The matrix may be Sparse. This function can only be used in rational mode, and assumes that every entry is an integer. (Any denominator encountered will be ignored, with unpredictable results.) By row and column manipulations, the argument is converted to a diagonal matrix of non-negative integers. Furthermore, each integer on the diagonal divides all the following integers. The set of such integers is an invariant of the matrix.

As with Normalize, you can omit any of the change of basis matrices.

If you do not require any conversion matrices then it is possible to greatly speed up Smith in most cases by working modulo a "pseudo-determinant", a multiple of the gcd of the determinants of all the maximal rank minors (see Kannan and Backem, SIAM Journal of Computing vol 8, no. 4, Nov 1979). Do this in Fermat with the command MSmith. Not to work modulo some number invites a horrendous explosion in the intermediate entries. On the other hand, for relatively small matrices or sparse matrices, it's faster to forgo the modding out. Fermat will compute the pseudo-determinant if the matrix is Sparse. If you already have a pseudo-determinant pd to use, use the syntax MSmith([x], pd). If the matrix is not Sparse, you must use the latter method. Pseudet may be helpful.

Hermite = Column reduce a matrix of integers. The matrix may be Sparse. This function can only be used in rational mode, and assumes that every entry is an integer. By column manipulations and row permutations, the argument is converted to a lower triangular matrix of integers. All diagonal entries are non-negative. This is often referred to as Hermite normal form.

If requested, Fermat will also return the integer change of basis (or conversion) matrices used in normalizing, exactly as in Smith.

Be aware that if the matrix is "large" and "dense" a horrendous explosion is possible in the intermediate entries, and in the entries of the conversion matrices.

Hampath = a pretty fast algorithm for finding a Hamiltonian path in a graph. Fermat does not really do Graph Theory; this is here because it's a sparse matrix application. The algorithm is from an exercise in a text book by Papadimitriou. Given the n X n adjacency matrix of a simple graph, the program constructs an n²+1 X n²+1 sparse matrix. Each entry is either 0, 1, or a polynomial variable x_i, i = 1, ..., n. The graph has a Hamiltonian path iff the term x₁ x₂ ... x_n appears in its determinant. We do not actually use any polynomial variables, nor do we compute the determinant, but simulate it to see if such a term would be present. We also throw a bit of randomness into the search. An exhaustive search on a reasonable graph of more than 15 or so nodes would be too time consuming if there is no such path in the graph. But if there is a Hamiltonian path, this method often proves it very quickly. Hampath returns 1 or 0 for path or no path.

The basic invocation is Hampath[w] where [w] is the adjacency matrix, symmetrical with each entry either 0 or 1. However, on any graph of more than 15 or so nodes, it is better to use the alternate form Hampath([w], k), where k is a positive integer. k controls the probabalistic search by terminating a search tree once more than k leaves in the tree are visited with no resolution, and starting over. A reasonable heuristic for k is around 2n². The command &W will display the total number of branches on the search tree and the leaf count.

Disclaimers: Hampath does not first check elementary properties of a graph that might easily decide the issue, like connectivity. That's up to the user. I make no claims about efficiency relative to other methods.

Redrowech = the reduced row echelon form, for elementary matrix equations of the form AX = B. (Other Fermat commands do column manipulations as well, which could be used to solve AX = B but take an extra step.) Invoke with Redrowech([a]), where all columns but the last in [a] represent the matrix A and the last represents B (i.e., Redrowech never pivots on the last column.) Altnately, Redrowech([a],[u],[v]) will return in [u] the transition matrix used in normalizing [a]. [v] is [u]^-1. As in other similar Fermat commands, you can also do Redrowech([a], ,[v]).

Minors: extract minors from sparse arrays. The syntax is e.g. [y] : = Minors([x],[r],[c]). [x] is an existing sparse array. [r] and [c] are existing ordinary arrays specifying the rows and columns to be extracted. The result is stored in [y], which will be a new sparse array of the right dimensions. [x] is untouched.

FFLU and FFLUC are for fraction-free LU factorization of matrices. See the two articles in the September 1997 SIGSAM Bulletin: "Fraction-free Algorithms for Linear and Polynomial Equations," by Nakos, Turner, and Williams; and "The Turing Factorization of a Rectangular Matrix," by Corless and Jeffrey.
    FFLU is invoked as: FFLU([x], [p], [l], [a], [b]). [x] is the n X m matrix to be factored. [p] is an n X n diagonal matrix consisting of the pivots used, [p] = diag(p_1, p_2, ... , p_(n-1), 1). [l] is the unit lower triangular matrix, the first factor. [a] (optional) is the n X n permutation matrix of row swaps. [b] (optional) is [a]^-1. At the end, [x] is in upper triangular form. Let [z] be a copy of the original [x]. If [f] and [g] are the matrices called f_1 and f_2 in the Corless and Jeffrey article, then at the end one has [f]*[a]*[z] = [l]*[g]*[x]. Note that [f] and [g] are not computed by FFLU; however it is obvious how to get them from [p]. Note also that [p] is not necessarily the diagonal of [x]: if columns of 0s are encountered along the way, [x] will be in row-echelon form and may have 0s on its main diagonal.
    FFLUC allows column swaps as well as row swaps. In this way, the size of the pivots can be reduced. FFLUC is invoked as: FFLUC([x], [p], [l], [a], [b], [c], [d]). As above, [x] is the n X m matrix to be factored. [p], [l], [a], and [b] are the same as above. At the end, [x] is in upper triangular form. [c] and [d] (optional) are permutation matrices coming from column swaps ([d] is [c]^-1). Let [z] be a copy of the original [x]. If [f] and [g] are the matrices called f_1 and f_2 in the Corless and Jeffrey article, then at the end one has [f]*[a]*[z]*[c] = [l]*[g]*[x].
     [p], [l], [a], etc. need not be existing matrices when the function is invoked. Matrices of those names with the right size will be created at the end.
    Saying that [a], [c], etc. are optional above means that they may be ommitted, for example FFLUC([x], [p], [l], [a], , [c]). Note the space to indicate no [b].

Sparse Access Loops: There is a need for a way to work efficiently with sparse arrays. For example, suppose you have a sparse array of 60000 rows and 50000 columns with only 10 or so entries in each row (this is quite realistic). Suppose you wanted to add up all the entries. Naively, one could write something like:

for  i  = 1,60000  do  for  j = 1,50000  do  sum  : =  sum + x[i,j]  od  od

But this will do 3,000,000,000 additions, almost all of which are adding 0! This is a preposterous waste of time. The solution is "sparse column access loops" for sparse arrays. The syntax is, continuing the example above,

for  i  = 1,60000  do  for  j = [x]i  do  sum  : =  sum + x[i,j]  od  od.

"for j = [x]i do" means find the i^th row of [x] and let j run down it - of course encountering only the entries actually there! So j takes on whatever the column indices are in which x[i,j] <> 0. [x] must be an existing sparse array, and i must have a value suitable for [x] at the start of the loop. More generally, one may use the syntax: for j = [x]i,k do ... . Here i and k both refer to rows of the sparse matrix [x]. At the start of the loop, all nonzero column coords in both rows are found. Then as the loop proceeds, j runs through those values in order. Any number of row indices is allowed. There is no analogous procedure for "sparse row loops" due to the way Fermat stores sparse matrices. If necessary, transpose the matrix.

Footnotes:

² This changed in MPW version 3.2. If you want to try it, have fun.

³ If you try setting a larger value, you don't get any kind of error message, but it won't work!