gcc/config/m68k/lb1sf68.asm

   1 /* libgcc routines for 68000 w/o floating-point hardware.
   2    Copyright (C) 1994, 1996, 1997, 1998 Free Software Foundation, Inc.
   3
   4 This file is part of GNU CC.
   5
   6 GNU CC is free software; you can redistribute it and/or modify it
   7 under the terms of the GNU General Public License as published by the
   8 Free Software Foundation; either version 2, or (at your option) any
   9 later version.
  10
  11 In addition to the permissions in the GNU General Public License, the
  12 Free Software Foundation gives you unlimited permission to link the
  13 compiled version of this file with other programs, and to distribute
  14 those programs without any restriction coming from the use of this
  15 file.  (The General Public License restrictions do apply in other
  16 respects; for example, they cover modification of the file, and
  17 distribution when not linked into another program.)
  18
  19 This file is distributed in the hope that it will be useful, but
  20 WITHOUT ANY WARRANTY; without even the implied warranty of
  21 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  22 General Public License for more details.
  23
  24 You should have received a copy of the GNU General Public License
  25 along with this program; see the file COPYING.  If not, write to
  26 the Free Software Foundation, 59 Temple Place - Suite 330,
  27 Boston, MA 02111-1307, USA.  */
  28
  29 /* As a special exception, if you link this library with files
  30    compiled with GCC to produce an executable, this does not cause
  31    the resulting executable to be covered by the GNU General Public License.
  32    This exception does not however invalidate any other reasons why
  33    the executable file might be covered by the GNU General Public License.  */
  34
  35 /* Use this one for any 680x0; assumes no floating point hardware.
  36    The trailing " '" appearing on some lines is for ANSI preprocessors.  Yuk.
  37    Some of this code comes from MINIX, via the folks at ericsson.
  38    D. V. Henkel-Wallace (gumby@cygnus.com) Fete Bastille, 1992
  39 */
  40
  41 /* These are predefined by new versions of GNU cpp.  */
  42
  43 #ifndef __USER_LABEL_PREFIX__
  44 #define __USER_LABEL_PREFIX__ _
  45 #endif
  46
  47 #ifndef __REGISTER_PREFIX__
  48 #define __REGISTER_PREFIX__
  49 #endif
  50
  51 #ifndef __IMMEDIATE_PREFIX__
  52 #define __IMMEDIATE_PREFIX__ #
  53 #endif
  54
  55 /* ANSI concatenation macros.  */
  56
  57 #define CONCAT1(a, b) CONCAT2(a, b)
  58 #define CONCAT2(a, b) a ## b
  59
  60 /* Use the right prefix for global labels.  */
  61
  62 #define SYM(x) CONCAT1 (__USER_LABEL_PREFIX__, x)
  63
  64 /* Use the right prefix for registers.  */
  65
  66 #define REG(x) CONCAT1 (__REGISTER_PREFIX__, x)
  67
  68 /* Use the right prefix for immediate values.  */
  69
  70 #define IMM(x) CONCAT1 (__IMMEDIATE_PREFIX__, x)
  71
  72 #define d0 REG (d0)
  73 #define d1 REG (d1)
  74 #define d2 REG (d2)
  75 #define d3 REG (d3)
  76 #define d4 REG (d4)
  77 #define d5 REG (d5)
  78 #define d6 REG (d6)
  79 #define d7 REG (d7)
  80 #define a0 REG (a0)
  81 #define a1 REG (a1)
  82 #define a2 REG (a2)
  83 #define a3 REG (a3)
  84 #define a4 REG (a4)
  85 #define a5 REG (a5)
  86 #define a6 REG (a6)
  87 #define fp REG (fp)
  88 #define sp REG (sp)
  89
  90 #ifdef L_floatex
  91
  92 | This is an attempt at a decent floating point (single, double and
  93 | extended double) code for the GNU C compiler. It should be easy to
  94 | adapt to other compilers (but beware of the local labels!).
  95
  96 | Starting date: 21 October, 1990
  97
  98 | It is convenient to introduce the notation (s,e,f) for a floating point
  99 | number, where s=sign, e=exponent, f=fraction. We will call a floating
 100 | point number fpn to abbreviate, independently of the precision.
 101 | Let MAX_EXP be in each case the maximum exponent (255 for floats, 1023
 102 | for doubles and 16383 for long doubles). We then have the following
 103 | different cases:
 104 |  1. Normalized fpns have 0 < e < MAX_EXP. They correspond to
 105 |     (-1)^s x 1.f x 2^(e-bias-1).
 106 |  2. Denormalized fpns have e=0. They correspond to numbers of the form
 107 |     (-1)^s x 0.f x 2^(-bias).
 108 |  3. +/-INFINITY have e=MAX_EXP, f=0.
 109 |  4. Quiet NaN (Not a Number) have all bits set.
 110 |  5. Signaling NaN (Not a Number) have s=0, e=MAX_EXP, f=1.
 111
 112 |=============================================================================
 113 |                                  exceptions
 114 |=============================================================================
 115
 116 | This is the floating point condition code register (_fpCCR):
 117 |
 118 | struct {
 119 |   short _exception_bits;
 120 |   short _trap_enable_bits;
 121 |   short _sticky_bits;
 122 |   short _rounding_mode;
 123 |   short _format;
 124 |   short _last_operation;
 125 |   union {
 126 |     float sf;
 127 |     double df;
 128 |   } _operand1;
 129 |   union {
 130 |     float sf;
 131 |     double df;
 132 |   } _operand2;
 133 | } _fpCCR;
 134
 135         .data
 136         .even
 137
 138         .globl  SYM (_fpCCR)
 139
 140 SYM (_fpCCR):
 141 __exception_bits:
 142         .word   0
 143 __trap_enable_bits:
 144         .word   0
 145 __sticky_bits:
 146         .word   0
 147 __rounding_mode:
 148         .word   ROUND_TO_NEAREST
 149 __format:
 150         .word   NIL
 151 __last_operation:
 152         .word   NOOP
 153 __operand1:
 154         .long   0
 155         .long   0
 156 __operand2:
 157         .long   0
 158         .long   0
 159
 160 | Offsets:
 161 EBITS  = __exception_bits - SYM (_fpCCR)
 162 TRAPE  = __trap_enable_bits - SYM (_fpCCR)
 163 STICK  = __sticky_bits - SYM (_fpCCR)
 164 ROUND  = __rounding_mode - SYM (_fpCCR)
 165 FORMT  = __format - SYM (_fpCCR)
 166 LASTO  = __last_operation - SYM (_fpCCR)
 167 OPER1  = __operand1 - SYM (_fpCCR)
 168 OPER2  = __operand2 - SYM (_fpCCR)
 169
 170 | The following exception types are supported:
 171 INEXACT_RESULT          = 0x0001
 172 UNDERFLOW               = 0x0002
 173 OVERFLOW                = 0x0004
 174 DIVIDE_BY_ZERO          = 0x0008
 175 INVALID_OPERATION       = 0x0010
 176
 177 | The allowed rounding modes are:
 178 UNKNOWN           = -1
 179 ROUND_TO_NEAREST  = 0 | round result to nearest representable value
 180 ROUND_TO_ZERO     = 1 | round result towards zero
 181 ROUND_TO_PLUS     = 2 | round result towards plus infinity
 182 ROUND_TO_MINUS    = 3 | round result towards minus infinity
 183
 184 | The allowed values of format are:
 185 NIL          = 0
 186 SINGLE_FLOAT = 1
 187 DOUBLE_FLOAT = 2
 188 LONG_FLOAT   = 3
 189
 190 | The allowed values for the last operation are:
 191 NOOP         = 0
 192 ADD          = 1
 193 MULTIPLY     = 2
 194 DIVIDE       = 3
 195 NEGATE       = 4
 196 COMPARE      = 5
 197 EXTENDSFDF   = 6
 198 TRUNCDFSF    = 7
 199
 200 |=============================================================================
 201 |                           __clear_sticky_bits
 202 |=============================================================================
 203
 204 | The sticky bits are normally not cleared (thus the name), whereas the
 205 | exception type and exception value reflect the last computation.
 206 | This routine is provided to clear them (you can also write to _fpCCR,
 207 | since it is globally visible).
 208
 209         .globl  SYM (__clear_sticky_bit)
 210
 211         .text
 212         .even
 213
 214 | void __clear_sticky_bits(void);
 215 SYM (__clear_sticky_bit):
 216         lea     SYM (_fpCCR),a0
 217 #ifndef __mcoldfire__
 218         movew   IMM (0),a0@(STICK)
 219 #else
 220         clr.w   a0@(STICK)
 221 #endif
 222         rts
 223
 224 |=============================================================================
 225 |                           $_exception_handler
 226 |=============================================================================
 227
 228         .globl  $_exception_handler
 229
 230         .text
 231         .even
 232
 233 | This is the common exit point if an exception occurs.
 234 | NOTE: it is NOT callable from C!
 235 | It expects the exception type in d7, the format (SINGLE_FLOAT,
 236 | DOUBLE_FLOAT or LONG_FLOAT) in d6, and the last operation code in d5.
 237 | It sets the corresponding exception and sticky bits, and the format.
 238 | Depending on the format if fills the corresponding slots for the
 239 | operands which produced the exception (all this information is provided
 240 | so if you write your own exception handlers you have enough information
 241 | to deal with the problem).
 242 | Then checks to see if the corresponding exception is trap-enabled,
 243 | in which case it pushes the address of _fpCCR and traps through
 244 | trap FPTRAP (15 for the moment).
 245
 246 FPTRAP = 15
 247
 248 $_exception_handler:
 249         lea     SYM (_fpCCR),a0
 250         movew   d7,a0@(EBITS)   | set __exception_bits
 251 #ifndef __mcoldfire__
 252         orw     d7,a0@(STICK)   | and __sticky_bits
 253 #else
 254         movew   a0@(STICK),d4
 255         orl     d7,d4
 256         movew   d4,a0@(STICK)
 257 #endif
 258         movew   d6,a0@(FORMT)   | and __format
 259         movew   d5,a0@(LASTO)   | and __last_operation
 260
 261 | Now put the operands in place:
 262 #ifndef __mcoldfire__
 263         cmpw    IMM (SINGLE_FLOAT),d6
 264 #else
 265         cmpl    IMM (SINGLE_FLOAT),d6
 266 #endif
 267         beq     1f
 268         movel   a6@(8),a0@(OPER1)
 269         movel   a6@(12),a0@(OPER1+4)
 270         movel   a6@(16),a0@(OPER2)
 271         movel   a6@(20),a0@(OPER2+4)
 272         bra     2f
 273 1:      movel   a6@(8),a0@(OPER1)
 274         movel   a6@(12),a0@(OPER2)
 275 2:
 276 | And check whether the exception is trap-enabled:
 277 #ifndef __mcoldfire__
 278         andw    a0@(TRAPE),d7   | is exception trap-enabled?
 279 #else
 280         clrl    d6
 281         movew   a0@(TRAPE),d6
 282         andl    d6,d7
 283 #endif
 284         beq     1f              | no, exit
 285         pea     SYM (_fpCCR)    | yes, push address of _fpCCR
 286         trap    IMM (FPTRAP)    | and trap
 287 #ifndef __mcoldfire__
 288 1:      moveml  sp@+,d2-d7      | restore data registers
 289 #else
 290 1:      moveml  sp@,d2-d7
 291         | XXX if frame pointer is ever removed, stack pointer must
 292         | be adjusted here.
 293 #endif
 294         unlk    a6              | and return
 295         rts
 296 #endif /* L_floatex */
 297
 298 #ifdef  L_mulsi3
 299         .text
 300         .proc
 301         .globl  SYM (__mulsi3)
 302 SYM (__mulsi3):
 303         movew   sp@(4), d0      /* x0 -> d0 */
 304         muluw   sp@(10), d0     /* x0*y1 */
 305         movew   sp@(6), d1      /* x1 -> d1 */
 306         muluw   sp@(8), d1      /* x1*y0 */
 307 #ifndef __mcoldfire__
 308         addw    d1, d0
 309 #else
 310         addl    d1, d0
 311 #endif
 312         swap    d0
 313         clrw    d0
 314         movew   sp@(6), d1      /* x1 -> d1 */
 315         muluw   sp@(10), d1     /* x1*y1 */
 316         addl    d1, d0
 317
 318         rts
 319 #endif /* L_mulsi3 */
 320
 321 #ifdef  L_udivsi3
 322         .text
 323         .proc
 324         .globl  SYM (__udivsi3)
 325 SYM (__udivsi3):
 326 #ifndef __mcoldfire__
 327         movel   d2, sp@-
 328         movel   sp@(12), d1     /* d1 = divisor */
 329         movel   sp@(8), d0      /* d0 = dividend */
 330
 331         cmpl    IMM (0x10000), d1 /* divisor >= 2 ^ 16 ?   */
 332         jcc     L3              /* then try next algorithm */
 333         movel   d0, d2
 334         clrw    d2
 335         swap    d2
 336         divu    d1, d2          /* high quotient in lower word */
 337         movew   d2, d0          /* save high quotient */
 338         swap    d0
 339         movew   sp@(10), d2     /* get low dividend + high rest */
 340         divu    d1, d2          /* low quotient */
 341         movew   d2, d0
 342         jra     L6
 343
 344 L3:     movel   d1, d2          /* use d2 as divisor backup */
 345 L4:     lsrl    IMM (1), d1     /* shift divisor */
 346         lsrl    IMM (1), d0     /* shift dividend */
 347         cmpl    IMM (0x10000), d1 /* still divisor >= 2 ^ 16 ?  */
 348         jcc     L4
 349         divu    d1, d0          /* now we have 16 bit divisor */
 350         andl    IMM (0xffff), d0 /* mask out divisor, ignore remainder */
 351
 352 /* Multiply the 16 bit tentative quotient with the 32 bit divisor.  Because of
 353    the operand ranges, this might give a 33 bit product.  If this product is
 354    greater than the dividend, the tentative quotient was too large. */
 355         movel   d2, d1
 356         mulu    d0, d1          /* low part, 32 bits */
 357         swap    d2
 358         mulu    d0, d2          /* high part, at most 17 bits */
 359         swap    d2              /* align high part with low part */
 360         tstw    d2              /* high part 17 bits? */
 361         jne     L5              /* if 17 bits, quotient was too large */
 362         addl    d2, d1          /* add parts */
 363         jcs     L5              /* if sum is 33 bits, quotient was too large */
 364         cmpl    sp@(8), d1      /* compare the sum with the dividend */
 365         jls     L6              /* if sum > dividend, quotient was too large */
 366 L5:     subql   IMM (1), d0     /* adjust quotient */
 367
 368 L6:     movel   sp@+, d2
 369         rts
 370
 371 #else /* __mcoldfire__ */
 372
 373 /* Coldfire implementation of non-restoring division algorithm from
 374    Hennessy & Patterson, Appendix A. */
 375         link    a6,IMM (-12)
 376         moveml  d2-d4,sp@
 377         movel   a6@(8),d0
 378         movel   a6@(12),d1
 379         clrl    d2              | clear p
 380         moveq   IMM (31),d4
 381 L1:     addl    d0,d0           | shift reg pair (p,a) one bit left
 382         addxl   d2,d2
 383         movl    d2,d3           | subtract b from p, store in tmp.
 384         subl    d1,d3
 385         jcs     L2              | if no carry,
 386         bset    IMM (0),d0      | set the low order bit of a to 1,
 387         movl    d3,d2           | and store tmp in p.
 388 L2:     subql   IMM (1),d4
 389         jcc     L1
 390         moveml  sp@,d2-d4       | restore data registers
 391         unlk    a6              | and return
 392         rts
 393 #endif /* __mcoldfire__ */
 394
 395 #endif /* L_udivsi3 */
 396
 397 #ifdef  L_divsi3
 398         .text
 399         .proc
 400         .globl  SYM (__divsi3)
 401 SYM (__divsi3):
 402         movel   d2, sp@-
 403
 404         moveq   IMM (1), d2     /* sign of result stored in d2 (=1 or =-1) */
 405         movel   sp@(12), d1     /* d1 = divisor */
 406         jpl     L1
 407         negl    d1
 408 #ifndef __mcoldfire__
 409         negb    d2              /* change sign because divisor <0  */
 410 #else
 411         negl    d2              /* change sign because divisor <0  */
 412 #endif
 413 L1:     movel   sp@(8), d0      /* d0 = dividend */
 414         jpl     L2
 415         negl    d0
 416 #ifndef __mcoldfire__
 417         negb    d2
 418 #else
 419         negl    d2
 420 #endif
 421
 422 L2:     movel   d1, sp@-
 423         movel   d0, sp@-
 424         jbsr    SYM (__udivsi3) /* divide abs(dividend) by abs(divisor) */
 425         addql   IMM (8), sp
 426
 427         tstb    d2
 428         jpl     L3
 429         negl    d0
 430
 431 L3:     movel   sp@+, d2
 432         rts
 433 #endif /* L_divsi3 */
 434
 435 #ifdef  L_umodsi3
 436         .text
 437         .proc
 438         .globl  SYM (__umodsi3)
 439 SYM (__umodsi3):
 440         movel   sp@(8), d1      /* d1 = divisor */
 441         movel   sp@(4), d0      /* d0 = dividend */
 442         movel   d1, sp@-
 443         movel   d0, sp@-
 444         jbsr    SYM (__udivsi3)
 445         addql   IMM (8), sp
 446         movel   sp@(8), d1      /* d1 = divisor */
 447 #ifndef __mcoldfire__
 448         movel   d1, sp@-
 449         movel   d0, sp@-
 450         jbsr    SYM (__mulsi3)  /* d0 = (a/b)*b */
 451         addql   IMM (8), sp
 452 #else
 453         mulsl   d1,d0
 454 #endif
 455         movel   sp@(4), d1      /* d1 = dividend */
 456         subl    d0, d1          /* d1 = a - (a/b)*b */
 457         movel   d1, d0
 458         rts
 459 #endif /* L_umodsi3 */
 460
 461 #ifdef  L_modsi3
 462         .text
 463         .proc
 464         .globl  SYM (__modsi3)
 465 SYM (__modsi3):
 466         movel   sp@(8), d1      /* d1 = divisor */
 467         movel   sp@(4), d0      /* d0 = dividend */
 468         movel   d1, sp@-
 469         movel   d0, sp@-
 470         jbsr    SYM (__divsi3)
 471         addql   IMM (8), sp
 472         movel   sp@(8), d1      /* d1 = divisor */
 473 #ifndef __mcoldfire__
 474         movel   d1, sp@-
 475         movel   d0, sp@-
 476         jbsr    SYM (__mulsi3)  /* d0 = (a/b)*b */
 477         addql   IMM (8), sp
 478 #else
 479         mulsl   d1,d0
 480 #endif
 481         movel   sp@(4), d1      /* d1 = dividend */
 482         subl    d0, d1          /* d1 = a - (a/b)*b */
 483         movel   d1, d0
 484         rts
 485 #endif /* L_modsi3 */
 486
 487
 488 #ifdef  L_double
 489
 490         .globl  SYM (_fpCCR)
 491         .globl  $_exception_handler
 492
 493 QUIET_NaN      = 0xffffffff
 494
 495 D_MAX_EXP      = 0x07ff
 496 D_BIAS         = 1022
 497 DBL_MAX_EXP    = D_MAX_EXP - D_BIAS
 498 DBL_MIN_EXP    = 1 - D_BIAS
 499 DBL_MANT_DIG   = 53
 500
 501 INEXACT_RESULT          = 0x0001
 502 UNDERFLOW               = 0x0002
 503 OVERFLOW                = 0x0004
 504 DIVIDE_BY_ZERO          = 0x0008
 505 INVALID_OPERATION       = 0x0010
 506
 507 DOUBLE_FLOAT = 2
 508
 509 NOOP         = 0
 510 ADD          = 1
 511 MULTIPLY     = 2
 512 DIVIDE       = 3
 513 NEGATE       = 4
 514 COMPARE      = 5
 515 EXTENDSFDF   = 6
 516 TRUNCDFSF    = 7
 517
 518 UNKNOWN           = -1
 519 ROUND_TO_NEAREST  = 0 | round result to nearest representable value
 520 ROUND_TO_ZERO     = 1 | round result towards zero
 521 ROUND_TO_PLUS     = 2 | round result towards plus infinity
 522 ROUND_TO_MINUS    = 3 | round result towards minus infinity
 523
 524 | Entry points:
 525
 526         .globl SYM (__adddf3)
 527         .globl SYM (__subdf3)
 528         .globl SYM (__muldf3)
 529         .globl SYM (__divdf3)
 530         .globl SYM (__negdf2)
 531         .globl SYM (__cmpdf2)
 532
 533         .text
 534         .even
 535
 536 | These are common routines to return and signal exceptions.
 537
 538 Ld$den:
 539 | Return and signal a denormalized number
 540         orl     d7,d0
 541         movew   IMM (INEXACT_RESULT+UNDERFLOW),d7
 542         moveq   IMM (DOUBLE_FLOAT),d6
 543         jmp     $_exception_handler
 544
 545 Ld$infty:
 546 Ld$overflow:
 547 | Return a properly signed INFINITY and set the exception flags
 548         movel   IMM (0x7ff00000),d0
 549         movel   IMM (0),d1
 550         orl     d7,d0
 551         movew   IMM (INEXACT_RESULT+OVERFLOW),d7
 552         moveq   IMM (DOUBLE_FLOAT),d6
 553         jmp     $_exception_handler
 554
 555 Ld$underflow:
 556 | Return 0 and set the exception flags
 557         movel   IMM (0),d0
 558         movel   d0,d1
 559         movew   IMM (INEXACT_RESULT+UNDERFLOW),d7
 560         moveq   IMM (DOUBLE_FLOAT),d6
 561         jmp     $_exception_handler
 562
 563 Ld$inop:
 564 | Return a quiet NaN and set the exception flags
 565         movel   IMM (QUIET_NaN),d0
 566         movel   d0,d1
 567         movew   IMM (INEXACT_RESULT+INVALID_OPERATION),d7
 568         moveq   IMM (DOUBLE_FLOAT),d6
 569         jmp     $_exception_handler
 570
 571 Ld$div$0:
 572 | Return a properly signed INFINITY and set the exception flags
 573         movel   IMM (0x7ff00000),d0
 574         movel   IMM (0),d1
 575         orl     d7,d0
 576         movew   IMM (INEXACT_RESULT+DIVIDE_BY_ZERO),d7
 577         moveq   IMM (DOUBLE_FLOAT),d6
 578         jmp     $_exception_handler
 579
 580 |=============================================================================
 581 |=============================================================================
 582 |                         double precision routines
 583 |=============================================================================
 584 |=============================================================================
 585
 586 | A double precision floating point number (double) has the format:
 587 |
 588 | struct _double {
 589 |  unsigned int sign      : 1;  /* sign bit */
 590 |  unsigned int exponent  : 11; /* exponent, shifted by 126 */
 591 |  unsigned int fraction  : 52; /* fraction */
 592 | } double;
 593 |
 594 | Thus sizeof(double) = 8 (64 bits).
 595 |
 596 | All the routines are callable from C programs, and return the result
 597 | in the register pair d0-d1. They also preserve all registers except
 598 | d0-d1 and a0-a1.
 599
 600 |=============================================================================
 601 |                              __subdf3
 602 |=============================================================================
 603
 604 | double __subdf3(double, double);
 605 SYM (__subdf3):
 606         bchg    IMM (31),sp@(12) | change sign of second operand
 607                                 | and fall through, so we always add
 608 |=============================================================================
 609 |                              __adddf3
 610 |=============================================================================
 611
 612 | double __adddf3(double, double);
 613 SYM (__adddf3):
 614 #ifndef __mcoldfire__
 615         link    a6,IMM (0)      | everything will be done in registers
 616         moveml  d2-d7,sp@-      | save all data registers and a2 (but d0-d1)
 617 #else
 618         link    a6,IMM (-24)
 619         moveml  d2-d7,sp@
 620 #endif
 621         movel   a6@(8),d0       | get first operand
 622         movel   a6@(12),d1      |
 623         movel   a6@(16),d2      | get second operand
 624         movel   a6@(20),d3      |
 625
 626         movel   d0,d7           | get d0's sign bit in d7 '
 627         addl    d1,d1           | check and clear sign bit of a, and gain one
 628         addxl   d0,d0           | bit of extra precision
 629         beq     Ladddf$b        | if zero return second operand
 630
 631         movel   d2,d6           | save sign in d6
 632         addl    d3,d3           | get rid of sign bit and gain one bit of
 633         addxl   d2,d2           | extra precision
 634         beq     Ladddf$a        | if zero return first operand
 635
 636         andl    IMM (0x80000000),d7 | isolate a's sign bit '
 637         swap    d6              | and also b's sign bit '
 638 #ifndef __mcoldfire__
 639         andw    IMM (0x8000),d6 |
 640         orw     d6,d7           | and combine them into d7, so that a's sign '
 641                                 | bit is in the high word and b's is in the '
 642                                 | low word, so d6 is free to be used
 643 #else
 644         andl    IMM (0x8000),d6
 645         orl     d6,d7
 646 #endif
 647         movel   d7,a0           | now save d7 into a0, so d7 is free to
 648                                 | be used also
 649
 650 | Get the exponents and check for denormalized and/or infinity.
 651
 652         movel   IMM (0x001fffff),d6 | mask for the fraction
 653         movel   IMM (0x00200000),d7 | mask to put hidden bit back
 654
 655         movel   d0,d4           |
 656         andl    d6,d0           | get fraction in d0
 657         notl    d6              | make d6 into mask for the exponent
 658         andl    d6,d4           | get exponent in d4
 659         beq     Ladddf$a$den    | branch if a is denormalized
 660         cmpl    d6,d4           | check for INFINITY or NaN
 661         beq     Ladddf$nf       |
 662         orl     d7,d0           | and put hidden bit back
 663 Ladddf$1:
 664         swap    d4              | shift right exponent so that it starts
 665 #ifndef __mcoldfire__
 666         lsrw    IMM (5),d4      | in bit 0 and not bit 20
 667 #else
 668         lsrl    IMM (5),d4      | in bit 0 and not bit 20
 669 #endif
 670 | Now we have a's exponent in d4 and fraction in d0-d1 '
 671         movel   d2,d5           | save b to get exponent
 672         andl    d6,d5           | get exponent in d5
 673         beq     Ladddf$b$den    | branch if b is denormalized
 674         cmpl    d6,d5           | check for INFINITY or NaN
 675         beq     Ladddf$nf
 676         notl    d6              | make d6 into mask for the fraction again
 677         andl    d6,d2           | and get fraction in d2
 678         orl     d7,d2           | and put hidden bit back
 679 Ladddf$2:
 680         swap    d5              | shift right exponent so that it starts
 681 #ifndef __mcoldfire__
 682         lsrw    IMM (5),d5      | in bit 0 and not bit 20
 683 #else
 684         lsrl    IMM (5),d5      | in bit 0 and not bit 20
 685 #endif
 686
 687 | Now we have b's exponent in d5 and fraction in d2-d3. '
 688
 689 | The situation now is as follows: the signs are combined in a0, the
 690 | numbers are in d0-d1 (a) and d2-d3 (b), and the exponents in d4 (a)
 691 | and d5 (b). To do the rounding correctly we need to keep all the
 692 | bits until the end, so we need to use d0-d1-d2-d3 for the first number
 693 | and d4-d5-d6-d7 for the second. To do this we store (temporarily) the
 694 | exponents in a2-a3.
 695
 696 #ifndef __mcoldfire__
 697         moveml  a2-a3,sp@-      | save the address registers
 698 #else
 699         movel   a2,sp@-
 700         movel   a3,sp@-
 701         movel   a4,sp@-
 702 #endif
 703
 704         movel   d4,a2           | save the exponents
 705         movel   d5,a3           |
 706
 707         movel   IMM (0),d7      | and move the numbers around
 708         movel   d7,d6           |
 709         movel   d3,d5           |
 710         movel   d2,d4           |
 711         movel   d7,d3           |
 712         movel   d7,d2           |
 713
 714 | Here we shift the numbers until the exponents are the same, and put
 715 | the largest exponent in a2.
 716 #ifndef __mcoldfire__
 717         exg     d4,a2           | get exponents back
 718         exg     d5,a3           |
 719         cmpw    d4,d5           | compare the exponents
 720 #else
 721         movel   d4,a4           | get exponents back
 722         movel   a2,d4
 723         movel   a4,a2
 724         movel   d5,a4
 725         movel   a3,d5
 726         movel   a4,a3
 727         cmpl    d4,d5           | compare the exponents
 728 #endif
 729         beq     Ladddf$3        | if equal don't shift '
 730         bhi     9f              | branch if second exponent is higher
 731
 732 | Here we have a's exponent larger than b's, so we have to shift b. We do
 733 | this by using as counter d2:
 734 1:      movew   d4,d2           | move largest exponent to d2
 735 #ifndef __mcoldfire__
 736         subw    d5,d2           | and subtract second exponent
 737         exg     d4,a2           | get back the longs we saved
 738         exg     d5,a3           |
 739 #else
 740         subl    d5,d2           | and subtract second exponent
 741         movel   d4,a4           | get back the longs we saved
 742         movel   a2,d4
 743         movel   a4,a2
 744         movel   d5,a4
 745         movel   a3,d5
 746         movel   a4,a3
 747 #endif
 748 | if difference is too large we don't shift (actually, we can just exit) '
 749 #ifndef __mcoldfire__
 750         cmpw    IMM (DBL_MANT_DIG+2),d2
 751 #else
 752         cmpl    IMM (DBL_MANT_DIG+2),d2
 753 #endif
 754         bge     Ladddf$b$small
 755 #ifndef __mcoldfire__
 756         cmpw    IMM (32),d2     | if difference >= 32, shift by longs
 757 #else
 758         cmpl    IMM (32),d2     | if difference >= 32, shift by longs
 759 #endif
 760         bge     5f
 761 2:
 762 #ifndef __mcoldfire__
 763         cmpw    IMM (16),d2     | if difference >= 16, shift by words
 764 #else
 765         cmpl    IMM (16),d2     | if difference >= 16, shift by words
 766 #endif
 767         bge     6f
 768         bra     3f              | enter dbra loop
 769
 770 4:
 771 #ifndef __mcoldfire__
 772         lsrl    IMM (1),d4
 773         roxrl   IMM (1),d5
 774         roxrl   IMM (1),d6
 775         roxrl   IMM (1),d7
 776 #else
 777         lsrl    IMM (1),d7
 778         btst    IMM (0),d6
 779         beq     10f
 780         bset    IMM (31),d7
 781 10:     lsrl    IMM (1),d6
 782         btst    IMM (0),d5
 783         beq     11f
 784         bset    IMM (31),d6
 785 11:     lsrl    IMM (1),d5
 786         btst    IMM (0),d4
 787         beq     12f
 788         bset    IMM (31),d5
 789 12:     lsrl    IMM (1),d4
 790 #endif
 791 3:
 792 #ifndef __mcoldfire__
 793         dbra    d2,4b
 794 #else
 795         subql   IMM (1),d2
 796         bpl     4b
 797 #endif
 798         movel   IMM (0),d2
 799         movel   d2,d3
 800         bra     Ladddf$4
 801 5:
 802         movel   d6,d7
 803         movel   d5,d6
 804         movel   d4,d5
 805         movel   IMM (0),d4
 806 #ifndef __mcoldfire__
 807         subw    IMM (32),d2
 808 #else
 809         subl    IMM (32),d2
 810 #endif
 811         bra     2b
 812 6:
 813         movew   d6,d7
 814         swap    d7
 815         movew   d5,d6
 816         swap    d6
 817         movew   d4,d5
 818         swap    d5
 819         movew   IMM (0),d4
 820         swap    d4
 821 #ifndef __mcoldfire__
 822         subw    IMM (16),d2
 823 #else
 824         subl    IMM (16),d2
 825 #endif
 826         bra     3b
 827
 828 9:
 829 #ifndef __mcoldfire__
 830         exg     d4,d5
 831         movew   d4,d6
 832         subw    d5,d6           | keep d5 (largest exponent) in d4
 833         exg     d4,a2
 834         exg     d5,a3
 835 #else
 836         movel   d5,d6
 837         movel   d4,d5
 838         movel   d6,d4
 839         subl    d5,d6
 840         movel   d4,a4
 841         movel   a2,d4
 842         movel   a4,a2
 843         movel   d5,a4
 844         movel   a3,d5
 845         movel   a4,a3
 846 #endif
 847 | if difference is too large we don't shift (actually, we can just exit) '
 848 #ifndef __mcoldfire__
 849         cmpw    IMM (DBL_MANT_DIG+2),d6
 850 #else
 851         cmpl    IMM (DBL_MANT_DIG+2),d6
 852 #endif
 853         bge     Ladddf$a$small
 854 #ifndef __mcoldfire__
 855         cmpw    IMM (32),d6     | if difference >= 32, shift by longs
 856 #else
 857         cmpl    IMM (32),d6     | if difference >= 32, shift by longs
 858 #endif
 859         bge     5f
 860 2:
 861 #ifndef __mcoldfire__
 862         cmpw    IMM (16),d6     | if difference >= 16, shift by words
 863 #else
 864         cmpl    IMM (16),d6     | if difference >= 16, shift by words
 865 #endif
 866         bge     6f
 867         bra     3f              | enter dbra loop
 868
 869 4:
 870 #ifndef __mcoldfire__
 871         lsrl    IMM (1),d0
 872         roxrl   IMM (1),d1
 873         roxrl   IMM (1),d2
 874         roxrl   IMM (1),d3
 875 #else
 876         lsrl    IMM (1),d3
 877         btst    IMM (0),d2
 878         beq     10f
 879         bset    IMM (31),d3
 880 10:     lsrl    IMM (1),d2
 881         btst    IMM (0),d1
 882         beq     11f
 883         bset    IMM (31),d2
 884 11:     lsrl    IMM (1),d1
 885         btst    IMM (0),d0
 886         beq     12f
 887         bset    IMM (31),d1
 888 12:     lsrl    IMM (1),d0
 889 #endif
 890 3:
 891 #ifndef __mcoldfire__
 892         dbra    d6,4b
 893 #else
 894         subql   IMM (1),d6
 895         bpl     4b
 896 #endif
 897         movel   IMM (0),d7
 898         movel   d7,d6
 899         bra     Ladddf$4
 900 5:
 901         movel   d2,d3
 902         movel   d1,d2
 903         movel   d0,d1
 904         movel   IMM (0),d0
 905 #ifndef __mcoldfire__
 906         subw    IMM (32),d6
 907 #else
 908         subl    IMM (32),d6
 909 #endif
 910         bra     2b
 911 6:
 912         movew   d2,d3
 913         swap    d3
 914         movew   d1,d2
 915         swap    d2
 916         movew   d0,d1
 917         swap    d1
 918         movew   IMM (0),d0
 919         swap    d0
 920 #ifndef __mcoldfire__
 921         subw    IMM (16),d6
 922 #else
 923         subl    IMM (16),d6
 924 #endif
 925         bra     3b
 926 Ladddf$3:
 927 #ifndef __mcoldfire__
 928         exg     d4,a2
 929         exg     d5,a3
 930 #else
 931         movel   d4,a4
 932         movel   a2,d4
 933         movel   a4,a2
 934         movel   d5,a4
 935         movel   a3,d5
 936         movel   a4,a3
 937 #endif
 938 Ladddf$4:
 939 | Now we have the numbers in d0--d3 and d4--d7, the exponent in a2, and
 940 | the signs in a4.
 941
 942 | Here we have to decide whether to add or subtract the numbers:
 943 #ifndef __mcoldfire__
 944         exg     d7,a0           | get the signs
 945         exg     d6,a3           | a3 is free to be used
 946 #else
 947         movel   d7,a4
 948         movel   a0,d7
 949         movel   a4,a0
 950         movel   d6,a4
 951         movel   a3,d6
 952         movel   a4,a3
 953 #endif
 954         movel   d7,d6           |
 955         movew   IMM (0),d7      | get a's sign in d7 '
 956         swap    d6              |
 957         movew   IMM (0),d6      | and b's sign in d6 '
 958         eorl    d7,d6           | compare the signs
 959         bmi     Lsubdf$0        | if the signs are different we have
 960                                 | to subtract
 961 #ifndef __mcoldfire__
 962         exg     d7,a0           | else we add the numbers
 963         exg     d6,a3           |
 964 #else
 965         movel   d7,a4
 966         movel   a0,d7
 967         movel   a4,a0
 968         movel   d6,a4
 969         movel   a3,d6
 970         movel   a4,a3
 971 #endif
 972         addl    d7,d3           |
 973         addxl   d6,d2           |
 974         addxl   d5,d1           |
 975         addxl   d4,d0           |
 976
 977         movel   a2,d4           | return exponent to d4
 978         movel   a0,d7           |
 979         andl    IMM (0x80000000),d7 | d7 now has the sign
 980
 981 #ifndef __mcoldfire__
 982         moveml  sp@+,a2-a3
 983 #else
 984         movel   sp@+,a4
 985         movel   sp@+,a3
 986         movel   sp@+,a2
 987 #endif
 988
 989 | Before rounding normalize so bit #DBL_MANT_DIG is set (we will consider
 990 | the case of denormalized numbers in the rounding routine itself).
 991 | As in the addition (not in the subtraction!) we could have set
 992 | one more bit we check this:
 993         btst    IMM (DBL_MANT_DIG+1),d0
 994         beq     1f
 995 #ifndef __mcoldfire__
 996         lsrl    IMM (1),d0
 997         roxrl   IMM (1),d1
 998         roxrl   IMM (1),d2
 999         roxrl   IMM (1),d3
1000         addw    IMM (1),d4
1001 #else
1002         lsrl    IMM (1),d3
1003         btst    IMM (0),d2
1004         beq     10f
1005         bset    IMM (31),d3
1006 10:     lsrl    IMM (1),d2
1007         btst    IMM (0),d1
1008         beq     11f
1009         bset    IMM (31),d2
1010 11:     lsrl    IMM (1),d1
1011         btst    IMM (0),d0
1012         beq     12f
1013         bset    IMM (31),d1
1014 12:     lsrl    IMM (1),d0
1015         addl    IMM (1),d4
1016 #endif
1017 1:
1018         lea     Ladddf$5,a0     | to return from rounding routine
1019         lea     SYM (_fpCCR),a1 | check the rounding mode
1020 #ifdef __mcoldfire__
1021         clrl    d6
1022 #endif
1023         movew   a1@(6),d6       | rounding mode in d6
1024         beq     Lround$to$nearest
1025 #ifndef __mcoldfire__
1026         cmpw    IMM (ROUND_TO_PLUS),d6
1027 #else
1028         cmpl    IMM (ROUND_TO_PLUS),d6
1029 #endif
1030         bhi     Lround$to$minus
1031         blt     Lround$to$zero
1032         bra     Lround$to$plus
1033 Ladddf$5:
1034 | Put back the exponent and check for overflow
1035 #ifndef __mcoldfire__
1036         cmpw    IMM (0x7ff),d4  | is the exponent big?
1037 #else
1038         cmpl    IMM (0x7ff),d4  | is the exponent big?
1039 #endif
1040         bge     1f
1041         bclr    IMM (DBL_MANT_DIG-1),d0
1042 #ifndef __mcoldfire__
1043         lslw    IMM (4),d4      | put exponent back into position
1044 #else
1045         lsll    IMM (4),d4      | put exponent back into position
1046 #endif
1047         swap    d0              |
1048 #ifndef __mcoldfire__
1049         orw     d4,d0           |
1050 #else
1051         orl     d4,d0           |
1052 #endif
1053         swap    d0              |
1054         bra     Ladddf$ret
1055 1:
1056         movew   IMM (ADD),d5
1057         bra     Ld$overflow
1058
1059 Lsubdf$0:
1060 | Here we do the subtraction.
1061 #ifndef __mcoldfire__
1062         exg     d7,a0           | put sign back in a0
1063         exg     d6,a3           |
1064 #else
1065         movel   d7,a4
1066         movel   a0,d7
1067         movel   a4,a0
1068         movel   d6,a4
1069         movel   a3,d6
1070         movel   a4,a3
1071 #endif
1072         subl    d7,d3           |
1073         subxl   d6,d2           |
1074         subxl   d5,d1           |
1075         subxl   d4,d0           |
1076         beq     Ladddf$ret$1    | if zero just exit
1077         bpl     1f              | if positive skip the following
1078         movel   a0,d7           |
1079         bchg    IMM (31),d7     | change sign bit in d7
1080         movel   d7,a0           |
1081         negl    d3              |
1082         negxl   d2              |
1083         negxl   d1              | and negate result
1084         negxl   d0              |
1085 1:
1086         movel   a2,d4           | return exponent to d4
1087         movel   a0,d7
1088         andl    IMM (0x80000000),d7 | isolate sign bit
1089 #ifndef __mcoldfire__
1090         moveml  sp@+,a2-a3      |
1091 #else
1092         movel   sp@+,a4
1093         movel   sp@+,a3
1094         movel   sp@+,a2
1095 #endif
1096
1097 | Before rounding normalize so bit #DBL_MANT_DIG is set (we will consider
1098 | the case of denormalized numbers in the rounding routine itself).
1099 | As in the addition (not in the subtraction!) we could have set
1100 | one more bit we check this:
1101         btst    IMM (DBL_MANT_DIG+1),d0
1102         beq     1f
1103 #ifndef __mcoldfire__
1104         lsrl    IMM (1),d0
1105         roxrl   IMM (1),d1
1106         roxrl   IMM (1),d2
1107         roxrl   IMM (1),d3
1108         addw    IMM (1),d4
1109 #else
1110         lsrl    IMM (1),d3
1111         btst    IMM (0),d2
1112         beq     10f
1113         bset    IMM (31),d3
1114 10:     lsrl    IMM (1),d2
1115         btst    IMM (0),d1
1116         beq     11f
1117         bset    IMM (31),d2
1118 11:     lsrl    IMM (1),d1
1119         btst    IMM (0),d0
1120         beq     12f
1121         bset    IMM (31),d1
1122 12:     lsrl    IMM (1),d0
1123         addl    IMM (1),d4
1124 #endif
1125 1:
1126         lea     Lsubdf$1,a0     | to return from rounding routine
1127         lea     SYM (_fpCCR),a1 | check the rounding mode
1128 #ifdef __mcoldfire__
1129         clrl    d6
1130 #endif
1131         movew   a1@(6),d6       | rounding mode in d6
1132         beq     Lround$to$nearest
1133 #ifndef __mcoldfire__
1134         cmpw    IMM (ROUND_TO_PLUS),d6
1135 #else
1136         cmpl    IMM (ROUND_TO_PLUS),d6
1137 #endif
1138         bhi     Lround$to$minus
1139         blt     Lround$to$zero
1140         bra     Lround$to$plus
1141 Lsubdf$1:
1142 | Put back the exponent and sign (we don't have overflow). '
1143         bclr    IMM (DBL_MANT_DIG-1),d0
1144 #ifndef __mcoldfire__
1145         lslw    IMM (4),d4      | put exponent back into position
1146 #else
1147         lsll    IMM (4),d4      | put exponent back into position
1148 #endif
1149         swap    d0              |
1150 #ifndef __mcoldfire__
1151         orw     d4,d0           |
1152 #else
1153         orl     d4,d0           |
1154 #endif
1155         swap    d0              |
1156         bra     Ladddf$ret
1157
1158 | If one of the numbers was too small (difference of exponents >=
1159 | DBL_MANT_DIG+1) we return the other (and now we don't have to '
1160 | check for finiteness or zero).
1161 Ladddf$a$small:
1162 #ifndef __mcoldfire__
1163         moveml  sp@+,a2-a3
1164 #else
1165         movel   sp@+,a4
1166         movel   sp@+,a3
1167         movel   sp@+,a2
1168 #endif
1169         movel   a6@(16),d0
1170         movel   a6@(20),d1
1171         lea     SYM (_fpCCR),a0
1172         movew   IMM (0),a0@
1173 #ifndef __mcoldfire__
1174         moveml  sp@+,d2-d7      | restore data registers
1175 #else
1176         moveml  sp@,d2-d7
1177         | XXX if frame pointer is ever removed, stack pointer must
1178         | be adjusted here.
1179 #endif
1180         unlk    a6              | and return
1181         rts
1182
1183 Ladddf$b$small:
1184 #ifndef __mcoldfire__
1185         moveml  sp@+,a2-a3
1186 #else
1187         movel   sp@+,a4
1188         movel   sp@+,a3
1189         movel   sp@+,a2
1190 #endif
1191         movel   a6@(8),d0
1192         movel   a6@(12),d1
1193         lea     SYM (_fpCCR),a0
1194         movew   IMM (0),a0@
1195 #ifndef __mcoldfire__
1196         moveml  sp@+,d2-d7      | restore data registers
1197 #else
1198         moveml  sp@,d2-d7
1199         | XXX if frame pointer is ever removed, stack pointer must
1200         | be adjusted here.
1201 #endif
1202         unlk    a6              | and return
1203         rts
1204
1205 Ladddf$a$den:
1206         movel   d7,d4           | d7 contains 0x00200000
1207         bra     Ladddf$1
1208
1209 Ladddf$b$den:
1210         movel   d7,d5           | d7 contains 0x00200000
1211         notl    d6
1212         bra     Ladddf$2
1213
1214 Ladddf$b:
1215 | Return b (if a is zero)
1216         movel   d2,d0
1217         movel   d3,d1
1218         bra     1f
1219 Ladddf$a:
1220         movel   a6@(8),d0
1221         movel   a6@(12),d1
1222 1:
1223         movew   IMM (ADD),d5
1224 | Check for NaN and +/-INFINITY.
1225         movel   d0,d7                   |
1226         andl    IMM (0x80000000),d7     |
1227         bclr    IMM (31),d0             |
1228         cmpl    IMM (0x7ff00000),d0     |
1229         bge     2f                      |
1230         movel   d0,d0                   | check for zero, since we don't  '
1231         bne     Ladddf$ret              | want to return -0 by mistake
1232         bclr    IMM (31),d7             |
1233         bra     Ladddf$ret              |
1234 2:
1235         andl    IMM (0x000fffff),d0     | check for NaN (nonzero fraction)
1236         orl     d1,d0                   |
1237         bne     Ld$inop                 |
1238         bra     Ld$infty                |
1239
1240 Ladddf$ret$1:
1241 #ifndef __mcoldfire__
1242         moveml  sp@+,a2-a3      | restore regs and exit
1243 #else
1244         movel   sp@+,a4
1245         movel   sp@+,a3
1246         movel   sp@+,a2
1247 #endif
1248
1249 Ladddf$ret:
1250 | Normal exit.
1251         lea     SYM (_fpCCR),a0
1252         movew   IMM (0),a0@
1253         orl     d7,d0           | put sign bit back
1254 #ifndef __mcoldfire__
1255         moveml  sp@+,d2-d7
1256 #else
1257         moveml  sp@,d2-d7
1258         | XXX if frame pointer is ever removed, stack pointer must
1259         | be adjusted here.
1260 #endif
1261         unlk    a6
1262         rts
1263
1264 Ladddf$ret$den:
1265 | Return a denormalized number.
1266 #ifndef __mcoldfire__
1267         lsrl    IMM (1),d0      | shift right once more
1268         roxrl   IMM (1),d1      |
1269 #else
1270         lsrl    IMM (1),d1
1271         btst    IMM (0),d0
1272         beq     10f
1273         bset    IMM (31),d1
1274 10:     lsrl    IMM (1),d0
1275 #endif
1276         bra     Ladddf$ret
1277
1278 Ladddf$nf:
1279         movew   IMM (ADD),d5
1280 | This could be faster but it is not worth the effort, since it is not
1281 | executed very often. We sacrifice speed for clarity here.
1282         movel   a6@(8),d0       | get the numbers back (remember that we
1283         movel   a6@(12),d1      | did some processing already)
1284         movel   a6@(16),d2      |
1285         movel   a6@(20),d3      |
1286         movel   IMM (0x7ff00000),d4 | useful constant (INFINITY)
1287         movel   d0,d7           | save sign bits
1288         movel   d2,d6           |
1289         bclr    IMM (31),d0     | clear sign bits
1290         bclr    IMM (31),d2     |
1291 | We know that one of them is either NaN of +/-INFINITY
1292 | Check for NaN (if either one is NaN return NaN)
1293         cmpl    d4,d0           | check first a (d0)
1294         bhi     Ld$inop         | if d0 > 0x7ff00000 or equal and
1295         bne     2f
1296         tstl    d1              | d1 > 0, a is NaN
1297         bne     Ld$inop         |
1298 2:      cmpl    d4,d2           | check now b (d1)
1299         bhi     Ld$inop         |
1300         bne     3f
1301         tstl    d3              |
1302         bne     Ld$inop         |
1303 3:
1304 | Now comes the check for +/-INFINITY. We know that both are (maybe not
1305 | finite) numbers, but we have to check if both are infinite whether we
1306 | are adding or subtracting them.
1307         eorl    d7,d6           | to check sign bits
1308         bmi     1f
1309         andl    IMM (0x80000000),d7 | get (common) sign bit
1310         bra     Ld$infty
1311 1:
1312 | We know one (or both) are infinite, so we test for equality between the
1313 | two numbers (if they are equal they have to be infinite both, so we
1314 | return NaN).
1315         cmpl    d2,d0           | are both infinite?
1316         bne     1f              | if d0 <> d2 they are not equal
1317         cmpl    d3,d1           | if d0 == d2 test d3 and d1
1318         beq     Ld$inop         | if equal return NaN
1319 1:
1320         andl    IMM (0x80000000),d7 | get a's sign bit '
1321         cmpl    d4,d0           | test now for infinity
1322         beq     Ld$infty        | if a is INFINITY return with this sign
1323         bchg    IMM (31),d7     | else we know b is INFINITY and has
1324         bra     Ld$infty        | the opposite sign
1325
1326 |=============================================================================
1327 |                              __muldf3
1328 |=============================================================================
1329
1330 | double __muldf3(double, double);
1331 SYM (__muldf3):
1332 #ifndef __mcoldfire__
1333         link    a6,IMM (0)
1334         moveml  d2-d7,sp@-
1335 #else
1336         link    a6,IMM (-24)
1337         moveml  d2-d7,sp@
1338 #endif
1339         movel   a6@(8),d0               | get a into d0-d1
1340         movel   a6@(12),d1              |
1341         movel   a6@(16),d2              | and b into d2-d3
1342         movel   a6@(20),d3              |
1343         movel   d0,d7                   | d7 will hold the sign of the product
1344         eorl    d2,d7                   |
1345         andl    IMM (0x80000000),d7     |
1346         movel   d7,a0                   | save sign bit into a0
1347         movel   IMM (0x7ff00000),d7     | useful constant (+INFINITY)
1348         movel   d7,d6                   | another (mask for fraction)
1349         notl    d6                      |
1350         bclr    IMM (31),d0             | get rid of a's sign bit '
1351         movel   d0,d4                   |
1352         orl     d1,d4                   |
1353         beq     Lmuldf$a$0              | branch if a is zero
1354         movel   d0,d4                   |
1355         bclr    IMM (31),d2             | get rid of b's sign bit '
1356         movel   d2,d5                   |
1357         orl     d3,d5                   |
1358         beq     Lmuldf$b$0              | branch if b is zero
1359         movel   d2,d5                   |
1360         cmpl    d7,d0                   | is a big?
1361         bhi     Lmuldf$inop             | if a is NaN return NaN
1362         beq     Lmuldf$a$nf             | we still have to check d1 and b ...
1363         cmpl    d7,d2                   | now compare b with INFINITY
1364         bhi     Lmuldf$inop             | is b NaN?
1365         beq     Lmuldf$b$nf             | we still have to check d3 ...
1366 | Here we have both numbers finite and nonzero (and with no sign bit).
1367 | Now we get the exponents into d4 and d5.
1368         andl    d7,d4                   | isolate exponent in d4
1369         beq     Lmuldf$a$den            | if exponent zero, have denormalized
1370         andl    d6,d0                   | isolate fraction
1371         orl     IMM (0x00100000),d0     | and put hidden bit back
1372         swap    d4                      | I like exponents in the first byte
1373 #ifndef __mcoldfire__
1374         lsrw    IMM (4),d4              |
1375 #else
1376         lsrl    IMM (4),d4              |
1377 #endif
1378 Lmuldf$1:
1379         andl    d7,d5                   |
1380         beq     Lmuldf$b$den            |
1381         andl    d6,d2                   |
1382         orl     IMM (0x00100000),d2     | and put hidden bit back
1383         swap    d5                      |
1384 #ifndef __mcoldfire__
1385         lsrw    IMM (4),d5              |
1386 #else
1387         lsrl    IMM (4),d5              |
1388 #endif
1389 Lmuldf$2:                               |
1390 #ifndef __mcoldfire__
1391         addw    d5,d4                   | add exponents
1392         subw    IMM (D_BIAS+1),d4       | and subtract bias (plus one)
1393 #else
1394         addl    d5,d4                   | add exponents
1395         subl    IMM (D_BIAS+1),d4       | and subtract bias (plus one)
1396 #endif
1397
1398 | We are now ready to do the multiplication. The situation is as follows:
1399 | both a and b have bit 52 ( bit 20 of d0 and d2) set (even if they were
1400 | denormalized to start with!), which means that in the product bit 104
1401 | (which will correspond to bit 8 of the fourth long) is set.
1402
1403 | Here we have to do the product.
1404 | To do it we have to juggle the registers back and forth, as there are not
1405 | enough to keep everything in them. So we use the address registers to keep
1406 | some intermediate data.
1407
1408 #ifndef __mcoldfire__
1409         moveml  a2-a3,sp@-      | save a2 and a3 for temporary use
1410 #else
1411         movel   a2,sp@-
1412         movel   a3,sp@-
1413         movel   a4,sp@-
1414 #endif
1415         movel   IMM (0),a2      | a2 is a null register
1416         movel   d4,a3           | and a3 will preserve the exponent
1417
1418 | First, shift d2-d3 so bit 20 becomes bit 31:
1419 #ifndef __mcoldfire__
1420         rorl    IMM (5),d2      | rotate d2 5 places right
1421         swap    d2              | and swap it
1422         rorl    IMM (5),d3      | do the same thing with d3
1423         swap    d3              |
1424         movew   d3,d6           | get the rightmost 11 bits of d3
1425         andw    IMM (0x07ff),d6 |
1426         orw     d6,d2           | and put them into d2
1427         andw    IMM (0xf800),d3 | clear those bits in d3
1428 #else
1429         moveq   IMM (11),d7     | left shift d2 11 bits
1430         lsll    d7,d2
1431         movel   d3,d6           | get a copy of d3
1432         lsll    d7,d3           | left shift d3 11 bits
1433         andl    IMM (0xffe00000),d6 | get the top 11 bits of d3
1434         moveq   IMM (21),d7     | right shift them 21 bits
1435         lsrl    d7,d6
1436         orl     d6,d2           | stick them at the end of d2
1437 #endif
1438
1439         movel   d2,d6           | move b into d6-d7
1440         movel   d3,d7           | move a into d4-d5
1441         movel   d0,d4           | and clear d0-d1-d2-d3 (to put result)
1442         movel   d1,d5           |
1443         movel   IMM (0),d3      |
1444         movel   d3,d2           |
1445         movel   d3,d1           |
1446         movel   d3,d0           |
1447
1448 | We use a1 as counter:
1449         movel   IMM (DBL_MANT_DIG-1),a1
1450 #ifndef __mcoldfire__
1451         exg     d7,a1
1452 #else
1453         movel   d7,a4
1454         movel   a1,d7
1455         movel   a4,a1
1456 #endif
1457
1458 1:
1459 #ifndef __mcoldfire__
1460         exg     d7,a1           | put counter back in a1
1461 #else
1462         movel   d7,a4
1463         movel   a1,d7
1464         movel   a4,a1
1465 #endif
1466         addl    d3,d3           | shift sum once left
1467         addxl   d2,d2           |
1468         addxl   d1,d1           |
1469         addxl   d0,d0           |
1470         addl    d7,d7           |
1471         addxl   d6,d6           |
1472         bcc     2f              | if bit clear skip the following
1473 #ifndef __mcoldfire__
1474         exg     d7,a2           |
1475 #else
1476         movel   d7,a4
1477         movel   a2,d7
1478         movel   a4,a2
1479 #endif
1480         addl    d5,d3           | else add a to the sum
1481         addxl   d4,d2           |
1482         addxl   d7,d1           |
1483         addxl   d7,d0           |
1484 #ifndef __mcoldfire__
1485         exg     d7,a2           |
1486 #else
1487         movel   d7,a4
1488         movel   a2,d7
1489         movel   a4,a2
1490 #endif
1491 2:
1492 #ifndef __mcoldfire__
1493         exg     d7,a1           | put counter in d7
1494         dbf     d7,1b           | decrement and branch
1495 #else
1496         movel   d7,a4
1497         movel   a1,d7
1498         movel   a4,a1
1499         subql   IMM (1),d7
1500         bpl     1b
1501 #endif
1502
1503         movel   a3,d4           | restore exponent
1504 #ifndef __mcoldfire__
1505         moveml  sp@+,a2-a3
1506 #else
1507         movel   sp@+,a4
1508         movel   sp@+,a3
1509         movel   sp@+,a2
1510 #endif
1511
1512 | Now we have the product in d0-d1-d2-d3, with bit 8 of d0 set. The
1513 | first thing to do now is to normalize it so bit 8 becomes bit
1514 | DBL_MANT_DIG-32 (to do the rounding); later we will shift right.
1515         swap    d0
1516         swap    d1
1517         movew   d1,d0
1518         swap    d2
1519         movew   d2,d1
1520         swap    d3
1521         movew   d3,d2
1522         movew   IMM (0),d3
1523 #ifndef __mcoldfire__
1524         lsrl    IMM (1),d0
1525         roxrl   IMM (1),d1
1526         roxrl   IMM (1),d2
1527         roxrl   IMM (1),d3
1528         lsrl    IMM (1),d0
1529         roxrl   IMM (1),d1
1530         roxrl   IMM (1),d2
1531         roxrl   IMM (1),d3
1532         lsrl    IMM (1),d0
1533         roxrl   IMM (1),d1
1534         roxrl   IMM (1),d2
1535         roxrl   IMM (1),d3
1536 #else
1537         moveq   IMM (29),d6
1538         lsrl    IMM (3),d3
1539         movel   d2,d7
1540         lsll    d6,d7
1541         orl     d7,d3
1542         lsrl    IMM (3),d2
1543         movel   d1,d7
1544         lsll    d6,d7
1545         orl     d7,d2
1546         lsrl    IMM (3),d1
1547         movel   d0,d7
1548         lsll    d6,d7
1549         orl     d7,d1
1550         lsrl    IMM (3),d0
1551 #endif
1552
1553 | Now round, check for over- and underflow, and exit.
1554         movel   a0,d7           | get sign bit back into d7
1555         movew   IMM (MULTIPLY),d5
1556
1557         btst    IMM (DBL_MANT_DIG+1-32),d0
1558         beq     Lround$exit
1559 #ifndef __mcoldfire__
1560         lsrl    IMM (1),d0
1561         roxrl   IMM (1),d1
1562         addw    IMM (1),d4
1563 #else
1564         lsrl    IMM (1),d1
1565         btst    IMM (0),d0
1566         beq     10f
1567         bset    IMM (31),d1
1568 10:     lsrl    IMM (1),d0
1569         addl    IMM (1),d4
1570 #endif
1571         bra     Lround$exit
1572
1573 Lmuldf$inop:
1574         movew   IMM (MULTIPLY),d5
1575         bra     Ld$inop
1576
1577 Lmuldf$b$nf:
1578         movew   IMM (MULTIPLY),d5
1579         movel   a0,d7           | get sign bit back into d7
1580         tstl    d3              | we know d2 == 0x7ff00000, so check d3
1581         bne     Ld$inop         | if d3 <> 0 b is NaN
1582         bra     Ld$overflow     | else we have overflow (since a is finite)
1583
1584 Lmuldf$a$nf:
1585         movew   IMM (MULTIPLY),d5
1586         movel   a0,d7           | get sign bit back into d7
1587         tstl    d1              | we know d0 == 0x7ff00000, so check d1
1588         bne     Ld$inop         | if d1 <> 0 a is NaN
1589         bra     Ld$overflow     | else signal overflow
1590
1591 | If either number is zero return zero, unless the other is +/-INFINITY or
1592 | NaN, in which case we return NaN.
1593 Lmuldf$b$0:
1594         movew   IMM (MULTIPLY),d5
1595 #ifndef __mcoldfire__
1596         exg     d2,d0           | put b (==0) into d0-d1
1597         exg     d3,d1           | and a (with sign bit cleared) into d2-d3
1598 #else
1599         movel   d2,d7
1600         movel   d0,d2
1601         movel   d7,d0
1602         movel   d3,d7
1603         movel   d1,d3
1604         movel   d7,d1
1605 #endif
1606         bra     1f
1607 Lmuldf$a$0:
1608         movel   a6@(16),d2      | put b into d2-d3 again
1609         movel   a6@(20),d3      |
1610         bclr    IMM (31),d2     | clear sign bit
1611 1:      cmpl    IMM (0x7ff00000),d2 | check for non-finiteness
1612         bge     Ld$inop         | in case NaN or +/-INFINITY return NaN
1613         lea     SYM (_fpCCR),a0
1614         movew   IMM (0),a0@
1615 #ifndef __mcoldfire__
1616         moveml  sp@+,d2-d7
1617 #else
1618         moveml  sp@,d2-d7
1619         | XXX if frame pointer is ever removed, stack pointer must
1620         | be adjusted here.
1621 #endif
1622         unlk    a6
1623         rts
1624
1625 | If a number is denormalized we put an exponent of 1 but do not put the
1626 | hidden bit back into the fraction; instead we shift left until bit 21
1627 | (the hidden bit) is set, adjusting the exponent accordingly. We do this
1628 | to ensure that the product of the fractions is close to 1.
1629 Lmuldf$a$den:
1630         movel   IMM (1),d4
1631         andl    d6,d0
1632 1:      addl    d1,d1           | shift a left until bit 20 is set
1633         addxl   d0,d0           |
1634 #ifndef __mcoldfire__
1635         subw    IMM (1),d4      | and adjust exponent
1636 #else
1637         subl    IMM (1),d4      | and adjust exponent
1638 #endif
1639         btst    IMM (20),d0     |
1640         bne     Lmuldf$1        |
1641         bra     1b
1642
1643 Lmuldf$b$den:
1644         movel   IMM (1),d5
1645         andl    d6,d2
1646 1:      addl    d3,d3           | shift b left until bit 20 is set
1647         addxl   d2,d2           |
1648 #ifndef __mcoldfire__
1649         subw    IMM (1),d5      | and adjust exponent
1650 #else
1651         subql   IMM (1),d5      | and adjust exponent
1652 #endif
1653         btst    IMM (20),d2     |
1654         bne     Lmuldf$2        |
1655         bra     1b
1656
1657
1658 |=============================================================================
1659 |                              __divdf3
1660 |=============================================================================
1661
1662 | double __divdf3(double, double);
1663 SYM (__divdf3):
1664 #ifndef __mcoldfire__
1665         link    a6,IMM (0)
1666         moveml  d2-d7,sp@-
1667 #else
1668         link    a6,IMM (-24)
1669         moveml  d2-d7,sp@
1670 #endif
1671         movel   a6@(8),d0       | get a into d0-d1
1672         movel   a6@(12),d1      |
1673         movel   a6@(16),d2      | and b into d2-d3
1674         movel   a6@(20),d3      |
1675         movel   d0,d7           | d7 will hold the sign of the result
1676         eorl    d2,d7           |
1677         andl    IMM (0x80000000),d7
1678         movel   d7,a0           | save sign into a0
1679         movel   IMM (0x7ff00000),d7 | useful constant (+INFINITY)
1680         movel   d7,d6           | another (mask for fraction)
1681         notl    d6              |
1682         bclr    IMM (31),d0     | get rid of a's sign bit '
1683         movel   d0,d4           |
1684         orl     d1,d4           |
1685         beq     Ldivdf$a$0      | branch if a is zero
1686         movel   d0,d4           |
1687         bclr    IMM (31),d2     | get rid of b's sign bit '
1688         movel   d2,d5           |
1689         orl     d3,d5           |
1690         beq     Ldivdf$b$0      | branch if b is zero
1691         movel   d2,d5
1692         cmpl    d7,d0           | is a big?
1693         bhi     Ldivdf$inop     | if a is NaN return NaN
1694         beq     Ldivdf$a$nf     | if d0 == 0x7ff00000 we check d1
1695         cmpl    d7,d2           | now compare b with INFINITY
1696         bhi     Ldivdf$inop     | if b is NaN return NaN
1697         beq     Ldivdf$b$nf     | if d2 == 0x7ff00000 we check d3
1698 | Here we have both numbers finite and nonzero (and with no sign bit).
1699 | Now we get the exponents into d4 and d5 and normalize the numbers to
1700 | ensure that the ratio of the fractions is around 1. We do this by
1701 | making sure that both numbers have bit #DBL_MANT_DIG-32-1 (hidden bit)
1702 | set, even if they were denormalized to start with.
1703 | Thus, the result will satisfy: 2 > result > 1/2.
1704         andl    d7,d4           | and isolate exponent in d4
1705         beq     Ldivdf$a$den    | if exponent is zero we have a denormalized
1706         andl    d6,d0           | and isolate fraction
1707         orl     IMM (0x00100000),d0 | and put hidden bit back
1708         swap    d4              | I like exponents in the first byte
1709 #ifndef __mcoldfire__
1710         lsrw    IMM (4),d4      |
1711 #else
1712         lsrl    IMM (4),d4      |
1713 #endif
1714 Ldivdf$1:                       |
1715         andl    d7,d5           |
1716         beq     Ldivdf$b$den    |
1717         andl    d6,d2           |
1718         orl     IMM (0x00100000),d2
1719         swap    d5              |
1720 #ifndef __mcoldfire__
1721         lsrw    IMM (4),d5      |
1722 #else
1723         lsrl    IMM (4),d5      |
1724 #endif
1725 Ldivdf$2:                       |
1726 #ifndef __mcoldfire__
1727         subw    d5,d4           | subtract exponents
1728         addw    IMM (D_BIAS),d4 | and add bias
1729 #else
1730         subl    d5,d4           | subtract exponents
1731         addl    IMM (D_BIAS),d4 | and add bias
1732 #endif
1733
1734 | We are now ready to do the division. We have prepared things in such a way
1735 | that the ratio of the fractions will be less than 2 but greater than 1/2.
1736 | At this point the registers in use are:
1737 | d0-d1 hold a (first operand, bit DBL_MANT_DIG-32=0, bit
1738 | DBL_MANT_DIG-1-32=1)
1739 | d2-d3 hold b (second operand, bit DBL_MANT_DIG-32=1)
1740 | d4    holds the difference of the exponents, corrected by the bias
1741 | a0    holds the sign of the ratio
1742
1743 | To do the rounding correctly we need to keep information about the
1744 | nonsignificant bits. One way to do this would be to do the division
1745 | using four registers; another is to use two registers (as originally
1746 | I did), but use a sticky bit to preserve information about the
1747 | fractional part. Note that we can keep that info in a1, which is not
1748 | used.
1749         movel   IMM (0),d6      | d6-d7 will hold the result
1750         movel   d6,d7           |
1751         movel   IMM (0),a1      | and a1 will hold the sticky bit
1752
1753         movel   IMM (DBL_MANT_DIG-32+1),d5
1754
1755 1:      cmpl    d0,d2           | is a < b?
1756         bhi     3f              | if b > a skip the following
1757         beq     4f              | if d0==d2 check d1 and d3
1758 2:      subl    d3,d1           |
1759         subxl   d2,d0           | a <-- a - b
1760         bset    d5,d6           | set the corresponding bit in d6
1761 3:      addl    d1,d1           | shift a by 1
1762         addxl   d0,d0           |
1763 #ifndef __mcoldfire__
1764         dbra    d5,1b           | and branch back
1765 #else
1766         subql   IMM (1), d5
1767         bpl     1b
1768 #endif
1769         bra     5f
1770 4:      cmpl    d1,d3           | here d0==d2, so check d1 and d3
1771         bhi     3b              | if d1 > d2 skip the subtraction
1772         bra     2b              | else go do it
1773 5:
1774 | Here we have to start setting the bits in the second long.
1775         movel   IMM (31),d5     | again d5 is counter
1776
1777 1:      cmpl    d0,d2           | is a < b?
1778         bhi     3f              | if b > a skip the following
1779         beq     4f              | if d0==d2 check d1 and d3
1780 2:      subl    d3,d1           |
1781         subxl   d2,d0           | a <-- a - b
1782         bset    d5,d7           | set the corresponding bit in d7
1783 3:      addl    d1,d1           | shift a by 1
1784         addxl   d0,d0           |
1785 #ifndef __mcoldfire__
1786         dbra    d5,1b           | and branch back
1787 #else
1788         subql   IMM (1), d5
1789         bpl     1b
1790 #endif
1791         bra     5f
1792 4:      cmpl    d1,d3           | here d0==d2, so check d1 and d3
1793         bhi     3b              | if d1 > d2 skip the subtraction
1794         bra     2b              | else go do it
1795 5:
1796 | Now go ahead checking until we hit a one, which we store in d2.
1797         movel   IMM (DBL_MANT_DIG),d5
1798 1:      cmpl    d2,d0           | is a < b?
1799         bhi     4f              | if b < a, exit
1800         beq     3f              | if d0==d2 check d1 and d3
1801 2:      addl    d1,d1           | shift a by 1
1802         addxl   d0,d0           |
1803 #ifndef __mcoldfire__
1804         dbra    d5,1b           | and branch back
1805 #else
1806         subql   IMM (1), d5
1807         bpl     1b
1808 #endif
1809         movel   IMM (0),d2      | here no sticky bit was found
1810         movel   d2,d3
1811         bra     5f
1812 3:      cmpl    d1,d3           | here d0==d2, so check d1 and d3
1813         bhi     2b              | if d1 > d2 go back
1814 4:
1815 | Here put the sticky bit in d2-d3 (in the position which actually corresponds
1816 | to it; if you don't do this the algorithm loses in some cases). '
1817         movel   IMM (0),d2
1818         movel   d2,d3
1819 #ifndef __mcoldfire__
1820         subw    IMM (DBL_MANT_DIG),d5
1821         addw    IMM (63),d5
1822         cmpw    IMM (31),d5
1823 #else
1824         subl    IMM (DBL_MANT_DIG),d5
1825         addl    IMM (63),d5
1826         cmpl    IMM (31),d5
1827 #endif
1828         bhi     2f
1829 1:      bset    d5,d3
1830         bra     5f
1831 #ifndef __mcoldfire__
1832         subw    IMM (32),d5
1833 #else
1834         subl    IMM (32),d5
1835 #endif
1836 2:      bset    d5,d2
1837 5:
1838 | Finally we are finished! Move the longs in the address registers to
1839 | their final destination:
1840         movel   d6,d0
1841         movel   d7,d1
1842         movel   IMM (0),d3
1843
1844 | Here we have finished the division, with the result in d0-d1-d2-d3, with
1845 | 2^21 <= d6 < 2^23. Thus bit 23 is not set, but bit 22 could be set.
1846 | If it is not, then definitely bit 21 is set. Normalize so bit 22 is
1847 | not set:
1848         btst    IMM (DBL_MANT_DIG-32+1),d0
1849         beq     1f
1850 #ifndef __mcoldfire__
1851         lsrl    IMM (1),d0
1852         roxrl   IMM (1),d1
1853         roxrl   IMM (1),d2
1854         roxrl   IMM (1),d3
1855         addw    IMM (1),d4
1856 #else
1857         lsrl    IMM (1),d3
1858         btst    IMM (0),d2
1859         beq     10f
1860         bset    IMM (31),d3
1861 10:     lsrl    IMM (1),d2
1862         btst    IMM (0),d1
1863         beq     11f
1864         bset    IMM (31),d2
1865 11:     lsrl    IMM (1),d1
1866         btst    IMM (0),d0
1867         beq     12f
1868         bset    IMM (31),d1
1869 12:     lsrl    IMM (1),d0
1870         addl    IMM (1),d4
1871 #endif
1872 1:
1873 | Now round, check for over- and underflow, and exit.
1874         movel   a0,d7           | restore sign bit to d7
1875         movew   IMM (DIVIDE),d5
1876         bra     Lround$exit
1877
1878 Ldivdf$inop:
1879         movew   IMM (DIVIDE),d5
1880         bra     Ld$inop
1881
1882 Ldivdf$a$0:
1883 | If a is zero check to see whether b is zero also. In that case return
1884 | NaN; then check if b is NaN, and return NaN also in that case. Else
1885 | return zero.
1886         movew   IMM (DIVIDE),d5
1887         bclr    IMM (31),d2     |
1888         movel   d2,d4           |
1889         orl     d3,d4           |
1890         beq     Ld$inop         | if b is also zero return NaN
1891         cmpl    IMM (0x7ff00000),d2 | check for NaN
1892         bhi     Ld$inop         |
1893         blt     1f              |
1894         tstl    d3              |
1895         bne     Ld$inop         |
1896 1:      movel   IMM (0),d0      | else return zero
1897         movel   d0,d1           |
1898         lea     SYM (_fpCCR),a0 | clear exception flags
1899         movew   IMM (0),a0@     |
1900 #ifndef __mcoldfire__
1901         moveml  sp@+,d2-d7      |
1902 #else
1903         moveml  sp@,d2-d7       |
1904         | XXX if frame pointer is ever removed, stack pointer must
1905         | be adjusted here.
1906 #endif
1907         unlk    a6              |
1908         rts                     |
1909
1910 Ldivdf$b$0:
1911         movew   IMM (DIVIDE),d5
1912 | If we got here a is not zero. Check if a is NaN; in that case return NaN,
1913 | else return +/-INFINITY. Remember that a is in d0 with the sign bit
1914 | cleared already.
1915         movel   a0,d7           | put a's sign bit back in d7 '
1916         cmpl    IMM (0x7ff00000),d0 | compare d0 with INFINITY
1917         bhi     Ld$inop         | if larger it is NaN
1918         tstl    d1              |
1919         bne     Ld$inop         |
1920         bra     Ld$div$0        | else signal DIVIDE_BY_ZERO
1921
1922 Ldivdf$b$nf:
1923         movew   IMM (DIVIDE),d5
1924 | If d2 == 0x7ff00000 we have to check d3.
1925         tstl    d3              |
1926         bne     Ld$inop         | if d3 <> 0, b is NaN
1927         bra     Ld$underflow    | else b is +/-INFINITY, so signal underflow
1928
1929 Ldivdf$a$nf:
1930         movew   IMM (DIVIDE),d5
1931 | If d0 == 0x7ff00000 we have to check d1.
1932         tstl    d1              |
1933         bne     Ld$inop         | if d1 <> 0, a is NaN
1934 | If a is INFINITY we have to check b
1935         cmpl    d7,d2           | compare b with INFINITY
1936         bge     Ld$inop         | if b is NaN or INFINITY return NaN
1937         tstl    d3              |
1938         bne     Ld$inop         |
1939         bra     Ld$overflow     | else return overflow
1940
1941 | If a number is denormalized we put an exponent of 1 but do not put the
1942 | bit back into the fraction.
1943 Ldivdf$a$den:
1944         movel   IMM (1),d4
1945         andl    d6,d0
1946 1:      addl    d1,d1           | shift a left until bit 20 is set
1947         addxl   d0,d0
1948 #ifndef __mcoldfire__
1949         subw    IMM (1),d4      | and adjust exponent
1950 #else
1951         subl    IMM (1),d4      | and adjust exponent
1952 #endif
1953         btst    IMM (DBL_MANT_DIG-32-1),d0
1954         bne     Ldivdf$1
1955         bra     1b
1956
1957 Ldivdf$b$den:
1958         movel   IMM (1),d5
1959         andl    d6,d2
1960 1:      addl    d3,d3           | shift b left until bit 20 is set
1961         addxl   d2,d2
1962 #ifndef __mcoldfire__
1963         subw    IMM (1),d5      | and adjust exponent
1964 #else
1965         subql   IMM (1),d5      | and adjust exponent
1966 #endif
1967         btst    IMM (DBL_MANT_DIG-32-1),d2
1968         bne     Ldivdf$2
1969         bra     1b
1970
1971 Lround$exit:
1972 | This is a common exit point for __muldf3 and __divdf3. When they enter
1973 | this point the sign of the result is in d7, the result in d0-d1, normalized
1974 | so that 2^21 <= d0 < 2^22, and the exponent is in the lower byte of d4.
1975
1976 | First check for underlow in the exponent:
1977 #ifndef __mcoldfire__
1978         cmpw    IMM (-DBL_MANT_DIG-1),d4
1979 #else
1980         cmpl    IMM (-DBL_MANT_DIG-1),d4
1981 #endif
1982         blt     Ld$underflow
1983 | It could happen that the exponent is less than 1, in which case the
1984 | number is denormalized. In this case we shift right and adjust the
1985 | exponent until it becomes 1 or the fraction is zero (in the latter case
1986 | we signal underflow and return zero).
1987         movel   d7,a0           |
1988         movel   IMM (0),d6      | use d6-d7 to collect bits flushed right
1989         movel   d6,d7           | use d6-d7 to collect bits flushed right
1990 #ifndef __mcoldfire__
1991         cmpw    IMM (1),d4      | if the exponent is less than 1 we
1992 #else
1993         cmpl    IMM (1),d4      | if the exponent is less than 1 we
1994 #endif
1995         bge     2f              | have to shift right (denormalize)
1996 1:
1997 #ifndef __mcoldfire__
1998         addw    IMM (1),d4      | adjust the exponent
1999         lsrl    IMM (1),d0      | shift right once
2000         roxrl   IMM (1),d1      |
2001         roxrl   IMM (1),d2      |
2002         roxrl   IMM (1),d3      |
2003         roxrl   IMM (1),d6      |
2004         roxrl   IMM (1),d7      |
2005         cmpw    IMM (1),d4      | is the exponent 1 already?
2006 #else
2007         addl    IMM (1),d4      | adjust the exponent
2008         lsrl    IMM (1),d7
2009         btst    IMM (0),d6
2010         beq     13f
2011         bset    IMM (31),d7
2012 13:     lsrl    IMM (1),d6
2013         btst    IMM (0),d3
2014         beq     14f
2015         bset    IMM (31),d6
2016 14:     lsrl    IMM (1),d3
2017         btst    IMM (0),d2
2018         beq     10f
2019         bset    IMM (31),d3
2020 10:     lsrl    IMM (1),d2
2021         btst    IMM (0),d1
2022         beq     11f
2023         bset    IMM (31),d2
2024 11:     lsrl    IMM (1),d1
2025         btst    IMM (0),d0
2026         beq     12f
2027         bset    IMM (31),d1
2028 12:     lsrl    IMM (1),d0
2029         cmpl    IMM (1),d4      | is the exponent 1 already?
2030 #endif
2031         beq     2f              | if not loop back
2032         bra     1b              |
2033         bra     Ld$underflow    | safety check, shouldn't execute '
2034 2:      orl     d6,d2           | this is a trick so we don't lose  '
2035         orl     d7,d3           | the bits which were flushed right
2036         movel   a0,d7           | get back sign bit into d7
2037 | Now call the rounding routine (which takes care of denormalized numbers):
2038         lea     Lround$0,a0     | to return from rounding routine
2039         lea     SYM (_fpCCR),a1 | check the rounding mode
2040 #ifdef __mcoldfire__
2041         clrl    d6
2042 #endif
2043         movew   a1@(6),d6       | rounding mode in d6
2044         beq     Lround$to$nearest
2045 #ifndef __mcoldfire__
2046         cmpw    IMM (ROUND_TO_PLUS),d6
2047 #else
2048         cmpl    IMM (ROUND_TO_PLUS),d6
2049 #endif
2050         bhi     Lround$to$minus
2051         blt     Lround$to$zero
2052         bra     Lround$to$plus
2053 Lround$0:
2054 | Here we have a correctly rounded result (either normalized or denormalized).
2055
2056 | Here we should have either a normalized number or a denormalized one, and
2057 | the exponent is necessarily larger or equal to 1 (so we don't have to  '
2058 | check again for underflow!). We have to check for overflow or for a
2059 | denormalized number (which also signals underflow).
2060 | Check for overflow (i.e., exponent >= 0x7ff).
2061 #ifndef __mcoldfire__
2062         cmpw    IMM (0x07ff),d4
2063 #else
2064         cmpl    IMM (0x07ff),d4
2065 #endif
2066         bge     Ld$overflow
2067 | Now check for a denormalized number (exponent==0):
2068         movew   d4,d4
2069         beq     Ld$den
2070 1:
2071 | Put back the exponents and sign and return.
2072 #ifndef __mcoldfire__
2073         lslw    IMM (4),d4      | exponent back to fourth byte
2074 #else
2075         lsll    IMM (4),d4      | exponent back to fourth byte
2076 #endif
2077         bclr    IMM (DBL_MANT_DIG-32-1),d0
2078         swap    d0              | and put back exponent
2079 #ifndef __mcoldfire__
2080         orw     d4,d0           |
2081 #else
2082         orl     d4,d0           |
2083 #endif
2084         swap    d0              |
2085         orl     d7,d0           | and sign also
2086
2087         lea     SYM (_fpCCR),a0
2088         movew   IMM (0),a0@
2089 #ifndef __mcoldfire__
2090         moveml  sp@+,d2-d7
2091 #else
2092         moveml  sp@,d2-d7
2093         | XXX if frame pointer is ever removed, stack pointer must
2094         | be adjusted here.
2095 #endif
2096         unlk    a6
2097         rts
2098
2099 |=============================================================================
2100 |                              __negdf2
2101 |=============================================================================
2102
2103 | double __negdf2(double, double);
2104 SYM (__negdf2):
2105 #ifndef __mcoldfire__
2106         link    a6,IMM (0)
2107         moveml  d2-d7,sp@-
2108 #else
2109         link    a6,IMM (-24)
2110         moveml  d2-d7,sp@
2111 #endif
2112         movew   IMM (NEGATE),d5
2113         movel   a6@(8),d0       | get number to negate in d0-d1
2114         movel   a6@(12),d1      |
2115         bchg    IMM (31),d0     | negate
2116         movel   d0,d2           | make a positive copy (for the tests)
2117         bclr    IMM (31),d2     |
2118         movel   d2,d4           | check for zero
2119         orl     d1,d4           |
2120         beq     2f              | if zero (either sign) return +zero
2121         cmpl    IMM (0x7ff00000),d2 | compare to +INFINITY
2122         blt     1f              | if finite, return
2123         bhi     Ld$inop         | if larger (fraction not zero) is NaN
2124         tstl    d1              | if d2 == 0x7ff00000 check d1
2125         bne     Ld$inop         |
2126         movel   d0,d7           | else get sign and return INFINITY
2127         andl    IMM (0x80000000),d7
2128         bra     Ld$infty
2129 1:      lea     SYM (_fpCCR),a0
2130         movew   IMM (0),a0@
2131 #ifndef __mcoldfire__
2132         moveml  sp@+,d2-d7
2133 #else
2134         moveml  sp@,d2-d7
2135         | XXX if frame pointer is ever removed, stack pointer must
2136         | be adjusted here.
2137 #endif
2138         unlk    a6
2139         rts
2140 2:      bclr    IMM (31),d0
2141         bra     1b
2142
2143 |=============================================================================
2144 |                              __cmpdf2
2145 |=============================================================================
2146
2147 GREATER =  1
2148 LESS    = -1
2149 EQUAL   =  0
2150
2151 | int __cmpdf2(double, double);
2152 SYM (__cmpdf2):
2153 #ifndef __mcoldfire__
2154         link    a6,IMM (0)
2155         moveml  d2-d7,sp@-      | save registers
2156 #else
2157         link    a6,IMM (-24)
2158         moveml  d2-d7,sp@
2159 #endif
2160         movew   IMM (COMPARE),d5
2161         movel   a6@(8),d0       | get first operand
2162         movel   a6@(12),d1      |
2163         movel   a6@(16),d2      | get second operand
2164         movel   a6@(20),d3      |
2165 | First check if a and/or b are (+/-) zero and in that case clear
2166 | the sign bit.
2167         movel   d0,d6           | copy signs into d6 (a) and d7(b)
2168         bclr    IMM (31),d0     | and clear signs in d0 and d2
2169         movel   d2,d7           |
2170         bclr    IMM (31),d2     |
2171         cmpl    IMM (0x7fff0000),d0 | check for a == NaN
2172         bhi     Ld$inop         | if d0 > 0x7ff00000, a is NaN
2173         beq     Lcmpdf$a$nf     | if equal can be INFINITY, so check d1
2174         movel   d0,d4           | copy into d4 to test for zero
2175         orl     d1,d4           |
2176         beq     Lcmpdf$a$0      |
2177 Lcmpdf$0:
2178         cmpl    IMM (0x7fff0000),d2 | check for b == NaN
2179         bhi     Ld$inop         | if d2 > 0x7ff00000, b is NaN
2180         beq     Lcmpdf$b$nf     | if equal can be INFINITY, so check d3
2181         movel   d2,d4           |
2182         orl     d3,d4           |
2183         beq     Lcmpdf$b$0      |
2184 Lcmpdf$1:
2185 | Check the signs
2186         eorl    d6,d7
2187         bpl     1f
2188 | If the signs are not equal check if a >= 0
2189         tstl    d6
2190         bpl     Lcmpdf$a$gt$b   | if (a >= 0 && b < 0) => a > b
2191         bmi     Lcmpdf$b$gt$a   | if (a < 0 && b >= 0) => a < b
2192 1:
2193 | If the signs are equal check for < 0
2194         tstl    d6
2195         bpl     1f
2196 | If both are negative exchange them
2197 #ifndef __mcoldfire__
2198         exg     d0,d2
2199         exg     d1,d3
2200 #else
2201         movel   d0,d7
2202         movel   d2,d0
2203         movel   d7,d2
2204         movel   d1,d7
2205         movel   d3,d1
2206         movel   d7,d3
2207 #endif
2208 1:
2209 | Now that they are positive we just compare them as longs (does this also
2210 | work for denormalized numbers?).
2211         cmpl    d0,d2
2212         bhi     Lcmpdf$b$gt$a   | |b| > |a|
2213         bne     Lcmpdf$a$gt$b   | |b| < |a|
2214 | If we got here d0 == d2, so we compare d1 and d3.
2215         cmpl    d1,d3
2216         bhi     Lcmpdf$b$gt$a   | |b| > |a|
2217         bne     Lcmpdf$a$gt$b   | |b| < |a|
2218 | If we got here a == b.
2219         movel   IMM (EQUAL),d0
2220 #ifndef __mcoldfire__
2221         moveml  sp@+,d2-d7      | put back the registers
2222 #else
2223         moveml  sp@,d2-d7
2224         | XXX if frame pointer is ever removed, stack pointer must
2225         | be adjusted here.
2226 #endif
2227         unlk    a6
2228         rts
2229 Lcmpdf$a$gt$b:
2230         movel   IMM (GREATER),d0
2231 #ifndef __mcoldfire__
2232         moveml  sp@+,d2-d7      | put back the registers
2233 #else
2234         moveml  sp@,d2-d7
2235         | XXX if frame pointer is ever removed, stack pointer must
2236         | be adjusted here.
2237 #endif
2238         unlk    a6
2239         rts
2240 Lcmpdf$b$gt$a:
2241         movel   IMM (LESS),d0
2242 #ifndef __mcoldfire__
2243         moveml  sp@+,d2-d7      | put back the registers
2244 #else
2245         moveml  sp@,d2-d7
2246         | XXX if frame pointer is ever removed, stack pointer must
2247         | be adjusted here.
2248 #endif
2249         unlk    a6
2250         rts
2251
2252 Lcmpdf$a$0:
2253         bclr    IMM (31),d6
2254         bra     Lcmpdf$0
2255 Lcmpdf$b$0:
2256         bclr    IMM (31),d7
2257         bra     Lcmpdf$1
2258
2259 Lcmpdf$a$nf:
2260         tstl    d1
2261         bne     Ld$inop
2262         bra     Lcmpdf$0
2263
2264 Lcmpdf$b$nf:
2265         tstl    d3
2266         bne     Ld$inop
2267         bra     Lcmpdf$1
2268
2269 |=============================================================================
2270 |                           rounding routines
2271 |=============================================================================
2272
2273 | The rounding routines expect the number to be normalized in registers
2274 | d0-d1-d2-d3, with the exponent in register d4. They assume that the
2275 | exponent is larger or equal to 1. They return a properly normalized number
2276 | if possible, and a denormalized number otherwise. The exponent is returned
2277 | in d4.
2278
2279 Lround$to$nearest:
2280 | We now normalize as suggested by D. Knuth ("Seminumerical Algorithms"):
2281 | Here we assume that the exponent is not too small (this should be checked
2282 | before entering the rounding routine), but the number could be denormalized.
2283
2284 | Check for denormalized numbers:
2285 1:      btst    IMM (DBL_MANT_DIG-32),d0
2286         bne     2f              | if set the number is normalized
2287 | Normalize shifting left until bit #DBL_MANT_DIG-32 is set or the exponent
2288 | is one (remember that a denormalized number corresponds to an
2289 | exponent of -D_BIAS+1).
2290 #ifndef __mcoldfire__
2291         cmpw    IMM (1),d4      | remember that the exponent is at least one
2292 #else
2293         cmpl    IMM (1),d4      | remember that the exponent is at least one
2294 #endif
2295         beq     2f              | an exponent of one means denormalized
2296         addl    d3,d3           | else shift and adjust the exponent
2297         addxl   d2,d2           |
2298         addxl   d1,d1           |
2299         addxl   d0,d0           |
2300 #ifndef __mcoldfire__
2301         dbra    d4,1b           |
2302 #else
2303         subql   IMM (1), d4
2304         bpl     1b
2305 #endif
2306 2:
2307 | Now round: we do it as follows: after the shifting we can write the
2308 | fraction part as f + delta, where 1 < f < 2^25, and 0 <= delta <= 2.
2309 | If delta < 1, do nothing. If delta > 1, add 1 to f.
2310 | If delta == 1, we make sure the rounded number will be even (odd?)
2311 | (after shifting).
2312         btst    IMM (0),d1      | is delta < 1?
2313         beq     2f              | if so, do not do anything
2314         orl     d2,d3           | is delta == 1?
2315         bne     1f              | if so round to even
2316         movel   d1,d3           |
2317         andl    IMM (2),d3      | bit 1 is the last significant bit
2318         movel   IMM (0),d2      |
2319         addl    d3,d1           |
2320         addxl   d2,d0           |
2321         bra     2f              |
2322 1:      movel   IMM (1),d3      | else add 1
2323         movel   IMM (0),d2      |
2324         addl    d3,d1           |
2325         addxl   d2,d0
2326 | Shift right once (because we used bit #DBL_MANT_DIG-32!).
2327 2:
2328 #ifndef __mcoldfire__
2329         lsrl    IMM (1),d0
2330         roxrl   IMM (1),d1
2331 #else
2332         lsrl    IMM (1),d1
2333         btst    IMM (0),d0
2334         beq     10f
2335         bset    IMM (31),d1
2336 10:     lsrl    IMM (1),d0
2337 #endif
2338
2339 | Now check again bit #DBL_MANT_DIG-32 (rounding could have produced a
2340 | 'fraction overflow' ...).
2341         btst    IMM (DBL_MANT_DIG-32),d0
2342         beq     1f
2343 #ifndef __mcoldfire__
2344         lsrl    IMM (1),d0
2345         roxrl   IMM (1),d1
2346         addw    IMM (1),d4
2347 #else
2348         lsrl    IMM (1),d1
2349         btst    IMM (0),d0
2350         beq     10f
2351         bset    IMM (31),d1
2352 10:     lsrl    IMM (1),d0
2353         addl    IMM (1),d4
2354 #endif
2355 1:
2356 | If bit #DBL_MANT_DIG-32-1 is clear we have a denormalized number, so we
2357 | have to put the exponent to zero and return a denormalized number.
2358         btst    IMM (DBL_MANT_DIG-32-1),d0
2359         beq     1f
2360         jmp     a0@
2361 1:      movel   IMM (0),d4
2362         jmp     a0@
2363
2364 Lround$to$zero:
2365 Lround$to$plus:
2366 Lround$to$minus:
2367         jmp     a0@
2368 #endif /* L_double */
2369
2370 #ifdef  L_float
2371
2372         .globl  SYM (_fpCCR)
2373         .globl  $_exception_handler
2374
2375 QUIET_NaN    = 0xffffffff
2376 SIGNL_NaN    = 0x7f800001
2377 INFINITY     = 0x7f800000
2378
2379 F_MAX_EXP      = 0xff
2380 F_BIAS         = 126
2381 FLT_MAX_EXP    = F_MAX_EXP - F_BIAS
2382 FLT_MIN_EXP    = 1 - F_BIAS
2383 FLT_MANT_DIG   = 24
2384
2385 INEXACT_RESULT          = 0x0001
2386 UNDERFLOW               = 0x0002
2387 OVERFLOW                = 0x0004
2388 DIVIDE_BY_ZERO          = 0x0008
2389 INVALID_OPERATION       = 0x0010
2390
2391 SINGLE_FLOAT = 1
2392
2393 NOOP         = 0
2394 ADD          = 1
2395 MULTIPLY     = 2
2396 DIVIDE       = 3
2397 NEGATE       = 4
2398 COMPARE      = 5
2399 EXTENDSFDF   = 6
2400 TRUNCDFSF    = 7
2401
2402 UNKNOWN           = -1
2403 ROUND_TO_NEAREST  = 0 | round result to nearest representable value
2404 ROUND_TO_ZERO     = 1 | round result towards zero
2405 ROUND_TO_PLUS     = 2 | round result towards plus infinity
2406 ROUND_TO_MINUS    = 3 | round result towards minus infinity
2407
2408 | Entry points:
2409
2410         .globl SYM (__addsf3)
2411         .globl SYM (__subsf3)
2412         .globl SYM (__mulsf3)
2413         .globl SYM (__divsf3)
2414         .globl SYM (__negsf2)
2415         .globl SYM (__cmpsf2)
2416
2417 | These are common routines to return and signal exceptions.
2418
2419         .text
2420         .even
2421
2422 Lf$den:
2423 | Return and signal a denormalized number
2424         orl     d7,d0
2425         movew   IMM (INEXACT_RESULT+UNDERFLOW),d7
2426         moveq   IMM (SINGLE_FLOAT),d6
2427         jmp     $_exception_handler
2428
2429 Lf$infty:
2430 Lf$overflow:
2431 | Return a properly signed INFINITY and set the exception flags
2432         movel   IMM (INFINITY),d0
2433         orl     d7,d0
2434         movew   IMM (INEXACT_RESULT+OVERFLOW),d7
2435         moveq   IMM (SINGLE_FLOAT),d6
2436         jmp     $_exception_handler
2437
2438 Lf$underflow:
2439 | Return 0 and set the exception flags
2440         movel   IMM (0),d0
2441         movew   IMM (INEXACT_RESULT+UNDERFLOW),d7
2442         moveq   IMM (SINGLE_FLOAT),d6
2443         jmp     $_exception_handler
2444
2445 Lf$inop:
2446 | Return a quiet NaN and set the exception flags
2447         movel   IMM (QUIET_NaN),d0
2448         movew   IMM (INEXACT_RESULT+INVALID_OPERATION),d7
2449         moveq   IMM (SINGLE_FLOAT),d6
2450         jmp     $_exception_handler
2451
2452 Lf$div$0:
2453 | Return a properly signed INFINITY and set the exception flags
2454         movel   IMM (INFINITY),d0
2455         orl     d7,d0
2456         movew   IMM (INEXACT_RESULT+DIVIDE_BY_ZERO),d7
2457         moveq   IMM (SINGLE_FLOAT),d6
2458         jmp     $_exception_handler
2459
2460 |=============================================================================
2461 |=============================================================================
2462 |                         single precision routines
2463 |=============================================================================
2464 |=============================================================================
2465
2466 | A single precision floating point number (float) has the format:
2467 |
2468 | struct _float {
2469 |  unsigned int sign      : 1;  /* sign bit */
2470 |  unsigned int exponent  : 8;  /* exponent, shifted by 126 */
2471 |  unsigned int fraction  : 23; /* fraction */
2472 | } float;
2473 |
2474 | Thus sizeof(float) = 4 (32 bits).
2475 |
2476 | All the routines are callable from C programs, and return the result
2477 | in the single register d0. They also preserve all registers except
2478 | d0-d1 and a0-a1.
2479
2480 |=============================================================================
2481 |                              __subsf3
2482 |=============================================================================
2483
2484 | float __subsf3(float, float);
2485 SYM (__subsf3):
2486         bchg    IMM (31),sp@(8) | change sign of second operand
2487                                 | and fall through
2488 |=============================================================================
2489 |                              __addsf3
2490 |=============================================================================
2491
2492 | float __addsf3(float, float);
2493 SYM (__addsf3):
2494 #ifndef __mcoldfire__
2495         link    a6,IMM (0)      | everything will be done in registers
2496         moveml  d2-d7,sp@-      | save all data registers but d0-d1
2497 #else
2498         link    a6,IMM (-24)
2499         moveml  d2-d7,sp@
2500 #endif
2501         movel   a6@(8),d0       | get first operand
2502         movel   a6@(12),d1      | get second operand
2503         movel   d0,d6           | get d0's sign bit '
2504         addl    d0,d0           | check and clear sign bit of a
2505         beq     Laddsf$b        | if zero return second operand
2506         movel   d1,d7           | save b's sign bit '
2507         addl    d1,d1           | get rid of sign bit
2508         beq     Laddsf$a        | if zero return first operand
2509
2510         movel   d6,a0           | save signs in address registers
2511         movel   d7,a1           | so we can use d6 and d7
2512
2513 | Get the exponents and check for denormalized and/or infinity.
2514
2515         movel   IMM (0x00ffffff),d4     | mask to get fraction
2516         movel   IMM (0x01000000),d5     | mask to put hidden bit back
2517
2518         movel   d0,d6           | save a to get exponent
2519         andl    d4,d0           | get fraction in d0
2520         notl    d4              | make d4 into a mask for the exponent
2521         andl    d4,d6           | get exponent in d6
2522         beq     Laddsf$a$den    | branch if a is denormalized
2523         cmpl    d4,d6           | check for INFINITY or NaN
2524         beq     Laddsf$nf
2525         swap    d6              | put exponent into first word
2526         orl     d5,d0           | and put hidden bit back
2527 Laddsf$1:
2528 | Now we have a's exponent in d6 (second byte) and the mantissa in d0. '
2529         movel   d1,d7           | get exponent in d7
2530         andl    d4,d7           |
2531         beq     Laddsf$b$den    | branch if b is denormalized
2532         cmpl    d4,d7           | check for INFINITY or NaN
2533         beq     Laddsf$nf
2534         swap    d7              | put exponent into first word
2535         notl    d4              | make d4 into a mask for the fraction
2536         andl    d4,d1           | get fraction in d1
2537         orl     d5,d1           | and put hidden bit back
2538 Laddsf$2:
2539 | Now we have b's exponent in d7 (second byte) and the mantissa in d1. '
2540
2541 | Note that the hidden bit corresponds to bit #FLT_MANT_DIG-1, and we
2542 | shifted right once, so bit #FLT_MANT_DIG is set (so we have one extra
2543 | bit).
2544
2545         movel   d1,d2           | move b to d2, since we want to use
2546                                 | two registers to do the sum
2547         movel   IMM (0),d1      | and clear the new ones
2548         movel   d1,d3           |
2549
2550 | Here we shift the numbers in registers d0 and d1 so the exponents are the
2551 | same, and put the largest exponent in d6. Note that we are using two
2552 | registers for each number (see the discussion by D. Knuth in "Seminumerical
2553 | Algorithms").
2554 #ifndef __mcoldfire__
2555         cmpw    d6,d7           | compare exponents
2556 #else
2557         cmpl    d6,d7           | compare exponents
2558 #endif
2559         beq     Laddsf$3        | if equal don't shift '
2560         bhi     5f              | branch if second exponent largest
2561 1:
2562         subl    d6,d7           | keep the largest exponent
2563         negl    d7
2564 #ifndef __mcoldfire__
2565         lsrw    IMM (8),d7      | put difference in lower byte
2566 #else
2567         lsrl    IMM (8),d7      | put difference in lower byte
2568 #endif
2569 | if difference is too large we don't shift (actually, we can just exit) '
2570 #ifndef __mcoldfire__
2571         cmpw    IMM (FLT_MANT_DIG+2),d7
2572 #else
2573         cmpl    IMM (FLT_MANT_DIG+2),d7
2574 #endif
2575         bge     Laddsf$b$small
2576 #ifndef __mcoldfire__
2577         cmpw    IMM (16),d7     | if difference >= 16 swap
2578 #else
2579         cmpl    IMM (16),d7     | if difference >= 16 swap
2580 #endif
2581         bge     4f
2582 2:
2583 #ifndef __mcoldfire__
2584         subw    IMM (1),d7
2585 #else
2586         subql   IMM (1), d7
2587 #endif
2588 3:
2589 #ifndef __mcoldfire__
2590         lsrl    IMM (1),d2      | shift right second operand
2591         roxrl   IMM (1),d3
2592         dbra    d7,3b
2593 #else
2594         lsrl    IMM (1),d3
2595         btst    IMM (0),d2
2596         beq     10f
2597         bset    IMM (31),d3
2598 10:     lsrl    IMM (1),d2
2599         subql   IMM (1), d7
2600         bpl     3b
2601 #endif
2602         bra     Laddsf$3
2603 4:
2604         movew   d2,d3
2605         swap    d3
2606         movew   d3,d2
2607         swap    d2
2608 #ifndef __mcoldfire__
2609         subw    IMM (16),d7
2610 #else
2611         subl    IMM (16),d7
2612 #endif
2613         bne     2b              | if still more bits, go back to normal case
2614         bra     Laddsf$3
2615 5:
2616 #ifndef __mcoldfire__
2617         exg     d6,d7           | exchange the exponents
2618 #else
2619         eorl    d6,d7
2620         eorl    d7,d6
2621         eorl    d6,d7
2622 #endif
2623         subl    d6,d7           | keep the largest exponent
2624         negl    d7              |
2625 #ifndef __mcoldfire__
2626         lsrw    IMM (8),d7      | put difference in lower byte
2627 #else
2628         lsrl    IMM (8),d7      | put difference in lower byte
2629 #endif
2630 | if difference is too large we don't shift (and exit!) '
2631 #ifndef __mcoldfire__
2632         cmpw    IMM (FLT_MANT_DIG+2),d7
2633 #else
2634         cmpl    IMM (FLT_MANT_DIG+2),d7
2635 #endif
2636         bge     Laddsf$a$small
2637 #ifndef __mcoldfire__
2638         cmpw    IMM (16),d7     | if difference >= 16 swap
2639 #else
2640         cmpl    IMM (16),d7     | if difference >= 16 swap
2641 #endif
2642         bge     8f
2643 6:
2644 #ifndef __mcoldfire__
2645         subw    IMM (1),d7
2646 #else
2647         subl    IMM (1),d7
2648 #endif
2649 7:
2650 #ifndef __mcoldfire__
2651         lsrl    IMM (1),d0      | shift right first operand
2652         roxrl   IMM (1),d1
2653         dbra    d7,7b
2654 #else
2655         lsrl    IMM (1),d1
2656         btst    IMM (0),d0
2657         beq     10f
2658         bset    IMM (31),d1
2659 10:     lsrl    IMM (1),d0
2660         subql   IMM (1),d7
2661         bpl     7b
2662 #endif
2663         bra     Laddsf$3
2664 8:
2665         movew   d0,d1
2666         swap    d1
2667         movew   d1,d0
2668         swap    d0
2669 #ifndef __mcoldfire__
2670         subw    IMM (16),d7
2671 #else
2672         subl    IMM (16),d7
2673 #endif
2674         bne     6b              | if still more bits, go back to normal case
2675                                 | otherwise we fall through
2676
2677 | Now we have a in d0-d1, b in d2-d3, and the largest exponent in d6 (the
2678 | signs are stored in a0 and a1).
2679
2680 Laddsf$3:
2681 | Here we have to decide whether to add or subtract the numbers
2682 #ifndef __mcoldfire__
2683         exg     d6,a0           | get signs back
2684         exg     d7,a1           | and save the exponents
2685 #else
2686         movel   d6,d4
2687         movel   a0,d6
2688         movel   d4,a0
2689         movel   d7,d4
2690         movel   a1,d7
2691         movel   d4,a1
2692 #endif
2693         eorl    d6,d7           | combine sign bits
2694         bmi     Lsubsf$0        | if negative a and b have opposite
2695                                 | sign so we actually subtract the
2696                                 | numbers
2697
2698 | Here we have both positive or both negative
2699 #ifndef __mcoldfire__
2700         exg     d6,a0           | now we have the exponent in d6
2701 #else
2702         movel   d6,d4
2703         movel   a0,d6
2704         movel   d4,a0
2705 #endif
2706         movel   a0,d7           | and sign in d7
2707         andl    IMM (0x80000000),d7
2708 | Here we do the addition.
2709         addl    d3,d1
2710         addxl   d2,d0
2711 | Note: now we have d2, d3, d4 and d5 to play with!
2712
2713 | Put the exponent, in the first byte, in d2, to use the "standard" rounding
2714 | routines:
2715         movel   d6,d2
2716 #ifndef __mcoldfire__
2717         lsrw    IMM (8),d2
2718 #else
2719         lsrl    IMM (8),d2
2720 #endif
2721
2722 | Before rounding normalize so bit #FLT_MANT_DIG is set (we will consider
2723 | the case of denormalized numbers in the rounding routine itself).
2724 | As in the addition (not in the subtraction!) we could have set
2725 | one more bit we check this:
2726         btst    IMM (FLT_MANT_DIG+1),d0
2727         beq     1f
2728 #ifndef __mcoldfire__
2729         lsrl    IMM (1),d0
2730         roxrl   IMM (1),d1
2731 #else
2732         lsrl    IMM (1),d1
2733         btst    IMM (0),d0
2734         beq     10f
2735         bset    IMM (31),d1
2736 10:     lsrl    IMM (1),d0
2737 #endif
2738         addl    IMM (1),d2
2739 1:
2740         lea     Laddsf$4,a0     | to return from rounding routine
2741         lea     SYM (_fpCCR),a1 | check the rounding mode
2742 #ifdef __mcoldfire__
2743         clrl    d6
2744 #endif
2745         movew   a1@(6),d6       | rounding mode in d6
2746         beq     Lround$to$nearest
2747 #ifndef __mcoldfire__
2748         cmpw    IMM (ROUND_TO_PLUS),d6
2749 #else
2750         cmpl    IMM (ROUND_TO_PLUS),d6
2751 #endif
2752         bhi     Lround$to$minus
2753         blt     Lround$to$zero
2754         bra     Lround$to$plus
2755 Laddsf$4:
2756 | Put back the exponent, but check for overflow.
2757 #ifndef __mcoldfire__
2758         cmpw    IMM (0xff),d2
2759 #else
2760         cmpl    IMM (0xff),d2
2761 #endif
2762         bhi     1f
2763         bclr    IMM (FLT_MANT_DIG-1),d0
2764 #ifndef __mcoldfire__
2765         lslw    IMM (7),d2
2766 #else
2767         lsll    IMM (7),d2
2768 #endif
2769         swap    d2
2770         orl     d2,d0
2771         bra     Laddsf$ret
2772 1:
2773         movew   IMM (ADD),d5
2774         bra     Lf$overflow
2775
2776 Lsubsf$0:
2777 | We are here if a > 0 and b < 0 (sign bits cleared).
2778 | Here we do the subtraction.
2779         movel   d6,d7           | put sign in d7
2780         andl    IMM (0x80000000),d7
2781
2782         subl    d3,d1           | result in d0-d1
2783         subxl   d2,d0           |
2784         beq     Laddsf$ret      | if zero just exit
2785         bpl     1f              | if positive skip the following
2786         bchg    IMM (31),d7     | change sign bit in d7
2787         negl    d1
2788         negxl   d0
2789 1:
2790 #ifndef __mcoldfire__
2791         exg     d2,a0           | now we have the exponent in d2
2792         lsrw    IMM (8),d2      | put it in the first byte
2793 #else
2794         movel   d2,d4
2795         movel   a0,d2
2796         movel   d4,a0
2797         lsrl    IMM (8),d2      | put it in the first byte
2798 #endif
2799
2800 | Now d0-d1 is positive and the sign bit is in d7.
2801
2802 | Note that we do not have to normalize, since in the subtraction bit
2803 | #FLT_MANT_DIG+1 is never set, and denormalized numbers are handled by
2804 | the rounding routines themselves.
2805         lea     Lsubsf$1,a0     | to return from rounding routine
2806         lea     SYM (_fpCCR),a1 | check the rounding mode
2807 #ifdef __mcoldfire__
2808         clrl    d6
2809 #endif
2810         movew   a1@(6),d6       | rounding mode in d6
2811         beq     Lround$to$nearest
2812 #ifndef __mcoldfire__
2813         cmpw    IMM (ROUND_TO_PLUS),d6
2814 #else
2815         cmpl    IMM (ROUND_TO_PLUS),d6
2816 #endif
2817         bhi     Lround$to$minus
2818         blt     Lround$to$zero
2819         bra     Lround$to$plus
2820 Lsubsf$1:
2821 | Put back the exponent (we can't have overflow!). '
2822         bclr    IMM (FLT_MANT_DIG-1),d0
2823 #ifndef __mcoldfire__
2824         lslw    IMM (7),d2
2825 #else
2826         lsll    IMM (7),d2
2827 #endif
2828         swap    d2
2829         orl     d2,d0
2830         bra     Laddsf$ret
2831
2832 | If one of the numbers was too small (difference of exponents >=
2833 | FLT_MANT_DIG+2) we return the other (and now we don't have to '
2834 | check for finiteness or zero).
2835 Laddsf$a$small:
2836         movel   a6@(12),d0
2837         lea     SYM (_fpCCR),a0
2838         movew   IMM (0),a0@
2839 #ifndef __mcoldfire__
2840         moveml  sp@+,d2-d7      | restore data registers
2841 #else
2842         moveml  sp@,d2-d7
2843         | XXX if frame pointer is ever removed, stack pointer must
2844         | be adjusted here.
2845 #endif
2846         unlk    a6              | and return
2847         rts
2848
2849 Laddsf$b$small:
2850         movel   a6@(8),d0
2851         lea     SYM (_fpCCR),a0
2852         movew   IMM (0),a0@
2853 #ifndef __mcoldfire__
2854         moveml  sp@+,d2-d7      | restore data registers
2855 #else
2856         moveml  sp@,d2-d7
2857         | XXX if frame pointer is ever removed, stack pointer must
2858         | be adjusted here.
2859 #endif
2860         unlk    a6              | and return
2861         rts
2862
2863 | If the numbers are denormalized remember to put exponent equal to 1.
2864
2865 Laddsf$a$den:
2866         movel   d5,d6           | d5 contains 0x01000000
2867         swap    d6
2868         bra     Laddsf$1
2869
2870 Laddsf$b$den:
2871         movel   d5,d7
2872         swap    d7
2873         notl    d4              | make d4 into a mask for the fraction
2874                                 | (this was not executed after the jump)
2875         bra     Laddsf$2
2876
2877 | The rest is mainly code for the different results which can be
2878 | returned (checking always for +/-INFINITY and NaN).
2879
2880 Laddsf$b:
2881 | Return b (if a is zero).
2882         movel   a6@(12),d0
2883         bra     1f
2884 Laddsf$a:
2885 | Return a (if b is zero).
2886         movel   a6@(8),d0
2887 1:
2888         movew   IMM (ADD),d5
2889 | We have to check for NaN and +/-infty.
2890         movel   d0,d7
2891         andl    IMM (0x80000000),d7     | put sign in d7
2892         bclr    IMM (31),d0             | clear sign
2893         cmpl    IMM (INFINITY),d0       | check for infty or NaN
2894         bge     2f
2895         movel   d0,d0           | check for zero (we do this because we don't '
2896         bne     Laddsf$ret      | want to return -0 by mistake
2897         bclr    IMM (31),d7     | if zero be sure to clear sign
2898         bra     Laddsf$ret      | if everything OK just return
2899 2:
2900 | The value to be returned is either +/-infty or NaN
2901         andl    IMM (0x007fffff),d0     | check for NaN
2902         bne     Lf$inop                 | if mantissa not zero is NaN
2903         bra     Lf$infty
2904
2905 Laddsf$ret:
2906 | Normal exit (a and b nonzero, result is not NaN nor +/-infty).
2907 | We have to clear the exception flags (just the exception type).
2908         lea     SYM (_fpCCR),a0
2909         movew   IMM (0),a0@
2910         orl     d7,d0           | put sign bit
2911 #ifndef __mcoldfire__
2912         moveml  sp@+,d2-d7      | restore data registers
2913 #else
2914         moveml  sp@,d2-d7
2915         | XXX if frame pointer is ever removed, stack pointer must
2916         | be adjusted here.
2917 #endif
2918         unlk    a6              | and return
2919         rts
2920
2921 Laddsf$ret$den:
2922 | Return a denormalized number (for addition we don't signal underflow) '
2923         lsrl    IMM (1),d0      | remember to shift right back once
2924         bra     Laddsf$ret      | and return
2925
2926 | Note: when adding two floats of the same sign if either one is
2927 | NaN we return NaN without regard to whether the other is finite or
2928 | not. When subtracting them (i.e., when adding two numbers of
2929 | opposite signs) things are more complicated: if both are INFINITY
2930 | we return NaN, if only one is INFINITY and the other is NaN we return
2931 | NaN, but if it is finite we return INFINITY with the corresponding sign.
2932
2933 Laddsf$nf:
2934         movew   IMM (ADD),d5
2935 | This could be faster but it is not worth the effort, since it is not
2936 | executed very often. We sacrifice speed for clarity here.
2937         movel   a6@(8),d0       | get the numbers back (remember that we
2938         movel   a6@(12),d1      | did some processing already)
2939         movel   IMM (INFINITY),d4 | useful constant (INFINITY)
2940         movel   d0,d2           | save sign bits
2941         movel   d1,d3
2942         bclr    IMM (31),d0     | clear sign bits
2943         bclr    IMM (31),d1
2944 | We know that one of them is either NaN of +/-INFINITY
2945 | Check for NaN (if either one is NaN return NaN)
2946         cmpl    d4,d0           | check first a (d0)
2947         bhi     Lf$inop
2948         cmpl    d4,d1           | check now b (d1)
2949         bhi     Lf$inop
2950 | Now comes the check for +/-INFINITY. We know that both are (maybe not
2951 | finite) numbers, but we have to check if both are infinite whether we
2952 | are adding or subtracting them.
2953         eorl    d3,d2           | to check sign bits
2954         bmi     1f
2955         movel   d0,d7
2956         andl    IMM (0x80000000),d7     | get (common) sign bit
2957         bra     Lf$infty
2958 1:
2959 | We know one (or both) are infinite, so we test for equality between the
2960 | two numbers (if they are equal they have to be infinite both, so we
2961 | return NaN).
2962         cmpl    d1,d0           | are both infinite?
2963         beq     Lf$inop         | if so return NaN
2964
2965         movel   d0,d7
2966         andl    IMM (0x80000000),d7 | get a's sign bit '
2967         cmpl    d4,d0           | test now for infinity
2968         beq     Lf$infty        | if a is INFINITY return with this sign
2969         bchg    IMM (31),d7     | else we know b is INFINITY and has
2970         bra     Lf$infty        | the opposite sign
2971
2972 |=============================================================================
2973 |                             __mulsf3
2974 |=============================================================================
2975
2976 | float __mulsf3(float, float);
2977 SYM (__mulsf3):
2978 #ifndef __mcoldfire__
2979         link    a6,IMM (0)
2980         moveml  d2-d7,sp@-
2981 #else
2982         link    a6,IMM (-24)
2983         moveml  d2-d7,sp@
2984 #endif
2985         movel   a6@(8),d0       | get a into d0
2986         movel   a6@(12),d1      | and b into d1
2987         movel   d0,d7           | d7 will hold the sign of the product
2988         eorl    d1,d7           |
2989         andl    IMM (0x80000000),d7
2990         movel   IMM (INFINITY),d6       | useful constant (+INFINITY)
2991         movel   d6,d5                   | another (mask for fraction)
2992         notl    d5                      |
2993         movel   IMM (0x00800000),d4     | this is to put hidden bit back
2994         bclr    IMM (31),d0             | get rid of a's sign bit '
2995         movel   d0,d2                   |
2996         beq     Lmulsf$a$0              | branch if a is zero
2997         bclr    IMM (31),d1             | get rid of b's sign bit '
2998         movel   d1,d3           |
2999         beq     Lmulsf$b$0      | branch if b is zero
3000         cmpl    d6,d0           | is a big?
3001         bhi     Lmulsf$inop     | if a is NaN return NaN
3002         beq     Lmulsf$inf      | if a is INFINITY we have to check b
3003         cmpl    d6,d1           | now compare b with INFINITY
3004         bhi     Lmulsf$inop     | is b NaN?
3005         beq     Lmulsf$overflow | is b INFINITY?
3006 | Here we have both numbers finite and nonzero (and with no sign bit).
3007 | Now we get the exponents into d2 and d3.
3008         andl    d6,d2           | and isolate exponent in d2
3009         beq     Lmulsf$a$den    | if exponent is zero we have a denormalized
3010         andl    d5,d0           | and isolate fraction
3011         orl     d4,d0           | and put hidden bit back
3012         swap    d2              | I like exponents in the first byte
3013 #ifndef __mcoldfire__
3014         lsrw    IMM (7),d2      |
3015 #else
3016         lsrl    IMM (7),d2      |
3017 #endif
3018 Lmulsf$1:                       | number
3019         andl    d6,d3           |
3020         beq     Lmulsf$b$den    |
3021         andl    d5,d1           |
3022         orl     d4,d1           |
3023         swap    d3              |
3024 #ifndef __mcoldfire__
3025         lsrw    IMM (7),d3      |
3026 #else
3027         lsrl    IMM (7),d3      |
3028 #endif
3029 Lmulsf$2:                       |
3030 #ifndef __mcoldfire__
3031         addw    d3,d2           | add exponents
3032         subw    IMM (F_BIAS+1),d2 | and subtract bias (plus one)
3033 #else
3034         addl    d3,d2           | add exponents
3035         subl    IMM (F_BIAS+1),d2 | and subtract bias (plus one)
3036 #endif
3037
3038 | We are now ready to do the multiplication. The situation is as follows:
3039 | both a and b have bit FLT_MANT_DIG-1 set (even if they were
3040 | denormalized to start with!), which means that in the product
3041 | bit 2*(FLT_MANT_DIG-1) (that is, bit 2*FLT_MANT_DIG-2-32 of the
3042 | high long) is set.
3043
3044 | To do the multiplication let us move the number a little bit around ...
3045         movel   d1,d6           | second operand in d6
3046         movel   d0,d5           | first operand in d4-d5
3047         movel   IMM (0),d4
3048         movel   d4,d1           | the sums will go in d0-d1
3049         movel   d4,d0
3050
3051 | now bit FLT_MANT_DIG-1 becomes bit 31:
3052         lsll    IMM (31-FLT_MANT_DIG+1),d6
3053
3054 | Start the loop (we loop #FLT_MANT_DIG times):
3055         movew   IMM (FLT_MANT_DIG-1),d3
3056 1:      addl    d1,d1           | shift sum
3057         addxl   d0,d0
3058         lsll    IMM (1),d6      | get bit bn
3059         bcc     2f              | if not set skip sum
3060         addl    d5,d1           | add a
3061         addxl   d4,d0
3062 2:
3063 #ifndef __mcoldfire__
3064         dbf     d3,1b           | loop back
3065 #else
3066         subql   IMM (1),d3
3067         bpl     1b
3068 #endif
3069
3070 | Now we have the product in d0-d1, with bit (FLT_MANT_DIG - 1) + FLT_MANT_DIG
3071 | (mod 32) of d0 set. The first thing to do now is to normalize it so bit
3072 | FLT_MANT_DIG is set (to do the rounding).
3073 #ifndef __mcoldfire__
3074         rorl    IMM (6),d1
3075         swap    d1
3076         movew   d1,d3
3077         andw    IMM (0x03ff),d3
3078         andw    IMM (0xfd00),d1
3079 #else
3080         movel   d1,d3
3081         lsll    IMM (8),d1
3082         addl    d1,d1
3083         addl    d1,d1
3084         moveq   IMM (22),d5
3085         lsrl    d5,d3
3086         orl     d3,d1
3087         andl    IMM (0xfffffd00),d1
3088 #endif
3089         lsll    IMM (8),d0
3090         addl    d0,d0
3091         addl    d0,d0
3092 #ifndef __mcoldfire__
3093         orw     d3,d0
3094 #else
3095         orl     d3,d0
3096 #endif
3097
3098         movew   IMM (MULTIPLY),d5
3099
3100         btst    IMM (FLT_MANT_DIG+1),d0
3101         beq     Lround$exit
3102 #ifndef __mcoldfire__
3103         lsrl    IMM (1),d0
3104         roxrl   IMM (1),d1
3105         addw    IMM (1),d2
3106 #else
3107         lsrl    IMM (1),d1
3108         btst    IMM (0),d0
3109         beq     10f
3110         bset    IMM (31),d1
3111 10:     lsrl    IMM (1),d0
3112         addql   IMM (1),d2
3113 #endif
3114         bra     Lround$exit
3115
3116 Lmulsf$inop:
3117         movew   IMM (MULTIPLY),d5
3118         bra     Lf$inop
3119
3120 Lmulsf$overflow:
3121         movew   IMM (MULTIPLY),d5
3122         bra     Lf$overflow
3123
3124 Lmulsf$inf:
3125         movew   IMM (MULTIPLY),d5
3126 | If either is NaN return NaN; else both are (maybe infinite) numbers, so
3127 | return INFINITY with the correct sign (which is in d7).
3128         cmpl    d6,d1           | is b NaN?
3129         bhi     Lf$inop         | if so return NaN
3130         bra     Lf$overflow     | else return +/-INFINITY
3131
3132 | If either number is zero return zero, unless the other is +/-INFINITY,
3133 | or NaN, in which case we return NaN.
3134 Lmulsf$b$0:
3135 | Here d1 (==b) is zero.
3136         movel   d1,d0           | put b into d0 (just a zero)
3137         movel   a6@(8),d1       | get a again to check for non-finiteness
3138         bra     1f
3139 Lmulsf$a$0:
3140         movel   a6@(12),d1      | get b again to check for non-finiteness
3141 1:      bclr    IMM (31),d1     | clear sign bit
3142         cmpl    IMM (INFINITY),d1 | and check for a large exponent
3143         bge     Lf$inop         | if b is +/-INFINITY or NaN return NaN
3144         lea     SYM (_fpCCR),a0 | else return zero
3145         movew   IMM (0),a0@     |
3146 #ifndef __mcoldfire__
3147         moveml  sp@+,d2-d7      |
3148 #else
3149         moveml  sp@,d2-d7
3150         | XXX if frame pointer is ever removed, stack pointer must
3151         | be adjusted here.
3152 #endif
3153         unlk    a6              |
3154         rts                     |
3155
3156 | If a number is denormalized we put an exponent of 1 but do not put the
3157 | hidden bit back into the fraction; instead we shift left until bit 23
3158 | (the hidden bit) is set, adjusting the exponent accordingly. We do this
3159 | to ensure that the product of the fractions is close to 1.
3160 Lmulsf$a$den:
3161         movel   IMM (1),d2
3162         andl    d5,d0
3163 1:      addl    d0,d0           | shift a left (until bit 23 is set)
3164 #ifndef __mcoldfire__
3165         subw    IMM (1),d2      | and adjust exponent
3166 #else
3167         subql   IMM (1),d2      | and adjust exponent
3168 #endif
3169         btst    IMM (FLT_MANT_DIG-1),d0
3170         bne     Lmulsf$1        |
3171         bra     1b              | else loop back
3172
3173 Lmulsf$b$den:
3174         movel   IMM (1),d3
3175         andl    d5,d1
3176 1:      addl    d1,d1           | shift b left until bit 23 is set
3177 #ifndef __mcoldfire__
3178         subw    IMM (1),d3      | and adjust exponent
3179 #else
3180         subl    IMM (1),d3      | and adjust exponent
3181 #endif
3182         btst    IMM (FLT_MANT_DIG-1),d1
3183         bne     Lmulsf$2        |
3184         bra     1b              | else loop back
3185
3186 |=============================================================================
3187 |                             __divsf3
3188 |=============================================================================
3189
3190 | float __divsf3(float, float);
3191 SYM (__divsf3):
3192 #ifndef __mcoldfire__
3193         link    a6,IMM (0)
3194         moveml  d2-d7,sp@-
3195 #else
3196         link    a6,IMM (-24)
3197         moveml  d2-d7,sp@
3198 #endif
3199         movel   a6@(8),d0               | get a into d0
3200         movel   a6@(12),d1              | and b into d1
3201         movel   d0,d7                   | d7 will hold the sign of the result
3202         eorl    d1,d7                   |
3203         andl    IMM (0x80000000),d7     |
3204         movel   IMM (INFINITY),d6       | useful constant (+INFINITY)
3205         movel   d6,d5                   | another (mask for fraction)
3206         notl    d5                      |
3207         movel   IMM (0x00800000),d4     | this is to put hidden bit back
3208         bclr    IMM (31),d0             | get rid of a's sign bit '
3209         movel   d0,d2                   |
3210         beq     Ldivsf$a$0              | branch if a is zero
3211         bclr    IMM (31),d1             | get rid of b's sign bit '
3212         movel   d1,d3                   |
3213         beq     Ldivsf$b$0              | branch if b is zero
3214         cmpl    d6,d0                   | is a big?
3215         bhi     Ldivsf$inop             | if a is NaN return NaN
3216         beq     Ldivsf$inf              | if a is INFINITY we have to check b
3217         cmpl    d6,d1                   | now compare b with INFINITY
3218         bhi     Ldivsf$inop             | if b is NaN return NaN
3219         beq     Ldivsf$underflow
3220 | Here we have both numbers finite and nonzero (and with no sign bit).
3221 | Now we get the exponents into d2 and d3 and normalize the numbers to
3222 | ensure that the ratio of the fractions is close to 1. We do this by
3223 | making sure that bit #FLT_MANT_DIG-1 (hidden bit) is set.
3224         andl    d6,d2           | and isolate exponent in d2
3225         beq     Ldivsf$a$den    | if exponent is zero we have a denormalized
3226         andl    d5,d0           | and isolate fraction
3227         orl     d4,d0           | and put hidden bit back
3228         swap    d2              | I like exponents in the first byte