gcc/config/m68k/lb1sf68.asm

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