sysdeps/alpha/divqu.S

   1 /* Copyright (C) 2004-2024 Free Software Foundation, Inc.
   2    This file is part of the GNU C Library.
   3
   4    The GNU C Library is free software; you can redistribute it and/or
   5    modify it under the terms of the GNU Lesser General Public
   6    License as published by the Free Software Foundation; either
   7    version 2.1 of the License, or (at your option) any later version.
   8
   9    The GNU C Library is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12    Lesser General Public License for more details.
  13
  14    You should have received a copy of the GNU Lesser General Public
  15    License along with the GNU C Library.  If not, see
  16    <https://www.gnu.org/licenses/>.  */
  17
  18 #include "div_libc.h"
  19
  20
  21 /* 64-bit unsigned long divide.  These are not normal C functions.  Argument
  22    registers are t10 and t11, the result goes in t12.  Only t12 and AT may be
  23    clobbered.
  24
  25    Theory of operation here is that we can use the FPU divider for virtually
  26    all operands that we see: all dividend values between -2**53 and 2**53-1
  27    can be computed directly.  Note that divisor values need not be checked
  28    against that range because the rounded fp value will be close enough such
  29    that the quotient is < 1, which will properly be truncated to zero when we
  30    convert back to integer.
  31
  32    When the dividend is outside the range for which we can compute exact
  33    results, we use the fp quotient as an estimate from which we begin refining
  34    an exact integral value.  This reduces the number of iterations in the
  35    shift-and-subtract loop significantly.
  36
  37    The FPCR save/restore is due to the fact that the EV6 _will_ set FPCR_INE
  38    for cvttq/c even without /sui being set.  It will not, however, properly
  39    raise the exception, so we don't have to worry about FPCR_INED being clear
  40    and so dying by SIGFPE.  */
  41
  42         .text
  43         .align  4
  44         .globl  __divqu
  45         .type   __divqu, @funcnoplt
  46         .usepv  __divqu, no
  47
  48         cfi_startproc
  49         cfi_return_column (RA)
  50 __divqu:
  51         lda     sp, -FRAME(sp)
  52         cfi_def_cfa_offset (FRAME)
  53         CALL_MCOUNT
  54
  55         /* Get the fp divide insn issued as quickly as possible.  After
  56            that's done, we have at least 22 cycles until its results are
  57            ready -- all the time in the world to figure out how we're
  58            going to use the results.  */
  59         beq     Y, DIVBYZERO
  60
  61         stt     $f0, 0(sp)
  62         excb
  63         stt     $f1, 8(sp)
  64         stt     $f3, 48(sp)
  65         cfi_rel_offset ($f0, 0)
  66         cfi_rel_offset ($f1, 8)
  67         cfi_rel_offset ($f3, 48)
  68         mf_fpcr $f3
  69
  70         _ITOFT2 X, $f0, 16, Y, $f1, 24
  71         cvtqt   $f0, $f0
  72         cvtqt   $f1, $f1
  73
  74         blt     X, $x_is_neg
  75         divt/c  $f0, $f1, $f0
  76
  77         /* Check to see if Y was mis-converted as signed value.  */
  78         ldt     $f1, 8(sp)
  79         blt     Y, $y_is_neg
  80
  81         /* Check to see if X fit in the double as an exact value.  */
  82         srl     X, 53, AT
  83         bne     AT, $x_big
  84
  85         /* If we get here, we're expecting exact results from the division.
  86            Do nothing else besides convert and clean up.  */
  87         cvttq/c $f0, $f0
  88         excb
  89         mt_fpcr $f3
  90         _FTOIT  $f0, RV, 16
  91
  92         ldt     $f0, 0(sp)
  93         ldt     $f3, 48(sp)
  94         lda     sp, FRAME(sp)
  95         cfi_remember_state
  96         cfi_restore ($f0)
  97         cfi_restore ($f1)
  98         cfi_restore ($f3)
  99         cfi_def_cfa_offset (0)
 100         ret     $31, (RA), 1
 101
 102         .align  4
 103         cfi_restore_state
 104 $x_is_neg:
 105         /* If we get here, X is so big that bit 63 is set, which made the
 106            conversion come out negative.  Fix it up lest we not even get
 107            a good estimate.  */
 108         ldah    AT, 0x5f80              /* 2**64 as float.  */
 109         stt     $f2, 24(sp)
 110         cfi_rel_offset ($f2, 24)
 111         _ITOFS  AT, $f2, 16
 112
 113         .align  4
 114         addt    $f0, $f2, $f0
 115         unop
 116         divt/c  $f0, $f1, $f0
 117         unop
 118
 119         /* Ok, we've now the divide issued.  Continue with other checks.  */
 120         ldt     $f1, 8(sp)
 121         unop
 122         ldt     $f2, 24(sp)
 123         blt     Y, $y_is_neg
 124         cfi_restore ($f1)
 125         cfi_restore ($f2)
 126         cfi_remember_state      /* for y_is_neg */
 127
 128         .align  4
 129 $x_big:
 130         /* If we get here, X is large enough that we don't expect exact
 131            results, and neither X nor Y got mis-translated for the fp
 132            division.  Our task is to take the fp result, figure out how
 133            far it's off from the correct result and compute a fixup.  */
 134         stq     t0, 16(sp)
 135         stq     t1, 24(sp)
 136         stq     t2, 32(sp)
 137         stq     t3, 40(sp)
 138         cfi_rel_offset (t0, 16)
 139         cfi_rel_offset (t1, 24)
 140         cfi_rel_offset (t2, 32)
 141         cfi_rel_offset (t3, 40)
 142
 143 #define Q       RV              /* quotient */
 144 #define R       t0              /* remainder */
 145 #define SY      t1              /* scaled Y */
 146 #define S       t2              /* scalar */
 147 #define QY      t3              /* Q*Y */
 148
 149         cvttq/c $f0, $f0
 150         _FTOIT  $f0, Q, 8
 151         mulq    Q, Y, QY
 152
 153         .align  4
 154         stq     t4, 8(sp)
 155         excb
 156         ldt     $f0, 0(sp)
 157         mt_fpcr $f3
 158         cfi_rel_offset (t4, 8)
 159         cfi_restore ($f0)
 160
 161         subq    QY, X, R
 162         mov     Y, SY
 163         mov     1, S
 164         bgt     R, $q_high
 165
 166 $q_high_ret:
 167         subq    X, QY, R
 168         mov     Y, SY
 169         mov     1, S
 170         bgt     R, $q_low
 171
 172 $q_low_ret:
 173         ldq     t4, 8(sp)
 174         ldq     t0, 16(sp)
 175         ldq     t1, 24(sp)
 176         ldq     t2, 32(sp)
 177
 178         ldq     t3, 40(sp)
 179         ldt     $f3, 48(sp)
 180         lda     sp, FRAME(sp)
 181         cfi_remember_state
 182         cfi_restore (t0)
 183         cfi_restore (t1)
 184         cfi_restore (t2)
 185         cfi_restore (t3)
 186         cfi_restore (t4)
 187         cfi_restore ($f3)
 188         cfi_def_cfa_offset (0)
 189         ret     $31, (RA), 1
 190
 191         .align  4
 192         cfi_restore_state
 193         /* The quotient that we computed was too large.  We need to reduce
 194            it by S such that Y*S >= R.  Obviously the closer we get to the
 195            correct value the better, but overshooting high is ok, as we'll
 196            fix that up later.  */
 197 0:
 198         addq    SY, SY, SY
 199         addq    S, S, S
 200 $q_high:
 201         cmpult  SY, R, AT
 202         bne     AT, 0b
 203
 204         subq    Q, S, Q
 205         unop
 206         subq    QY, SY, QY
 207         br      $q_high_ret
 208
 209         .align  4
 210         /* The quotient that we computed was too small.  Divide Y by the
 211            current remainder (R) and add that to the existing quotient (Q).
 212            The expectation, of course, is that R is much smaller than X.  */
 213         /* Begin with a shift-up loop.  Compute S such that Y*S >= R.  We
 214            already have a copy of Y in SY and the value 1 in S.  */
 215 0:
 216         addq    SY, SY, SY
 217         addq    S, S, S
 218 $q_low:
 219         cmpult  SY, R, AT
 220         bne     AT, 0b
 221
 222         /* Shift-down and subtract loop.  Each iteration compares our scaled
 223            Y (SY) with the remainder (R); if SY <= R then X is divisible by
 224            Y's scalar (S) so add it to the quotient (Q).  */
 225 2:      addq    Q, S, t3
 226         srl     S, 1, S
 227         cmpule  SY, R, AT
 228         subq    R, SY, t4
 229
 230         cmovne  AT, t3, Q
 231         cmovne  AT, t4, R
 232         srl     SY, 1, SY
 233         bne     S, 2b
 234
 235         br      $q_low_ret
 236
 237         .align  4
 238         cfi_restore_state
 239 $y_is_neg:
 240         /* If we get here, Y is so big that bit 63 is set.  The results
 241            from the divide will be completely wrong.  Fortunately, the
 242            quotient must be either 0 or 1, so just compute it directly.  */
 243         cmpule  Y, X, RV
 244         excb
 245         mt_fpcr $f3
 246         ldt     $f0, 0(sp)
 247         ldt     $f3, 48(sp)
 248         lda     sp, FRAME(sp)
 249         cfi_restore ($f0)
 250         cfi_restore ($f3)
 251         cfi_def_cfa_offset (0)
 252         ret     $31, (RA), 1
 253
 254         cfi_endproc
 255         .size   __divqu, .-__divqu
 256
 257         DO_DIVBYZERO