Contents Previous Next Index
7 Numerical Programming in Maple
An important part of efficient scientific and mathematical programming is numerical computation. Maple provides many tools for computing with floating-point numbers, some for improving efficiency and some for improving accuracy.
7.1 In This Chapter
An Overview of Numeric Types in Maple
An Explanation of Floating-Point Numbers in Maple
Maple Commands for Numerical Computing
Efficient Numerical Programs
7.2 Numeric Types in Maple
Before discussing numerical computing in Maple, we will first introduce the various numeric data types used in Maple and briefly describe how they are represented. All of the real numbers discussed in this section will pass checks of type,numeric or type,extended_numeric.
Integers
The most basic numeric type in Maple is the integer. Small integers are represented directly as hardware integers (similar to the int type in C or integer type in Fortran), which allows for maximum efficiency of both CPU time used for arithmetic and memory used for storage. That is, the number can be stored in one machine word and arithmetic operations can be performed with one CPU operation. On 32-bit architectures, integers in the range −230−1 to 230−1 are stored in this way, while on 64-bit architectures, integers in the range −262−1 to 262−1. Integers stored in this way are referred to as immediate integers.
Larger integers are stored in DAGs of type INTPOS or INTNEG, which contain pointers to arrays of digits that can store integers up to magnitude 10228−218 on 32-bit architectures and 10235+232−18 on 64-bit architectures.
dismantle(2^80-1);
INTPOS(6): 1208925819614629174706175
dismantle(-2^101+6);
INTNEG(6): -2535301200456458802993406410746
The arithmetic for these large integers is computed using the GNU Multiple Precision Arithmetic (GMP) library. This library is quite efficient, but still several orders of magnitude slower than arithmetic on immediate integers since each arithmetic operation will require more than one CPU operation and the larger the integer, the more operations and memory will be needed for arithmetic.
CodeTools:-Usage(add(i,i=-2^15..2^16));
memory used=161.02KiB, alloc change=0 bytes, cpu time=4.00ms, real time=3.00ms, gc time=0ns
1610629120
CodeTools:-Usage(add(i,i=2^88-2^15..2^88+2^16));
memory used=8.92MiB, alloc change=0 bytes, cpu time=16.00ms, real time=16.00ms, gc time=0ns
30423923890487326980991212339200
CodeTools:-Usage(add(i,i=2^4097-2^15..2^4097+2^16));
memory used=109.34MiB, alloc change=-3.01MiB, cpu time=562.00ms, real time=325.00ms, gc time=499.32ms
205337297974639914340665500453995519859046771005206125061344145148458177846287911591047839494539869945402403881125903466576736813681708881124024847232456927296954601158258822773698492236321182857799702845030997587567848047013984302926637722288989174866047402930648066257455156822794719737995549112033208483004719072244728296595021115951707203531283316339427422299668078787669834908390333537991291917652425009429047159125028756356966402889004692279651401969448167419628544878495940166131118356838592642274676376886478814390080150121184593927849550735415253857623428878524437164429933525844071108151834935561910383653050597127603219544447136154079024731191638245863569443556792160565672892410184321563392936412391807603283785170896557622787216442839497627198271642314069273549773891778249817599015826649550121400193546930160852466240321147395833226174991830052909513217459249145986409972890607772408436946645869934972341681160127289004812013361613475609253983756470541201044624106898890791412083639760908985411624428440457165177561278545720011754678944275988097254471054453900542709445889615221249910582150846966313990266569518867263462971416395393555619749389086067511784172161221094396003584450254897514753869702834443706806664146972590080
Any transitions between GMP integers and immediate integers will be completely transparent and it is not possible to tell them apart in general without use low-level tools such as addressof. However, you can check if an integer is small enough to fit into a single machine word with types integer[4] and integer[8] for 4-byte and 8-byte words respectively.
Integers of all types pass a type,integer type check.
The Integer constructor is guaranteed to return an integer, an extended numeric symbol such as infinity or undefined, a complex number with integer parts, or return unevaluated.
Integer(-2^160);
−1461501637330902918203684832716283019655932542976
Integer(infinity);
∞
Integer(1/2);
Integer⁡12
The system dependent value for the largest immediate integer can be found with kernelopts(maximmediate), the maximum number of decimal digits in an integer can be found with kernelopts(maxdigits), and the version of the GMP library being used can be found with kernelopts(gmpversion).
Rationals
Exact rational numbers are stored in DAGs of type RATIONAL, which consist of a pair of integers. The first integer is the numerator and can be a POSINT or NEGINT. The second integer is the denominator and is a POSINT. Most low-level programming languages such as C or Fortran do not have an equivalent rational number type.
dismantle(1/2);
RATIONAL(3): 1/2 INTPOS(2): 1 INTPOS(2): 2
dismantle(-10/3);
RATIONAL(3): -10/3 INTNEG(2): -10 INTPOS(2): 3
Rational numbers can be constructed by using the division operator or the Fraction constructor. In either case, automatic simplification will occur to ensure that the denominator is positive and that the fraction is in lowest terms (the numerator and denominator do not have factors in common). This means that the Fraction constructor may return integers in some cases.
dismantle(Fraction(21,7));
INTPOS(2): 3
dismantle(Fraction(40,-14));
RATIONAL(3): -20/7 INTNEG(2): -20 INTPOS(2): 7
Rational number arithmetic is performed in the natural way using integer arithmetic and the igcd and ilcm operations to reduce to lowest terms.
Fraction(2^20+2^12,2^27-2^13) + Fraction(2^12-1,2^13);
68141057134209536
Fraction(2^20+2^12,2^27-2^13) * Fraction(23,187);
59116127242
Rational numbers of all types will pass a type,rational type check. Only rational numbers that are not also integers will pass a type,fraction type check. Additionally, type,extended_rational includes all rationals as well as the extended numeric symbols infinity, -infinity, and undefined.
type(1, fraction);
false
Like the Integer constructor, the Fraction constructor will return unevaluated if it cannot return a value of type extended_rational.
Fraction(x,1);
Fraction⁡x,1
Fraction(infinity);
Floating-Point Numbers
Floating-point numbers are stored in DAGs of type FLOAT.
In Maple, as in nearly every programming language, floating-point numbers can be constructed using and visually distinguished from integers with a decimal point symbol, '.'. The floating-point number 1. is often treated as equal to the exact integer 1.
evalb(1. = 1);
true
Maple floating-point numbers can also be constructed with the SFloat constructor (or the equivalent Float constructor) and can be checked with the nearly equivalent type,sfloat and type,float types. We will generally refer to these numbers as sfloats to when we need to distinguish them from hardware floating-point numbers (hfloats), introduced below.
Float(1);
1.
dismantle(SFloat(0.3333));
FLOAT(3): .3333 INTPOS(2): 3333 INTNEG(2): -4
type(.1, float);
type(.1, sfloat);
type(1, float);
A floating-point number represents a rational number with a fixed precision. That rational number can be recovered with convert/rational.
convert(.3333333333, rational, exact);
333333333310000000000
However, not every rational number can be represented exactly by a floating-point number. For example, the closest floating-point number to 13 is 0.3333333333.
convert(1/3, float);
0.3333333333
Also, unlike numeric types integer and rational, integer and float do not have compatible arithmetic. Floating-point arithmetic has a fixed finite precision, and does round off if the result of arithmetic does not fit into that precision.
9123456789 + 8123456789 <> convert( 9123456789. + 8123456789., rational, exact);
17246913578≠17246913580
123456 * 1234567 <> convert( 123456.*1234567., rational, exact);
152414703552≠152414703600
Unlike many other programming languages the precision of sfloat arithmetic can be changed. For this reason, sfloats are known as arbitrary precision floating-point numbers.
More information on sfloats and how they differ from the floating-point types in languages such as C and Fortran will be discussed in greater detail in More about Floating-Point Numbers in Maple.
Hardware Floating-Point Numbers
Floating-point numbers of the type used in languages such as C and Fortran can also be constructed in Maple; they are known as hardware floating-point numbers or hfloats. These types are stored as 8-byte double precision IEEE floating-point numbers contained in DAGs of type HFLOAT. Since the . notation is used to construct Maple sfloats, hfloats must be constructed with the HFloat constructor. Maple will display sfloats and hfloats the same way, using decimal notation.
HFloat(1);
dismantle(HFloat(0.3333));
HFLOAT(2): .3333
The advantage of hfloats over sfloats is that their arithmetic is computed directly using a single CPU operation for each arithmetic operation. Maple sfloats, however, offer much more flexibility and precision. In many ways the difference is analogous to the difference between immediate integers and GMP integers.
Hardware floats can be distinguished from sfloats with the type,hfloat type.
type(HFloat(1), float);
type(HFloat(1), sfloat);
type(HFloat(1), hfloat);
type(SFloat(1), hfloat);
For more information on hardware floats and how they differ from sfloats, see More about Floating-Point Numbers in Maple.
Extended Numeric Types
The special built-in symbols infinity (∞), and undefined can be used in numeric arithmetic in Maple. In general, operations involving ∞ simplify automatically to a signed infinity or a complex infinity. For details, refer to the type,infinity help page.
-1*infinity;
−∞
1/2*infinity;
1/infinity;
0
The undefined symbol is usually produced as the result of attempting to carry out an operation that cannot result in a number for the given operands. Almost every arithmetic operation involving undefined returns undefined. For details, refer to the type,undefined help page.
infinity-infinity;
undefined
undefined-undefined;
undefined+1;
Integer and rational numbers share exact undefined and infinite symbols while sfloat and hfloat numbers have their own versions of these, which are displayed differently but treated similarly.
Float(infinity);
Float⁡∞
HFloat(undefined);
Float⁡undefined
Complex Numbers
A complex number in Maple is a DAG of type COMPLEX, which consists of a pair of any of the two numeric types. They can be constructed in the natural way using the symbol I for the imaginary unit, or using the Complex constructor.
dismantle(1+I);
COMPLEX(3) INTPOS(2): 1 INTPOS(2): 1
dismantle(Complex(1/2,1/3));
COMPLEX(3) RATIONAL(3): 1/2 INTPOS(2): 1 INTPOS(2): 2 RATIONAL(3): 1/3 INTPOS(2): 1 INTPOS(2): 3
Automatic simplification will ensure that if one of the parts of a complex number is a float (or hfloat), then other will be made into a float (hfloat).
dismantle(Complex(1., 1/1001));
COMPLEX(3) FLOAT(3): 1. INTPOS(2): 1 INTPOS(2): 0 FLOAT(3): .9990009990e-3 INTPOS(2): 9990009990 INTNEG(2): -13
dismantle(Complex(HFloat(1.), 1/1001));
COMPLEX(3) HFLOAT(2): 1. HFLOAT(2): .000999000999
dismantle(Complex(HFloat(1.), 1.));
COMPLEX(3) HFLOAT(2): 1. HFLOAT(2): 1.
Complex numbers are not of type type,numeric but can be checked with type type,complex which can also be structured to check for the numeric subtypes of its two components.
type(1+I,numeric);
type(1+I,complex(integer));
Non-numeric Constants
Many Maple expressions represent constants, but are not considered to be of type numeric. This means that arithmetic performed on these constants will be more generic symbolic operations on DAGs of type SUM, PROD, NAME, or FUNCTION. Some examples of non-numeric constants are Pi (π), Catalan, sin⁡1, 5, and π+ⅇ2−1+5⁢Catalan.
type(Pi, numeric);
type(sqrt(5)-1, constant);
7.3 More about Floating-Point Numbers in Maple
To take full advantage of floating-point numbers and to avoid many common pitfalls in numerical computing, it is important to understand exactly what floating-point numbers are and how they are represented.
Representation of Floating-Point Numbers in Maple
The dismantle command shows that the two numbers 1 and 1. have different internal representations. 1 is simply stored as an integer while 1. is stored as a pair of integers.
dismantle(1);
INTPOS(2): 1
dismantle(1.);
FLOAT(3): 1. INTPOS(2): 1 INTPOS(2): 0
Similarly, the numbers 12 and 0.5 are also different even though they are both stored as pairs of integers.
dismantle(0.5);
FLOAT(3): .5 INTPOS(2): 5 INTNEG(2): -1
In Maple, the FLOAT DAG-type represents a floating-point number in the form S * 10^E where both S and E are integers. For 1., the significand (or mantissa) is S=1 and the exponent is E=0. In addition to being specified in decimal notation, floats of this form can be constructed by using scientific notation, or the Float constructor.
Float(2,0);
2.
2*1e0;
The advantage of using this significand-exponent representation is that fixed precision approximations of large and small numbers can be stored compactly and their arithmetic can be done efficiently. Storing the integer 10^50 needs at least 167 bits or 3 words on a 64-bit machine. The floating-point number 1e50 can be stored in less than 8 bits but in in practice uses 2 words (one for each integer).
dismantle(10^50);
INTPOS(8): 100000000000000000000000000000000000000000000000000
dismantle(1e50);
FLOAT(3): .1e51 INTPOS(2): 1 INTPOS(2): 50
Using two immediate integers, a float can store a much larger range of numbers than a rational number with two immediate integers. The range a rational can represent is about 1.⁢10-9..1.⁢109 while a float can represent a range of about 1.⁢10-1073741823..9.⁢101073741823. This is a much larger range for the same storage cost. Of course, that larger range means that floats of a fixed size can represent fewer numbers in that range. And since floating-point numbers are always of a fixed size, this means that arithmetic will always be of limited precision. That is, each operation will have to round the result to a number that can be represented as another floating-point number.
In Maple, the significand is limited to 10 decimal digits of precision by default but can be changed while the exponent is restricted to being a word-sized integer.
More information on the restrictions on the size of software floats in Maple can be found by using the Maple_floats command.
By contrast, hfloats, are represented in base 2, rather than base 10. So they represent numbers using the form S * 2^E, where the significand, S, is a 52-bit integer and the exponent, E, is a 10-bit integer. Thus, the range of numbers representable as hardware floats is 2.225073859⁢10-308..1.797693135⁢10308. Because the largest possible significand of a hardware float has about floor⁡log10⁡252=15 base-10 digits of precision, hardware floats can be converted to software floats without meaningful loss of precision when Digits is 15. Conversely, so long as their exponent is smaller than 307 and their significand had fewer than 15 digits sfloats can be converted to hfloats without loss of precision.
Precision and Accuracy
By default, 10-digit precision is used for floating-point arithmetic, which means that the arithmetic will be rounded to 10 digits. This means any single floating-point operation will be accurate to 10 digits.
For example, storing 10^50+1 requires 50 decimal digits so it will be rounded in floating-point arithmetic. By contrast, 10^50+10^41 can be stored with 10 digits so it will still be computed accurately.
.1e51 + 1.;
1.×1050
.1e51 + .1e42;
1.000000001×1050
The Digits environment variable can be used to change the working precision used by Maple. Larger values of Digits will allow more accurate computation, but at the cost of slower arithmetic.
Digits := 100:
1.00000000000000000000000000000000000000000000000001×1050
The maximum value for Digits is system dependent and can be found with the Maple_floats command.
Maple_floats(MAX_DIGITS);
38654705646
For the default value of Digits, the significand is an immediate integer and so arithmetic will be fast in general. It also means that some numerical function evaluations (such as sin in the following example) will be able to use the CPU's native hardware floating-point arithmetic to achieve the needed precision. However, raising Digits about the default value will lead to slower arithmetic and slower function evaluation.
Digits:=10:
CodeTools:-CPUTime(add(sin(1e-6*i),i=1..100000));
1.300,4995.884639
Digits:=22:
6.494,4995.884638682140998954
Reducing Digits below its default value does not usually lead to large improvements in performance.
Digits:=5:
1.202,4996.0
It is also important to note that changing Digits does not necessarily change the accuracy of sequences of multiple floating-point computations; it changes only the precision of the individual operations performed. The following example computes two additions using 10 digits of precision, but catastrophic cancellation leads to a mere one digit of accuracy in the final answer.
Digits := 10:
x := 1234567890.;
x≔1.234567890×109
y := -x+1;
y≔−1.234567889×109
z := 3.141592654;
z≔3.141592654
x+z+y<>z+1;
4.≠4.141592654
Ensuring accuracy requires careful study of the problem at hand. In this example, you need 19 digits of precision to get 10 digits of accuracy.
Digits := 19:
x+z+y=z+1;
4.141592654=4.141592654
Floating-Point Contagion
An important property of floating-point numbers in Maple, and nearly every other computing environment, is contagion. When numerical expressions are created involving both floating-point numbers and exact numbers, the floating property is contagious and causes the answer to become a floating-point number.
1. * 10;
10.
0. + 10;
As you can see, this contagion property can be used as a quick method to convert exact values to floating-point numbers. However, while floating-point contagion extends to all Maple structures of type numeric (except, in some cases, hfloats), it does not apply to non-numeric constants.
type(3/4, numeric);
4/3 + 0.;
1.333333333
1.*sqrt(3);
1.⁢3
The hfloat type is also contagious, but the precise behavior of the contagion is determined by the UseHardwareFloats environment variable. By default, hfloats are contagious for small values of Digits:
type(4/3 + HFloat(0.), hfloat);
type(1. + HFloat(0.), hfloat);
HFloat(1.1) * sin(4*Pi/7) -1;
1.10000000000000⁢sin⁡3⁢π7−1
For large values of Digits, hfloats in computations will be converted to sfloats so that the results are sfloats.
Digits := 20;
Digits≔20
type(1 + HFloat(0.), hfloat);
type(1 + HFloat(0.), sfloat);
If UseHardwareFloats=true then hfloats are completely contagious.
UseHardwareFloats := true;
UseHardwareFloats≔true
a := 10.^19+1;
a≔1.0000000000000000001×1019
b := a + HFloat(0.1);
b≔1.00000000000000×1019
type(b, hfloat);
If UseHardwareFloats=false then hfloats will always be converted to sfloats in computations, regardless of the setting of Digits. The HFloat constructor will still create hfloats, however.
UseHardwareFloats := false;
UseHardwareFloats≔false
Digits := 10;
Digits≔10
c := 1 + HFloat(0.1);
c≔1.100000000
type(c, hfloat);
type(HFloat(0.1), hfloat);
Table 7.1 summarizes the floating-point contagion rules.
UseHardwareFloats
deduced
Digits
any
1...15
16...
hfloat + hfloat
hfloat
sfloat
hfloat + sfloat
sfloat + sfloat
More on the Floating-Point Model
The software floating-point system is designed as a natural extension of the industry standard for hardware floating-point computation, known as IEEE 754. Thus, there are representations for infinity and undefined (what IEEE 754 calls a NaN, meaning Not a Number) as discussed in Extended Numeric Types.
The IEEE 754 standard defines five rounding algorithms. Two methods called nearest and simple round to nearest values, and the other three are directed roundings that round up or down (as needed) towards −∞, ∞, or 0. Maple implements all of these rounding modes and the desired mode can be selected by setting the Rounding environment variable.
Rounding;
nearest
1.4^10;
28.92546550
Rounding := 0;
Rounding≔0
28.92546549
Another important feature of this system is that the floating-point representation of zero, 0., retains its arithmetic sign in computations. That is, Maple distinguishes between +0. and -0. when necessary. In most situations, this difference is irrelevant, but when dealing with functions that have a discontinuity across the negative real axis, such as ln⁡x, preserving the sign of the imaginary part of a number on the negative real axis is important.
For more intricate applications, Maple implements extensions of the IEEE 754 notion of a numeric event, and provides facilities for monitoring events and their associated status flags. For more information about this system, refer to the numerics help page.
7.4 Maple Commands for Numerical Computing
In this section we will discuss some of the commands available in Maple for floating-point computation.
The evalf Command
The evalf command is the primary tool in Maple for performing floating-point calculations in software floating-point mode. You can use evalf to compute approximations of non-numeric constants.
evalf(Pi);
3.141592654
You can alter the number of digits of the approximation by changing the value of the environment variable Digits, or by specifying the number as an index to evalf (which leaves the value of Digits unchanged).
Digits := 20:
3.1415926535897932385
evalf[200](Pi);
3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303820
evalf(sqrt(2));
1.4142135623730950488
Remember that the Digits command specifies the precision in decimal digits, unlike hardware floating-point numbers which specify precision in binary digits.
All floating-point computations are performed in finite precision, with intermediate results generally being rounded to Digits precision. As such, it is possible for round-off errors to accumulate in long computations. Maple ensures that the result of any single floating-point arithmetic operation (+, -, *, /, or sqrt) is fully accurate. Further, many of the basic functions in Maple, such as the trigonometric functions and their inverses, the exponential and logarithmic functions, and some of the other standard special functions for mathematics, are accurate to within .6 units in the last place (ulps), meaning that if the Digits + 1st digit of the true result is a 4, Maple may round it up, or if it is a 6, Maple may round it down. Most mathematical functions in Maple, including numerical integration, achieve this accuracy on nearly all inputs.
It is possible to create software floats with different precisions. Changing the value of Digits will not change these numbers; it affects only the precision of subsequent operations on those numbers.
Digits := 50;
Digits≔50
a := evalf(Pi);
a≔3.1415926535897932384626433832795028841971693993751
a;
3.1415926535897932384626433832795028841971693993751
a+1;
4.141592654
evalf(a);
From this example, you can see that evalf can be used to create a lower precision float from one of higher precision. This can be used to round a result to a desired number of digits. However, evalf will not increase the precision of a low precision float.
evalf[100](1.0);
1.0
evalf[10000](a);
The evalf command also provides an interface to purely numerical computations of integrals, limits, and sums.
Some definite integrals have no closed-form solution in terms of standard mathematical functions. You can use evalf to obtain a numerical answer directly using numerical techniques.
r := Int(exp(x^3), x=0..1);
Int⁡exp⁡x3,x=0..1
value(r);
int⁡exp⁡x3,x=0..1
evalf(r);
1.341904418
In other cases, Maple can find an exact solution, but the form of the exact solution is almost incomprehensible. The function Beta in the following example is a special function that appears in mathematical literature.
q := Int( x^99 * (1-x)^199 / Beta(100, 200), x=0..1/5 );
Int⁡277216764217237649652225568421760372018697733716243247244202243820317088554933908000⁢x99⁢1−x199,x=0..15
value(q);
27852290545780521179255248650434305998403849800909690342170417622052715523897761906828166964420518416902474524718187972029459617663867797175746341349064425727501861101435750157352018112989492972548449785454954447636248495323512797804102876034481999911930417847858749936840755474537033615661445973112364349371450421100562106866977667955024449202371857434152360496874313577908566230689757503569126129150390625
evalf(q);
3.546007367×10−8
The two previous examples use the Int command rather than int for the integration. If you use int, Maple first tries to integrate the expression symbolically. Thus, when evaluating the following commands, Maple determines a symbolic answer and then converts it to a floating-point approximation, rather than performing direct numerical integration. In general, the symbolic computation is more difficult, and thus slower than the numerical computation.
evalf( int(x^99 * (1-x)^199 / Beta(100, 200), x=0..1/5) );
Similarly, evalf can be used on the inert forms Limit and Sum to compute using numerical algorithms for computing numeric limits and sums.
evalf(Limit(sin(erf(1)*x)/(erf(1)^2*x),x=0));
1.186660803
evalf( Sum(exp(x), x=RootOf(_Z^5+_Z+1)) );
4.791792042+0.⁢I
When Not to Use evalf
In general the symbolic commands in Maple are able to handle floating-point numbers in their input, but, by their nature floats are not as precise as rationals or symbolic constants. So, even if you want a numerical answer from a command, you should avoid calling evalf on the input.
The following command does not compute the expected answer of 0.1111111111.
limit(n*(evalf(1/3) - 1/(3+1/n)), n=infinity);
Float⁡-∞
It would have been computed correctly with non-float values in the input.
evalf( limit(n*(1/3 - 1/(3+1/n)), n=infinity) );
0.1111111111
Numeric Solvers
There are also a number of numerical algorithms available in Maple in commands other than evalf. One of the most important is fsolve which is short for floating-point solve. This command computes numerical solutions to equations or systems of equations. In general, it is much more efficient than calling evalf on the result of solve, especially if you are interested in only a single solution.
fsolve( exp(x) + 2*sin(x), x);