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