This routine takes the cube root of a 24-bit unsigned number and returns the cube root and the remainder. Note that:
number = root * root * root + remainder
- Inputs:
- NUMH = bits 23-16 of the number
- NUMM = bits 15-8 of the number
- NUML = bits 7-0 of the number
- Outputs:
- ROOT = the cube root of NUM
- REMH = bits 23-16 of the number
- REMM = bits 15-8 of the number
- REML = bits 7-0 of the number
This routine requires 3 bytes of temporary RAM storage: DIVL, DIVM, and DIVH.
;
LDA #0
STA ROOT
STA REML
STA REMM
STA REMH
STA DIVL
STA DIVM
STA DIVH
JSR L1
L1 JSR L2
L2 JSR L3
L3 LDX #3
L4 ASL NUML
ROL NUMM
ROL NUMH
ROL REML
ROL REMM
ROL REMH
DEX
BNE L4
;
; Note: carry is clear here!
;
LDA REML
SBC DIVL
TAX
LDA REMM
SBC DIVM
TAY
LDA REMH
SBC DIVH
BCC L9
STA REMH
STY REMM
STX REML
ROL ROOT
LDX #0
LDA ROOT
ASL
BCC L5
LDX #2
L5 ASL
BCC L6
INX
CLC
L6 LDY #2
L7 ADC DIVL
STA DIVL
TXA
ADC DIVM
STA DIVM
BCC L8
INC DIVH
L8 ASL DIVL
ROL DIVM
ROL DIVH
LDA ROOT
LDX #0
DEY
BNE L7
RTS
L9 LDX #2
L10 SEC
LDA DIVL
SBC ROOT
STA DIVL
BCS L12
LDA DIVM
BNE L11
DEC DIVH
L11 DEC DIVM
L12 ASL ROOT
DEX
BNE L10
;
; undo the ASL ROOT from the previous loop
; (new ROOT = old ROOT * 2, not old ROOT * 4)
;
ROR ROOT
ASL DIVL
ROL DIVM
ROL DIVH
ASL DIVL
ROL DIVM
ROL DIVH
RTS
How it works:
Calculating a cube root digit-by-digit is similar to calculating a square root digit-by-digit. There are two main differences.
First, for a square root, two digits are shifted in (brought down); for a cube root, three digits are shifted in
Second, for a square root
divisor = 4 * root so far + 1
since:
(2*x+y)*(2*x+y) = 4*x*x + 4*x*y + y*y
For a cube root:
divisor = 12 * root so far * root so far + 6 * root so far + 1
since:
(2*x+y)*(2*x+y)*(2*x+y) = 2*x*x*x + 12*x*x*y + 6*x*y*y + y*y*y
Thus, the algorithm for the cube root is:
cube_root (n)
dividend = n
root = 0
LOOP 8 TIMES
dividend << 3
NOTE: discard the lower 24 bits to produce the partial dividend
partial_dividend = dividend >> 24
divisor = 12 * root * root + 6 * root + 1
IF partial_dividend >= divisor
partial_dividend -= divisor
root = root * 2 + 1
ELSE
root = root * 2
remainder = partial_dividend
RETURN (root, remainder)
Rather than calculate the square of the root, the assembly routine calculates a new divisor from the old divisor and the root so far. First, note that the first time through the loop, the divisor will be 1. Second, note that since the carry is clear when the comparison (and subtraction) is performed, the divisor stored in DIV is:
12 * root * root + 6 * root
rather than
12 * root * root + 6 * root + 1
Let r = old root
Thus,
old divisor = 12 * r*r + 6 * r
There are two cases:
Case #1: when the IF clause is true
new root = 2 * old root + 1 = 2 * r + 1
new divisor = 12 * new root * new root + 6 * new root
= 12 * (2 * r + 1) * (2 * r + 1) + 6 * (2 * r + 1)
= 48 * r*r + 48 * r + 12 + 12 * r +6
= 48 * r*r + 60 * r + 18
The new divisor is calculated as follows:
divisor = 2 * (divisor + 4 * new root)
= 2 * (12 * r*r + 6 * r + 4 * (2 * r + 1))
= 2 * (12 * r*r + 6 * r + 8 * r + 4)
= 24 * r*r + 28 * r + 8
divisor = 2 * (divisor + new root)
= 2 * (24 * r*r + 28 * r + 8 + 2 * r + 1)
= 48 * r + 60 * r + 18
Case #2: when the IF clause is false (i.e. when the ELSE clause is true)
new root = 2 * old root = 2 * r
new divisor = 12 * new root * new root + 6 * new root
= 12 * (2 * r) * (2 * r) + 6 * (2 * r)
= 48 * r*r + 12 * r
The new divisor is calculated as follows:
divisor = divisor - old root
= 12 * r*r + 6 * r - r
= 12 * r*r + 5 * r
divisor = divisor - 2 * old root
= 12 * r*r + 5 * r - 2 * r
= 12 * r*r + 3 * r
divisor = divisor * 4
= 48 * r*r + 12 * r
