MACHINE CODE CALCULATOR
part 5 of 5
by Toni Baker
from ZX Computing, November 1986

Toni Baker rounds off the series with a look
at the function generator.


There is just one calculator instruction left to cover. It is the
function "series" otherwise known as the function generator or series
generator. Its code will be a value between 80h and 9F. The last five
bits of the code form a parameter, so that 86h means "series 6", 8C
means "series twelve", 99 means "series twenty-five" and so forth.
This is the single most powerful instruction in the whole of the
calculator set. It is the function with which SIN and EXP and others
were written in the ROM. With it we may create our own designer
functions, or implement mathematical functions which are not present
(and not otherwise possible) on the Spectrum. The "series" function is
the calculator's final frontier.

Essentially, what the "series" instruction does is to evaluate a
polynomial expression. "Polynomial" is simply a highbrow piece of
mathematical jargon - it means an expression which looks something
like this ...

a1 + a2X + a3X^2 + ... anX^(n-1)

I apologise to those of you not acquainted with mathematical notation.
The same expression written in Spectrum BASIC would be as follows ...

A(1) + A(2)*X + A(3)X^2 + ... + A(N)*X^(N-1)

The dots in the middle of the above simply mean "and so on until". If
you know the value of N you can fill in the missing chunk and put the
whole lot into a single BASIC expression. If you don't know the value
of N then you would have to use a FOR/NEXT loop from 1 to N in order
to work it all out. But what does it all mean? Let's take it apart and
find out.

Firstly the number 'n' - this is called the "degree" of the polynomial
(technically the degree is n-1, not n). The "series" instruction
specifies the number 'n' explicitly, so that, for example, "series 6"
tells us that n equals six; "series 8" would mean that n equals eight.
The parameter, which is part of the hex code (bits 4 to 0) actually
specifies the degree of the polynomial. You should note that since you
only have five bits in which to specify this parameter, its maximum
value (other than zero) is hex 1F, or thirty-one. Note also that if
you specify n to be zero, the ROM will mistakenly interpret this as
two hundred and fifty six.

A polynomial expression such as we have been discussing consists of a
number of terms, each separated by a "+" sign. The first term is a1,
and the last term is ai x^(i-1). In general, the i-th term will be ai
xi-1. The value of the whole polynomial is therefore the sum of all of
its terms. There are n such terms, and each of them contains the
variable 'x' (except the first one - this is because in the first term
x would have to be raised to the pawer of zero, and anything to the
power of zero is one). Because every term except the first one
contains the variable 'x' (and the first term is simply a number to be
added), it follows that the whole thing is simply a function of x. In
other words - you put in a value for x at one end, and you get a new
number out at the other. This is in common with all of the other
functions of the Spectrum.

The 'x' in this case is the value at the top of the calculator stack.
When the "series" instruction is encountered, the value 'x' is removed
from the calculator stack, the expression is evaluated, and the result
put back onto the calculator stack in place of the original 'x'.

All we need to know now are the value of the 'a's- You can think of
the 'a's as being a BASIC array A() dimensioned up to n, so that the
first term is simply A(1), and the last term is A(N) multiplied by
x^(N-1).

Before we look at how to specify the value of the 'A's, let's have a
look at how we may use such a polynomial to calculate functions.


Polynomials

Suppose that the value of 'x' is somewhere between minus one and plus
one (ie. that x is greater than minus one, and less than one) - we can
of course ensure this since we already know how to manipulate the
calculator stack. If x is zero then the value of the polynomial will
simply be a, since all of the other terms will be multiplied by zero.
If x is non-zero, but still between minus and plus one, then ABS(x^2)
will be less than ABS(x): ABS(x^3) will be less than ABS(x^2), and so
on. it follows, therefore, that if the 'A's are all roughly the same
size, then each term will be smaller than the last. Furthermore, if
the 'A's themselves also get smaller (ie. if ABS(A(2)) is less than
ABS(A(1)): ABS(A(3)) is less than ABS(A(2)); and so on) then each term
will be smaller still than the last. Indeed it is possible to ensure
that the last term is so small, by comparison to the first term, that
if any more terms were added it would be outside the limits of the
computer's accuracy, and would thus make no more difference. For
instance, if A(1) equals 1, and if the last term equals 2^-64 or less,
then the Spectrum would have to round its answer up by simply ignoring
the last term.

It is in theory possible to simulate any function whatsoever, provided
that x is between minus one and plus one, that each term is suitably
smaller than the last (although I can't go into the precise
mathematical details of this condition in an article such as this),
and that there are an infinite number of terms.

You see, with an infinite number of terms, we can make the polynomial
closer and closer to the desired function with each new term.
Unfortunately for us, we are only allowed a maximum of thirty one
terms altogether. We can get round this problem by ensuring that the
polynomial is close enough - it doesn't have to be exactly right. As
long as the maximum error is smaller than the accuracy possible on the
Spectrum, then the polynomial will calculate the function - at least
to the limits of Spectrum accuracy.

So how do we work out the 'A's? - Well, mathematicians may care to use
Taylor's or Maclaurin's Theorem (which I can't go into here), and
everyone else will have to look the values up in books. Almost any
A-level maths book (available from your local library) will tell you
what the correct series is to simulate SIN, or EXP, or whatever.

For instance, take the function 2^X. If X is between minus one and
plus one then the function can be simulated by the following
polynomial. 2^X is approximately equal to 1 + 0.69314718*X +
0.24022651*X^2 + 0.055504109*X^3 + 0.0096181291*X^4 + 0.0013333558*X^5
+ 0.0001540353*X^6 + 0.000015252734*X^7 ... (Incidentally -I didn't
work out the above numbers - I got them out of a book. I advise you to
do the same. Mathematicians are, of course, welcome to work as much
out for themselves as they want.)

How does this help us? Well - you see we are still restricted to only
using values of X between minus one and plus one (note: if you want to
use negative numbers on the Spectrum for the above formula you'll have
to use "X*X*X*X*X" instead of "X^5", and so on). This is no good - the
Spectrum can calculate EXP(X) for all values of X, not just small
ones. To find the way round this problem we'll need to do a bit of
maths. It's more difficult, but not, I hope, too difficult.

Firstly, note that EXP(X) is defined as e^x, where e = 2.7182818

LET            u =x / (LN 2)
therefore:     x =u * (LN 2)
therefore:   e^x =e^(u*(LN 2))
                 =(e^(LN 2))^u
                 =2^u

Now, since u is a floating point number, it must have an integer part,
and a fractional part. Thus:

LET            i =INT u
LET            f =u - i
therefore:     u =i + f
thus we      e^x =2^u
have             =2^(i+f)
                 =(2^i) * (2^f)

Now, since f is a fraction, between zero and one, we can use a
polynomial to simulate the function. The polynomial I gave earlier may
be used to calculate 2^f. All we need to do now is to multiply the
result by 2^i (where i s an integer). We can do this simply by adding
i to the exponent byte of the five-byte form of the number. This is
the procedure used by the ROM to calculate EXP for all values of x.

Well, I won't bore you to tears with any more maths, but I hope you
can see how ingenious little tricks like the above can be used to
ensure that polynomial approximations are only used on numbers between
minus one and plus one. Such ingenious little tricks are necessary,
because the polynomial formula won't work with numbers outside this
range (for instance - the value of SIN(x) is always between -1 and +1,
and yet, even if we used 256 terms, there would always be some large
value of x for which a256 x^256 was greater than one. Clearly this
would be the wrong answer. You must devise a means of ensuring that
polynomial approximations to functions are always fed with numbers in
range -1 to +1).


Instructions

Now - at last - we need to turn our attention to exactly how the
values for the 'A's are passed to the "series" instruction. They are
not specified explicitly - that would be too simple! Instead, we have
to work out a new array, which I'll call B(), also dimensioned up to
N. It is the values of the 'B's which are passed to the Spectrum.
Therefore, our first task, knowing the value of the 'A's, is to work
out the values of the 'B's. The BASIC program of Figure One
[MCCALC.TAP "part 5"] will do this task for us. Don't worry too much
about the algorithm used, just take my word for it - it works!. (The
mathematics required to prove that it works is beyond the scope of
this article).

The program of Figure One doesn't print the 'B's in decimal - instead
it prints them in the form required by the series instruction,
converting them first to five byte form, and then to the compressed
form used by the "stk_data" instruction. This form is also used by the
"series" instruction.

In other words, if you input the 'A's - which technically are the
coefficients of a power series (don't worry if you don't understand
the terms), then the program will miraculously transform the numbers
into a new set of numbers - those used by the "series" instruction.
The program leaves no work to you at all. The output of the program is
the complete series instruction (in hex) required to evaluate the
required polynomial, including the initial "series" byte. For extra
clarity, you could, if you so desired, add an extra BASIC line at 1285
PRINT,B(i) which would show you the 'B's in decimal as well as in hex.

If you now feed the bytes printed on the screen into a calculator
program, then hey presto - your function will be operational!

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Listing 1

1010 LET k=0
1020 INPUT "N = ";n
1030 PRINT "N = ";n
1040 DIM a(n)
1050 FOR i=1 TO n
1060 INPUT "A(";(i);") = ";a(i)
1070 PRINT "A(";(i);") = ";a(i)
1080 NEXT i
1090 GO SUB 1340
1100 PRINT 
1110 PRINT FN h$(n+128)
1120 FOR i=1 TO n
1130 LET k=b(i)
1140 LET q=PEEK 23627+256*PEEK 23628
1150 FOR c=3 TO 1 STEP -1
1160 IF PEEK (q+c+2)=0 THEN NEXT c
1170 LET d=PEEK (q+1)
1180 IF d<81 OR d>143 THEN GO TO 1210
1190 LET d=d-80+64*c
1200 GO TO 1240
1210 LET d=d-80
1220 IF d<0 THEN LET d=d+256
1230 PRINT FN h$(64*c);
1240 PRINT FN h$(d);
1250 FOR j=0 TO c
1260 PRINT FN h$(PEEK (q+j+2));
1270 NEXT j
1280 PRINT 
1290 NEXT i
1300 STOP 
1310 DEF FN h$(a)=FN k$(INT (a/16))+FN k$(a-16*INT (a/16))+" "
1320 DEF FN k$(a)="0123456789ABCDEF"(a+1)
1330 REM CALCULATES VALUES               SUBROUTINE
1340 DIM b(n)
1350 FOR a=1 TO n
1360 FOR b=1+n-a TO n STEP 2
1370 LET x=b-1
1380 GO SUB 1520
1390 LET j=y
1400 LET x=(a+b-n-1)/2
1410 GO SUB 1520
1420 LET j=j/y
1430 LET x=(b-a+n-1)/2
1440 GO SUB 1520
1450 LET j=j/y
1460 LET j=j*2^(1-b)
1470 LET b(a)=b(a)+j*a(b)
1480 NEXT b
1490 NEXT a
1500 RETURN 
1510 REM LET Y=X FACTORIAL
1520 LET y=1
1530 FOR i=1 TO x
1540 LET y=y*i
1550 NEXT i
1560 RETURN 
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


Gamma

We shall demonstrate this procedure now by introducing a brand new
function to Spectrum BASIC. It will be referred to in BASIC as FN G(),
but the function will be defined mathematically as the GAMMA FUNCTION.
It is very similar to the FACTORIAL function described last month, in
that GAMMA(x) equals FACTORIAL(x-1) for all positive integer values of
x, but whereas FACTORIAL only works for integers, GAMMA works for all
numbers. A graph drawn of Y = GAMMA(x) will form a smooth continuous
curve, and for all values of x it is true that GAMMA(x+1) equals
x*GAMMA(x).

Furthermore, I shall go through, in detail, the process by which the
function is being created. First of all we need a polynomial. I came
across such a polynomial approximation in the book "Mathematical
Methods in the Physical Sciences" by Merle C Potter. The polynomial is
only valid if x is between 0 and 1, but is more accurate if x is near
zero. For this reason I have chosen only to use values of x between 0
and 1/2. The polynomial calculates the value of GAMMA(x+1) if x is in
this range. The polynomial is as follows:

GAMMA(x+1) = 1 - 0.577191652x
+ 0.988205891x^2 - 0.897056937x^3
+ 0.918206857x^4 - 0.756704078x^5
+ 0.482199394x^6 - 0.193527818x^7
+ 0.035868343x^8

Now, as I have said, I want to make sure that the polynomial is only
used for values of x between 0 and 1/2, and yet I want the GAMMA
function, FN G(), to work for all values of x, so we must find a way
around this problem.

Firstly, we can make use of the rule GAMMA(x+1) equals x*GAMMA(x). We
can use the rule repeatedly (in a loop) to reach any value of x, so
long as the fractional part (the part to the right of the decimal
point) is no greater than a half. For these remaining numbers we can
use a different mathematical rule; for all x, GAMMA(x)*GAMMA(1-x) =
PI/SIN(x*PI)

The first step, then, is to convert the polynomial into a "series"
instruction. To do this we simply run the BASIC program of Figure One,
and input 9 (because there are nine terms), and then input the nine
numbers (1, -0.577191652, and so on). The program will print out the
required form of the "series 9" instructions, which you will find
incorporated into the main machine code program accompanying this
article. The rest of it is just a hard slog, isolating the various
cases, and dealing with them accordingly.

Incidentally, although I have used a loop to implement the GAMMA(x+1)
= x*GAMMA(x) rule, this is not the only way of doing things. Another
way would have been to call the GAMMA subroutine recursively (ie. from
within itself). I decided not to use this method because it would use
up more memory on the calculator stack (an additional five bytes for
each recursive call), but if you feel that that would not be too much
of a disadvantage you might like to re-write the routine using this
method.

Finally, the proof of the pudding is in the eating. The final test is
whether or not it works. One way we can check is to calculate
GAMMA(1/2) - the answer should come out to be SQR(PI). The best way to
check, however, is to plot a graph. The function should produce a
smooth continuous curve. The BASIC program of Figure Two [MCCALC.TAP
"part 5"] will plot such a graph, and I will leave it to you to draw
your own conclusions from it.

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Listing 2

2010 CLS 
2020 FOR i=1 TO 255
2030 LET x=i/12.5
2040 LET y=LN FN g(x)
2050 LET j=y*4+10
2060 IF j<176 THEN PLOT i,j
2070 NEXT i
2099 STOP 
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

Well that's it from me for this series. We have now covered everything
the calculator can do (or have we?). We have covered every single
calculator instruction possible (or have we?). This now concludes our
series. I hope you have enjoyed reading it. Good programming, and may
the force be with you.

(But - just a thought - what would happen if you used "invalid"
calculator codes? - codes between 3E and 7F? With code 43h for
instance, control could be directed into RAM, and enable the use of
user-defined calculator instructions - defined in machine code of
course. Interesting? Maybe one day the saga may continue ... But not
for the moment).


- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Gamma Function
BASIC
-----
DEF FN g(x)=USR gamma_fn
The number following the keyword USR should be the address
of the label GAMMA_FN in the machine code program.

MACHINE CODE
------------
EA60 EF       GAMMA     RST  #28           x
EA61 C5                 store M5           (M5 contains x)
EA62 31                 duplicate          x,x
     31                 duplicate          x,x,x
     27                 int                x,x,INT x
EA65 03                 subtract           x,x-INT x
EA66 31                 duplicate          x,x-INT x,x-INT x
     A2                 const half         x,x-INT x,x-INT x,1/2
     03                 subtract           x,x-INT x,x-INT x-1/2
EA69 37                 gt zero            x,x-INT x,x-INT x>1/2?
EA6A 0066               jump true,G_HALF   x,x-INT x
EA6C 89                 series 9
EA6D E412EAAF03         .00014011071
EA72 E7C62C2805         -.0015119361
EA77 EA0DCEC143         .0086552512
EA7C EC8C3556A3         -.034230555
EA81 ED5A25B37E         0.10651722
EA86 EF86344D75         -0.26211779
EA8B F018F3752D         0.59746487
EA90 F0EA14CE2B         -0.91437996
EA95 F17FDCCB17         1.9989256

                                           x,GAMMA(x+1-INT x)
EA9A 01                 exchange           GAMMA(x+1-INT x),x
     27                 int                GAMMA(x+1-INT x),INT x
     A1                 const one          GAMMA(x+1-INT x),INT x,1
EA9D 03                 subtract           GAMMA(x+1-INT x),INT x-1
EA9E 31                 duplicate          GAMMA(x+1-INT x),INT x-1,INT x-1
     30                 eq zero            GAMMA(x+1-INT x),INT x-1,INT x=1?
     0017               jump true,G_1_TO_2 Jump if INT x = 1;
                                           ie. if x is between one and two

                                           GAMMA(x+1-INT x),INT x-1
EAA2 31                 duplicate          GAMMA(x+1-INT x),INT x-1,INT x-1
     36                 lt zero            GAMMA(x+1-INT x),INT x-1,INT x<1?
     0017               jump true,G_LT_1   Jump if INT x < 1; ie. if x<1

The following deals with the remaining case;
ie. if x is greater than or equal to two.

EAA6 38                 end_calc           GAMMA(x+1-INT x),INT x-1
EAA7 CDD52D             CALL FP_TO_A       ;A= INT x-1
EAAA 47                 LD   B,A           ;B= INT x-1
EAAB EF                 RST  #28           GAMMA(x+1-INT x)
EAAC A1                 const one          GAMMA(x+1-INT x),1

At each pass around the following loop, the topmost item on the calculator
stack is assumed to be the "product-so-far", which I shall refer to as Pb.
At any stage this will actually be the product x(x-1)(x-2)... to INT x-1-B
terms. In addition, memory five will contain x-INT x+B+1.

                                           GAMMA(x+1-INT x),Pb
EAAD E5       G_LOOP_1  recall M5          GAMMA(x+1-INT x),Pb,x-INT x+B+1
EAAE A1                 const one          GAMMA(x+1-INT x),Pb,x-INT x+B+1,1
EAAF 03                 subtract           GAMMA(x+1-INT x),Pb,x-INT x+B
EAB0 C5                 store M5
EAB1 04                 multiply           GAMMA(x+1-INT x),P(b-1)
EAB2 35FA               djnz G_LOOP_1      GAMMA(x+1-INT x),P0

At this stage the topmost item on the calculator stack is x(x-1)(x-2)(x-3)...K,
where K is less than two (but greater than or equal to one).

EAB4 04                 multiply           GAMMA(x)
     38                 end_calc
EAB6 C9                 RET                ;Return with correct value
                                           ;on calculator stack

The following applies for the case INT x=1. Note that since INT x=1 it follows
that GAMMA(x+1-INT x) actually equals GAMMA(x), and that INT x-1 equals zero.

EAB7 EF                 RST  #28           GAMMA(x),0
EAB8 02       G_1_TO_2  delete             GAMMA(x)
EAB9 38                 end_calc
EABA C9                 RET                ;Return with correct value
                                           ;on calculator stack

The following deals with the case where x is less than one.

EABB EF                 RST  #28           GAMMA(x+1-INT x),INT x-1
EABC 1B       G_LT_1    negate             GAMMA(x+1-INT x),1-INT x
EABD 38                 end_calc
EABE CDD52D             CALL FP_TO_A       ;A= 1-INT x
EAC1 47                 LD   B,A           ;B= 1-INT x
EAC2 EF                 RST  #28           GAMMA(x+1-INT x)
EAC3 A1                 const one          GAMMA(x+1-INT x),1

At each pass around the following loop, the topmost item on the calculator
stack is assumed to be the "product-so-far", which I shall refer to as Pb.
At any stage this will actually be the product (x+1)(x+2)(x+3)... to 1-INT x-B
terms. In addition, memory five will contain x+1-INT x-B.

                                           GAMMA(x+1-INT x),Pb
EAC4 E5       G_LOOP_2  recall M5          GAMMA(x+1-INT x),Pb,x+1-INT x-B
EAC5 31                 duplicate          GAMMA(x+1-INT x),Pb,x+1-INT x-B,x+1-INT x-B
     A1                 const one          GAMMA(x+1-INT x),Pb,x+1-INT x-B,x+1-INT x-B,1
     0F                 add                GAMMA(x+1-INT x),Pb,x+1-INT x-B,x+2-INT x-B
EAC8 C5                 store M5
EAC9 02                 delete             GAMMA(x+1-INT x),Pb,x+1-INT x-B
EACA 04                 multiply           GAMMA(x+1-INT x),P(b-1)
EACB 35F8               djnz G_LOOP_2      GAMMA(x+1-INT x),P0

At this stage the topmost item on the calculator stack is (x+1)(x+2)(x+3)...K,
where K is less than one, but greater than or equal to zero.

EACD 05                 divide             GAMMA(x)
EACE 38                 end_calc
EACF C9                 RET                ;Return with correct value
                                           ;on calculator stack

Finally, note that numbers such that x-INT x>1/2 are dealt with separately,
for greater accuracy. These are numbers such as: -1.5 to -1; -0.5 to 0;
0.5 to 1; 1.5 to 2; and so on.

EAD0 EF                 RST  #28           x,x-INT x
EAD1 02       G_HALF    delete             x
EAD2 31                 duplicate          x,x
     A1                 const one          x,x,1
     01                 exchange           x,1,x
EAD5 03                 subtract           x,1-x
EAD6 38                 end_calc           x,1-x
EAD7 CD60EA             CALL GAMMA         ;Call the GAMMA s/rtn recursively
EADA EF                 RST  #28           x,GAMMA(1-x)
EADB 01                 exchange           GAMMA(1-x),x
     A3                 const PI/2         GAMMA(1-x),x,PI/2
     38                 end_calc
EADE 34                 INC  (HL)          GAMMA(1-x),x,PI
EADF EF                 RST  #28
EAE0 C4                 store M4           (M4 contains PI)
     04                 multiply           GAMMA(1-x),x*PI
     1F                 sin                GAMMA(1-x),SIN(x*PI)
EAE3 04                 multiply           GAMMA(1-x)*SIN(x*PI)
EAE4 E4                 recall M4          GAMMA(1-x)*SIN(x*PI),PI
     01                 exchange           PI,GAMMA(1-x)*SIN(x*PI)
     05                 divide             PI/(GAMMA(1-x)*SIN(x*PI))
EAE7 38                 end_calc           GAMMA(x)
EAE8 C9                 RET

It can be shown algebraically that PI/(GAMMA(1-x)*SIN(x*PI)) = GAMMA(x). The
proof of this is beyond the scope of this article, however I have made use of
this fact to ensure maximum accuracy. It means that the "series" instruction
is only ever passed values of x between zero and half. The last part of the
program links the machine code to the BASIC DEF FN g(x) statement.

EAE9 2A0B5C   GAMMA_FN  LD   HL,(DEFADD)   ;Point HL to DEF FN record
EAEC 23                 INC  HL
EAED 23                 INC  HL            ;HL points to 5-byte form of 'x'
EAEE CDB433             CALL STACK_NUM     ;Place x on calculator stack
EAF1 CD60EA             CALL GAMMA         ;Calculate GAMMA(x)
EAF4 C1                 POP  BC            ;Delete STACK_BC address from
                                           ;machine stack
EAF5 C9                 RET
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
