\" geqn -d$$ 5_Algorithms | groff -pt -ms | lpr .fp 4 C .EQ delim $$ gfont R .EN .bp .EH '\fIScientific Computation''Algorithms\fR %' .OH '\fIScientific Computation''Algorithms\fR %' .sp 5 .ps 14 .ce .B 5. Algorithms .R .NL .sp 3 .NH Introduction .PP \fBalgorithms\fR are processes or procedures leading to desired results. The term \fIalgorithm\fR is derived from the name of the ninth-century Persian mathematician Mohammed A-Khow\*^arizm\*^i (Latin, Algorismus), who first described step-by-step procedures for performing arithmetic operations on Arabic decimal numbers (addition, subtraction, multiplication and division). The term \fIalgebra\fR, referring to systems for manipulating symbols according to equivalence rules, has the same root. The idea of formal algorithmic procedures dates to much earlier times, sophisticated examples being Greek geometric constructions and Mayan time reckoning. .PP We employ algorithms for most of our activities. Such procedures many times are comprised of a number of component processes, called \fBsteps\fR. Assembly, operation and repair manuals list step-by-step algorithms, as do income tax forms and cooking recipe books. The procedures for sequential (linear) tasks, such as assembly-line production, and repetitive (looping) tasks, such as "rote" memorization are relatively easy to describe. More complicated processes, such as creative thinking may also be procedural in nature, but the path of activity may be more difficult to discern in terms of component steps. Many mathematical procedures, such as derivations and proofs, are algorithmic in nature, and problem-solving algorithms are at the heart of computer science. .PP \fBHeuristics\fR are unproven strategies based more on intuition than derivation. They are the only choice for situations which don't lend themselves to deterministic analysis. We use heuristics in games of chance stock market investing, and highly complex processes. In practice, heuristics can sometimes outperform well-defined algorithms, even in simple applications. .PP Computer programs consist of sets of data structures and instructions. Translating an algorithm into a computer program is called \fBimplementing\fR the algorithm. Programming languages are often designed to implement specific algorithmic tasks. For example, FORTRAN implements certain mathematical structures and statements in a natural way. Similarly COBOL is designed to implement business transaction structures and statements. While there may be a variety of languages from which to select for the implementation of a given algorithm, programming effort is reduced by selecting programming languages that incorporate data and control structures appropriate to the structure of the algorithm to be implemented. .PP \fIAlgorithmetic analysis\fR consists in verifying that a given algorithm does what is intended and in deducing the expenditure of resources required to execute it. Computer resource expenditures may be measured in \fItime\fR units (e.g. CPU cycles) or \fIspace\fR units (e.g. memory words), or as the integrated product of time and space units. Just as there is usually more than one way to solve a given problem, there are often alternative algorithms which produce the same results from the same input. \fIEfficiency analysis\fR compares alternative algorithms and implementations in terms of absolute space and time resource usage and relative tradeoffs between space and time. Comparisons are usually expressed in terms of required resources as functions of input parameters representing the \fIsize\fR of the problem to be solved. For example, analysis of procedures to sort lists of items into some specified order would be expressed in terms of the number of program steps or number of storage locations needed to perform the sorting as a function of the number of items to be sorted. .PP The \fBorder\fR of an algorithm expresses the dependence of space and/or time on the size of the task. Rather than explore alternative algorithms for a given task in any great detail (a favorite pastime of computer scientists), we will survey a variety of tasks, exploring a few variations and including some known efficient algorithms. .NH Simple Algorithms .SH Balancing a Bank Account .PP Let's begin with a simple (or perhaps not so simple) common algorithm, namely reconciling, or "balancing" a bank account statement. The idea is to start with a given amount of money in an account (the "balance") at a given time, and compute the amount at a later time by addition of deposits and subtraction of withdrawals from the given amount. A mathematical statement of the balancing process is .EQ F~=~I~+~sum T, .EN where I is the initial balance, F the final balance and T is the transaction amount (positive for deposits, negative for withdrawals). .PP We need to start somewhere, so let's assume the first step is "open a bank account". Actually this step presumes some number of preceding steps, such as "first obtain some money", or "first get a job", or "first register with Social Security, or "first establish residency", etc. Before we begin, it would be well to know some arithmetic too. This illustrates the need to establish boundaries between what is given and what is expected. The list of given input and expected output items needs to be carefully defined and well understood. We will describe our expectation with the statement, "given a bank statement and the ability to scan and interpret the statement and perform addition and subtraction, reconcile the transactions with the stated initial and final account amounts." .PP Note that some implementation is implicit in the above expected "abilities". \fIScanning\fR ability could be implemented with human eyes or optical scanners, \fIinterpretation\fR requires communication (interfacing) between the scanner and the interpreter, interpretation and arithmetic \fIoperations\fR could involve human reasoning or hardware computer operations. It might be useful to look ahead and decide how we are going to implement the calculational aspects of the algorithm (pencil and paper, pocket calculator, supercomputer, etc), but the algorithm for the expected balancing task should (and in this case will) be independent of how it is implemented. .PP Given a list of (signed, floating point) numbers the essential act of the balancing algorithm will be addition of numbers representing deposits and subtraction of numbers representing withdrawals. These results only produce the change in amount in the account, however. To reconcile the bank statement, it is necessary to know the initial amount in the account, how many transactions are to be computed and the declared balance after the transactions. .PP A realistic account balancing algorithm would assume there are at least two numbers in an input list, the first number being the total amount in the account before the transactions begin and the last (which could be the next item) being the final amount in the account after any transactions have taken place. Assigning an initial balance equal to the initial input value, a repetitive, \fIlooping\fR process employing addition and subtraction operations would be an appropriate algorithm. Looping would continue until the last entry in the input list is encountered thus producing a final computed balance. Note that determining whether a given entry is the last in the list may not be trivial. The computed balance up to the final entry could be compared with the final entry in the input, and the results reported. .PP Implementation of the algorithm requires additional considerations. Some format for the data would need to be agreed on, such as positive or negative (or possibly parenthesized) decimal numbers representing, say, dollars and cents in the balances and transactions. (Processing transactions in non-decimal monetary systems would obviously be more difficult if the available integer arithmetic operations deal only with decimal integers.) How the entries are distinguished (separated by spaces, tabs commas, lines, etc.) would need to be established. Some error checking on the input could be performed to confirm the input is in the correct format and within reasonable range, and some action taken if the input is suspicious. Humans can sometimes do this almost automatically; computers may be more accurate but less suspicious. .PP Invoking arithmetic addition and subtraction algorithms during the account balancing algorithm illustrates a \fInested\fR algorithmic structure. Program implementations of such algorithms ordinarily employ \fBsubroutines\fR to perform repetitive sub-tasks on variable data. .SH Integer Multiplication .PP Continuing with familiar processes, let's explore some algorithms for multiplying two positive decimal integers. .PP One possible way to proceed is to locate a "multiplication table" of all possible products of integers and look up the product corresponding to the given multiplicands. Needless to say, this becomes impractical for arbitrarily large numbers. Another possibility, based on the definition of multiplication as multiple additions, could proceed by adding one of the integers n times, where n is the other integer. This becomes inefficient for large integers. Note that these are \fInot incorrect\fR, merely \fIimpractical\fR algorithms. .PP A better way to proceed is to break the problem down into smaller problems from which the final result can be accumulated. \*QLong-hand\*U multiplication using a pencil and paper is based on the recognition that an integer may be expanded as products of digits times powers of the base (the base is 10 for decimal integers), and the distributive and associative laws of arithmetic. These properties allow the multiplication of arbitrarily large integers to be reduced to the sum of partial products obtained from the multiplication of two digits. In this way multiplication of decimal integers requires looking up products in a table of products of single digits (containing ${n(n~+~1)} over 2~=~55$ unique entries for $n~=~10$ digits, 0 through 9). For example, .EQ ~~~~~~~~~~~~~~~~~~~~123 times 456~mark =~(1 times 100~+~2 times 10~+~3 times 1) times 456 .EN .EQ lineup =~(3 times 1~+~2 times 10~+~1 times 100) times 456 .EN .EQ lineup =~3 times 456~+~2 times 4560~+~1 times 45600 .EN .EQ lineup =~3 times (4 times 100~+~5 times 10~+~6 times 1) .EN .EQ lineup ~~+~2 times (4 times 1000~+~5 times 100~+~6 times 10~+~0 times 1) .EN .EQ lineup ~~+~1 times (4 times 10000~+~5 times 1000~+~6 times 100~+~0 times 10~+~0 times 1) .EN .EQ lineup =~(12 times 100~+~15 times 10~+~18 times 1) .EN .EQ lineup ~~+~(8 times 1000~+~10 times 100~+~12 times 10~+~0 times 1) .EN .EQ lineup ~~+~(4 times 10000~+~5 times 1000~+~6 times 100~+~0 times 10~+~0 times 1) .EN .EQ lineup =~1368~+~9120~+~45600 .EN .EQ lineup =~56088 .EN .PP The rearrangements don't affect the result (distributivity and associativity properties over addition) and renders the process more familiar to the long-hand multiplication, algorithm which scans one of the numbers from right to left. Also note that it doesn't matter which of the original integers is expanded. Multiplying out the powers of 10 in the third line is equivalent to shifting the second number to the left by one place for each power of 10, which accounts for the shifting operation of sub-products in some long-hand multiplication schemes; filling, or "padding" the spaces produced at the right end with zeros would be more explicit. .PP Generalizing the process outlined in the example yields an algorithm for multiplying arbitrary decimal integers: .sp .nr LL -0.5i .KS .B1 .sp .ce \fBInteger Multiplication Algorithm\fR .sp .5 .LP Purpose: To obtain the product of two arbitrary decimal integer numbers. .LP Procedure: .RS .IP 1) 3 Scan the digits of one of the numbers from right (least) to left (greatest). For each digit of the number scanned, repeatedly: .RS .IP a) 3 Obtain a sub-product by multiplying the digit by the other number. .IP b) 3 Increase the other number tenfold. .RE .IP 2) 3 Sum the sub-products to obtain the result. .IP 3) 3 Assign the sign of the product to be the product of the signs of the numbers. .RE .sp .B2 .KE .nr LL +0.5i .PP This algorithm illustrates a number of general features of algorithms. Notice the use of nested operations in step 1 (steps a and b); this happens with repetitive processes. Note also that the multiplication of an integer by a digit in step 1a) may involve \*Qcarrying\*U digits when a subproduct exceeds one digit. Thus, multiplication of a multidigit number by a single digit number involves yet another process, requiring another (nested) algorithm.\** .FS Actually this would involve a similar integer expansion in products of digits and powers of the base and accumulation of subproducts. .FE It is not uncommon for one algorithm to apply sub-algorithms in its steps. In fact, reducing a process to its most fundamental steps could involve many layers of subprocesses. .PP Digital computers can do the equivalent of long-hand integer multiplication efficiently by working in a binary base representation, which has only two digits, 0 and 1 (called \fIbi\fRnary digi\fIts\fR, or \fIbits\fR). Arbitrary integers are represented as binary numbers with a string of binary digits (0 or 1), representing a sum of binary digits (0 or 1) times powers of the base (2), similar to decimal numbers. Scanning the digits of one number (step 1) is effected by repeatedly reading the rightmost digit (bit) and shifting the number one place (bit) to the right. Multiplication of a given number by a binary digit (step 1a) simply produces 0 (if the binary digit is 0), or reproduces the given number (if the binary digit = 1). Increasing the other number by powers of the base (step 1b) is effected by shifting to the left (and filling in a zero). Step 2 becomes binary addition. Binary integer multiplication is very efficient because binary operations (the subprocesses) are very basic, efficient process in digital computers, and therefore among the fastest operations. .NH Structured Programming .PP The main point of the discussion so far is to illustrate how much care and planning are needed to formalize even ordinary processes. Whenever we analyze a procedure to the point of describing it in terms in terms of simpler processes, we find there must be some assumptions made and some attempt to communicate in the language of the simpler processors. Complex processes can often be broken down into simpler tasks such that the entire process is a composite of the smaller parts. \fBStructured programming\fR refers to the organization of the programming data and instructions to perform complex tasks in terms of simpler structures. .PP The order in which multiple-step processes are carried out is referred to as \fBcontrol structure\fR. There are three basic control structures, \fIsequential, branched and looping\fR. The terminology connotes a graphical representation of the path of action as a function of time, with intersections and vertices representing the steps, connected by lines indicating the order of operations. Fig. 4.1 shows an example of the graph of a typical control structure. .KS .PS line from 5.487,9.762 to 5.487,8.762 line from 5.487,8.762 to 5.050,8.325 line from 5.050,8.325 to 5.487,7.825 line from 5.487,7.825 to 5.925,8.325 line from 5.925,8.325 to 5.487,8.762 line from 5.925,8.325 to 6.362,8.325 line from 6.362,8.325 to 6.362,9.387 line from 6.362,9.387 to 5.487,9.387 line from 5.487,7.825 to 5.487,7.075 box with .sw at (2.99,6.01) width 5.06 height 4.00 "\s14\fRs\fP" at 4.800,8.301 ljust "\s14\fRs\fP" at 5.175,7.801 ljust "\s10\fR2\fP" at 5.237,7.665 ljust "\s12\fRFig. 4.1. Test statements are indicated by t, other execution statements by s.\fP" at 3.362,5.421 ljust "\s12\fR1\fP" at 3.737,9.421 ljust "\s12\fR2\fP" at 3.737,9.171 ljust "\s12\fRs\fP" at 3.612,9.046 ljust "\s12\fR1\fP" at 3.737,8.921 ljust "\s12\fRs\fP" at 3.300,8.796 ljust "\s12\fR2\fP" at 3.362,8.671 ljust "\s14\fRF\fP" at 5.737,8.614 ljust "\s10\fR2\fP" at 5.300,8.728 ljust "\s14\fRt\fP" at 5.237,8.801 ljust "\s14\fRt\fP" at 5.175,9.426 ljust "\s14\fRT\fP" at 5.112,8.551 ljust "\s10\fR1\fP" at 4.862,8.228 ljust "\s10\fR1\fP" at 5.237,9.353 ljust "\s12\fRwhile(t )\fP" at 3.300,9.546 ljust "\s12\fRif (t )\fP" at 3.487,9.358 ljust .PE .KE .sp 2 .PP Algorithm control structures are expressed in English by keywords, such as \fIdo\fR something \fInext\fR, \fIif\fR something \fIthen\fR, and \fIwhile\fR something \fIrepeat\fR. More generally, any action verb can appear in place of \fIdo\fR, as can a period (.) or semicolon (;) in place of \fInext\fR; \fIcheck\fR might be used in place of \fIif\fR, \fIuntil\fR in place of \fIwhile\fR, etc. .PP There are \fIstorage structure\fR analogues to \fIcontrol structures\fR. Linear sequential data structures can be organized into homogeneous \fIarrays\fR and heterogeneous \fIrecords\fR, called \fIstructures\fR in C. Alternative data structures are called \fIunions\fR in C. Looping control structures are analogous to \fIrepeated subroutine (function) calls\fR. .NH Recursion .PP Another control structure which is natural to a certain class of algorithms involves \fInested\fR operations related to the mathematical construct of functions of functions, called \fBfunctional composition\fR. Functional composition may be represented with the notation f(g(h(...))). If the composed functions have the same functional form, the total is called a \fBrecursive function\fR,\** $f(f(f~...~f(x)))~=~f sup n (x)$. .FS Note that n refers to the nesting depth, not to a power as in algebra, or to an order as with differential equations. .FE In programming languages, functions which refer to ("call") themselves are also called \fIrecursive functions\fR. Structured data types can be defined that refer to themselves as well. .PP Not all programming languages permit recursive calls and data structures, as recursion implementation requires some mechanism to store intermediate results that can be passed back to the calling functions and referencing data structures. Temporary storage identified by the language compiler and created dynamically during execution is called an "activation record" and may be implemented during execution with a "run-time stack". .PP When a problem can be solved by reducing the problem to smaller problems of the same type, a "divide and conquer" strategy is suggested, which may be implemented with recursive function calls. To avoid reduction \fIad infinitum\fR, the recursion must stop at some point. Thus recursive programs always \fItest\fR some condition (usually related to the size of the input) and "exit" (recursively return) when the test is passed. .PP It should be noted that most tasks are \fInot\fR recursive in nature. Furthermore, any recursive process or structure can be cast, or "derecursed" into an equivalent non-recursive process or structure. This justifies including programming languages which do not support recursive structures in the class of "general purpose" languages. The power of recursion lies in the ease in implementing naturally recursive structures with those languages which support recursive subroutine calls and data structures. .NH Algorithm Efficiency .PP The most efficient algorithms are independent of the size of the problem; these are order(1) algorithms, because the resources needed are a function of $n sup 0 ~=~1$. .PP Looping operations require resources proportional to the number of repetitions. If the number of repetitions is proportional to the size of the problem, n, the efficiency of looping portion of the algorithm is of order(n). .PP Recursive algorithms that solve successively smaller problems require space/time resources proportional to $log(n)$, where n is a measure of the size of the original problem.\** Such algorithms may be shown to be of order $log(n)$. .FS Note that changing the base of the logarithm only changes the proportionally constant, not the basic functional dependence on size, or order. .FE When the results of smaller problems are recombined recursively to obtain the solution to the original problem, the total number of operations, or \fIwork\fR becomes proportional to $n~log(n)$, and the algorithm is of order($n~log(n)$). .PP Alternative algorithms which accomplish the same task can trade space for time, so a fair comparison between algorithms should consider both space and time requirements. Nevertheless, \fIit is common practice to measure efficiency in execution time\fR, which is roughly proportional to the number of program steps, discounting startup, or "overhead" charges. Order(1) algorithms are thus commonly called "constant time" algorithms. .PP As an example of these efficiency concepts, consider the game "I'm thinking of a number between one and N", where N is some maximum value, say one ten, one hundred or one thousand. The person guessing the number may employ a looping iterative approach, perhaps starting with the minimum possible value (one) and incrementing by unity until the guess matches the value selected by the person thinking of the number. This algorithm obviously takes longer on average to match the answer the greater N is, and is of order N. An alternative strategy is to divide the guessing range in half and determine which half the answer lies in. This requires some cooperation on the part of the person posing the problem, but if obtained, divides the range in half at each guess, and the answer is obtained within M guesses for any case, where $2 sup M ~=~N$, producing a lg(N) order process (lg symbolizes the logarithm in base 2), typical of divide-and-conquer algorithms. .NH Recursive Arithmetic Algorithms .PP Algorithms related to mathematics are among the most fundamental. We will concentrate on algorithms which deal with numbers, as distinct from (semantic) algorithms which deal with symbols. Examples include mathematical functions, such as the sine and logarithm functions, which are supplied in mathematical libraries supporting programming languages. Routines to compute these functions accurately and efficiently for all legal arguments are usually not trivial. Here we consider some processes which gain power when cast in recursive form. .SH Integer Multiplication .PP \fIArithmetic integer multiplication\fR is defined as repeated, or multiple addition: .EQ m~times~n~==~n~times~n~times~n~times~...~(m~"times")~==~m~times~m~times~m~times~...~(n~"times")~=~prod from i=1 to m n~=~prod from i=1 to n m .EN For comparison purposes, we will show an \fIiterative\fR C program function for integer multiplication followed by a \fIrecursive\fR C program function. .sp 2 .KS .ps -2 .nf \f4 .eo product (m,n) /* Iterative positive integer multiplication */ int m, n; { int prod; for (prod = 0; n > 0; --n) prod += m; return (prod); } product (m,n) /* Recursive positive integer multiplication */ int m, n; { if (n == 1) return (m); else if (n > 0) return (m + product (m, --n)); } .ec .R .fi .ps +2 .KE .PP Note the trade for explicit storage in the variable \f4prod\fR in the iteration program for implicit storage of the variable \f4n\fR in the recursive program stored (in the activation record). Recursion is somewhat contrived in this case since integer multiplication is iterative by nature. However, recursion is based on the property: .EQ m times n~mark =~prod from {i~=~1} to n m~=~m prod from {i~=~1} to n 1 .EN .EQ lineup =~m(1~+~prod from {i~=~1} to {n~-~1} 1 ) .EN .EQ lineup =~m~+~m(1~+~prod from {i~=~1} to {n~-~2} 1 ) .EN .EQ lineup =~.~.~. .EN .PP This is a rather trivial programming example because most languages support the arithmetic multiplication operations (symbolized by, say, the star operator, \(**) provided by the CPU hardware. .SH Factorials .PP Integer factorials are defined as the multiple product: .EQ n!~==~n~times~(n~-~1)~times~(n~-~2)~times~...~times~1 .EN This example is more practical because, on the one hand, most general purpose compilers don't support a primitive factorial function, but, on the other hand, the factorial function is easy to implement. Note the naturally recursive nature of the factorial function: .EQ n !~==~n (n~-~1) ! .EN .EQ 0 !~==~1 .EN Thus, the alternative definitions of the factorial function lend themselves readily to iterative and iterative and recursive implementations. .sp .KS .ps -2 .nf \f4 .eo long factorial (n) /* Iterative positive integer factorial function */ int n; { long fact; for (fact = 1; n > 0; --n) fact *= n; return (fact); } long factorial (n) /* Recursive positive integer factorial function */ int n; { if (n == 0) return (1); else if (n > 0) return (n * factorial (n - 1)); } .ec .R .fi .ps +2 .KE .PP Note the similarity to the integer multiplication implementations. Also note the differences. One might be tempted to decrement n with \f4--n\fR in place of \f4n-1\fR in the last line of the recursive factorial function, but that would be a mistake in this case because \f4n\fR appears as an operand in two places in the last expression, and since the definition of the C language does not guarantee the order of evaluation of operands in expressions, the value stored in the variable \f4n\fR could change when it should not. one only wants a different value of \f4n\fR store in each of the activation records for each time the \f4factorial\fR function calls itself. .PP Although the factorial function is naturally recursive, it is more efficient to implement the iterative algorithm in this case, due to the large overhead in creating the multiple activation records for the recursive calls. .SH Euclid's Algorithm .PP The earliest identified mathematical algorithm beyond arithmetic, attributed to Euclid about 300 BC, but known to the Chinese about the same time (Yun, 1987), deals with finding the \fIgreatest common divisor\fR (gcd) of two integers, which has the practical application of "reducing" (through division by the gcd) an integer fraction (ratio of two integers) to its "lowest terms" (smallest integers). .PP An inefficient, "brute-force" method would test every integer, from unity to the smaller of the given numbers, keeping track of the largest discovered common divisor along the way. The amount of work in such a scheme in the worst case (where the largest common even divisor is unity) is clearly proportional to the magnitude of the smaller integer. Some efficiency might be gained by starting with the smaller of the two given integers, if they differ, and test that integer and every smaller positive integer until one is found that evenly divides both of the given integers, but it can be seen that in the worst case the number of computations needed to reach a conclusion would again be proportional to the magnitude of the input integer. .PP An implementation of the testing phase of the algorithm would need some method of determining whether the quotient of two integers has a remainder or not ("evenness" means there is no remainder). This can be accomplished using a remainder, or \fImodulo\fR (mod) function, available as a standard binary operator in most programming languages that support integer arithmetic. .PP Euclid's method, which is much more clever than a simple search, is based on the fact that the gcd that will evenly divide two integers will also divide their \fIdifference\fR. Repeated differences equate to integer division, with the modulo of the two integers equating to the remainder. If m and n are the two given integers and $m~>=~n$, then either $n~=~0$, in which case m is the greatest common divisor (because a divides itself evenly and 0 is divided by any value without remainder), or $n~>~0$, in which case n divides n evenly and m divides n with remainder $(m~mod~n)~=~r$, with $n~>=~r$. Replacing m with n and n with r yields a similar situation: either $r~=~0$, in which case n is the gcd, or $r~>~0$, etc. A C program implementing this method is given below. .sp .KS .ps -2 .nf \f4 .eo /* Input two positive integers and output their greatest common divisor */ main() { int a, b; printf ("Input two positive integers: "); scanf ("%d%d", &a, &b); if (a > 0 && b > 0) printf ("%d is the greatest common divisor of %d and %d.\n", gcd(a,b), a, b); else printf ("\nIntegers must be positive!\n"); } /* Euclid algorithm for greatest common divisor of two positive integers */ int gcd(a,b) int a, b; { if (b == 0) return a; else return gcd(b, a % b); /* Interchanges a and b for the case a < b */ } .ec .R .fi .ps +2 .KE .PP Euclid's algorithm works more efficiently than simple searching because it successively solves much smaller-sized problems than does the simple decrement method. The analysis of Euclid's recursive algorithm is not trivial but has been shown to involve approximately $left ( {12~ln (2)} over {pi sup 2} right ) ln (n)$ steps on average, where n is the larger of the two input (decimal) integers. This is much more efficient than brute-force decremental testing, which requires n/2 steps on average (an order(n) algorithm), as can be illustrated by comparing the ratios of the order formulas for pairs of values of n, say 1 and 1000. .SH Random Number Generators .PP Numbers which appear not to have any identifiable relationship are said to be \fBrandom\fR. The key word is \fIappear\fR, for not only is it difficult to devise criteria for even testing the lack of a functional relationship, in a practical sense any sequence of random numbers generated with a computer must follow some deterministic (reproducible) algorithm producing numbers from a given initial value ("seed") in some limited range. Hence such numbers are referred to as "pseudo-random" (or false random). Furthermore, if at any point in the internal generation of the sequence some prior value is regenerated, the sequence enters an unavoidable repeating cycle.\** .FS The internal representation of numbers is distinguished from the external presentation of the numbers obtained according to some decoding scheme. If there is a many-to-one mapping from the internal numbers to the external numbers, several different internal numbers can produce the same output value, thereby producing the desirable property that the occurrence of a given random number in a sequence of random numbers occurs with some probability (or fraction). .FE Any computer with a finite word size must thus be limited in the number of possible patterns which can be expressed before a cycle is initiated. The object of random number generating algorithms is to produce as many patterns as possible in an unordered fashion before the deterministic cycle recurs, however "unordered" may be defined. .PP The study of random number generating algorithms and implementations is very broad, ranging from mathematical number theory to cryptography to chaos theory. Needless to say, it is not a closed subject, necessitating more art than science in practical implementations. Yet random numbers are used in science from statistics to stochastic simulation, so they merit some attention. Historically, the "Monte Carlo" method of statistical analysis was invented by John Von Neuman and Nickolas Metropolis during the second world war to simulate neutron flows in thermonuclear devices. It now finds use in a wide variety of simulation applications. .PP An often used algorithm for generating pseudo-random numbers is based on the so-called \fIlinear congruential formula\fR, $a~<-~(a~*~b~+~1)~mod~m$, which for integer values of a, and integer constants b and m, replaces a with a linear function of a mapped onto [0, m-1] by way of the modulus function. The parameter m should be large (but limited to the maximum integer that can be computed) to ensure that many different numbers are generated. The parameter b should be of medium size, and by empirical exploration is recommended to end with x21, with x even (remember, this is art, not science). The first value of a is called the random number generator "seed". .PP The limitation of finite word sizes in digital computers limits the number of different patterns that can be represented, and consequently the number of random paths through the possibilities. Eventually a previously computed value must be encountered, and the path is retraced. Thus finite-state machines ultimately produce cycles of generated values, called "pseudo" random numbers. Hopefully, the cycle length is long enough to be useful in simulating random paths before repetition occurs. There are no known guaranteed strategies to prevent short cycles, but obviously longer word length is necessary, if not sufficient, to produce long cycles. .PP There are connections between random number generators and mathematical number theory and physical dynamics. A cycle length of one represents mathematically a repeating fractional decimal, such as .111..., and physical equilibrium (fixed point). An n-cycle represents a repeating decimal of a certain length, such as .1210...1210...1210, and physical oscillating systems. An infinite cycle never repeats and represents mathematical irrational numbers, such as $pi$, $sqrt 2$, etc, and physical chaotic motion. .sp 2 .KS .ps -2 .nf \f4 .eo static int a; /* Print n random real decimal proper fractions */ main() { int n; extern float random_real(); printf ("Input seed and number of random numbers desired: "); scanf ("%d%d", &a, &n); while (n-- > 0) printf ("%f\n", random_real()); } /* Return a linear congruential random decimal fraction between 0 and 1 */ float random_real() { int b = 2221, m = 65535; return ((a = (a * b + 1) % m) / (float) m); /* Division by m normalizes */ } .ec .R .fi .ps +2 .KE .PP The use of the C external static storage class for the random variable a allows it to be stored between calls to random_real(). Before the value is returned the random integer a is normalized by m and cast to a float to map it into the range (0,1). The value of m is the maximum value of a 16-bit unsigned integer ($2 sup 16~-~1$, starting with 0), and is implementation dependent. Reference 2 at the end of the chapter gives an extensive discussion on random number generators. .NH Sorting .PP Searching for a given item among a collection of items is a very common process to many tasks. The search is facilitated more efficiently if the items are stored in some order, such as words in a dictionary. Using the divide-and-conquer strategy of repeatedly dividing an ordered list of items to be searched into two halves ("binary search") and determining which half the item searched for resides, it can be shown that the maximum work to find an item is order(n log n). .PP Here we will discuss the process of arranging lists of items having some order property into strict order, called \fBsorting\fR. Sorting can be used to prepare a list for efficient multiple searches, and has intrinsic value as well. The sorting process generally involves two basic operations, \fIcomparison\fR and \fIexchange\fR (at least if the sorted items are to occupy the same space as the original unsorted list). It may be useful to practice with a hand of playing cards to understand the algorithms to be discussed. The cards should be face down and comparison tests made by turning cards (temporarily) face up. This aids algorithmic thinking by preventing efficiency improving exchanges obvious to a more global scanning and processing device like the human brain. (While a human may be able to beat a computer in sorting relatively short lists using locally optimized temporary processes, the computer programmed with systematic efficient sorting algorithms will eventually win on sufficiently long lists.) We will use arrays of positive integers for our items to be sorted, to concentrate more on the process than the data structure. The principles can be extended to more complicated items and structures as needed. .SH Selection Sorting .PP There are many elementary algorithms for sorting items, from which we will choose one that is relatively easy to understand and program. In the process called \fIselection sorting\fR, a list of n items to be sorted into increasing order is scanned for the minimum-value item, (or maximum-value item if it is to be sorted in decreasing order), which is then exchanged with the first item. After the smallest (largest) item has been placed in the first position, the search/replace process is repeated on the remaining list of n - 1 items, which are put successively into the remaining positions. This continues until the last item is encountered (which does not need to be sorted). A C function to implement this algorithm is given below. .sp 2 .KS .ps -2 .nf \f4 .eo # define EXCHANGE(x,y) { float temp; temp = *x; *x = *y; *y = temp; } select_sort(a, n) /* Selection sort array a[1] ... a[n] into increasing order */ float a[]; long n; { long i, j, min; for (i = 1; i < n; ++i) { for (min = i, j = i + 1; j <= n; j++) /* find minimum element beyond i */ if (a[j] < a[min]) min = j; EXCHANGE(&a[i], &a[min]); /* move minimum element to position i */ } } .ec .R .fi .ps +2 .KE .PP Note the use of the C macro feature to define the exchange function. Pointers (addresses) are used to reference the values to be exchanged because only copies of the arguments are passed to functions ("call by value"). If only the values were passed to the exchange function, they would be unchanged on return. This is a conscious design feature to protect the integrity of data passed as arguments to functions. .PP The double-nested \f4for\fR loop over n items requires a total of n(n-1) iterations on n items, so this algorithm is order($n sup 2$). .SH Merge Sorting .PP More efficient sorting algorithms require only order(n log n) operations. We present one here, which superficially doesn't require many steps, but in fact passes much of the work on to the compiler through recursive calls. Merge sorting involves two subprocesses. The list of items to be sorted is divided into two parts (roughly halves), each of which is sorted and then the two parts are \fImerged\fR, or interleaved together in sorted order. This is repeated recursively until only two items are to be sorted. The division and recursive merging is easily implemented through a double recursive call to the merge function. .sp 2 .KS .ps -2 .nf \f4 .eo float *merge_sort (a, n) /* recursive array merge sort */ float *a; long n; { if (n > 1) merge (merge_sort(a, n/2), n/2, merge_sort(a+n/2, n-n/2), n-n/2); return a; } merge (a, m, b, n) /* merge contiguous sorted arrays a and b into sorted a */ float *a, *b; long m, n; { long i=0, j=0, k=0; float *c; c = (float *) calloc ((unsigned) (m+n), sizeof(a)); /* c is size of a + b */ while (i < m && j < n) /* move smaller value from a or b to c */ c[k++] = a[i] <= b[j] ? a[i++] : b[j++]; while (i < m) /* copy any remaining a values into c */ c[k++] = a[i++]; for (i = 0; i < k; ++i) /* copy sorted part of c back to a */ a[i] = c[i]; /* (any remaining b values are contiguous to a) */ } .ec .R .fi .ps +2 .KE .sp .PP Temporary storage space is needed by the merge function to store the result of the merged sublists. This is obtained during execution "dynamically" using the \f4calloc\fR memory allocation function. To be considerate, temporary allocated memory should be freed after use with the \f4cfree(pointer)\fR function. The program is deceptively simple in that attending to the index values and their ranges is tricky. .PP Derecursing the merge sort into an equivalent iterative algorithm would be a tedious and can be acomplicated process. Merge sorting nicely illustrates the advantage of being able to employ recursive function calls. .SH References .IP 1. \fIAlgorithms\fR, by Robert Sedgewick, Addison-Wesley Publishing Company, Reading, MA, 1983. .IP 2. \fIThe Art of Computer Programming. Volume 2: Seminumerical Algorithms\fR, 2nd. ed., Donald E. Knuth, Addison-Wesley Publishing Company, Reading, MA 1981. .IP 3. \fIAlgorithmics: The Spirit of Computing\fR, by David Harel, Addison-Wesley Publishing Company, Reading, MA, 1987. .bp .SH Exercises .IP \ 1. Describe an algorithm to perform long-hand (pencil and paper) integer addition. .IP \ 2. Describe an algorithm to perform long-hand (pencil and paper) integer multiplication. .IP \ 3. Describe an heuristic for guessing a number in a given range in the least number of trials. .IP \ 4. Describe an heuristic to never loose at tic-tac-toe. .IP \ 5. Describe an algorithm to solve cryptograms, like .TS center; r. SEND +MORE _ MONEY .TE .IP \ 6. Write, compile and execute a computer program to balance a bank account. .IP \ 7. How does the algorithm for integer multiplication in the section on Simple Algorithms differ from that in the section on Recursive Arithmetic Algorithms? .IP \ 8. Write a computer program to determine the greatest common divisor of two positive integers by the brute-force method and exercise it to determine its efficiency (how execution time depends on input size). (The \f4/bin/time\fR command may be useful on UNIX systems.) .IP \ 9. Exercise Euclid's method for finding the greatest common divisor of two positive integers to determine its efficiency (how execution time depends on input size). Compare the efficiency for Euclid's algorithm with the that of the brute-force algorithm (exercise 2). .IP 10. To see an odd case, run the random_real function using the parameter values b = 19, m = 381 and the seed (initial a) value 0. .IP 11. What is the expected value of the average of many random decimal fractions between 0 and 1? How should the numbers be distributed? Write programs to test the quality of the random numbers generated by the random_real function. .IP 12. Insertion sorting repeatedly scans the unsorted portion of a list to find the smallest (or largest) element and inserts it into the boundary between the the sorted and unsorted portions of the list. Develop the insertion sorting algorithm and implement an insertion_sort(a,n) function, where a is a floating array of n elements. How does the observed number of comparisons in insertion sorting depend on the number of items in the list (i.e. what is the order of the algorithm), in the case of an inverted-order (worst-case) input list? .IP 13. Write driver programs for the select_sort() and merge_sort() functions and test their relative efficiencies by executing them on random-valued arrays of various lengths (number of items). .OH '''' .EH ''''