REDUCE
UpNextPrevPrevTailTail

9.12 Polynomial Coefficient Arithmetic

REDUCE allows for a variety of numerical domains for the numerical coefficients of polynomials used in calculations. The default mode is integer arithmetic, although the possibility of using real coefficients has been discussed elsewhere. Rational coefficients have also been available by using integer coefficients in both the numerator and denominator of an expression, using the on div option to print the coefficients as rationals. However, REDUCE includes several other coefficient options in its basic version which we shall describe in this section. All such coefficient modes are supported in a table-driven manner so that it is straightforward to extend the range of possibilities. A description of how to do this is given in [BHPS86].

9.12.1 Rational Coefficients in Polynomials

Instead of treating rational numbers as the numerator and denominator of a rational expression, it is also possible to use them as polynomial coefficients directly. This is accomplished by turning on the switch rational.

Example: With rational off, the input expression a/2 would be converted into a rational expression, whose numerator was a and denominator 2. With rational on, the same input would become a rational expression with numerator 1/2*a and denominator 1. Thus the latter can be used in operations that require polynomial input whereas the former could not.

9.12.2 Real Coefficients in Polynomials

The switch rounded permits the use of arbitrary sized real coefficients in polynomial expressions. The actual precision of these coefficients can be set by the operator precision. For example, precision 50; sets the precision to fifty decimal digits. The default precision is system dependent and can be found by precision 0;. In this mode, denominators are automatically made monic, and an appropriate adjustment is made to the numerator.

Example: With ROUNDED on, the input expression a/2 would be converted into a rational expression whose numerator is 0.5*a and denominator 1.

Internally, REDUCE uses floating point numbers up to the precision supported by the underlying machine hardware, and so-called bigfloats for higher precision or whenever necessary to represent numbers whose value cannot be represented in floating point. The internal precision is two decimal digits greater than the external precision to guard against roundoff inaccuracies. Bigfloats represent the fraction and exponent parts of a floating-point number by means of (arbitrary precision) integers, which is a more precise representation in many cases than the machine floating point arithmetic, but not as efficient. If a case arises where use of the machine arithmetic leads to problems, a user can force REDUCE to use the bigfloat representation at all precisions by turning on the switch roundbf. In rare cases, this switch is turned on by the system, and the user informed by the message 

        ROUNDBF turned on to increase accuracy

Rounded numbers are normally printed to the specified precision. However, if the user wishes to print such numbers with less precision, the printing precision can be set by the command print_precision. For example, print_precision 5; will cause such numbers to be printed with five digits maximum.

Under normal circumstances when rounded is on, REDUCE converts the number 1.0 to the integer 1. If this is not desired, the switch noconvert can be turned on.

Numbers that are stored internally as bigfloats are normally printed with a space between every five digits to improve readability. If this feature is not required, it can be suppressed by turning off the switch bfspace.

Further information on the bigfloat arithmetic may be found in T. Sasaki, “Manual for Arbitrary Precision Real Arithmetic System in REDUCE”, Department of Computer Science, University of Utah, Technical Note No. TR-8 (1979).

When a real number is input, it is normally truncated to the precision in effect at the time the number is read. If it is desired to keep the full precision of all numbers input, the switch adjprec (for adjust precision) can be turned on. While on, adjprec will automatically increase the precision, when necessary, to match that of any integer or real input, and a message printed to inform the user of the precision increase.

When rounded is on, rational numbers are normally converted to rounded representation. However, if a user wishes to keep such numbers in a rational form until used in an operation that returns a real number, the switch roundall can be turned off. This switch is normally on.

Results from rounded calculations are returned in rounded form with two exceptions: if the result is recognized as 0 or 1 to the current precision, the integer result is returned.

9.12.3 Modular Number Coefficients in Polynomials

REDUCE includes facilities for manipulating polynomials whose coefficients are computed modulo a given base. To use this option, two commands must be used; setmod \(\langle \)integer\(\rangle \), to set the prime modulus, and on modular to cause the actual modular calculations to occur. For example, with setmod 3; and on modular;, the polynomial (a+2*b)^3 would become a^3+2*b^3.

The argument of setmod is evaluated algebraically, except that non-modular (integer) arithmetic is used. Thus the sequence 

        setmod 3; on modular; setmod 7;

will correctly set the modulus to 7.

Modular numbers are by default represented by integers in the interval [0,p-1] where p is the current modulus. Sometimes it is more convenient to use an equivalent symmetric representation in the interval [-p/2+1,p/2], or more precisely [-floor((p-1)/2), ceiling((p-1)/2)], especially if the modular numbers map objects that include negative quantities. The switch balanced_mod allows you to select the symmetric representation for output.

Users should note that the modular calculations are on the polynomial coefficients only. It is not currently possible to reduce the exponents since no check for a prime modulus is made (which would allow \(x^{p-1}\) to be reduced to 1 mod p). Note also that any division by a number not co-prime with the modulus will result in the error “Invalid modular division”.

9.12.4 Complex Number Coefficients in Polynomials

Although REDUCE routinely treats the square of the variable i as equivalent to \(-1\), this is not sufficient to reduce expressions involving i to lowest terms, or to factor such expressions over the complex numbers. For example, in the default case, 

        factorize(a^2+1);

gives the result 

        {{a**2+1,1}}

and 

        (a^2+b^2)/(a+i*b)

is not reduced further. However, if the switch complex is turned on, full complex arithmetic is then carried out. In other words, the above factorization will give the result 

        {{a + i,1},{a - i,1}}

and the quotient will be reduced to a-i*b.

The switch complex may be combined with rounded to give complex real numbers; the appropriate arithmetic is performed in this case. Similarly, combining complex with rational performs polynomial arithmetic with complex rational numbers.

Complex conjugation is used to remove complex numbers from denominators of expressions. To do this if complex is off, you must turn the switch rationalize on.

9.12.5 Algebraic Numbers as Coefficients in Polynomial

The ARNUM package1 provides facilities for handling algebraic numbers as polynomial coefficients in REDUCE calculations. It includes facilities for introducing indeterminates to represent algebraic numbers, for calculating splitting fields, and for factoring and finding greatest common divisors in such domains.

Algebraic numbers are the solutions of an irreducible polynomial over some ground domain. The algebraic number \(i\) (imaginary unit), for example, would be defined by the polynomial \(i^2 + 1\). The arithmetic of algebraic number \(s\) can be viewed as a polynomial arithmetic modulo the defining polynomial.

Given a defining polynomial for an algebraic number \(a\)

\[ a^n ~ + ~ {p _{n-1}} {a ^ {n -1}} ~ + ~ ... ~ + ~ {p_0} \]

All algebraic numbers which can be built up from \(a\) are then of the form:

\[ {r_{n-1}} {a ^{n-1}} ~+~ {r_{n-2}} {a ^{n-2}} ~+~ ... ~+~ {r_0} \]
where the \(r_j\)’s are rational numbers.

The operation of addition is defined by

\begin{multline} ({r_{n-1}} {a ^{n-1}} ~+~ {r_{n-2}} {a ^{n-2}} ~+~ ...) ~ + ~ ({s_{n-1}} {a ^{n-1}} ~+~ {s_{n-2}} {a ^{n-2}} ~+~ ...) ~ = \\ ({r_{n-1}+s_{n-1}}) {a ^{n-1}} ~+~ ({r_{n-2}+s_{n-2}}) {a ^{n-2}} ~+~ ... \end{multline}

Multiplication of two algebraic numbers can be performed by normal polynomial multiplication followed by a reduction of the result with the help of the defining polynomial.

\begin{multline} ({r_{n-1}} {a ^{n-1}} + {r_{n-2}} {a ^{n-2}} + ...) ~ \times ~ ({s_{n-1}} {a ^{n-1}} + {s_{n-2}} {a ^{n-2}} + ...) = \\ {r_{n-1}} {s ^{n-1}}{a^{2n-2}} + ... ~ \mathbin {\mathbf {modulo}} ~ (a^n ~ + ~ {p _{n-1}} {a ^ {n -1}} ~ + ~ ... ~ + ~ {p_0}) \\ = ~~~{q_{n-1}} a^{n-1} ~ + ~ {q _{n-2}} {a ^ {n -2}} ~ + ~ ... \end{multline}

Division of two algebraic numbers r and s yields another algebraic number q.

\( \frac {r}{s} = q\) or \( r = q s \).

The last equation written out explicitly reads

\begin{align*} ({r_{n-1}} & {a^{n-1}} + {r_{n-2}} {a^{n-2}} + \ldots ) \\ =\ & ({q_{n-1}} {a^{n-1}} + {q_{n-2}} {a^{n-2}} + \ldots ) \times ({s_{n-1}} {a^{n-1}} + {s_{n-2}} {a^{n-2}} + \ldots ) \\ & \mathbin {\mathbf {modulo}} (a^n + {p _{n-1}} {a^{n -1}} + \ldots ) \\ =\ & ({t_{n-1}} {a^{n-1}} + {t_{n-2}} {a^{n-2}} + \ldots ) \end{align*}

The \(t_i\) are linear in the \(q_j\). Equating equal powers of \(a\) yields a linear system for the quotient coefficients \(q_j\).

With this, all field operations for the algebraic numbers are available. The translation into algorithms is straightforward. For an implementation we have to decide on a data structure for an algebraic number. We have chosen the representation REDUCE normally uses for polynomials, the so-called standard form. Since our polynomials have in general rational coefficients, we must allow for a rational number domain inside the algebraic number.

\(<\) algebraic number \(>\) ::=
  :ar: . \(<\) univariate polynomial over the rationals \(>\)
\(<\) univariate polynomial over the rationals \(>\) ::=
  \(<\) variable \(>\) .** \(<\) ldeg \(>\) .* \(<\) rational \(>\) .+ \(<\) reductum \(>\)
\(<\) ldeg \(>\) ::= integer

\(<\) rational \(>\) ::=
  :rn: . \(<\) integer numerator \(>\) . \(<\) integer denominator \(>\) : integer
\(<\) reductum \(>\) ::= \(<\) univariate polynomial \(>\) : \(<\) rational \(>\) : nil

This representation allows us to use the REDUCE functions for adding and multiplying polynomials on the tail of the tagged algebraic number. Also, the routines for solving linear equations can easily be used for the calculation of quotients. We are still left with the problem of introducing a particular algebraic number. In the current version this is done by giving the defining polynomial to the statement defpoly. The algebraic number sqrt(2), for example, can be introduced by 

   defpoly sqrt2**2 - 2;

This statement associates a simplification function for the translation of the variable in the defining polynomial into its tagged internal form and also generates a power reduction rule used by the operations times and quotient for the reduction of their result modulo the defining polynomial. A basis for the representation of an algebraic number is also set up by the statement. In the working version, the basis is a list of powers of the indeterminate of the defining polynomial up to one less then its degree. Experiments with integral bases, however, have been very encouraging, and these bases might be available in a later version. If the defining polynomial is not monic, it will be made so by an appropriate substitution.

Example 1 

defpoly sqrt2**2-2;

1/(sqrt2+1);

sqrt2 - 1

(x**2+2*sqrt2*x+2)/(x+sqrt2);

x + sqrt2

on gcd;

(x**3+(sqrt2-2)*x**2-(2*sqrt2+3)*x-3*sqrt2)/(x**2-2);

  2
(x  - 2*x - 3)/(x - sqrt2)

off gcd;

sqrt(x**2-2*sqrt2*x*y+2*y**2);

abs(x - sqrt2*y)

Until now we have dealt with only a single algebraic number. In practice this is not sufficient as very often several algebraic numbers appear in an expression. There are two possibilities for handling this: one can use multivariate extensions [Dav81] or one can construct a defining polynomial that contains all specified extensions. This package implements the latter case (the so called primitive representation). The algorithm we use for the construction of the primitive element is the same as given by Trager [Tra76]. In the implementation, multiple extensions can be given as a list of equations to the statement defpoly, which, among other things, adds the new extension to the previously defined one. All algebraic numbers are then expressed in terms of the primitive element.

Example 2 

defpoly sqrt2**2-2,cbrt5**3-5;

*** defining polynomial for primitive element:

  6       4        3        2
a1  - 6*a1  - 10*a1  + 12*a1  - 60*a1 + 17

sqrt2;

  48     5     45     4    320     3    780     2
------*a1  + ------*a1  - ------*a1  - ------*a1
 1187         1187         1187         1187

    735         1820
 + ------*a1 - ------
    1187        1187

sqrt2**2;

   2

We can provide factorization of polynomials over the algebraic number domain by using Trager’s algorithm. The polynomial to be factored is first mapped to a polynomial over the integers by computing the norm of the polynomial, which is the resultant with respect to the primitive element of the polynomial and the defining polynomial. After factoring over the integers, the factors over the algebraic number field are recovered by GCD calculations.

Example 3 

   defpoly a**2-5;

   on factor;

   x**2 + x - 1;

   (x + (1/2*a + 1/2))*(x - (1/2*a - 1/2))

We have also incorporated a function split_field for the calculation of a primitive element of minimal degree for which a given polynomial splits into linear factors. The algorithm as described in Trager’s article is essentially a repeated primitive element calculation.

Example 4 

   split_field(x**3-3*x+7);

   *** Splitting field is generated by:

     6        4        2
   a2  - 18*a2  + 81*a2  + 1215



            4          2
   {1/126*a2  - 5/42*a2  - 1/2*a2 + 2/7,


              4          2
    - (1/63*a2  - 5/21*a2  + 4/7),


           4          2
   1/126*a2  - 5/42*a2  + 1/2*a2 + 2/7}


   for each j in ws product (x-j);

    3
   x  - 3*x + 7

A more complete description can be found in [BHPS86].

UpNextPrevPrevTailFront

Hosted by Download REDUCE Powered by MathJax