X-Git-Url: https://git.gag.com/?a=blobdiff_plain;f=device%2Flib%2F_fsdiv.c;h=7bf1c48f010293a75054c269f0f63a06bcb0ac53;hb=a8a7fc8a514bc79114c6d76d6a6aeb885387478c;hp=a89d698b5f0d4efe93297d9ea775c69ee1a1c34d;hpb=5a964d04cfff5accede4e488135633bd1a780576;p=fw%2Fsdcc diff --git a/device/lib/_fsdiv.c b/device/lib/_fsdiv.c index a89d698b..7bf1c48f 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. @@ -16,8 +254,6 @@ /* (c)2000/2001: hacked a little by johan.knol@iduna.nl for sdcc */ -#include - 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