|
57
|
1 ;FBASIC, Floating point BASIC.
|
|
|
2
|
|
|
3 ;This is not a BASIC interpreter, but a simple RPN calculator
|
|
|
4 ;to test the floating point routines. As such it is not a finished
|
|
|
5 ;application.
|
|
|
6 ;Written in 1996 by Lennart Benschoo.
|
|
|
7 ;
|
|
|
8 ;2014-07-26: Added welcome message, a few more comments.
|
|
|
9
|
|
|
10 ;Configuration info, change this for different apps.
|
|
|
11 ROM equ 0 ;Flag to indicate that BASIC is in ROM
|
|
|
12 ROMSTART equ $8000 ;First ROM address.
|
|
|
13 RAMSTART equ $400 ;First RAM address.
|
|
|
14 RAMTOP equ $8000 ;Last RAM address +1.
|
|
|
15
|
|
|
16 PROGORG equ ROM*ROMSTART+(1-ROM)*RAMSTART
|
|
|
17
|
|
|
18 ;First the O.S. vectors in the zero page.
|
|
|
19 org $0000
|
|
|
20 * First the I/O routine vectors.
|
|
|
21 getchar rmb 3 ;Jump to getchar routine.
|
|
|
22 putchar rmb 3 ;Jump to putchar routine.
|
|
|
23 getline rmb 3 ;Jump to getline routine.
|
|
|
24 putline rmb 3 ;Jump to putline routine.
|
|
|
25 putcr rmb 3 ;Jump to putcr routine.
|
|
|
26 getpoll rmb 3 ;Jump to getpoll routine.
|
|
|
27 xopenin rmb 3 ;Jump to xopenin routine.
|
|
|
28 xopenout rmb 3 ;Jump to xopenout routine.
|
|
|
29 xabortin rmb 3 ;Jump to xabortin routine.
|
|
|
30 xclosein rmb 3 ;Jump to xclosein routine.
|
|
|
31 xcloseout rmb 3 ;Jump to xcloseout routine.
|
|
|
32 delay rmb 3 ;Jump to delay routine.
|
|
|
33
|
|
|
34 timer equ *+6 ;3-byte timer.
|
|
|
35 linebuf equ $200 ;Input line buffer
|
|
|
36 xerrvec equ $280+6*3 ;Error vector.
|
|
|
37
|
|
|
38 * Now BASIC's own zero-page allocations.
|
|
|
39 org $40
|
|
|
40 startprog rmb 2 ;Start of BASIC program.
|
|
|
41 endprog rmb 2 ;End of BASIC program.
|
|
|
42 endvar rmb 2 ;End of variable area.
|
|
|
43 fpsp rmb 2 ;Floating point stack pointer (grows up).
|
|
|
44 endmem rmb 2
|
|
|
45
|
|
|
46 intbuf rmb 4 ;Buffer to store integer.
|
|
|
47 intbuf2 rmb 4
|
|
|
48 bcdbuf rmb 5 ;Buffer for BCD conversion.
|
|
|
49
|
|
|
50 endstr rmb 2 ;End address of string.
|
|
|
51 dpl rmb 1 ;Decimal point location.
|
|
|
52
|
|
|
53 * The BASIC interpreter starts here.
|
|
|
54 org PROGORG
|
|
|
55 cold jmp docold
|
|
|
56 warm bra noboot
|
|
|
57
|
|
|
58 * Cold start routine.
|
|
|
59 docold ldx #RAMTOP
|
|
|
60 stx endmem
|
|
|
61 tfr x,s
|
|
|
62 ldu #FREEMEM
|
|
|
63 stu startprog
|
|
|
64 clr ,u+
|
|
|
65 clr ,u+
|
|
|
66 stu endprog
|
|
|
67 ldx PROGEND
|
|
|
68 leax 1,x
|
|
|
69 beq noboot ;Test for autoboot program.
|
|
|
70 ldx #PROGEND
|
|
|
71 stx startprog
|
|
|
72 jmp dorun
|
|
|
73 noboot jsr doclear
|
|
|
74 ;; Print a welcome message first.
|
|
|
75 ldx #nbmesg
|
|
|
76 ldb #nbmend-nbmesg
|
|
|
77 jsr putline
|
|
|
78 jsr putcr
|
|
|
79 ldd #$4000
|
|
|
80 std fpsp
|
|
|
81 ldu fpsp
|
|
|
82 ;; Main loop. This is a simple RPN calculator that treat
|
|
|
83 nbloop ldx #$5000
|
|
|
84 ldb #20
|
|
|
85 jsr getline
|
|
|
86 clr b,x
|
|
|
87 cmpb #1
|
|
|
88 lbne donum ; All commands are single-character, everything
|
|
|
89 ; else is treated as a number. Also the single-character lines that
|
|
|
90 ; are not commands are later parsed as numbers.
|
|
|
91 ldb ,x
|
|
|
92 cmpb #'+' ; Add
|
|
|
93 bne nb1
|
|
|
94 jsr fpadd
|
|
|
95 lbra doprint
|
|
|
96 nb1 cmpb #'-' ; Subtract
|
|
|
97 bne nb2
|
|
|
98 jsr fpsub
|
|
|
99 lbra doprint
|
|
|
100 nb2 cmpb #'*' ; Multiply
|
|
|
101 bne nb3
|
|
|
102 jsr fpmul
|
|
|
103 lbra doprint
|
|
|
104 nb3 cmpb #'/' ; Divide
|
|
|
105 bne nb4
|
|
|
106 jsr fpdiv
|
|
|
107 lbra doprint
|
|
|
108 nb4 cmpb #'q' ; Square root.
|
|
|
109 bne nb5
|
|
|
110 jsr fpsqrt
|
|
|
111 lbra doprint
|
|
|
112 nb5 cmpb #'i' ; Round to -Inf INT() in BASIC.
|
|
|
113 bne nb6
|
|
|
114 jsr fpfloor
|
|
|
115 lbra doprint
|
|
|
116 nb6 cmpb #'s' ; SIN() function
|
|
|
117 bne nb7
|
|
|
118 jsr fpsin
|
|
|
119 lbra doprint
|
|
|
120 nb7 cmpb #'=' ; Compare top two numbers Show < = or >
|
|
|
121 bne nb8
|
|
|
122 jsr fpcmp
|
|
|
123 beq nbeq
|
|
|
124 bcc nbgt
|
|
|
125 ldb #'<'
|
|
|
126 bra nbcmp
|
|
|
127 nbeq ldb #'='
|
|
|
128 bra nbcmp
|
|
|
129 nbgt ldb #'>'
|
|
|
130 nbcmp leau -10,u
|
|
|
131 jsr putchar
|
|
|
132 jsr putcr
|
|
|
133 bra nbloop
|
|
|
134 nb8 cmpb #'c' ; COS() function.
|
|
|
135 bne nb9
|
|
|
136 jsr fpcos
|
|
|
137 bra doprint
|
|
|
138 nb9 cmpb #'t' ; TAN() function.
|
|
|
139 bne nb10
|
|
|
140 jsr fptan
|
|
|
141 bra doprint
|
|
|
142 nb10 cmpb #'a' ; ATAN() function.
|
|
|
143 bne nb11
|
|
|
144 jsr fpatan
|
|
|
145 bra doprint
|
|
|
146 nb11 cmpb #'e' ; EXP() function.
|
|
|
147 bne nb12
|
|
|
148 jsr fpexp
|
|
|
149 bra doprint
|
|
|
150 nb12 cmpb #'l' ; LN() function.
|
|
|
151 bne nb13
|
|
|
152 jsr fln
|
|
|
153 bra doprint
|
|
|
154 nb13 cmpb #'d' ; Duplicate top number on stack.
|
|
|
155 bne nb14
|
|
|
156 jsr fpdup
|
|
|
157 bra doprint ; Exchange top two numbers on stack.
|
|
|
158 nb14 cmpb #'x'
|
|
|
159 bne nb15
|
|
|
160 jsr fpexg
|
|
|
161 bra doprint
|
|
|
162 nb15 cmpb #'r' ; Drop top from stack.
|
|
|
163 bne nb16
|
|
|
164 leau -5,u
|
|
|
165 bra doprint
|
|
|
166 nb16
|
|
|
167 donum ldy #$5000
|
|
|
168 jsr scannum
|
|
|
169 lbra nbloop
|
|
|
170 doprint ldy #$5000
|
|
|
171 jsr fpdup
|
|
|
172 jsr fpscient
|
|
|
173 ldx #$5000
|
|
|
174 ldb ,x+
|
|
|
175 jsr putline
|
|
|
176 jsr putcr
|
|
|
177 lbra nbloop
|
|
|
178 nbmesg fcc "Welcome to RPN calculator"
|
|
|
179 nbmend
|
|
|
180
|
|
|
181 doclear rts
|
|
|
182 dorun swi
|
|
|
183 makefree rts
|
|
|
184
|
|
|
185 * Floating point primitives.
|
|
|
186
|
|
|
187 * U is the floating point stack pointer and points to the first free
|
|
|
188 * location. Each number occupies 5 bytes,
|
|
|
189 * Format: byte0: binary exponent $00-$FF $80 means number in range 1..2.
|
|
|
190 * byte1-byte4 binary fraction between 1.0 and 2.0, msb would
|
|
|
191 * always be set, but replaced by sign.
|
|
|
192 * Special case: all bytes zero, number=0.
|
|
|
193
|
|
|
194 * Exchange top two numbers on stack.
|
|
|
195 fpexg ldx -2,u
|
|
|
196 ldd -7,u
|
|
|
197 stx -7,u
|
|
|
198 std -2,u
|
|
|
199 ldx -4,u
|
|
|
200 ldd -9,u
|
|
|
201 stx -9,u
|
|
|
202 std -4,u
|
|
|
203 lda -5,u
|
|
|
204 ldb -10,u
|
|
|
205 sta -10,u
|
|
|
206 stb -5,u
|
|
|
207 rts
|
|
|
208
|
|
|
209 fpdup leax -5,u
|
|
|
210 * Load fp number from address X and push onto stack.
|
|
|
211 fplod ldd ,x++
|
|
|
212 std ,u++
|
|
|
213 ldd ,x++
|
|
|
214 std ,u++
|
|
|
215 lda ,x+
|
|
|
216 sta ,u+
|
|
|
217 fckmem tfr s,d
|
|
|
218 stu fpsp
|
|
|
219 subd fpsp
|
|
|
220 subd #40
|
|
|
221 lbcs makefree ;Test for sufficient free space.
|
|
|
222 rts
|
|
|
223
|
|
|
224 * Pop fp number from stack and store into address X.
|
|
|
225 fpsto lda -5,u
|
|
|
226 sta ,x+
|
|
|
227 ldd -4,u
|
|
|
228 std ,x++
|
|
|
229 ldd -2,u
|
|
|
230 std ,x++
|
|
|
231 leau -5,u
|
|
|
232 rts
|
|
|
233
|
|
|
234 * Compare magnitude (second-top).
|
|
|
235 fpcmpmag lda -10,u
|
|
|
236 cmpa -5,u ;Compare exponents.
|
|
|
237 bne cmpend
|
|
|
238 ldd -4,u
|
|
|
239 anda #$7F ;Eliminate sign bit.
|
|
|
240 std ,--s
|
|
|
241 ldd -9,u
|
|
|
242 anda #$7F ;Eliminate sign bit.
|
|
|
243 subd ,s++ ;Compare msb of mantissa.
|
|
|
244 bne cmpend
|
|
|
245 ldd -7,u
|
|
|
246 subd -2,u
|
|
|
247 bne cmpend
|
|
|
248 cmpend rts
|
|
|
249
|
|
|
250 * Test a top number for 0.
|
|
|
251 fptest0 tst -5,u
|
|
|
252 bne cmpend
|
|
|
253 ldd -4,u
|
|
|
254 bne cmpend
|
|
|
255 ldd -2,u
|
|
|
256 rts
|
|
|
257
|
|
|
258 * Floating point subtraction.
|
|
|
259 fpsub jsr fpneg
|
|
|
260
|
|
|
261 * Floating point addition.
|
|
|
262 fpadd bsr fpcmpmag ;First compare magnitudes.
|
|
|
263 bcc fpadd1
|
|
|
264 jsr fpexg ;Put the biggest one second.
|
|
|
265 fpadd1 bsr fptest0
|
|
|
266 beq fpaddend ;Done if smallest number is 0.
|
|
|
267 lda -10,u
|
|
|
268 suba -5,u ;Determine exponent difference.
|
|
|
269 cmpa #32
|
|
|
270 bhi fpaddend ;Done if difference too big.
|
|
|
271 ldb -9,u
|
|
|
272 andb #$80
|
|
|
273 stb ,-s ;Store sign of biggest number.
|
|
|
274 eorb -4,u
|
|
|
275 stb ,-s ;Store difference of signs.
|
|
|
276 ldb -9,u
|
|
|
277 orb #$80
|
|
|
278 stb -9,u
|
|
|
279 ldb -4,u
|
|
|
280 orb #$80
|
|
|
281 stb -4,u ;Put the hidden msbs back in.
|
|
|
282 clr ,u ;Make extra mantissa byte.
|
|
|
283 tsta
|
|
|
284 beq fpadd2b ;Skip the alignment phase.
|
|
|
285 fpalign lsr -4,u
|
|
|
286 ror -3,u
|
|
|
287 ror -2,u
|
|
|
288 ror -1,u ;Shift the smaller number right to align
|
|
|
289 ror ,u
|
|
|
290 deca
|
|
|
291 bne fpalign
|
|
|
292 fpadd2b tst ,s+
|
|
|
293 bmi dosub ;Did signs differ? Then subtract.
|
|
|
294 ldd -7,u ;Add the mantissas.
|
|
|
295 addd -2,u
|
|
|
296 std -7,u
|
|
|
297 ldd -9,u
|
|
|
298 adcb -3,u
|
|
|
299 adca -4,u
|
|
|
300 std -9,u
|
|
|
301 bcc fpadd2
|
|
|
302 fpadd2a inc -10,u ;Sum overflowed, inc exp, shift mant.
|
|
|
303 lbeq fpovf ;If exponent overflowed, too bad.
|
|
|
304 ror -9,u
|
|
|
305 ror -8,u
|
|
|
306 ror -7,u
|
|
|
307 ror -6,u
|
|
|
308 ror ,u
|
|
|
309 fpadd2 tst ,u
|
|
|
310 bpl fpadd3 ;test msb of extra mantissa byte.
|
|
|
311 ldd -7,u ;Add 1 to mantissa if this is set
|
|
|
312 addd #1
|
|
|
313 std -7,u
|
|
|
314 bcc fpadd3
|
|
|
315 ldd -9,u
|
|
|
316 clr ,u
|
|
|
317 addd #1
|
|
|
318 std -9,u
|
|
|
319 bcs fpadd2a
|
|
|
320 fpadd3 ldb -9,u
|
|
|
321 andb #$7F
|
|
|
322 eorb ,s+
|
|
|
323 stb -9,u ;Put original sign back in.
|
|
|
324 fpaddend leau -5,u
|
|
|
325 rts
|
|
|
326 dosub ldb ,u
|
|
|
327 negb
|
|
|
328 stb ,u
|
|
|
329 ldd -7,u ;Signs differed, so sbutract.
|
|
|
330 sbcb -1,u
|
|
|
331 sbca -2,u
|
|
|
332 std -7,u
|
|
|
333 ldd -9,u
|
|
|
334 sbcb -3,u
|
|
|
335 sbca -4,u
|
|
|
336 std -9,u
|
|
|
337 bmi fpadd2 ;Number still normalized, then done.
|
|
|
338 ldd -9,u
|
|
|
339 bne fpnorm
|
|
|
340 ldd -7,u
|
|
|
341 bne fpnorm
|
|
|
342 tst ,u
|
|
|
343 beq fpundf ;If mantissa exactly zero, underflow.
|
|
|
344 fpnorm tst -10,u ;dec exp, shift mant left
|
|
|
345 beq fpundf ;Underflow, put a zero in.
|
|
|
346 dec -10,u
|
|
|
347 asl ,u
|
|
|
348 rol -6,u
|
|
|
349 rol -7,u
|
|
|
350 rol -8,u
|
|
|
351 rol -9,u
|
|
|
352 bpl fpnorm ;Until number is normalized.
|
|
|
353 bra fpadd2
|
|
|
354
|
|
|
355 fpundf clr -10,u ;Underflow, substitute zero.
|
|
|
356 clr -9,u
|
|
|
357 clr -8,u
|
|
|
358 clr -7,u
|
|
|
359 clr -6,u
|
|
|
360 leas 1,s ;Discard the sign on stack.
|
|
|
361 bra fpaddend
|
|
|
362
|
|
|
363 * Compare Floating Point Numbers, flags as with unsigned comparison.
|
|
|
364 fpcmp lda -9,u
|
|
|
365 anda #$80
|
|
|
366 sta ,-s
|
|
|
367 lda -4,u
|
|
|
368 anda #$80
|
|
|
369 suba ,s+ ;Subtract the signs, subtraction is reversed.
|
|
|
370 bne fpcmpend
|
|
|
371 tst -9,u
|
|
|
372 bmi fpcmpneg ;Are numbers negative?
|
|
|
373 jmp fpcmpmag
|
|
|
374 fpcmpneg jsr fpcmpmag
|
|
|
375 beq fpcmpend
|
|
|
376 tfr cc,a
|
|
|
377 eora #$1
|
|
|
378 tfr a,cc ;Reverse the carry flag.
|
|
|
379 fpcmpend rts
|
|
|
380
|
|
|
381 * Multiply floating point numbers.
|
|
|
382 fpmul lda -9,u
|
|
|
383 eora -4,u
|
|
|
384 anda #$80
|
|
|
385 sta ,-s ;Sign difference to stack.
|
|
|
386 jsr fptest0 ;Test one operand for 0
|
|
|
387 beq fpundf
|
|
|
388 ldd -7,u
|
|
|
389 bne fpmula
|
|
|
390 ldd -9,u
|
|
|
391 bne fpmula ;And the other one.
|
|
|
392 ldb -10,u
|
|
|
393 beq fpundf
|
|
|
394 fpmula ldb -9,u
|
|
|
395 orb #$80
|
|
|
396 stb -9,u
|
|
|
397 ldb -4,u
|
|
|
398 orb #$80
|
|
|
399 stb -4,u ;Put hidden msb back in.
|
|
|
400 lda -10,u
|
|
|
401 suba #$80 ;Make unbiased signed num of exponents.
|
|
|
402 sta ,-s
|
|
|
403 lda -5,u
|
|
|
404 suba #$80
|
|
|
405 adda ,s+ ;add exponents.
|
|
|
406 bvc fpmul1 ;Check over/underflow
|
|
|
407 lbmi fpovf
|
|
|
408 bra fpundf
|
|
|
409 fpmul1 adda #$80 ;Make exponent biased again.
|
|
|
410 sta -10,u ;Store result exponent.
|
|
|
411 * Now perform multiplication of mantissas to 40-bit product.
|
|
|
412 * 0,u--4,u product. 5,u--9,u added term
|
|
|
413 * Having a mul instruction is nice, but using it for an efficient
|
|
|
414 * multiprecision multiplicaton is hard. This routine has 13 mul instructions.
|
|
|
415 lda -1,u
|
|
|
416 ldb -8,u
|
|
|
417 mul ;b4*a2
|
|
|
418 sta 4,u
|
|
|
419 lda -1,u
|
|
|
420 ldb -9,u
|
|
|
421 mul ;b4*a1
|
|
|
422 addb 4,u
|
|
|
423 adca #0
|
|
|
424 std 3,u
|
|
|
425 lda -2,u
|
|
|
426 ldb -7,u
|
|
|
427 mul ;b3*a3
|
|
|
428 sta 9,u
|
|
|
429 lda -2,u
|
|
|
430 ldb -8,u
|
|
|
431 mul ;b3*a2
|
|
|
432 addb 9,u
|
|
|
433 adca #0
|
|
|
434 std 8,u
|
|
|
435 lda -2,u
|
|
|
436 ldb -9,u
|
|
|
437 mul ;b3*a1
|
|
|
438 addb 8,u
|
|
|
439 adca #0
|
|
|
440 std 7,u
|
|
|
441 ldd 8,u
|
|
|
442 addd 3,u
|
|
|
443 std 3,u
|
|
|
444 ldb 7,u
|
|
|
445 adcb #0
|
|
|
446 stb 2,u ;Add b4*a and b3*a partial products.
|
|
|
447 lda -3,u
|
|
|
448 ldb -6,u
|
|
|
449 mul ;b2*a4
|
|
|
450 sta 9,u
|
|
|
451 lda -3,u
|
|
|
452 ldb -7,u
|
|
|
453 mul ;b2*a3
|
|
|
454 addb 9,u
|
|
|
455 adca #0
|
|
|
456 std 8,u
|
|
|
457 lda -3,u
|
|
|
458 ldb -8,u
|
|
|
459 mul ;b2*a2
|
|
|
460 addb 8,u
|
|
|
461 adca #0
|
|
|
462 std 7,u
|
|
|
463 lda -3,u
|
|
|
464 ldb -9,u ;b2*a1
|
|
|
465 mul
|
|
|
466 addb 7,u
|
|
|
467 adca #0
|
|
|
468 std 6,u
|
|
|
469 ldd 8,u
|
|
|
470 addd 3,u
|
|
|
471 std 3,u
|
|
|
472 ldd 6,u
|
|
|
473 adcb 2,u
|
|
|
474 adca #0
|
|
|
475 std 1,u ;Add b2*a partial product in.
|
|
|
476 lda -4,u
|
|
|
477 ldb -6,u
|
|
|
478 mul ;b1*a4
|
|
|
479 std 8,u
|
|
|
480 lda -4,u
|
|
|
481 ldb -7,u
|
|
|
482 mul ;b1*a3
|
|
|
483 addb 8,u
|
|
|
484 adca #0
|
|
|
485 std 7,u
|
|
|
486 lda -4,u
|
|
|
487 ldb -8,u
|
|
|
488 mul ;b1*a2
|
|
|
489 addb 7,u
|
|
|
490 adca #0
|
|
|
491 std 6,u
|
|
|
492 lda -4,u
|
|
|
493 ldb -9,u
|
|
|
494 mul ;b1*a1
|
|
|
495 addb 6,u
|
|
|
496 adca #0
|
|
|
497 std 5,u
|
|
|
498 ldd 8,u
|
|
|
499 addd 3,u
|
|
|
500 std -6,u
|
|
|
501 ldd 6,u
|
|
|
502 adcb 2,u
|
|
|
503 adca 1,u
|
|
|
504 std -8,u
|
|
|
505 ldb 5,u
|
|
|
506 adcb #0
|
|
|
507 stb -9,u ;Add product term b1*a in, result to dest.
|
|
|
508 bmi fpmul2
|
|
|
509 asl -5,u
|
|
|
510 rol -6,u
|
|
|
511 rol -7,u
|
|
|
512 rol -8,u
|
|
|
513 rol -9,u ;Normalize by shifting mantissa left.
|
|
|
514 bra fpmul3
|
|
|
515 fpmul2 inc -10,u ;increment exponent.
|
|
|
516 lbeq fpovf ;Test for overflow.
|
|
|
517 fpmul3 tst -5,u
|
|
|
518 lbpl fpadd3
|
|
|
519 ldd -7,u ;Add 1 if msb of 5th nibble is set.
|
|
|
520 addd #1
|
|
|
521 std -7,u
|
|
|
522 lbcc fpadd3
|
|
|
523 ldd -9,u
|
|
|
524 addd #1
|
|
|
525 std -9,u
|
|
|
526 bcs fpmul4 ;It could overflow.
|
|
|
527 lbra fpadd3
|
|
|
528 fpmul4 clr -5,u
|
|
|
529 bra fpmul2
|
|
|
530
|
|
|
531 * Divide floating point numbers.
|
|
|
532 fpdiv lda -9,u
|
|
|
533 eora -4,u
|
|
|
534 anda #$80
|
|
|
535 sta ,-s ;Sign difference to stack.
|
|
|
536 jsr fptest0 ;Test divisor for 0
|
|
|
537 lbeq fpovf
|
|
|
538 ldd -7,u
|
|
|
539 bne fpdiva
|
|
|
540 ldd -9,u
|
|
|
541 bne fpdiva ;And the other one.
|
|
|
542 ldb -10,u
|
|
|
543 lbeq fpundf
|
|
|
544 fpdiva ldb -9,u
|
|
|
545 orb #$80
|
|
|
546 stb -9,u
|
|
|
547 ldb -4,u
|
|
|
548 orb #$80
|
|
|
549 stb -4,u ;Put hidden msb back in.
|
|
|
550 lda -5,u
|
|
|
551 suba #$80 ;Make unbiased signed difference of exponents.
|
|
|
552 sta ,-s
|
|
|
553 lda -10,u
|
|
|
554 suba #$80
|
|
|
555 suba ,s+ ;subtract exponents.
|
|
|
556 bvc fpdiv1 ;Check over/underflow
|
|
|
557 lbmi fpovf
|
|
|
558 lbra fpundf
|
|
|
559 fpdiv1 adda #$80 ;Make exponent biased again.
|
|
|
560 sta -10,u ;Store result exponent.
|
|
|
561 * Now start the division of mantissas. Temprorary 34-bit quotient in 0,u--4,u
|
|
|
562 * -5,u is extra byte of dividend.
|
|
|
563 lda #34
|
|
|
564 sta ,-s
|
|
|
565 clr ,u
|
|
|
566 clr 1,u
|
|
|
567 clr 2,u
|
|
|
568 clr 3,u
|
|
|
569 clr 4,u
|
|
|
570 clr -5,u
|
|
|
571 fpdivloop asl 4,u ;Shift quotient left.
|
|
|
572 rol 3,u
|
|
|
573 rol 2,u
|
|
|
574 rol 1,u
|
|
|
575 rol ,u
|
|
|
576 ldd -7,u ;Perform trial subtraction.
|
|
|
577 subd -2,u
|
|
|
578 std -7,u
|
|
|
579 ldd -9,u
|
|
|
580 sbcb -3,u
|
|
|
581 sbca -4,u
|
|
|
582 std -9,u
|
|
|
583 ldb -5,u
|
|
|
584 sbcb #0
|
|
|
585 bcc fpdiv2
|
|
|
586 ldd -7,u ;Undo the trial subtraction.
|
|
|
587 addd -2,u
|
|
|
588 std -7,u
|
|
|
589 ldd -9,u
|
|
|
590 adcb -3,u
|
|
|
591 adca -4,u
|
|
|
592 std -9,u
|
|
|
593 bra fpdiv4
|
|
|
594 fpdiv2 stb -5,u ;Store new msb of quotient.
|
|
|
595 lda 4,u ;Add 1 to quotient.
|
|
|
596 adda #$40
|
|
|
597 sta 4,u
|
|
|
598 fpdiv4 asl -6,u ;Shift dividend left.
|
|
|
599 rol -7,u
|
|
|
600 rol -8,u
|
|
|
601 rol -9,u
|
|
|
602 rol -5,u
|
|
|
603 dec ,s
|
|
|
604 bne fpdivloop
|
|
|
605 leas 1,s
|
|
|
606 ldd 3,u
|
|
|
607 std -6,u
|
|
|
608 ldd 1,u
|
|
|
609 std -8,u
|
|
|
610 ldb ,u
|
|
|
611 stb -9,u ;Move quotient to final location.
|
|
|
612 bmi fpdiv3
|
|
|
613 fpdiv5 asl -5,u
|
|
|
614 rol -6,u
|
|
|
615 rol -7,u
|
|
|
616 rol -8,u
|
|
|
617 rol -9,u ;Normalize by shifting mantissa left.
|
|
|
618 ldb -10,u ;decrement exponent.
|
|
|
619 lbeq fpundf ;Test for underflow.
|
|
|
620 decb
|
|
|
621 stb -10,u
|
|
|
622 fpdiv3 tst -5,u
|
|
|
623 lbpl fpadd3
|
|
|
624 ldd -7,u ;Add 1 if msb of 5th nibble is set.
|
|
|
625 addd #1
|
|
|
626 std -7,u
|
|
|
627 lbcc fpadd3
|
|
|
628 ldd -9,u
|
|
|
629 addd #1
|
|
|
630 std -9,u
|
|
|
631 lbcs fpmul4 ;This addition could overflow.
|
|
|
632 lbra fpadd3
|
|
|
633
|
|
|
634 * Floating point negation.
|
|
|
635 fpneg jsr fptest0
|
|
|
636 beq fpnegend ;Do nothing if number equals zero.
|
|
|
637 lda -4,u
|
|
|
638 eora #$80
|
|
|
639 sta -4,u ;Invert the sign bit.
|
|
|
640 fpnegend rts
|
|
|
641
|
|
|
642 * Convert unsigned double number at X to float.
|
|
|
643 ufloat leau 5,u ;Make room for extra number on stack.
|
|
|
644 ldd ,x
|
|
|
645 std -4,u
|
|
|
646 ldd 2,x
|
|
|
647 clr -5,u
|
|
|
648 uf16 std -2,u ;Transfer integer to FP number.
|
|
|
649 jsr fptest0
|
|
|
650 beq ufzero
|
|
|
651 ldb #$9f ;Number is not zero.
|
|
|
652 stb -5,u
|
|
|
653 tst -4,u
|
|
|
654 bmi ufdone
|
|
|
655 ufloop dec -5,u ;Decrement exponent.
|
|
|
656 asl -1,u
|
|
|
657 rol -2,u
|
|
|
658 rol -3,u
|
|
|
659 rol -4,u ;Shift mantissa.
|
|
|
660 bpl ufloop ;until normalized.
|
|
|
661 ufdone ldb -4,u
|
|
|
662 andb #$7f
|
|
|
663 stb -4,u ;Remove the hidden msb.
|
|
|
664 ufend jmp fckmem ;Check that fp stack does not overflow
|
|
|
665 ufzero clr -5,u ;Make exponent zero as well.
|
|
|
666 bra ufend
|
|
|
667
|
|
|
668 * Convert unsigned 16-bit integer in D to floating point.
|
|
|
669 unint2fp clr ,-s
|
|
|
670 bra i2fp2
|
|
|
671 * Convert signed 16-bit integer in D to floating point.
|
|
|
672 int2fp sta ,-s ;Store sign byte.
|
|
|
673 bpl i2fp2
|
|
|
674 comb
|
|
|
675 coma
|
|
|
676 addd #1 ;Negate D if negative.
|
|
|
677 i2fp2 leau 5,u
|
|
|
678 clr -4,u
|
|
|
679 clr -5,u
|
|
|
680 clr -3,u ;Clear msb
|
|
|
681 jsr uf16
|
|
|
682 tst ,s+
|
|
|
683 bmi fpneg
|
|
|
684 rts ;Negate number if it was negative.
|
|
|
685
|
|
|
686 * Convert float to unsigned 32-bit integer at X.
|
|
|
687 * A is nonzero if number was not integer or zero.
|
|
|
688 uint ldd -4,u
|
|
|
689 ora #$80 ;Put the hidden msb back in.
|
|
|
690 std ,x
|
|
|
691 ldd -2,u
|
|
|
692 std 2,x ;Transfer mantissa.
|
|
|
693 clra
|
|
|
694 ldb -5,u
|
|
|
695 cmpb #$80 ;If less than 1, it's 0
|
|
|
696 blo uizero
|
|
|
697 cmpb #$9f
|
|
|
698 lbhi intrange ;2^32 or higher, that's too bad.
|
|
|
699 beq uidone
|
|
|
700 uiloop lsr ,x
|
|
|
701 ror 1,x
|
|
|
702 ror 2,x
|
|
|
703 ror 3,x ;Adjust integer by shifting to right
|
|
|
704 adca #0 ;Add any shifted out bit into A.
|
|
|
705 incb
|
|
|
706 cmpb #$9f
|
|
|
707 blo uiloop
|
|
|
708 uidone leau -5,u
|
|
|
709 rts
|
|
|
710 uizero inca ; Indicate non-integer.
|
|
|
711 clr ,x ; Number is zero
|
|
|
712 clr 1,x
|
|
|
713 clr 2,x
|
|
|
714 clr 3,x
|
|
|
715 leau -5,u
|
|
|
716 rts
|
|
|
717
|
|
|
718 * Convert fp number to signed or unsigned 16-bit number in D.
|
|
|
719 * Acceptable values are -65535..65535.
|
|
|
720 fp2uint ldb -5,u
|
|
|
721 stb ,-s ;Store sign.
|
|
|
722 ldx #intbuf
|
|
|
723 bsr uint
|
|
|
724 ldx ,x
|
|
|
725 lbne intrange ;Integer must be in 16-bit range.
|
|
|
726 ldd intbuf+2
|
|
|
727 tst ,s+
|
|
|
728 bpl fp2iend
|
|
|
729 comb
|
|
|
730 coma
|
|
|
731 addd #1 ;Negate number if negative.
|
|
|
732 fp2iend rts
|
|
|
733 * Convert fp number to signed 16-bit number in D.
|
|
|
734 fp2int ldb -5,u
|
|
|
735 stb ,-s ;Store sign of FP number.
|
|
|
736 bsr fp2uint
|
|
|
737 pshs d
|
|
|
738 eora ,s+
|
|
|
739 lbmi intrange ;Compare sign to what it should be.
|
|
|
740 puls d,pc
|
|
|
741
|
|
|
742 * Scan a number at address Y and convert to integer or floating point
|
|
|
743 scannum jsr skipspace
|
|
|
744 clr ,-s ;Store sign on stack.
|
|
|
745 cmpb #'-' ;Test for minus sign.
|
|
|
746 bne sn1
|
|
|
747 inc ,s ;Set sign on stack
|
|
|
748 ldb ,y+
|
|
|
749 sn1 jsr scanint ;First scan the number as an integer.
|
|
|
750 ldx #intbuf
|
|
|
751 jsr ufloat ;Convert to float.
|
|
|
752 ldb -1,y
|
|
|
753 sn1loop cmpb #'.'
|
|
|
754 bne sn1c
|
|
|
755 tst dpl ;If dpl already set, accept no other point.
|
|
|
756 bne sn1d
|
|
|
757 inc dpl
|
|
|
758 ldb ,y+
|
|
|
759 bra sn1loop
|
|
|
760 sn1c subb #'0'
|
|
|
761 blo sn1d
|
|
|
762 cmpb #9
|
|
|
763 bhi sn1d
|
|
|
764 clra
|
|
|
765 jsr int2fp ;Convert digit to fp
|
|
|
766 jsr fpexg
|
|
|
767 ldx #fpten
|
|
|
768 jsr fplod
|
|
|
769 jsr fpmul ;Multiply original number by 10.
|
|
|
770 jsr fpadd ;Add digit to it.
|
|
|
771 tst dpl
|
|
|
772 beq sn1k
|
|
|
773 inc dpl ;Adjust dpl (one more digit after .)
|
|
|
774 sn1k ldb ,y+
|
|
|
775 bra sn1loop
|
|
|
776 sn1d tst ,s+
|
|
|
777 beq sn1a
|
|
|
778 jsr fpneg ;Negate the number if negative.
|
|
|
779 sn1a clr ,-s
|
|
|
780 clr ,-s ;Prepare exponent part on stack.
|
|
|
781 ldb -1,y
|
|
|
782 cmpb #'e'
|
|
|
783 beq sn1e
|
|
|
784 cmpb #'E'
|
|
|
785 bne sn1f ;Test for exponent part.
|
|
|
786 sn1e ldb ,y+
|
|
|
787 clr ,-s ;Prepare exponent sign on stack.
|
|
|
788 cmpb #'+'
|
|
|
789 beq sn1g
|
|
|
790 cmpb #'-'
|
|
|
791 bne sn1h
|
|
|
792 inc ,s ;Set sign to negative.
|
|
|
793 sn1g ldb ,y+
|
|
|
794 sn1h lda dpl
|
|
|
795 pshs a
|
|
|
796 clr dpl
|
|
|
797 inc dpl
|
|
|
798 jsr scanint ;Scan the exponent part.
|
|
|
799 puls a
|
|
|
800 sta dpl ;Restore dpl.
|
|
|
801 lda intbuf
|
|
|
802 ora intbuf+1
|
|
|
803 ora intbuf+2
|
|
|
804 lbne fpovf ;Exponent may not be greater than 255.
|
|
|
805 ldb intbuf+3
|
|
|
806 lbmi fpovf ;Not even greater than 127.
|
|
|
807 tst ,s+
|
|
|
808 beq sn1i
|
|
|
809 negb
|
|
|
810 sn1i sex
|
|
|
811 std ,s
|
|
|
812 sn1f ldb dpl
|
|
|
813 beq sn1j
|
|
|
814 decb
|
|
|
815 sn1j negb
|
|
|
816 sex
|
|
|
817 addd ,s++ ;Add exponent part as well
|
|
|
818 pshs d
|
|
|
819 ldx #fpten
|
|
|
820 jsr fplod
|
|
|
821 puls d
|
|
|
822 jsr fpipower
|
|
|
823 jsr fpmul
|
|
|
824 sn1b rts
|
|
|
825
|
|
|
826 * Scan integer number below 1e9 at address Y, first digit in B.
|
|
|
827 scanint clr dpl
|
|
|
828 scanint1 clr intbuf
|
|
|
829 clr intbuf+1
|
|
|
830 clr intbuf+2
|
|
|
831 clr intbuf+3 ;Initialize number
|
|
|
832 snloop cmpb #'.'
|
|
|
833 bne sn2a ;Test for decimal point.
|
|
|
834 tst dpl
|
|
|
835 bne sndone ;Done if second point found.
|
|
|
836 inc dpl ;Set dpl to indicate decimal point.
|
|
|
837 bra sn3
|
|
|
838 sn2a subb #'0'
|
|
|
839 blo sndone
|
|
|
840 cmpb #9
|
|
|
841 bhi sndone ;Check that character is a digit.
|
|
|
842 tst dpl
|
|
|
843 beq sn2b
|
|
|
844 inc dpl ;Incremend deecimal point loc if set.
|
|
|
845 sn2b pshs b
|
|
|
846 ldd intbuf+2
|
|
|
847 aslb
|
|
|
848 rola
|
|
|
849 std intbuf+2
|
|
|
850 std intbuf2+2
|
|
|
851 ldd intbuf
|
|
|
852 rolb
|
|
|
853 rola
|
|
|
854 std intbuf
|
|
|
855 std intbuf2
|
|
|
856 asl intbuf+3
|
|
|
857 rol intbuf+2
|
|
|
858 rol intbuf+1
|
|
|
859 rol intbuf
|
|
|
860 asl intbuf+3
|
|
|
861 rol intbuf+2
|
|
|
862 rol intbuf+1
|
|
|
863 rol intbuf
|
|
|
864 ldd intbuf+2
|
|
|
865 addd intbuf2+2
|
|
|
866 std intbuf +2
|
|
|
867 ldd intbuf
|
|
|
868 adcb intbuf2+1
|
|
|
869 adca intbuf2
|
|
|
870 std intbuf ;Multiply the integer by 10
|
|
|
871 ldd intbuf+2
|
|
|
872 addb ,s+ ;Add the digit in.
|
|
|
873 adca #0
|
|
|
874 std intbuf+2
|
|
|
875 bcc sn2
|
|
|
876 ldd intbuf
|
|
|
877 addd #1
|
|
|
878 std intbuf
|
|
|
879 sn2 ldd intbuf
|
|
|
880 cmpd #$5f5
|
|
|
881 blo sn3
|
|
|
882 bhi snovf
|
|
|
883 ldd intbuf+2 ;note $5f5e100 is 100 million
|
|
|
884 cmpd #$e100 ;Compare result to 100 million
|
|
|
885 bhs snovf
|
|
|
886 sn3 ldb ,y+ ;get next digit.
|
|
|
887 bra snloop
|
|
|
888 snovf ldb ,y+ ;get next digit.
|
|
|
889 sndone ldb -1,y
|
|
|
890 rts
|
|
|
891
|
|
|
892 *Convert integer at X to BCD.
|
|
|
893 int2bcd clr bcdbuf
|
|
|
894 clr bcdbuf+1
|
|
|
895 clr bcdbuf+2
|
|
|
896 clr bcdbuf+3
|
|
|
897 clr bcdbuf+4
|
|
|
898 ldb #4
|
|
|
899 tstzero tst ,x+
|
|
|
900 bne bcd1
|
|
|
901 decb
|
|
|
902 bne tstzero ;Skip bytes that are zero.
|
|
|
903 bra sndone ;Done if number already zero.
|
|
|
904 bcd1 stb ,-s ;Store number of bytes.
|
|
|
905 leax -1,x
|
|
|
906 bcdloop ldb #8
|
|
|
907 bcdloop1 rol ,x ;Get next bit of binary nunber
|
|
|
908 lda bcdbuf+4
|
|
|
909 adca bcdbuf+4
|
|
|
910 daa
|
|
|
911 sta bcdbuf+4
|
|
|
912 lda bcdbuf+3
|
|
|
913 adca bcdbuf+3
|
|
|
914 daa
|
|
|
915 sta bcdbuf+3
|
|
|
916 lda bcdbuf+2
|
|
|
917 adca bcdbuf+2
|
|
|
918 daa
|
|
|
919 sta bcdbuf+2
|
|
|
920 lda bcdbuf+1
|
|
|
921 adca bcdbuf+1
|
|
|
922 daa
|
|
|
923 sta bcdbuf+1
|
|
|
924 lda bcdbuf
|
|
|
925 adca bcdbuf
|
|
|
926 daa
|
|
|
927 sta bcdbuf ;Add BCD number to itself plus the extra bit.
|
|
|
928 decb
|
|
|
929 bne bcdloop1
|
|
|
930 leax 1,x
|
|
|
931 dec ,s
|
|
|
932 bne bcdloop
|
|
|
933 leas 1,s ;Remove counter from stack.
|
|
|
934 rts
|
|
|
935
|
|
|
936 * Raise fp number to an integer power contained in D.
|
|
|
937 fpipower sta ,-s ;Store sign of exponent.
|
|
|
938 bpl fppow1 ;Is exponent negative.
|
|
|
939 coma
|
|
|
940 comb
|
|
|
941 addd #1 ;Take absolute value of exponent.
|
|
|
942 fppow1 std ,--s ;Store the exponent.
|
|
|
943 ldx #fpone
|
|
|
944 jsr fplod ;Start with number one.
|
|
|
945 fppowloop lsr ,s
|
|
|
946 ror 1,s ;Divide exponent by 2.
|
|
|
947 bcc fppow2 ;Test if it was odd.
|
|
|
948 leax -10,u
|
|
|
949 jsr fplod
|
|
|
950 jsr fpmul ;Multiply result by factor.
|
|
|
951 fppow2 ldd ,s
|
|
|
952 beq fppowdone ;Is exponent zero?
|
|
|
953 leax -10,u
|
|
|
954 jsr fplod
|
|
|
955 jsr fpdup
|
|
|
956 jsr fpmul ;Sqaure the factor.
|
|
|
957 leax -15,u
|
|
|
958 jsr fpsto ;Store it in its place on stack.
|
|
|
959 bra fppowloop
|
|
|
960 fppowdone leas 2,s ;Remove exponent.
|
|
|
961 tst ,s+
|
|
|
962 bpl fppow3 ;Was exponent negative?
|
|
|
963 ldx #fpone
|
|
|
964 jsr fplod
|
|
|
965 jsr fpexg
|
|
|
966 jsr fpdiv :compute 1/result.
|
|
|
967 fppow3 jsr fpexg
|
|
|
968 leau -5,u ;Remove factor from stack.
|
|
|
969 rts
|
|
|
970
|
|
|
971
|
|
|
972 * Convert fp number to string at address Y in scientific notation.
|
|
|
973 fpscient ldb #15
|
|
|
974 stb ,y+ ;Store the string length.
|
|
|
975 lda #' '
|
|
|
976 ldb -4,u
|
|
|
977 bpl fpsc1
|
|
|
978 lda #'-'
|
|
|
979 fpsc1 sta ,y+ ;Store - or space depending on sign.
|
|
|
980 andb #$7f
|
|
|
981 stb -4,u ;Make number positive.
|
|
|
982 clr ,-s ;Store decimal exponent (default 0)
|
|
|
983 jsr fptest0
|
|
|
984 beq fpsc2 ;Test for zero
|
|
|
985 lda -5,u
|
|
|
986 suba #$80
|
|
|
987 suba #$1D ;Adjust exponent.
|
|
|
988 bvc fpsc11a
|
|
|
989 lda #-128
|
|
|
990 fpsc11a sta ,-s ;store it to recover sign later.
|
|
|
991 bpl posexp
|
|
|
992 nega ;Take absolute value.
|
|
|
993 posexp ldb #5
|
|
|
994 mul
|
|
|
995 lsra
|
|
|
996 rorb
|
|
|
997 lsra
|
|
|
998 rorb
|
|
|
999 lsra
|
|
|
1000 rorb
|
|
|
1001 lsra
|
|
|
1002 rorb ;multiply by 5/16 approx 10log 2
|
|
|
1003 cmpb #37
|
|
|
1004 bls expmax
|
|
|
1005 ldb #37 ;Maximum decimal exponent=37
|
|
|
1006 expmax tst ,s+
|
|
|
1007 bpl posexp1
|
|
|
1008 negb
|
|
|
1009 posexp1 stb ,s ;Store approximate decimal exponent.
|
|
|
1010 negb
|
|
|
1011 sex ;Approximate (negated) decimal exponent in D.
|
|
|
1012 pshs d
|
|
|
1013 ldx #fpten
|
|
|
1014 jsr fplod
|
|
|
1015 puls d
|
|
|
1016 jsr fpipower ;Take 10^-exp
|
|
|
1017 jsr fpmul
|
|
|
1018 fpsc1a ldx #fplolim
|
|
|
1019 jsr fplod
|
|
|
1020 jsr fpcmpmag ;Compare number to 100 million
|
|
|
1021 leau -5,u
|
|
|
1022 bhs fpsc1c
|
|
|
1023 dec ,s ;Decrement approximate exponent.
|
|
|
1024 ldx #fpten
|
|
|
1025 jsr fplod
|
|
|
1026 jsr fpmul ;Multiply by ten.
|
|
|
1027 bra fpsc1a
|
|
|
1028 fpsc1c ldx #fphilim
|
|
|
1029 jsr fplod
|
|
|
1030 jsr fpcmpmag ;Compare number to 1 billion
|
|
|
1031 leau -5,u
|
|
|
1032 blo fpsc1d
|
|
|
1033 inc ,s ;Increment approximate exponent.
|
|
|
1034 ldx #fpten
|
|
|
1035 jsr fplod
|
|
|
1036 jsr fpdiv ;Divide by ten.
|
|
|
1037 bra fpsc1c
|
|
|
1038 fpsc1d ldb ,s
|
|
|
1039 addb #8
|
|
|
1040 stb ,s ;Adjust decimal exponent (8 decimals)
|
|
|
1041 ldx #fphalf
|
|
|
1042 jsr fplod
|
|
|
1043 jsr fpadd ;Add 0.5 for the final round to integer.
|
|
|
1044 * Number is either zero or between 100 million and 1 billion.
|
|
|
1045 fpsc2 ldx #intbuf
|
|
|
1046 jsr uint ;Convert decimal mantissa to integer.
|
|
|
1047 jsr int2bcd ;Convert to bcd.
|
|
|
1048 ldb bcdbuf
|
|
|
1049 addb #'0'
|
|
|
1050 stb ,y+ ;Store digit before decimal point
|
|
|
1051 ldb #'.'
|
|
|
1052 stb ,y+ ;Store decimal point.
|
|
|
1053 lda #4
|
|
|
1054 sta ,-s
|
|
|
1055 ldx #bcdbuf+1
|
|
|
1056 fpscloop lda ,x+
|
|
|
1057 tfr a,b
|
|
|
1058 lsrb
|
|
|
1059 lsrb
|
|
|
1060 lsrb
|
|
|
1061 lsrb
|
|
|
1062 addb #'0'
|
|
|
1063 stb ,y+
|
|
|
1064 anda #$0f
|
|
|
1065 adda #'0
|
|
|
1066 sta ,y+
|
|
|
1067 dec ,s ;Convert the other 8 digits to ASCII
|
|
|
1068 bne fpscloop
|
|
|
1069 leas 1,s ;Remove loop counter.
|
|
|
1070 ldb #'E'
|
|
|
1071 stb ,y+ ;Store the E character.
|
|
|
1072 lda #'+'
|
|
|
1073 ldb ,s+ ;Get decimal exponent.
|
|
|
1074 bpl fpsc3 ;Test sign of exponent.
|
|
|
1075 lda #'-'
|
|
|
1076 negb ;Take absolute value of exponent.
|
|
|
1077 fpsc3 sta ,y+ ;Store sign of exponent.
|
|
|
1078 stb intbuf+3
|
|
|
1079 clr intbuf+2
|
|
|
1080 clr intbuf+1
|
|
|
1081 clr intbuf
|
|
|
1082 ldx #intbuf
|
|
|
1083 jsr int2bcd ;Convert decimal exponent to bcd.
|
|
|
1084 lda bcdbuf+4
|
|
|
1085 tfr a,b
|
|
|
1086 lsrb
|
|
|
1087 lsrb
|
|
|
1088 lsrb
|
|
|
1089 lsrb
|
|
|
1090 addb #'0'
|
|
|
1091 stb ,y+ ;Convert first exp digit to ascii
|
|
|
1092 anda #$0f
|
|
|
1093 adda #'0'
|
|
|
1094 sta ,y+ ;And the second one.
|
|
|
1095 rts
|
|
|
1096
|
|
|
1097
|
|
|
1098 include "floatnum.inc"
|
|
|
1099
|
|
|
1100 fpovf swi
|
|
|
1101 intrange swi
|
|
|
1102 inval swi
|
|
|
1103
|
|
|
1104 * This routine takes the square root of an FP number.
|
|
|
1105 * Uses Newton's algorithm.
|
|
|
1106 fpsqrt tst -4,u
|
|
|
1107 lbmi inval ;Negative arguments are invalid.
|
|
|
1108 jsr fptest0
|
|
|
1109 beq sqdone ;Sqaure root of 0 is 0.
|
|
|
1110 jsr fpdup
|
|
|
1111 ldb -5,u
|
|
|
1112 subb #$80 ;Unbias the exponent.
|
|
|
1113 bpl sq1
|
|
|
1114 addb #1
|
|
|
1115 sq1 asrb ;Divide exponent by 2.
|
|
|
1116 addb #$80 ;Make it biased again.
|
|
|
1117 stb -5,u ;This is the initial guess for the root.
|
|
|
1118 ldb #4 ;Do the loop 4 times.
|
|
|
1119 stb ,-s
|
|
|
1120 sqloop leax -10,u
|
|
|
1121 jsr fplod
|
|
|
1122 leax -10,u
|
|
|
1123 jsr fplod
|
|
|
1124 jsr fpdiv ;Divide argument by guess.
|
|
|
1125 jsr fpadd ;Add to guess.
|
|
|
1126 dec -5,u ;Divide this by two, giving new guess.
|
|
|
1127 dec ,s
|
|
|
1128 bne sqloop
|
|
|
1129 leas 1,s
|
|
|
1130 jsr fpexg
|
|
|
1131 leau -5,u ;Remove argument, leave final guess.
|
|
|
1132 sqdone rts
|
|
|
1133
|
|
|
1134 * Compute the floor of an fp number (result is still fp.
|
|
|
1135 fpfloor ldb -5,u
|
|
|
1136 cmpb #$9f
|
|
|
1137 bhs sqdone ;If abs value >=2^31, then already integer.
|
|
|
1138 ldb -4,u
|
|
|
1139 stb ,-s ;Stroe sign of number
|
|
|
1140 andb #$7f
|
|
|
1141 stb -4,u ;Take absolute value of number.
|
|
|
1142 ldx #intbuf
|
|
|
1143 jsr uint ;Convert to int (truncation)
|
|
|
1144 sta ,-s ;Store number of fraction bits.
|
|
|
1145 ldx #intbuf
|
|
|
1146 jsr ufloat ;Convert back to float
|
|
|
1147 ldd ,s++
|
|
|
1148 tstb
|
|
|
1149 bpl sqdone
|
|
|
1150 sta ,-s
|
|
|
1151 jsr fpneg ;Negate number if it was negative
|
|
|
1152 lda ,s+
|
|
|
1153 beq sqdone
|
|
|
1154 ldx #fpone
|
|
|
1155 jsr fplod
|
|
|
1156 jmp fpsub ;Subtract 1 if negative & not integer.
|
|
|
1157
|
|
|
1158 * Floating point modulo operation (floored modulo).
|
|
|
1159 * Integer part of quotient is still left in intbuf
|
|
|
1160 fpmod leax -10,u
|
|
|
1161 jsr fplod
|
|
|
1162 leax -10,u
|
|
|
1163 jsr fplod
|
|
|
1164 jsr fpdiv ;Perform division.
|
|
|
1165 jsr fpfloor
|
|
|
1166 jsr fpmul ;Multiply Quotient and Divisor
|
|
|
1167 leax -10,u
|
|
|
1168 jmp fpsub ;Dividend - quotient*divisor = modulus.
|
|
|
1169
|
|
|
1170
|
|
|
1171 * Now the transcendental functions follow.
|
|
|
1172 * They use approximation polynomials as defined in the
|
|
|
1173 * Handbook of Mathematical Functions by Abramowitz & Stegun.
|
|
|
1174
|
|
|
1175 * Compute polynomial, number of terms in B, coefficients start at Y
|
|
|
1176 fppoly stb ,-s
|
|
|
1177 ldx #fpzero
|
|
|
1178 jsr fplod ;Start with zero.
|
|
|
1179 polyloop leax ,y
|
|
|
1180 jsr fplod
|
|
|
1181 jsr fpadd ;Add next coefficient.
|
|
|
1182 leay 5,y
|
|
|
1183 leax -10,u
|
|
|
1184 jsr fplod
|
|
|
1185 jsr fpmul ;Multiply by x.
|
|
|
1186 dec ,s
|
|
|
1187 bne polyloop
|
|
|
1188 leas 1,s
|
|
|
1189 jsr fpexg
|
|
|
1190 leau -5,u ;Remove x from stack.
|
|
|
1191 rts
|
|
|
1192
|
|
|
1193 add1 ldx #fpone
|
|
|
1194 jsr fplod
|
|
|
1195 jsr fpadd
|
|
|
1196 rts
|
|
|
1197
|
|
|
1198 halfpi ldx #fpi
|
|
|
1199 jsr fplod
|
|
|
1200 dec -5,u
|
|
|
1201 rts
|
|
|
1202
|
|
|
1203 * sin(x)
|
|
|
1204 fpsin ldx #fpi
|
|
|
1205 jsr fplod
|
|
|
1206 inc -5,u ;Load 2*pi
|
|
|
1207 jsr fpmod ;Modulo 2pi
|
|
|
1208 bsr halfpi
|
|
|
1209 jsr fpcmp ;Compare x to pi/2
|
|
|
1210 bls sin2
|
|
|
1211 inc -5,u ;Change pi/2 to pi
|
|
|
1212 jsr fpsub
|
|
|
1213 jsr fpneg ;x := pi-x if x>pi/2
|
|
|
1214 bsr halfpi
|
|
|
1215 jsr fpneg
|
|
|
1216 jsr fpcmp ;Compare x to -pi/2
|
|
|
1217 bhs sin2
|
|
|
1218 inc -5,u ;Change -pi/2 to -pi
|
|
|
1219 jsr fpsub
|
|
|
1220 jsr fpneg
|
|
|
1221 bra sin3
|
|
|
1222 sin2 leau -5,u ;Drop the compare limit pi/2 or -pi/2
|
|
|
1223 sin3 jsr fpdup
|
|
|
1224 jsr fpdup
|
|
|
1225 jsr fpmul ;On stack: x, x*x
|
|
|
1226 ldy #sincoeff
|
|
|
1227 ldb #5
|
|
|
1228 jsr fppoly ;Do the sine polynomial with x*x as argument
|
|
|
1229 jsr add1 ;Add 1 to the result.
|
|
|
1230 jmp fpmul ;multiply the polynomial result with x.
|
|
|
1231 * cos(x)
|
|
|
1232 fpcos jsr halfpi
|
|
|
1233 jsr fpsub
|
|
|
1234 jsr fpneg
|
|
|
1235 bra fpsin ;Compute sin(pi/2-x)
|
|
|
1236
|
|
|
1237 * tan(x)
|
|
|
1238 fptan jsr fpdup
|
|
|
1239 jsr fpsin
|
|
|
1240 jsr fpexg
|
|
|
1241 jsr fpcos
|
|
|
1242 jmp fpdiv ;Compute sin(x)/cos(x)
|
|
|
1243
|
|
|
1244 * atan(x)
|
|
|
1245 fpatan clr ,-s ;Make flag on stack
|
|
|
1246 ldb -5,u
|
|
|
1247 cmpb #$80 ;Compare magnitude to 1.
|
|
|
1248 blo atn1
|
|
|
1249 inc ,s ;Set flag on stack.
|
|
|
1250 ldx #fpone ;if x>1 then compute 1/x
|
|
|
1251 jsr fplod
|
|
|
1252 jsr fpexg
|
|
|
1253 jsr fpdiv
|
|
|
1254 atn1 jsr fpdup
|
|
|
1255 jsr fpdup
|
|
|
1256 jsr fpmul ;On stack: x, x*x
|
|
|
1257 ldb #8
|
|
|
1258 ldy #atancoeff
|
|
|
1259 jsr fppoly ;Doe the arctan polynomyal, x*x as argument.
|
|
|
1260 jsr add1 ;Add 1 to result
|
|
|
1261 jsr fpmul ;multiply result by x.
|
|
|
1262 tst ,s+
|
|
|
1263 beq atndone
|
|
|
1264 jsr halfpi
|
|
|
1265 jsr fpsub
|
|
|
1266 jsr fpneg ;Compute pi/2 - result when x was >1
|
|
|
1267 atndone rts
|
|
|
1268
|
|
|
1269 * exp(x)
|
|
|
1270 fpexp ldb -4,u
|
|
|
1271 stb ,-s ;Store sign of x.
|
|
|
1272 andb #$7f
|
|
|
1273 stb -4,u ;Take absolute value.
|
|
|
1274 ldx #fln2
|
|
|
1275 jsr fplod
|
|
|
1276 jsr fpmod ;modulo ln2.
|
|
|
1277 ldb #7
|
|
|
1278 ldy #expcoeff
|
|
|
1279 jsr fppoly ;Do the exp(-x) polynomial.
|
|
|
1280 jsr add1
|
|
|
1281 tst ,s+
|
|
|
1282 bpl exppos
|
|
|
1283 ldb -5,u ;Number was negative.
|
|
|
1284 subb intbuf+3 ;Subtract the integer quotient of the modln2
|
|
|
1285 bcs expund
|
|
|
1286 lda intbuf
|
|
|
1287 ora intbuf+1
|
|
|
1288 ora intbuf+2
|
|
|
1289 bne expund ;Underflow also if quotient >255
|
|
|
1290 stb -5,u ;Store exponent.
|
|
|
1291 rts
|
|
|
1292 exppos ldx #fpone
|
|
|
1293 jsr fplod
|
|
|
1294 jsr fpexg
|
|
|
1295 jsr fpdiv ;x was postitive, compute 1/exp(-x)
|
|
|
1296 ldb intbuf
|
|
|
1297 orb intbuf+1
|
|
|
1298 orb intbuf+2 ;Check int part is less than 255
|
|
|
1299 lbne fpovf
|
|
|
1300 ldb -5,u
|
|
|
1301 addb intbuf+3 ;Add integer part to exponent.
|
|
|
1302 lbcs fpovf ;Check for overflow.
|
|
|
1303 stb -5,u
|
|
|
1304 rts
|
|
|
1305 expund leau -5,u
|
|
|
1306 ldx #fpzero
|
|
|
1307 jmp fplod ;underflow, result is zero.
|
|
|
1308
|
|
|
1309 * ln(x) Natural logarithm
|
|
|
1310 fln jsr fptest0
|
|
|
1311 lbeq inval ;Don't accept zero as argument.
|
|
|
1312 tst -4,u
|
|
|
1313 lbmi inval ;No negative numbers either.
|
|
|
1314 ldb -5,u
|
|
|
1315 stb ,-s ;Save the binary exponent.
|
|
|
1316 ldb #$80
|
|
|
1317 stb -5,u ;Replace exponent with 1.
|
|
|
1318 ldx #fpone ;Argument is now in range 1..2
|
|
|
1319 jsr fplod
|
|
|
1320 jsr fpsub ;Subtract 1.
|
|
|
1321 ldy #lncoeff
|
|
|
1322 ldb #8
|
|
|
1323 jsr fppoly ;Do the ln(1+x) polynomial.
|
|
|
1324 ldb ,s+ ;Get original exponent.
|
|
|
1325 subb #$80 ;Unbias it.
|
|
|
1326 sex
|
|
|
1327 jsr int2fp ;Convert to fp.
|
|
|
1328 ldx #fln2
|
|
|
1329 jsr fplod
|
|
|
1330 jsr fpmul ;Multiply it by ln2.
|
|
|
1331 jmp fpadd ;Add that to result.
|
|
|
1332
|
|
|
1333 skipspace ldb ,y+
|
|
|
1334 cmpb #' '
|
|
|
1335 beq skipspace
|
|
|
1336 rts
|
|
|
1337
|
|
|
1338 PROGEND fdb $FFFF ;Indicate there is no AUTOBOOT app.
|
|
|
1339 ;Flag can be overwritten by it.
|
|
|
1340 FREEMEM equ ROM*RAMSTART+(1-ROM)*(PROGEND+2)
|
|
|
1341
|
|
|
1342
|
|
|
1343 end
|
|
|
1344
|
