git.gag.com Git - debian/cc1111/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, ...
  46         cjne    r7, #0, 00003$
  47         // if dividend is also zero, return NaN
  48         cjne    r4, #0, 00002$
  49         ljmp    fs_return_nan
  50 00002$:
  51         // but dividend is non-zero, return infinity
  52         ljmp    fs_return_inf
  53 00003$:
  54         // if dividend is zero, return zero
  55         cjne    r4, #0, 00004$
  56         ljmp    fs_return_zero
  57 00004$:
  58         // if divisor is infinity, ...
  59         mov     a, exp_b
  60         cjne    a, #0xFF, 00006$
  61         // and dividend is also infinity, return NaN
  62         mov     a, exp_a
  63         cjne    a, #0xFF, 00005$
  64         ljmp    fs_return_nan
  65 00005$:
  66         // but dividend is not infinity, return zero
  67         ljmp    fs_return_zero
  68 00006$:
  69         // if dividend is infinity, return infinity
  70         mov     a, exp_a
  71         cjne    a, #0xFF, 00007$
  72         ljmp    fs_return_inf
  73 00007$:
  74
  75         // subtract exponents
  76         clr     c
  77         subb    a, exp_b
  78         // if no carry then no underflow
  79         jnc     00008$
  80         add     a, #127
  81         jc      00009$
  82         ljmp    fs_return_zero
  83
  84 00008$:
  85         add     a, #128
  86         dec     a
  87         jnc     00009$
  88         ljmp    fs_return_inf
  89
  90 00009$:
  91         mov     exp_a, a
  92
  93         // need extra bits on a's mantissa
  94 #ifdef FLOAT_FULL_ACCURACY
  95         clr     c
  96         mov     a, r5
  97         subb    a, r2
  98         mov     a, r6
  99         subb    a, r3
 100         mov     a, r7
 101         subb    a, r4
 102         jc      00010$
 103         dec     exp_a
 104         clr     c
 105         mov     a, r2
 106         rlc     a
 107         mov     r1, a
 108         mov     a, r3
 109         rlc     a
 110         mov     r2, a
 111         mov     a, r4
 112         rlc     a
 113         mov     r3, a
 114         clr     a
 115         rlc     a
 116         mov     r4, a
 117         sjmp    00011$
 118 00010$:
 119 #endif
 120         clr     a
 121         xch     a, r4
 122         xch     a, r3
 123         xch     a, r2
 124         mov     r1, a
 125 00011$:
 126
 127         // begin long division
 128         push    exp_a
 129 #ifdef FLOAT_FULL_ACCURACY
 130         mov     b, #25
 131 #else
 132         mov     b, #24
 133 #endif
 134 00012$:
 135         // compare
 136         clr     c
 137         mov     a, r1
 138         subb    a, r5
 139         mov     a, r2
 140         subb    a, r6
 141         mov     a, r3
 142         subb    a, r7
 143         mov     a, r4
 144         subb    a, #0           // carry==0 if mant1 >= mant2
 145
 146 #ifdef FLOAT_FULL_ACCURACY
 147         djnz    b, 00013$
 148         sjmp    00015$
 149 00013$:
 150 #endif
 151         jc      00014$
 152         // subtract
 153         mov     a, r1
 154         subb    a, r5
 155         mov     r1, a
 156         mov     a, r2
 157         subb    a, r6
 158         mov     r2, a
 159         mov     a, r3
 160         subb    a, r7
 161         mov     r3, a
 162         mov     a, r4
 163         subb    a, #0
 164         mov     r4, a
 165         clr     c
 166
 167 00014$:
 168         // shift result
 169         cpl     c
 170         mov     a, r0
 171         rlc     a
 172         mov     r0, a
 173         mov     a, dpl
 174         rlc     a
 175         mov     dpl, a
 176         mov     a, dph
 177         rlc     a
 178         mov     dph, a
 179
 180         // shift partial remainder
 181         clr     c
 182         mov     a, r1
 183         rlc     a
 184         mov     r1, a
 185         mov     a, r2
 186         rlc     a
 187         mov     r2, a
 188         mov     a, r3
 189         rlc     a
 190         mov     r3, a
 191         mov     a, r4
 192         rlc     a
 193         mov     r4, a
 194
 195 #ifdef FLOAT_FULL_ACCURACY
 196         sjmp    00012$
 197 00015$:
 198 #else
 199         djnz    b, 00012$
 200 #endif
 201
 202         // now we've got a division result, so all we need to do
 203         // is round off properly, normalize and output a float
 204
 205 #ifdef FLOAT_FULL_ACCURACY
 206         cpl     c
 207         clr     a
 208         mov     r1, a
 209         addc    a, r0
 210         mov     r2, a
 211         clr     a
 212         addc    a, dpl
 213         mov     r3, a
 214         clr     a
 215         addc    a, dph
 216         mov     r4, a
 217         pop     exp_a
 218         jnc     00016$
 219         inc     exp_a
 220         // incrementing exp_a without checking carry is dangerous
 221         mov     r4, #0x80
 222 00016$:
 223 #else
 224         mov     r1, #0
 225         mov     a, r0
 226         mov     r2, a
 227         mov     r3, dpl
 228         mov     r4, dph
 229         pop     exp_a
 230 #endif
 231
 232         lcall   fs_normalize_a
 233         ljmp    fs_zerocheck_return
 234         __endasm;
 235 }
 236
 237 #else
 238
 239 /*
 240 ** libgcc support for software floating point.
 241 ** Copyright (C) 1991 by Pipeline Associates, Inc.  All rights reserved.
 242 ** Permission is granted to do *anything* you want with this file,
 243 ** commercial or otherwise, provided this message remains intact.  So there!
 244 ** I would appreciate receiving any updates/patches/changes that anyone
 245 ** makes, and am willing to be the repository for said changes (am I
 246 ** making a big mistake?).
 247 **
 248 ** Pat Wood
 249 ** Pipeline Associates, Inc.
 250 ** pipeline!phw@motown.com or
 251 ** sun!pipeline!phw or
 252 ** uunet!motown!pipeline!phw
 253 */
 254
 255 /* (c)2000/2001: hacked a little by johan.knol@iduna.nl for sdcc */
 256
 257 union float_long
 258   {
 259     float f;
 260     long l;
 261   };
 262
 263 /* divide two floats */
 264 float __fsdiv (float a1, float a2)
 265 {
 266   volatile union float_long fl1, fl2;
 267   volatile long result;
 268   volatile unsigned long mask;
 269   volatile long mant1, mant2;
 270   volatile int exp;
 271   char sign;
 272
 273   fl1.f = a1;
 274   fl2.f = a2;
 275
 276   /* subtract exponents */
 277   exp = EXP (fl1.l) ;
 278   exp -= EXP (fl2.l);
 279   exp += EXCESS;
 280
 281   /* compute sign */
 282   sign = SIGN (fl1.l) ^ SIGN (fl2.l);
 283
 284   /* divide by zero??? */
 285   if (!fl2.l)
 286     {/* return NaN or -NaN */
 287       fl2.l = 0x7FC00000;
 288       return (fl2.f);
 289     }
 290
 291   /* numerator zero??? */
 292   if (!fl1.l)
 293     return (0);
 294
 295   /* now get mantissas */
 296   mant1 = MANT (fl1.l);
 297   mant2 = MANT (fl2.l);
 298
 299   /* this assures we have 25 bits of precision in the end */
 300   if (mant1 < mant2)
 301     {
 302       mant1 <<= 1;
 303       exp--;
 304     }
 305
 306   /* now we perform repeated subtraction of fl2.l from fl1.l */
 307   mask = 0x1000000;
 308   result = 0;
 309   while (mask)
 310     {
 311       if (mant1 >= mant2)
 312         {
 313           result |= mask;
 314           mant1 -= mant2;
 315         }
 316       mant1 <<= 1;
 317       mask >>= 1;
 318     }
 319
 320   /* round */
 321   result += 1;
 322
 323   /* normalize down */
 324   exp++;
 325   result >>= 1;
 326
 327   result &= ~HIDDEN;
 328
 329   /* pack up and go home */
 330   if (exp >= 0x100)
 331     fl1.l = (sign ? SIGNBIT : 0) | __INFINITY;
 332   else if (exp < 0)
 333     fl1.l = 0;
 334   else
 335     fl1.l = PACK (sign ? SIGNBIT : 0 , exp, result);
 336   return (fl1.f);
 337 }
 338
 339 #endif