Mercurial > hg > index.cgi

comparison src/fps.s @ 80:bb50ac9fdf37

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Checkpoint with very basic integer and floating point arithmetic, untested This commit has implementations for floating point add, subtract, multiply, and divide, along with 32 bit signed integer equivalents. These can probably be optimized and they are untested.
author William Astle <lost@l-w.ca>
date Sat, 07 Oct 2023 02:56:59 -0600
parents
children 9a4e2364a966
comparison
equal deleted inserted replaced
79:df86e6d64ce2 80:bb50ac9fdf37
1 *pragmapush list
2 *pragma list
3 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
4 ; Single precision floating point arithmetic package
5 ;
6 ; Floating point values are stored in 6 byte packets (single precision) organized as follows:
7 ; Offset Length Contents
8 ; 0 1 8 bit binary exponent with a bias of 128; 0 means the number is 0
9 ; 1 4 32 bit significand
10 ; 5 1 sign flag; zero for positive, 0xff for negative
11 ;
12 ; Binary operateions take pointers to their arguments in X and U. Y contains the result location. In all cases, it is
13 ; safe to specify either one of the arguments as the result location. Unary operations take their operand pointer in
14 ; X and their result pointer in Y. In all cases, there is an in place version of the unary operation that operates on
15 ; its input and only needs the single pointer in X.
16 ;
17 ; On the topic of optimization: these routines copy their operands to accumulators in the direct page. This saves one
18 ; cycle per instruction which has a nonzero offset. In addition, this means that the input operands can remain
19 ; unmodified by the actual operations. For instance, addition requires denormalizing one operand which would otherwise
20 ; have to be done in place.
21 ;
22 ; NOTE: the input pointers X and U may be clobbered.
23 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
24 ; Convert 32 bit unsigned integer at (X) to single precision floating point at (U)
25 fps_fromuint32 clr fpa0+fps.sign ; set result sign to positive
26 bra fps_fromint32a ; go do conversion
27 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
28 ; Convert 32 bit signed integer at (X) to single precision floating point at value accumulator in (Y)
29 fps_fromint32 ldb ,x ; set sign based on signed integer
30 sex
31 sta fpa0+fps.sign ; set result sign to the two's complement sign
32 bpl fps_fromint32a ; brif positive - no need to mess with bits
33 jsr int32_neg ; make the bits positive
34 fps_fromint32a ldb valtype_float ; set output value type to floating point
35 stb val.type,y
36 ldd ,x ; copy integer bits to fpa0
37 std fpa0+fps.sig
38 ldd 2,x
39 std fpa0+fps.sig+2
40 ldb #0xa0 ; put binary point to the right of the significand
41 stb fpa0+fps.exp
42 clr fpa0extra ; clear extra precision
43 jmp fps_add10 ; go normalize the result and return
44 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
45 ; Unary negation - negate (X) to (Y)
46 fps_neg ldd 2,x ; copy to output and keep exponent in A
47 std 2,y
48 ldd 4,x
49 std 4,y
50 ldd ,x
51 std ,y
52 tsta ; is the number zero?
53 beq fps_neg0 ; brif so - do nothing
54 com fps.sign,y ; flip the sign
55 fps_neg0 rts
56 fps_negx lda fps.exp,x ; is the number zero?
57 beq fps_negx0 ; brif so - do nothing
58 com fps.sign,x ; flip the sign
59 fps_negx0 rts
60 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
61 ; Copy (X) to fpa0 and (U) to fpa1
62 fps_copyinputs ldd ,x ; copy (X) to fpa0
63 std fpa0
64 ldd 2,x
65 std fpa0+2
66 ldd 4,x
67 std fpa0+4
68 ldd ,u ; copy (U) to fpa0
69 std fpa1
70 ldd 2,u
71 std fpa1+2
72 ldd 4,u
73 std fpa1+4
74 rts
75 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
76 ; Subtraction (X) - (U) to (Y)
77 fps_sub bsr fps_copyinputs ; copy input operands
78 com fpa1+fps.sign ; negate the subtrahend (don't need to handle zero here)
79 bra fps_add0 ; go handle it as addition
80 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
81 ; Addition (X) + (U) to (Y)
82 fps_add bsr fps_copyinputs ; copy input operands
83 fps_add0 lda fpa1+fps.exp ; is the second operand 0?
84 bne fps_add1 ; brif not
85 ldd fpa0 ; copy first operand to output
86 std ,y
87 ldd fpa0+2
88 std 2,y
89 ldd fpa0+4
90 std 4,y
91 rts
92 fps_add1 lda fpa0+fps.exp ; get exponent of first operand
93 bne fps_add2 ; brif not zero
94 ldd fpa1 ; copy second operand to output
95 std ,y
96 ldd fpa1+2
97 std 2,y
98 ldd fpa1+4
99 std 4,y
100 rts
101 fps_add2 clr fpa0extra ; clear extra precision bits
102 lda fpa0+fps.exp ; get first operand exponent
103 suba fpa1+fps.exp ; get difference in exponents
104 beq fps_add8 ; brif we don't need to denormalize - they're the same
105 blo fps_add3 ; brif we need to denormalize the first operand
106 ldx #fpa1 ; point to second operand to denormalize
107 bra fps_add4 ; go denormalize
108 fps_add3 ldb fpa1+fps.exp ; get exponent of second operand (result exponent)
109 stb fpa0+fps.exp ; set result exponent
110 ldx #fpa0 ; point to first operand to denormalize
111 nega ; get number of bits to shift as positive number
112 fps_add4 suba #8 ; is there 8 bits left?
113 blo fps_add5 ; brif not
114 ldb fps.sig+3,x ; shift significand 8 bits right
115 stb fpa0extra ; save extra precision bits
116 ldb fps.sig+2,x
117 stb fps.sig+3,x
118 ldb fps.sig+1,x
119 stb fps.sig+2,x
120 ldb fps.sig,x
121 stb fps.sig+1,x
122 clr fps.sig,x
123 bra fps_add4 ; see if we have another byte to shift
124 fps_add5 adda #8 ; adjust for extra subtract above
125 bra fps_add7 ; do the bit shifting
126 fps_add6 lsr fps.sig,x ; shift one bit right
127 ror fps.sig+1,x
128 ror fps.sig+2,x
129 ror fps.sig+3,x
130 ror fpa0extra
131 deca ; done all of the bits?
132 fps_add7 bne fps_add6 ; brif not
133 fps_add8 lda fpa0+fps.sign ; compare the signs of the operands
134 eora fpa1+fps.sign ; set A if signs differ
135 bne fps_add9 ; brif they differ - do subtraction
136 ldd fpa0+fps.sig+2 ; add low word of significands
137 addd fpa1+fps.sig+2
138 std fpa0+fps.sig+2
139 ldd fpa0+fps.sig ; and the high word
140 adcb fpa1+fps.sig+1
141 adca fpa1+fps.sig
142 std fpa0+fps.sig
143 bcc fps_add9 ; brif no carry
144 ror fpa0+fps.sig ; shift carry into significand
145 ror fpa0+fps.sig+1
146 ror fpa0+fps.sig+2
147 ror fpa0+fps.sig+3
148 ror fpa0extra ; and the extra bits (for rounding)
149 inc fpa0+fps.exp ; bump exponent to account for bit shift
150 bne fps_add14 ; go check for round-off if not overflow
151 OVERROR ldb #err_ov ; raise overflow
152 jmp ERROR
153 fps_add9 ldd fpa0+fps.sig+2 ; subtract low word
154 subd fpa1+fps.sig+2
155 std fpa0+fps.sig+2
156 ldd fpa0+fps.sig ; and now the high word
157 sbcb fpa1+fps.sig+1
158 sbca fpa1+fps.sig
159 std fpa0+fps.sig
160 bcc fps_add10 ; brif we didn't carry
161 com fpa0+fps.sign ; flip result sign - other number was bigger
162 com fpa0+fps.sig+3 ; negate two's complement result to be just the magnitude
163 com fpa0+fps.sig+2
164 com fpa0+fps.sig+1
165 com fpa0+fps.sig
166 ldx fpa0+fps.sig+2 ; add 1 to complete negation
167 leax 1,x
168 stx fpa0+fps.sig+2
169 bne fps_add10 ; brif carry doesn't propagate
170 ldx fpa0+fps.sig ; propagate carry
171 leax 1,x
172 stx fpa0+fps.sig ; NOTE: this cannot carry because magnitude got smaller
173 fps_add10 clra ; initialize exponent offset
174 fps_add11 ldb fpa0+fps.sig ; do we have nonzero bits in high byte of significand?
175 bne fps_add13 ; brif so
176 ldb fpa0+fps.sig+1 ; shift left 8 bits
177 stb fpa0+fps.sig
178 ldb fpa0+fps.sig+2
179 stb fpa0+fps.sig+1
180 ldb fpa0+fps.sig+3
181 stb fpa0+fps.sig+2
182 ldb fpa0extra ; and extra precision bits
183 stb fpa0+fps.sig+3
184 clr fpa0extra
185 addb #8 ; account for number of bits shifted
186 cmpb #40 ; done 40 bits?
187 blo fps_add11 ; brif not - see if we have more bits to shift
188 clr fpa0+fps.exp ; number underflowed to zero - set it so
189 clr fpa0+fps.sign
190 clr fpa0+fps.sig
191 clr fpa0+fps.sig+1
192 clr fpa0+fps.sig+2
193 clr fpa0+fps.sig+3
194 bra fps_add16 ; go return result
195 fps_add12 inca ; account for a bit shift
196 lsl fpa0extra ; shift significand and extra bits left
197 rol fpa0+fps.sig+3
198 rol fpa0+fps.sig+2
199 rol fpa0+fps.sig+1
200 rol fpa0+fps.sig
201 fps_add13 bpl fps_add12 ; brif we haven't normalized yet
202 fps_add14 ldb fpa0extra ; do we need to round?
203 bpl fps_add16 ; brif not
204 ldx fpa0+fps.sig+2 ; add one to significand
205 leax 1,x
206 stx fpa0+fps.sig+2
207 bne fps_add16 ; brif no carry
208 ldx fpa0+fps.sig ; bump the upper word of significand
209 leax 1,x
210 stx fpa0+fps.sig
211 bne fps_add16
212 fps_add15 inc fpa0+fps.exp ; bump exponent
213 beq OVERROR ; brif it overflowed
214 lda #0x80 ; set high bit of significand (rest is zero)
215 sta fpa0+fps.sig
216 fps_add16 ldd fpa0 ; copy result to destination
217 std ,y
218 ldd fpa0+2
219 std 2,y
220 ldd fpa0+4
221 std 4,y
222 rts
223 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
224 ; Single precision multiplication (X) × (U) to (Y)
225 fps_mul lda fps.exp,x ; is first operand zero?
226 beq fps_mul0 ; brif so - return zero
227 lda fps.exp,u ; is second operand zero?
228 bne fps_mul1 ; brif not - have to do the multiply
229 fps_mul0 ldd zero ; return zero result
230 std ,y
231 std 2,y
232 std 4,y
233 rts
234 fps_mul1 jsr fps_copyinputs ; copy input operands
235 lda fpa0+fps.sign ; calculate sign of result - xor of the two signs
236 eora fpa1+fps.sign
237 sta fpa0+fps.sign ; save result sign
238 lda fpa0+fps.exp ; get exponent of first value
239 adda fpa1+fps.exp ; calculate new exponent; also cancels the bias
240 rora ; set V if C and N differ
241 rola
242 bvc fps_mul3 ; brif maybe an overflow
243 adda #0x80 ; add back the bias
244 sta fpa0+fps.sign ; save new sign
245 beq fps_mul0 ; brif we underflowed - zero out sign
246 ; This does a shift-and-add multiplication algorithm. This is slower than an equivalent using MUL but is smaller.
247 ; The high order bytes will be left in fpa0 with the low order bytes in fpa0extra and fpa0extra[1-3]. The low order
248 ; bytes are kept for the odd case where extra precision is useful. Uses fpa0extra[4-7] as temporaries.
249 fps_mul2 ldd zero ; zero out temporary bytes
250 std fpa0extra4
251 std fpa0extra6
252 ldb fpa0+fps.sig+3 ; multiply by low byte of fpa0
253 bsr fps_mul4
254 ldb fpa0extra ; move low bites
255 stb fpa0extra3
256 ldb fpa0+fps.sig+2 ; multiply by next higher byte
257 bsr fps_mul4
258 ldb fpa0extra ; move low bits
259 stb fpa0extra2
260 ldb fpa0+fps.sig+1 ; and again for the next higher byte
261 bsr fps_mul4
262 ldb fpa0extra
263 stb fpa0extra1
264 ldb fpa0+fps.sig ; and the high order byte
265 bsr fps_mul4
266 ldd fpa0extra4 ; copy high order product bits to result
267 std fpa0+fps.sig
268 ldd fpa0extra6
269 std fpa0+fps.sig+2
270 jmp fps_add10 ; go normalize the result and return
271 fps_mul3 bpl fps_mul0 ; brif we underflowed - return 0
272 jmp OVERROR ; raise overflow
273 fps_mul4 bne fps_mul5 ; brif not multiply by zero
274 lda fpa0+fps.sig+3 ; shift 8 bits right
275 sta fpa0extra
276 lda fpa0+fps.sig+2
277 sta fpa0+fps.sig+3
278 lda fpa0+fps.sig+1
279 sta fpa0+fps.sig+2
280 lda fpa0+fps.sig
281 sta fpa0+fps.sig+1
282 clr fpa0+fps.sig
283 fps_mul8 rts
284 fps_mul5 coma ; set C
285 fps_mul6 lda fpa0extra4 ; get high byte of result bytes
286 rorb ; is multiplier bit set?
287 beq fps_mul8 ; brif 8 shifts done (C set above makes sure of that)
288 bcc fps_mul7 ; brif bit not set - don't do addition
289 lda fpa0extra7 ; add multiplier (fpa1) to result
290 adda fpa1+fps.sig+3
291 sta fpa0extra7
292 lda fpa0extra6
293 adca fpa1+fps.sig+2
294 sta fpa0extra6
295 lda fpa0extra5
296 adca fpa1+fps.sig+1
297 sta fpa0extra5
298 lda fpa0extra4
299 adca fpa1+fps.sig
300 fps_mul7 rora ; rotate carry in (from add or 0)
301 sta fpa0extra4
302 ror fpa0extra5
303 ror fpa0extra6
304 ror fpa0extra7
305 ror fpa0extra ; and into the extra precision bits
306 clra ; clear carry - so shift above will terminate
307 bra fps_mul6 ; go do another bit
308 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
309 ; Single precision division (X) ÷ (U) -> (Y)
310 ;
311 ; This is basically the same algorithm used in the Color Basic ROM
312 fps_div lda fps.exp,u ; is divisor 0?
313 bne fps_div0
314 DIV0ERROR ldb #err_div0 ; raise division by zero
315 jmp ERROR
316 fps_div0 lda fps.exp,x ; is dividend 0?
317 lbeq fps_mul0 ; brif so - return 0
318 jsr fps_copyinputs ; copy input values
319 lda fpa0+fps.sign ; calculate result sign - xor of operand signs
320 eora fpa1+fps.sign
321 sta fpa0+fps.sign ; save result sign
322 lda fpa1+fps.exp ; get divisor exponent
323 nega ; negate for subtraction
324 adda fpa0+fps.exp ; subtract it from dividend exponent
325 rora ; set V if C and N differ
326 rola
327 bvc fps_mul3 ; brif overflow or underflow
328 adda #0x80 ; add back the bias
329 sta fpa0+fps.sign ; save new sign
330 lbeq fps_mul0 ; brif we underflowed - zero out sign
331 inca ; bump exponent - why?? related to the bias stuff above?
332 lbeq OVERROR ; brif it overflows
333 sta fpa0+fps.exp ; save result exponent
334 ldx #fpa0extra4 ; point to temporary storage bytes for quotient
335 ldb #4 ; counter for 4 significand bytes and one extra partial byte
336 stb fpa0extra1 ; save counter since we need both accumulators
337 ldb #1 ; shift counter flag and quotient byte
338 fps_div1 lda fpa1+fps.sig ; set C if fpa0 significand <= fpa1 significand
339 cmpa fpa0+fps.sig
340 bne fps_div2
341 lda fpa1+fps.sig+1
342 cmpa fpa0+fps.sig+1
343 bne fps_div2
344 lda fpa1+fps.sig+2
345 cmpa fpa0+fps.sig+2
346 bne fps_div2
347 lda fpa1+fps.sig+3
348 cmpa fpa0+fps.sig+3
349 bne fps_div2
350 coma ; set C if values are the same
351 fps_div2 tfr cc,a ; save C for later, C clear if fpa1 > fpa0
352 rolb ; shift carry into quotient
353 bcc fps_div3 ; brif carry clear - we haven't done 8 bits yet
354 stb ,x+ ; save quotient byte after a full set of bits
355 dec fpa0extra1 ; have we done all the bytes?
356 bmi fps_div7 ; brif all bytes plus extra precision - done all
357 beq fps_div6 ; brif all main sigificand bytes - do a couple extra bits
358 ldb #1 ; reset the bit counter flag
359 fps_div3 tfr cc,a ; get back original carry from compare
360 bcs fps_div5 ; brif it "went"
361 fps_div4 lsl fpa0+fps.sig+3 ; shift dividend left
362 rol fpa0+fps.sig+2
363 rol fpa0+fps.sig+1
364 rol fpa0+fps.sig
365 bcs fps_div2 ; brif it carries - next bit "goes"
366 bmi fps_div1 ; check magnitudes of next bit
367 bra fps_div2 ; carry clear - check another bit
368 fps_div5 lda fpa0+fps.sig+3 ; subtract divisor from dividend bits
369 suba fpa1+fps.sig+3
370 sta fpa0+fps.sig+3
371 lda fpa0+fps.sig+2
372 sbca fpa1+fps.sig+2
373 sta fpa0+fps.sig+2
374 lda fpa0+fps.sig+1
375 sbca fpa1+fps.sig+1
376 sta fpa0+fps.sig+1
377 lda fpa0+fps.sig
378 sbca fpa1+fps.sig
379 sta fpa0+fps.sig
380 bra fps_div4 ; now do the bit shift to line things up
381 fps_div6 ldb #0x40 ; only do two bits of extra precision byte
382 bra fps_div2 ; go handle these bitsd
383 fps_div7 rorb ; get extra quotient bits to bit 7,6 and bit 5 set
384 rorb
385 rorb
386 stb fpa0extra ; save extra precision bits
387 ldd fpa0extra4 ; copy result bits to fpa0
388 std fpa0+fps.sig
389 ldd fpa0extra6
390 std fpa0+fps.sig+2
391 jmp fps_add10 ; go normalize the result
392 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
393 ; Pack single precision number at (X) to (U)
394 fps_pack lda fps.sign,x ; get sign of number (will be 0 for 0)
395 ora #0x7f ; make sure low bits are set for merging
396 anda fps.sig,x ; merge with high bits of significand
397 ldb fps.sig+1,x ; get upper mid bits of significand
398 std fps.sig,u
399 ldd fps.sig+2,x
400 std fps.sig+2,u
401 lda fps.exp,x
402 sta fps.exp,u
403 rts
404 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
405 ; Unpack single precision number at (X) to (U)
406 fps_unpack lda fps.exp,x ; get exponent of value
407 beq fps_unpack0 ; brif value is zero
408 sta fps.exp,u
409 ldb fps.sig,x ; get high byte of significand
410 sex ; make sign value in A
411 sta fps.sign,u ; set sign in result
412 ldd fps.sig,x ; get high word of sifnificand
413 ora #0x80 ; make sure high bit is set
414 std fps.sig,u ; save high word in result
415 ldd fps.sig+2,x ; copy middle bytes of significand
416 std fps.sig+2,u
417 rts
418 fps_unpack0 sta fps.exp,u ; zero out destination
419 sta fps.sig,u
420 sta fps.sig+1,u
421 sta fps.sig+2,u
422 sta fps.sig+3,u
423 sta fps.sign,u
424 rts
425 *pragmapop list