* device/include/float.h: added __INFINITY
[fw/sdcc] / device / lib / _fsdiv.c
index a89d698b5f0d4efe93297d9ea775c69ee1a1c34d..7bf1c48f010293a75054c269f0f63a06bcb0ac53 100644 (file)
@@ -1,3 +1,241 @@
+/* Floating point library in optimized assembly for 8051
+ * Copyright (c) 2004, Paul Stoffregen, paul@pjrc.com
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Library General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This library is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+ */
+
+
+#define SDCC_FLOAT_LIB
+#include <float.h>
+
+
+#ifdef FLOAT_ASM_MCS51
+
+// float __fsdiv (float a, float b) __reentrant
+static void dummy(void) __naked
+{
+       __asm
+       .globl  ___fsdiv
+___fsdiv:
+       // extract the two inputs, placing them into:
+       //      sign     exponent   mantiassa
+       //      ----     --------   ---------
+       //  a:  sign_a   exp_a      r4/r3/r2
+       //  b:  sign_b   exp_b      r7/r6/r5
+
+       lcall   fsgetargs
+
+       // compute final sign bit
+       jnb     sign_b, 00001$
+       cpl     sign_a
+00001$:
+
+       // if divisor is zero, ...
+       cjne    r7, #0, 00003$
+       // if dividend is also zero, return NaN
+       cjne    r4, #0, 00002$
+       ljmp    fs_return_nan
+00002$:
+       // but dividend is non-zero, return infinity
+       ljmp    fs_return_inf
+00003$:
+       // if dividend is zero, return zero
+       cjne    r4, #0, 00004$
+       ljmp    fs_return_zero
+00004$:
+       // if divisor is infinity, ...
+       mov     a, exp_b
+       cjne    a, #0xFF, 00006$
+       // and dividend is also infinity, return NaN
+       mov     a, exp_a
+       cjne    a, #0xFF, 00005$
+       ljmp    fs_return_nan
+00005$:
+       // but dividend is not infinity, return zero
+       ljmp    fs_return_zero
+00006$:
+       // if dividend is infinity, return infinity
+       mov     a, exp_a
+       cjne    a, #0xFF, 00007$
+       ljmp    fs_return_inf
+00007$:
+
+       // subtract exponents
+       clr     c
+       subb    a, exp_b
+       // if no carry then no underflow
+       jnc     00008$
+       add     a, #127
+       jc      00009$
+       ljmp    fs_return_zero
+
+00008$:
+       add     a, #128
+       dec     a
+       jnc     00009$
+       ljmp    fs_return_inf
+
+00009$:
+       mov     exp_a, a
+
+       // need extra bits on a's mantissa
+#ifdef FLOAT_FULL_ACCURACY
+       clr     c
+       mov     a, r5
+       subb    a, r2
+       mov     a, r6
+       subb    a, r3
+       mov     a, r7
+       subb    a, r4
+       jc      00010$
+       dec     exp_a
+       clr     c
+       mov     a, r2
+       rlc     a
+       mov     r1, a
+       mov     a, r3
+       rlc     a
+       mov     r2, a
+       mov     a, r4
+       rlc     a
+       mov     r3, a
+       clr     a
+       rlc     a
+       mov     r4, a
+       sjmp    00011$
+00010$:
+#endif
+       clr     a
+       xch     a, r4
+       xch     a, r3
+       xch     a, r2
+       mov     r1, a
+00011$:
+
+       // begin long division
+       push    exp_a
+#ifdef FLOAT_FULL_ACCURACY
+       mov     b, #25
+#else
+       mov     b, #24
+#endif
+00012$:
+       // compare
+       clr     c
+       mov     a, r1
+       subb    a, r5
+       mov     a, r2
+       subb    a, r6
+       mov     a, r3
+       subb    a, r7
+       mov     a, r4
+       subb    a, #0           // carry==0 if mant1 >= mant2
+
+#ifdef FLOAT_FULL_ACCURACY
+       djnz    b, 00013$
+       sjmp    00015$
+00013$:
+#endif
+       jc      00014$
+       // subtract
+       mov     a, r1
+       subb    a, r5
+       mov     r1, a
+       mov     a, r2
+       subb    a, r6
+       mov     r2, a
+       mov     a, r3
+       subb    a, r7
+       mov     r3, a
+       mov     a, r4
+       subb    a, #0
+       mov     r4, a
+       clr     c
+
+00014$:
+       // shift result
+       cpl     c
+       mov     a, r0
+       rlc     a
+       mov     r0, a
+       mov     a, dpl
+       rlc     a
+       mov     dpl, a
+       mov     a, dph
+       rlc     a
+       mov     dph, a
+
+       // shift partial remainder
+       clr     c
+       mov     a, r1
+       rlc     a
+       mov     r1, a
+       mov     a, r2
+       rlc     a
+       mov     r2, a
+       mov     a, r3
+       rlc     a
+       mov     r3, a
+       mov     a, r4
+       rlc     a
+       mov     r4, a
+
+#ifdef FLOAT_FULL_ACCURACY
+       sjmp    00012$
+00015$:
+#else
+       djnz    b, 00012$
+#endif
+
+       // now we've got a division result, so all we need to do
+       // is round off properly, normalize and output a float
+
+#ifdef FLOAT_FULL_ACCURACY
+       cpl     c
+       clr     a
+       mov     r1, a
+       addc    a, r0
+       mov     r2, a
+       clr     a
+       addc    a, dpl
+       mov     r3, a
+       clr     a
+       addc    a, dph
+       mov     r4, a
+       pop     exp_a
+       jnc     00016$
+       inc     exp_a
+       // incrementing exp_a without checking carry is dangerous
+       mov     r4, #0x80
+00016$:
+#else
+       mov     r1, #0
+       mov     a, r0
+       mov     r2, a
+       mov     r3, dpl
+       mov     r4, dph
+       pop     exp_a
+#endif
+
+       lcall   fs_normalize_a
+       ljmp    fs_zerocheck_return
+       __endasm;
+}
+
+#else
+
 /*
 ** libgcc support for software floating point.
 ** Copyright (C) 1991 by Pipeline Associates, Inc.  All rights reserved.
 
 /* (c)2000/2001: hacked a little by johan.knol@iduna.nl for sdcc */
 
-#include <float.h>
-
 union float_long
   {
     float f;
@@ -31,7 +267,7 @@ float __fsdiv (float a1, float a2)
   volatile long result;
   volatile unsigned long mask;
   volatile long mant1, mant2;
-  volatile int exp ;
+  volatile int exp;
   char sign;
 
   fl1.f = a1;
@@ -47,8 +283,10 @@ float __fsdiv (float a1, float a2)
 
   /* divide by zero??? */
   if (!fl2.l)
-    /* return NaN or -NaN */
-    return (-1.0);
+    {/* return NaN or -NaN */
+      fl2.l = 0x7FC00000;
+      return (fl2.f);
+    }
 
   /* numerator zero??? */
   if (!fl1.l)
@@ -89,7 +327,13 @@ float __fsdiv (float a1, float a2)
   result &= ~HIDDEN;
 
   /* pack up and go home */
-  fl1.l = PACK (sign ? 1ul<<31 : 0, (unsigned long) exp, result);
+  if (exp >= 0x100)
+    fl1.l = (sign ? SIGNBIT : 0) | __INFINITY;
+  else if (exp < 0)
+    fl1.l = 0;
+  else
+    fl1.l = PACK (sign ? SIGNBIT : 0 , exp, result);
   return (fl1.f);
 }
 
+#endif