| blob | 10dc2cd88a4abc11f7967b902e00f7a66e036764 |
1 {
2 $Id$
3 This file is part of the Free Pascal run time library.
4 Copyright (c) 1999-2000 by Carl-Eric Codere,
5 member of the Free Pascal development team
7 See the file COPYING.FPC, included in this distribution,
8 for details about the copyright.
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
14 **********************************************************************}
15 {*************************************************************************}
16 { lowmath.inc }
17 { Ported to FPK-Pascal by Carl Eric Codere }
18 { Terms of use: This source code is freeware. }
19 {*************************************************************************}
20 { This inc. implements low-level mathemtical routines for the motorola }
21 { 68000 family of processors. }
22 {*************************************************************************}
23 { single floating point routines taken from GCC 2.5.2 Atari compiler }
24 { library source. }
25 { Original credits: }
26 { written by Kai-Uwe Bloem (I5110401@dbstu1.bitnet). }
27 { Based on a 80x86 floating point packet from comp.os.minix, }
28 { written by P.Housel }
29 { Patched by Olaf Flebbe (flebbe@tat.physik.uni-tuebingen.de) }
30 { Revision by michal 05-93 (ntomczak@vm.ucs.ualberta.ca) }
31 {*************************************************************************}
32 {--------------------------------------------------------------------}
33 { LEFT TO DO: }
34 { o Add support for FPU if present. }
35 { o Verify if single comparison works in all cases. }
36 { o Add support for NANs in SINGLE_CMP }
37 { o Add comp (80-bit) multiplication,addition,substract,division, }
38 { shift. }
39 { o Add stack checking for the routines which use the stack. }
40 { (This will probably have to be done in the code generator). }
41 {--------------------------------------------------------------------}
45 Procedure Single_Norm;[alias : 'FPC_SINGLE_NORM'];Assembler;
46 {--------------------------------------------}
47 { Low-level routine to normalize single e }
48 { IEEE floating point values. Never called }
49 { directly. }
50 { On Exit: }
51 { d0 = result. }
52 { Registers destroyed: d0,d1 }
53 {--------------------------------------------}
54 Asm
55 tst.l d4 { rounding and u.mant == 0 ? }
56 bne @normlab1
57 tst.b d1
58 beq @retzok
59 @normlab1:
60 clr.b d2 { "sticky byte" }
61 @normlab3:
62 move.l #$ff000000,d5
63 @normlab4:
64 tst.w d0 { divide (shift) }
65 ble @normlab2 { denormalized number }
66 move.l d4,d3
67 and.l d5,d3 { or until no bits above 23 }
68 beq @normlab5
69 @normlab2:
70 addq.w #1,d0 { increment exponent }
71 lsr.l #1,d4
72 or.b d1,d2 { set "sticky" }
73 roxr.b #1,d1 { shift into rounding bits }
74 bra @normlab4
75 @normlab5:
76 and.b #1,d2
77 or.b d2,d1 { make least sig bit "sticky" }
78 asr.l #1,d5 { #0xff800000 -> d5 }
79 @normlab6:
80 move.l d4,d3 { multiply (shift) until }
81 and.l d5,d3 { one in "implied" position }
82 bne @normlab7
83 subq.w #1,d0 { decrement exponent }
84 beq @normlab7 { too small. store as denormalized number }
85 add.b d1,d1 { some doubt about this one * }
86 addx.l d4,d4
87 bra @normlab6
88 @normlab7:
89 tst.b d1 { check rounding bits }
90 bge @normlab9 { round down - no action neccessary }
91 neg.b d1
92 bvc @normlab8 { round up }
93 move.w d4,d1 { tie case - round to even }
94 { dont need rounding bits any more }
95 and.w #1,d1 { check if even }
96 beq @normlab9 { mantissa is even - no action necessary }
97 { fall through }
98 @normlab8:
99 clr.w d1 { zero rounding bits }
100 add.l #1,d4
101 tst.w d0
102 bne @normlab10 { renormalize if number was denormalized }
103 add.w #1,d0 { correct exponent for denormalized numbers }
104 bra @normlab3
105 @normlab10:
106 move.l d4,d3 { check for rounding overflow }
107 asl.l #1,d5 { #0xff000000 -> d5 }
108 and.l d5,d3
109 bne @normlab4 { go back and renormalize }
110 @normlab9:
111 tst.l d4 { check if normalization caused an underflow }
112 beq @retz
113 tst.w d0 { check for exponent overflow or underflow }
114 blt @retz
115 cmp.w #255,d0
116 bge @oflow
118 lsl.w #8,d0 { re-position exponent - one bit too high }
119 lsl.w #1,d2 { get X bit }
120 roxr.w #1,d0 { shift it into sign position }
121 swap d0 { map to upper word }
122 clr.w d0
123 and.l #$7fffff,d4 { top mantissa bits }
124 or.l d4,d0 { insert exponent and sign }
125 movem.l (sp)+,d2-d5
126 rts
128 @retz:
129 { handling underflow should be done here... }
130 { by default simply return 0 as retzok... }
131 @retzok:
132 moveq.l #0,d0
133 lsl.w #1,d2
134 roxr.l #1,d0 { sign of 0 is the same as of d2 }
135 movem.l (sp)+,d2-d5
136 rts
138 @oflow:
139 move.l #$7f800000,d0 { +infinity as proposed by IEEE }
141 tst.w d2 { transfer sign }
142 bge @ofl_clear { (mjr++) }
143 bset #31,d0 { }
144 @ofl_clear:
145 or.b #2,ccr { set overflow flag. }
146 movem.l (sp)+,d2-d5
147 rts
148 end;
151 Procedure Single_AddSub; Assembler;
152 {--------------------------------------------}
153 { Low-level routine to add/subtract single }
154 { IEEE floating point values. Never called }
155 { directly. }
156 { On Exit: }
157 { d0 = result -- from normalize routine }
158 { Flags : V set if overflow. }
159 { on underflow d0 = 0 }
160 { Registers destroyed: d0,d1 }
161 {--------------------------------------------}
162 Asm
163 {--------------------------------------------}
164 { On Entry: }
165 { d1-d0 = single values to subtract. }
166 {--------------------------------------------}
167 XDEF SINGLE_SUB
168 eor.l #$80000000,d0 { reverse sign of v }
169 {--------------------------------------------}
170 { On Entry: }
171 { d0, d1 = single values to add. }
172 {--------------------------------------------}
173 XDEF SINGLE_ADD
174 movem.l d2-d5,-(sp) { save registers }
175 move.l d0,d4 { d4 = d0 = v }
176 move.l d1,d5 { d5 = d1 = u }
178 move.l #$7fffff,d3
179 move.l d5,d0 { d0 = u.exp }
180 move.l d5,d2 { d2.h = u.sign }
181 swap d0
182 move.w d0,d2 { d2 = u.sign }
183 and.l d3,d5 { remove exponent from u.mantissa }
185 move.l d4,d1 { d1 = v.exp }
186 and.l d3,d4 { remove exponent from v.mantissa }
187 swap d1
188 eor.w d1,d2 { d2 = u.sign ^ v.sign (in bit 15)}
189 clr.b d2 { we will use the lowest byte as a flag }
190 moveq.l #15,d3
191 bclr d3,d1 { kill sign bit u.exp }
192 bclr d3,d0 { kill sign bit u.exp }
193 btst d3,d2 { same sign for u and v? }
194 beq @slabel1
195 cmp.l d0,d1 { different signs - maybe x - x ? }
196 seq d2 { set 'cancellation' flag }
197 @slabel1:
198 lsr.w #7,d0 { keep here exponents only }
199 lsr.w #7,d1
200 {--------------------------------------------------------------------}
201 { Now perform testing of NaN and infinities }
202 {--------------------------------------------------------------------}
203 moveq.l #-1,d3
204 cmp.b d3,d0
205 beq @alabel1
206 cmp.b d3,d1
207 bne @nospec
208 bra @alabel2
209 {--------------------------------------------------------------------}
210 { u is special. }
211 {--------------------------------------------------------------------}
212 @alabel1:
213 tst.b d2
214 bne @retnan { cancellation of specials -> NaN }
215 tst.l d5
216 bne @retnan { arith with Nan gives always NaN }
218 addq.w #4,a0 { here is an infinity }
219 cmp.b d3,d1
220 bne @alabel3 { skip check for NaN if v not special }
221 {--------------------------------------------------------------------}
222 { v is special. }
223 {--------------------------------------------------------------------}
224 @alabel2:
225 tst.l d4
226 bne @retnan
227 @alabel3:
228 move.l (a0),d0
229 bra @return
230 {--------------------------------------------------------------------}
231 { Return a quiet nan }
232 {--------------------------------------------------------------------}
233 @retnan:
234 moveq.l #-1,d0
235 lsr.l #1,d0 { 0x7fffffff -> d0 }
236 bra @return
237 { Ok, no inifinty or NaN involved.. }
238 @nospec:
239 tst.b d2
240 beq @alabel4
241 moveq.l #0,d0 { x - x hence we always return +0 }
242 @return:
243 movem.l (sp)+,d2-d5
244 rts
246 @alabel4:
247 moveq.l #23,d3
248 bset d3,d5 { restore implied leading "1" }
249 tst.w d0 { check for zero exponent - no leading "1" }
250 bne @alabel5
251 bclr d3,d5 { remove it }
252 addq.w #1,d0 { "normalize" exponent }
253 @alabel5:
254 bset d3,d4 { restore implied leading "1" }
255 tst.w d1 { check for zero exponent - no leading "1" }
256 bne @alabel6
257 bclr d3,d4 { remove it }
258 addq.w #1,d1 { "normalize" exponent }
259 @alabel6:
260 moveq.l #0,d3 { (put initial zero rounding bits in d3) }
261 neg.w d1 { d1 = u.exp - v.exp }
262 add.w d0,d1
263 beq @alabel8 { exponents are equal - no shifting neccessary }
264 bgt @alabel7 { not equal but no exchange neccessary }
265 exg d4,d5 { exchange u and v }
266 sub.w d1,d0 { d0 = u.exp - (u.exp - v.exp) = v.exp }
267 neg.w d1
268 tst.w d2 { d2.h = u.sign ^ (u.sign ^ v.sign) = v.sign }
269 bpl @alabel7
270 bchg #31,d2
271 @alabel7:
272 cmp.w #26,d1 { is u so much bigger that v is not }
273 bge @alabel9 { significant ? }
274 {--------------------------------------------------------------------}
275 { shift mantissa left two digits, to allow cancellation of }
276 { most significant digit, while gaining an additional digit for }
277 { rounding. }
278 {--------------------------------------------------------------------}
279 moveq.l #1,d3
280 @alabel10:
281 add.l d5,d5
282 subq.w #1,d0 { decrement exponent }
283 subq.w #1,d1 { done shifting altogether ? }
284 dbeq d3,@alabel10 { loop if still can shift u.mant more }
285 moveq.l #0,d3
287 cmp.w #16,d1 { see if fast rotate possible }
288 blt @alabel11
289 or.w d4,d3 { set rounding bits }
290 clr.w d4
291 swap d4
292 subq.w #8,d1
293 subq.w #8,d1
294 bra @alabel11
296 @alabel12:
297 move.b d4,d2
298 and.b #1,d2
299 or.b d2,d3
300 lsr.l #1,d4 { shift v.mant right the rest of the way }
301 @alabel11:
302 dbra d1,@alabel12 { loop }
304 @alabel8:
305 tst.w d2 { are the signs equal ? }
306 bpl @alabel13 { yes, no negate necessary }
309 tst.w d3 { negate rounding bits and v.mant }
310 beq @alabel14
311 addq.l #1,d4
312 @alabel14:
313 neg.l d4
315 @alabel13:
316 add.l d4,d5 { u.mant = u.mant + v.mant }
317 bcs @alabel9 { needn not negate }
318 tst.w d2 { opposite signs ? }
319 bpl @alabel9 { do not need to negate result }
321 neg.l d5
322 not.l d2 { switch sign }
323 @alabel9:
324 move.l d5,d4 { move result for normalization }
325 clr.l d1
326 tst.l d3
327 beq @alabel15
328 moveq.l #-1,d1
329 @alabel15:
330 swap d2 { put sign into d2 (exponent is in d0) }
331 jmp FPC_SINGLE_NORM { leave registers on stack for norm_sf }
332 end;
335 Procedure Single_Mul;Assembler;
336 {--------------------------------------------}
337 { Low-level routine to multiply two single }
338 { IEEE floating point values. Never called }
339 { directly. }
340 { Om Entry: }
341 { d0,d1 = values to multiply }
342 { On Exit: }
343 { d0 = result. }
344 { Registers destroyed: d0,d1 }
345 { stack space used (and restored): 8 bytes. }
346 {--------------------------------------------}
347 Asm
348 XDEF SINGLE_MUL
349 movem.l d2-d5,-(sp)
350 move.l d0,d4 { d4 = v }
351 move.l d1,d5 { d5 = u }
353 move.l #$7fffff,d3
354 move.l d5,d0 { d0 = u.exp }
355 and.l d3,d5 { remove exponent from u.mantissa }
356 swap d0
357 move.w d0,d2 { d2 = u.sign }
359 move.l d4,d1 { d1 = v.exp }
360 and.l d3,d4 { remove exponent from v.mantissa }
361 swap d1
362 eor.w d1,d2 { d2 = u.sign ^ v.sign (in bit 15)}
364 moveq.l #15,d3
365 bclr d3,d0 { kill sign bit }
366 bclr d3,d1 { kill sign bit }
367 tst.l d0 { test if one of factors is 0 }
368 beq @mlabel1
369 tst.l d1
370 @mlabel1:
371 seq d2 { 'one of factors is 0' flag in the lowest byte }
372 lsr.w #7,d0 { keep here exponents only }
373 lsr.w #7,d1
375 {--------------------------------------------------------------------}
376 { Now perform testing of NaN and infinities }
377 {--------------------------------------------------------------------}
378 moveq.l #-1,d3
379 cmp.b d3,d0
380 beq @mlabel2
381 cmp.b d3,d1
382 bne @mnospec
383 bra @mlabel3
384 {--------------------------------------------------------------------}
385 { first operand is special }
386 {--------------------------------------------------------------------}
387 @mlabel2:
388 tst.l d5 { is it NaN? }
389 bne @mretnan
390 @mlabel3:
391 tst.b d2 { 0 times special or special times 0 ? }
392 bne @mretnan { yes -> NaN }
393 cmp.b d3,d1 { is the other special ? }
394 beq @mlabel4 { maybe it is NaN }
395 {--------------------------------------------------------------------}
396 { Return infiny with correct sign }
397 {--------------------------------------------------------------------}
398 @mretinf:
399 move.l #$ff000000,d0 { we will return #0xff800000 or #0x7f800000 }
400 lsl.w #1,d2
401 roxr.l #1,d0 { shift in high bit as given by d2 }
402 @mreturn:
403 movem.l (sp)+,d2-d5
404 rts
406 {--------------------------------------------------------------------}
407 { v is special. }
408 {--------------------------------------------------------------------}
409 @mlabel4:
410 tst.l d4 { is this NaN? }
411 beq @mretinf { we know that the other is not zero }
412 @mretnan:
413 moveq.l #-1,d0
414 lsr.l #1,d0 { 0x7fffffff -> d0 }
415 bra @mreturn
416 {--------------------------------------------------------------------}
417 { End of NaN and Inf }
418 {--------------------------------------------------------------------}
419 @mnospec:
420 tst.b d2 { not needed - but we can waste two instr. }
421 bne @mretzz { return signed 0 if one of factors is 0 }
422 moveq.l #23,d3
423 bset d3,d5 { restore implied leading "1" }
424 subq.w #8,sp { multiplication accumulator }
425 tst.w d0 { check for zero exponent - no leading "1" }
426 bne @mlabel5
427 bclr d3,d5 { remove it }
428 addq.w #1,d0 { "normalize" exponent }
429 @mlabel5:
430 tst.l d5
431 beq @mretz { multiplying zero }
433 moveq.l #23,d3
434 bset d3,d4 { restore implied leading "1" }
435 tst.w d1 { check for zero exponent - no leading "1" }
436 bne @mlabel6
437 bclr d3,d4 { remove it }
438 addq.w #1,d1 { "normalize" exponent }
439 @mlabel6:
440 tst.l d4
441 beq @mretz { multiply by zero }
443 add.w d1,d0 { add exponents, }
444 sub.w #BIAS4+16-8,d0 { remove excess bias, acnt for repositioning }
446 clr.l (sp) { initialize 64-bit product to zero }
447 clr.l 4(sp)
448 {--------------------------------------------------------------------}
449 { see Knuth, Seminumerical Algorithms, section 4.3. algorithm M }
450 {--------------------------------------------------------------------}
451 move.w d4,d3
452 mulu.w d5,d3 { mulitply with bigit from multiplier }
453 move.l d3,4(sp) { store into result }
455 move.l d4,d3
456 swap d3
457 mulu.w d5,d3
458 add.l d3,2(sp) { add to result }
460 swap d5 { [TOP 8 BITS SHOULD BE ZERO !] }
462 move.w d4,d3
463 mulu.w d5,d3 { mulitply with bigit from multiplier }
464 add.l d3,2(sp) { store into result (no carry can occur here) }
466 move.l d4,d3
467 swap d3
468 mulu.w d5,d3
469 add.l d3,(sp) { add to result }
470 { [TOP 16 BITS SHOULD BE ZERO !] }
471 movem.l 2(sp),d4-d5 { get the 48 valid mantissa bits }
472 clr.w d5 { (pad to 64) }
474 move.l #$0000ffff,d3
475 @mlabel7:
476 cmp.l d3,d4 { multiply (shift) until }
477 bhi @mlabel8 { 1 in upper 16 result bits }
478 cmp.w #9,d0 { give up for denormalized numbers }
479 ble @mlabel8
480 swap d4 { (we''re getting here only when multiplying }
481 swap d5 { with a denormalized number; there''s an }
482 move.w d5,d4 { eventual loss of 4 bits in the rounding }
483 clr.w d5 { byte -- what a pity 8-) }
484 subq.w #8,d0 { decrement exponent }
485 subq.w #8,d0
486 bra @mlabel7
487 @mlabel8:
488 move.l d5,d1 { get rounding bits }
489 rol.l #8,d1
490 move.l d1,d3 { see if sticky bit should be set }
491 and.l #$ffffff00,d3
492 beq @mlabel9
493 or.b #1,d1 { set "sticky bit" if any low-order set }
494 @mlabel9:
495 addq.w #8,sp { remove accumulator from stack }
496 jmp FPC_SINGLE_NORM{ (result in d4) }
498 @mretz:
499 addq.w #8,sp { release accumulator space }
500 @mretzz:
501 moveq.l #0,d0 { save zero as result }
502 lsl.w #1,d2 { and set it sign as for d2 }
503 roxr.l #1,d0
504 movem.l (sp)+,d2-d5
505 rts { no normalizing neccessary }
506 end;
509 Procedure Single_Div;Assembler;
510 {--------------------------------------------}
511 { Low-level routine to dividr two single }
512 { IEEE floating point values. Never called }
513 { directly. }
514 { Om Entry: }
515 { d1/d0 = u/v = operation to perform. }
516 { On Exit: }
517 { d0 = result. }
518 { Registers destroyed: d0,d1 }
519 { stack space used (and restored): 8 bytes. }
520 {--------------------------------------------}
521 ASM
522 XDEF SINGLE_DIV
523 { u = d1 = dividend }
524 { v = d0 = divisor }
525 tst.l d0 { check if divisor is 0 }
526 bne @dno_exception
528 move.l #$7f800000,d0
529 btst #31,d1 { transfer sign of dividend }
530 beq @dclear
531 bset #31,d0
532 @dclear:
533 rts
535 @dno_exception:
536 move.l d1,d4 { d4 = u, d5 = v }
537 move.l d0,d5
538 movem.l d2-d5,-(sp) { save registers }
540 move.l #$7fffff,d3
541 move.l d4,d0 { d0 = u.exp }
542 and.l d3,d4 { remove exponent from u.mantissa }
543 swap d0
544 move.w d0,d2 { d2 = u.sign }
546 move.l d5,d1 { d1 = v.exp }
547 and.l d3,d5 { remove exponent from v.mantissa }
548 swap d1
549 eor.w d1,d2 { d2 = u.sign ^ v.sign (in bit 15) }
551 moveq.l #15,d3
552 bclr d3,d0 { kill sign bit }
553 bclr d3,d1 { kill sign bit }
554 lsr.w #7,d0
555 lsr.w #7,d1
557 moveq.l #-1,d3
558 cmp.b d3,d0 { comparison with #0xff }
559 beq @dlabel1 { u == NaN ;; u== Inf }
560 cmp.b d3,d1
561 beq @dlabel2 { v == NaN ;; v == Inf }
562 tst.b d0
563 bne @dlabel4 { u not zero nor denorm }
564 tst.l d4
565 beq @dlabel3 { 0/ ? }
567 @dlabel4:
568 tst.w d1
569 bne @dnospec
571 tst.l d5
572 bne @dnospec
573 bra @dretinf { x/0 -> +/- Inf }
575 @dlabel1:
576 tst.l d4 { u == NaN ? }
577 bne @dretnan { NaN/ x }
578 cmp.b d3,d1
579 beq @dretnan { Inf/Inf or Inf/NaN }
580 { bra dretinf ; Inf/x ; x != Inf && x != NaN }
581 {--------------------------------------------------------------------}
582 { Return infinity with correct sign. }
583 {--------------------------------------------------------------------}
584 @dretinf:
585 move.l #$ff000000,d0
586 lsl.w #1,d2
587 roxr.l #1,d0 { shift in high bit as given by d2 }
588 @dreturn:
589 movem.l (sp)+,d2-d5
590 rts
592 @dlabel2:
593 tst.l d5
594 bne @dretnan { x/NaN }
595 { bra dretzero ; x/Inf -> +/- 0 }
596 {--------------------------------------------------------------------}
597 { Return correct signed zero. }
598 {--------------------------------------------------------------------}
599 @dretzero:
600 moveq.l #0,d0 { zero destination }
601 lsl.w #1,d2 { set X bit accordingly }
602 roxr.l #1,d0
603 bra @dreturn
605 @dlabel3:
606 tst.w d1
607 bne @dretzero { 0/x ->+/- 0 }
608 tst.l d4
609 bne @dretzero { 0/x }
610 { bra dretnan 0/0 }
611 {--------------------------------------------------------------------}
612 { Return NotANumber }
613 {--------------------------------------------------------------------}
614 @dretnan:
615 move.l d3,d0 { d3 contains 0xffffffff }
616 lsr.l #1,d0
617 bra @dreturn
618 {--------------------------------------------------------------------}
619 { End of Special Handling }
620 {--------------------------------------------------------------------}
621 @dnospec:
622 moveq.l #23,d3
623 bset d3,d4 { restore implied leading "1" }
624 tst.w d0 { check for zero exponent - no leading "1" }
625 bne @dlabel5
626 bclr d3,d4 { remove it }
627 add.w #1,d0 { "normalize" exponent }
628 @dlabel5:
629 tst.l d4
630 beq @dretzero { dividing zero }
632 bset d3,d5 { restore implied leading "1" }
633 tst.w d1 { check for zero exponent - no leading "1"}
634 bne @dlabel6
635 bclr d3,d5 { remove it }
636 add.w #1,d1 { "normalize" exponent }
637 @dlabel6:
639 sub.w d1,d0 { subtract exponents, }
640 add.w #BIAS4-8+1,d0 { add bias back in, account for shift }
641 add.w #34,d0 { add loop offset, +2 for extra rounding bits}
642 { for denormalized numbers (2 implied by dbra)}
643 move.w #27,d1 { bit number for "implied" pos (+4 for rounding)}
644 moveq.l #-1,d3 { zero quotient (for speed a one''s complement) }
645 sub.l d5,d4 { initial subtraction, u = u - v }
646 @dlabel7:
647 btst d1,d3 { divide until 1 in implied position }
648 beq @dlabel9
650 add.l d4,d4
651 bcs @dlabel8 { if carry is set, add, else subtract }
653 addx.l d3,d3 { shift quotient and set bit zero }
654 sub.l d5,d4 { subtract, u = u - v }
655 dbra d0,@dlabel7 { give up if result is denormalized }
656 bra @dlabel9
657 @dlabel8:
658 addx.l d3,d3 { shift quotient and clear bit zero }
659 add.l d5,d4 { add (restore), u = u + v }
660 dbra d0,@dlabel7 { give up if result is denormalized }
661 @dlabel9:
662 subq.w #2,d0 { remove rounding offset for denormalized nums }
663 not.l d3 { invert quotient to get it right }
665 clr.l d1 { zero rounding bits }
666 tst.l d4 { check for exact result }
667 beq @dlabel10
668 moveq.l #-1,d1 { prevent tie case }
669 @dlabel10:
670 move.l d3,d4 { save quotient mantissa }
671 jmp FPC_SINGLE_NORM{ (registers on stack removed by norm_sf) }
672 end;
675 Procedure Single_Cmp; Assembler;
676 {--------------------------------------------}
677 { Low-level routine to compare single two }
678 { single point values.. }
679 { Never called directly. }
680 { On Entry: }
681 { d1 and d0 Values to compare }
682 { d0 = first operand }
683 { On Exit: }
684 { Flags according to result }
685 { Registers destroyed: d0,d1 }
686 {--------------------------------------------}
687 Asm
688 XDEF SINGLE_CMP
689 tst.l d0 { check sign bit }
690 bpl @cmplab1
691 neg.l d0 { negate }
692 bchg #31,d0 { toggle sign bit }
693 @cmplab1:
694 tst.l d1 { check sign bit }
695 bpl @cmplab2
696 neg.l d1 { negate }
697 bchg #31,d1 { toggle sign bit }
698 @cmplab2:
699 cmp.l d0,d1 { compare... }
700 rts
701 end;
705 Procedure LongMul;Assembler;
706 {--------------------------------------------}
707 { Low-level routine to multiply two signed }
708 { 32-bit values. Never called directly. }
709 { On entry: d1,d0 = 32-bit signed values to }
710 { multiply. }
711 { On Exit: }
712 { d0 = result. }
713 { Registers destroyed: d0,d1 }
714 { stack space used and restored: 10 bytes }
715 {--------------------------------------------}
716 Asm
717 XDEF LONGMUL
718 cmp.b #2,Test68000 { Are we on a 68020+ cpu }
719 blt @Lmulcontinue
720 muls.l d1,d0 { yes, then directly mul... }
721 rts { return... result in d0 }
722 @Lmulcontinue:
723 move.l d2,a0 { save registers }
724 move.l d3,a1
726 move.l d0,-(sp)
727 move.l d1,-(sp)
729 movem.w (sp)+,d0-d3 { u = d0-d1, v = d2-d3 }
731 move.w d0,-(sp) { sign flag }
732 bpl @LMul1 { is u negative ? }
733 neg.w d1 { yes, force it positive }
734 negx.w d0
735 @LMul1:
736 tst.w d2 { is v negative ? }
737 bpl @LMul2
738 neg.w d3 { yes, force it positive ... }
739 negx.w d2
740 not.w (sp) { ... and modify flag word }
741 @LMul2:
742 ext.l d0 { u.h <> 0 ? }
743 beq @LMul3
744 mulu.w d3,d0 { r = v.l * u.h }
745 @LMul3:
746 tst.w d2 { v.h <> 0 ? }
747 beq @LMul4
748 mulu.w d1,d2 { r += v.h * u.l }
749 add.w d2,d0
750 @LMul4:
751 swap d0
752 clr.w d0
753 mulu.w d3,d1 { r += v.l * u.l }
754 add.l d1,d0
755 move.l a1,d3
756 move.l a0,d2
757 tst.w (sp)+ { should the result be negated ? }
758 bpl @LMul5 { no, just return }
759 neg.l d0 { else r = -r }
760 @LMul5:
761 rts
762 end;
766 Procedure Long2Single;Assembler;
767 {--------------------------------------------}
768 { Low-level routine to convert a longint }
769 { to a single floating point value. }
770 { On entry: d0 = longint value to convert. }
771 { On Exit: }
772 { d0 = single IEEE value }
773 { Registers destroyed: d0,d1 }
774 { stack space used and restored: 8 bytes }
775 {--------------------------------------------}
776 Asm
777 XDEF LONG2SINGLE
778 movem.l d2-d5,-(sp) { save registers to make norm_sf happy}
780 move.l d0,d4 { prepare result mantissa }
781 move.w #BIAS4+32-8,d0 { radix point after 32 bits }
782 move.l d4,d2 { set sign flag }
783 bge @l2slabel1 { nonnegative }
784 neg.l d4 { take absolute value }
785 @l2slabel1:
786 swap d2 { follow SINGLE_NORM conventions }
787 clr.w d1 { set rounding = 0 }
788 jmp FPC_SINGLE_NORM
789 end;
792 Procedure LongDiv; [alias : 'FPC_LONGDIV'];Assembler;
793 {--------------------------------------------}
794 { Low-level routine to do signed long }
795 { division. }
796 { On entry: d0/d1 operation to perform }
797 { On Exit: }
798 { d0 = quotient }
799 { d1 = remainder }
800 { Registers destroyed: d0,d1,d6 }
801 { stack space used and restored: 10 bytes }
802 {--------------------------------------------}
803 asm
804 XDEF LONGDIV
805 cmp.b #2,Test68000 { can we use divs ? }
806 blt @continue
807 tst.l d1
808 beq @zerodiv2
809 move.l d1,d6
810 clr.l d1 { clr }
811 tst.l d0 { check sign of d0 }
812 bpl @posdiv
813 move.l #$ffffffff,d1{ sign extend into d1 }
814 @posdiv:
815 divsl.l d6,d1:d0
816 rts
817 @continue:
819 move.l d2,a0 { save registers }
820 move.l d3,a1
821 move.l d4,-(sp) { divisor = d1 = d4 }
822 move.l d5,-(sp) { divident = d0 = d5 }
824 move.l d1,d4 { save divisor }
825 move.l d0,d5 { save dividend }
827 clr.w -(sp) { sign flag }
829 clr.l d0 { prepare result }
830 move.l d4,d2 { get divisor }
831 beq @zerodiv { divisor = 0 ? }
832 bpl @LDiv1 { divisor < 0 ? }
833 neg.l d2 { negate it }
834 not.w (sp) { remember sign }
835 @LDiv1:
836 move.l d5,d1 { get dividend }
837 bpl @LDiv2 { dividend < 0 ? }
838 neg.l d1 { negate it }
839 not.w (sp) { remember sign }
840 @LDiv2:
841 {;== case 1) divident < divisor}
842 cmp.l d2,d1 { is divident smaller then divisor ? }
843 bcs @LDiv7 { yes, return immediately }
844 {;== case 2) divisor has <= 16 significant bits}
845 move.l d4,d6 { put divisor in d6 register }
846 lsr.l #8,d6 { rotate into low word }
847 lsr.l #8,d6
848 tst.l d6
849 bne @LDiv3 { divisor has only 16 bits }
850 move.w d1,d3 { save dividend }
851 clr.w d1 { divide dvd.h by dvs }
852 swap d1
853 beq @LDiv4 { (no division necessary if dividend zero)}
854 divu d2,d1
855 @LDiv4:
856 move.w d1,d0 { save quotient.h }
857 swap d0
858 move.w d3,d1 { (d0.h = remainder of prev divu) }
859 divu d2,d1 { divide dvd.l by dvs }
860 move.w d1,d0 { save quotient.l }
861 clr.w d1 { get remainder }
862 swap d1
863 bra @LDiv7 { and return }
864 {;== case 3) divisor > 16 bits (corollary is dividend > 16 bits, see case 1)}
865 @LDiv3:
866 moveq.l #31,d3 { loop count }
867 @LDiv5:
868 add.l d1,d1 { shift divident ... }
869 addx.l d0,d0 { ... into d0 }
870 cmp.l d2,d0 { compare with divisor }
871 bcs @LDiv6
872 sub.l d2,d0 { big enough, subtract }
873 addq.w #1,d1 { and note bit into result }
874 @LDiv6:
875 dbra d3,@LDiv5
876 exg d0,d1 { put quotient and remainder in their registers}
877 @LDiv7:
878 tst.l d5 { must the remainder be corrected ? }
879 bpl @LDiv8
880 neg.l d1 { yes, apply sign }
881 { the following line would be correct if modulus is defined as in algebra}
882 {; add.l sp@(6),d1 ; algebraic correction: modulus can only be >= 0}
883 @LDiv8:
884 tst.w (sp)+ { result should be negative ? }
885 bpl @LDiv9
886 neg.l d0 { yes, negate it }
887 @LDiv9:
888 move.l a1,d3
889 move.l a0,d2
890 move.l (sp)+,d5
891 move.l (sp)+,d4
892 rts { en exit : remainder = d1, quotient = d0 }
893 @zerodiv:
894 move.l a1,d3 { restore stack... }
895 move.l a0,d2
896 move.w (sp)+,d1 { remove sign word }
897 move.l (sp)+,d5
898 move.l (sp)+,d4
899 @zerodiv2:
900 move.l #200,d0
901 jsr FPC_HALT_ERROR { RunError(200) }
902 rts { this should never occur... }
903 end;
906 Procedure LongMod;[alias : 'FPC_LONGMOD'];Assembler;
907 { see longdiv for info on calling convention }
908 asm
909 XDEF LONGMOD
910 jsr FPC_LONGDIV
911 move.l d1,d0 { return the remainder in d0 }
912 rts
913 end;
915 {
916 $Log$
917 Revision 1.1 2002/02/19 08:25:41 sasu
918 Initial revision
920 Revision 1.1 2000/07/13 06:30:55 michael
921 + Initial import
923 Revision 1.6 2000/01/07 16:41:42 daniel
924 * copyright 2000
926 Revision 1.5 2000/01/07 16:32:28 daniel
927 * copyright 2000 added
929 Revision 1.4 1998/10/13 08:00:04 pierre
930 * some bugs related to FPC_ prefix fixed
931 * problems with pbyte sometimes defined and sometimes not for rttip.inc solved
933 Revision 1.3 1998/07/01 14:28:32 carl
934 * LONGDIV bugfixed with signed and modulo
935 * LONGDIV bugfix when divisor is less then 16 bits
937 Revision 1.1.1.1 1998/03/25 11:18:44 root
938 * Restored version
940 Revision 1.5 1998/01/26 12:01:47 michael
941 + Added log at the end
945 Working file: rtl/m68k/lowmath.inc
946 description:
947 ----------------------------
948 revision 1.4
949 date: 1998/01/05 00:33:30; author: carl; state: Exp; lines: +144 -156
950 * Several minor bugfixes
951 ----------------------------
952 revision 1.3
953 date: 1997/12/14 19:03:36; author: carl; state: Exp; lines: +11 -12
954 + uses now correct register conventions
955 ----------------------------
956 revision 1.2
957 date: 1997/12/01 12:37:20; author: michael; state: Exp; lines: +14 -0
958 + added copyright reference in header.
959 ----------------------------
960 revision 1.1
961 date: 1997/11/27 13:55:47; author: carl; state: Exp;
962 Internal RTL math functions for m68k.
963 =============================================================================
964 }