This chapter explains how to use Calc and its many features, in a step-by-step, tutorial way. You are encouraged to run Calc and work along with the examples as you read (see section Starting Calc). If you are already familiar with advanced calculators, you may wish to skip on to the rest of this manual.
This tutorial describes the standard user interface of Calc only. The "Quick Mode" and "Keypad Mode" interfaces are fairly self-explanatory. See section Embedded Mode, for a description of the "Embedded Mode" interface.
The easiest way to read this tutorial on-line is to have two windows on your Emacs screen, one with Calc and one with the Info system. (If you have a printed copy of the manual you can use that instead.) Press M-# c to turn Calc on or to switch into the Calc window, and press M-# i to start the Info system or to switch into its window.
This tutorial is designed to be done in sequence. But the rest of this manual does not assume you have gone through the tutorial. The tutorial does not cover everything in the Calculator, but it touches on most general areas.
The Calc Summary at the end of the reference manual includes some blank space for your own use. You may wish to keep notes there as you learn Calc.
In this section, we learn how RPN and algebraic-style calculations work, how to undo and redo an operation done by mistake, and how to control various modes of the Calculator.
The central component of an RPN calculator is the stack. A calculator stack is like a stack of dishes. New dishes (numbers) are added at the top of the stack, and numbers are normally only removed from the top of the stack.
In an operation like 2+3, the 2 and 3 are called the operands and the + is the operator. In an RPN calculator you always enter the operands first, then the operator. Each time you type a number, Calc adds or pushes it onto the top of the Stack. When you press an operator key like +, Calc pops the appropriate number of operands from the stack and pushes back the result.
Thus we could add the numbers 2 and 3 in an RPN calculator by typing: 2 RET 3 RET +. (The RET key, Return, corresponds to the ENTER key on traditional RPN calculators.) Try this now if you wish; type M-# c to switch into the Calc window (you can type M-# c again or M-# o to switch back to the Tutorial window). The first four keystrokes "push" the numbers 2 and 3 onto the stack. The + key "pops" the top two numbers from the stack, adds them, and pushes the result (5) back onto the stack. Here's how the stack will look at various points throughout the calculation:
. 1: 2 2: 2 1: 5 .
. 1: 3 .
.
M-# c 2 RET 3 RET + DEL
The `.' symbol is a marker that represents the top of the stack. Note that the "top" of the stack is really shown at the bottom of the Stack window. This may seem backwards, but it turns out to be less distracting in regular use.
The numbers `1:' and `2:' on the left are stack level numbers. Old RPN calculators always had four stack levels called x, y, z, and t. Calc's stack can grow as large as you like, so it uses numbers instead of letters. Some stack-manipulation commands accept a numeric argument that says which stack level to work on. Normal commands like + always work on the top few levels of the stack.
The Stack buffer is just an Emacs buffer, and you can move around in it using the regular Emacs motion commands. But no matter where the cursor is, even if you have scrolled the `.' marker out of view, most Calc commands always move the cursor back down to level 1 before doing anything. It is possible to move the `.' marker upwards through the stack, temporarily "hiding" some numbers from commands like +. This is called stack truncation and we will not cover it in this tutorial; see section Truncating the Stack, if you are interested.
You don't really need the second RET in 2 RET 3 RET +. That's because if you type any operator name or other non-numeric key when you are entering a number, the Calculator automatically enters that number and then does the requested command. Thus 2 RET 3 + will work just as well.
Examples in this tutorial will often omit RET even when the stack displays shown would only happen if you did press RET:
1: 2 2: 2 1: 5
. 1: 3 .
.
2 RET 3 +
Here, after pressing 3 the stack would really show `1: 2' with `Calc: 3' in the minibuffer. In these situations, you can press the optional RET to see the stack as the figure shows.
(*) Exercise 1. (This tutorial will include exercises at various points. Try them if you wish. Answers to all the exercises are located at the end of the Tutorial chapter. Each exercise will include a cross-reference to its particular answer. If you are reading with the Emacs Info system, press f and the exercise number to go to the answer, then the letter l to return to where you were.)
Here's the first exercise: What will the keystrokes 1 RET 2 RET 3 RET 4 + * - compute? (`*' is the symbol for multiplication.) Figure it out by hand, then try it with Calc to see if you're right. See section RPN Tutorial Exercise 1. (*)
(*) Exercise 2. Compute @c{$(2\times4) + (7\times9.4) + {5\over4}$} 2*4 + 7*9.5 + 5/4 using the stack. See section RPN Tutorial Exercise 2. (*)
The DEL key is called Backspace on some keyboards. It is whatever key you would use to correct a simple typing error when regularly using Emacs. The DEL key pops and throws away the top value on the stack. (You can still get that value back from the Trail if you should need it later on.) There are many places in this tutorial where we assume you have used DEL to erase the results of the previous example at the beginning of a new example. In the few places where it is really important to use DEL to clear away old results, the text will remind you to do so.
(It won't hurt to let things accumulate on the stack, except that whenever you give a display-mode-changing command Calc will have to spend a long time reformatting such a large stack.)
Since the - key is also an operator (it subtracts the top two stack elements), how does one enter a negative number? Calc uses the _ (underscore) key to act like the minus sign in a number. So, typing -5 RET won't work because the - key will try to do a subtraction, but _5 RET works just fine.
You can also press n, which means "change sign." It changes the number at the top of the stack (or the number being entered) from positive to negative or vice-versa: 5 n RET.
If you press RET when you're not entering a number, the effect is to duplicate the top number on the stack. Consider this calculation:
1: 3 2: 3 1: 9 2: 9 1: 81
. 1: 3 . 1: 9 .
. .
3 RET RET * RET *
(Of course, an easier way to do this would be 3 RET 4 ^, to raise 3 to the fourth power.)
The space-bar key (denoted SPC here) performs the same function as RET; you could replace all three occurrences of RET in the above example with SPC and the effect would be the same.
Another stack manipulation key is TAB. This exchanges the top two stack entries. Suppose you have computed 2 RET 3 + to get 5, and then you realize what you really wanted to compute was 20 / (2+3).
1: 5 2: 5 2: 20 1: 4
. 1: 20 1: 5 .
. .
2 RET 3 + 20 TAB /
Planning ahead, the calculation would have gone like this:
1: 20 2: 20 3: 20 2: 20 1: 4
. 1: 2 2: 2 1: 5 .
. 1: 3 .
.
20 RET 2 RET 3 + /
A related stack command is M-TAB (hold META and type TAB). It rotates the top three elements of the stack upward, bringing the object in level 3 to the top.
1: 10 2: 10 3: 10 3: 20 3: 30
. 1: 20 2: 20 2: 30 2: 10
. 1: 30 1: 10 1: 20
. . .
10 RET 20 RET 30 RET M-TAB M-TAB
(*) Exercise 3. Suppose the numbers 10, 20, and 30 are on the stack. Figure out how to add one to the number in level 2 without affecting the rest of the stack. Also figure out how to add one to the number in level 3. See section RPN Tutorial Exercise 3. (*)
Operations like +, -, *, /, and ^ pop two arguments from the stack and push a result. Operations like n and Q (square root) pop a single number and push the result. You can think of them as simply operating on the top element of the stack.
1: 3 1: 9 2: 9 1: 25 1: 5
. . 1: 16 . .
.
3 RET RET * 4 RET RET * + Q
(Note that capital Q means to hold down the Shift key while typing q. Remember, plain unshifted q is the Quit command.)
Here we've used the Pythagorean Theorem to determine the hypotenuse of a right triangle. Calc actually has a built-in command for that called f h, but let's suppose we can't remember the necessary keystrokes. We can still enter it by its full name using M-x notation:
1: 3 2: 3 1: 5
. 1: 4 .
.
3 RET 4 RET M-x calc-hypot
All Calculator commands begin with the word `calc-'. Since it gets tiring to type this, Calc provides an x key which is just like the regular Emacs M-x key except that it types the `calc-' prefix for you:
1: 3 2: 3 1: 5
. 1: 4 .
.
3 RET 4 RET x hypot
What happens if you take the square root of a negative number?
1: 4 1: -4 1: (0, 2)
. . .
4 RET n Q
The notation (a, b) represents a complex number. Complex numbers are more traditionally written @c{$a + b i$} a + b i; Calc can display in this format, too, but for now we'll stick to the (a, b) notation.
If you don't know how complex numbers work, you can safely ignore this feature. Complex numbers only arise from operations that would be errors in a calculator that didn't have complex numbers. (For example, taking the square root or logarithm of a negative number produces a complex result.)
Complex numbers are entered in the notation shown. The ( and , and ) keys manipulate "incomplete complex numbers."
1: ( ... 2: ( ... 1: (2, ... 1: (2, ... 1: (2, 3)
. 1: 2 . 3 .
. .
( 2 , 3 )
You can perform calculations while entering parts of incomplete objects. However, an incomplete object cannot actually participate in a calculation:
1: ( ... 2: ( ... 3: ( ... 1: ( ... 1: ( ...
. 1: 2 2: 2 5 5
. 1: 3 . .
.
(error)
( 2 RET 3 + +
Adding 5 to an incomplete object makes no sense, so the last command produces an error message and leaves the stack the same.
Incomplete objects can't participate in arithmetic, but they can be moved around by the regular stack commands.
2: 2 3: 2 3: 3 1: ( ... 1: (2, 3)
1: 3 2: 3 2: ( ... 2 .
. 1: ( ... 1: 2 3
. . .
2 RET 3 RET ( M-TAB M-TAB )
Note that the , (comma) key did not have to be used here. When you press ) all the stack entries between the incomplete entry and the top are collected, so there's never really a reason to use the comma. It's up to you.
(*) Exercise 4. To enter the complex number (2, 3), your friend Joe typed ( 2 , SPC 3 ). What happened? (Joe thought of a clever way to correct his mistake in only two keystrokes, but it didn't quite work. Try it to find out why.) See section RPN Tutorial Exercise 4. (*)
Vectors are entered the same way as complex numbers, but with square brackets in place of parentheses. We'll meet vectors again later in the tutorial.
Any Emacs command can be given a numeric prefix argument by typing a series of META-digits beforehand. If META is awkward for you, you can instead type C-u followed by the necessary digits. Numeric prefix arguments can be negative, as in M-- M-3 M-5 or C-u - 3 5. Calc commands use numeric prefix arguments in a variety of ways. For example, a numeric prefix on the + operator adds any number of stack entries at once:
1: 10 2: 10 3: 10 3: 10 1: 60
. 1: 20 2: 20 2: 20 .
. 1: 30 1: 30
. .
10 RET 20 RET 30 RET C-u 3 +
For stack manipulation commands like RET, a positive numeric prefix argument operates on the top n stack entries at once. A negative argument operates on the entry in level n only. An argument of zero operates on the entire stack. In this example, we copy the second-to-top element of the stack:
1: 10 2: 10 3: 10 3: 10 4: 10
. 1: 20 2: 20 2: 20 3: 20
. 1: 30 1: 30 2: 30
. . 1: 20
.
10 RET 20 RET 30 RET C-u -2 RET
Another common idiom is M-0 DEL, which clears the stack. (The M-0 numeric prefix tells DEL to operate on the entire stack.)
If you are not used to RPN notation, you may prefer to operate the Calculator in "algebraic mode," which is closer to the way non-RPN calculators work. In algebraic mode, you enter formulas in traditional 2+3 notation.
You don't really need any special "mode" to enter algebraic formulas. You can enter a formula at any time by pressing the apostrophe (') key. Answer the prompt with the desired formula, then press RET. The formula is evaluated and the result is pushed onto the RPN stack. If you don't want to think in RPN at all, you can enter your whole computation as a formula, read the result from the stack, then press DEL to delete it from the stack.
Try pressing the apostrophe key, then 2+3+4, then RET. The result should be the number 9.
Algebraic formulas use the operators `+', `-', `*', `/', and `^'. You can use parentheses to make the order of evaluation clear. In the absence of parentheses, `^' is evaluated first, then `*', then `/', then finally `+' and `-'. For example, the expression
2 + 3*4*5 / 6*7^8 - 9
is equivalent to
2 + ((3*4*5) / (6*(7^8)) - 9
or, in large mathematical notation,
The result of this expression will be the number -6.99999826533.
Calc's order of evaluation is the same as for most computer languages, except that `*' binds more strongly than `/', as the above example shows. As in normal mathematical notation, the `*' symbol can often be omitted: `2 a' is the same as `2*a'.
Operators at the same level are evaluated from left to right, except that `^' is evaluated from right to left. Thus, `2-3-4' is equivalent to `(2-3)-4' or -5, whereas `2^3^4' is equivalent to `2^(3^4)' (a very large integer; try it!).
If you tire of typing the apostrophe all the time, there is an "algebraic mode" you can select in which Calc automatically senses when you are about to type an algebraic expression. To enter this mode, press the two letters m a. (An `Alg' indicator should appear in the Calc window's mode line.)
Press m a, then 2+3+4 with no apostrophe, then RET.
In algebraic mode, when you press any key that would normally begin entering a number (such as a digit, a decimal point, or the _ key), or if you press ( or [, Calc automatically begins an algebraic entry.
Functions which do not have operator symbols like `+' and `*'
must be entered in formulas using function-call notation. For example,
the function name corresponding to the square-root key Q is
sqrt. To compute a square root in a formula, you would use
the notation `sqrt(x)'.
Press the apostrophe, then type sqrt(5*2) - 3. The result should be 0.16227766017.
Note that if the formula begins with a function name, you need to use the apostrophe even if you are in algebraic mode. If you type arcsin out of the blue, the a r will be taken as an Algebraic Rewrite command, and the csin will be taken as the name of the rewrite rule to use!
Some people prefer to enter complex numbers and vectors in algebraic form because they find RPN entry with incomplete objects to be too distracting, even though they otherwise use Calc as an RPN calculator.
Still in algebraic mode, type:
1: (2, 3) 2: (2, 3) 1: (8, -1) 2: (8, -1) 1: (9, -1)
. 1: (1, -2) . 1: 1 .
. .
(2,3) RET (1,-2) RET * 1 RET +
Algebraic mode allows us to enter complex numbers without pressing an apostrophe first, but it also means we need to press RET after every entry, even for a simple number like 1.
(You can type C-u m a to enable a special "incomplete algebraic mode" in which the ( and [ keys use algebraic entry even though regular numeric keys still use RPN numeric entry. There is also a "total algebraic mode," started by typing m t, in which all normal keys begin algebraic entry. You must then use the META key to type Calc commands: M-m t to get back out of total algebraic mode, M-q to quit, etc. Total algebraic mode is not supported under Emacs 19.)
If you're still in algebraic mode, press m a again to turn it off.
Actual non-RPN calculators use a mixture of algebraic and RPN styles. In general, operators of two numbers (like + and *) use algebraic form, but operators of one number (like n and Q) use RPN form. Also, a non-RPN calculator allows you to see the intermediate results of a calculation as you go along. You can accomplish this in Calc by performing your calculation as a series of algebraic entries, using the $ sign to tie them together. In an algebraic formula, $ represents the number on the top of the stack. Here, we perform the calculation @c{$\sqrt{2\times4+1}$} sqrt(2*4+1), which on a traditional calculator would be done by pressing 2 * 4 + 1 = and then the square-root key.
1: 8 1: 9 1: 3
. . .
' 2*4 RET $+1 RET Q
Notice that we didn't need to press an apostrophe for the $+1, because the dollar sign always begins an algebraic entry.
(*) Exercise 1. How could you get the same effect as pressing Q but using an algebraic entry instead? How about if the Q key on your keyboard were broken? See section Algebraic Entry Tutorial Exercise 1. (*)
The notations $$, $$$, and so on stand for higher stack entries. For example, ' $$+$ RET is just like typing +.
Algebraic formulas can include variables. To store in a variable, press s s, then type the variable name, then press RET. (There are actually two flavors of store command: s s stores a number in a variable but also leaves the number on the stack, while s t removes a number from the stack and stores it in the variable.) A variable name should consist of one or more letters or digits, beginning with a letter.
1: 17 . 1: a + a^2 1: 306
. . .
17 s t a RET ' a+a^2 RET =
The = key evaluates a formula by replacing all its variables by the values that were stored in them.
For RPN calculations, you can recall a variable's value on the stack either by entering its name as a formula and pressing =, or by using the s r command.
1: 17 2: 17 3: 17 2: 17 1: 306
. 1: 17 2: 17 1: 289 .
. 1: 2 .
.
s r a RET ' a RET = 2 ^ +
If you press a single digit for a variable name (as in s t 3, you
get one of ten quick variables q0 through q9.
They are "quick" simply because you don't have to type the letter
q or the RET after their names. In fact, you can type
simply s 3 as a shorthand for s s 3, and likewise for
t 3 and r 3.
Any variables in an algebraic formula for which you have not stored values are left alone, even when you evaluate the formula.
1: 2 a + 2 b 1: 34 + 2 b
. .
' 2a+2b RET =
Calls to function names which are undefined in Calc are also left alone, as are calls for which the value is undefined.
1: 2 + log10(0) + log10(x) + log10(5, 6) + foo(3)
.
' log10(100) + log10(0) + log10(x) + log10(5,6) + foo(3) RET
In this example, the first call to log10 works, but the other
calls are not evaluated. In the second call, the logarithm is
undefined for that value of the argument; in the third, the argument
is symbolic, and in the fourth, there are too many arguments. In the
fifth case, there is no function called foo. You will see a
"Wrong number of arguments" message referring to `log10(5,6)'.
Press the w ("why") key to see any other messages that may
have arisen from the last calculation. In this case you will get
"logarithm of zero," then "number expected: x". Calc
automatically displays the first message only if the message is
sufficiently important; for example, Calc considers "wrong number
of arguments" and "logarithm of zero" to be important enough to
report automatically, while a message like "number expected: x"
will only show up if you explicitly press the w key.
(*) Exercise 2. Joe entered the formula `2 x y',
stored 5 in x, pressed =, and got the expected result,
`10 y'. He then tried the same for the formula `2 x (1+y)',
expecting `10 (1+y)', but it didn't work. Why not?
See section Algebraic Entry Tutorial Exercise 2. (*)
(*) Exercise 3. What result would you expect 1 RET 0 / to give? What if you then type 0 *? See section Algebraic Entry Tutorial Exercise 3. (*)
One interesting way to work with variables is to use the evaluates-to (`=>') operator. It works like this: Enter a formula algebraically in the usual way, but follow the formula with an `=>' symbol. (There is also an s = command which builds an `=>' formula using the stack.) On the stack, you will see two copies of the formula with an `=>' between them. The lefthand formula is exactly like you typed it; the righthand formula has been evaluated as if by typing =.
2: 2 + 3 => 5 2: 2 + 3 => 5
1: 2 a + 2 b => 34 + 2 b 1: 2 a + 2 b => 20 + 2 b
. .
' 2+3 => RET ' 2a+2b RET s = 10 s t a RET
Notice that the instant we stored a new value in a, all
`=>' operators already on the stack that referred to a
were updated to use the new value. With `=>', you can push a
set of formulas on the stack, then change the variables experimentally
to see the effects on the formulas' values.
You can also "unstore" a variable when you are through with it:
2: 2 + 5 => 5
1: 2 a + 2 b => 2 a + 2 b
.
s u a RET
We will encounter formulas involving variables and functions again when we discuss the algebra and calculus features of the Calculator.
If you make a mistake, you can usually correct it by pressing shift-U, the "undo" command. First, clear the stack (M-0 DEL) and exit and restart Calc (M-# M-# M-# M-#) to make sure things start off with a clean slate. Now:
1: 2 2: 2 1: 8 2: 2 1: 6
. 1: 3 . 1: 3 .
. .
2 RET 3 ^ U *
You can undo any number of times. Calc keeps a complete record of all you have done since you last opened the Calc window. After the above example, you could type:
1: 6 2: 2 1: 2 . .
. 1: 3 .
.
(error)
U U U U
You can also type D to "redo" a command that you have undone mistakenly.
. 1: 2 2: 2 1: 6 1: 6
. 1: 3 . .
.
(error)
D D D D
It was not possible to redo past the 6, since that was placed there by something other than an undo command.
You can think of undo and redo as a sort of "time machine." Press U to go backward in time, D to go forward. If you go backward and do something (like *) then, as any science fiction reader knows, you have changed your future and you cannot go forward again. Thus, the inability to redo past the 6 even though there was an earlier undo command.
You can always recall an earlier result using the Trail. We've ignored the trail so far, but it has been faithfully recording everything we did since we loaded the Calculator. If the Trail is not displayed, press t d now to turn it on.
Let's try grabbing an earlier result. The 8 we computed was undone by a U command, and was lost even to Redo when we pressed *, but it's still there in the trail. There should be a little `>' arrow (the trail pointer) resting on the last trail entry. If there isn't, press t ] to reset the trail pointer. Now, press t p to move the arrow onto the line containing 8, and press t y to "yank" that number back onto the stack.
If you press t ] again, you will see that even our Yank command went into the trail.
Let's go further back in time. Earlier in the tutorial we computed a huge integer using the formula `2^3^4'. We don't remember what it was, but the first digits were "241". Press t r (which stands for trail-search-reverse), then type 241. The trail cursor will jump back to the next previous occurrence of the string "241" in the trail. This is just a regular Emacs incremental search; you can now press C-s or C-r to continue the search forwards or backwards as you like.
To finish the search, press RET. This halts the incremental search and leaves the trail pointer at the thing we found. Now we can type t y to yank that number onto the stack. If we hadn't remembered the "241", we could simply have searched for 2^3^4, then pressed RET t n to halt and then move to the next item.
You may have noticed that all the trail-related commands begin with the letter t. (The store-and-recall commands, on the other hand, all began with s.) Calc has so many commands that there aren't enough keys for all of them, so various commands are grouped into two-letter sequences where the first letter is called the prefix key. If you type a prefix key by accident, you can press C-g to cancel it. (In fact, you can press C-g to cancel almost anything in Emacs.) To get help on a prefix key, press that key followed by ?. Some prefixes have several lines of help, so you need to press ? repeatedly to see them all. This may not work under Lucid Emacs, but you can also type h h to see all the help at once.
Try pressing t ? now. You will see a line of the form,
trail/time: Display; Fwd, Back; Next, Prev, Here, [, ]; Yank: [MORE] t-
The word "trail" indicates that the t prefix key contains trail-related commands. Each entry on the line shows one command, with a single capital letter showing which letter you press to get that command. We have used t n, t p, t ], and t y so far. The `[MORE]' means you can press ? again to see more t-prefix comands. Notice that the commands are roughly divided (by semicolons) into related groups.
When you are in the help display for a prefix key, the prefix is still active. If you press another key, like y for example, it will be interpreted as a t y command. If all you wanted was to look at the help messages, press C-g afterwards to cancel the prefix.
One more way to correct an error is by editing the stack entries. The actual Stack buffer is marked read-only and must not be edited directly, but you can press ` (the backquote or accent grave) to edit a stack entry.
Try entering `3.141439' now. If this is supposed to represent pi, it's got several errors. Press ` to edit this number. Now use the normal Emacs cursor motion and editing keys to change the second 4 to a 5, and to transpose the 3 and the 9. When you press RET, the number on the stack will be replaced by your new number. This works for formulas, vectors, and all other types of values you can put on the stack. The ` key also works during entry of a number or algebraic formula.
Calc has many types of modes that affect the way it interprets your commands or the way it displays data. We have already seen one mode, namely algebraic mode. There are many others, too; we'll try some of the most common ones here.
Perhaps the most fundamental mode in Calc is the current precision. Notice the `12' on the Calc window's mode line:
--%%-Calc: 12 Deg (Calculator)----All------
Most of the symbols there are Emacs things you don't need to worry about, but the `12' and the `Deg' are mode indicators. The `12' means that calculations should always be carried to 12 significant figures. That is why, when we type 1 RET 7 /, we get 0.142857142857 with exactly 12 digits, not counting leading and trailing zeros.
You can set the precision to anything you like by pressing p, then entering a suitable number. Try pressing p 30 RET, then doing 1 RET 7 / again:
1: 0.142857142857
2: 0.142857142857142857142857142857
.
Although the precision can be set arbitrarily high, Calc always has to have some value for the current precision. After all, the true value 1/7 is an infinitely repeating decimal; Calc has to stop somewhere.
Of course, calculations are slower the more digits you request. Press p 12 now to set the precision back down to the default.
Calculations always use the current precision. For example, even though we have a 30-digit value for 1/7 on the stack, if we use it in a calculation in 12-digit mode it will be rounded down to 12 digits before it is used. Try it; press RET to duplicate the number, then 1 +. Notice that the RET key didn't round the number, because it doesn't do any calculation. But the instant we pressed +, the number was rounded down.
1: 0.142857142857
2: 0.142857142857142857142857142857
3: 1.14285714286
.
In fact, since we added a digit on the left, we had to lose one digit on the right from even the 12-digit value of 1/7.
How did we get more than 12 digits when we computed `2^3^4'? The answer is that Calc makes a distinction between integers and floating-point numbers, or floats. An integer is a number that does not contain a decimal point. There is no such thing as an "infinitely repeating fraction integer," so Calc doesn't have to limit itself. If you asked for `2^10000' (don't try this!), you would have to wait a long time but you would eventually get an exact answer. If you ask for `2.^10000', you will quickly get an answer which is correct only to 12 places. The decimal point tells Calc that it should use floating-point arithmetic to get the answer, not exact integer arithmetic.
You can use the F (calc-floor) command to convert a
floating-point value to an integer, and c f (calc-float)
to convert an integer to floating-point form.
Let's try entering that last calculation:
1: 2. 2: 2. 1: 1.99506311689e3010
. 1: 10000 .
.
2.0 RET 10000 RET ^
Notice the letter `e' in there. It represents "times ten to the power of," and is used by Calc automatically whenever writing the number out fully would introduce more extra zeros than you probably want to see. You can enter numbers in this notation, too.
1: 2. 2: 2. 1: 1.99506311678e3010
. 1: 10000. .
.
2.0 RET 1e4 RET ^
Hey, the answer is different! Look closely at the middle columns of the two examples. In the first, the stack contained the exact integer 10000, but in the second it contained a floating-point value with a decimal point. When you raise a number to an integer power, Calc uses repeated squaring and multiplication to get the answer. When you use a floating-point power, Calc uses logarithms and exponentials. As you can see, a slight error crept in during one of these methods. Which one should we trust? Let's raise the precision a bit and find out:
. 1: 2. 2: 2. 1: 1.995063116880828e3010
. 1: 10000. .
.
p 16 RET 2. RET 1e4 ^ p 12 RET
Presumably, it doesn't matter whether we do this higher-precision calculation using an integer or floating-point power, since we have added enough "guard digits" to trust the first 12 digits no matter what. And the verdict is... Integer powers were more accurate; in fact, the result was only off by one unit in the last place.
Calc does many of its internal calculations to a slightly higher precision, but it doesn't always bump the precision up enough. In each case, Calc added about two digits of precision during its calculation and then rounded back down to 12 digits afterward. In one case, it was enough; in the the other, it wasn't. If you really need x digits of precision, it never hurts to do the calculation with a few extra guard digits.
What if we want guard digits but don't want to look at them? We can set the float format. Calc supports four major formats for floating-point numbers, called normal, fixed-point, scientific notation, and engineering notation. You get them by pressing d n, d f, d s, and d e, respectively. In each case, you can supply a numeric prefix argument which says how many digits should be displayed. As an example, let's put a few numbers onto the stack and try some different display modes. First, use M-0 DEL to clear the stack, then enter the four numbers shown here:
4: 12345 4: 12345 4: 12345 4: 12345 4: 12345
3: 12345. 3: 12300. 3: 1.2345e4 3: 1.23e4 3: 12345.000
2: 123.45 2: 123. 2: 1.2345e2 2: 1.23e2 2: 123.450
1: 12.345 1: 12.3 1: 1.2345e1 1: 1.23e1 1: 12.345
. . . . .
d n M-3 d n d s M-3 d s M-3 d f
Notice that when we typed M-3 d n, the numbers were rounded down to three significant digits, but then when we typed d s all five significant figures reappeared. The float format does not affect how numbers are stored, it only affects how they are displayed. Only the current precision governs the actual rounding of numbers in the Calculator's memory.
Engineering notation, not shown here, is like scientific notation except the exponent (the power-of-ten part) is always adjusted to be a multiple of three (as in "kilo," "micro," etc.). As a result there will be one, two, or three digits before the decimal point.
Whenever you change a display-related mode, Calc redraws everything in the stack. This may be slow if there are many things on the stack, so Calc allows you to type shift-H before any mode command to prevent it from updating the stack. Anything Calc displays after the mode-changing command will appear in the new format.
4: 12345 4: 12345 4: 12345 4: 12345 4: 12345
3: 12345.000 3: 12345.000 3: 12345.000 3: 1.2345e4 3: 12345.
2: 123.450 2: 123.450 2: 1.2345e1 2: 1.2345e1 2: 123.45
1: 12.345 1: 1.2345e1 1: 1.2345e2 1: 1.2345e2 1: 12.345
. . . . .
H d s DEL U TAB d SPC d n
Here the H d s command changes to scientific notation but without updating the screen. Deleting the top stack entry and undoing it back causes it to show up in the new format; swapping the top two stack entries reformats both entries. The d SPC command refreshes the whole stack. The d n command changes back to the normal float format; since it doesn't have an H prefix, it also updates all the stack entries to be in d n format.
Notice that the integer 12345 was not affected by any of the float formats. Integers are integers, and are always displayed exactly.
Large integers have their own problems. Let's look back at the result of 2^3^4.
2417851639229258349412352
Quick--how many digits does this have? Try typing d g:
2,417,851,639,229,258,349,412,352
Now how many digits does this have? It's much easier to tell! We can actually group digits into clumps of any size. Some people prefer M-5 d g:
24178,51639,22925,83494,12352
Let's see what happens to floating-point numbers when they are grouped. First, type p 25 RET to make sure we have enough precision to get ourselves into trouble. Now, type 1e13 /:
24,17851,63922.9258349412352
The integer part is grouped but the fractional part isn't. Now try M-- M-5 d g (that's meta-minus-sign, meta-five):
24,17851,63922.92583,49412,352
If you find it hard to tell the decimal point from the commas, try changing the grouping character to a space with d , SPC:
24 17851 63922.92583 49412 352
Type d , , to restore the normal grouping character, then d g again to turn grouping off. Also, press p 12 to restore the default precision.
Press U enough times to get the original big integer back. (Notice that U does not undo each mode-setting command; if you want to undo a mode-setting command, you have to do it yourself.) Now, type d r 16 RET:
16#200000000000000000000
The number is now displayed in hexadecimal, or "base-16" form. Suddenly it looks pretty simple; this should be no surprise, since we got this number by computing a power of two, and 16 is a power of 2. In fact, we can use d r 2 RET to see it in actual binary form:
2#1000000000000000000000000000000000000000000000000000000 ...
We don't have enough space here to show all the zeros! They won't fit on a typical screen, either, so you will have to use horizontal scrolling to see them all. Press < and > to scroll the stack window left and right by half its width. Another way to view something large is to press ` (back-quote) to edit the top of stack in a separate window. (Press M-# M-# when you are done.)
You can enter non-decimal numbers using the # symbol, too. Let's see what the hexadecimal number `5FE' looks like in binary. Type 16#5FE (the letters can be typed in upper or lower case; they will always appear in upper case). It will also help to turn grouping on with d g:
2#101,1111,1110
Notice that d g groups by fours by default if the display radix is binary or hexadecimal, but by threes if it is decimal, octal, or any other radix.
Now let's see that number in decimal; type d r 10:
1,534
Numbers are not stored with any particular radix attached. They're just numbers; they can be entered in any radix, and are always displayed in whatever radix you've chosen with d r. The current radix applies to integers, fractions, and floats.
(*) Exercise 1. Your friend Joe tried to enter one-third as `3#0.1' in d r 3 mode with a precision of 12. He got `3#0.0222222...' (with 25 2's) in the display. When he multiplied that by three, he got `3#0.222222...' instead of the expected `3#1'. Next, Joe entered `3#0.2' and, to his great relief, saw `3#0.2' on the screen. But when he typed 2 /, he got `3#0.10000001' (some zeros omitted). What's going on here? See section Modes Tutorial Exercise 1. (*)
(*) Exercise 2. Scientific notation works in non-decimal modes in the natural way (the exponent is a power of the radix instead of a power of ten, although the exponent itself is always written in decimal). Thus `8#1.23e3 = 8#1230.0'. Suppose we have the hexadecimal number `f.e8f' times 16 to the 15th power: We write `16#f.e8fe15'. What is wrong with this picture? What could we write instead that would work better? See section Modes Tutorial Exercise 2. (*)
The m prefix key has another set of modes, relating to the way Calc interprets your inputs and does computations. Whereas d-prefix modes generally affect the way things look, m-prefix modes affect the way they are actually computed.
The most popular m-prefix mode is the angular mode. Notice the `Deg' indicator in the mode line. This means that if you use a command that interprets a number as an angle, it will assume the angle is measured in degrees. For example,
1: 45 1: 0.707106781187 1: 0.500000000001 1: 0.5
. . . .
45 S 2 ^ c 1
The shift-S command computes the sine of an angle. The sine of 45 degrees is @c{$\sqrt{2}/2$} sqrt(2)/2; squaring this yields 2/4 = 0.5. However, there has been a slight roundoff error because the representation of @c{$\sqrt{2}/2$} sqrt(2)/2 wasn't exact. The c 1 command is a handy way to clean up numbers in this case; it temporarily reduces the precision by one digit while it re-rounds the number on the top of the stack.
(*) Exercise 3. Your friend Joe computed the sine of 45 degrees as shown above, then, hoping to avoid an inexact result, he increased the precision to 16 digits before squaring. What happened? See section Modes Tutorial Exercise 3. (*)
To do this calculation in radians, we would type m r first. (The indicator changes to `Rad'.) 45 degrees corresponds to pi/4 radians. To get @c{$\pi$} pi, press the P key. (Once again, this is a shifted capital P. Remember, unshifted p sets the precision.)
1: 3.14159265359 1: 0.785398163398 1: 0.707106781187
. . .
P 4 / m r S
Likewise, inverse trigonometric functions generate results in either radians or degrees, depending on the current angular mode.
1: 0.707106781187 1: 0.785398163398 1: 45.
. . .
.5 Q m r I S m d U I S
Here we compute the Inverse Sine of @c{$\sqrt{0.5}$} sqrt(0.5), first in radians, then in degrees.
Use c d and c r to convert a number from radians to degrees and vice-versa.
1: 45 1: 0.785398163397 1: 45.
. . .
45 c r c d
Another interesting mode is fraction mode. Normally, dividing two integers produces a floating-point result if the quotient can't be expressed as an exact integer. Fraction mode causes integer division to produce a fraction, i.e., a rational number, instead.
2: 12 1: 1.33333333333 1: 4:3
1: 9 . .
.
12 RET 9 / m f U / m f
In the first case, we get an approximate floating-point result. In the second case, we get an exact fractional result (four-thirds).
You can enter a fraction at any time using : notation. (Calc uses : instead of / as the fraction separator because / is already used to divide the top two stack elements.) Calculations involving fractions will always produce exact fractional results; fraction mode only says what to do when dividing two integers.
(*) Exercise 4. If fractional arithmetic is exact, why would you ever use floating-point numbers instead? See section Modes Tutorial Exercise 4. (*)
Typing m f doesn't change any existing values in the stack. In the above example, we had to Undo the division and do it over again when we changed to fraction mode. But if you use the evaluates-to operator you can get commands like m f to recompute for you.
1: 12 / 9 => 1.33333333333 1: 12 / 9 => 1.333 1: 12 / 9 => 4:3
. . .
' 12/9 => RET p 4 RET m f
In this example, the righthand side of the `=>' operator on the stack is recomputed when we change the precision, then again when we change to fraction mode. All `=>' expressions on the stack are recomputed every time you change any mode that might affect their values.
In this section, we explore the arithmetic and scientific functions available in the Calculator.
The standard arithmetic commands are +, -, *, /, and ^. Each normally takes two numbers from the top of the stack and pushes back a result. The n and & keys perform change-sign and reciprocal operations, respectively.
1: 5 1: 0.2 1: 5. 1: -5. 1: 5.
. . . . .
5 & & n n
You can apply a "binary operator" like + across any number of stack entries by giving it a numeric prefix. You can also apply it pairwise to several stack elements along with the top one if you use a negative prefix.
3: 2 1: 9 3: 2 4: 2 3: 12
2: 3 . 2: 3 3: 3 2: 13
1: 4 1: 4 2: 4 1: 14
. . 1: 10 .
.
2 RET 3 RET 4 M-3 + U 10 M-- M-3 +
You can apply a "unary operator" like & to the top n stack entries with a numeric prefix, too.
3: 2 3: 0.5 3: 0.5
2: 3 2: 0.333333333333 2: 3.
1: 4 1: 0.25 1: 4.
. . .
2 RET 3 RET 4 M-3 & M-2 &
Notice that the results here are left in floating-point form. We can convert them back to integers by pressing F, the "floor" function. This function rounds down to the next lower integer. There is also R, which rounds to the nearest integer.
7: 2. 7: 2 7: 2
6: 2.4 6: 2 6: 2
5: 2.5 5: 2 5: 3
4: 2.6 4: 2 4: 3
3: -2. 3: -2 3: -2
2: -2.4 2: -3 2: -2
1: -2.6 1: -3 1: -3
. . .
M-7 F U M-7 R
Since dividing-and-flooring (i.e., "integer quotient") is such a common operation, Calc provides a special command for that purpose, the backslash \. Another common arithmetic operator is %, which computes the remainder that would arise from a \ operation, i.e., the "modulo" of two numbers. For example,
2: 1234 1: 12 2: 1234 1: 34
1: 100 . 1: 100 .
. .
1234 RET 100 \ U %
These commands actually work for any real numbers, not just integers.
2: 3.1415 1: 3 2: 3.1415 1: 0.1415
1: 1 . 1: 1 .
. .
3.1415 RET 1 \ U %
(*) Exercise 1. The \ command would appear to be a frill, since you could always do the same thing with / F. Think of a situation where this is not true---/ F would be inadequate. Now think of a way you could get around the problem if Calc didn't provide a \ command. See section Arithmetic Tutorial Exercise 1. (*)
We've already seen the Q (square root) and S (sine) commands. Other commands along those lines are C (cosine), T (tangent), E (e^x) and L (natural logarithm). These can be modified by the I (inverse) and H (hyperbolic) prefix keys.
Let's compute the sine and cosine of an angle, and verify the identity @c{$\sin^2x + \cos^2x = 1$} sin(x)^2 + cos(x)^2 = 1. We'll arbitrarily pick -64 degrees as a good value for x. With the angular mode set to degrees (type m d), do:
2: -64 2: -64 2: -0.89879 2: -0.89879 1: 1.
1: -64 1: -0.89879 1: -64 1: 0.43837 .
. . . .
64 n RET RET S TAB C f h
(For brevity, we're showing only five digits of the results here. You can of course do these calculations to any precision you like.)
Remember, f h is the calc-hypot, or square-root of sum
of squares, command.
Another identity is @c{$\displaystyle\tan x = {\sin x \over \cos x}$} tan(x) = sin(x) / cos(x).
2: -0.89879 1: -2.0503 1: -64.
1: 0.43837 . .
.
U / I T
A physical interpretation of this calculation is that if you move 0.89879 units downward and 0.43837 units to the right, your direction of motion is -64 degrees from horizontal. Suppose we move in the opposite direction, up and to the left:
2: -0.89879 2: 0.89879 1: -2.0503 1: -64.
1: 0.43837 1: -0.43837 . .
. .
U U M-2 n / I T
How can the angle be the same? The answer is that the / operation
loses information about the signs of its inputs. Because the quotient
is negative, we know exactly one of the inputs was negative, but we
can't tell which one. There is an f T [arctan2] function which
computes the inverse tangent of the quotient of a pair of numbers.
Since you feed it the two original numbers, it has enough information
to give you a full 360-degree answer.
2: 0.89879 1: 116. 3: 116. 2: 116. 1: 180.
1: -0.43837 . 2: -0.89879 1: -64. .
. 1: 0.43837 .
.
U U f T M-RET M-2 n f T -
The resulting angles differ by 180 degrees; in other words, they point in opposite directions, just as we would expect.
The META-RET we used in the third step is the "last-arguments" command. It is sort of like Undo, except that it restores the arguments of the last command to the stack without removing the command's result. It is useful in situations like this one, where we need to do several operations on the same inputs. We could have accomplished the same thing by using M-2 RET to duplicate the top two stack elements right after the U U, then a pair of M-TAB commands to cycle the 116 up around the duplicates.
A similar identity is supposed to hold for hyperbolic sines and cosines, except that it is the difference cosh(x)^2 - sinh(x)^2 that always equals one. Let's try to verify this identity.
2: -64 2: -64 2: -64 2: 9.7192e54 2: 9.7192e54
1: -64 1: -3.1175e27 1: 9.7192e54 1: -64 1: 9.7192e54
. . . . .
64 n RET RET H C 2 ^ TAB H S 2 ^
Something's obviously wrong, because when we subtract these numbers the answer will clearly be zero! But if you think about it, if these numbers did differ by one, it would be in the 55th decimal place. The difference we seek has been lost entirely to roundoff error.
We could verify this hypothesis by doing the actual calculation with, say, 60 decimal places of precision. This will be slow, but not enormously so. Try it if you wish; sure enough, the answer is 0.99999, reasonably close to 1.
Of course, a more reasonable way to verify the identity is to use a more reasonable value for x!
Some Calculator commands use the Hyperbolic prefix for other purposes. The logarithm and exponential functions, for example, work to the base e normally but use base-10 instead if you use the Hyperbolic prefix.
1: 1000 1: 6.9077 1: 1000 1: 3
. . . .
1000 L U H L
First, we mistakenly compute a natural logarithm. Then we undo and compute a common logarithm instead.
The B key computes a general base-b logarithm for any value of b.
2: 1000 1: 3 1: 1000. 2: 1000. 1: 6.9077
1: 10 . . 1: 2.71828 .
. .
1000 RET 10 B H E H P B
Here we first use B to compute the base-10 logarithm, then use the "hyperbolic" exponential as a cheap hack to recover the number 1000, then use B again to compute the natural logarithm. Note that P with the hyperbolic prefix pushes the constant e onto the stack.
You may have noticed that both times we took the base-10 logarithm of 1000, we got an exact integer result. Calc always tries to give an exact rational result for calculations involving rational numbers where possible. But when we used H E, the result was a floating-point number for no apparent reason. In fact, if we had computed 10 RET 3 ^ we would have gotten an exact integer 1000. But the H E command is rigged to generate a floating-point result all of the time so that 1000 H E will not waste time computing a thousand-digit integer when all you probably wanted was `1e1000'.
(*) Exercise 2. Find a pair of integer inputs to the B command for which Calc could find an exact rational result but doesn't. See section Arithmetic Tutorial Exercise 2. (*)
The Calculator also has a set of functions relating to combinatorics and statistics. You may be familiar with the factorial function, which computes the product of all the integers up to a given number.
1: 100 1: 93326215443... 1: 100. 1: 9.3326e157
. . . .
100 ! U c f !
Recall, the c f command converts the integer or fraction at the top of the stack to floating-point format. If you take the factorial of a floating-point number, you get a floating-point result accurate to the current precision. But if you give ! an exact integer, you get an exact integer result (158 digits long in this case).
If you take the factorial of a non-integer, Calc uses a generalized factorial function defined in terms of Euler's Gamma function gamma(n) (which is itself available as the f g command).
3: 4. 3: 24. 1: 5.5 1: 52.342777847
2: 4.5 2: 52.3427777847 . .
1: 5. 1: 120.
. .
M-3 ! M-0 DEL 5.5 f g
Here we verify the identity @c{$n! = \Gamma(n+1)$} n! = gamma(n+1).
The binomial coefficient n-choose-m@c{ or $\displaystyle {n \choose m}$} is defined by n! / m! (n-m)! for all reals n and m. The intermediate results in this formula can become quite large even if the final result is small; the k c command computes a binomial coefficient in a way that avoids large intermediate values.
The k prefix key defines several common functions out of combinatorics and number theory. Here we compute the binomial coefficient 30-choose-20, then determine its prime factorization.
2: 30 1: 30045015 1: [3, 3, 5, 7, 11, 13, 23, 29]
1: 20 . .
.
30 RET 20 k c k f
You can verify these prime factors by using v u to "unpack" this vector into 8 separate stack entries, then M-8 * to multiply them back together. The result is the original number, 30045015.
Suppose a program you are writing needs a hash table with at least 10000 entries. It's best to use a prime number as the actual size of a hash table. Calc can compute the next prime number after 10000:
1: 10000 1: 10007 1: 9973
. . .
10000 k n I k n
Just for kicks we've also computed the next prime less than 10000.
See section Financial Functions, for a description of the Calculator
commands that deal with business and financial calculations (functions
like pv, rate, and sln).
See section Binary Number Functions, to read about the commands for operating
on binary numbers (like and, xor, and lsh).
A vector is a list of numbers or other Calc data objects. Calc provides a large set of commands that operate on vectors. Some are familiar operations from vector analysis. Others simply treat a vector as a list of objects.
If you add two vectors, the result is a vector of the sums of the elements, taken pairwise.
1: [1, 2, 3] 2: [1, 2, 3] 1: [8, 8, 3]
. 1: [7, 6, 0] .
.
[1,2,3] s 1 [7 6 0] s 2 +
Note that we can separate the vector elements with either commas or spaces. This is true whether we are using incomplete vectors or algebraic entry. The s 1 and s 2 commands save these vectors so we can easily reuse them later.
If you multiply two vectors, the result is the sum of the products of the elements taken pairwise. This is called the dot product of the vectors.
2: [1, 2, 3] 1: 19
1: [7, 6, 0] .
.
r 1 r 2 *
The dot product of two vectors is equal to the product of their lengths times the cosine of the angle between them. (Here the vector is interpreted as a line from the origin (0,0,0) to the specified point in three-dimensional space.) The A (absolute value) command can be used to compute the length of a vector.
3: 19 3: 19 1: 0.550782 1: 56.579
2: [1, 2, 3] 2: 3.741657 . .
1: [7, 6, 0] 1: 9.219544
. .
M-RET M-2 A * / I C
First we recall the arguments to the dot product command, then we compute the absolute values of the top two stack entries to obtain the lengths of the vectors, then we divide the dot product by the product of the lengths to get the cosine of the angle. The inverse cosine finds that the angle between the vectors is about 56 degrees.
The cross product of two vectors is a vector whose length is the product of the lengths of the inputs times the sine of the angle between them, and whose direction is perpendicular to both input vectors. Unlike the dot product, the cross product is defined only for three-dimensional vectors. Let's double-check our computation of the angle using the cross product.
2: [1, 2, 3] 3: [-18, 21, -8] 1: [-0.52, 0.61, -0.23] 1: 56.579
1: [7, 6, 0] 2: [1, 2, 3] . .
. 1: [7, 6, 0]
.
r 1 r 2 V C s 3 M-RET M-2 A * / A I S
First we recall the original vectors and compute their cross product, which we also store for later reference. Now we divide the vector by the product of the lengths of the original vectors. The length of this vector should be the sine of the angle; sure enough, it is!
Vector-related commands generally begin with the v prefix key. Some are uppercase letters and some are lowercase. To make it easier to type these commands, the shift-V prefix key acts the same as the v key. (See section General Mode Commands, for a way to make all prefix keys have this property.)
If we take the dot product of two perpendicular vectors we expect to get zero, since the cosine of 90 degrees is zero. Let's check that the cross product is indeed perpendicular to both inputs:
2: [1, 2, 3] 1: 0 2: [7, 6, 0] 1: 0
1: [-18, 21, -8] . 1: [-18, 21, -8] .
. .
r 1 r 3 * DEL r 2 r 3 *
(*) Exercise 1. Given a vector on the top of the stack, what keystrokes would you use to normalize the vector, i.e., to reduce its length to one without changing its direction? See section Vector Tutorial Exercise 1. (*)
(*) Exercise 2. Suppose a certain particle can be at any of several positions along a ruler. You have a list of those positions in the form of a vector, and another list of the probabilities for the particle to be at the corresponding positions. Find the average position of the particle. See section Vector Tutorial Exercise 2. (*)
A matrix is just a vector of vectors, all the same length. This means you can enter a matrix using nested brackets. You can also use the semicolon character to enter a matrix. We'll show both methods here:
1: [ [ 1, 2, 3 ] 1: [ [ 1, 2, 3 ]
[ 4, 5, 6 ] ] [ 4, 5, 6 ] ]
. .
[[1 2 3] [4 5 6]] ' [1 2 3; 4 5 6] RET
We'll be using this matrix again, so type s 4 to save it now.
Note that semicolons work with incomplete vectors, but they work better in algebraic entry. That's why we use the apostrophe in the second example.
When two matrices are multiplied, the lefthand matrix must have the same number of columns as the righthand matrix has rows. Row i, column j of the result is effectively the dot product of row i of the left matrix by column j of the right matrix.
If we try to duplicate this matrix and multiply it by itself, the dimensions are wrong and the multiplication cannot take place:
1: [ [ 1, 2, 3 ] * [ [ 1, 2, 3 ]
[ 4, 5, 6 ] ] [ 4, 5, 6 ] ]
.
RET *
Though rather hard to read, this is a formula which shows the product of two matrices. The `*' function, having invalid arguments, has been left in symbolic form.
We can multiply the matrices if we transpose one of them first.
2: [ [ 1, 2, 3 ] 1: [ [ 14, 32 ] 1: [ [ 17, 22, 27 ]
[ 4, 5, 6 ] ] [ 32, 77 ] ] [ 22, 29, 36 ]
1: [ [ 1, 4 ] . [ 27, 36, 45 ] ]
[ 2, 5 ] .
[ 3, 6 ] ]
.
U v t * U TAB *
Matrix multiplication is not commutative; indeed, switching the order of the operands can even change the dimensions of the result matrix, as happened here!
If you multiply a plain vector by a matrix, it is treated as a single row or column depending on which side of the matrix it is on. The result is a plain vector which should also be interpreted as a row or column as appropriate.
2: [ [ 1, 2, 3 ] 1: [14, 32]
[ 4, 5, 6 ] ] .
1: [1, 2, 3]
.
r 4 r 1 *
Multiplying in the other order wouldn't work because the number of rows in the matrix is different from the number of elements in the vector.
(*) Exercise 1. Use `*' to sum along the rows of the above @c{$2\times3$} 2x3 matrix to get [6, 15]. Now use `*' to sum along the columns to get [5, 7, 9]. See section Matrix Tutorial Exercise 1. (*)
An identity matrix is a square matrix with ones along the diagonal and zeros elsewhere. It has the property that multiplication by an identity matrix, on the left or on the right, always produces the original matrix.
1: [ [ 1, 2, 3 ] 2: [ [ 1, 2, 3 ] 1: [ [ 1, 2, 3 ]
[ 4, 5, 6 ] ] [ 4, 5, 6 ] ] [ 4, 5, 6 ] ]
. 1: [ [ 1, 0, 0 ] .
[ 0, 1, 0 ]
[ 0, 0, 1 ] ]
.
r 4 v i 3 RET *
If a matrix is square, it is often possible to find its inverse, that is, a matrix which, when multiplied by the original matrix, yields an identity matrix. The & (reciprocal) key also computes the inverse of a matrix.
1: [ [ 1, 2, 3 ] 1: [ [ -2.4, 1.2, -0.2 ]
[ 4, 5, 6 ] [ 2.8, -1.4, 0.4 ]
[ 7, 6, 0 ] ] [ -0.73333, 0.53333, -0.2 ] ]
. .
r 4 r 2 | s 5 &
The vertical bar | concatenates numbers, vectors, and matrices together. Here we have used it to add a new row onto our matrix to make it square.
We can multiply these two matrices in either order to get an identity.
1: [ [ 1., 0., 0. ] 1: [ [ 1., 0., 0. ]
[ 0., 1., 0. ] [ 0., 1., 0. ]
[ 0., 0., 1. ] ] [ 0., 0., 1. ] ]
. .
M-RET * U TAB *
Matrix inverses are related to systems of linear equations in algebra. Suppose we had the following set of equations:
This can be cast into the matrix equation,
We can solve this system of equations by multiplying both sides by the inverse of the matrix. Calc can do this all in one step:
2: [6, 2, 3] 1: [-12.6, 15.2, -3.93333]
1: [ [ 1, 2, 3 ] .
[ 4, 5, 6 ]
[ 7, 6, 0 ] ]
.
[6,2,3] r 5 /
The result is the [a, b, c] vector that solves the equations. (Dividing by a square matrix is equivalent to multiplying by its inverse.)
Let's verify this solution:
2: [ [ 1, 2, 3 ] 1: [6., 2., 3.]
[ 4, 5, 6 ] .
[ 7, 6, 0 ] ]
1: [-12.6, 15.2, -3.93333]
.
r 5 TAB *
Note that we had to be careful about the order in which we multiplied the matrix and vector. If we multiplied in the other order, Calc would assume the vector was a row vector in order to make the dimensions come out right, and the answer would be incorrect. If you don't feel safe letting Calc take either interpretation of your vectors, use explicit @c{$N\times1$} Nx1 or @c{$1\times N$} 1xN matrices instead. In this case, you would enter the original column vector as `[[6], [2], [3]]' or `[6; 2; 3]'.
(*) Exercise 2. Algebraic entry allows you to make vectors and matrices that include variables. Solve the following system of equations to get expressions for x and y in terms of a and b.
See section Matrix Tutorial Exercise 2. (*)
(*) Exercise 3. A system of equations is "over-determined" if it has more equations than variables. It is often the case that there are no values for the variables that will satisfy all the equations at once, but it is still useful to find a set of values which "nearly" satisfy all the equations. In terms of matrix equations, you can't solve A X = B directly because the matrix A is not square for an over-determined system. Matrix inversion works only for square matrices. One common trick is to multiply both sides on the left by the transpose of A: Now @c{$A^T A$} trn(A)*A is a square matrix so a solution is possible. It turns out that the X vector you compute in this way will be a "least-squares" solution, which can be regarded as the "closest" solution to the set of equations. Use Calc to solve the following over-determined system:
See section Matrix Tutorial Exercise 3. (*)
Although Calc has a number of features for manipulating vectors and matrices as mathematical objects, you can also treat vectors as simple lists of values. For example, we saw that the k f command returns a vector which is a list of the prime factors of a number.
You can pack and unpack stack entries into vectors:
3: 10 1: [10, 20, 30] 3: 10
2: 20 . 2: 20
1: 30 1: 30
. .
M-3 v p v u
You can also build vectors out of consecutive integers, or out of many copies of a given value:
1: [1, 2, 3, 4] 2: [1, 2, 3, 4] 2: [1, 2, 3, 4]
. 1: 17 1: [17, 17, 17, 17]
. .
v x 4 RET 17 v b 4 RET
You can apply an operator to every element of a vector using the map command.
1: [17, 34, 51, 68] 1: [289, 1156, 2601, 4624] 1: [17, 34, 51, 68]
. . .
V M * 2 V M ^ V M Q
In the first step, we multiply the vector of integers by the vector of 17's elementwise. In the second step, we raise each element to the power two. (The general rule is that both operands must be vectors of the same length, or else one must be a vector and the other a plain number.) In the final step, we take the square root of each element.
(*) Exercise 1. Compute a vector of powers of two from @c{$2^{-4}$} 2^-4 to 2^4. See section List Tutorial Exercise 1. (*)
You can also reduce a binary operator across a vector. For example, reducing `*' computes the product of all the elements in the vector:
1: 123123 1: [3, 7, 11, 13, 41] 1: 123123
. . .
123123 k f V R *
In this example, we decompose 123123 into its prime factors, then multiply those factors together again to yield the original number.
We could compute a dot product "by hand" using mapping and reduction:
2: [1, 2, 3] 1: [7, 12, 0] 1: 19
1: [7, 6, 0] . .
.
r 1 r 2 V M * V R +
Recalling two vectors from the previous section, we compute the sum of pairwise products of the elements to get the same answer for the dot product as before.
A slight variant of vector reduction is the accumulate operation, V U. This produces a vector of the intermediate results from a corresponding reduction. Here we compute a table of factorials:
1: [1, 2, 3, 4, 5, 6] 1: [1, 2, 6, 24, 120, 720]
. .
v x 6 RET V U *
Calc allows vectors to grow as large as you like, although it gets rather slow if vectors have more than about a hundred elements. Actually, most of the time is spent formatting these large vectors for display, not calculating on them. Try the following experiment (if your computer is very fast you may need to substitute a larger vector size).
1: [1, 2, 3, 4, ... 1: [2, 3, 4, 5, ...
. .
v x 500 RET 1 V M +
Now press v . (the letter v, then a period) and try the experiment again. In v . mode, long vectors are displayed "abbreviated" like this:
1: [1, 2, 3, ..., 500] 1: [2, 3, 4, ..., 501]
. .
v x 500 RET 1 V M +
(where now the `...' is actually part of the Calc display). You will find both operations are now much faster. But notice that even in v . mode, the full vectors are still shown in the Trail. Type t . to cause the trail to abbreviate as well, and try the experiment one more time. Operations on long vectors are now quite fast! (But of course if you use t . you will lose the ability to get old vectors back using the t y command.)
An easy way to view a full vector when v . mode is active is to press ` (back-quote) to edit the vector; editing always works with the full, unabbreviated value.
As a larger example, let's try to fit a straight line to some data, using the method of least squares. (Calc has a built-in command for least-squares curve fitting, but we'll do it by hand here just to practice working with vectors.) Suppose we have the following list of values in a file we have loaded into Emacs:
x y -- --- 1.34 0.234 1.41 0.298 1.49 0.402 1.56 0.412 1.64 0.466 1.73 0.473 1.82 0.601 1.91 0.519 2.01 0.603 2.11 0.637 2.22 0.645 2.33 0.705 2.45 0.917 2.58 1.009 2.71 0.971 2.85 1.062 3.00 1.148 3.15 1.157 3.32 1.354
If you are reading this tutorial in printed form, you will find it easiest to press M-# i to enter the on-line Info version of the manual and find this table there. (Press g, then type List Tutorial, to jump straight to this section.)
Position the cursor at the upper-left corner of this table, just
to the left of the 1.34. Press C-@ to set the mark.
(On your system this may be C-2, C-SPC, or NUL.)
Now position the cursor to the lower-right, just after the 1.354.
You have now defined this region as an Emacs "rectangle." Still
in the Info buffer, type M-# r. This command
(calc-grab-rectangle) will pop you back into the Calculator, with
the contents of the rectangle you specified in the form of a matrix.
1: [ [ 1.34, 0.234 ]
[ 1.41, 0.298 ]
...
(You may wish to use v . mode to abbreviate the display of this large matrix.)
We want to treat this as a pair of lists. The first step is to transpose this matrix into a pair of rows. Remember, a matrix is just a vector of vectors. So we can unpack the matrix into a pair of row vectors on the stack.
1: [ [ 1.34, 1.41, 1.49, ... ] 2: [1.34, 1.41, 1.49, ... ]
[ 0.234, 0.298, 0.402, ... ] ] 1: [0.234, 0.298, 0.402, ... ]
. .
v t v u
Let's store these in quick variables 1 and 2, respectively.
1: [1.34, 1.41, 1.49, ... ] .
.
t 2 t 1
(Recall that t 2 is a variant of s 2 that removes the stored value from the stack.)
In a least squares fit, the slope m is given by the formula
where @c{$\sum x$}
sum(x) represents the sum of all the values of x.
While there is an actual sum function in Calc, it's easier to
sum a vector using a simple reduction. First, let's compute the four
different sums that this formula uses.
1: 41.63 1: 98.0003
. .
r 1 V R + t 3 r 1 2 V M ^ V R + t 4
1: 13.613 1: 33.36554
. .
r 2 V R + t 5 r 1 r 2 V M * V R + t 6
Finally, we also need N, the number of data points. This is just the length of either of our lists.
1: 19
.
r 1 v l t 7
(That's v followed by a lower-case l.)
Now we grind through the formula:
1: 633.94526 2: 633.94526 1: 67.23607
. 1: 566.70919 .
.
r 7 r 6 * r 3 r 5 * -
2: 67.23607 3: 67.23607 2: 67.23607 1: 0.52141679
1: 1862.0057 2: 1862.0057 1: 128.9488 .
. 1: 1733.0569 .
.
r 7 r 4 * r 3 2 ^ - / t 8
That gives us the slope m. The y-intercept b can now be found with the simple formula,
1: 13.613 2: 13.613 1: -8.09358 1: -0.425978
. 1: 21.70658 . .
.
r 5 r 8 r 3 * - r 7 / t 9
Let's "plot" this straight line approximation, @c{$y \approx m x + b$} m x + b, and compare it with the original data.
1: [0.699, 0.735, ... ] 1: [0.273, 0.309, ... ]
. .
r 1 r 8 * r 9 + s 0
Notice that multiplying a vector by a constant, and adding a constant to a vector, can be done without mapping commands since these are common operations from vector algebra. As far as Calc is concerned, we've just been doing geometry in 19-dimensional space!
We can subtract this vector from our original y vector to get a feel for the error of our fit. Let's find the maximum error:
1: [0.0387, 0.0112, ... ] 1: [0.0387, 0.0112, ... ] 1: 0.0897
. . .
r 2 - V M A V R X
First we compute a vector of differences, then we take the absolute
values of these differences, then we reduce the max function
across the vector. (The max function is on the two-key sequence
f x; because it is so common to use max in a vector
operation, the letters X and N are also accepted for
max and min in this context. In general, you answer
the V M or V R prompt with the actual key sequence that
invokes the function you want. You could have typed V R f x or
even V R x max RET if you had preferred.)
If your system has the GNUPLOT program, you can see graphs of your data and your straight line to see how well they match. (If you have GNUPLOT 3.0, the following instructions will work regardless of the kind of display you have. Some GNUPLOT 2.0, non-X-windows systems may require additional steps to view the graphs.)
Let's start by plotting the original data. Recall the "x" and "y" vectors onto the stack and press g f. This "fast" graphing command does everything you need to do for simple, straightforward plotting of data.
2: [1.34, 1.41, 1.49, ... ]
1: [0.234, 0.298, 0.402, ... ]
.
r 1 r 2 g f
If all goes well, you will shortly get a new window containing a graph of the data. (If not, contact your GNUPLOT or Calc installer to find out what went wrong.) In the X window system, this will be a separate graphics window. For other kinds of displays, the default is to display the graph in Emacs itself using rough character graphics. Press q when you are done viewing the character graphics.
Next, let's add the line we got from our least-squares fit:
2: [1.34, 1.41, 1.49, ... ]
1: [0.273, 0.309, 0.351, ... ]
.
DEL r 0 g a g p
It's not very useful to get symbols to mark the data points on this second curve; you can type g S g p to remove them. Type g q when you are done to remove the X graphics window and terminate GNUPLOT.
(*) Exercise 2. An earlier exercise showed how to do least squares fitting to a general system of equations. Our 19 data points are really 19 equations of the form y_i = m x_i + b for different pairs of (x_i,y_i). Use the matrix-transpose method to solve for m and b, duplicating the above result. See section List Tutorial Exercise 2. (*)
(*) Exercise 3. If the input data do not form a
rectangle, you can use M-# g (calc-grab-region)
to grab the data the way Emacs normally works with regions--it reads
left-to-right, top-to-bottom, treating line breaks the same as spaces.
Use this command to find the geometric mean of the following numbers.
(The geometric mean is the nth root of the product of n numbers.)
2.3 6 22 15.1 7 15 14 7.5 2.5
The M-# g command accepts numbers separated by spaces or commas, with or without surrounding vector brackets. See section List Tutorial Exercise 3. (*)
1: [1, 2, 3, 4, 5, 6, 7] 1: [0, 1, 2, 3, 4, 5, 6]
. .
v x 7 RET 1 -
1: [1, -6, 15, -20, 15, -6, 1] 1: 0
. .
V M ' (-1)^$ choose(6,$) RET V R +
The V M ' command prompts you to enter any algebraic expression to define the function to map over the vector. The symbol `$' inside this expression represents the argument to the function. The Calculator applies this formula to each element of the vector, substituting each element's value for the `$' sign(s) in turn.
To define a two-argument function, use `$$' for the first argument and `$' for the second: V M ' $$-$ RET is equivalent to V M -. This is analogous to regular algebraic entry, where `$$' would refer to the next-to-top stack entry and `$' would refer to the top stack entry, and ' $$-$ RET would act exactly like -.
Notice that the V M ' command has recorded two things in the trail: The result, as usual, and also a funny-looking thing marked `oper' that represents the operator function you typed in. The function is enclosed in `< >' brackets, and the argument is denoted by a `#' sign. If there were several arguments, they would be shown as `#1', `#2', and so on. (For example, V M ' $$-$ will put the function `<#1 - #2>' on the trail.) This object is a "nameless function"; you can use nameless `< >' notation to answer the V M ' prompt if you like. Nameless function notation has the interesting, occasionally useful property that a nameless function is not actually evaluated until it is used. For example, V M ' $+random(2.0) evaluates `random(2.0)' once and adds that random number to all elements of the vector, but V M ' <#+random(2.0)> evaluates the `random(2.0)' separately for each vector element.
Another group of operators that are often useful with V M are the relational operators: a =, for example, compares two numbers and gives the result 1 if they are equal, or 0 if not. Similarly, a < checks for one number being less than another.
Other useful vector operations include v v, to reverse a vector end-for-end; V S, to sort the elements of a vector into increasing order; and v r and v c, to extract one row or column of a matrix, or (in both cases) to extract one element of a plain vector. With a negative argument, v r and v c instead delete one row, column, or vector element.
(*) Exercise 4. The kth divisor function is the sum of the kth powers of all the divisors of an integer n. Figure out a method for computing the divisor function for reasonably small values of n. As a test, the 0th and 1st divisor functions of 30 are 8 and 72, respectively. See section List Tutorial Exercise 4. (*)
(*) Exercise 5. The k f command produces a list of prime factors for a number. Sometimes it is important to know that a number is square-free, i.e., that no prime occurs more than once in its list of prime factors. Find a sequence of keystrokes to tell if a number is square-free; your method should leave 1 on the stack if it is, or 0 if it isn't. See section List Tutorial Exercise 5. (*)
(*) Exercise 6. Build a list of lists that looks like the following diagram. (You may wish to use the v / command to enable multi-line display of vectors.)
1: [ [1],
[1, 2],
[1, 2, 3],
[1, 2, 3, 4],
[1, 2, 3, 4, 5],
[1, 2, 3, 4, 5, 6] ]
See section List Tutorial Exercise 6. (*)
(*) Exercise 7. Build the following list of lists.
1: [ [0],
[1, 2],
[3, 4, 5],
[6, 7, 8, 9],
[10, 11, 12, 13, 14],
[15, 16, 17, 18, 19, 20] ]
See section List Tutorial Exercise 7. (*)
(*) Exercise 8. Compute a list of values of Bessel's J1 function `besJ(1,x)' for x from 0 to 5 in steps of 0.25. Find the value of x (from among the above set of values) for which `besJ(1,x)' is a maximum. Use an "automatic" method, i.e., just reading along the list by hand to find the largest value is not allowed! (There is an a X command which does this kind of thing automatically; see section Numerical Solutions.) See section List Tutorial Exercise 8. (*)
(*) Exercise 9. You are given an integer in the range 0 <= N < 10^m for m=12 (i.e., an integer of less than twelve digits). Convert this integer into a vector of m digits, each in the range from 0 to 9. In vector-of-digits notation, add one to this integer to produce a vector of m+1 digits (since there could be a carry out of the most significant digit). Convert this vector back into a regular integer. A good integer to try is 25129925999. See section List Tutorial Exercise 9. (*)
(*) Exercise 10. Your friend Joe tried to use V R a = to test if all numbers in a list were equal. What happened? How would you do this test? See section List Tutorial Exercise 10. (*)
(*) Exercise 11. The area of a circle of radius one is @c{$\pi$} pi. The area of the @c{$2\times2$} 2x2 square that encloses that circle is 4. So if we throw N darts at random points in the square, about @c{$\pi/4$} pi/4 of them will land inside the circle. This gives us an entertaining way to estimate the value of @c{$\pi$} pi. The k r command picks a random number between zero and the value on the stack. We could get a random floating-point number between -1 and 1 by typing 2.0 k r 1 -. Build a vector of 100 random (x,y) points in this square, then use vector mapping and reduction to count how many points lie inside the unit circle. Hint: Use the v b command. See section List Tutorial Exercise 11. (*)
(*) Exercise 12. The matchstick problem provides another way to calculate @c{$\pi$} pi. Say you have an infinite field of vertical lines with a spacing of one inch. Toss a one-inch matchstick onto the field. The probability that the matchstick will land crossing a line turns out to be @c{$2/\pi$} 2/pi. Toss 100 matchsticks to estimate pi. (If you want still more fun, the probability that the GCD (k g) of two large integers is one turns out to be @c{$6/\pi^2$} 6/pi^2. That provides yet another way to estimate @c{$\pi$} pi.) See section List Tutorial Exercise 12. (*)
(*) Exercise 13. An algebraic entry of a string in double-quote marks, `"hello"', creates a vector of the numerical (ASCII) codes of the characters (here, [104, 101, 108, 108, 111]). Sometimes it is convenient to compute a hash code of a string, which is just an integer that represents the value of that string. Two equal strings have the same hash code; two different strings probably have different hash codes. (For example, Calc has over 400 function names, but Emacs can quickly find the definition for any given name because it has sorted the functions into "buckets" by their hash codes. Sometimes a few names will hash into the same bucket, but it is easier to search among a few names than among all the names.) One popular hash function is computed as follows: First set h = 0. Then, for each character from the string in turn, set h = 3h + c_i where c_i is the character's ASCII code. If we have 511 buckets, we then take the hash code modulo 511 to get the bucket number. Develop a simple command or commands for converting string vectors into hash codes. The hash code for `"Testing, 1, 2, 3"' is 1960915098, which modulo 511 is 121. See section List Tutorial Exercise 13. (*)
(*) Exercise 14. The H V R and H V U
commands do nested function evaluations. H V U takes a starting
value and a number of steps n from the stack; it then applies the
function you give to the starting value 0, 1, 2, up to n times
and returns a vector of the results. Use this command to create a
"random walk" of 50 steps. Start with the two-dimensional point
(0,0); then take one step a random distance between -1 and 1
in both x and y; then take another step, and so on. Use the
g f command to display this random walk. Now modify your random
walk to walk a unit distance, but in a random direction, at each step.
(Hint: The sincos function returns a vector of the cosine and
sine of an angle.) See section List Tutorial Exercise 14. (*)
Calc understands a variety of data types as well as simple numbers. In this section, we'll experiment with each of these types in turn.
The numbers we've been using so far have mainly been either integers or floats. We saw that floats are usually a good approximation to the mathematical concept of real numbers, but they are only approximations and are susceptible to roundoff error. Calc also supports fractions, which can exactly represent any rational number.
1: 3628800 2: 3628800 1: 518400:7 1: 518414:7 1: 7:518414
. 1: 49 . . .
.
10 ! 49 RET : 2 + &
The : command divides two integers to get a fraction; / would normally divide integers to get a floating-point result. Notice we had to type RET between the 49 and the : since the : would otherwise be interpreted as part of a fraction beginning with 49.
You can convert between floating-point and fractional format using c f and c F:
1: 1.35027217629e-5 1: 7:518414
. .
c f c F
The c F command replaces a floating-point number with the "simplest" fraction whose floating-point representation is the same, to within the current precision.
1: 3.14159265359 1: 1146408:364913 1: 3.1416 1: 355:113
. . . .
P c F DEL p 5 RET P c F
(*) Exercise 1. A calculation has produced the result 1.26508260337. You suspect it is the square root of the product of @c{$\pi$} pi and some rational number. Is it? (Be sure to allow for roundoff error!) See section Types Tutorial Exercise 1. (*)
Complex numbers can be stored in both rectangular and polar form.
1: -9 1: (0, 3) 1: (3; 90.) 1: (6; 90.) 1: (2.4495; 45.)
. . . . .
9 n Q c p 2 * Q
The square root of -9 is by default rendered in rectangular form (0 + 3i), but we can convert it to polar form (3 with a phase angle of 90 degrees). All the usual arithmetic and scientific operations are defined on both types of complex numbers.
Another generalized kind of number is infinity. Infinity
isn't really a number, but it can sometimes be treated like one.
Calc uses the symbol inf to represent positive infinity,
i.e., a value greater than any real number. Naturally, you can
also write `-inf' for minus infinity, a value less than any
real number. The word inf can only be input using
algebraic entry.
2: inf 2: -inf 2: -inf 2: -inf 1: nan
1: -17 1: -inf 1: -inf 1: inf .
. . . .
' inf RET 17 n * RET 72 + A +
Since infinity is infinitely large, multiplying it by any finite
number (like -17) has no effect, except that since -17
is negative, it changes a plus infinity to a minus infinity.
("A huge positive number, multiplied by -17, yields a huge
negative number.") Adding any finite number to infinity also
leaves it unchanged. Taking an absolute value gives us plus
infinity again. Finally, we add this plus infinity to the minus
infinity we had earlier. If you work it out, you might expect
the answer to be -72 for this. But the 72 has been completely
lost next to the infinities; by the time we compute `inf - inf'
the finite difference between them, if any, is indetectable.
So we say the result is indeterminate, which Calc writes
with the symbol nan (for Not A Number).
Dividing by zero is normally treated as an error, but you can get Calc to write an answer in terms of infinity by pressing m i to turn on "infinite mode."
3: nan 2: nan 2: nan 2: nan 1: nan
2: 1 1: 1 / 0 1: uinf 1: uinf .
1: 0 . . .
.
1 RET 0 / m i U / 17 n * +
Dividing by zero normally is left unevaluated, but after m i
it instead gives an infinite result. The answer is actually
uinf, "undirected infinity." If you look at a graph of
1 / x around x = 0, you'll see that it goes toward
plus infinity as you approach zero from above, but toward minus
infinity as you approach from below. Since we said only 1 / 0,
Calc knows that the answer is infinite but not in which direction.
That's what uinf means. Notice that multiplying uinf
by a negative number still leaves plain uinf; there's no
point in saying `-uinf' because the sign of uinf is
unknown anyway. Finally, we add uinf to our nan,
yielding nan again. It's easy to see that, because
nan means "totally unknown" while uinf means
"unknown sign but known to be infinite," the more mysterious
nan wins out when it is combined with uinf, or, for
that matter, with anything else.
(*) Exercise 2. Predict what Calc will answer for each of these formulas: `inf / inf', `exp(inf)', `exp(-inf)', `sqrt(-inf)', `sqrt(uinf)', `abs(uinf)', `ln(0)'. See section Types Tutorial Exercise 2. (*)
(*) Exercise 3. We saw that `inf - inf = nan',
which stands for an unknown value. Can nan stand for
a complex number? Can it stand for infinity?
See section Types Tutorial Exercise 3. (*)
HMS forms represent a value in terms of hours, minutes, and seconds.
1: 2@ 30' 0" 1: 3@ 30' 0" 2: 3@ 30' 0" 1: 2.
. . 1: 1@ 45' 0." .
.
2@ 30' RET 1 + RET 2 / /
HMS forms can also be used to hold angles in degrees, minutes, and seconds.
1: 0.5 1: 26.56505 1: 26@ 33' 54.18" 1: 0.44721
. . . .
0.5 I T c h S
First we convert the inverse tangent of 0.5 to degrees-minutes-seconds form, then we take the sine of that angle. Note that the trigonometric functions will accept HMS forms directly as input.
(*) Exercise 4. The Beatles' Abbey Road is 47 minutes and 26 seconds long, and contains 17 songs. What is the average length of a song on Abbey Road? If the Extended Disco Version of Abbey Road added 20 seconds to the length of each song, how long would the album be? See section Types Tutorial Exercise 4. (*)
A date form represents a date, or a date and time. Dates must be entered using algebraic entry. Date forms are surrounded by `< >' symbols; most standard formats for dates are recognized.
2: <Sun Jan 13, 1991> 1: 2.25
1: <6:00pm Thu Jan 10, 1991> .
.
' <13 Jan 1991>, <1/10/91, 6pm> RET -
In this example, we enter two dates, then subtract to find the number of days between them. It is also possible to add an HMS form or a number (of days) to a date form to get another date form.
1: <4:45:59pm Mon Jan 14, 1991> 1: <2:50:59am Thu Jan 17, 1991>
. .
t N 2 + 10@ 5' +
The t N ("now") command pushes the current date and time on the stack; then we add two days, ten hours and five minutes to the date and time. Other date-and-time related commands include t J, which does Julian day conversions, t W, which finds the beginning of the week in which a date form lies, and t I, which increments a date by one or several months. See section Date Arithmetic, for more.
(*) Exercise 5. How many days until the next Friday the 13th? See section Types Tutorial Exercise 5. (*)
(*) Exercise 6. How many leap years will there be between now and the year 10001 A.D.? See section Types Tutorial Exercise 6. (*)
An error form represents a mean value with an attached standard deviation, or error estimate. Suppose our measurements indicate that a certain telephone pole is about 30 meters away, with an estimated error of 1 meter, and 8 meters tall, with an estimated error of 0.2 meters. What is the slope of a line from here to the top of the pole, and what is the equivalent angle in degrees?
1: 8 +/- 0.2 2: 8 +/- 0.2 1: 0.266 +/- 0.011 1: 14.93 +/- 0.594
. 1: 30 +/- 1 . .
.
8 p .2 RET 30 p 1 / I T
This means that the angle is about 15 degrees, and, assuming our original error estimates were valid standard deviations, there is about a 60% chance that the result is correct within 0.59 degrees.
(*) Exercise 7. The volume of a torus (a donut shape) is 2 pi^2 R r^2 where R is the radius of the circle that defines the center of the tube and r is the radius of the tube itself. Suppose R is 20 cm and r is 4 cm, each known to within 5 percent. What is the volume and the relative uncertainty of the volume? See section Types Tutorial Exercise 7. (*)
An interval form represents a range of values. While an error form is best for making statistical estimates, intervals give you exact bounds on an answer. Suppose we additionally know that our telephone pole is definitely between 28 and 31 meters away, and that it is between 7.7 and 8.1 meters tall.
1: [7.7 .. 8.1] 2: [7.7 .. 8.1] 1: [0.24 .. 0.28] 1: [13.9 .. 16.1]
. 1: [28 .. 31] . .
.
[ 7.7 .. 8.1 ] [ 28 .. 31 ] / I T
If our bounds were correct, then the angle to the top of the pole is sure to lie in the range shown.
The square brackets around these intervals indicate that the endpoints themselves are allowable values. In other words, the distance to the telephone pole is between 28 and 31, inclusive. You can also make an interval that is exclusive of its endpoints by writing parentheses instead of square brackets. You can even make an interval which is inclusive ("closed") on one end and exclusive ("open") on the other.
1: [1 .. 10) 1: (0.1 .. 1] 2: (0.1 .. 1] 1: (0.2 .. 3)
. . 1: [2 .. 3) .
.
[ 1 .. 10 ) & [ 2 .. 3 ) *
The Calculator automatically keeps track of which end values should be open and which should be closed. You can also make infinite or semi-infinite intervals by using `-inf' or `inf' for one or both endpoints.
(*) Exercise 8. What answer would you expect from `1 / (0 .. 10)'? What about `1 / (-10 .. 0)'? What about `1 / [0 .. 10]' (where the interval actually includes zero)? What about `1 / (-10 .. 10)'? See section Types Tutorial Exercise 8. (*)
(*) Exercise 9. Two easy ways of squaring a number are RET * and 2 ^. Normally these produce the same answer. Would you expect this still to hold true for interval forms? If not, which of these will result in a larger interval? See section Types Tutorial Exercise 9. (*)
A modulo form is used for performing arithmetic modulo M. For example, arithmetic involving time is generally done modulo 12 or 24 hours.
1: 17 mod 24 1: 3 mod 24 1: 21 mod 24 1: 9 mod 24
. . . .
17 M 24 RET 10 + n 5 /
In this last step, Calc has found a new number which, when multiplied by 5 modulo 24, produces the original number, 21. If M is prime it is always possible to find such a number. For non-prime M like 24, it is only sometimes possible.
1: 10 mod 24 1: 16 mod 24 1: 1000000... 1: 16
. . . .
10 M 24 RET 100 ^ 10 RET 100 ^ 24 %
These two calculations get the same answer, but the first one is much more efficient because it avoids the huge intermediate value that arises in the second one.
(*) Exercise 10. A theorem of Pierre de Fermat says that @c{\w{$x^{n-1} \bmod n = 1$}} x^(n-1) mod n = 1 if n is a prime number and x is an integer less than n. If n is not a prime number, this will not be true for most values of x. Thus we can test informally if a number is prime by trying this formula for several values of x. Use this test to tell whether the following numbers are prime: 811749613, 15485863. See section Types Tutorial Exercise 10. (*)
It is possible to use HMS forms as parts of error forms, intervals,
modulo forms, or as the phase part of a polar complex number.
For example, the calc-time command pushes the current time
of day on the stack as an HMS/modulo form.
1: 17@ 34' 45" mod 24@ 0' 0" 1: 6@ 22' 15" mod 24@ 0' 0"
. .
x time RET n
This calculation tells me it is six hours and 22 minutes until midnight.
(*) Exercise 11. A rule of thumb is that one year is about @c{$\pi \times 10^7$} pi * 10^7 seconds. What time will it be that many seconds from right now? See section Types Tutorial Exercise 11. (*)
(*) Exercise 12. You are preparing to order packaging for the CD release of the Extended Disco Version of Abbey Road. You are told that the songs will actually be anywhere from 20 to 60 seconds longer than the originals. One CD can hold about 75 minutes of music. Should you order single or double packages? See section Types Tutorial Exercise 12. (*)
Another kind of data the Calculator can manipulate is numbers with units. This isn't strictly a new data type; it's simply an application of algebraic expressions, where we use variables with suggestive names like `cm' and `in' to represent units like centimeters and inches.
1: 2 in 1: 5.08 cm 1: 0.027778 fath 1: 0.0508 m
. . . .
' 2in RET u c cm RET u c fath RET u b
We enter the quantity "2 inches" (actually an algebraic expression which means two times the variable `in'), then we convert it first to centimeters, then to fathoms, then finally to "base" units, which in this case means meters.
1: 9 acre 1: 3 sqrt(acre) 1: 190.84 m 1: 190.84 m + 30 cm
. . . .
' 9 acre RET Q u s ' $+30 cm RET
1: 191.14 m 1: 36536.3046 m^2 1: 365363046 cm^2
. . .
u s 2 ^ u c cgs
Since units expressions are really just formulas, taking the square
root of `acre' is undefined. After all, acre might be an
algebraic variable that you will someday assign a value. We use the
"units-simplify" command to simplify the expression with variables
being interpreted as unit names.
In the final step, we have converted not to a particular unit, but to a units system. The "cgs" system uses centimeters instead of meters as its standard unit of length.
There is a wide variety of units defined in the Calculator.
1: 55 mph 1: 88.5139 kph 1: 88.5139 km / hr 1: 8.201407e-8 c
. . . .
' 55 mph RET u c kph RET u c km/hr RET u c c RET
We express a speed first in miles per hour, then in kilometers per hour, then again using a slightly more explicit notation, then finally in terms of fractions of the speed of light.
Temperature conversions are a bit more tricky. There are two ways to interpret "20 degrees Fahrenheit"---it could mean an actual temperature, or it could mean a change in temperature. For normal units there is no difference, but temperature units have an offset as well as a scale factor and so there must be two explicit commands for them.
1: 20 degF 1: 11.1111 degC 1: -20:3 degC 1: -6.666 degC
. . . .
' 20 degF RET u c degC RET U u t degC RET c f
First we convert a change of 20 degrees Fahrenheit into an equivalent change in degrees Celsius (or Centigrade). Then, we convert the absolute temperature 20 degrees Fahrenheit into Celsius. Since this comes out as an exact fraction, we then convert to floating-point for easier comparison with the other result.
For simple unit conversions, you can put a plain number on the stack. Then u c and u t will prompt for both old and new units. When you use this method, you're responsible for remembering which numbers are in which units:
1: 55 1: 88.5139 1: 8.201407e-8
. . .
55 u c mph RET kph RET u c km/hr RET c RET
To see a complete list of built-in units, type u v. Press M-# c again to re-enter the Calculator when you're done looking at the units table.
(*) Exercise 13. How many seconds are there really in a year? See section Types Tutorial Exercise 13. (*)
(*) Exercise 14. Supercomputer designs are limited by the speed of light (and of electricity, which is nearly as fast). Suppose a computer has a 4.1 ns (nanosecond) clock cycle, and its cabinet is one meter across. Is speed of light going to be a significant factor in its design? See section Types Tutorial Exercise 14. (*)
(*) Exercise 15. Sam the Slug normally travels about five yards in an hour. He has obtained a supply of Power Pills; each Power Pill he eats doubles his speed. How many Power Pills can he swallow and still travel legally on most US highways? See section Types Tutorial Exercise 15. (*)
This section shows how to use Calc's algebra facilities to solve equations, do simple calculus problems, and manipulate algebraic formulas.
If you enter a formula in algebraic mode that refers to variables, the formula itself is pushed onto the stack. You can manipulate formulas as regular data objects.
1: 2 x^2 - 6 1: 6 - 2 x^2 1: (6 - 2 x^2) (3 x^2 + y)
. . .
' 2x^2-6 RET n ' 3x^2+y RET *
(*) Exercise 1. Do ' x RET Q 2 ^ and ' x RET 2 ^ Q both wind up with the same result (`x')? Why or why not? See section Algebra Tutorial Exercise 1. (*)
There are also commands for doing common algebraic operations on formulas. Continuing with the formula from the last example,
1: 18 x^2 + 6 y - 6 x^4 - 2 x^2 y 1: (18 - 2 y) x^2 - 6 x^4 + 6 y
. .
a x a c x RET
First we "expand" using the distributive law, then we "collect" terms involving like powers of x.
Let's find the value of this expression when x is 2 and y is one-half.
1: 17 x^2 - 6 x^4 + 3 1: -25
. .
1:2 s l y RET 2 s l x RET
The s l command means "let"; it takes a number from the top of the stack and temporarily assigns it as the value of the variable you specify. It then evaluates (as if by the = key) the next expression on the stack. After this command, the variable goes back to its original value, if any.
(An earlier exercise in this tutorial involved storing a value in the
variable x; if this value is still there, you will have to
unstore it with s u x RET before the above example will work
properly.)
Let's find the maximum value of our original expression when y is one-half and x ranges over all possible values. We can do this by taking the derivative with respect to x and examining values of x for which the derivative is zero. If the second derivative of the function at that value of x is negative, the function has a local maximum there.
1: 17 x^2 - 6 x^4 + 3 1: 34 x - 24 x^3
. .
U DEL s 1 a d x RET s 2
Well, the derivative is clearly zero when x is zero. To find the other root(s), let's divide through by x and then solve:
1: (34 x - 24 x^3) / x 1: 34 x / x - 24 x^3 / x 1: 34 - 24 x^2
. . .
' x RET / a x a s
1: 34 - 24 x^2 = 0 1: x = 1.19023
. .
0 a = s 3 a S x RET
Notice the use of a s to "simplify" the formula. When the default algebraic simplifications don't do enough, you can use a s to tell Calc to spend more time on the job.
Now we compute the second derivative and plug in our values of x:
1: 1.19023 2: 1.19023 2: 1.19023
. 1: 34 x - 24 x^3 1: 34 - 72 x^2
. .
a . r 2 a d x RET s 4
(The a . command extracts just the righthand side of an equation. Another method would have been to use v u to unpack the equation `x = 1.19' to `x' and `1.19', then use M-- M-2 DEL to delete the `x'.)
2: 34 - 72 x^2 1: -68. 2: 34 - 72 x^2 1: 34
1: 1.19023 . 1: 0 .
. .
TAB s l x RET U DEL 0 s l x RET
The first of these second derivatives is negative, so we know the function has a maximum value at x = 1.19023. (The function also has a local minimum at x = 0.)
When we solved for x, we got only one value even though 34 - 24 x^2 = 0 is a quadratic equation that ought to have two solutions. The reason is that a S normally returns a single "principal" solution. If it needs to come up with an arbitrary sign (as occurs in the quadratic formula) it picks +. If it needs an arbitrary integer, it picks zero. We can get a full solution by pressing H (the Hyperbolic flag) before a S.
1: 34 - 24 x^2 = 0 1: x = 1.19023 s1 1: x = -1.19023
. . .
r 3 H a S x RET s 5 1 n s l s1 RET
Calc has invented the variable `s1' to represent an unknown sign; it is supposed to be either +1 or -1. Here we have used the "let" command to evaluate the expression when the sign is negative. If we plugged this into our second derivative we would get the same, negative, answer, so x = -1.19023 is also a maximum.
To find the actual maximum value, we must plug our two values of x into the original formula.
2: 17 x^2 - 6 x^4 + 3 1: 24.08333 s1^2 - 12.04166 s1^4 + 3
1: x = 1.19023 s1 .
.
r 1 r 5 s l RET
(Here we see another way to use s l; if its input is an equation with a variable on the lefthand side, then s l treats the equation like an assignment to that variable if you don't give a variable name.)
It's clear that this will have the same value for either sign of
s1, but let's work it out anyway, just for the exercise:
2: [-1, 1] 1: [15.04166, 15.04166]
1: 24.08333 s1^2 ... .
.
[ 1 n , 1 ] TAB V M $ RET
Here we have used a vector mapping operation to evaluate the function at several values of `s1' at once. V M $ is like V M ' except that it takes the formula from the top of the stack. The formula is interpreted as a function to apply across the vector at the next-to-top stack level. Since a formula on the stack can't contain `$' signs, Calc assumes the variables in the formula stand for different arguments. It prompts you for an argument list, giving the list of all variables in the formula in alphabetical order as the default list. In this case the default is `(s1)', which is just what we want so we simply press RET at the prompt.
If there had been several different values, we could have used V R X to find the global maximum.
Calc has a built-in a P command that solves an equation using H a S and returns a vector of all the solutions. It simply automates the job we just did by hand. Applied to our original cubic polynomial, it would produce the vector of solutions [1.19023, -1.19023, 0]. (There is also an a X command which finds a local maximum of a function. It uses a numerical search method rather than examining the derivatives, and thus requires you to provide some kind of initial guess to show it where to look.)
(*) Exercise 2. Given a vector of the roots of a polynomial (such as the output of an a P command), what sequence of commands would you use to reconstruct the original polynomial? (The answer will be unique to within a constant multiple; choose the solution where the leading coefficient is one.) See section Algebra Tutorial Exercise 2. (*)
The m s command enables "symbolic mode," in which formulas like `sqrt(5)' that can't be evaluated exactly are left in symbolic form rather than giving a floating-point approximate answer. Fraction mode (m f) is also useful when doing algebra.
2: 34 x - 24 x^3 2: 34 x - 24 x^3
1: 34 x - 24 x^3 1: [sqrt(51) / 6, sqrt(51) / -6, 0]
. .
r 2 RET m s m f a P x RET
One more mode that makes reading formulas easier is "Big mode."
3
2: 34 x - 24 x
____ ____
V 51 V 51
1: [-----, -----, 0]
6 -6
.
d B
Here things like powers, square roots, and quotients and fractions are displayed in a two-dimensional pictorial form. Calc has other language modes as well, such as C mode, FORTRAN mode, and TeX mode.
2: 34*x - 24*pow(x, 3) 2: 34*x - 24*x**3
1: {sqrt(51) / 6, sqrt(51) / -6, 0} 1: /sqrt(51) / 6, sqrt(51) / -6, 0/
. .
d C d F
3: 34 x - 24 x^3
2: [{\sqrt{51} \over 6}, {\sqrt{51} \over -6}, 0]
1: {2 \over 3} \sqrt{5}
.
d T ' 2 \sqrt{5} \over 3 RET
As you can see, language modes affect both entry and display of formulas. They affect such things as the names used for built-in functions, the set of arithmetic operators and their precedences, and notations for vectors and matrices.
Notice that `sqrt(51)' may cause problems with older implementations of C and FORTRAN, which would require something more like `sqrt(51.0)'. It is always wise to check over the formulas produced by the various language modes to make sure they are fully correct.
Type m s, m f, and d N to reset these modes. (You may prefer to remain in Big mode, but all the examples in the tutorial are shown in normal mode.)
What is the area under the portion of this curve from x = 1 to 2? This is simply the integral of the function:
1: 17 x^2 - 6 x^4 + 3 1: 5.6666 x^3 - 1.2 x^5 + 3 x
. .
r 1 a i x
We want to evaluate this at our two values for x and subtract. One way to do it is again with vector mapping and reduction:
2: [2, 1] 1: [12.93333, 7.46666] 1: 5.46666 1: 5.6666 x^3 ... . . [ 2 , 1 ] TAB V M $ RET V R -
(*) Exercise 3. Find the integral from 1 to y of @c{$x \sin \pi x$} x sin(pi x) (where the sine is calculated in radians). Find the values of the integral for integers y from 1 to 5. See section Algebra Tutorial Exercise 3. (*)
Calc's integrator can do many simple integrals symbolically, but many others are beyond its capabilities. Suppose we wish to find the area under the curve @c{$\sin x \ln x$} sin(x) ln(x) over the same range of x. If you entered this formula and typed a i x RET (don't bother to try this), Calc would work for a long time but would be unable to find a solution. In fact, there is no closed-form solution to this integral. Now what do we do?
One approach would be to do the integral numerically. It is not hard to do this by hand using vector mapping and reduction. It is rather slow, though, since the sine and logarithm functions take a long time. We can save some time by reducing the working precision.
3: 10 1: [1, 1.1, 1.2, ... , 1.8, 1.9]
2: 1 .
1: 0.1
.
10 RET 1 RET .1 RET C-u v x
(Note that we have used the extended version of v x; we could also have used plain v x as follows: v x 10 RET 9 + .1 *.)
2: [1, 1.1, ... ] 1: [0., 0.084941, 0.16993, ... ]
1: sin(x) ln(x) .
.
' sin(x) ln(x) RET s 1 m r p 5 RET V M $ RET
1: 3.4195 0.34195
. .
V R + 0.1 *
(If you got wildly different results, did you remember to switch to radians mode?)
Here we have divided the curve into ten segments of equal width; approximating these segments as rectangular boxes (i.e., assuming the curve is nearly flat at that resolution), we compute the areas of the boxes (height times width), then sum the areas. (It is faster to sum first, then multiply by the width, since the width is the same for every box.)
The true value of this integral turns out to be about 0.374, so we're not doing too well. Let's try another approach.
1: sin(x) ln(x) 1: 0.84147 x - 0.84147 + 0.11957 (x - 1)^2 - ...
. .
r 1 a t x=1 RET 4 RET
Here we have computed the Taylor series expansion of the function about the point x=1. We can now integrate this polynomial approximation, since polynomials are easy to integrate.
1: 0.42074 x^2 + ... 1: [-0.0446, -0.42073] 1: 0.3761
. . .
a i x RET [ 2 , 1 ] TAB V M $ RET V R -
Better! By increasing the precision and/or asking for more terms
in the Taylor series, we can get a result as accurate as we like.
(Taylor series converge better away from singularities in the
function such as the one at ln(0), so it would also help to
expand the series about the points x=2 or x=1.5 instead
of x=1.)
(*) Exercise 4. Our first method approximated the curve by stairsteps of width 0.1; the total area was then the sum of the areas of the rectangles under these stairsteps. Our second method approximated the function by a polynomial, which turned out to be a better approximation than stairsteps. A third method is Simpson's rule, which is like the stairstep method except that the steps are not required to be flat. Simpson's rule boils down to the formula,
where n (which must be even) is the number of slices and h is the width of each slice. These are 10 and 0.1 in our example. For reference, here is the corresponding formula for the stairstep method:
Compute the integral from 1 to 2 of @c{$\sin x \ln x$} sin(x) ln(x) using Simpson's rule with 10 slices. See section Algebra Tutorial Exercise 4. (*)
Calc has a built-in a I command for doing numerical integration. It uses Romberg's method, which is a more sophisticated cousin of Simpson's rule. In particular, it knows how to keep refining the result until the current precision is satisfied.
Aside from the commands we've seen so far, Calc also provides a large set of commands for operating on parts of formulas. You indicate the desired sub-formula by placing the cursor on any part of the formula before giving a selection command. Selections won't be covered in the tutorial; see section Selecting Sub-Formulas, for details and examples.
No matter how many built-in commands Calc provided for doing algebra, there would always be something you wanted to do that Calc didn't have in its repertoire. So Calc also provides a rewrite rule system that you can use to define your own algebraic manipulations.
Suppose we want to simplify this trigonometric formula:
1: 1 / cos(x) - sin(x) tan(x)
.
' 1/cos(x) - sin(x) tan(x) RET s 1
If we were simplifying this by hand, we'd probably replace the `tan' with a `sin/cos' first, then combine over a common denominator. There is no Calc command to do the former; the a n algebra command will do the latter but we'll do both with rewrite rules just for practice.
Rewrite rules are written with the `:=' symbol.
1: 1 / cos(x) - sin(x)^2 / cos(x)
.
a r tan(a) := sin(a)/cos(a) RET
(The "assignment operator" `:=' has several uses in Calc. All by itself the formula `tan(a) := sin(a)/cos(a)' doesn't do anything, but when it is given to the a r command, that command interprets it as a rewrite rule.)
The lefthand side, `tan(a)', is called the pattern of the rewrite rule. Calc searches the formula on the stack for parts that match the pattern. Variables in a rewrite pattern are called meta-variables, and when matching the pattern each meta-variable can match any sub-formula. Here, the meta-variable `a' matched the actual variable `x'.
When the pattern part of a rewrite rule matches a part of the formula, that part is replaced by the righthand side with all the meta-variables substituted with the things they matched. So the result is `sin(x) / cos(x)'. Calc's normal algebraic simplifications then mix this in with the rest of the original formula.
To merge over a common denominator, we can use another simple rule:
1: (1 - sin(x)^2) / cos(x)
.
a r a/x + b/x := (a+b)/x RET
This rule points out several interesting features of rewrite patterns. First, if a meta-variable appears several times in a pattern, it must match the same thing everywhere. This rule detects common denominators because the same meta-variable `x' is used in both of the denominators.
Second, meta-variable names are independent from variables in the target formula. Notice that the meta-variable `x' here matches the subformula `cos(x)'; Calc never confuses the two meanings of `x'.
And third, rewrite patterns know a little bit about the algebraic properties of formulas. The pattern called for a sum of two quotients; Calc was able to match a difference of two quotients by matching `a = 1', `b = -sin(x)^2', and `x = cos(x)'.
We could just as easily have written `a/x - b/x := (a-b)/x' for
the rule. It would have worked just the same in all cases. (If we
really wanted the rule to apply only to `+' or only to `-',
we could have used the plain symbol. See section Algebraic Properties of Rewrite Rules, for some examples of this.)
One more rewrite will complete the job. We want to use the identity `sin(x)^2 + cos(x)^2 = 1', but of course we must first rearrange the identity in a way that matches our formula. The obvious rule would be `1 - sin(x)^2 := cos(x)^2', but a little thought shows that the rule `sin(x)^2 := 1 - cos(x)^2' will also work. The latter rule has a more general pattern so it will work in many other situations, too.
1: (1 + cos(x)^2 - 1) / cos(x) 1: cos(x)
. .
a r sin(x)^2 := 1 - cos(x)^2 RET a s
You may ask, what's the point of using the most general rule if you have to type it in every time anyway? The answer is that Calc allows you to store a rewrite rule in a variable, then give the variable name in the a r command. In fact, this is the preferred way to use rewrites. For one, if you need a rule once you'll most likely need it again later. Also, if the rule doesn't work quite right you can simply Undo, edit the variable, and run the rule again without having to retype it.
' tan(x) := sin(x)/cos(x) RET s t tsc RET
' a/x + b/x := (a+b)/x RET s t merge RET
' sin(x)^2 := 1 - cos(x)^2 RET s t sinsqr RET
1: 1 / cos(x) - sin(x) tan(x) 1: cos(x)
. .
r 1 a r tsc RET a r merge RET a r sinsqr RET a s
To edit a variable, type s e and the variable name, use regular Emacs editing commands as necessary, then type M-# M-# or C-c C-c to store the edited value back into the variable. You can also use s e to create a new variable if you wish.
Notice that the first time you use each rule, Calc puts up a "compiling" message briefly. The pattern matcher converts rules into a special optimized pattern-matching language rather than using them directly. This allows a r to apply even rather complicated rules very efficiently. If the rule is stored in a variable, Calc compiles it only once and stores the compiled form along with the variable. That's another good reason to store your rules in variables rather than entering them on the fly.
(*) Exercise 1. Type m s to get symbolic mode, then enter the formula `(2 + sqrt(2)) / (1 + sqrt(2))'. Using a rewrite rule, simplify this formula by multiplying both sides by the conjugate `1 - sqrt(2)'. The result will have to be expanded by the distributive law; do this with another rewrite. See section Rewrites Tutorial Exercise 1. (*)
The a r command can also accept a vector of rewrite rules, or a variable containing a vector of rules.
1: [tsc, merge, sinsqr] 1: [tan(x) := sin(x) / cos(x), ... ]
. .
' [tsc,merge,sinsqr] RET =
1: 1 / cos(x) - sin(x) tan(x) 1: cos(x)
. .
s t trig RET r 1 a r trig RET a s
Calc tries all the rules you give against all parts of the formula, repeating until no further change is possible. (The exact order in which things are tried is rather complex, but for simple rules like the ones we've used here the order doesn't really matter. See section Nested Formulas with Rewrite Rules.)
Calc actually repeats only up to 100 times, just in case your rule set has gotten into an infinite loop. You can give a numeric prefix argument to a r to specify any limit. In particular, M-1 a r does only one rewrite at a time.
1: 1 / cos(x) - sin(x)^2 / cos(x) 1: (1 - sin(x)^2) / cos(x)
. .
r 1 M-1 a r trig RET M-1 a r trig RET
You can type M-0 a r if you want no limit at all on the number of rewrites that occur.
Rewrite rules can also be conditional. Simply follow the rule with a `::' symbol and the desired condition. For example,
1: exp(2 pi i) + exp(3 pi i) + exp(4 pi i)
.
' exp(2 pi i) + exp(3 pi i) + exp(4 pi i) RET
1: 1 + exp(3 pi i) + 1
.
a r exp(k pi i) := 1 :: k % 2 = 0 RET
(Recall, `k % 2' is the remainder from dividing `k' by 2, which will be zero only when `k' is an even integer.)
An interesting point is that the variables `pi' and `i' were matched literally rather than acting as meta-variables. This is because they are special-constant variables. The special constants `e', `phi', and so on also match literally. A common error with rewrite rules is to write, say, `f(a,b,c,d,e) := g(a+b+c+d+e)', expecting to match any `f' with five arguments but in fact matching only when the fifth argument is literally `e'!
Rewrite rules provide an interesting way to define your own functions. Suppose we want to define `fib(n)' to produce the nth Fibonacci number. The first two Fibonacci numbers are each 1; later numbers are formed by summing the two preceding numbers in the sequence. This is easy to express in a set of three rules:
' [fib(1) := 1, fib(2) := 1, fib(n) := fib(n-1) + fib(n-2)] RET s t fib
1: fib(7) 1: 13
. .
' fib(7) RET a r fib RET
One thing that is guaranteed about the order that rewrites are tried is that, for any given subformula, earlier rules in the rule set will be tried for that subformula before later ones. So even though the first and third rules both match `fib(1)', we know the first will be used preferentially.
This rule set has one dangerous bug: Suppose we apply it to the formula `fib(x)'? (Don't actually try this.) The third rule will match `fib(x)' and replace it with `fib(x-1) + fib(x-2)'. Each of these will then be replaced to get `fib(x-2) + 2 fib(x-3) + fib(x-4)', and so on, expanding forever. What we really want is to apply the third rule only when `n' is an integer greater than two. Type s e fib RET, then edit the third rule to:
fib(n) := fib(n-1) + fib(n-2) :: integer(n) :: n > 2
Now:
1: fib(6) + fib(x) + fib(0) 1: 8 + fib(x) + fib(0)
. .
' fib(6)+fib(x)+fib(0) RET a r fib RET
We've created a new function, fib, and a new command,
a r fib RET, which means "evaluate all fib calls in
this formula." To make things easier still, we can tell Calc to
apply these rules automatically by storing them in the special
variable EvalRules.
1: [fib(1) := ...] . 1: [8, 13]
. .
s r fib RET s t EvalRules RET ' [fib(6), fib(7)] RET
It turns out that this rule set has the problem that it does far more work than it needs to when `n' is large. Consider the first few steps of the computation of `fib(6)':
fib(6) = fib(5) + fib(4) = fib(4) + fib(3) + fib(3) + fib(2) = fib(3) + fib(2) + fib(2) + fib(1) + fib(2) + fib(1) + 1 = ...
Note that `fib(3)' appears three times here. Unless Calc's
algebraic simplifier notices the multiple `fib(3)'s and combines
them (and, as it happens, it doesn't), this rule set does lots of
needless recomputation. To cure the problem, type s e EvalRules
to edit the rules (or just s E, a shorthand command for editing
EvalRules) and add another condition:
fib(n) := fib(n-1) + fib(n-2) :: integer(n) :: n > 2 :: remember
If a `:: remember' condition appears anywhere in a rule, then if
that rule succeeds Calc will add another rule that describes that match
to the front of the rule set. (Remembering works in any rule set, but
for technical reasons it is most effective in EvalRules.) For
example, if the rule rewrites `fib(7)' to something that evaluates
to 13, then the rule `fib(7) := 13' will be added to the rule set.
Type ' fib(8) RET to compute the eighth Fibonacci number, then type s E again to see what has happened to the rule set.
With the remember feature, our rule set can now compute
`fib(n)' in just n steps. In the process it builds
up a table of all Fibonacci numbers up to n. After we have
computed the result for a particular n, we can get it back
(and the results for all smaller n) later in just one step.
All Calc operations will run somewhat slower whenever EvalRules
contains any rules. You should type s u EvalRules RET now to
un-store the variable.
(*) Exercise 2. Sometimes it is possible to reformulate
a problem to reduce the amount of recursion necessary to solve it.
Create a rule that, in about n simple steps and without recourse
to the remember option, replaces `fib(n, 1, 1)' with
`fib(1, x, y)' where x and y are the
nth and n+1st Fibonacci numbers, respectively. This rule is
rather clunky to use, so add a couple more rules to make the "user
interface" the same as for our first version: enter `fib(n)',
get back a plain number. See section Rewrites Tutorial Exercise 2. (*)
There are many more things that rewrites can do. For example, there are `&&&' and `|||' pattern operators that create "and" and "or" combinations of rules. As one really simple example, we could combine our first two Fibonacci rules thusly:
[fib(1 ||| 2) := 1, fib(n) := ... ]
That means "fib of something matching either 1 or 2 rewrites
to 1."
You can also make meta-variables optional by enclosing them in opt.
For example, the pattern `a + b x' matches `2 + 3 x' but not
`2 + x' or `3 x' or `x'. The pattern `opt(a) + opt(b) x'
matches all of these forms, filling in a default of zero for `a'
and one for `b'.
(*) Exercise 3. Your friend Joe had `2 + 3 x' on the stack and tried to use the rule `opt(a) + opt(b) x := f(a, b, x)'. What happened? See section Rewrites Tutorial Exercise 3. (*)
(*) Exercise 4. Starting with a positive integer a, divide a by two if it is even, otherwise compute 3 a + 1. Now repeat this step over and over. A famous unproved conjecture is that for any starting a, the sequence always eventually reaches 1. Given the formula `seq(a, 0)', write a set of rules that convert this into `seq(1, n)' where n is the number of steps it took the sequence to reach the value 1. Now enhance the rules to accept `seq(a)' as a starting configuration, and to stop with just the number n by itself. Now make the result be a vector of values in the sequence, from a to 1. (The formula `x|y' appends the vectors x and y.) For example, rewriting `seq(6)' should yield the vector [6, 3, 10, 5, 16, 8, 4, 2, 1]. See section Rewrites Tutorial Exercise 4. (*)
(*) Exercise 5. Define, using rewrite rules, a function `nterms(x)' that returns the number of terms in the sum x, or 1 if x is not a sum. (A sum for our purposes is one or more non-sum terms separated by `+' or `-' signs, so that 2 - 3 (x + y) + x y is a sum of three terms.) See section Rewrites Tutorial Exercise 5. (*)
(*) Exercise 6. Calc considers the form 0^0 to be "indeterminate," and leaves it unevaluated (assuming infinite mode is not enabled). Some people prefer to define 0^0 = 1, so that the identity x^0 = 1 can safely be used for all x. Find a way to make Calc follow this convention. What happens if you now type m i to turn on infinite mode? See section Rewrites Tutorial Exercise 6. (*)
(*) Exercise 7. A Taylor series for a function is an infinite series that exactly equals the value of that function at values of x near zero.
The a t command produces a truncated Taylor series which is obtained by dropping all the terms higher than, say, x^2. Calc represents the truncated Taylor series as a polynomial in x. Mathematicians often write a truncated series using a "big-O" notation that records what was the lowest term that was truncated.
The meaning of O(x^3) is "a quantity which is negligibly small if x^3 is considered negligibly small as x goes to zero."
The exercise is to create rewrite rules that simplify sums and products of power series represented as `polynomial + O(var^n)'. For example, given `1 - x^2 / 2 + O(x^3)' and `x - x^3 / 6 + O(x^4)' on the stack, we want to be able to type * and get the result `x - 2:3 x^3 + O(x^4)'. Don't worry if the terms of the sum are rearranged or if a s needs to be typed after rewriting. (This one is rather tricky; the solution at the end of this chapter uses 6 rewrite rules. Hint: The `constant(x)' condition tests whether `x' is a number.) See section Rewrites Tutorial Exercise 7. (*)
See section Rewrite Rules, for the whole story on rewrite rules.
The Calculator is written entirely in Emacs Lisp, a highly extensible language. If you know Lisp, you can program the Calculator to do anything you like. Rewrite rules also work as a powerful programming system. But Lisp and rewrite rules take a while to master, and often all you want to do is define a new function or repeat a command a few times. Calc has features that allow you to do these things easily.
(Note that the programming commands relating to user-defined keys are not yet supported under Lucid Emacs 19.)
One very limited form of programming is defining your own functions. Calc's Z F command allows you to define a function name and key sequence to correspond to any formula. Programming commands use the shift-Z prefix; the user commands they create use the lower case z prefix.
1: 1 + x + x^2 / 2 + x^3 / 6 1: 1 + x + x^2 / 2 + x^3 / 6
. .
' 1 + x + x^2/2! + x^3/3! RET Z F e myexp RET RET RET y
This polynomial is a Taylor series approximation to `exp(x)'.
The Z F command asks a number of questions. The above answers
say that the key sequence for our function should be z e; the
M-x equivalent should be calc-myexp; the name of the
function in algebraic formulas should also be myexp; the
default argument list `(x)' is acceptable; and finally y
answers the question "leave it in symbolic form for non-constant
arguments?"
1: 1.3495 2: 1.3495 3: 1.3495
. 1: 1.34986 2: 1.34986
. 1: myexp(a + 1)
.
.3 z e .3 E ' a+1 RET z e
First we call our new exp approximation with 0.3 as an
argument, and compare it with the true exp function. Then
we note that, as requested, if we try to give z e an
argument that isn't a plain number, it leaves the myexp
function call in symbolic form. If we had answered n to the
final question, `myexp(a + 1)' would have evaluated by plugging
in `a + 1' for `x' in the defining formula.
(*) Exercise 1. The "sine integral" function
Si(x) is defined as the integral of `sin(t)/t' for
t = 0 to x in radians. (It was invented because this
integral has no solution in terms of basic functions; if you give it
to Calc's a i command, it will ponder it for a long time and then
give up.) We can use the numerical integration command, however,
which in algebraic notation is written like `ninteg(f(t), t, 0, x)'
with any integrand `f(t)'. Define a z s command and
Si function that implement this. You will need to edit the
default argument list a bit. As a test, `Si(1)' should return
0.946083. (Hint: ninteg will run a lot faster if you reduce
the precision to, say, six digits beforehand.)
See section Programming Tutorial Exercise 1. (*)
The simplest way to do real "programming" of Emacs is to define a keyboard macro. A keyboard macro is simply a sequence of keystrokes which Emacs has stored away and can play back on demand. For example, if you find yourself typing H a S x RET often, you may wish to program a keyboard macro to type this for you.
1: y = sqrt(x) 1: x = y^2
. .
' y=sqrt(x) RET C-x ( H a S x RET C-x )
1: y = cos(x) 1: x = s1 arccos(y) + 2 pi n1
. .
' y=cos(x) RET X
When you type C-x (, Emacs begins recording. But it is also still ready to execute your keystrokes, so you're really "training" Emacs by walking it through the procedure once. When you type C-x ), the macro is recorded. You can now type X to re-execute the same keystrokes.
You can give a name to your macro by typing Z K.
1: . 1: y = x^4 1: x = s2 sqrt(s1 sqrt(y))
. .
Z K x RET ' y=x^4 RET z x
Notice that we use shift-Z to define the command, and lower-case z to call it up.
Keyboard macros can call other macros.
1: abs(x) 1: x = s1 y 1: 2 / x 1: x = 2 / y
. . . .
' abs(x) RET C-x ( ' y RET a = z x C-x ) ' 2/x RET X
(*) Exercise 2. Define a keyboard macro to negate the item in level 3 of the stack, without disturbing the rest of the stack. See section Programming Tutorial Exercise 2. (*)
(*) Exercise 3. Define keyboard macros to compute the following functions:
See section Programming Tutorial Exercise 3. (*)
(*) Exercise 4. Define a keyboard macro to compute the average (mean) value of a list of numbers. See section Programming Tutorial Exercise 4. (*)
In many programs, some of the steps must execute several times. Calc has looping commands that allow this. Loops are useful inside keyboard macros, but actually work at any time.
1: x^6 2: x^6 1: 360 x^2
. 1: 4 .
.
' x^6 RET 4 Z < a d x RET Z >
Here we have computed the fourth derivative of x^6 by enclosing a derivative command in a "repeat loop" structure. This structure pops a repeat count from the stack, then executes the body of the loop that many times.
If you make a mistake while entering the body of the loop, type Z C-g to cancel the loop command.
3: 1 2: 10946
2: 1 1: 17711
1: 20 .
.
1 RET RET 20 Z < TAB C-j + Z >
The numbers in levels 2 and 1 should be the 21st and 22nd Fibonacci numbers, respectively. (To see what's going on, try a few repetitions of the loop body by hand; C-j, also on the Line-Feed or LFD key if you have one, makes a copy of the number in level 2.)
A fascinating property of the Fibonacci numbers is that the nth
Fibonacci number can be found directly by computing @c{$\phi^n / \sqrt{5}$}
phi^n / sqrt(5)
and then rounding to the nearest integer, where @c{$\phi$ ("phi")}
phi, the
"golden ratio," is @c{$(1 + \sqrt{5}) / 2$}
(1 + sqrt(5)) / 2. (For convenience, this constant is available
from the phi variable, or the I H P command.)
1: 1.61803 1: 24476.0000409 1: 10945.9999817 1: 10946
. . . .
I H P 21 ^ 5 Q / R
(*) Exercise 5. The continued fraction representation of @c{$\phi$} phi is @c{$1 + 1/(1 + 1/(1 + 1/( \ldots )))$} 1 + 1/(1 + 1/(1 + 1/( ... ))). We can compute an approximate value by carrying this however far and then replacing the innermost @c{$1/( \ldots )$} 1/( ... ) by 1. Approximate phi using a twenty-term continued fraction. See section Programming Tutorial Exercise 5. (*)
(*) Exercise 6. Linear recurrences like the one for Fibonacci numbers can be expressed in terms of matrices. Given a vector [a, b] determine a matrix which, when multiplied by this vector, produces the vector [b, c], where a, b and c are three successive Fibonacci numbers. Now write a program that, given an integer n, computes the nth Fibonacci number using matrix arithmetic. See section Programming Tutorial Exercise 6. (*)
A more sophisticated kind of loop is the for loop. Suppose we wish to compute the 20th "harmonic" number, which is equal to the sum of the reciprocals of the integers from 1 to 20.
3: 0 1: 3.597739
2: 1 .
1: 20
.
0 RET 1 RET 20 Z ( & + 1 Z )
The "for" loop pops two numbers, the lower and upper limits, then repeats the body of the loop as an internal counter increases from the lower limit to the upper one. Just before executing the loop body, it pushes the current loop counter. When the loop body finishes, it pops the "step," i.e., the amount by which to increment the loop counter. As you can see, our loop always uses a step of one.
This harmonic number function uses the stack to hold the running total as well as for the various loop housekeeping functions. If you find this disorienting, you can sum in a variable instead:
1: 0 2: 1 . 1: 3.597739
. 1: 20 .
.
0 t 7 1 RET 20 Z ( & s + 7 1 Z ) r 7
The s + command adds the top-of-stack into the value in a variable (and removes that value from the stack).
It's worth noting that many jobs that call for a "for" loop can also be done more easily by Calc's high-level operations. Two other ways to compute harmonic numbers are to use vector mapping and reduction (v x 20, then V M &, then V R +), or to use the summation command a +. Both of these are probably easier than using loops. However, there are some situations where loops really are the way to go:
(*) Exercise 7. Use a "for" loop to find the first harmonic number which is greater than 4.0. See section Programming Tutorial Exercise 7. (*)
Of course, if we're going to be using variables in our programs, we have to worry about the programs clobbering values that the caller was keeping in those same variables. This is easy to fix, though:
. 1: 0.6667 1: 0.6667 3: 0.6667
. . 2: 3.597739
1: 0.6667
.
Z ` p 4 RET 2 RET 3 / s 7 s s a RET Z ' r 7 s r a RET
When we type Z ` (that's a back-quote character), Calc saves its mode settings and the contents of the ten "quick variables" for later reference. When we type Z ' (that's an apostrophe now), Calc restores those saved values. Thus the p 4 and s 7 commands have no effect outside this sequence. Wrapping this around the body of a keyboard macro ensures that it doesn't interfere with what the user of the macro was doing. Notice that the contents of the stack, and the values of named variables, survive past the Z ' command.
The Bernoulli numbers are a sequence with the interesting property that all of the odd Bernoulli numbers are zero, and the even ones, while difficult to compute, can be roughly approximated by the formula @c{$\displaystyle{2 n! \over (2 \pi)^n}$} 2 n! / (2 pi)^n. Let's write a keyboard macro to compute (approximate) Bernoulli numbers. (Calc has a command, k b, to compute exact Bernoulli numbers, but this command is very slow for large n since the higher Bernoulli numbers are very large fractions.)
1: 10 1: 0.0756823
. .
10 C-x ( RET 2 % Z [ DEL 0 Z : ' 2 $! / (2 pi)^$ RET = Z ] C-x )
You can read Z [ as "then," Z : as "else," and Z ] as "end-if." There is no need for an explicit "if" command. For the purposes of Z [, the condition is "true" if the value it pops from the stack is a nonzero number, or "false" if it pops zero or something that is not a number (like a formula). Here we take our integer argument modulo 2; this will be nonzero if we're asking for an odd Bernoulli number.
The actual tenth Bernoulli number is 5/66.
3: 0.0756823 1: 0 1: 0.25305 1: 0 1: 1.16659
2: 5:66 . . . .
1: 0.0757575
.
10 k b RET c f M-0 DEL 11 X DEL 12 X DEL 13 X DEL 14 X
Just to exercise loops a bit more, let's compute a table of even Bernoulli numbers.
3: [] 1: [0.10132, 0.03079, 0.02340, 0.033197, ...]
2: 2 .
1: 30
.
[ ] 2 RET 30 Z ( X | 2 Z )
The vertical-bar | is the vector-concatenation command. When we execute it, the list we are building will be in stack level 2 (initially this is an empty list), and the next Bernoulli number will be in level 1. The effect is to append the Bernoulli number onto the end of the list. (To create a table of exact fractional Bernoulli numbers, just replace X with k b in the above sequence of keystrokes.)
With loops and conditionals, you can program essentially anything in Calc. One other command that makes looping easier is Z /, which takes a condition from the stack and breaks out of the enclosing loop if the condition is true (non-zero). You can use this to make "while" and "until" style loops.
If you make a mistake when entering a keyboard macro, you can edit it using Z E. First, you must attach it to a key with Z K. One technique is to enter a throwaway dummy definition for the macro, then enter the real one in the edit command.
1: 3 1: 3 Keyboard Macro Editor.
. . Original keys: 1 RET 2 +
type "1\r"
type "2"
calc-plus
C-x ( 1 RET 2 + C-x ) Z K h RET Z E h
This shows the screen display assuming you have the `macedit' keyboard macro editing package installed, which is usually the case since a copy of `macedit' comes bundled with Calc.
A keyboard macro is stored as a pure keystroke sequence. The `macedit' package (invoked by Z E) scans along the macro and tries to decode it back into human-readable steps. If a key or keys are simply shorthand for some command with a M-x name, that name is shown. Anything that doesn't correspond to a M-x command is written as a `type' command.
Let's edit in a new definition, for computing harmonic numbers. First, erase the three lines of the old definition. Then, type in the new definition (or use Emacs M-w and C-y commands to copy it from this page of the Info file; you can skip typing the comments that begin with `#').
calc-kbd-push # Save local values (Z `) type "0" # Push a zero calc-store-into # Store it in variable 1 type "1" type "1" # Initial value for loop calc-roll-down # This is the TAB key; swap initial & final calc-kbd-for # Begin "for" loop... calc-inv # Take reciprocal calc-store-plus # Add to accumulator type "1" type "1" # Loop step is 1 calc-kbd-end-for # End "for" loop calc-recall # Now recall final accumulated value type "1" calc-kbd-pop # Restore values (Z ')
Press M-# M-# to finish editing and return to the Calculator.
1: 20 1: 3.597739
. .
20 z h
If you don't know how to write a particular command in `macedit'
format, you can always write it as keystrokes in a type command.
There is also a keys command which interprets the rest of the
line as standard Emacs keystroke names. In fact, `macedit' defines
a handy read-kbd-macro command which reads the current region
of the current buffer as a sequence of keystroke names, and defines that
sequence on the X (and C-x e) key. Because this is so
useful, Calc puts this command on the M-# m key. Try reading in
this macro in the following form: Press C-@ (or C-SPC) at
one end of the text below, then type M-# m at the other.
Z ` 0 t 1
1 TAB
Z ( & s + 1 1 Z )
r 1
Z '
(*) Exercise 8. A general algorithm for solving equations numerically is Newton's Method. Given the equation f(x) = 0 for any function f, and an initial guess x_0 which is reasonably close to the desired solution, apply this formula over and over:
where f'(x) is the derivative of f. The x
values will quickly converge to a solution, i.e., eventually
new_x and x will be equal to within the limits
of the current precision. Write a program which takes a formula
involving the variable x, and an initial guess x_0,
on the stack, and produces a value of x for which the formula
is zero. Use it to find a solution of @c{$\sin(\cos x) = 0.5$}
sin(cos(x)) = 0.5
near x = 4.5. (Use angles measured in radians.) Note that
the built-in a R (calc-find-root) command uses Newton's
method when it is able. See section Programming Tutorial Exercise 8. (*)
(*) Exercise 9. The digamma function @c{$\psi(z)$ ("psi")} psi(z) is defined as the derivative of @c{$\ln \Gamma(z)$} ln(gamma(z)). For large values of z, it can be approximated by the infinite sum
where @c{$\sum$}
sum represents the sum over n from 1 to infinity
(or to some limit high enough to give the desired accuracy), and
the bern function produces (exact) Bernoulli numbers.
While this sum is not guaranteed to converge, in practice it is safe.
An interesting mathematical constant is Euler's gamma, which is equal
to about 0.5772. One way to compute it is by the formula,
gamma = -psi(1). Unfortunately, 1 isn't a large enough argument
for the above formula to work (5 is a much safer value for z).
Fortunately, we can compute @c{$\psi(1)$}
psi(1) from @c{$\psi(5)$}
psi(5) using
the recurrence @c{$\psi(z+1) = \psi(z) + {1 \over z}$}
psi(z+1) = psi(z) + 1/z. Your task: Develop
a program to compute @c{$\psi(z)$}
psi(z); it should "pump up" z
if necessary to be greater than 5, then use the above summation
formula. Use looping commands to compute the sum. Use your function
to compute @c{$\gamma$}
gamma to twelve decimal places. (Calc has a built-in command
for Euler's constant, I P, which you can use to check your answer.)
See section Programming Tutorial Exercise 9. (*)
(*) Exercise 10. Given a polynomial in x and a number m on the stack, where the polynomial is of degree m or less (i.e., does not have any terms higher than x^m), write a program to convert the polynomial into a list-of-coefficients notation. For example, 5 x^4 + (x + 1)^2 with m = 6 should produce the list [1, 2, 1, 0, 5, 0, 0]. Also develop a way to convert from this form back to the standard algebraic form. See section Programming Tutorial Exercise 10. (*)
(*) Exercise 11. The Stirling numbers of the first kind are defined by the recurrences,
This can be implemented using a recursive program in Calc; the
program must invoke itself in order to calculate the two righthand
terms in the general formula. Since it always invokes itself with
"simpler" arguments, it's easy to see that it must eventually finish
the computation. Recursion is a little difficult with Emacs keyboard
macros since the macro is executed before its definition is complete.
So here's the recommended strategy: Create a "dummy macro" and assign
it to a key with, e.g., Z K s. Now enter the true definition,
using the z s command to call itself recursively, then assign it
to the same key with Z K s. Now the z s command will run
the complete recursive program. (Another way is to use Z E
or M-# m (read-kbd-macro) to read the whole macro at once,
thus avoiding the "training" phase.) The task: Write a program
that computes Stirling numbers of the first kind, given n and
m on the stack. Test it with small inputs like
s(4,2). (There is a built-in command for Stirling numbers,
k s, which you can use to check your answers.)
See section Programming Tutorial Exercise 11. (*)
The programming commands we've seen in this part of the tutorial are low-level, general-purpose operations. Often you will find that a higher-level function, such as vector mapping or rewrite rules, will do the job much more easily than a detailed, step-by-step program can:
(*) Exercise 12. Write another program for computing Stirling numbers of the first kind, this time using rewrite rules. Once again, n and m should be taken from the stack. See section Programming Tutorial Exercise 12. (*)
This ends the tutorial section of the Calc manual. Now you know enough about Calc to use it effectively for many kinds of calculations. But Calc has many features that were not even touched upon in this tutorial. The rest of this manual tells the whole story.
This section includes answers to all the exercises in the Calc tutorial.
1 RET 2 RET 3 RET 4 + * -
The result is @c{$1 - (2 \times (3 + 4)) = -13$} 1 - (2 * (3 + 4)) = -13.
2*4 + 7*9.5 + 5/4 = 75.75
After computing the intermediate term @c{$2\times4 = 8$} 2*4 = 8, you can leave that result on the stack while you compute the second term. With both of these results waiting on the stack you can then compute the final term, then press + + to add everything up.
2: 2 1: 8 3: 8 2: 8
1: 4 . 2: 7 1: 66.5
. 1: 9.5 .
.
2 RET 4 * 7 RET 9.5 *
4: 8 3: 8 2: 8 1: 75.75
3: 66.5 2: 66.5 1: 67.75 .
2: 5 1: 1.25 .
1: 4 .
.
5 RET 4 / + +
Alternatively, you could add the first two terms before going on with the third term.
2: 8 1: 74.5 3: 74.5 2: 74.5 1: 75.75
1: 66.5 . 2: 5 1: 1.25 .
. 1: 4 .
.
... + 5 RET 4 / +
On an old-style RPN calculator this second method would have the advantage of using only three stack levels. But since Calc's stack can grow arbitrarily large this isn't really an issue. Which method you choose is purely a matter of taste.
The TAB key provides a way to operate on the number in level 2.
3: 10 3: 10 4: 10 3: 10 3: 10
2: 20 2: 30 3: 30 2: 30 2: 21
1: 30 1: 20 2: 20 1: 21 1: 30
. . 1: 1 . .
.
TAB 1 + TAB
Similarly, M-TAB gives you access to the number in level 3.
3: 10 3: 21 3: 21 3: 30 3: 11
2: 21 2: 30 2: 30 2: 11 2: 21
1: 30 1: 10 1: 11 1: 21 1: 30
. . . . .
M-TAB 1 + M-TAB M-TAB
Either ( 2 , 3 ) or ( 2 SPC 3 ) would have worked, but using both the comma and the space at once yields:
1: ( ... 2: ( ... 1: (2, ... 2: (2, ... 2: (2, ...
. 1: 2 . 1: (2, ... 1: (2, 3)
. . .
( 2 , SPC 3 )
Joe probably tried to type TAB DEL to swap the extra incomplete object to the top of the stack and delete it. But a feature of Calc is that DEL on an incomplete object deletes just one component out of that object, so he had to press DEL twice to finish the job.
2: (2, ... 2: (2, 3) 2: (2, 3) 1: (2, 3)
1: (2, 3) 1: (2, ... 1: ( ... .
. . .
TAB DEL DEL
(As it turns out, deleting the second-to-top stack entry happens often enough that Calc provides a special key, M-DEL, to do just that. M-DEL is just like TAB DEL, except that it doesn't exhibit the "feature" that tripped poor Joe.)
Type ' sqrt($) RET.
If the Q key is broken, you could use ' $^0.5 RET. Or, RPN style, 0.5 ^.
(Actually, `$^1:2', using the fraction one-half as the power, is a closer equivalent, since `9^0.5' yields 3.0 whereas `sqrt(9)' and `9^1:2' yield the exact integer 3.)
In the formula `2 x (1+y)', `x' was interpreted as a function name with `1+y' as its argument. Assigning a value to a variable has no relation to a function by the same name. Joe needed to use an explicit `*' symbol here: `2 x*(1+y)'.
The result from 1 RET 0 / will be the formula 1 / 0. The "function" `/' cannot be evaluated when its second argument is zero, so it is left in symbolic form. When you now type 0 *, the result will be zero because Calc uses the general rule that "zero times anything is zero."
The m i command enables an infinite mode in which 1 / 0 results in a special symbol that represents "infinity." If you multiply infinity by zero, Calc uses another special new symbol to show that the answer is "indeterminate." See section Infinities, for further discussion of infinite and indeterminate values.
Calc always stores its numbers in decimal, so even though one-third has an exact base-3 representation (`3#0.1'), it is still stored as 0.3333333 (chopped off after 12 or however many decimal digits) inside the calculator's memory. When this inexact number is converted back to base 3 for display, it may still be slightly inexact. When we multiply this number by 3, we get 0.999999, also an inexact value.
When Calc displays a number in base 3, it has to decide how many digits
to show. If the current precision is 12 (decimal) digits, that corresponds
to `12 / log10(3) = 25.15' base-3 digits. Because 25.15 is not an
exact integer, Calc shows only 25 digits, with the result that stored
numbers carry a little bit of extra information that may not show up on
the screen. When Joe entered `3#0.2', the stored number 0.666666
happened to round to a pleasing value when it lost that last 0.15 of a
digit, but it was still inexact in Calc's memory. When he divided by 2,
he still got the dreaded inexact value 0.333333. (Actually, he divided
0.666667 by 2 to get 0.333334, which is why he got something a little
higher than 3#0.1 instead of a little lower.)
If Joe didn't want to be bothered with all this, he could have typed M-24 d n to display with one less digit than the default. (If you give d n a negative argument, it uses default-minus-that, so M-- d n would be an easier way to get the same effect.) Those inexact results would still be lurking there, but they would now be rounded to nice, natural-looking values for display purposes. (Remember, `0.022222' in base 3 is like `0.099999' in base 10; rounding off one digit will round the number up to `0.1'.) Depending on the nature of your work, this hiding of the inexactness may be a benefit or a danger. With the d n command, Calc gives you the choice.
Incidentally, another consequence of all this is that if you type M-30 d n to display more digits than are "really there," you'll see garbage digits at the end of the number. (In decimal display mode, with decimally-stored numbers, these garbage digits are always zero so they vanish and you don't notice them.) Because Calc rounds off that 0.15 digit, there is the danger that two numbers could be slightly different internally but still look the same. If you feel uneasy about this, set the d n precision to be a little higher than normal; you'll get ugly garbage digits, but you'll always be able to tell two distinct numbers apart.
An interesting side note is that most computers store their floating-point numbers in binary, and convert to decimal for display. Thus everyday programs have the same problem: Decimal 0.1 cannot be represented exactly in binary (try it: 0.1 d 2), so `0.1 * 10' comes out as an inexact approximation to 1 on some machines (though they generally arrange to hide it from you by rounding off one digit as we did above). Because Calc works in decimal instead of binary, you can be sure that numbers that look exact are exact as long as you stay in decimal display mode.
It's not hard to show that any number that can be represented exactly in binary, octal, or hexadecimal is also exact in decimal, so the kinds of problems we saw in this exercise are likely to be severe only when you use a relatively unusual radix like 3.
If the radix is 15 or higher, we can't use the letter `e' to mark the exponent because `e' is interpreted as a digit. When Calc needs to display scientific notation in a high radix, it writes `16#F.E8F*16.^15'. You can enter a number like this as an algebraic entry. Also, pressing e without any digits before it normally types 1e, but in a high radix it types 16.^ and puts you in algebraic entry: 16#f.e8f RET e 15 RET * is another way to enter this number.
The reason Calc puts a decimal point in the `16.^' is to prevent huge integers from being generated if the exponent is large (consider `16#1.23*16^1000', where we compute `16^1000' as a giant exact integer and then throw away most of the digits when we multiply it by the floating-point `16#1.23'). While this wouldn't normally matter for display purposes, it could give you a nasty surprise if you copied that number into a file and later moved it back into Calc.
The answer he got was 0.5000000000006399.
The problem is not that the square operation is inexact, but that the sine of 45 that was already on the stack was accurate to only 12 places. Arbitrary-precision calculations still only give answers as good as their inputs.
The real problem is that there is no 12-digit number which, when squared, comes out to 0.5 exactly. The f [ and f ] commands decrease or increase a number by one unit in the last place (according to the current precision). They are useful for determining facts like this.
1: 0.707106781187 1: 0.500000000001
. .
45 S 2 ^
1: 0.707106781187 1: 0.707106781186 1: 0.499999999999
. . .
U DEL f [ 2 ^
A high-precision calculation must be carried out in high precision all the way. The only number in the original problem which was known exactly was the quantity 45 degrees, so the precision must be raised before anything is done after the number 45 has been entered in order for the higher precision to be meaningful.
Many calculations involve real-world quantities, like the width and height of a piece of wood or the volume of a jar. Such quantities can't be measured exactly anyway, and if the data that is input to a calculation is inexact, doing exact arithmetic on it is a waste of time.
Fractions become unwieldy after too many calculations have been done with them. For example, the sum of the reciprocals of the integers from 1 to 10 is 7381:2520. The sum from 1 to 30 is 9304682830147:2329089562800. After a point it will take a long time to add even one more term to this sum, but a floating-point calculation of the sum will not have this problem.
Also, rational numbers cannot express the results of all calculations. There is no fractional form for the square root of two, so if you type 2 Q, Calc has no choice but to give you a floating-point answer.
Dividing two integers that are larger than the current precision may give a floating-point result that is inaccurate even when rounded down to an integer. Consider 123456789 / 2 when the current precision is 6 digits. The true answer is 61728394.5, but with a precision of 6 this will be rounded to @c{$12345700.0/2.0 = 61728500.0$} 12345700. / 2. = 61728500.. The result, when converted to an integer, will be off by 106.
Here are two solutions: Raise the precision enough that the floating-point round-off error is strictly to the right of the decimal point. Or, convert to fraction mode so that 123456789 / 2 produces the exact fraction 123456789:2, which can be rounded down by the F command without ever switching to floating-point format.
27 RET 9 B could give the exact result 3:2, but it does a floating-point calculation instead and produces 1.5.
Calc will find an exact result for a logarithm if the result is an integer or the reciprocal of an integer. But there is no efficient way to search the space of all possible rational numbers for an exact answer, so Calc doesn't try.
Duplicate the vector, compute its length, then divide the vector by its length: RET A /.
1: [1, 2, 3] 2: [1, 2, 3] 1: [0.27, 0.53, 0.80] 1: 1.
. 1: 3.74165738677 . .
.
r 1 RET A / A
The final A command shows that the normalized vector does indeed have unit length.
The average position is equal to the sum of the products of the positions times their corresponding probabilities. This is the definition of the dot product operation. So all you need to do is to put the two vectors on the stack and press *.
The trick is to multiply by a vector of ones. Use r 4 [1 1 1] * to get the row sum. Similarly, use [1 1] r 4 * to get the column sum.
Just enter the righthand side vector, then divide by the lefthand side matrix as usual.
1: [6, 10] 2: [6, 10] 1: [6 - 4 a / (b - a), 4 / (b - a) ]
. 1: [ [ 1, a ] .
[ 1, b ] ]
.
' [6 10] RET ' [1 a; 1 b] RET /
This can be made more readable using d B to enable "big" display mode:
4 a 4
1: [6 - -----, -----]
b - a b - a
Type d N to return to "normal" display mode afterwards.
To solve @c{$A^T A \, X = A^T B$} trn(A) * A * X = trn(A) * B, first we compute A2 = trn(A) * A and @c{$B' = A^T B$} B2 = trn(A) * B; now, we have a system @c{$A' X = B'$} A2 * X = B2 which we can solve using Calc's `/' command.
The first step is to enter the coefficient matrix. We'll store it in quick variable number 7 for later reference. Next, we compute the B2 vector.
1: [ [ 1, 2, 3 ] 2: [ [ 1, 4, 7, 2 ] 1: [57, 84, 96]
[ 4, 5, 6 ] [ 2, 5, 6, 4 ] .
[ 7, 6, 0 ] [ 3, 6, 0, 6 ] ]
[ 2, 4, 6 ] ] 1: [6, 2, 3, 11]
. .
' [1 2 3; 4 5 6; 7 6 0; 2 4 6] RET s 7 v t [6 2 3 11] *
Now we compute the matrix @c{$A'$} A2 and divide.
2: [57, 84, 96] 1: [-11.64, 14.08, -3.64]
1: [ [ 70, 72, 39 ] .
[ 72, 81, 60 ]
[ 39, 60, 81 ] ]
.
r 7 v t r 7 * /
(The actual computed answer will be slightly inexact due to round-off error.)
Notice that the answers are similar to those for the @c{$3\times3$} 3x3 system solved in the text. That's because the fourth equation that was added to the system is almost identical to the first one multiplied by two. (If it were identical, we would have gotten the exact same answer since the @c{$4\times3$} 4x3 system would be equivalent to the original @c{$3\times3$} 3x3 system.)
Since the first and fourth equations aren't quite equivalent, they can't both be satisfied at once. Let's plug our answers back into the original system of equations to see how well they match.
2: [-11.64, 14.08, -3.64] 1: [5.6, 2., 3., 11.2]
1: [ [ 1, 2, 3 ] .
[ 4, 5, 6 ]
[ 7, 6, 0 ]
[ 2, 4, 6 ] ]
.
r 7 TAB *
This is reasonably close to our original B vector, [6, 2, 3, 11].
We can use v x to build a vector of integers. This needs to be adjusted to get the range of integers we desire. Mapping `-' across the vector will accomplish this, although it turns out the plain `-' key will work just as well.
2: 2 2: 2
1: [1, 2, 3, 4, 5, 6, 7, 8, 9] 1: [-4, -3, -2, -1, 0, 1, 2, 3, 4]
. .
2 v x 9 RET 5 V M - or 5 -
Now we use V M ^ to map the exponentiation operator across the vector.
1: [0.0625, 0.125, 0.25, 0.5, 1, 2, 4, 8, 16]
.
V M ^
Given x and y vectors in quick variables 1 and 2 as before, the first job is to form the matrix that describes the problem.
Thus we want a @c{$19\times2$} 19x2 matrix with our x vector as one column and ones as the other column. So, first we build the column of ones, then we combine the two columns to form our A matrix.
2: [1.34, 1.41, 1.49, ... ] 1: [ [ 1.34, 1 ]
1: [1, 1, 1, ...] [ 1.41, 1 ]
. [ 1.49, 1 ]
...
r 1 1 v b 19 RET M-2 v p v t s 3
Now we compute @c{$A^T y$} trn(A) * y and @c{$A^T A$} trn(A) * A and divide.
1: [33.36554, 13.613] 2: [33.36554, 13.613]
. 1: [ [ 98.0003, 41.63 ]
[ 41.63, 19 ] ]
.
v t r 2 * r 3 v t r 3 *
(Hey, those numbers look familiar!)
1: [0.52141679, -0.425978]
.
/
Since we were solving equations of the form @c{$m \times x + b \times 1 = y$} m*x + b*1 = y, these numbers should be m and b, respectively. Sure enough, they agree exactly with the result computed using V M and V R!
The moral of this story: V M and V R will probably solve your problem, but there is often an easier way using the higher-level arithmetic functions!
In fact, there is a built-in a F command that does least-squares fits. See section Curve Fitting.
Move to one end of the list and press C-@ (or C-SPC or whatever) to set the mark, then move to the other end of the list and type M-# g.
1: [2.3, 6, 22, 15.1, 7, 15, 14, 7.5, 2.5]
.
To make things interesting, let's assume we don't know at a glance how many numbers are in this list. Then we could type:
2: [2.3, 6, 22, ... ] 2: [2.3, 6, 22, ... ]
1: [2.3, 6, 22, ... ] 1: 126356422.5
. .
RET V R *
2: 126356422.5 2: 126356422.5 1: 7.94652913734
1: [2.3, 6, 22, ... ] 1: 9 .
. .
TAB v l I ^
(The I ^ command computes the nth root of a number. You could also type & ^ to take the reciprocal of 9 and then raise the number to that power.)
A number j is a divisor of n if @c{$n \mathbin{\hbox{\code{\%}}} j = 0$} `n % j = 0'. The first step is to get a vector that identifies the divisors.
2: 30 2: [0, 0, 0, 2, ...] 1: [1, 1, 1, 0, ...]
1: [1, 2, 3, 4, ...] 1: 0 .
. .
30 RET v x 30 RET s 1 V M % 0 V M a = s 2
This vector has 1's marking divisors of 30 and 0's marking non-divisors.
The zeroth divisor function is just the total number of divisors. The first divisor function is the sum of the divisors.
1: 8 3: 8 2: 8 2: 8
2: [1, 2, 3, 4, ...] 1: [1, 2, 3, 0, ...] 1: 72
1: [1, 1, 1, 0, ...] . .
.
V R + r 1 r 2 V M * V R +
Once again, the last two steps just compute a dot product for which a simple * would have worked equally well.
The obvious first step is to obtain the list of factors with k f. This list will always be in sorted order, so if there are duplicates they will be right next to each other. A suitable method is to compare the list with a copy of itself shifted over by one.
1: [3, 7, 7, 7, 19] 2: [3, 7, 7, 7, 19] 2: [3, 7, 7, 7, 19, 0]
. 1: [3, 7, 7, 7, 19, 0] 1: [0, 3, 7, 7, 7, 19]
. .
19551 k f RET 0 | TAB 0 TAB |
1: [0, 0, 1, 1, 0, 0] 1: 2 1: 0
. . .
V M a = V R + 0 a =
Note that we have to arrange for both vectors to have the same length so that the mapping operation works; no prime factor will ever be zero, so adding zeros on the left and right is safe. From then on the job is pretty straightforward.
Incidentally, Calc provides the @c{\dfn{M\"obius} $\mu$} Moebius mu function which is zero if and only if its argument is square-free. It would be a much more convenient way to do the above test in practice.
First use v x 6 RET to get a list of integers, then V M v x to get a list of lists of integers!
Here's one solution. First, compute the triangular list from the previous exercise and type 1 - to subtract one from all the elements.
1: [ [0],
[0, 1],
[0, 1, 2],
...
1 -
The numbers down the lefthand edge of the list we desire are called the "triangular numbers" (now you know why!). The nth triangular number is the sum of the integers from 1 to n, and can be computed directly by the formula @c{$n (n+1) \over 2$} n * (n+1) / 2.
2: [ [0], [0, 1], ... ] 2: [ [0], [0, 1], ... ]
1: [0, 1, 2, 3, 4, 5] 1: [0, 1, 3, 6, 10, 15]
. .
v x 6 RET 1 - V M ' $ ($+1)/2 RET
Adding this list to the above list of lists produces the desired result:
1: [ [0],
[1, 2],
[3, 4, 5],
[6, 7, 8, 9],
[10, 11, 12, 13, 14],
[15, 16, 17, 18, 19, 20] ]
.
V M +
If we did not know the formula for triangular numbers, we could have computed them using a V U + command. We could also have gotten them the hard way by mapping a reduction across the original triangular list.
2: [ [0], [0, 1], ... ] 2: [ [0], [0, 1], ... ]
1: [ [0], [0, 1], ... ] 1: [0, 1, 3, 6, 10, 15]
. .
RET V M V R +
(This means "map a V R + command across the vector," and since each element of the main vector is itself a small vector, V R + computes the sum of its elements.)
The first step is to build a list of values of x.
1: [1, 2, 3, ..., 21] 1: [0, 1, 2, ..., 20] 1: [0, 0.25, 0.5, ..., 5]
. . .
v x 21 RET 1 - 4 / s 1
Next, we compute the Bessel function values.
1: [0., 0.124, 0.242, ..., -0.328]
.
V M ' besJ(1,$) RET
(Another way to do this would be 1 TAB V M f j.)
A way to isolate the maximum value is to compute the maximum using V R X, then compare all the Bessel values with that maximum.
2: [0., 0.124, 0.242, ... ] 1: [0, 0, 0, ... ] 2: [0, 0, 0, ... ]
1: 0.5801562 . 1: 1
. .
RET V R X V M a = RET V R + DEL
It's a good idea to verify, as in the last step above, that only one value is equal to the maximum. (After all, a plot of @c{$\sin x$} sin(x) might have many points all equal to the maximum value, 1.)
The vector we have now has a single 1 in the position that indicates the maximum value of x. Now it is a simple matter to convert this back into the corresponding value itself.
2: [0, 0, 0, ... ] 1: [0, 0., 0., ... ] 1: 1.75
1: [0, 0.25, 0.5, ... ] . .
.
r 1 V M * V R +
If a = had produced more than one 1 value, this method
would have given the sum of all maximum x values; not very
useful! In this case we could have used v m (calc-mask-vector)
instead. This command deletes all elements of a "data" vector that
correspond to zeros in a "mask" vector, leaving us with, in this
example, a vector of maximum x values.
The built-in a X command maximizes a function using more efficient methods. Just for illustration, let's use a X to maximize `besJ(1,x)' over this same interval.
2: besJ(1, x) 1: [1.84115, 0.581865]
1: [0 .. 5] .
.
' besJ(1,x), [0..5] RET a X x RET
The output from a X is a vector containing the value of x that maximizes the function, and the function's value at that maximum. As you can see, our simple search got quite close to the right answer.
Step one is to convert our integer into vector notation.
1: 25129925999 3: 25129925999
. 2: 10
1: [11, 10, 9, ..., 1, 0]
.
25129925999 RET 10 RET 12 RET v x 12 RET -
1: 25129925999 1: [0, 2, 25, 251, 2512, ... ]
2: [100000000000, ... ] .
.
V M ^ s 1 V M \
(Recall, the \ command computes an integer quotient.)
1: [0, 2, 5, 1, 2, 9, 9, 2, 5, 9, 9, 9]
.
10 V M % s 2
Next we must increment this number. This involves adding one to the last digit, plus handling carries. There is a carry to the left out of a digit if that digit is a nine and all the digits to the right of it are nines.
1: [0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 1] 1: [1, 1, 1, 0, 0, 1, ... ]
. .
9 V M a = v v
1: [1, 1, 1, 0, 0, 0, ... ] 1: [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1]
. .
V U * v v 1 |
Accumulating * across a vector of ones and zeros will preserve only the initial run of ones. These are the carries into all digits except the rightmost digit. Concatenating a one on the right takes care of aligning the carries properly, and also adding one to the rightmost digit.
2: [0, 0, 0, 0, ... ] 1: [0, 0, 2, 5, 1, 2, 9, 9, 2, 6, 0, 0, 0]
1: [0, 0, 2, 5, ... ] .
.
0 r 2 | V M + 10 V M %
Here we have concatenated 0 to the left of the original number; this takes care of shifting the carries by one with respect to the digits that generated them.
Finally, we must convert this list back into an integer.
3: [0, 0, 2, 5, ... ] 2: [0, 0, 2, 5, ... ]
2: 1000000000000 1: [1000000000000, 100000000000, ... ]
1: [100000000000, ... ] .
.
10 RET 12 ^ r 1 |
1: [0, 0, 20000000000, 5000000000, ... ] 1: 25129926000
. .
V M * V R +
Another way to do this final step would be to reduce the formula `10 $$ + $' across the vector of digits.
1: [0, 0, 2, 5, ... ] 1: 25129926000
. .
V R ' 10 $$ + $ RET
For the list [a, b, c, d], the result is ((a = b) = c) = d, which will compare a and b to produce a 1 or 0, which is then compared with c to produce another 1 or 0, which is then compared with d. This is not at all what Joe wanted.
Here's a more correct method:
1: [7, 7, 7, 8, 7] 2: [7, 7, 7, 8, 7]
. 1: 7
.
' [7,7,7,8,7] RET RET v r 1 RET
1: [1, 1, 1, 0, 1] 1: 0
. .
V M a = V R *
The circle of unit radius consists of those points (x,y) for which x^2 + y^2 < 1. We start by generating a vector of x^2 and a vector of y^2.
We can make this go a bit faster by using the v . and t . commands.
2: [2., 2., ..., 2.] 2: [2., 2., ..., 2.]
1: [2., 2., ..., 2.] 1: [1.16, 1.98, ..., 0.81]
. .
v . t . 2. v b 100 RET RET V M k r
2: [2., 2., ..., 2.] 1: [0.026, 0.96, ..., 0.036]
1: [0.026, 0.96, ..., 0.036] 2: [0.53, 0.81, ..., 0.094]
. .
1 - 2 V M ^ TAB V M k r 1 - 2 V M ^
Now we sum the x^2 and y^2 values, compare with 1 to get a vector of 1/0 truth values, then sum the truth values.
1: [0.56, 1.78, ..., 0.13] 1: [1, 0, ..., 1] 1: 84
. . .
+ 1 V M a < V R +
The ratio 84/100 should approximate the ratio @c{$\pi/4$} pi/4.
1: 0.84 1: 3.36 2: 3.36 1: 1.0695
. . 1: 3.14159 .
100 / 4 * P /
Our estimate, 3.36, is off by about 7%. We could get a better estimate by taking more points (say, 1000), but it's clear that this method is not very efficient!
(Naturally, since this example uses random numbers your own answer will be slightly different from the one shown here!)
If you typed v . and t . before, type them again to return to full-sized display of vectors.
This problem can be made a lot easier by taking advantage of some symmetries. First of all, after some thought it's clear that the y axis can be ignored altogether. Just pick a random x component for one end of the match, pick a random direction @c{$\theta$} theta, and see if x and @c{$x + \cos \theta$} x + cos(theta) (which is the x coordinate of the other endpoint) cross a line. The lines are at integer coordinates, so this happens when the two numbers surround an integer.
Since the two endpoints are equivalent, we may as well choose the leftmost of the two endpoints as x. Then theta is an angle pointing to the right, in the range -90 to 90 degrees. (We could use radians, but it would feel like cheating to refer to @c{$\pi/2$} pi/2 radians while trying to estimate @c{$\pi$} pi!)
In fact, since the field of lines is infinite we can choose the coordinates 0 and 1 for the lines on either side of the leftmost endpoint. The rightmost endpoint will be between 0 and 1 if the match does not cross a line, or between 1 and 2 if it does. So: Pick random x and @c{$\theta$} theta, compute @c{$x + \cos \theta$} x + cos(theta), and count how many of the results are greater than one. Simple!
We can make this go a bit faster by using the v . and t . commands.
1: [0.52, 0.71, ..., 0.72] 2: [0.52, 0.71, ..., 0.72]
. 1: [78.4, 64.5, ..., -42.9]
.
v . t . 1. v b 100 RET V M k r 180. v b 100 RET V M k r 90 -
(The next step may be slow, depending on the speed of your computer.)
2: [0.52, 0.71, ..., 0.72] 1: [0.72, 1.14, ..., 1.45]
1: [0.20, 0.43, ..., 0.73] .
.
m d V M C +
1: [0, 1, ..., 1] 1: 0.64 1: 3.125
. . .
1 V M a > V R + 100 / 2 TAB /
Let's try the third method, too. We'll use random integers up to one million. The k r command with an integer argument picks a random integer.
2: [1000000, 1000000, ..., 1000000] 2: [78489, 527587, ..., 814975]
1: [1000000, 1000000, ..., 1000000] 1: [324014, 358783, ..., 955450]
. .
1000000 v b 100 RET RET V M k r TAB V M k r
1: [1, 1, ..., 25] 1: [1, 1, ..., 0] 1: 0.56
. . .
V M k g 1 V M a = V R + 100 /
1: 10.714 1: 3.273
. .
6 TAB / Q
For a proof of this property of the GCD function, see section 4.5.2, exercise 10, of Knuth's Art of Computer Programming, volume II.
If you typed v . and t . before, type them again to return to full-sized display of vectors.
First, we put the string on the stack as a vector of ASCII codes.
1: [84, 101, 115, ..., 51]
.
"Testing, 1, 2, 3 RET
Note that the " key, like $, initiates algebraic entry so there was no need to type an apostrophe. Also, Calc didn't mind that we omitted the closing ". (The same goes for all closing delimiters like ) and ] at the end of a formula.
We'll show two different approaches here. In the first, we note that if the input vector is [a, b, c, d], then the hash code is 3 (3 (3a + b) + c) + d = 27a + 9b + 3c + d. In other words, it's a sum of descending powers of three times the ASCII codes.
2: [84, 101, 115, ..., 51] 2: [84, 101, 115, ..., 51]
1: 16 1: [15, 14, 13, ..., 0]
. .
RET v l v x 16 RET -
2: [84, 101, 115, ..., 51] 1: 1960915098 1: 121
1: [14348907, ..., 1] . .
.
3 TAB V M ^ * 511 %
Once again, * elegantly summarizes most of the computation. But there's an even more elegant approach: Reduce the formula 3 $$ + $ across the vector. Recall that this represents a function of two arguments that computes its first argument times three plus its second argument.
1: [84, 101, 115, ..., 51] 1: 1960915098
. .
"Testing, 1, 2, 3 RET V R ' 3$$+$ RET
If you did the decimal arithmetic exercise, this will be familiar. Basically, we're turning a base-3 vector of digits into an integer, except that our "digits" are much larger than real digits.
Instead of typing 511 % again to reduce the result, we can be cleverer still and notice that rather than computing a huge integer and taking the modulo at the end, we can take the modulo at each step without affecting the result. While this means there are more arithmetic operations, the numbers we operate on remain small so the operations are faster.
1: [84, 101, 115, ..., 51] 1: 121
. .
"Testing, 1, 2, 3 RET V R ' (3$$+$)%511 RET
Why does this work? Think about a two-step computation: 3 (3a + b) + c. Taking a result modulo 511 basically means subtracting off enough 511's to put the result in the desired range. So the result when we take the modulo after every step is,
for some suitable integers m and n. Expanding out by the distributive law yields
The m term in the latter formula is redundant because any contribution it makes could just as easily be made by the n term. So we can take it out to get an equivalent formula with n' = 3m + n,
which is just the formula for taking the modulo only at the end of the calculation. Therefore the two methods are essentially the same.
Later in the tutorial we will encounter modulo forms, which basically automate the idea of reducing every intermediate result modulo some value M.
We want to use H V U to nest a function which adds a random step to an (x,y) coordinate. The function is a bit long, but otherwise the problem is quite straightforward.
2: [0, 0] 1: [ [ 0, 0 ]
1: 50 [ 0.4288, -0.1695 ]
. [ -0.4787, -0.9027 ]
...
[0,0] 50 H V U ' <# + [random(2.0)-1, random(2.0)-1]> RET
Just as the text recommended, we used `< >' nameless function
notation to keep the two random calls from being evaluated
before nesting even begins.
We now have a vector of [x, y] sub-vectors, which by Calc's rules acts like a matrix. We can transpose this matrix and unpack to get a pair of vectors, x and y, suitable for graphing.
2: [ 0, 0.4288, -0.4787, ... ]
1: [ 0, -0.1696, -0.9027, ... ]
.
v t v u g f
Incidentally, because the x and y are completely independent in this case, we could have done two separate commands to create our x and y vectors of numbers directly.
To make a random walk of unit steps, we note that sincos of
a random direction exactly gives us an [x, y] step of unit
length; in fact, the new nesting function is even briefer, though
we might want to lower the precision a bit for it.
2: [0, 0] 1: [ [ 0, 0 ]
1: 50 [ 0.1318, 0.9912 ]
. [ -0.5965, 0.3061 ]
...
[0,0] 50 m d p 6 RET H V U ' <# + sincos(random(360.0))> RET
Another v t v u g f sequence will graph this new random walk.
An interesting twist on these random walk functions would be to use complex numbers instead of 2-vectors to represent points on the plane. In the first example, we'd use something like `random + random*(0,1)', and in the second we could use polar complex numbers with random phase angles. (This exercise was first suggested in this form by Randal Schwartz.)
If the number is the square root of @c{$\pi$} pi times a rational number, then its square, divided by @c{$\pi$} pi, should be a rational number.
1: 1.26508260337 1: 0.509433962268 1: 2486645810:4881193627
. . .
2 ^ P / c F
Technically speaking this is a rational number, but not one that is likely to have arisen in the original problem. More likely, it just happens to be the fraction which most closely represents some irrational number to within 12 digits.
But perhaps our result was not quite exact. Let's reduce the precision slightly and try again:
1: 0.509433962268 1: 27:53
. .
U p 10 RET c F
Aha! It's unlikely that an irrational number would equal a fraction this simple to within ten digits, so our original number was probably sqrt(27 pi / 53).
Notice that we didn't need to re-round the number when we reduced the precision. Remember, arithmetic operations always round their inputs to the current precision before they begin.
`inf / inf = nan'. Perhaps `1' is the "obvious" answer. But if `17 inf = inf', then `17 inf / inf = inf / inf = 17', too.
`exp(inf) = inf'. It's tempting to say that the exponential of infinity must be "bigger" than "regular" infinity, but as far as Calc is concerned all infinities are as just as big. In other words, as x goes to infinity, e^x also goes to infinity, but the fact the e^x grows much faster than x is not relevant here.
`exp(-inf) = 0'. Here we have a finite answer even though the input is infinite.
`sqrt(-inf) = (0, 1) inf'. Remember that (0, 1)
represents the imaginary number i. Here's a derivation:
`sqrt(-inf) = sqrt((-1) * inf) = sqrt(-1) * sqrt(inf)'.
The first part is, by definition, i; the second is inf
because, once again, all infinities are the same size.
`sqrt(uinf) = uinf'. In fact, we do know something about the
direction because sqrt is defined to return a value in the
right half of the complex plane. But Calc has no notation for this,
so it settles for the conservative answer uinf.
`abs(uinf) = inf'. No matter which direction x points, `abs(x)' always points along the positive real axis.
`ln(0) = -inf'. Here we have an infinite answer to a finite input. As in the 1 / 0 case, Calc will only use infinities here if you have turned on "infinite" mode. Otherwise, it will treat `ln(0)' as an error.
We can make `inf - inf' be any real number we like, say,
a, just by claiming that we added a to the first
infinity but not to the second. This is just as true for complex
values of a, so nan can stand for a complex number.
(And, similarly, uinf can stand for an infinity that points
in any direction in the complex plane, such as `(0, 1) inf').
In fact, we can multiply the first inf by two. Surely
`2 inf - inf = inf', but also `2 inf - inf = inf - inf = nan'.
So nan can even stand for infinity. Obviously it's just
as easy to make it stand for minus infinity as for plus infinity.
The moral of this story is that "infinity" is a slippery fish
indeed, and Calc tries to handle it by having a very simple model
for infinities (only the direction counts, not the "size"); but
Calc is careful to write nan any time this simple model is
unable to tell what the true answer is.
2: 0@ 47' 26" 1: 0@ 2' 47.411765"
1: 17 .
.
0@ 47' 26" RET 17 /
The average song length is two minutes and 47.4 seconds.
2: 0@ 2' 47.411765" 1: 0@ 3' 7.411765" 1: 0@ 53' 6.000005"
1: 0@ 0' 20" . .
.
20" + 17 *
The album would be 53 minutes and 6 seconds long.
Let's suppose it's January 14, 1991. The easiest thing to do is to keep trying 13ths of months until Calc reports a Friday. We can do this by manually entering dates, or by using t I:
1: <Wed Feb 13, 1991> 1: <Wed Mar 13, 1991> 1: <Sat Apr 13, 1991>
. . .
' <2/13> RET DEL ' <3/13> RET t I
(Calc assumes the current year if you don't say otherwise.)
This is getting tedious--we can keep advancing the date by typing t I over and over again, but let's automate the job by using vector mapping. The t I command actually takes a second "how-many-months" argument, which defaults to one. This argument is exactly what we want to map over:
2: <Sat Apr 13, 1991> 1: [<Mon May 13, 1991>, <Thu Jun 13, 1991>,
1: [1, 2, 3, 4, 5, 6] <Sat Jul 13, 1991>, <Tue Aug 13, 1991>,
. <Fri Sep 13, 1991>, <Sun Oct 13, 1991>]
.
v x 6 RET V M t I
1: 242
.
' <sep 13> - <jan 14> RET
And the answer to our original question: 242 days to go.
The full rule for leap years is that they occur in every year divisible by four, except that they don't occur in years divisible by 100, except that they do in years divisible by 400. We could work out the answer by carefully counting the years divisible by four and the exceptions, but there is a much simpler way that works even if we don't know the leap year rule.
Let's assume the present year is 1991. Years have 365 days, except that leap years (whenever they occur) have 366 days. So let's count the number of days between now and then, and compare that to the number of years times 365. The number of extra days we find must be equal to the number of leap years there were.
1: <Mon Jan 1, 10001> 2: <Mon Jan 1, 10001> 1: 2925593
. 1: <Tue Jan 1, 1991> .
.
' <jan 1 10001> RET ' <jan 1 1991> RET -
3: 2925593 2: 2925593 2: 2925593 1: 1943
2: 10001 1: 8010 1: 2923650 .
1: 1991 . .
.
10001 RET 1991 - 365 * -
There will be 1943 leap years before the year 10001. (Assuming, of course, that the algorithm for computing leap years remains unchanged for that long. See section Date Forms, for some interesting background information in that regard.)
The relative errors must be converted to absolute errors so that `+/-' notation may be used.
1: 1. 2: 1.
. 1: 0.2
.
20 RET .05 * 4 RET .05 *
Now we simply chug through the formula.
1: 19.7392088022 1: 394.78 +/- 19.739 1: 6316.5 +/- 706.21
. . .
2 P 2 ^ * 20 p 1 * 4 p .2 RET 2 ^ *
It turns out the v u command will unpack an error form as well as a vector. This saves us some retyping of numbers.
3: 6316.5 +/- 706.21 2: 6316.5 +/- 706.21
2: 6316.5 1: 0.1118
1: 706.21 .
.
RET v u TAB /
Thus the volume is 6316 cubic centimeters, within about 11 percent.
The first answer is pretty simple: `1 / (0 .. 10) = (0.1 .. inf)'. Since a number in the interval `(0 .. 10)' can get arbitrarily close to zero, its reciprocal can get arbitrarily large, so the answer is an interval that effectively means, "any number greater than 0.1" but with no upper bound.
The second answer, similarly, is `1 / (-10 .. 0) = (-inf .. -0.1)'.
Calc normally treats division by zero as an error, so that the formula `1 / 0' is left unsimplified. Our third problem, `1 / [0 .. 10]', also (potentially) divides by zero because zero is now a member of the interval. So Calc leaves this one unevaluated, too.
If you turn on "infinite" mode by pressing m i, you will instead get the answer `[0.1 .. inf]', which includes infinity as a possible value.
The fourth calculation, `1 / (-10 .. 10)', has the same problem. Zero is buried inside the interval, but it's still a possible value. It's not hard to see that the actual result of `1 / (-10 .. 10)' will be either greater than 0.1, or less than -0.1. Thus the interval goes from minus infinity to plus infinity, with a "hole" in it from -0.1 to 0.1. Calc doesn't have any way to represent this, so it just reports `[-inf .. inf]' as the answer. It may be disappointing to hear "the answer lies somewhere between minus infinity and plus infinity, inclusive," but that's the best that interval arithmetic can do in this case.
1: [-3 .. 3] 2: [-3 .. 3] 2: [0 .. 9]
. 1: [0 .. 9] 1: [-9 .. 9]
. .
[ 3 n .. 3 ] RET 2 ^ TAB RET *
In the first case the result says, "if a number is between -3 and 3, its square is between 0 and 9." The second case says, "the product of two numbers each between -3 and 3 is between -9 and 9."
An interval form is not a number; it is a symbol that can stand for many different numbers. Two identical-looking interval forms can stand for different numbers.
The same issue arises when you try to square an error form.
Testing the first number, we might arbitrarily choose 17 for x.
1: 17 mod 811749613 2: 17 mod 811749613 1: 533694123 mod 811749613
. 811749612 .
.
17 M 811749613 RET 811749612 ^
Since 533694123 is (considerably) different from 1, the number 811749613 must not be prime.
It's awkward to type the number in twice as we did above. There are various ways to avoid this, and algebraic entry is one. In fact, using a vector mapping operation we can perform several tests at once. Let's use this method to test the second number.
2: [17, 42, 100000] 1: [1 mod 15485863, 1 mod ... ]
1: 15485863 .
.
[17 42 100000] 15485863 RET V M ' ($$ mod $)^($-1) RET
The result is three ones (modulo n), so it's very probable that 15485863 is prime. (In fact, this number is the millionth prime.)
Note that the functions `($$^($-1)) mod $' or `$$^($-1) % $' would have been hopelessly inefficient, since they would have calculated the power using full integer arithmetic.
Calc has a k p command that does primality testing. For small numbers it does an exact test; for large numbers it uses a variant of the Fermat test we used here. You can use k p repeatedly to prove that a large integer is prime with any desired probability.
There are several ways to insert a calculated number into an HMS form. One way to convert a number of seconds to an HMS form is simply to multiply the number by an HMS form representing one second:
1: 31415926.5359 2: 31415926.5359 1: 8726@ 38' 46.5359"
. 1: 0@ 0' 1" .
.
P 1e7 * 0@ 0' 1" *
2: 8726@ 38' 46.5359" 1: 6@ 6' 2.5359" mod 24@ 0' 0"
1: 15@ 27' 16" mod 24@ 0' 0" .
.
x time RET +
It will be just after six in the morning.
The algebraic hms function can also be used to build an
HMS form:
1: hms(0, 0, 10000000. pi) 1: 8726@ 38' 46.5359"
. .
' hms(0, 0, 1e7 pi) RET =
The = key is necessary to evaluate the symbol `pi' to the actual number 3.14159...
As we recall, there are 17 songs of about 2 minutes and 47 seconds each.
2: 0@ 2' 47" 1: [0@ 3' 7" .. 0@ 3' 47"]
1: [0@ 0' 20" .. 0@ 1' 0"] .
.
[ 0@ 20" .. 0@ 1' ] +
1: [0@ 52' 59." .. 1@ 4' 19."]
.
17 *
No matter how long it is, the album will fit nicely on one CD.
Type ' 1 yr RET u c s RET. The answer is 31557600 seconds.
How long will it take for a signal to get from one end of the computer to the other?
1: m / c 1: 3.3356 ns
. .
' 1 m / c RET u c ns RET
(Recall, `c' is a "unit" corresponding to the speed of light.)
1: 3.3356 ns 1: 0.81356 ns / ns 1: 0.81356
2: 4.1 ns . .
.
' 4.1 ns RET / u s
Thus a signal could take up to 81 percent of a clock cycle just to go from one place to another inside the computer, assuming the signal could actually attain the full speed of light. Pretty tight!
The speed limit is 55 miles per hour on most highways. We want to find the ratio of Sam's speed to the US speed limit.
1: 55 mph 2: 55 mph 3: 11 hr mph / yd
. 1: 5 yd / hr .
.
' 55 mph RET ' 5 yd/hr RET /
The u s command cancels out these units to get a plain number. Now we take the logarithm base two to find the final answer, assuming that each successive pill doubles his speed.
1: 19360. 2: 19360. 1: 14.24
. 1: 2 .
.
u s 2 B
Thus Sam can take up to 14 pills without a worry.
The result `sqrt(x)^2' is simplified back to x by the Calculator, but `sqrt(x^2)' is not. (Consider what happens if x = -4.) If x is real, this formula could be simplified to `abs(x)', but for general complex arguments even that is not safe. (See section Declarations, for a way to tell Calc that x is known to be real.)
Suppose our roots are [a, b, c]. We want a polynomial which is zero when x is any of these values. The trivial polynomial x-a is zero when x=a, so the product (x-a)(x-b)(x-c) will do the job. We can use a c x to write this in a more familiar form.
1: 34 x - 24 x^3 1: [1.19023, -1.19023, 0]
. .
r 2 a P x RET
1: [x - 1.19023, x + 1.19023, x] 1: (x - 1.19023) (x + 1.19023) x
. .
V M ' x-$ RET V R *
1: x^3 - 1.41666 x 1: 34 x - 24 x^3
. .
a c x RET 24 n * a x
Sure enough, our answer (multiplied by a suitable constant) is the same as the original polynomial.
1: x sin(pi x) 1: (sin(pi x) - pi x cos(pi x)) / pi^2
. .
' x sin(pi x) RET m r a i x RET
1: [y, 1]
2: (sin(pi x) - pi x cos(pi x)) / pi^2
.
' [y,1] RET TAB
1: [(sin(pi y) - pi y cos(pi y)) / pi^2, (sin(pi) - pi cos(pi)) / pi^2]
.
V M $ RET
1: (sin(pi y) - pi y cos(pi y)) / pi^2 + (pi cos(pi) - sin(pi)) / pi^2
.
V R -
1: (sin(3.14159 y) - 3.14159 y cos(3.14159 y)) / 9.8696 - 0.3183
.
=
1: [0., -0.95493, 0.63662, -1.5915, 1.2732]
.
v x 5 RET TAB V M $ RET
The hard part is that V R + is no longer sufficient to add up all the contributions from the slices, since the slices have varying coefficients. So first we must come up with a vector of these coefficients. Here's one way:
2: -1 2: 3 1: [4, 2, ..., 4]
1: [1, 2, ..., 9] 1: [-1, 1, ..., -1] .
. .
1 n v x 9 RET V M ^ 3 TAB -
1: [4, 2, ..., 4, 1] 1: [1, 4, 2, ..., 4, 1]
. .
1 | 1 TAB |
Now we compute the function values. Note that for this method we need eleven values, including both endpoints of the desired interval.
2: [1, 4, 2, ..., 4, 1]
1: [1, 1.1, 1.2, ... , 1.8, 1.9, 2.]
.
11 RET 1 RET .1 RET C-u v x
2: [1, 4, 2, ..., 4, 1]
1: [0., 0.084941, 0.16993, ... ]
.
' sin(x) ln(x) RET m r p 5 RET V M $ RET
Once again this calls for V M * V R +; a simple * does the same thing.
1: 11.22 1: 1.122 1: 0.374
. . .
* .1 * 3 /
Wow! That's even better than the result from the Taylor series method.
We'll use Big mode to make the formulas more readable.
___
2 + V 2
1: (2 + sqrt(2)) / (1 + sqrt(2)) 1: --------
. ___
1 + V 2
.
' (2+sqrt(2)) / (1+sqrt(2)) RET d B
Multiplying by the conjugate helps because (a+b) (a-b) = a^2 - b^2.
___ ___
1: (2 + V 2 ) (V 2 - 1)
.
a r a/(b+c) := a*(b-c) / (b^2-c^2) RET
___ ___
1: 2 + V 2 - 2 1: V 2
. .
a r a*(b+c) := a*b + a*c a s
(We could have used a x instead of a rewrite rule for the second step.)
The multiply-by-conjugate rule turns out to be useful in many
different circumstances, such as when the denominator involves
sines and cosines or the imaginary constant i.
Here is the rule set:
[ fib(n) := fib(n, 1, 1) :: integer(n) :: n >= 1, fib(1, x, y) := x, fib(n, x, y) := fib(n-1, y, x+y) ]
The first rule turns a one-argument fib that people like to write
into a three-argument fib that makes computation easier. The
second rule converts back from three-argument form once the computation
is done. The third rule does the computation itself. It basically
says that if x and y are two consecutive Fibonacci numbers,
then y and x+y are the next (overlapping) pair of Fibonacci
numbers.
Notice that because the number n was "validated" by the conditions on the first rule, there is no need to put conditions on the other rules because the rule set would never get that far unless the input were valid. That further speeds computation, since no extra conditions need to be checked at every step.
Actually, a user with a nasty sense of humor could enter a bad
three-argument fib call directly, say, `fib(0, 1, 1)',
which would get the rules into an infinite loop. One thing that would
help keep this from happening by accident would be to use something like
`ZzFib' instead of fib as the name of the three-argument
function.
He got an infinite loop. First, Calc did as expected and rewrote `2 + 3 x' to `f(2, 3, x)'. Then it looked for ways to apply the rule again, and found that `f(2, 3, x)' looks like `a + b x' with `a = 0' and `b = 1', so it rewrote to `f(0, 1, f(2, 3, x))'. It then wrapped another `f(0, 1, ...)' around that, and so on, ad infinitum. Joe should have used M-1 a r to make sure the rule applied only once.
(Actually, even the first step didn't work as he expected. What Calc really gives for M-1 a r in this situation is `f(3 x, 1, 2)', treating 2 as the "variable," and `3 x' as a constant being added to it. While this may seem odd, it's just as valid a solution as the "obvious" one. One way to fix this would be to add the condition `:: variable(x)' to the rule, to make sure the thing that matches `x' is indeed a variable, or to change `x' to `quote(x)' on the lefthand side, so that the rule matches the actual variable `x' rather than letting `x' stand for something else.)
Here is a suitable set of rules to solve the first part of the problem:
[ seq(n, c) := seq(n/2, c+1) :: n%2 = 0, seq(n, c) := seq(3n+1, c+1) :: n%2 = 1 :: n > 1 ]
Given the initial formula `seq(6, 0)', application of these rules produces the following sequence of formulas:
seq( 3, 1) seq(10, 2) seq( 5, 3) seq(16, 4) seq( 8, 5) seq( 4, 6) seq( 2, 7) seq( 1, 8)
whereupon neither of the rules match, and rewriting stops.
We can pretty this up a bit with a couple more rules:
[ seq(n) := seq(n, 0), seq(1, c) := c, ... ]
Now, given `seq(6)' as the starting configuration, we get 8 as the result.
The change to return a vector is quite simple:
[ seq(n) := seq(n, []) :: integer(n) :: n > 0, seq(1, v) := v | 1, seq(n, v) := seq(n/2, v | n) :: n%2 = 0, seq(n, v) := seq(3n+1, v | n) :: n%2 = 1 ]
Given `seq(6)', the result is `[6, 3, 10, 5, 16, 8, 4, 2, 1]'.
Notice that the n > 1 guard is no longer necessary on the last rule since the n = 1 case is now detected by another rule. But a guard has been added to the initial rule to make sure the initial value is suitable before the computation begins.
While still a good idea, this guard is not as vitally important as it
was for the fib function, since calling, say, `seq(x, [])'
will not get into an infinite loop. Calc will not be able to prove
the symbol `x' is either even or odd, so none of the rules will
apply and the rewrites will stop right away.
If x is the sum a + b, then `nterms(x)' must be `nterms(a)' plus `nterms(b)'. If x is not a sum, then `nterms(x)' = 1.
[ nterms(a + b) := nterms(a) + nterms(b), nterms(x) := 1 ]
Here we have taken advantage of the fact that earlier rules always match before later rules; `nterms(x)' will only be tried if we already know that `x' is not a sum.
Just put the rule `0^0 := 1' into EvalRules. For example,
before making this definition we have:
2: [-2, -1, 0, 1, 2] 1: [1, 1, 0^0, 1, 1]
1: 0 .
.
v x 5 RET 3 - 0 V M ^
But then:
2: [-2, -1, 0, 1, 2] 1: [1, 1, 1, 1, 1]
1: 0 .
.
U ' 0^0:=1 RET s t EvalRules RET V M ^
Perhaps more surprisingly, this rule still works with infinite mode
turned on. Calc tries EvalRules before any built-in rules for
a function. This allows you to override the default behavior of any
Calc feature: Even though Calc now wants to evaluate 0^0 to
nan, your rule gets there first and evaluates it to 1 instead.
Just for kicks, try adding the rule 2+3 := 6 to EvalRules.
What happens? (Be sure to remove this rule afterward, or you might get
a nasty surprise when you use Calc to balance your checkbook!)
Here is a rule set that will do the job:
[ a*(b + c) := a*b + a*c,
opt(a) O(x^n) + opt(b) O(x^m) := O(x^n) :: n <= m
:: constant(a) :: constant(b),
opt(a) O(x^n) + opt(b) x^m := O(x^n) :: n <= m
:: constant(a) :: constant(b),
a O(x^n) := O(x^n) :: constant(a),
x^opt(m) O(x^n) := O(x^(n+m)),
O(x^n) O(x^m) := O(x^(n+m)) ]
If we really want the + and * keys to operate naturally
on power series, we should put these rules in EvalRules. For
testing purposes, it is better to put them in a different variable,
say, O, first.
The first rule just expands products of sums so that the rest of the
rules can assume they have an expanded-out polynomial to work with.
Note that this rule does not mention `O' at all, so it will
apply to any product-of-sum it encounters--this rule may surprise
you if you put it into EvalRules!
In the second rule, the sum of two O's is changed to the smaller O. The optional constant coefficients are there mostly so that `O(x^2) - O(x^3)' and `O(x^3) - O(x^2)' are handled as well as `O(x^2) + O(x^3)'.
The third rule absorbs higher powers of `x' into O's.
The fourth rule says that a constant times a negligible quantity is still negligible. (This rule will also match `O(x^3) / 4', with `a = 1/4'.)
The fifth rule rewrites, for example, `x^2 O(x^3)' to `O(x^5)'. (It is easy to see that if one of these forms is negligible, the other is, too.) Notice the `x^opt(m)' to pick up terms like `x O(x^3)'. Optional powers will match `x' as `x^1' but not 1 as `x^0'. This turns out to be exactly what we want here.
The sixth rule is the corresponding rule for products of two O's.
Another way to solve this problem would be to create a new "data type"
that represents truncated power series. We might represent these as
function calls `series(coefs, x)' where coefs is
a vector of coefficients for x^0, x^1, x^2, and so
on. Rules would exist for sums and products of such series
objects, and as an optional convenience could also know how to combine a
series object with a normal polynomial. (With this, and with a
rule that rewrites `O(x^n)' to the equivalent series form,
you could still enter power series in exactly the same notation as
before.) Operations on such objects would probably be more efficient,
although the objects would be a bit harder to read.
Some other symbolic math programs provide a power series data type
similar to this. Mathematica, for example, has an object that looks
like `PowerSeries[x, x0, coefs, nmin,
nmax, den]', where x0 is the point about which the
power series is taken (we've been assuming this was always zero),
and nmin, nmax, and den allow pseudo-power-series
with fractional or negative powers. Also, the PowerSeries
objects have a special display format that makes them look like
`2 x^2 + O(x^4)' when they are printed out. (See section Compositions,
for a way to do this in Calc, although for something as involved as
this it would probably be better to write the formatting routine
in Lisp.)
Just enter the formula `ninteg(sin(t)/t, t, 0, x)', type
Z F, and answer the questions. Since this formula contains two
variables, the default argument list will be `(t x)'. We want to
change this to `(x)' since t is really a dummy variable
to be used within ninteg.
The exact keystrokes are Z F s Si RET RET C-b C-b DEL DEL RET y. (The C-b C-b DEL DEL are what fix the argument list.)
One way is to move the number to the top of the stack, operate on it, then move it back: C-x ( M-TAB n M-TAB M-TAB C-x ).
Another way is to negate the top three stack entries, then negate again the top two stack entries: C-x ( M-3 n M-2 n C-x ).
Finally, it turns out that a negative prefix argument causes a command like n to operate on the specified stack entry only, which is just what we want: C-x ( M-- 3 n C-x ).
Just for kicks, let's also do it algebraically: C-x ( ' -$$$, $$, $ RET C-x ).
Each of these functions can be computed using the stack, or using algebraic entry, whichever way you prefer:
Computing @c{$\displaystyle{\sin x \over x}$} sin(x) / x:
Using the stack: C-x ( RET S TAB / C-x ).
Using algebraic entry: C-x ( ' sin($)/$ RET C-x ).
Computing the logarithm:
Using the stack: C-x ( TAB B C-x )
Using algebraic entry: C-x ( ' log($,$$) RET C-x ).
Computing the vector of integers:
Using the stack: C-x ( 1 RET 1 C-u v x C-x ). (Recall that C-u v x takes the vector size, starting value, and increment from the stack.)
Alternatively: C-x ( ~ v x C-x ). (The ~ key pops a number from the stack and uses it as the prefix argument for the next command.)
Using algebraic entry: C-x ( ' index($) RET C-x ).
Here's one way: C-x ( RET V R + TAB v l / C-x ).
2: 1 1: 1.61803398502 2: 1.61803398502
1: 20 . 1: 1.61803398875
. .
1 RET 20 Z < & 1 + Z > I H P
This answer is quite accurate.
Here is the matrix:
[ [ 0, 1 ] * [a, b] = [b, a + b] [ 1, 1 ] ]
Thus `[0, 1; 1, 1]^n * [1, 1]' computes Fibonacci numbers n+1 and n+2. Here's one program that does the job:
C-x ( ' [0, 1; 1, 1] ^ ($-1) * [1, 1] RET v u DEL C-x )
This program is quite efficient because Calc knows how to raise a matrix (or other value) to the power n in only @c{$\log_2 n$} log(n,2) steps. For example, this program can compute the 1000th Fibonacci number (a 209-digit integer!) in about 10 steps; even though the Z < ... Z > solution had much simpler steps, it would have required so many steps that it would not have been practical.
The trick here is to compute the harmonic numbers differently, so that the loop counter itself accumulates the sum of reciprocals. We use a separate variable to hold the integer counter.
1: 1 2: 1 1: .
. 1: 4
.
1 t 1 1 RET 4 Z ( t 2 r 1 1 + s 1 & Z )
The body of the loop goes as follows: First save the harmonic sum so far in variable 2. Then delete it from the stack; the for loop itself will take care of remembering it for us. Next, recall the count from variable 1, add one to it, and feed its reciprocal to the for loop to use as the step value. The for loop will increase the "loop counter" by that amount and keep going until the loop counter exceeds 4.
2: 31 3: 31
1: 3.99498713092 2: 3.99498713092
. 1: 4.02724519544
.
r 1 r 2 RET 31 & +
Thus we find that the 30th harmonic number is 3.99, and the 31st harmonic number is 4.02.
The first step is to compute the derivative f'(x) and thus the formula @c{$\displaystyle{x - {f(x) \over f'(x)}}$} x - f(x)/f'(x).
(Because this definition is long, it will be repeated in concise form below. You can use M-# m to load it from there. While you are entering a Z ` Z ' body in a macro, Calc simply collects keystrokes without executing them. In the following diagrams we'll pretend Calc actually executed the keystrokes as you typed them, just for purposes of illustration.)
2: sin(cos(x)) - 0.5 3: 4.5
1: 4.5 2: sin(cos(x)) - 0.5
. 1: -(sin(x) cos(cos(x)))
.
' sin(cos(x))-0.5 RET 4.5 m r C-x ( Z ` TAB RET a d x RET
2: 4.5
1: x + (sin(cos(x)) - 0.5) / sin(x) cos(cos(x))
.
/ ' x RET TAB - t 1
Now, we enter the loop. We'll use a repeat loop with a 20-repetition limit just in case the method fails to converge for some reason. (Normally, the Z / command will stop the loop before all 20 repetitions are done.)
1: 4.5 3: 4.5 2: 4.5
. 2: x + (sin(cos(x)) ... 1: 5.24196456928
1: 4.5 .
.
20 Z < RET r 1 TAB s l x RET
This is the new guess for x. Now we compare it with the old one to see if we've converged.
3: 5.24196 2: 5.24196 1: 5.24196 1: 5.26345856348
2: 5.24196 1: 0 . .
1: 4.5 .
.
RET M-TAB a = Z / Z > Z ' C-x )
The loop converges in just a few steps to this value. To check the result, we can simply substitute it back into the equation.
2: 5.26345856348
1: 0.499999999997
.
RET ' sin(cos($)) RET
Let's test the new definition again:
2: x^2 - 9 1: 3.
1: 1 .
.
' x^2-9 RET 1 X
Once again, here's the full Newton's Method definition:
C-x ( Z ` TAB RET a d x RET / ' x RET TAB - t 1
20 Z < RET r 1 TAB s l x RET
RET M-TAB a = Z /
Z >
Z '
C-x )
It turns out that Calc has a built-in command for applying a formula repeatedly until it converges to a number. See section Nesting and Fixed Points, to see how to use it.
Also, of course, a R is a built-in command that uses Newton's method (among others) to look for numerical solutions to any equation. See section Root Finding.
The first step is to adjust z to be greater than 5. A simple "for" loop will do the job here. If z is less than 5, we reduce the problem using @c{$\psi(z) = \psi(z+1) - 1/z$} psi(z) = psi(z+1) - 1/z. We go on to compute @c{$\psi(z+1)$} psi(z+1), and remember to add back a factor of -1/z when we're done. This step is repeated until z > 5.
(Because this definition is long, it will be repeated in concise form below. You can use M-# m to load it from there. While you are entering a Z ` Z ' body in a macro, Calc simply collects keystrokes without executing them. In the following diagrams we'll pretend Calc actually executed the keystrokes as you typed them, just for purposes of illustration.)
1: 1. 1: 1.
. .
1.0 RET C-x ( Z ` s 1 0 t 2
Here, variable 1 holds z and variable 2 holds the adjustment factor. If z < 5, we use a loop to increase it.
(By the way, we started with `1.0' instead of the integer 1 because otherwise the calculation below will try to do exact fractional arithmetic, and will never converge because fractions compare equal only if they are exactly equal, not just equal to within the current precision.)
3: 1. 2: 1. 1: 6.
2: 1. 1: 1 .
1: 5 .
.
RET 5 a < Z [ 5 Z ( & s + 2 1 s + 1 1 Z ) r 1 Z ]
Now we compute the initial part of the sum: @c{$\ln z - {1 \over 2z}$} ln(z) - 1/2z minus the adjustment factor.
2: 1.79175946923 2: 1.7084261359 1: -0.57490719743
1: 0.0833333333333 1: 2.28333333333 .
. .
L r 1 2 * & - r 2 -
Now we evaluate the series. We'll use another "for" loop counting up the value of 2 n. (Calc does have a summation command, a +, but we'll use loops just to get more practice with them.)
3: -0.5749 3: -0.5749 4: -0.5749 2: -0.5749
2: 2 2: 1:6 3: 1:6 1: 2.3148e-3
1: 40 1: 2 2: 2 .
. . 1: 36.
.
2 RET 40 Z ( RET k b TAB RET r 1 TAB ^ * /
3: -0.5749 3: -0.5772 2: -0.5772 1: -0.577215664892
2: -0.5749 2: -0.5772 1: 0 .
1: 2.3148e-3 1: -0.5749 .
. .
TAB RET M-TAB - RET M-TAB a = Z / 2 Z ) Z ' C-x )
This is the value of @c{$-\gamma$} - gamma, with a slight bit of roundoff error. To get a full 12 digits, let's use a higher precision:
2: -0.577215664892 2: -0.577215664892
1: 1. 1: -0.577215664901532
1. RET p 16 RET X
Here's the complete sequence of keystrokes:
C-x ( Z ` s 1 0 t 2
RET 5 a < Z [ 5 Z ( & s + 2 1 s + 1 1 Z ) r 1 Z ]
L r 1 2 * & - r 2 -
2 RET 40 Z ( RET k b TAB RET r 1 TAB ^ * /
TAB RET M-TAB - RET M-TAB a = Z /
2 Z )
Z '
C-x )
Taking the derivative of a term of the form x^n will produce a term like @c{$n x^{n-1}$} n x^(n-1). Taking the derivative of a constant produces zero. From this it is easy to see that the nth derivative of a polynomial, evaluated at x = 0, will equal the coefficient on the x^n term times n!.
(Because this definition is long, it will be repeated in concise form below. You can use M-# m to load it from there. While you are entering a Z ` Z ' body in a macro, Calc simply collects keystrokes without executing them. In the following diagrams we'll pretend Calc actually executed the keystrokes as you typed them, just for purposes of illustration.)
2: 5 x^4 + (x + 1)^2 3: 5 x^4 + (x + 1)^2
1: 6 2: 0
. 1: 6
.
' 5 x^4 + (x+1)^2 RET 6 C-x ( Z ` [ ] t 1 0 TAB
Variable 1 will accumulate the vector of coefficients.
2: 0 3: 0 2: 5 x^4 + ...
1: 5 x^4 + ... 2: 5 x^4 + ... 1: 1
. 1: 1 .
.
Z ( TAB RET 0 s l x RET M-TAB ! / s | 1
Note that s | 1 appends the top-of-stack value to the vector in a variable; it is completely analogous to s + 1. We could have written instead, r 1 TAB | t 1.
1: 20 x^3 + 2 x + 2 1: 0 1: [1, 2, 1, 0, 5, 0, 0]
. . .
a d x RET 1 Z ) DEL r 1 Z ' C-x )
To convert back, a simple method is just to map the coefficients against a table of powers of x.
2: [1, 2, 1, 0, 5, 0, 0] 2: [1, 2, 1, 0, 5, 0, 0]
1: 6 1: [0, 1, 2, 3, 4, 5, 6]
. .
6 RET 1 + 0 RET 1 C-u v x
2: [1, 2, 1, 0, 5, 0, 0] 2: 1 + 2 x + x^2 + 5 x^4
1: [1, x, x^2, x^3, ... ] .
.
' x RET TAB V M ^ *
Once again, here are the whole polynomial to/from vector programs:
C-x ( Z ` [ ] t 1 0 TAB
Z ( TAB RET 0 s l x RET M-TAB ! / s | 1
a d x RET
1 Z ) r 1
Z '
C-x )
C-x ( 1 + 0 RET 1 C-u v x ' x RET TAB V M ^ * C-x )
First we define a dummy program to go on the z s key. The true z s key is supposed to take two numbers from the stack and return one number, so DEL as a dummy definition will make sure the stack comes out right.
2: 4 1: 4 2: 4
1: 2 . 1: 2
. .
4 RET 2 C-x ( DEL C-x ) Z K s RET 2
The last step replaces the 2 that was eaten during the creation of the dummy z s command. Now we move on to the real definition. The recurrence needs to be rewritten slightly, to the form s(n,m) = s(n-1,m-1) - (n-1) s(n-1,m).
(Because this definition is long, it will be repeated in concise form below. You can use M-# m to load it from there.)
2: 4 4: 4 3: 4 2: 4
1: 2 3: 2 2: 2 1: 2
. 2: 4 1: 0 .
1: 2 .
.
C-x ( M-2 RET a = Z [ DEL DEL 1 Z :
4: 4 2: 4 2: 3 4: 3 4: 3 3: 3
3: 2 1: 2 1: 2 3: 2 3: 2 2: 2
2: 2 . . 2: 3 2: 3 1: 3
1: 0 1: 2 1: 1 .
. . .
RET 0 a = Z [ DEL DEL 0 Z : TAB 1 - TAB M-2 RET 1 - z s
(Note that the value 3 that our dummy z s produces is not correct; it is merely a placeholder that will do just as well for now.)
3: 3 4: 3 3: 3 2: 3 1: -6
2: 3 3: 3 2: 3 1: 9 .
1: 2 2: 3 1: 3 .
. 1: 2 .
.
M-TAB M-TAB TAB RET M-TAB z s * -
1: -6 2: 4 1: 11 2: 11
. 1: 2 . 1: 11
. .
Z ] Z ] C-x ) Z K s RET DEL 4 RET 2 z s M-RET k s
Even though the result that we got during the definition was highly bogus, once the definition is complete the z s command gets the right answers.
Here's the full program once again:
C-x ( M-2 RET a =
Z [ DEL DEL 1
Z : RET 0 a =
Z [ DEL DEL 0
Z : TAB 1 - TAB M-2 RET 1 - z s
M-TAB M-TAB TAB RET M-TAB z s * -
Z ]
Z ]
C-x )
You can read this definition using M-# m (read-kbd-macro)
followed by Z K s, without having to make a dummy definition
first, because read-kbd-macro doesn't need to execute the
definition as it reads it in. For this reason, M-# m is often
the easiest way to create recursive programs in Calc.
This turns out to be a much easier way to solve the problem. Let's denote Stirling numbers as calls of the function `s'.
First, we store the rewrite rules corresponding to the definition of Stirling numbers in a convenient variable:
s e StirlingRules RET [ s(n,n) := 1 :: n >= 0, s(n,0) := 0 :: n > 0, s(n,m) := s(n-1,m-1) - (n-1) s(n-1,m) :: n >= m :: m >= 1 ] C-c C-c
Now, it's just a matter of applying the rules:
2: 4 1: s(4, 2) 1: 11
1: 2 . .
.
4 RET 2 C-x ( ' s($$,$) RET a r StirlingRules RET C-x )
As in the case of the fib rules, it would be useful to put these
rules in EvalRules and to add a `:: remember' condition to
the last rule.