4 REM Table-based trig functions
6 REM © 1995-1998 Straylight
9 REM ----- Licensing note ----------------------------------------------------
11 REM This file is part of Straylight's Sapphire library.
13 REM Sapphire is free software; you can redistribute it and/or modify
14 REM it under the terms of the GNU General Public License as published by
15 REM the Free Software Foundation; either version 2, or (at your option)
18 REM Sapphire is distributed in the hope that it will be useful,
19 REM but WITHOUT ANY WARRANTY; without even the implied warranty of
20 REM MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 REM GNU General Public License for more details.
23 REM You should have received a copy of the GNU General Public License
24 REM along with Sapphire. If not, write to the Free Software Foundation,
25 REM 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
27 REM --- How the table works ---
29 REM We represent a number in 16.16 fixed point. We set up a table of
30 REM arctan values for values of x between 0 and 1. All the others may be
31 REM calculated from this, hopefully. Certainly, this is enough for
32 REM calculating polar angles, which is the really interesting bit.
34 REM For niceness purposes, the values are actually stored in degrees, rather
35 REM than radians. Division by 360 isn't a major difficulty, so it's not
36 REM hard to convert back again.
38 REM --- Set up some variables ---
40 ON ERROR ERROR 0,REPORT$+" ["+STR$(ERL)+"]"
41 bas%=TRUE :REM Generate BAS code
44 REM --- BAS library ---
51 REM --- A little bit of test code ---
54 PROCbas_aofInit(&10000)
63 REM --- General parameters ---
65 scalebits% = 16 :REM Avoiding tweaking this: it's part of the interface
66 scale% = 1 << scalebits%
68 REM --- Parameters for the tan table --
70 atn_tblbits% = 7 :REM TWEAK ME!
71 atn_entries% = 1 << atn_tblbits%
72 atn_shift% = scalebits% - atn_tblbits%
73 atn_revshift% = 32 - atn_shift%
74 atn_interpoff% = 1 << (atn_shift% - 1)
76 REM --- Parameters for the sin table ---
78 sin_tblbits% = 6 :REM TWEAK ME! sin_tblbits may be at most 8
79 sin_entries% = 1 << sin_tblbits%
80 sin_shift% = 30 - sin_tblbits%
81 sin_modmask% = (sin_entries% - 1) << sin_shift%
82 sin_interpoff% = 1 << (sin_shift% - 1)
84 REM --- Start the assembly ---
86 FOR o%=o_base% TO o_lim% STEP 2 :REM Standard two-pass assembly loop
87 IF bas%=0 THEN P%=code% :REM Start generating code in buffer
93 FNarea("Sapphire$$Code","CODE,READONLY")
95 FNimport("div_unsigned")
96 FNimport("div_u64x32")
100 [ OPT o% :REM Start assembly
104 ; On entry R0 == x, in 16.16 fixed point form
106 ; On exit R0 == arctan x, in degrees, in 16.16 fixed point
108 ; Use Calculates arctan x, hopefully fairly swiftly. The
109 ; accuracy of the result is open to doubt, although
110 ; it's usually good to about 3 significant figures.
111 ; It uses a small lookup table and linear interpolation
112 ; to calculate the result.
114 .arctan STMFD R13!,{R1-R3,R14} ;Save some work registers
116 ; --- Set some initial things up ---
118 MOV R3,#0 ;Clear out my flags register
120 ; --- Now sort out the argument ---
122 CMP R0,#0 ;Is the argument negative?
123 ORRLT R3,R3,#1 ;Yes -- set the `negative' bit
124 RSBLT R0,R0,#0 ;And make the argument positive
126 CMP R0,#scale% ;Is the argument bigger than 1?
127 BLE atn_skip ;No; skip ahead a bit
129 ORR R3,R3,#2 ;Set the `reciprocal' bit
131 IF scalebits% < 16 THEN
133 MOV R1,R0 ;Stand by to do some division
134 MOV R0,#1 << (scalebits% * 2)
139 MOV R2,R0 ;Stand by to do some division
141 MOV R1,#1 << (scalebits% * 2 - 32)
147 ; --- Look up the arctan value ---
149 .atn_skip ADR R14,atn_table ;Load the address of the table
150 MOV R1,R0,LSL #atn_revshift% ;Get bottom bits in R1
151 MOV R1,R1,LSR #atn_revshift% ;Shift back to bit 0
152 MOV R0,R0,LSR #atn_shift% ;Get the table index
153 ADD R14,R14,R0,LSL #2 ;Get the address of the item
154 LDR R2,[R14,#0] ;Locate the lower entry and load
155 LDR R0,[R14,#4] ;Locate the upper entry and load
156 SUB R0,R0,R2 ;Get the difference between them
157 MUL R0,R1,R0 ;Multiply this by the x offset
158 ADD R0,R0,#atn_interpoff% ;Round to nearest
159 ADD R0,R2,R0,LSR #atn_shift% ;Do interpolation
161 ; --- Now postprocess this result nicely ---
163 TST R3,#2 ;Did I have to reciprocal it?
164 RSBNE R0,R0,#90*scale% ;Yes -- work out the correct one
165 TST R3,#1 ;Did I have to negate the value?
166 RSBNE R0,R0,#0 ;Yes -- negate the result
168 LDMFD R13!,{R1-R3,PC}^ ;Return to caller
172 ; On entry R0 == x coordinate
175 ; On exit R0 == angle in degrees, in 16.16 form
177 ; Use Calculates the angle a vector makes with the +ve x axis.
178 ; The angle is given in degrees, rather than radians,
179 ; although this isn't really overly significant, it just
180 ; makes it slightly easier to work with, because it's
183 ; This routine uses the arctan table and linear
184 ; interpolation, so it's fairly quick, but the accuracy
185 ; of its results is questionable.
187 .pol STMFD R13!,{R1-R3,R14} ;Save some registers away
188 MOV R3,#0 ;Clear out my status flags
190 ; --- Preprocess the arguments ---
192 CMP R0,#0 ;Is x less than zero?
193 RSBLT R0,R0,#0 ;Yes -- make it positive
194 ORRLT R3,R3,#1 ;And remember that we did this
196 CMP R1,#0 ;Is y less than zero?
197 RSBLT R1,R1,#0 ;Yes -- make it positive
198 ORRLT R3,R3,#2 ;And remember we did this
200 CMP R0,R1 ;Is x bigger than y?
201 EORGT R0,R0,R1 ;Yes -- then swap them round
204 ORRGT R3,R3,#4 ;And remember we did this
206 ; --- Handle silly case where point is at (0,0) ---
208 CMP R1,#0 ;Are we about to divide by 0?
209 LDMEQFD R13!,{R1-R3,PC}^ ;Yes -- return now -- R0 is 0
211 ; --- Now work out arctan(R0/R1) ---
213 MOV R2,R1 ;Get the divisor
214 MOV R1,R0,LSR #32 - scalebits% ;Sort out dividend
215 MOV R0,R0,LSL #scalebits%
216 BL div_u64x32 ;Work out the ratio nicely
218 ADR R14,atn_table ;Load the address of the table
219 MOV R1,R0,LSL #atn_revshift% ;Get bottom bits in R1
220 MOV R1,R1,LSR #atn_revshift% ;Shift back to bit 0
221 MOV R0,R0,LSR #atn_shift% ;Get correct table index
222 ADD R14,R14,R0,LSL #2 ;Get the address of the item
223 LDR R2,[R14,#0] ;Locate the lower entry and load
224 LDR R0,[R14,#4] ;Locate the upper entry and load
225 SUB R0,R0,R2 ;Get the difference between them
226 MUL R0,R1,R0 ;Multiply this by the x offset
227 ADD R0,R0,#atn_interpoff% ;Just round to nearest
228 ADD R0,R2,R0,LSR #atn_shift% ;And add to the base
230 ; --- arctan is now in R0 -- process it for angle ---
232 TST R3,#4 ;Did we swap x and y round?
233 RSBEQ R0,R0,#90*scale% ;No -- compensate for this
234 TST R3,#1 ;Was x negative?
235 RSBNE R0,R0,#180*scale% ;Yes -- push it to left side
236 TST R3,#2 ;Was y negative?
237 RSBNE R0,R0,#360*scale% ;Yes -- push it to bottom
239 LDMFD R13!,{R1-R3,PC}^ ;Return to caller
241 ; --- arctan table ---
246 FOR i%=0 TO atn_entries% :REM Go through each index value
247 a=DEG ATN(i%/atn_entries%) :REM Calculate the arctan value
248 [opt o%:dcd a*scale%:] :REM Store it in the table nicely
255 ; On entry R0 == angle in degrees, in 16.16 form
257 ; On exit R0 == sin of angle, in 16.16 form
259 ; Use Calculates a sin of an angle with a degree of swiftness
260 ; and a lot less accuracy.
262 .cos RSB R0,R0,#90*scale%
264 .sin STMFD R13!,{R1-R3,R14} ;Save a job load of registers
266 ; --- How it all works ---
268 ; Having angles in degrees imposes inconvenient units
269 ; on use, forcing us to do all sorts of horrible things
270 ; with division by 360.
272 ; Hence, we translate the degree units into nice things
273 ; which represent a right angle as a power of two,
276 MOV R1,R0,LSR #26 ;Keep the top half somewhere
277 MOV R0,R0,LSL #6 ;Shift up so as not to lose
278 MOV R2,#45 ;Divide by 360 to translate units
279 BL div_u64x32 ;Now we have Internal Angle Units
280 MOV R0,R0,LSL #32 - 9 - scalebits%
281 ;Shift extra revolutions off
283 ; --- Now fix it into the first quadrant ---
285 AND R3,R0,#&C0000000 ;Get quadrant number in R3
286 BIC R0,R0,#&C0000000 ;And angle within quadrant here
287 TST R3,#&40000000 ;Is it on the left hand side?
288 RSBNE R0,R0,#&40000000 ;Yes; reflect across the circle
290 ; --- Look up the angle in the table ---
292 BIC R2,R0,#sin_modmask% ;Get remainder for interp
293 MOV R0,R0,LSR #sin_shift% ;Convert to a nice index
294 ADR R14,sin_table ;Find the sin table
295 ADD R14,R14,R0,LSL #2 ;Find the item to read
296 LDR R1,[R14,#0] ;Load the lower entry
297 LDR R0,[R14,#4] ;Load the upper entry too
298 SUB R0,R0,R1 ;Get the difference out
299 MUL R0,R2,R0 ;Multiply this by the x offset
300 ADD R0,R0,#sin_interpoff% ;Add on for rounding
301 ADD R0,R1,R0,LSR #sin_shift% ;Do linear interp
303 ; --- Now postprocess as usual ---
305 TST R3,#&80000000 ;Are we down at the bottom there?
306 RSBNE R0,R0,#0 ;Yes -- negate the sin value
308 LDMFD R13!,{R1-R3,PC}^
315 FOR i%=0 TO sin_entries% :REM Go through each index value
316 index=RAD(i%*90/sin_entries%) :REM Convert to a value in radians
317 [opt o%:dcd SIN(index)*scale%:] :REM Fill in the sin value nicely
322 FNexportAs("arctan","fxp_atan")
323 FNexportAs("pol","fxp_pol")
324 FNexportAs("sin","fxp_sin")
325 FNexportAs("cos","fxp_cos")
332 ; On entry R0 == dividend
335 ; On exit R0 == quotient
338 ; Use Divides on number by another in a vaguely slow way
340 .divide CMP R1,#0 ;Is this a division by zero?
341 MOVEQ PC,#0 ;Yes -- cause an exception
342 STMFD R13!,{R2,R14} ;Save some working registers
344 ; --- Multiply up the divisor ---
346 MOV R14,R1 ;Take a copy of the divisor
347 CMP R14,R0,LSR #1 ;Is this almost big enough?
348 .a MOVLS R14,R14,LSL #1 ;No -- shift the divisor up
349 CMP R14,R0,LSR #1 ;Is this almost big enough?
350 BLS a ;No -- loop round again
352 ; --- Now we can start dividing properly ---
354 MOV R2,#0 ;Start the quotient at 0
355 .a CMP R0,R14 ;Can we divide here?
356 SUBCS R0,R0,R14 ;Yes -- subtract the divisor
357 ADC R2,R2,R2 ;Add with carry nicely
358 MOV R14,R14,LSR #1 ;Shift divisor back down
359 CMP R14,R1 ;Have we finished this loop?
360 BCS a ;And go round again
362 ; --- Set up the return values ---
364 MOV R1,R0 ;Return remainder in R1
365 MOV R0,R2 ;Return quotient in R2
366 LDMFD R13!,{R2,PC}^ ;Return to caller
370 ; On entry R0,R1 = dividend
373 ; On exit R0 = quotient
376 ; Use Does long division.
378 .div_u64x32 STMFD R13!,{R14}
398 .testpol stmfd r13!,{r14}
403 .testatn stmfd r13!,{r14}
408 .testsin stmfd r13!,{r14}
418 NEXT :REM End current pass
424 REM --- Test the divide routine ---
426 REPEAT :REM Go into a nice loop
427 INPUT a$, A, B :REM Read the divide arguments
431 WHEN "sin": CALL testsin: r=SINRAD(A)
432 WHEN "atn": CALL testatn: r=DEGATN(A)
433 WHEN "pol": CALL testpol: r=DEGFNpol(A, B)
434 OTHERWISE: PRINT "wtf?"
436 PRINT !result/scale%'r
437 UNTIL FALSE :REM Keep on looping until bored
445 IF x=0 THEN a=PI/2 ELSE a=ATN(ABS(y/x))
447 IF x<0 THEN a=PI+a ELSE a=2*PI-a