Skip to content

Instantly share code, notes, and snippets.

@jac18281828
Created April 26, 2026 19:26
Show Gist options
  • Select an option

  • Save jac18281828/1f7d040cfa3747541b0599d5b1ba8531 to your computer and use it in GitHub Desktop.

Select an option

Save jac18281828/1f7d040cfa3747541b0599d5b1ba8531 to your computer and use it in GitHub Desktop.
ATAN MMIX
% =====================================================================
% test_atan.mms -- Test harness for the Atan library.
% =====================================================================
%
% Build: mmixal test_atan.mms
% Run: mmix test_atan
%
% Output looks like:
%
% atan tests:
% [00] PASS x=0000000000000000 got=0000000000000000 exp=0000000000000000
% [01] PASS x=3FF0000000000000 got=3FE921FB54442D18 exp=3FE921FB54442D18
% ...
% summary: 9/9 passed
%
% The test passes if either:
% (a) the bit pattern matches the expected reference exactly, OR
% (b) |computed - expected| <= 1e-12 (absolute tolerance).
%
% Reference values were computed in IEEE 754 double precision externally.
% =====================================================================
% Pull in the library. This places :Atan and its data into the program.
INCLUDE atan.mms
% --- Test data ----------------------------------------------------------
LOC Data_Segment+#1000 % keep clear of library data
TestBase GREG @
NumTests OCTA 9
% Inputs (IEEE 754 double bit patterns)
Inputs OCTA #0000000000000000 % 0.0
OCTA #3FF0000000000000 % 1.0
OCTA #BFF0000000000000 % -1.0
OCTA #3FE0000000000000 % 0.5
OCTA #4000000000000000 % 2.0
OCTA #C000000000000000 % -2.0
OCTA #4024000000000000 % 10.0
OCTA #7FF0000000000000 % +Inf
OCTA #FFF0000000000000 % -Inf
% Expected outputs (atan applied)
Expects OCTA #0000000000000000 % atan(0) = 0
OCTA #3FE921FB54442D18 % atan(1) = pi/4
OCTA #BFE921FB54442D18 % atan(-1) = -pi/4
OCTA #3FDDAC670561BB4F % atan(0.5)
OCTA #3FF1B6E192EBBE44 % atan(2)
OCTA #BFF1B6E192EBBE44 % atan(-2)
OCTA #3FF789BD2C160054 % atan(10)
OCTA #3FF921FB54442D18 % atan(+Inf) = pi/2
OCTA #BFF921FB54442D18 % atan(-Inf) = -pi/2
Tol OCTA #3D719799812DEA11 % 1e-12 in IEEE double
% --- Strings (NUL-terminated; each followed by an explicit 0 byte) ------
MsgHdr BYTE "atan tests:",#a,0
MsgOpen BYTE "[",0
MsgPass BYTE "] PASS x=",0
MsgFail BYTE "] FAIL x=",0
MsgGot BYTE " got=",0
MsgExp BYTE " exp=",0
MsgNL BYTE #a,0
MsgSum1 BYTE "summary: ",0
MsgSum2 BYTE "/",0
MsgSum3 BYTE " passed",#a,0
% Scratch buffer for formatting numbers (large enough for 16 hex digits + NUL)
Buf BYTE 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
BYTE 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
% --- Globals for main loop (survive PUSHJ calls) ------------------------
% These GREGs allocate global registers usable from any code below.
idx GREG 0 % current test index
total GREG 0 % NumTests cached
passes GREG 0 % pass counter
% --- Code ---------------------------------------------------------------
LOC #1000
GREG @
% PrintHex2: print a NUL-terminated 16-digit lowercase-hex of the value
% in $0 to stdout. Uses Buf.
% PUSHJ $X,PrintHex with x in $X+1.
PrintHex IS @
% $0 = value to print
GETA $1,Buf % $1 = buffer base
SETL $2,16 % digit count
SETL $3,0 % i
:PHLoop CMP $4,$3,$2
BNN $4,:PHEmit
SUBU $5,$2,$3 % position from left? we'll fill MSB-first
SUBU $5,$5,1 % index 15..0
% extract nibble (3-i)*4 from top: easier to extract from the value
% directly: we want digit at position $3 from MSB, i.e. shift by
% (15 - $3) * 4 to the right and mask 0xF.
SUBU $6,$2,$3
SUBU $6,$6,1 % $6 = 15 - i
SLU $6,$6,2 % *4 (bits)
SRU $7,$0,$6 % shifted value
AND $7,$7,#F % nibble
CMP $4,$7,10
BN $4,:PHDigit
% hex letter A..F (uppercase for consistency with constants)
ADDU $7,$7,'A'-10
JMP :PHStore
:PHDigit ADDU $7,$7,'0'
:PHStore STBU $7,$1,$3
ADDU $3,$3,1
JMP :PHLoop
:PHEmit SETL $4,0
STBU $4,$1,16 % NUL terminate
SET $255,$1
TRAP 0,Fputs,StdOut
POP 0,0
% PrintInt: print a small unsigned int in $0 as decimal, NUL-terminated.
PrintInt IS @
GETA $1,Buf
SETL $2,0 % length
SETL $3,$0 % n
BNZ $3,:PILoop
SETL $4,'0'
STBU $4,$1,0
SETL $4,0
STBU $4,$1,1
JMP :PIEmit
:PILoop BZ $3,:PIRev
DIV $4,$3,10
GET $5,rR % remainder
ADDU $5,$5,'0'
STBU $5,$1,$2
ADDU $2,$2,1
SET $3,$4
JMP :PILoop
:PIRev % digits stored low-to-high; reverse in place
SETL $3,0 % i
SUBU $4,$2,1 % j
:PIRevL CMP $5,$3,$4
BNN $5,:PIEnd
LDBU $6,$1,$3
LDBU $7,$1,$4
STBU $7,$1,$3
STBU $6,$1,$4
ADDU $3,$3,1
SUBU $4,$4,1
JMP :PIRevL
:PIEnd SETL $5,0
STBU $5,$1,$2 % NUL terminate
:PIEmit SET $255,$1
TRAP 0,Fputs,StdOut
POP 0,0
% PrintStr: print NUL-terminated string whose address is in $0.
PrintStr IS @
SET $255,$0
TRAP 0,Fputs,StdOut
POP 0,0
% =====================================================================
% Main
% =====================================================================
Main IS @
% Header
GETA $255,MsgHdr
TRAP 0,Fputs,StdOut
% Initialize loop state in globals
LDOU total,NumTests
SETL idx,0
SETL passes,0
:Loop CMP $0,idx,total
BNN $0,:Done
% Compute byte offset = idx * 8
SLU $1,idx,3
% Load x and expected
GETA $2,Inputs
LDOU $10,$2,$1 % x
GETA $2,Expects
LDOU $11,$2,$1 % expected
% Call Atan: argument goes in $X+1, result returns in $X.
% Use $20 as the PUSHJ target ("hole"); put x in $21 first.
SET $21,$10
PUSHJ $20,:Atan
SET $12,$20 % computed result
% Pass test if bit-equal OR |computed - expected| <= Tol.
CMP $0,$12,$11
BZ $0,:Pass
% Compute |diff| as float, then reinterpret as bits and mask off sign.
FSUB $0,$12,$11 % diff
LDOU $1,:Atan_AbsMask
AND $0,$0,$1 % |diff| bits
LDOU $1,Tol
FCMP $2,$0,$1
BNP $2,:Pass
JMP :Fail
:Pass ADDU passes,passes,1
% Print "[" idx "] PASS x=..."
GETA $0,MsgOpen
PUSHJ $30,PrintStr
SET $0,idx
PUSHJ $30,PrintInt
GETA $0,MsgPass
PUSHJ $30,PrintStr
SET $0,$10
PUSHJ $30,PrintHex
GETA $0,MsgGot
PUSHJ $30,PrintStr
SET $0,$12
PUSHJ $30,PrintHex
GETA $0,MsgExp
PUSHJ $30,PrintStr
SET $0,$11
PUSHJ $30,PrintHex
GETA $0,MsgNL
PUSHJ $30,PrintStr
JMP :Step
:Fail GETA $0,MsgOpen
PUSHJ $30,PrintStr
SET $0,idx
PUSHJ $30,PrintInt
GETA $0,MsgFail
PUSHJ $30,PrintStr
SET $0,$10
PUSHJ $30,PrintHex
GETA $0,MsgGot
PUSHJ $30,PrintStr
SET $0,$12
PUSHJ $30,PrintHex
GETA $0,MsgExp
PUSHJ $30,PrintStr
SET $0,$11
PUSHJ $30,PrintHex
GETA $0,MsgNL
PUSHJ $30,PrintStr
:Step ADDU idx,idx,1
JMP :Loop
:Done GETA $0,MsgSum1
PUSHJ $30,PrintStr
SET $0,passes
PUSHJ $30,PrintInt
GETA $0,MsgSum2
PUSHJ $30,PrintStr
SET $0,total
PUSHJ $30,PrintInt
GETA $0,MsgSum3
PUSHJ $30,PrintStr
TRAP 0,Halt,0
% =====================================================================
% atan.mms -- Arctangent for MMIX (IEEE 754 double precision)
% =====================================================================
%
% Public symbol: :Atan
%
% Calling convention (standard MMIX PUSHJ):
%
% SET $X+1,<x bits> ; or load x into $X+1
% PUSHJ $X,:Atan
% ; result in $X (IEEE 754 double, bit pattern)
%
% Equivalently: place input in the register one above the PUSHJ target,
% call, and the result lands in the PUSHJ target register.
%
% Behavior:
% atan(NaN) = NaN
% atan(+Inf) = +pi/2
% atan(-Inf) = -pi/2
% atan(+/-0) = +/-0
% atan(x) = arctangent in (-pi/2, pi/2) for finite x
%
% Algorithm:
% 1. If x is NaN, return NaN.
% 2. s := signbit(x); u := |x|.
% 3. If u is +Inf, return s ? -pi/2 : +pi/2.
% 4. If u > 1, set rec=1 and u := 1/u.
% 5. If u > tan(pi/8), set shf=1 and u := (u-1)/(u+1).
% 6. p := u + u^3 * Q(u^2), where Q is a degree-5 polynomial (Horner).
% 7. If shf: p := p + pi/4.
% 8. If rec: p := pi/2 - p.
% 9. If s: p := -p (toggle sign bit).
%
% Local register usage inside the routine: $0..$7. The PUSHJ register-stack
% mechanism guarantees these are fresh and do not clobber the caller.
% =====================================================================
% --- Constants (data segment) ---------------------------------------
LOC Data_Segment
GREG @
:Atan_PiOver2 OCTA #3FF921FB54442D18
:Atan_PiOver4 OCTA #3FE921FB54442D18
:Atan_TanPi8 OCTA #3FDA827999FCEF32
:Atan_One OCTA #3FF0000000000000
:Atan_AbsMask OCTA #7FFFFFFFFFFFFFFF
:Atan_SignBit OCTA #8000000000000000
:Atan_InfBits OCTA #7FF0000000000000
% Polynomial Q(t) = c3 + c5 t + c7 t^2 + c9 t^3 + c11 t^4 + c13 t^5
% used as: atan(u) ~= u + u^3 * Q(u^2) on |u| <= tan(pi/8)
% (After both range reductions, |u| stays well below tan(pi/8) so this
% Maclaurin-style truncation gives a few-ulp result; for production use
% one would substitute a true minimax polynomial of equal degree.)
:Atan_C3 OCTA #BFD5555555555555
:Atan_C5 OCTA #3FC9999999999999
:Atan_C7 OCTA #BFC2492492492492
:Atan_C9 OCTA #3FBC71C71C71C71C
:Atan_C11 OCTA #BFB745D1745D1746
:Atan_C13 OCTA #3FB3B13B13B13B14
% --- Code ---------------------------------------------------------
LOC #200
GREG @
:Atan IS @
% Step 1: NaN test (FUN sets result to 1 if unordered)
FUN $1,$0,$0
BNZ $1,:Atan_Ret % NaN -> return x unchanged
% Step 2: split into sign + magnitude
SRU $1,$0,63 % $1 = sign bit (0 or 1)
LDOU $2,:Atan_AbsMask
AND $0,$0,$2 % $0 = |x| (bit pattern of magnitude)
% Step 3: Infinity test (after abs, +Inf == 0x7FF0000000000000)
LDOU $3,:Atan_InfBits
CMP $4,$0,$3
BNZ $4,:Atan_NotInf
LDOU $0,:Atan_PiOver2
JMP :Atan_ApplySign
:Atan_NotInf IS @
% Reduction flags
SET $5,0 % rec
SET $6,0 % shf
% Step 4: reciprocal reduction if |x| > 1
LDOU $3,:Atan_One
FCMP $4,$0,$3
BNP $4,:Atan_NoRec
SET $5,1
FDIV $0,$3,$0 % u := 1/u
:Atan_NoRec IS @
% Step 5: shift reduction if u > tan(pi/8)
LDOU $3,:Atan_TanPi8
FCMP $4,$0,$3
BNP $4,:Atan_NoShf
SET $6,1
LDOU $7,:Atan_One
FSUB $3,$0,$7 % u - 1
FADD $4,$0,$7 % u + 1
FDIV $0,$3,$4 % u := (u-1)/(u+1)
:Atan_NoShf IS @
% Step 6: Horner evaluation of Q(t), t = u^2
FMUL $2,$0,$0 % t = u*u
LDOU $3,:Atan_C13
FMUL $3,$3,$2
LDOU $4,:Atan_C11
FADD $3,$3,$4
FMUL $3,$3,$2
LDOU $4,:Atan_C9
FADD $3,$3,$4
FMUL $3,$3,$2
LDOU $4,:Atan_C7
FADD $3,$3,$4
FMUL $3,$3,$2
LDOU $4,:Atan_C5
FADD $3,$3,$4
FMUL $3,$3,$2
LDOU $4,:Atan_C3
FADD $3,$3,$4 % $3 = Q(t)
FMUL $3,$3,$2 % Q*t = Q*u^2
FMUL $3,$3,$0 % Q*u^3
FADD $0,$0,$3 % p = u + Q*u^3
% Step 7: undo shift
BZ $6,:Atan_NoUnShf
LDOU $4,:Atan_PiOver4
FADD $0,$0,$4
:Atan_NoUnShf IS @
% Step 8: undo reciprocal
BZ $5,:Atan_NoUnRec
LDOU $4,:Atan_PiOver2
FSUB $0,$4,$0
:Atan_NoUnRec IS @
% Step 9: reapply sign by toggling the sign bit
:Atan_ApplySign IS @
BZ $1,:Atan_Ret
LDOU $4,:Atan_SignBit
XOR $0,$0,$4
:Atan_Ret IS @
POP 1,0
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment