Z80:Advanced Math

Integer Math

Addition/Subtraction

To do addition or subtraction on an integer of arbitrary size, we use ADC.

Here we add two 32-bit integers:

   LD     HL, (word1)        ; Get least-significant word of word1
       LD     DE, (word2)        ; Get least-significant word of word2
       ADD    HL, DE              ; Add them
       LD     (result), HL        ; Store least-significant result
   
       LD     HL, (dword1 + 2)    ; Get most-significant word of word1
       LD     DE, (dword2 + 2)    ; Get most-significant word of word2
       ADC    HL, DE              ; Add them with the carry from the previous addition
       LD     (result + 2), HL    ; Store most-significant result
       RET
   
   word1:    .DB    $B3, $90, $12, $32    ; Each word is stored with the least-significant
   word2:    .DB    $F1, $15, $06, $B8    ; bytes first. You could just as easily have stored
   result:    .DB    $00, $00, $00, $00    ; them in big-endian, but because of how registers are
                                        ; loaded from RAM, it wouldn't work.

But now, what if you wanted to add truly arbitrary sized integers? We might do something like this:

   ld hl,Num1
        ld de,Num2
        ld bc,Size
   AddLoop:
   ;Inputs:
   ;     HL points to Num1 (stored little-endian)
   ;     DE points to Num2 (stored little-endian)
   ;     BC is the size of the integers
   ;Outputs:
   ;     BC is 0
   ;     HL points to the byte after Num1
   ;     DE points tot he byte after Num2
   ;     c flag is the carry
   ;     z flag is set
   ;     p/v flag reset
   ;     Num1 is overwritten with the result
          ld a,(de)
          adc a,(hl)
          ld (hl),a
          inc de
          cpi         ;increment HL, decrement BC, set parity flag of BC becomes 0
          jp pe,AddLoop
        ret
   Num1:
        .db $B3,$90,$12,$32,$47
   Num2:
        .db $F1,$15,$06,$B8,0

Similarly, one can perform arbitrary subtractions.

Multiplication

For multiplication, we could just add them in a DJNZ loop. But that is slow for multiplying big numbers. There is a much more efficient algorithm that you probably learned as a child in school that involves multiplying multi-digit numbers by hand. Similarly, we will do that here by using the binary representation of numbers and the astute reader will realise that binary digits are only 1 or 0 and multiplying by 1 or 0 is rather trivial!

To understand this algorithm, it is suggested that the reader try this by hand with some simple 4-bit numbers, then follow the algorithm or try to develop it in pseudocode.

Here is the first example:

   .module    DE_Times_A
   
   ;Note that we choose DE and A because of the readily available instructions that make this even more efficient.
   
   DE_Times_A:           ; HL = DE × A
       LD     HL, 0      ; Use HL to store the product
       LD     B, 8       ; Eight bits to check
   _loop:
       RRCA              ; Check least-significant bit of accumulator
       JR     NC, _skip  ; If zero, skip addition
       ADD    HL, DE
   _skip:
       SLA    E          ; Shift DE one bit left
       RL     D
   
   ;In binary, a shift left is multiplying by 2, just like a shift left in decimal is like multiplying by 10.
   ;We essentially added a 0 to the end of the number
   
       DJNZ   _loop
       RET

Or alternatively for a slightly faster and smaller code:

   .module    DE_Times_A
   DE_Times_A:      ; HL = DE × A
       ld hl,0      ; Use HL to store the product
       ld b,8       ; Eight bits to check
   _loop:
       add hl,hl
       rlca         ; Check most-significant bit of accumulator
       jr nc,_skip  ; If zero, skip addition
       add hl,de
   _skip:
       djnz _loop
       ret

And you might realise that multiplying an 8-bit number by a 16-bit number can yield up to a 24-bit number, so to take care of that overflow by turning rlca into adc a,a which is the same size and speed. Note that adc a,a can be replaced with rla for the same effect. Some overflow may occur from add hl,de, so we take care of that, too

   .module    DE_Times_A
   DE_Times_A:      ; AHL = DE × A
       ld hl,0      ; Use HL to store the product
       ld b,8       ; Eight bits to check
       ld c,h
   _loop:
       add hl,hl
       adc a,a      ; Check most-significant bit of accumulator
       jr nc,_skip  ; If zero, skip addition
       add hl,de
       adc a,c      ;add in overflow as needed.
   _skip:
       djnz _loop
       ret

Division

Division is also just like base-10 division as well, but instead of checking if the next digit is from 0 to 9, we only need to check if it is 0 or 1 which amounts to the question of "is our remainder smaller than the divisor (1) or not (0) ?"

Again, it is recommended that the reader try this with small examples to understand and possibly develop their own code.

   .module    Div_HL_D
   Div_HL_D:            ; HL = HL ÷ D, A = remainder
       XOR    A         ; Clear upper eight bits of AHL
       LD     B, 16     ; Sixteen bits in dividend
   _loop:
       ADD    HL, HL    ; Do a SLA HL. If the upper bit was 1, the c flag is set
       RLA              ; This moves the upper bits of the dividend into A
       JR     C, _overflow
       CP     D         ; Check if we can subtract the divisor
       JR     C, _skip  ; Carry means D > A
   _overflow:
       SUB    D         ; Do subtraction for real this time
       INC    L         ; Set the next bit of the quotient (currently bit 0)
   _skip:
       DJNZ   _loop
       RET

Fixed Point Math

Fixed Point refers to a number format that includes a fixed size integer part, and a fixed size fractional part. For example, 8.8 fixed point refers to a format where the upper 8 bits form the integer part, the lower 8 bits form the fractional part. Fixed point math often provides a desired balance between the numerical precision of floating point and the computing speed of integer maths.

To elaborate a little more, to represent 3.0 as a fixed point number stored in HL, we might have HL=0300, where the decimal point is between H and L (03.00). Then HL=0341 would be equivalent to 3+65/256=3.25390625. If we were to square this, we get HL=0A96 where we have lost some bits of precision.

Addition and Subtraction

As with integer arithmetic, addition and subtraction are trivial on the Z80. To add H.L+D.E, just use add hl,de, and for subtraction, 'or a \ sbc hl,de' will sufice.

Multiplication

H.L*D.E is essentially HL*DE/65536 If we view HL as a 16 bit integer, then H.L is essentially HL/256. For multiplication, the math for H.L*D.E is then: (HL/256)(DE/256) = (HL*DE)/65536. This means we can use a 16-bit signed multiplication routine that returns 32 bits and we now have a 16.16 fixed point result. To keep output as an 8.8 fixed point number, we need to grab the middle bytes around the decimal. If output was DE.HL, then E.H is a our 8.8 result. Overflow and underflow should be handled by the programmer depending on the task. note: Really, only a routine to return the lower 24-bits of the 16x16 multiplication is needed.

Allowing speed to take priority, here is a 16-bit unsigned integer multiplication:

   BC_Times_DE:
   ;  BC*DE->BHLA
   	ld a,b
   	ld hl,0
   	ld b,h
   	add a,a \ jr nc,$+5 \ ld h,d \ ld l,e
   	add hl,hl \ rla \ jr nc,$+4 \ add hl,de \ adc a,b
   	add hl,hl \ rla \ jr nc,$+4 \ add hl,de \ adc a,b
   	add hl,hl \ rla \ jr nc,$+4 \ add hl,de \ adc a,b
   	add hl,hl \ rla \ jr nc,$+4 \ add hl,de \ adc a,b
   	add hl,hl \ rla \ jr nc,$+4 \ add hl,de \ adc a,b
   	add hl,hl \ rla \ jr nc,$+4 \ add hl,de \ adc a,b
   	add hl,hl \ rla \ jr nc,$+4 \ add hl,de \ adc a,b
   	push hl
   	ld h,b
   	ld l,b
   	ld b,a
   	ld a,c
   	ld c,h
   	add a,a \ jr nc,$+5 \ ld h,d \ ld l,e
   	add hl,hl \ rla \ jr nc,$+4 \ add hl,de \ adc a,c
   	add hl,hl \ rla \ jr nc,$+4 \ add hl,de \ adc a,c
   	add hl,hl \ rla \ jr nc,$+4 \ add hl,de \ adc a,c
   	add hl,hl \ rla \ jr nc,$+4 \ add hl,de \ adc a,c
   	add hl,hl \ rla \ jr nc,$+4 \ add hl,de \ adc a,c
   	add hl,hl \ rla \ jr nc,$+4 \ add hl,de \ adc a,c
   	add hl,hl \ rla \ jr nc,$+4 \ add hl,de \ adc a,c    
   	pop de
   	ld c,a
   	ld a,l
   	ld l,h
   	ld h,c
   	add hl,de
   	ret nc
   	inc b
   	ret

Division

For division, we should notice that H.L/D.E evaluates as: (HL/256)/(DE/256)=HL/DE We could use a standard signed integer division routine, however, unlike with multiplication, the result can only be used a 16.0 integer. We need 8 more bits of accuracy for 16.8 output (of which we only use the lower 8-bits of the integer part). Since we are going to discard the upper 8 bits of the acquired result, we can optimise the first 8 iterations to not store the bits anywhere.

   DE_Div_BC_88:
   ;Inputs:
   ;     DE,BC are 8.8 Fixed Point numbers
   ;Outputs:
   ;     DE is the 8.8 Fixed Point result (rounded to the least significant bit)
       ld a,d
       xor b
       push af
       call m,negDE
       call div8sub
       pop af
       ret p
   negDE:
       xor a
       sub e
       ld e,a
       sbc a,a
       sub d
       ld d,a
       ret
   div8sub:
       ld a,8
       ld hl,0
   Loop1:
       rl d
       adc hl,hl
       sbc hl,bc
       jr nc,$+3
   	add hl,bc
       dec a
       jr nz,Loop1
       ld d,e
       ld e,a
       ld a,16
  	jp $+6
   DivLoop:
       add hl,bc
       dec a
       ret z
       sla e
       rl d
       adc hl,hl
   	jr c,overflow
       sbc hl,bc
       jr c,DivLoop
       inc e
       jp DivLoop+1
   overflow:
   	or a
   	sbc hl,bc
   	inc e
   	jp DivLoop

Inverse (1/x)

If we want to invert a number >1, we know that its result will be <1, so for 8.8 fixed point, this means there will only be 8 bits in the result. We can optimise out most of the cycles of the division, making it much faster. In the DE_Div_BC_88 routine, we initialise HL=0 and rotate a bit of DE into HL at each iteration. If DE=100h=1.0, then no bits will get rotated into HL until the 8th iteration. Then after that, HL=1, but we know BC>256, so we are guaranteed to have nothing happen for another 8 iterations until HL=256. At this point, we have only 8 more iterations to compute, and DE is 0, so we no longer need to even rotate DE into HL. So now at this point, we can just initialise HL=256, then rotate in a 0 each iteration. This can be done with add hl,hl, then we compare to BC to see if the remainder (here, HL). However, since we now only need two register pairs, lets use B for the iteration counter, and DE is the denominator. If we compare HL to DE, we can use sbc hl,de. If the C flag is set, it means HL<DE, so we need to rotate a 0 into our accumulator and then add DE to HL again. We could use sbc hl,de \ jr nc,$+3 \ add hl,de \ ccf \ rla, but if we note that we have a cpl instruction to invert the bits of a:

   DEgt1_Inv:
   ;1/D.E->B.C
   ; 470 to 518 t-states
   	ld hl,256
   	ld b,8
   InvLoop:
   	add hl,hl
   	sbc hl,de
   	jr nc,$+3
   	add hl,de
   	rla
   	djnz InvLoop
   	cpl
   	ld c,a
   	ret

And to unroll, saving 102 t-states:

   DEgt1_Inv:
   ;1/D.E->B.C
   ; 343 to 397 t-states
   	ld hl,512
   	ld b,l
   InvLoop:
   	or a       \ sbc hl,de \ jr nc,$+3 \ add hl,de \ adc a,a
   	add hl,hl \ sbc hl,de \ jr nc,$+3 \ add hl,de \ adc a,a
   	add hl,hl \ sbc hl,de \ jr nc,$+3 \ add hl,de \ adc a,a
   	add hl,hl \ sbc hl,de \ jr nc,$+3 \ add hl,de \ adc a,a
   	add hl,hl \ sbc hl,de \ jr nc,$+3 \ add hl,de \ adc a,a
   	add hl,hl \ sbc hl,de \ jr nc,$+3 \ add hl,de \ adc a,a
   	add hl,hl \ sbc hl,de \ jr nc,$+3 \ add hl,de \ adc a,a
   	add hl,hl \ sbc hl,de \ adc a,a
   	cpl
   	ld c,a
   	ret

Note that this actually works as an HL_Div_DE_88 routine, by removing the ld hl,256 and making sure HL<DE.

Square Root

Using the integer square root of HL, sqrt(H.L) would be returned in 4.4 fixed point format. The upper 4 bits of the integer part will be 0 regardless, so we can essentially say that we have 8.4 Fixed Point, so we need to modify the algorithm to iterate 4 more times to get 8.8 fixed point output. The following is adapted from a speed-optimised routine used in a floating point library (which required 16 bits instead of 12):

   SqrtHL_prec12:
   ;input: HL
   ;Output: HL
   ;183 bytes
   	xor a
   	ld b,a
   
   	ld e,l
   	ld l,h
   	ld h,a
   
   	add hl,hl
   	add hl,hl
   	cp h
   	jr nc,$+5
   	dec h
   	ld a,4
   
   	add hl,hl
   	add hl,hl
   	ld c,a
   	sub h
   	jr nc,$+6
   	cpl
   	ld h,a
   	inc c
   	inc c
   
   	ld a,c
   	add hl,hl
   	add hl,hl
   	add a,a
   	ld c,a
   	sub h
   	jr nc,$+6
   	cpl
   	ld h,a
   	inc c
   	inc c
   
   	ld a,c
   	add hl,hl
   	add hl,hl
   	add a,a
   	ld c,a
   	sub h
   	jr nc,$+6
   	cpl
   	ld h,a
   	inc c
   	inc c
   
   	ld a,c
   	ld l,e
   
   	add hl,hl
   	add hl,hl
   	add a,a
   	ld c,a
   	sub h
   	jr nc,$+6
   	cpl
   	ld h,a
   	inc c
   	inc c
   
   	ld a,c
   	add hl,hl
   	add hl,hl
   	add a,a
   	ld c,a
   	sub h
   	jr nc,$+6
   	cpl
   	ld h,a
   	inc c
   	inc c
   
   	ld a,c
   	add a,a \ ld c,a
   	add hl,hl
   	add hl,hl
   	jr nc,$+6
   	sub h \ jp $+6
   	sub h
   	jr nc,$+6
   	inc c \ inc c
   	cpl
   	ld h,a
   
   
   	ld a,l
   	ld l,h
   	add a,a
   	ld h,a
   	adc hl,hl
   	adc hl,hl
   	sll c \ rl b
   	sbc hl,bc
   	jr nc,$+3
   	add hl,bc
   	sbc a,a \ add a,a \ inc a \ add a,c \ ld c,a
   
   ;iteration 9
   	add hl,hl \ add hl,hl
   	sll c \ rl b
   	sbc hl,bc
   	jr nc,$+3
   	add hl,bc
   	sbc a,a \ add a,a \ inc a \ add a,c \ ld c,a
   
   	add hl,hl \ add hl,hl
   	sll c \ rl b
   	sbc hl,bc
   	jr nc,$+3
   	add hl,bc
   	sbc a,a \ add a,a \ inc a \ add a,c \ ld c,a
   
   	add hl,hl \ add hl,hl
   	sll c \ rl b
   	sbc hl,bc
   	jr nc,$+3
   	add hl,bc
   	sbc a,a \ add a,a \ inc a \ add a,c \ ld c,a
   
   	add hl,hl \ add hl,hl
   	sll c \ rl b
   	sbc hl,bc
   	jr nc,$+3
   	add hl,bc
   	sbc a,a \ add a,a \ inc a \ add a,c \ ld c,a
   ;12th iteration completed
   ; output in BC
       srl b \ rr c
       ld h,b
       ld l,c
       ret

Exponentiation and Logarithms

Exponentiation

Here, exponentiation refers to <math>2^{x}</math>. Observe that <math>2^{39/16}=2^{2+7/16} = 2^{2+0/2+1/4+1/8+1/16}=2^{2}2^{0/2}2^{1/4}2^{1/8}2^{1/16}</math> The observant reader will see that taking away the trivial integer part of the power, the fractional part can be represented as a sequence of square roots with 2s wherever there are 1s in the binary expansion. So to compute <math>2^{x}</math> where x is on [0,1), in pseudo code we have:

   If x==0
     return 1.0
   repeat until carry is set
     x>>=1    ;shift X right
   y=sqrt(2)
   while x>0
     x>>=1
     if (carry is set)
       y*=2
     y=sqrt(y)

Converting that to assembly:

   ;Inputs:
   ;     HL is the 8.8 fixed point number 'x' for 2^x
   ;Outputs:
   ;     DEHL is the 24.8 fixed point result. If there was overflow exceeding 2^24, then this value is set to the max.
        ld a,l
        or a
        push hl     ;save H for later, H is the integer part of the power
        ld hl,1
        jr z,integerpow
        scf      ;set the carry flag so that a bit is rotated into a. This will act as our counter.
   ;wait until we come across the lowest bit. Also note that we 
        rra
        jr nc,$-1
        ld hl,2*256
   powloop:
        push af
        call FPSqrtHL    ;returns in HL
        pop af
        srl a
        jr z,integerpow
        jr nc,powloop
        add hl,hl
        jp powloop
   integerpow:
        pop bc
   ;Now b is the integer part for 2^x that we need to multiply HL by.
        ld de,0
        ld a,b
        or a
        ret z
   
        add hl,hl
        rl e \ rl d \ jr c,wayoverflow
        djnz $-7
        ret
   wayoverflow:
        ld hl,-1
        ld d,h
        ld e,l
        ret

Note that in the spirit of optimisation, that we can significantly speed up the fixed point sqrtHL routine since HL is never greater than 4.0.

Logarithms

Here, logarithm refers to log,,2,, also written as lg(). Using logarithms will allow us to extend the exponential routine to other bases, such as 3^n without much difficulty. Here is a fast, fixed point log_2 routine, optimised for size:

   Log_2_88_size:
   ;Inputs:
   ;     HL is an unsigned 8.8 fixed point number.
   ;Outputs:
   ;     HL is the signed 8.8 fixed point value of log base 2 of the input.
   ;Example:
   ;     pass HL = 3.0, returns 1.58203125 (actual is ~1.584962501...)
   ;averages about 39 t-states slower than original
   ;62 bytes
        ex de,hl
        ld hl,0
        ld a,d
        ld c,8
        or a
        jr z,DE_lessthan_1
        srl d
        jr z,logloop-1
        inc l
        rr e
        jr $-7
   DE_lessthan_1:
        ld a,e
        dec hl
        or a
        ret z
        inc l
        dec l
        add a,a
        jr nc,$-2
        ld e,a
   
        inc d
   logloop:
        add hl,hl
        push hl
        ld h,d
        ld l,e
        ld a,e
        ld b,8
   
        add hl,hl
        rla
        jr nc,$+5
          add hl,de
          adc a,0
        djnz $-7
   
        ld e,h
        ld d,a
        pop hl
        rr a         ;this is NOT supposed to be rra, we need the z flag affected
        jr z,$+7
          srl d
          rr e
          inc l
        dec c
        jr nz,logloop
        ret

Arctangent and Natural Logarithm

Although these are not the fastest methods, the approach to the problem of computing these functions efficiently is interesting. As well, they can work faster than CORDIC or AGM methods. Given efficient routines for computing multiplication and division, the following are accurate to better than 8 bits of accuracy:

   atan(x) : x(9+2x^2)/(9+5x^2) on [-1,1]
   ln(x+1) : x(8+x)/(8+5x) on [0,1]

For a more complete analysis of how these formulas were derived, see this document. What is convenient to note is that the constant multiplications can be reduced to a few shift and add instuctions. So to put the functions into code:

   ArcTan_88:
   ;Inputs: DE is positive
       ld a,d
       or a
       jr z,ArcTan+2
       inc a
       jr z,ArcTan+2
       ld b,d
       ld c,e
       ld de,256
       call DE_Div_BC_88
       call ArcTan
       ld de,402     ;about pi/2
       ex de,hl
       or a
       sbc hl,de
       ret
   ArcTan:
   ;Input: DE
   ;Output: HL
       push de    ;save X for later
       ld b,d
       ld c,e
       call BC_Times_DE ;x^2 computed
       ld d,h
       ld e,l
       add hl,hl
       add hl,hl
       add hl,de   ;5x^2 is now in HL
       ld a,9
       add a,h
       ld h,a      ;9+5x^2->HL
       ld a,9
       sla e
       rl d
       add a,d
       ld d,a      ;9+2x^2->DE
       ld b,h
       ld c,l
       call DE_Div_BC_88   ;->DE     
       pop bc
       jp BC_Times_DE

For the following ln() routine, the approximating function is employed taking advantage of lots of optimisations, putting it between the computing speed of multiplication and division:

   lognat:
   ;Input:  H.L needs to be on [1,2]
   ;Output: H.L if z flag is set, else if nz, no result
   ;avg speed: 677 t-states
   	dec h
      	dec h
   	jr nz,$+8
   	inc l \ dec l
   	ret nz
   	ld l,177
   	ret
   	inc h
   	ret nz
   	ld b,h
   	ld c,l
   	ld e,l
   	ld d,8
   	add hl,hl
   	add hl,hl
   	add hl,de
   	ex de,hl
   
   	add hl,hl \ sbc hl,de \ jr nc,$+3 \ add hl,de \ adc a,a
   	add hl,hl \ sbc hl,de \ jr nc,$+3 \ add hl,de \ adc a,a
   	add hl,hl \ sbc hl,de \ jr nc,$+3 \ add hl,de \ adc a,a
   	add hl,hl \ sbc hl,de \ jr nc,$+3 \ add hl,de \ adc a,a
   	add hl,hl \ sbc hl,de \ jr nc,$+3 \ add hl,de \ adc a,a
   	add hl,hl \ sbc hl,de \ jr nc,$+3 \ add hl,de \ adc a,a
   	add hl,hl \ sbc hl,de \ jr nc,$+3 \ add hl,de \ adc a,a
   	add hl,hl \ sbc hl,de \ adc a,a
   
   	ld h,a \ ld l,b
   	sla h \ jr c,$+3 \ ld l,c
   	add hl,hl \ jr c,$+3 \ add hl,bc
   	add hl,hl \ jr c,$+3 \ add hl,bc
   	add hl,hl \ jr c,$+3 \ add hl,bc
   	add hl,hl \ jr c,$+3 \ add hl,bc
   	add hl,hl \ jr c,$+3 \ add hl,bc
   	add hl,hl \ jr c,$+3 \ add hl,bc
   	add hl,hl \ jr c,$+3 \ add hl,bc
   	rl l
   	ld a,h
   	adc a,b
   	ld l,a
   	ld h,b
   	cp a
   	ret

But that only works on [1,2], so we have to reduce speed a little to extend it to all positive inputs (0,128):

   lognat:
   ;Input:  H.L needs to be on (0,128.0)
   ;Output: H.L if c flag set
   ;	 returns nc if input is negative (HL not modified)
   ;Error:
   ;	The error on the outputs is as follows:
   ;		20592 inputs are exact
   ;		12075 inputs are off by 1/256
   ;		100 inputs are off by 2/256
   ;	So all 32767 inputs are within 2/256, with average error being <1/683 which is smaller than 1/256.
   ;Size: 177 bytes
   ;Speed: average speed is less than 1250 t-states
   
   	ld a,h \ or l \ jr nz,$+5
   	ld h,80h \ ret
   	dec h
   	dec h
   	jr nz,$+9
   	inc l \ dec l
   	jr nz,normalizeln-1
   	ld l,177
   	ret
   	inc h
   	jr nz,normalizeln
   	ld b,h
   	ld c,l
   	ld e,l
   	ld d,8
   	add hl,hl
   	add hl,hl
   	add hl,de
   	ex de,hl
   	call HL_Div_DE
   	ld h,a \ ld l,b
   	sla h \ jr c,$+3 \ ld l,c
   	add hl,hl \ jr c,$+3 \ add hl,bc
   	add hl,hl \ jr c,$+3 \ add hl,bc
   	add hl,hl \ jr c,$+3 \ add hl,bc
   	add hl,hl \ jr c,$+3 \ add hl,bc
   	add hl,hl \ jr c,$+3 \ add hl,bc
   	add hl,hl \ jr c,$+3 \ add hl,bc
   	add hl,hl \ jr c,$+3 \ add hl,bc
   	rl l
   	ld a,h
   	adc a,b
   	ld h,b
   	ld l,a
   	scf
   	ret
   HL_Div_DE:
   	add hl,hl \ sbc hl,de \ jr nc,$+3 \ add hl,de \ adc a,a
   	add hl,hl \ sbc hl,de \ jr nc,$+3 \ add hl,de \ adc a,a
   	add hl,hl \ sbc hl,de \ jr nc,$+3 \ add hl,de \ adc a,a
   	add hl,hl \ sbc hl,de \ jr nc,$+3 \ add hl,de \ adc a,a
   	add hl,hl \ sbc hl,de \ jr nc,$+3 \ add hl,de \ adc a,a
   	add hl,hl \ sbc hl,de \ jr nc,$+3 \ add hl,de \ adc a,a
   	add hl,hl \ sbc hl,de \ jr nc,$+3 \ add hl,de \ adc a,a
   	add hl,hl \ sbc hl,de \ adc a,a \ ret
   
   	inc h
   normalizeln:
   	xor a
   	inc h \ ret m
   	ld d,a \ ld e,a
   	ld a,l
   	jr z,toosmall
   	inc e \ srl h \ rra \ jr nz,$-4
   	rla \ rl h
   	dec e
   stepin:
   	ld l,a
   	push de
   	call lognat
   	pop de
   	;now multiply DE by 355, then divide by 2 (rounding)
   	ld b,d \ ld c,e \ ld a,d
   	ex de,hl
   	add hl,hl
   	add hl,hl	;4
   	add hl,bc	;5
   	add hl,hl	;10
   	add hl,bc	;11
   	add hl,hl	;22
   	add hl,hl
   	add hl,hl
   	add hl,hl
   	add hl,bc
   	add hl,hl
   	add hl,bc
   	sra h \ rr l
   	adc hl,de
   	scf
   	ret
   toosmall:
   	dec d
   	dec e \ add a,a \ jr nc,$-2
   	inc h
   	jp stepin

Now it averages about 1250 t-states, still faster than some implementations of division.

If you pass HL=768, in fixed point that is 3.0, so output is 282 which is 1.1015625 in fixed point.

Sine and Cosine

Unlike atan() and ln(), these have quickly converging Taylor series. From these, we can derive:

   sin(x) : x-85x^3/512+x^5/128 is within 8 bits of accuracy on [-pi/2,pi/2]
   cos(x) : 1-x^2/2+5x^4/128 is within 8 bits of accuracy on [-pi/2,pi/2] for almost all values

For sine, noting that 85=17*5 will help with fast constant multiplication:

   sine_88:
   ;Inputs: de
       push de
       sra d \ rr e
       ld b,d \ ld c,e
       call BC_Times_DE
       push hl     ;x^2/4
       sra h \ rr l
       ex de,hl
       ld b,d \ ld c,e
       call BC_Times_DE
       sra h \ rr l
       inc h
       ex (sp),hl    ;x^4/128+1 is on stack, HL=x^2/4
       xor a
       ld d,a
       ld b,h
       ld c,l
       add hl,hl \ rla
       add hl,hl \ rla
       add hl,bc \ adc a,d
       ld b,h
       ld c,l
       add hl,hl \ rla
       add hl,hl \ rla
       add hl,hl \ rla
       add hl,hl \ rla
       add hl,bc \ adc a,d
       ld e,l
       ld l,h
       ld h,a
       rl e
       adc hl,hl
       rl e
       jr nc,$+3
       inc hl
   
       pop de
       ex hl,de
       or a
       sbc hl,de
       ex de,hl
       pop bc
       jp BC_Times_DE

Table-Based Fixed 8.8

When speed is priority, using table-based routines are often the fastest way to achieve a task.

Logarithms

Because we are working in binary, lg(), which is log,,2,,() is the easiest and most practical to compute. From this, we can derive other bases, too. Now to use some identities, note that lg(x/2)=lg(x)-lg(2) = lg(x)-1. So then we have lg(x)=1+lg(x/2). So what we will do is make a routine that works on [1,2] and extend that to all other values. If our input is less than 1, we have lg(x)=lg(x*2)-1, so we double our input, and subtract from an accumulator. If input is >1, we will divide it by 2 and add to the accumulator. Then, we only have 256 values to worry about in our table. But, if our input is now even, we can divide by 2 without losing any accuracy (no bits lost) and just increment the accumulator. We can keep doing this until we finally have an odd number (which means we need to take care of input=0 as a separate case). Then we have a table of 128 elements instead--all of the odd numbers. And now we can build our algorithm:

   lg_88:
   ;Input: HL is a fixed point number
   ;Output: lg(H.L)->H.L
   ;Speed: Avg: 340
    ld de,lgLUT
    ld b,0
    ld a,h
    or a
    ret m     ;If input is negative, this algorithm won't suffice (this is only for the Reals, not Complex)
    ld a,l
    jr z,$+8
    inc b \ srl h \ rra \ jr nz,$-4   ;input >2, so we rotate, using b as the accumulator
    or a \ jr nz,$+6
    ld hl,8000h \ ret     ;If input was 0, use this as -inf
    rra \ inc b \ jr nc,$-2      ;input<1, so we rotate
   ;A is the element to look up in the LUT
    ld l,a
    ld c,h
    dec b
    add hl,hl
    add hl,de
    ld e,(hl)
    inc hl
    ld d,(hl)
    ex de,hl
    add hl,bc
    ret
   lglut:
   .dw $F800
   .dw $F996
   .dw $FA52
   .dw $FACF
   .dw $FB2C
   .dw $FB76
   .dw $FBB3
   .dw $FBE8
   .dw $FC16
   .dw $FC3F
   .dw $FC64
   .dw $FC86
   .dw $FCA5
   .dw $FCC1
   .dw $FCDC
   .dw $FCF4
   .dw $FD0B
   .dw $FD21
   .dw $FD36
   .dw $FD49
   .dw $FD5C
   .dw $FD6D
   .dw $FD7E
   .dw $FD8E
   .dw $FD9D
   .dw $FDAC
   .dw $FDBA
   .dw $FDC8
   .dw $FDD5
   .dw $FDE2
   .dw $FDEE
   .dw $FDFA
   .dw $FE06
   .dw $FE11
   .dw $FE1C
   .dw $FE26
   .dw $FE31
   .dw $FE3B
   .dw $FE44
   .dw $FE4E
   .dw $FE57
   .dw $FE60
   .dw $FE69
   .dw $FE71
   .dw $FE7A
   .dw $FE82
   .dw $FE8A
   .dw $FE92
   .dw $FE9A
   .dw $FEA1
   .dw $FEA9
   .dw $FEB0
   .dw $FEB7
   .dw $FEBE
   .dw $FEC5
   .dw $FECB
   .dw $FED2
   .dw $FED8
   .dw $FEDF
   .dw $FEE5
   .dw $FEEB
   .dw $FEF1
   .dw $FEF7
   .dw $FEFD
   .dw $FF03
   .dw $FF09
   .dw $FF0E
   .dw $FF14
   .dw $FF19
   .dw $FF1E
   .dw $FF24
   .dw $FF29
   .dw $FF2E
   .dw $FF33
   .dw $FF38
   .dw $FF3D
   .dw $FF42
   .dw $FF47
   .dw $FF4B
   .dw $FF50
   .dw $FF55
   .dw $FF59
   .dw $FF5E
   .dw $FF62
   .dw $FF67
   .dw $FF6B
   .dw $FF6F
   .dw $FF74
   .dw $FF78
   .dw $FF7C
   .dw $FF80
   .dw $FF84
   .dw $FF88
   .dw $FF8C
   .dw $FF90
   .dw $FF94
   .dw $FF98
   .dw $FF9B
   .dw $FF9F
   .dw $FFA3
   .dw $FFA7
   .dw $FFAA
   .dw $FFAE
   .dw $FFB2
   .dw $FFB5
   .dw $FFB9
   .dw $FFBC
   .dw $FFC0
   .dw $FFC3
   .dw $FFC6
   .dw $FFCA
   .dw $FFCD
   .dw $FFD0
   .dw $FFD4
   .dw $FFD7
   .dw $FFDA
   .dw $FFDD
   .dw $FFE0
   .dw $FFE4
   .dw $FFE7
   .dw $FFEA
   .dw $FFED
   .dw $FFF0
   .dw $FFF3
   .dw $FFF6
   .dw $FFF9
   .dw $FFFC
   .dw $FFFF

Adding a little overhead code, we can get an ln() routine. Note that to convert, we have ln(x)=lg(x)/lg(e), but 1/lg(e) is a constant, so we can use constant multiplication. (it is approximately 355/512 for a good estimate):

   ln_88:
   ;Input: HL is a fixed point number
   ;Output: ln(H.L)->H.L
   ;Speed: Avg: 340+(325 worst case)
    call lg_88
    ;now signed multiply HL by 355, then divide by 2 (rounding)
       ld de,0
       bit 7,h
       jr z,$+9
       dec e \ xor a \ sub l \ ld l,a
       sbc a,a \ sub h \ ld h,a
       ld b,h
       ld c,l
       xor a
       add hl,hl
       add hl,hl \ rla
       add hl,bc \ adc a,d
       add hl,hl \ rla
       add hl,bc \ adc a,d
       add hl,hl \ rla
       add hl,hl \ rla
       add hl,hl \ rla
       add hl,hl \ rla
       add hl,bc \ adc a,d
       add hl,hl \ rla
       add hl,bc \ adc a,d
       sra a \ rr h
       ld l,h
       ld h,a
       inc e
       ret nz
       xor a \ sub l \ ld l,a
       sbc a,a \ sub h \ ld h,a
       ret

Arctan

Arctan also can be computed on the whole real line just from the values on [0,1] by the identities: arctan(-x)=-arctan(x) and arctan(x)=pi/2-arctan(1/x). Unlike ln() where we need an LUT (Look Up Table) of 16-bit numbers, we can just use 8-bit numbers. Plus, since arctan is pretty close to a linear function for a while, a large chunk of the LUT can be omitted by checking if the input is within a certain range:

   arctan_88:
   ;Input:
   ;	D.E
   ;Output: atan(D.E)->D.E
   	push de
   	ld a,d
   	or a
   	jp p,$+5
   	neg
   	ld d,a
   	dec a
   	jr nz,checkneedinv
   	inc e \ dec e \ jr nz,checkneedinv
   	pop af \ rla \ ld de,201 \ ret nc \ ld de,-201 \ ret
   checkneedinv:
   	inc a
   	call nz,DEgt1_Inv
   ;0.E is the value to atan
   	ld hl,adjustatan
   	push hl
   	ld a,e
   	cp 46 \ ret c
   	dec a \ cp 42h \ ret c
   	dec a \ cp 4Eh \ ret c
   	dec a \ cp 57h \ ret c
   	dec a \ cp 5Eh \ ret c
   	dec a \ cp 64h \ ret c
   	dec a \ cp 6Ah \ ret c
   	dec a \ cp 6Fh \ ret c
   	sub 6Fh \ ld e,a
   	ld hl,atanlut
   	add hl,de
   	ld a,(hl)
   	ret
   adjustatan:
   	ld e,a
   	pop bc
   	ld a,b
   	or a
   	jp p,$+5
   	neg
   	jr z,$+9
   	ld hl,402
   	or a
   	sbc hl,de
   	ex de,hl
   	rl b
   	ret nc
   	xor a
   	sub e
   	ld e,a
   	sbc a,a
   	sub d
   	ld d,a
   	ret
   DEgt1_Inv:
   ;Works if DE>1
   	ld hl,256
   	ld b,8
   InvLoop:
   	add hl,hl
   	sbc hl,de
   	jr nc,$+3
   	add hl,de
   	adc a,a
   	djnz InvLoop
       cpl
   	ld e,a
       ld d,b
       ret
   atanlut:
   .db $6F
   .db $6F
   .db $70
   .db $71
   .db $72
   .db $73
   .db $73
   .db $74
   .db $75
   .db $76
   .db $77
   .db $77
   .db $78
   .db $79
   .db $7A
   .db $7B
   .db $7B
   .db $7C
   .db $7D
   .db $7E
   .db $7F
   .db $7F
   .db $80
   .db $81
   .db $82
   .db $82
   .db $83
   .db $84
   .db $85
   .db $85
   .db $86
   .db $87
   .db $88
   .db $88
   .db $89
   .db $8A
   .db $8B
   .db $8B
   .db $8C
   .db $8D
   .db $8E
   .db $8E
   .db $8F
   .db $90
   .db $90
   .db $91
   .db $92
   .db $93
   .db $93
   .db $94
   .db $95
   .db $95
   .db $96
   .db $97
   .db $97
   .db $98
   .db $99
   .db $9A
   .db $9A
   .db $9B
   .db $9C
   .db $9C
   .db $9D
   .db $9E
   .db $9E
   .db $9F
   .db $A0
   .db $A0
   .db $A1
   .db $A2
   .db $A2
   .db $A3
   .db $A3
   .db $A4
   .db $A5
   .db $A5
   .db $A6
   .db $A7
   .db $A7
   .db $A8
   .db $A9
   .db $A9
   .db $AA
   .db $AA
   .db $AB
   .db $AC
   .db $AC
   .db $AD
   .db $AD
   .db $AE
   .db $AF
   .db $AF
   .db $B0
   .db $B0
   .db $B1
   .db $B2
   .db $B2
   .db $B3
   .db $B3
   .db $B4
   .db $B5
   .db $B5
   .db $B6
   .db $B6
   .db $B7
   .db $B7
   .db $B8
   .db $B9
   .db $B9
   .db $BA
   .db $BA
   .db $BB
   .db $BB
   .db $BC
   .db $BC
   .db $BD
   .db $BE
   .db $BE
   .db $BF
   .db $BF
   .db $C0
   .db $C0
   .db $C1
   .db $C1
   .db $C2
   .db $C2
   .db $C3
   .db $C3
   .db $C4
   .db $C4
   .db $C5
   .db $C6
   .db $C6
   .db $C7
   .db $C7
   .db $C8
   .db $C8
   .db $C9

Random Numbers

16-bit Pseudo-Random Number Generator

The previous routine in this place was 1592cc and 280 bytes (compare to 291cc, 58 bytes). It worked well enough in terms of randomness, but there were some issues that merited a replacement. The following was tested and passed all tests at the CAcert Research Lab (see zrng2 by Zeda Thomas). The data that was output from this routine was from calling this routine to generate two bytes at a time to generate 16MiB of data.

Usage Proper use of this routine is to use the upper bits. The lower bits are not nearly as "random" as the upper bits. With this said, if you want to genreate a number on a range from [0,N), treat the output as a 0.16 fixed-point number, and multiply by N and take the integer part. In assembly, this translated to multiplying the output by N and taking the overflow (bits[16:]).

   #ifndef smc     ;\
   #define smc     ; |Remove or comment for Apps or if you *really* don't want SMC
   #end            ;/
   #ifndef smc             ;\
   #define seed1_0 8000h   ; |If you are not using SMC, you'll need to define
   #define seed1_1 8002h   ; |these vars yourself.
   #define seed2_0 8004h   ; |
   #define seed2_1 8006h   ; |
   #endif                  ;/
   rand16:
   ;Inputs:
   ;   seed1
   ;   seed2>0
   ;Output:
   ;   HL
   ;Destroys:
   ;   A,BC,DE. Could all be used as 'less random' options, too.
   ;Notes:
   ;   Tested and passes all CAcert Labs tests.
   ;   Has a period of 18,446,744,069,414,584,320 (roughly 18.4 quintillion)
   ;   This would take over 11.3 million years at 15MHz to finish a full period.
   ;291cc, +32cc if not using smc
   ;58 bytes, +2 if not using smc
   ;;lcg
   #ifdef smc
   seed1_0=$+1
       ld hl,12345
   seed1_1=$+1
       ld de,6789
   #else
       ld hl,(seed1_0)
       ld de,(seed1_1)
   #endif
       ld b,h
       ld c,l
       add hl,hl \ rl e \ rl d
       add hl,hl \ rl e \ rl d
       inc l
       add hl,bc
       ld (seed1_0),hl
       ld hl,(seed1_1)
       adc hl,de
       ld (seed1_1),hl
       ex de,hl
   ;;lfsr
   #ifdef smc
   seed2_0=$+1
       ld hl,9876
   seed2_1=$+1
       ld bc,54321
   #else
       ld hl,(seed2_0)
       ld bc,(seed2_1)
   #endif
       add hl,hl
       rl c
       rl b
       ld (seed2_1),bc
       sbc a,a
       and %11000101
       xor l
       ld l,a
       ld (seed2_0),hl
       ex de,hl
       add hl,bc
       ret

Conclusion

For more info, see the TI-83 Asm in 28 Days, Day 15.

When multiplying by a constant, there are often more efficient ways than using these subprograms. See here for more information.

For ways to perform more accurate math (Floating Point), there are many prebuilt functions that operate on the OP operators. They can be found here]. We also have Floating Point routines that cover other formats that may be of more use in other situations.

Sources: ASM in 28 Days, Day 15 FP_Math