Mercurial > hg > Members > kono > os9 > sbc09

comparison basic/fbasic.asm @ 57:2088fd998865

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
sbc09 directry clean up
author Shinji KONO <kono@ie.u-ryukyu.ac.jp>
date Mon, 23 Jul 2018 16:07:12 +0900
parents
children
comparison
equal deleted inserted replaced
56:4fa2bdb0c457 57:2088fd998865
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