6.4 The Standard Library (2)--Mathematical Functions

Every implementation of Modula-2, and of every computer language that is to be used in a scientific or academic environment, must also provide a number of standard mathematical functions. A name and contents of the module that contains these were suggested by Wirth--he called it MathLib0. However, in some versions, vendors called it either MathLib or MathLib1, and in a few the procedure names started with an uppercase letter. The ISO standard version is called RealMath. It makes available the two real constants:

CONST
  pi   = 3.1415926535897932384626433832795028841972;
  exp1 = 2.7182818284590452353602874713526624977572;

to as many decimal places as the implementation allows. In addition, the following sections detail the basic functions that are exported by RealMath, with their parameter lists, some comments and examples. (A few have been used before .)

6.4.1 Square Root

PROCEDURE sqrt (x : REAL) : REAL;

This function procedure returns the square root of a positive real number. Naturally, there will be an error generated if one attempts to take the square root of a negative number.

Problem:

Given two sides of a right triangle, compute the hypotenuse.

Refinement:

Ask the user for two numbers, the sides' lengths
  Read each into a real variable
Compute the hypotenuse using Pythagoras' Theorem
  hyp = sqrt (a * a + b * b)
Print out the result

Data Table:

Variables:

Two reals to hold the side lengths, one for the hypotenuse

Imports

WriteString, WriteLn, ReadReal, WriteReal, sqrt

Code:

MODULE Pythagoras;

(* Written by R.J. Sutcliffe *)
(* using ISO Modula-2 *)
(* to illustrate RealMath.sqrt *)
(* last revision 1993 03 01 *)

FROM STextIO IMPORT
  WriteString, WriteLn, SkipLine, ReadChar;
FROM SRealIO IMPORT
  ReadReal, WriteFixed;
FROM SIOResult IMPORT
  ReadResult, ReadResults;
FROM RealMath IMPORT
  sqrt;

VAR
  side1, side2, side3 : REAL;
  key : CHAR;

PROCEDURE GetReal (VAR numToGet : REAL);
VAR
  tempResult : ReadResults;
  
BEGIN
  REPEAT
    WriteString ("Please type in a real number ===> ");
    ReadReal (numToGet);
    tempResult := ReadResult ();
    SkipLine; (* swallow line marker *)
    WriteLn;
  UNTIL tempResult = allRight;
END GetReal;

BEGIN   (* main *)
  WriteString ("What is the first side length? ");
  WriteLn;
  GetReal (side1);
  WriteString ("What is the second side length? ");
  WriteLn;
  GetReal (side2);
  side3 := sqrt (side1 * side1 + side2 * side2);
  WriteString ("The hypotenuse is ");
  WriteFixed (side3, 2, 0);
  WriteString (" units long.");
  WriteLn;
  WriteString ("Press any key to continue");
  ReadChar (key);
END Pythagoras.

Results of a run:

What is the first side length? 
Please type in a real number ===> 20.0
What is the second side length? 
Please type in a real number ===> 21.0
The hypotenuse is 29.00 units long.

6.4.2 Exponential and Logarithmic Functions

It was observed in section 4.9 that an amount A placed at compound interest rate for time years would grow to A(1 + rate)time. Notice that if the interest is compounded n times a year, the amount will be higher than if it is compounded annually, even though the interest rate at each compounding period must be divided by n. For instance, the first example in that section found that $1000 at 6% compounded for 10 years would grow to 1790.85. If the interest is computed monthly, the rate at each application of interest is .06/12, or .005, and the number of applications becomes 10*12 or 120. Modifying the formula above yields A(1 + rate/n)n*time and in this case the $1000 grows to $1819.40. If the compounding is done daily, the result is 1822.20, not much of a difference from monthly compounding. One might ask whether continuing to increase the number of compounding periods indefinitely would yield an indefinitely large (infinite) amount, or if there is some limit beyond which the amount will not grow.

The latter turns out to be the case. To see that this is so, consider a principal amount of $1.00 placed at 100% for a year, and increase the number of compounding periods indefinitely. In mathematical terms, this computes:

As the number of periods grows and the interest rate per period shrinks, this converges to a definite limit at about 2.7182818. (It has a non-repeating, non-terminating decimal representation.) This number, denoted e arises naturally in a variety of situations in mathematics, and mathematical libraries provide the exponential function to compute y = ex. The inverse function, to compute the exponent, or logarithm of x given the number y is the logarithm function x = ln(y). These were mentioned in section 4.5 in connection with writing the function procedure APowerB and are found in RealMath as:

PROCEDURE exp (numb : REAL) : REAL;

and

PROCEDURE ln (numb : REAL) : REAL;

In addition, RealMath exports the related function procedure:

PROCEDURE power (base, exponent: REAL): REAL;

NOTE: Some non-standard mathematical library modules also export some or all of the following related function procedures:

For Base 10 Logarithms:

    PROCEDURE log (numb : REAL) : REAL;
    PROCEDURE TenToX (numb : REAL) : REAL;

For other power and magnitude operations:

    PROCEDURE ipower (numb1 : REAL; numb2 : INTEGER) : REAL;
        (* Both return numb1 to the numb2 power *)
    PROCEDURE Magnitude (numb : REAL) : INTEGER;
        (* returns the order of magnitude of numb, namely the largest integer less than or equal to the scale factor or log10 of numb *)

One application of the logarithmic and exponential functions is to compute radioactive (and other) decay processes. Under normal conditions, a quantity of radioactive material decays over time according to the formula

A = A0 ekt

where A0 is the amount of the substance present at time zero, A is the amount at the time being examined, t is the elapsed time in appropriate units, and k is a constant that is a property of the substance.

In the standard literature, one often finds the constant k expressed indirectly as the half-life, that is, the time it would take for half of any given quantity of the substance to decay.

Problem:

A lab is gathering data from experiments done on radioactive samples and determines experimentally the amount of a radioactive substance present at time zero and also at some subsequent time. Write a program to calculate the half-life of the substance from this data. (This is often one way of identifying an unknown radioactive material.)

Discussion:

The formula A = A0 ekt may be rewritten as

and, upon taking natural logarithms on both sides and solving for k, one obtains

In the case where half of the material is supposed to have decayed, the right hand side of the latter formula becomes

or, solving for t,

With these variations on the initial formula, all the tools are at hand to write the code to do the computation.

Code:

MODULE HalfLife;

(* Written by R.J. Sutcliffe *)
(* using ISO Modula-2 *)
(* to illustrate RealMath.ln *)
(* last revision 1994 08 30 *)

FROM STextIO IMPORT
  WriteString, WriteLn, ReadChar, SkipLine;
FROM SRealIO IMPORT
  ReadReal, WriteFixed;
FROM SIOResult IMPORT
  ReadResult, ReadResults;
FROM RealMath IMPORT
  ln;

PROCEDURE GetReal (VAR numToGet : REAL);
VAR
  tempResult : ReadResults;
  
BEGIN
  REPEAT
    WriteString ("Please type in a real number ===> ");
    ReadReal (numToGet);
    tempResult := ReadResult ();
    SkipLine; (* swallow line marker *)
    WriteLn;
  UNTIL tempResult = allRight;
END GetReal;

VAR
  initialAmount, laterAmount, timePassed,
  constant, halfLife : REAL;
  key: CHAR;

BEGIN
  WriteString ("What was the initial amount? ");
  GetReal (initialAmount);
  WriteString ("How much time elapsed til the second reading? ");
  GetReal (timePassed);
  WriteString ("And, how much material was left then? ");
  GetReal (laterAmount);
  constant := ln (laterAmount / initialAmount) / timePassed;
  halfLife := ln ( 0.5) / constant;
  WriteString ("The half life of this material is ");
  WriteFixed (halfLife, 6, 10);
  WriteLn;
  WriteString ("Press a key to continue ==>");
  ReadChar (key);
END HalfLife.

Trial Run:

What was the initial amount? Please type in a real number ===> 100.0
How much time elapsed til the second reading? Please type in a real number ===> 10.0
And, how much material was left then? Please type in a real number ===> 25.0
The half life of this material is   5.000000
Press a key to continue ==>

Naturally, the units for the half life will be the same as those of the elapsed time given by the user.

Exponential growth works in much the same way, except that the constant k is positive instead of negative.

6.4.3 Trigonometric Functions

It is discovered in elementary geometry that in two similar triangles, the sides are in proportion.

In particular, for right triangles, one may relate these fixed ratios to a particular acute angle, such as angle B in figure 6.8. When this is done, and the right triangle labelled as in figure 6.9, the trigonometric ratios are defined as follows:

The symbols sin, cos, and tan are themselves abbreviations for sine, cosine, and tangent, respectively. Auxiliary to these three are their inverse functions arcsin, arccos, and arctan for producing an angle given one of the fixed ratios. Wirth suggested that MathLib0 provide only three of these functions; the minimum necessary for work in trigonometry, but the ISO library RealMath supplies all six. They are:

PROCEDURE sin (x : REAL) : REAL;
PROCEDURE cos (x : REAL) : REAL;
PROCEDURE tan (x : REAL) : REAL;
PROCEDURE arcsin (x : REAL) : REAL;
PROCEDURE arccos (x : REAL) : REAL;
PROCEDURE arctan (x : REAL) : REAL;

NOTES: 1. The first three require the angle to be in radians, and return the sine, cosine, and tangent, respectively, of the angle supplied. The last three takes the sine, cosine, and tangent of an angle and returns the principal value of the angle measure (in radians). For arcsine and arctangent, the values returned are in the range -/2 < ø < /2. For arccosine, the values returned are in the range 0 < ø < . (Recall that one radian is 180/ degrees.)

2. Many non-standard implementations are much less generous in their supply of trigonometric functions than this, and may omit as many as three of these.

3. While angles larger than 2 (about 6.28) will work correctly, specific implementations may have a maximum argument for trigonometric functions that is much smaller than MAX (REAL).

Here is another useful procedure:

  PROCEDURE degToRad (x : REAL) : REAL;
    (* converts degrees to radians *)  
  BEGIN
    RETURN (x * pi / 180.0); 
  END degToRad;  

Basic trigonometric identities and formulas can be employed to extend the scope of the available mathematical functions.

Example:

Write a module that computes the area of any triangle given two adjacent sides and an included angle.

Discussion:

Consider the triangle ABC in figure 6.10 below. If the base b (here AC) and the altitude h (here BD) is known, the area of the triangle is given by

S = .5bh

However, since BD is the side opposite angle C, in the right triangle BCD, and sin (C) = BD/BC,

h = BC sin (C) = a sin (C).

That is, provided the given data consists of two sides a and b of the triangle and an included angle C, its area can be computed by using the formula

s = .5ab sin (C)

which reduces to the original formula when C is a right angle, because the sine of 90° is one.

Code:

MODULE TriArea;

(* Written by R.J. Sutcliffe *)
(* using ISO Modula-2 *)
(* to illustrate RealMath.sin *)
(* last revision 1993 03 01 *)

FROM STextIO IMPORT
  WriteString, WriteLn, ReadChar, SkipLine;
FROM SRealIO IMPORT
  ReadReal, WriteFixed;
FROM SIOResult IMPORT
  ReadResult, ReadResults;
FROM RealMath IMPORT
  sin, pi;

PROCEDURE GetReal (VAR numToGet : REAL);
(* prompts for a real number, reads it, loops until a correct one is typed, swallows the end-of-line state and returns the number read *)

VAR
  tempResult : ReadResults;
  
BEGIN
  REPEAT
    WriteString ("Please type in a real number ===> ");
    ReadReal (numToGet);
    tempResult := ReadResult ();
    SkipLine; (* swallow line marker *)
    WriteLn;
  UNTIL tempResult = allRight;
END GetReal;

PROCEDURE degToRad (x : REAL) : REAL;
  (* converts degrees to radians *)  
BEGIN
  RETURN (x * pi / 180.0); 
END degToRad; 

VAR
  angleC, sideA, sideB, area: REAL;
  key : CHAR;

BEGIN  (* main *)
  (* obtain triangle data *)
  WriteString ("What is the first side length? ");
  GetReal (sideA);
  WriteString ("What is the second side length? ");
  GetReal (sideB);
  WriteString ("Now, what is the included angle ");
  WriteString ("in degrees? ");
  GetReal (angleC);
  (* do calculation *)
  angleC := degToRad (angleC);
  area := 0.5 * sideA * sideB * sin (angleC);
  (* inform user of result *)
  WriteString ("The area is ");
  WriteFixed (area, 5, 0);
  WriteString (" square units ");
  WriteLn;
  WriteString ("Press a key to continue ==>");
  ReadChar (key);
END TriArea.

First run:

What is the first side length? Please type in a real number ===> 10.0
What is the second side length? Please type in a real number ===> 8.0
Now, what is the included angle in degrees? Please type in a real number ===> 30.0
The area is 20.00000 square units 

Second run:

What is the first side length? Please type in a real number ===> 10.0
What is the second side length? Please type in a real number ===> 20.0
Now, what is the included angle in degrees? Please type in a real number ===> 90.0
The area is 100.0000 square units 

NOTE: Some nonstandard versions of the math library module use atan instead of arctan and may not export asin, or acos. Others provide the hyperbolic functions sinh, cosh and tanh.

6.4.4 Conversions

RealMath also offers the useful conversion function:

PROCEDURE round (x: REAL): INTEGER;

which rounds off a real to the nearest integer.

Two function procedures that may be in the traditionally named MathLib0 that are of more specialized interest are the following integer/real conversions.

PROCEDURE real (m : INTEGER) : REAL;
PROCEDURE entier (x : REAL) : INTEGER;

The first of these is essentially the same as FLOAT except that it only operates on the type INTEGER (and assignment compatible cardinals.) This is important in older versions of Modula-2 where FLOAT works only on CARDINAL (not on either one as in the ISO standard.)

The second one is sometimes called the greatest integer function. It takes a real argument, and returns the greatest integer less than or equal to the real. Note that this is not the same as TRUNC even in those versions where both can return integers. Compare the following:

entier (5.7) produces 5 and TRUNC (5.7) also produces 5, but

entier (-4.3) produces -5 while TRUNC (-4.3) yields -4

That is, for positive numbers, the result is the same, but for negative ones, it will be different, because in those cases entier gives the nearest integer less than the argument and TRUNC simply "hacks off" the decimal fractional portion. Notice that an order of magnitude function would be written using entier rather than TRUNC.

  PROCEDURE Magnitude (num : REAL) : INTEGER;
(* uses non-ISO functions *)

  BEGIN
    RETURN (entier (ln (num) / ln (10.0) ))
  END Magnitude;

This procedure returns -6 when given 4.5E-6 and 2 when given 3.8E2, having computed the base ten logarithms as approximately -5.346 and 2.579 respectively. Notice that a base ten logarithm of a number (or one in any other base) is computed by dividing the natural logarithm of the number by the natural logarithm of the base, for if

x = log10y then 10y = x,

so that, taking natural logarithms on both sides yields

y ln(10) = ln (x)

and therefore

y = ln(10)/ln(x)

as used in the procedure.

6.4.5 Summary of RealMath

DEFINITION MODULE RealMath;

CONST
  pi   = 3.1415926535897932384626433832795028841972;
  exp1 = 2.7182818284590452353602874713526624977572;

PROCEDURE sqrt (x: REAL): REAL;
PROCEDURE exp (x: REAL): REAL;
PROCEDURE ln (x: REAL): REAL;
PROCEDURE sin (x: REAL): REAL;
PROCEDURE cos (x: REAL): REAL;
PROCEDURE tan (x: REAL): REAL;
PROCEDURE arcsin (x: REAL): REAL;
PROCEDURE arccos (x: REAL): REAL;
PROCEDURE arctan (x: REAL): REAL;
PROCEDURE power (base, exponent: REAL): REAL;
PROCEDURE round (x: REAL): INTEGER;
END RealMath.

6.4.6 Other Mathematical functions

A wide variety of other function procedures and error handling may be provided in some auxiliary modules associated with RealMath, or, in non-standard versions, added to MathLib0.

The ISO standard libraries, and some non-standard versions as well, include a second module that is identical to RealMath but that acts on and returns long types.

DEFINITION MODULE LongMath;

CONST
  pi   = 3.1415926535897932384626433832795028841972;
  exp1 = 2.7182818284590452353602874713526624977572;

PROCEDURE sqrt (x: LONGREAL): LONGREAL;
PROCEDURE exp (x: LONGREAL): LONGREAL;
PROCEDURE ln (x: LONGREAL): LONGREAL;
PROCEDURE sin (x: LONGREAL): LONGREAL;
PROCEDURE cos (x: LONGREAL): LONGREAL;
PROCEDURE tan (x: LONGREAL): LONGREAL;
PROCEDURE arcsin (x: LONGREAL): LONGREAL;
PROCEDURE arccos (x: LONGREAL): LONGREAL;
PROCEDURE arctan (x: LONGREAL): LONGREAL;
PROCEDURE power (base, exponent: LONGREAL): LONGREAL;
PROCEDURE round (x: LONGREAL): INTEGER;

END LongMath.

This second module (along with the built-in type LONGREAL itself) are provided because many systems have two or more real types of different precisions. ISO Modula-2 defines the precision of LONGREAL to be equal to or greater than that of REAL. Thus, if there is only one underlying type in the actual system being used, the programmer may use either or both of the Modula-2 logical types to refer to this actual type. Both RealMath and LongMath also include an error enquiry function not listed here but the use of this will be postponed to a later chapter.

It should be noted that the name and the contents of both modules in non-standard versions, can vary widely from one implementation to another. For further information, see the Appendix on standard module definitions or consult the manuals that are available with the system.


Contents