1COMMENT
2
3                  REDUCE INTERACTIVE LESSON NUMBER 4
4
5                         David R. Stoutemyer
6                         University of Hawaii
7
8
9COMMENT This is lesson 4 of 7 REDUCE lessons.  As before, please
10refrain from using variables beginning with the letters F through H
11during the lesson.
12
13In theory, assignments and LET statements are sufficient to accomplish
14anything that any other practical computing mechanism is capable of
15doing.  However, it is more convenient for some purposes to use
16function procedures which can employ branch selection and iteration as
17do most traditional programming languages.  As a trivial example, if
18we invariably wanted to replace cotangents with the corresponding
19tangents, we could input:;
20
21algebraic procedure cotan(x); 1/tan(x);
22
23COMMENT As an example of the use of this function, we have;
24
25cotan(log(f));
26
27pause;
28
29COMMENT Note:
30
31   1.  The procedure definition automatically declares the procedure
32       name as an operator.
33   2.  A procedure can be executed any time after its definition,
34       until it is cleared.
35   3.  Any parameters are dummy variables that are distinct from any
36       other variables with the same name outside the procedure
37       definition, and the corresponding arguments can be arbitrary
38       expressions.
39   4.  The value returned by a procedure is the value of the
40       expression following the procedure statement.
41   5.  The function COT is already defined in REDUCE and should not be
42       redefined.
43
44We can replace this definition with a different one:;
45
46algebraic procedure cotan(y); cos(y)/sin(y);
47
48g1 := cotan(log(f));
49
50COMMENT In place of the word ALGEBRAIC, we can optionally use the word
51INTEGER when a function always returns an integer value, or we can
52optionally use the word REAL when a function always returns a
53floating-point value.  (ALGEBRAIC can also be omitted, since it is the
54default procedure type.)
55
56Try writing a procedure definition for the sine in terms of the
57cosine, then type G1.;
58
59pause;
60
61COMMENT Here is a more complicated function which introduces the
62notion of a conditional expression:;
63
64algebraic procedure sumcheck(aj, j, m, n, s);
65   COMMENT J is an indeterminate and the other parameters are
66      expressions.  This function returns the global variable named
67      PROVED if the function can inductively verify that S equals the
68      sum of AJ for J going from M through N, returning the global
69      variable named UNPROVED otherwise.  For the best chance of
70      proving a correct sum, the function should be executed under the
71      influence of ON EXP, ON MCD, and any other user-supplied
72      simplification rules relevant to the expression classes of AJ
73      and S;
74   if sub(j=m,aj) - sub(n=m,s) neq 0 or
75      s + sub(j=n+1,aj) - sub(n=n+1,s) neq 0 then unproved
76    else proved;
77
78on exp, mcd;
79
80clear x, j, n;
81
82sumcheck(j, j, 1, n, n*(n+1)/2);
83
84sumcheck(x^j, j, 0, n, (x^(n+1)-1)/(x-1));
85
86COMMENT Within procedures of this sort a global variable is any
87variable which is not one of the parameters, and a global variable has
88the value, if any, which is current for that name at the point from
89where the procedure is used.;
90
91pause;
92
93COMMENT Conditional expressions have the form
94
95   IF condition THEN expression1 ELSE expression2.
96
97There are generally several equivalent ways of writing a conditional
98expression.  For example, the body of the above procedure could have
99been written
100
101   IF SUB(J=M,AJ) - SUB(N=M,S) = 0 AND
102      S + SUB(J=N+1,AJ) - SUB(N=N+1,S) = 0 THEN PROVED
103    ELSE UNPROVED.
104
105Note how we compare a difference with 0, rather than comparing two
106nonzero expressions, for reasons explained in lesson 3.
107
108As an exercise, write a procedure analogous to SUMCHECK for proving
109closed-form product formulas, then test it on the valid formula that
110COS(N*X) equals the product of COS(J*X)/COS(J*X-X) for J ranging from
1111 through N.  You do not need to include prefatory comments describing
112parameters and the returned value until you learn how to use a text
113editor.;
114
115pause;
116
117COMMENT Most REDUCE statements are also expressions because they have
118a value.  The value is usually 0 if nothing else makes sense, but I
119will mention the value only if it is useful.
120
121The value of an assignment statement is the assigned value.  Thus a
122multiple assignment, performed right to left, can be achieved by a
123sequence of the form
124
125    variable1 := variable2 := ... := variableN := expression.
126
127Moreover, assignments can be inserted within ordinary expressions such
128as X*(Y:=5).  Such assignments must usually be parenthesized because
129of the low precedence of the assignment operator, and excessive use of
130this construct tends to make programs confusing.;
131
132pause;
133
134COMMENT REDUCE treats as a single expression any sequence of
135statements preceded by the pair of adjacent characters << and followed
136by the pair >>.  The value of such a group expression is the value of
137the last statement in the group.
138
139Group expressions facilitate the implementation of tasks that are most
140easily stated as a sequence of operations.  However, such sequences
141often utilize temporary variables to count, hold intermediate results,
142etc., and it is hazardous to use global variables for that purpose.
143If a top-level REDUCE statement or another function directly or
144indirectly uses that variable name, then its value or its virgin
145indeterminate status there might be damaged by our use as a temporary
146variable.  In large programs or programs which rely on the work of
147others, such interference has a non-negligible probability, even if
148all programmers agree to the convention that all such temporary
149variables should begin with the function name as a prefix and all
150programmers attempt to comply with the convention.  For this reason,
151REDUCE provides another expression-valued sequence called a
152BEGIN-block, which permits the declaration of local variables that are
153distinct from any other variables outside the block having the same
154name.  Another advantage of using local variables for temporary
155variables is that the perhaps large amount of storage occupied by
156their values can be reclaimed after leaving their block.;
157
158pause;
159
160COMMENT A BEGIN-block consists of the word BEGIN, followed by optional
161declarations, followed by a sequence of statements, followed by the
162word END.  Within BEGIN-blocks, it is often convenient to return
163control and possibly a value from someplace other than the end of the
164block.  Control and a value may be returned via a RETURN-statement of
165the form
166
167          RETURN expression
168or
169          RETURN,
170
1710 being returned in the latter case.  A BEGIN-block does not return
172the value of the last statement.  If a value is to be returned then
173RETURN must be used.  These features and others are illustrated by the
174following function:;
175
176pause;
177
178algebraic procedure limit(ex, indet, pnt);
179   begin COMMENT This function uses up through 4 iterations of L'Hospital's
180      rule to attempt determination of the limit of expression EX as
181      indeterminate INDET approaches expression PNT.  This function is
182      intended for the case where SUB(INDET=PNT, EX) yields 0/0,
183      provoking a zero-divide message.  This function returns the
184      global variable named UNDEFINED when the limit is 0 dividing an
185      expression which did not simplify to 0, and this function
186      returns the global variable named UNKNOWN when it cannot
187      determine the limit.  Otherwise this function returns an
188      expression which is the limit.  For best results, this function
189      should be executed under the influence of ON EXP, ON MCD, and
190      any user-supplied simplification rules appropriate to the
191      expression classes of EX and PNT;
192   integer iteration;
193   scalar n, d, nlim, dlim;
194   iteration := 0;
195   n := num(ex);
196   d := den(ex);
197   nlim := sub(indet=pnt, n);
198   dlim := sub(indet=pnt, d);
199   while nlim=0 and dlim=0 and iteration<5 do <<
200      n := df(n, indet);
201      d := df(d, indet);
202      nlim := sub(indet=pnt, n);
203      dlim := sub(indet=pnt, d);
204      iteration := iteration + 1 >>;
205   return (if nlim=0 then
206              if dlim=0 then unknown
207              else 0
208           else if dlim=0 then undefined
209           else nlim/dlim)
210   end;
211
212% Examples follow...
213pause;
214
215g1 := (e^x-1)/x;
216
217% Evaluation at 0 causes a zero denominator error at top level but
218% continue anyway.
219sub(x=0, g1);
220
221limit(g1, x, 0);
222
223g1:= ((1-x)/log(x))^2;
224
225% Evaluation at 1 causes a zero denominator error at top level but
226% continue anyway.
227sub(x=1, g1);
228
229limit(g1, x, 1);
230
231COMMENT  Note:
232
233   1.  The idea behind L'Hospital's rule is that as long as the
234       numerator and denominator are both zero at the limit point, we
235       can replace them by their derivatives without altering the
236       limit of the quotient.
237   2.  Assignments within groups and BEGIN-blocks do not automatically
238       cause output.
239   3.  Local variables are declared INTEGER, REAL, or SCALAR, the
240       latter corresponding to the same most general class denoted by
241       ALGEBRAIC in a procedure statement.  All local variables are
242       initialized to zero, so they cannot serve as indeterminates.
243       Moreover, if we attempted to overcome this by clearing them, we
244       would clear all variables with their names.
245   4.  We do not declare the attributes of parameters.
246   5.  The NUM and DEN functions respectively extract the numerator
247       and denominator of their arguments.  (With OFF MCD, the
248       denominator of 1+1/X would be 1.)
249   6.  The WHILE-loop has the general form
250
251          WHILE condition DO statement.
252
253       REDUCE also has a "GO TO" statement, and using commas rather
254       than semicolons to prevent termination of this comment, the
255       above general form of a WHILE-loop is equivalent to
256
257          BEGIN  GO TO TEST,
258       LOOP: statement,
259       TEST: IF condition THEN GO TO LOOP,
260          RETURN 0
261          END.
262
263       A GOTO statement is permitted only within a block, and the GOTO
264       statement cannot refer to a label outside the same block or to
265       a label inside a block that the GOTO statement is not also
266       within.  Actually, 99.99% of REDUCE BEGIN-blocks are less
267       confusing if written entirely without GOTOs, and I mention them
268       primarily to explain WHILE-loops in terms of a more primitive
269       notion.;
270
271pause;
272
273COMMENT
274   7.  The LIMIT function provides a good illustration of nested
275       conditional expressions.  Proceeding sequentially through such
276       nests, each ELSE clause is matched with the nearest preceding
277       unmatched THEN clause in the group or block.  In order to help
278       reveal their structure, I have consistently indented nested
279       conditional statements, continuations of multi-line statements
280       and loop-bodies according to one of the many staunchly defended
281       indentation styles.  (If you have an instructor, I also urge
282       you to humor him by adopting his style for the duration of the
283       course.)
284   8.  C and Java programmers take note: "IF ... THEN ... ELSE ..." is
285       regarded as one expression, and semicolons are used to separate
286       rather than terminate statements.  Moreover, BEGIN and END are
287       brackets rather than statements, so a semicolon is never needed
288       immediately after BEGIN, and a semicolon is necessary
289       immediately preceding END only if the END is intended as a
290       labeled destination for a GOTO.  Within conditional
291       expressions, an inappropriate semicolon after an END, a >>, or
292       an ELSE-clause is likely to be one of your most prevalent
293       mistakes.;
294
295pause;
296
297COMMENT The next exercise is based on the above LIMIT function:
298
299For the sum of positive expressions AJ for J ranging from some finite
300initial value to infinity, the infinite series converges if the limit
301of the ratio SUB(J=J+1,AJ)/AJ is less than 1 as J approaches infinity.
302The series diverges if this limit exceeds 1, and the test is
303inconclusive if the limit is 1.  To convert the problem to the form
304required by the above LIMIT program, we can replace J by 1/!*FOO in
305the ratio, then take the limit as the indeterminate !*FOO approaches
306zero.  (Since an indeterminate is necessary here, I picked the weird
307name !*FOO to make the chance of conflict negligible.)
308
309After writing such a function to perform the ratio test, test it on
310the examples AJ=J/2^J, AJ=1/J^2, AJ=2^J/J^10, and AJ=1/J.  (The first
311two converge and the second two diverge.);
312
313pause;
314
315COMMENT Groups or blocks can be used wherever any arbitrary expression
316is allowed, including the right-hand side of a LET rule.
317
318The need for loops with an integer index variable running from a given
319initial value through a given final value by a given increment is so
320prevalent that REDUCE offers a convenient special way of accomplishing
321it via a FOR-loop, which has the general form
322
323   FOR index := initial STEP increment UNTIL final DO statement.
324
325Except for the use of commas as statement separators, this construct
326is equivalent to
327
328   BEGIN INTEGER index,
329   index := initial,
330   IF increment>0 THEN WHILE index <= final DO <<
331      statement,
332      index := index + increment >>
333   ELSE WHILE index >= final DO <<
334      statement,
335      index := index + increment >>,
336   RETURN 0
337   END;
338
339pause;
340
341COMMENT Note:
342
343   1.  The index variable is automatically declared local to the FOR-
344       loop.
345   2.  "initial", "increment", and "final" must have integer values.
346   3.  FORTRAN programmers take note: the body of the loop is not
347       automatically executed at least once.
348   4.  An abbreviation for "STEP 1 UNTIL" is ":".
349   5.  Since the WHILE-loop and the FOR-loop have implied BEGIN-
350       blocks, a RETURN statement within their bodies cannot transfer
351       control further than the point following the loops.
352
353Another frequent need is to produce output from within a group or
354block, because such output is not automatically produced.  This can be
355done using the WRITE-statement, which has the form
356
357   WRITE expression1, expression2, ..., expressionN.
358
359Beginning a new line with expression1, the expressions are printed
360immediately adjacent to each other, split over line boundaries if
361necessary.  The value of the WRITE-statement is the value of its last
362expression, and any of the expressions can be a character-string of
363the form "character1 character2 ... characterM".
364
365Inserting the word "WRITE" on a separate line before an assignment is
366convenient for debugging, because the word is then easily deleted
367afterward.  These features and others are illustrated by the following
368equation solver:;
369
370pause;
371
372operator solvefor, soln;
373
374for all x, lhs, rhs let solvefor(x, lhs, rhs) = solvefor(x, lhs-rhs);
375
376COMMENT LHS and RHS are expressions such that P=NUM(LHS-RHS) is a
377polynomial of degree at most 2 in the indeterminate or functional form
378X.  Otherwise an error message is printed.  As a convenience, RHS can
379be omitted if it is 0.  If P is quadratic in X, the two values of X
380which satisfy P=0 are stored as the values of the functional forms
381SOLN(1) and SOLN(2).  If P is a first-degree polynomial in X, SOLN(1)
382is set to the one solution.  If P simplifies to 0, SOLN(1) is set to
383the identifier ARBITRARY.  If P is an expression which does not
384simplify to zero but does not contain X, SOLN(1) is set to the
385identifier NONE.  In all other cases, SOLN(1) is set to the identifier
386UNKNOWN.  The function then returns the number of SOLN forms which
387were set.  This function prints a well deserved warning message if the
388denominator of LHS-RHS contains X.  If LHS-RHS is not polynomial in X,
389it is wise to execute this function under the influence of ON GCD.;
390
391pause;
392
393for all x, lhsmrhs let solvefor(x, lhsmrhs) =
394   begin integer hipow;  scalar temp, cflist, cf0, cf1, cf2;
395   if lhsmrhs = 0 then <<
396      soln(1) := arbitrary;
397      return 1 >>;
398   cflist :=  coeff(lhsmrhs, x);
399   hipow := hipow!*;
400   if hipow = 0 then <<
401      soln(1) := none;
402      return 1 >>;
403   if hipow > 2 then <<
404      soln(1) := unknown;
405      return 1 >>;
406   if hipow = 1 then <<
407      soln(1) := first(cflist)/second(cflist);
408      if df(sub(x=!*foo, soln(1)), !*foo) neq 0 then
409         soln(1) := unknown;
410      return 1 >>;
411   cf0 := first(cflist)/third(cflist);
412   cf1 := -second(cflist)/third(cflist)/2;
413   if df(sub(x=!*foo, cf0), !*foo) neq 0
414         or df(sub(x=!*foo, cf1), !*foo) neq 0  then <<
415      soln(1) := unknown;
416      return 1 >>;
417   temp := (cf1^2 - cf0)^(1/2);
418   soln(1) := cf1 + temp;
419   soln(2) := cf1 - temp;
420   return 2
421   end;
422
423COMMENT And some examples:;
424
425pause;
426
427for k:=1:solvefor(x, a*x^2, -b*x-c) do write soln(k) := soln(k);
428
429for k:=1:solvefor(log(x), 5*log(x)-7) do write soln(k) := soln(k);
430
431for k:=1:solvefor(x, x, x) do write soln(k) := soln(k);
432
433for k:= 1:solvefor(x, 5) do write soln(k) := soln(k);
434
435for k:=1:solvefor(x, x^3+x+1) do write soln(k) := soln(k);
436
437for k:=1:solvefor(x, x*e^x, 1) do write soln(k) := soln(k);
438
439g1 := x/(e^x-1);
440
441% Results in 'invalid as polynomial' error, continue anyway:;
442for k:=1:solvefor(x, g1) do write soln(k) := soln(k);
443
444sub(x=soln(1), g1);
445
446limit(g1, x, soln(1));
447
448pause;
449
450COMMENT Here we have used LET rules to permit the user the convenience
451of omitting default arguments.  (Function definitions have to have a
452fixed number of parameters.)
453
454Array elements are designated by the same syntax as matrix elements,
455namely as functional forms having integer arguments.  Here are some
456desiderata that may help you decide which of these alternatives is
457most appropriate for a particular application:
458
459   1.  The lower bound of each array subscript is 0 vs. 1 for
460       matrices vs. unrestricted for functional forms.
461   2.  The upper bound of each array subscript must have a specific
462       integer value at the time the array is declared, as must the
463       upper bounds of matrix subscripts when a matrix is first
464       referred to, on the left side of a matrix assignment.  In
465       contrast, functional forms never require a commitment to a
466       specific upper bound.
467   3.  An array can have any fixed number of subscripts, a matrix must
468       have exactly 2, and a functional form can have a varying
469       arbitrary number.
470   4.  Matrix operations, such as transpose and inverse, are built-in
471       only for matrices.
472   5.  For most implementations, access to array elements requires
473       time approximately proportional to the number of subscripts,
474       whereas access to matrix elements takes time approximately
475       proportional to the sum of the two subscript values, whereas
476       access to functional forms takes average time approximately
477       proportional to the number of bound functional forms having
478       that name.
479   6.  Only functional forms permit the effect of a subscripted
480       indeterminate such as having an answer be "A(M,N) + B(3,4)".
481   7.  Only functional forms can be used alone in the LHS of LET
482       substitutions.;
483
484pause;
485
486COMMENT
487   8.  All arrays, matrices, and operators are global regardless of
488       where they are declared, so declaring them within a BEGIN-block
489       does not afford the protection and automatic storage recovery
490       of local variables.  Moreover, clearing them within a
491       BEGIN-block will clear them globally, and normal functions
492       cannot return an array or a matrix value.  Furthermore, REDUCE
493       parameters are referenced by value, which means that an
494       assignment to a parameter has no effect on the corresponding
495       argument.  Thus, matrix or array results cannot be transmitted
496       back to an argument either.
497
498   9.  It is often advantageous to use two or more of these
499       alternatives to represent a set of quantities at different
500       times in the same program.  For example, to get the general
501       form of the inverse of a 3-by-3 matrix, we could write
502
503          MATRIX AA(3,3),
504          OPERATOR A,
505          FOR J:=1:3 DO
506             FOR K:=1:3 DO AA(J,K) := A(J,K),
507          AA^-1.
508
509       As another example, we might use an array to receive some
510       polynomial coefficients, then transfer the values to a matrix
511       for inversion.;
512
513pause;
514
515COMMENT The COEFF function is the remaining new feature in our
516SOLVEFOR example.  The first argument is a polynomial expression in
517the indeterminate or functional form which is the second argument.
518The polynomial coefficients of the integer powers of the indeterminate
519are returned as a LIST, with the independent coefficient first.  The
520highest and lowest non-zero powers are placed in the variables HIPOW!*
521and LOWPOW!* respectively.
522
523A LIST is a kind of data structure, just as matrices and arrays are.
524It is represented as a comma-separated sequence of elements enclosed
525in braces.  The elements can be accessed with the functions FIRST,
526SECOND, THIRD, PART(i) which returns the i-th element, and REST, which
527returns a list of all but the first element.  For example:;
528
529clear x;
530
531coeff(x^5+2, x);
532
533lowpow!*;
534
535hipow!*;
536
537pause;
538
539COMMENT COEFF does not check to make sure that the coefficients do not
540contain its second argument within a functional form, so that is the
541reason we differentiated.  The reason we first substituted the
542indeterminate !*FOO for the second argument is that differentiation
543does not work with respect to a functional form.
544
545The last exercise is to rewrite the last rule so that we can solve
546equations which simplify to the form
547
548   a*x^(m+2*l) + b*x^(m+l) + c*x^m = 0,  where m >= 0 and l >= 1.
549
550The solutions are
551
552   0,  with multiplicity m,
553   x1*E^(2*j*I*pi/l),
554   x2*E^(2*j*I*pi/l),  with j = 0, 1, ..., l-1,
555
556where x1 and x2 are the solutions to the quadratic equation
557
558   a*x^2 + b*x + c = 0.
559
560As a convenience to the user, you might also wish to have a global
561switch named SOLVEPRINT, such that when it is nonzero, the solutions
562are automatically printed.
563
564This is the end of lesson 4.  When you are ready to run lesson 5,
565start a new REDUCE session.
566
567;end;
568