git.gag.com Git - fw/sdcc/blob - device/lib/_fsdiv.c

   1 /* Floating point library in optimized assembly for 8051
   2  * Copyright (c) 2004, Paul Stoffregen, paul@pjrc.com
   3  *
   4  * This program is free software; you can redistribute it and/or
   5  * modify it under the terms of the GNU Library General Public License
   6  * as published by the Free Software Foundation; either version 2
   7  * of the License, or (at your option) any later version.
   8  *
   9  * This 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
  12  * GNU General Public License for more details.
  13  *
  14  * You should have received a copy of the GNU General Public License
  15  * along with this program; if not, write to the Free Software
  16  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
  17  */
  18
  19
  20 #define SDCC_FLOAT_LIB
  21 #include <float.h>
  22
  23
  24 #ifdef FLOAT_ASM_MCS51
  25
  26 // float __fsdiv (float a, float b) reentrant
  27 static void dummy(void) _naked
  28 {
  29         _asm
  30         .globl  ___fsdiv
  31 ___fsdiv:
  32         // extract the two inputs, placing them into:
  33         //      sign     exponent   mantiassa
  34         //      ----     --------   ---------
  35         //  a:  sign_a   exp_a      r4/r3/r2
  36         //  b:  sign_b   exp_b      r7/r6/r5
  37
  38         lcall   fsgetargs
  39
  40         // compute final sign bit
  41         jnb     sign_b, 00001$
  42         cpl     sign_a
  43 00001$:
  44
  45         // if divisor is zero, return infinity
  46         cjne    r7, #0, 00002$
  47         ljmp    fs_return_inf
  48 00002$:
  49         // if dividend is zero, return zero
  50         cjne    r4, #0, 00003$
  51         ljmp    fs_return_zero
  52 00003$:
  53
  54         // subtract exponents
  55         clr     c
  56         mov     a, exp_a
  57         subb    a, exp_b
  58         add     a, #127
  59         mov     exp_a, a
  60
  61         // need extra bits on a's mantissa
  62 #ifdef FLOAT_FULL_ACCURACY
  63         clr     c
  64         mov     a, r5
  65         subb    a, r2
  66         mov     a, r6
  67         subb    a, r3
  68         mov     a, r7
  69         subb    a, r4
  70         jc      00005$
  71         dec     exp_a
  72         clr     c
  73         mov     a, r2
  74         rlc     a
  75         mov     r1, a
  76         mov     a, r3
  77         rlc     a
  78         mov     r2, a
  79         mov     a, r4
  80         rlc     a
  81         mov     r3, a
  82         clr     a
  83         rlc     a
  84         mov     r4, a
  85         sjmp    00006$
  86 00005$:
  87 #endif
  88         clr     a
  89         xch     a, r4
  90         xch     a, r3
  91         xch     a, r2
  92         mov     r1, a
  93 00006$:
  94
  95         // begin long division
  96         push    exp_a
  97 #ifdef FLOAT_FULL_ACCURACY
  98         mov     b, #25
  99 #else
 100         mov     b, #24
 101 #endif
 102 00010$:
 103         // compare
 104         clr     c
 105         mov     a, r1
 106         subb    a, r5
 107         mov     a, r2
 108         subb    a, r6
 109         mov     a, r3
 110         subb    a, r7
 111         mov     a, r4
 112         subb    a, #0           // carry==0 if mant1 >= mant2
 113
 114 #ifdef FLOAT_FULL_ACCURACY
 115         djnz    b, 00011$
 116         sjmp    00015$
 117 00011$:
 118 #endif
 119         jc      00012$
 120         // subtract
 121         mov     a, r1
 122         subb    a, r5
 123         mov     r1, a
 124         mov     a, r2
 125         subb    a, r6
 126         mov     r2, a
 127         mov     a, r3
 128         subb    a, r7
 129         mov     r3, a
 130         mov     a, r4
 131         subb    a, #0
 132         mov     r4, a
 133         clr     c
 134
 135 00012$:
 136         // shift result
 137         cpl     c
 138         mov     a, r0
 139         rlc     a
 140         mov     r0, a
 141         mov     a, dpl
 142         rlc     a
 143         mov     dpl, a
 144         mov     a, dph
 145         rlc     a
 146         mov     dph, a
 147
 148         // shift partial remainder
 149         clr     c
 150         mov     a, r1
 151         rlc     a
 152         mov     r1, a
 153         mov     a, r2
 154         rlc     a
 155         mov     r2, a
 156         mov     a, r3
 157         rlc     a
 158         mov     r3, a
 159         mov     a, r4
 160         rlc     a
 161         mov     r4, a
 162
 163 #ifdef FLOAT_FULL_ACCURACY
 164         sjmp    00010$
 165 00015$:
 166 #else
 167         djnz    b, 00010$
 168 #endif
 169
 170         // now we've got a division result, so all we need to do
 171         // is round off properly, normalize and output a float
 172
 173 #ifdef FLOAT_FULL_ACCURACY
 174         cpl     c
 175         clr     a
 176         mov     r1, a
 177         addc    a, r0
 178         mov     r2, a
 179         clr     a
 180         addc    a, dpl
 181         mov     r3, a
 182         clr     a
 183         addc    a, dph
 184         mov     r4, a
 185         pop     exp_a
 186         jnc     00016$
 187         inc     exp_a
 188         mov     r4, #0x80
 189 00016$:
 190 #else
 191         mov     r1, #0
 192         mov     a, r0
 193         mov     r2, a
 194         mov     r3, dpl
 195         mov     r4, dph
 196         pop     exp_a
 197 #endif
 198
 199         lcall   fs_normalize_a
 200         ljmp    fs_zerocheck_return
 201         _endasm;
 202 }
 203
 204 #else
 205
 206
 207
 208
 209 /*
 210 ** libgcc support for software floating point.
 211 ** Copyright (C) 1991 by Pipeline Associates, Inc.  All rights reserved.
 212 ** Permission is granted to do *anything* you want with this file,
 213 ** commercial or otherwise, provided this message remains intact.  So there!
 214 ** I would appreciate receiving any updates/patches/changes that anyone
 215 ** makes, and am willing to be the repository for said changes (am I
 216 ** making a big mistake?).
 217 **
 218 ** Pat Wood
 219 ** Pipeline Associates, Inc.
 220 ** pipeline!phw@motown.com or
 221 ** sun!pipeline!phw or
 222 ** uunet!motown!pipeline!phw
 223 */
 224
 225 /* (c)2000/2001: hacked a little by johan.knol@iduna.nl for sdcc */
 226
 227
 228 union float_long
 229   {
 230     float f;
 231     long l;
 232   };
 233
 234 /* divide two floats */
 235 float __fsdiv (float a1, float a2)
 236 {
 237   volatile union float_long fl1, fl2;
 238   volatile long result;
 239   volatile unsigned long mask;
 240   volatile long mant1, mant2;
 241   volatile int exp ;
 242   char sign;
 243
 244   fl1.f = a1;
 245   fl2.f = a2;
 246
 247   /* subtract exponents */
 248   exp = EXP (fl1.l) ;
 249   exp -= EXP (fl2.l);
 250   exp += EXCESS;
 251
 252   /* compute sign */
 253   sign = SIGN (fl1.l) ^ SIGN (fl2.l);
 254
 255   /* divide by zero??? */
 256   if (!fl2.l)
 257     /* return NaN or -NaN */
 258     return (-1.0);
 259
 260   /* numerator zero??? */
 261   if (!fl1.l)
 262     return (0);
 263
 264   /* now get mantissas */
 265   mant1 = MANT (fl1.l);
 266   mant2 = MANT (fl2.l);
 267
 268   /* this assures we have 25 bits of precision in the end */
 269   if (mant1 < mant2)
 270     {
 271       mant1 <<= 1;
 272       exp--;
 273     }
 274
 275   /* now we perform repeated subtraction of fl2.l from fl1.l */
 276   mask = 0x1000000;
 277   result = 0;
 278   while (mask)
 279     {
 280       if (mant1 >= mant2)
 281         {
 282           result |= mask;
 283           mant1 -= mant2;
 284         }
 285       mant1 <<= 1;
 286       mask >>= 1;
 287     }
 288
 289   /* round */
 290   result += 1;
 291
 292   /* normalize down */
 293   exp++;
 294   result >>= 1;
 295
 296   result &= ~HIDDEN;
 297
 298   /* pack up and go home */
 299   fl1.l = PACK (sign ? 1ul<<31 : 0, (unsigned long) exp, result);
 300   return (fl1.f);
 301 }
 302
 303 #endif
 304