X-Git-Url: https://git.gag.com/?a=blobdiff_plain;f=device%2Flib%2F_fsdiv.c;h=5a1a21a744387cdceb947e01f97f4fd461f0dea8;hb=025f549fbb89aa070e3f3baeb6c8ef4ca72d8210;hp=e56e4ddad7541ca17480c984fc9fa5825396a425;hpb=b09af35f2f1cde7649d3ac4a6f5d2af6d97895a0;p=fw%2Fsdcc diff --git a/device/lib/_fsdiv.c b/device/lib/_fsdiv.c index e56e4dda..5a1a21a7 100644 --- a/device/lib/_fsdiv.c +++ b/device/lib/_fsdiv.c @@ -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 + + +#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. @@ -6,61 +244,15 @@ ** I would appreciate receiving any updates/patches/changes that anyone ** makes, and am willing to be the repository for said changes (am I ** making a big mistake?). - -Warning! Only single-precision is actually implemented. This file -won't really be much use until double-precision is supported. - -However, once that is done, this file might eventually become a -replacement for libgcc1.c. It might also make possible -cross-compilation for an IEEE target machine from a non-IEEE -host such as a VAX. - -If you'd like to work on completing this, please talk to rms@gnu.ai.mit.edu. - - ** ** Pat Wood ** Pipeline Associates, Inc. ** pipeline!phw@motown.com or ** sun!pipeline!phw or ** uunet!motown!pipeline!phw -** -** 05/01/91 -- V1.0 -- first release to gcc mailing lists -** 05/04/91 -- V1.1 -- added float and double prototypes and return values -** -- fixed problems with adding and subtracting zero -** -- fixed rounding in truncdfsf2 -** -- fixed SWAP define and tested on 386 */ -/* -** The following are routines that replace the libgcc soft floating point -** routines that are called automatically when -msoft-float is selected. -** The support single and double precision IEEE format, with provisions -** for byte-swapped machines (tested on 386). Some of the double-precision -** routines work at full precision, but most of the hard ones simply punt -** and call the single precision routines, producing a loss of accuracy. -** long long support is not assumed or included. -** Overall accuracy is close to IEEE (actually 68882) for single-precision -** arithmetic. I think there may still be a 1 in 1000 chance of a bit -** being rounded the wrong way during a multiply. I'm not fussy enough to -** bother with it, but if anyone is, knock yourself out. -** -** Efficiency has only been addressed where it was obvious that something -** would make a big difference. Anyone who wants to do this right for -** best speed should go in and rewrite in assembler. -** -** I have tested this only on a 68030 workstation and 386/ix integrated -** in with -msoft-float. -*/ - -/* the following deal with IEEE single-precision numbers */ -#define EXCESS 126 -#define SIGNBIT ((unsigned long)0x80000000) -#define HIDDEN (unsigned long)(1 << 23) -#define SIGN(fp) ((fp >> (8*sizeof(fp)-1)) & 1) -#define EXP(fp) (((fp) >> 23) & (unsigned int)0x00FF) -#define MANT(fp) (((fp) & (unsigned long)0x007FFFFF) | HIDDEN) -#define PACK(s,e,m) ((s) | ((e) << 23) | (m)) +/* (c)2000/2001: hacked a little by johan.knol@iduna.nl for sdcc */ union float_long { @@ -69,15 +261,14 @@ union float_long }; /* divide two floats */ -float -__fsdiv (float a1, float a2) +float __fsdiv (float a1, float a2) { volatile union float_long fl1, fl2; volatile long result; volatile unsigned long mask; volatile long mant1, mant2; - volatile int exp ; - short sign; + volatile int exp; + char sign; fl1.f = a1; fl2.f = a2; @@ -92,8 +283,10 @@ __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) @@ -134,6 +327,13 @@ __fsdiv (float a1, float a2) result &= ~HIDDEN; /* pack up and go home */ - fl1.l = PACK (sign, (unsigned long) exp, result); + if (exp >= 0x100) + fl1.l = (sign ? SIGNBIT : 0) | 0x7F800000; + else if (exp < 0) + fl1.l = 0; + else + fl1.l = PACK (sign ? SIGNBIT : 0 , exp, result); return (fl1.f); } + +#endif