Return to floatlib.c CVS log

Up to [Apple XNU] / GNUtools / cc

Annotation of GNUtools/cc/floatlib.c, revision 1.1

1.1     ! root        1: /*
        !             2: ** libgcc support for software floating point.
        !             3: ** Copyright (C) 1991 by Pipeline Associates, Inc.  All rights reserved.
        !             4: ** Permission is granted to do *anything* you want with this file,
        !             5: ** commercial or otherwise, provided this message remains intact.  So there!
        !             6: ** I would appreciate receiving any updates/patches/changes that anyone
        !             7: ** makes, and am willing to be the repository for said changes (am I
        !             8: ** making a big mistake?).
        !             9: 
        !            10: Warning! Only single-precision is actually implemented.  This file
        !            11: won't really be much use until double-precision is supported.
        !            12: 
        !            13: However, once that is done, this file might eventually become a
        !            14: replacement for libgcc1.c.  It might also make possible
        !            15: cross-compilation for an IEEE target machine from a non-IEEE
        !            16: host such as a VAX.
        !            17: 
        !            18: If you'd like to work on completing this, please talk to [email protected].
        !            19: 
        !            20: 
        !            21: **
        !            22: ** Pat Wood
        !            23: ** Pipeline Associates, Inc.
        !            24: ** [email protected] or
        !            25: ** sun!pipeline!phw or
        !            26: ** uunet!motown!pipeline!phw
        !            27: **
        !            28: ** 05/01/91 -- V1.0 -- first release to gcc mailing lists
        !            29: ** 05/04/91 -- V1.1 -- added float and double prototypes and return values
        !            30: **                  -- fixed problems with adding and subtracting zero
        !            31: **                  -- fixed rounding in truncdfsf2
        !            32: **                  -- fixed SWAP define and tested on 386
        !            33: */
        !            34: 
        !            35: /*
        !            36: ** The following are routines that replace the libgcc soft floating point
        !            37: ** routines that are called automatically when -msoft-float is selected.
        !            38: ** The support single and double precision IEEE format, with provisions
        !            39: ** for byte-swapped machines (tested on 386).  Some of the double-precision
        !            40: ** routines work at full precision, but most of the hard ones simply punt
        !            41: ** and call the single precision routines, producing a loss of accuracy.
        !            42: ** long long support is not assumed or included.
        !            43: ** Overall accuracy is close to IEEE (actually 68882) for single-precision
        !            44: ** arithmetic.  I think there may still be a 1 in 1000 chance of a bit
        !            45: ** being rounded the wrong way during a multiply.  I'm not fussy enough to
        !            46: ** bother with it, but if anyone is, knock yourself out.
        !            47: **
        !            48: ** Efficiency has only been addressed where it was obvious that something
        !            49: ** would make a big difference.  Anyone who wants to do this right for
        !            50: ** best speed should go in and rewrite in assembler.
        !            51: **
        !            52: ** I have tested this only on a 68030 workstation and 386/ix integrated
        !            53: ** in with -msoft-float.
        !            54: */
        !            55: 
        !            56: /* the following deal with IEEE single-precision numbers */
        !            57: #define EXCESS         126
        !            58: #define SIGNBIT                0x80000000
        !            59: #define HIDDEN         (1 << 23)
        !            60: #define SIGN(fp)       ((fp) & SIGNBIT)
        !            61: #define EXP(fp)                (((fp) >> 23) & 0xFF)
        !            62: #define MANT(fp)       (((fp) & 0x7FFFFF) | HIDDEN)
        !            63: #define PACK(s,e,m)    ((s) | ((e) << 23) | (m))
        !            64: 
        !            65: /* the following deal with IEEE double-precision numbers */
        !            66: #define EXCESSD                1022
        !            67: #define HIDDEND                (1 << 20)
        !            68: #define EXPD(fp)       (((fp.l.upper) >> 20) & 0x7FF)
        !            69: #define SIGND(fp)      ((fp.l.upper) & SIGNBIT)
        !            70: #define MANTD(fp)      (((((fp.l.upper) & 0xFFFFF) | HIDDEND) << 10) | \
        !            71:                                (fp.l.lower >> 22))
        !            72: 
        !            73: /* define SWAP for 386/960 reverse-byte-order brain-damaged CPUs */
        !            74: union double_long
        !            75:   {
        !            76:     double d;
        !            77: #ifdef SWAP
        !            78:     struct {
        !            79:       unsigned long lower;
        !            80:       long upper;
        !            81:     } l;
        !            82: #else
        !            83:     struct {
        !            84:       long upper;
        !            85:       unsigned long lower;
        !            86:     } l;
        !            87: #endif
        !            88:   };
        !            89: 
        !            90: union float_long
        !            91:   {
        !            92:     float f;
        !            93:     long l;
        !            94:   };
        !            95: 
        !            96: /* add two floats */
        !            97: float
        !            98: __addsf3 (float a1, float a2)
        !            99: {
        !           100:   register long mant1, mant2;
        !           101:   register union float_long fl1, fl2;
        !           102:   register int exp1, exp2;
        !           103:   int sign = 0;
        !           104: 
        !           105:   fl1.f = a1;
        !           106:   fl2.f = a2;
        !           107: 
        !           108:   /* check for zero args */
        !           109:   if (!fl1.l)
        !           110:     return (fl2.f);
        !           111:   if (!fl2.l)
        !           112:     return (fl1.f);
        !           113: 
        !           114:   exp1 = EXP (fl1.l);
        !           115:   exp2 = EXP (fl2.l);
        !           116: 
        !           117:   if (exp1 > exp2 + 25)
        !           118:     return (fl1.l);
        !           119:   if (exp2 > exp1 + 25)
        !           120:     return (fl2.l);
        !           121: 
        !           122:   /* do everything in excess precision so's we can round later */
        !           123:   mant1 = MANT (fl1.l) << 6;
        !           124:   mant2 = MANT (fl2.l) << 6;
        !           125: 
        !           126:   if (SIGN (fl1.l))
        !           127:     mant1 = -mant1;
        !           128:   if (SIGN (fl2.l))
        !           129:     mant2 = -mant2;
        !           130: 
        !           131:   if (exp1 > exp2)
        !           132:     {
        !           133:       mant2 >>= exp1 - exp2;
        !           134:     }
        !           135:   else
        !           136:     {
        !           137:       mant1 >>= exp2 - exp1;
        !           138:       exp1 = exp2;
        !           139:     }
        !           140:   mant1 += mant2;
        !           141: 
        !           142:   if (mant1 < 0)
        !           143:     {
        !           144:       mant1 = -mant1;
        !           145:       sign = SIGNBIT;
        !           146:     }
        !           147:   else if (!mant1)
        !           148:     return (0);
        !           149: 
        !           150:   /* normalize up */
        !           151:   while (!(mant1 & 0xE0000000))
        !           152:     {
        !           153:       mant1 <<= 1;
        !           154:       exp1--;
        !           155:     }
        !           156: 
        !           157:   /* normalize down? */
        !           158:   if (mant1 & (1 << 30))
        !           159:     {
        !           160:       mant1 >>= 1;
        !           161:       exp1++;
        !           162:     }
        !           163: 
        !           164:   /* round to even */
        !           165:   mant1 += (mant1 & 0x40) ? 0x20 : 0x1F;
        !           166: 
        !           167:   /* normalize down? */
        !           168:   if (mant1 & (1 << 30))
        !           169:     {
        !           170:       mant1 >>= 1;
        !           171:       exp1++;
        !           172:     }
        !           173: 
        !           174:   /* lose extra precision */
        !           175:   mant1 >>= 6;
        !           176: 
        !           177:   /* turn off hidden bit */
        !           178:   mant1 &= ~HIDDEN;
        !           179: 
        !           180:   /* pack up and go home */
        !           181:   fl1.l = PACK (sign, exp1, mant1);
        !           182:   return (fl1.f);
        !           183: }
        !           184: 
        !           185: /* subtract two floats */
        !           186: float
        !           187: __subsf3 (float a1, float a2)
        !           188: {
        !           189:   register union float_long fl1, fl2;
        !           190: 
        !           191:   fl1.f = a1;
        !           192:   fl2.f = a2;
        !           193: 
        !           194:   /* check for zero args */
        !           195:   if (!fl2.l)
        !           196:     return (fl1.f);
        !           197:   if (!fl1.l)
        !           198:     return (-fl2.f);
        !           199: 
        !           200:   /* twiddle sign bit and add */
        !           201:   fl2.l ^= SIGNBIT;
        !           202:   return __addsf3 (a1, fl2.f);
        !           203: }
        !           204: 
        !           205: /* compare two floats */
        !           206: long
        !           207: __cmpsf2 (float a1, float a2)
        !           208: {
        !           209:   register union float_long fl1, fl2;
        !           210: 
        !           211:   fl1.f = a1;
        !           212:   fl2.f = a2;
        !           213: 
        !           214:   if (SIGN (fl1.l) && SIGN (fl2.l))
        !           215:     {
        !           216:       fl1.l ^= SIGNBIT;
        !           217:       fl2.l ^= SIGNBIT;
        !           218:     }
        !           219:   if (fl1.l < fl2.l)
        !           220:     return (-1);
        !           221:   if (fl1.l > fl2.l)
        !           222:     return (1);
        !           223:   return (0);
        !           224: }
        !           225: 
        !           226: /* multiply two floats */
        !           227: float
        !           228: __mulsf3 (float a1, float a2)
        !           229: {
        !           230:   register union float_long fl1, fl2;
        !           231:   register unsigned long result;
        !           232:   register int exp;
        !           233:   int sign;
        !           234: 
        !           235:   fl1.f = a1;
        !           236:   fl2.f = a2;
        !           237: 
        !           238:   if (!fl1.l || !fl2.l)
        !           239:     return (0);
        !           240: 
        !           241:   /* compute sign and exponent */
        !           242:   sign = SIGN (fl1.l) ^ SIGN (fl2.l);
        !           243:   exp = EXP (fl1.l) - EXCESS;
        !           244:   exp += EXP (fl2.l);
        !           245: 
        !           246:   fl1.l = MANT (fl1.l);
        !           247:   fl2.l = MANT (fl2.l);
        !           248: 
        !           249:   /* the multiply is done as one 16x16 multiply and two 16x8 multiples */
        !           250:   result = (fl1.l >> 8) * (fl2.l >> 8);
        !           251:   result += ((fl1.l & 0xFF) * (fl2.l >> 8)) >> 8;
        !           252:   result += ((fl2.l & 0xFF) * (fl1.l >> 8)) >> 8;
        !           253: 
        !           254:   if (result & 0x80000000)
        !           255:     {
        !           256:       /* round */
        !           257:       result += 0x80;
        !           258:       result >>= 8;
        !           259:     }
        !           260:   else
        !           261:     {
        !           262:       /* round */
        !           263:       result += 0x40;
        !           264:       result >>= 7;
        !           265:       exp--;
        !           266:     }
        !           267: 
        !           268:   result &= ~HIDDEN;
        !           269: 
        !           270:   /* pack up and go home */
        !           271:   fl1.l = PACK (sign, exp, result);
        !           272:   return (fl1.f);
        !           273: }
        !           274: 
        !           275: /* divide two floats */
        !           276: float
        !           277: __divsf3 (float a1, float a2)
        !           278: {
        !           279:   register union float_long fl1, fl2;
        !           280:   register int result;
        !           281:   register int mask;
        !           282:   register int exp, sign;
        !           283: 
        !           284:   fl1.f = a1;
        !           285:   fl2.f = a2;
        !           286: 
        !           287:   /* subtract exponents */
        !           288:   exp = EXP (fl1.l) - EXP (fl2.l) + EXCESS;
        !           289: 
        !           290:   /* compute sign */
        !           291:   sign = SIGN (fl1.l) ^ SIGN (fl2.l);
        !           292: 
        !           293:   /* divide by zero??? */
        !           294:   if (!fl2.l)
        !           295:     /* return NaN or -NaN */
        !           296:     return (sign ? 0xFFFFFFFF : 0x7FFFFFFF);
        !           297: 
        !           298:   /* numerator zero??? */
        !           299:   if (!fl1.l)
        !           300:     return (0);
        !           301: 
        !           302:   /* now get mantissas */
        !           303:   fl1.l = MANT (fl1.l);
        !           304:   fl2.l = MANT (fl2.l);
        !           305: 
        !           306:   /* this assures we have 25 bits of precision in the end */
        !           307:   if (fl1.l < fl2.l)
        !           308:     {
        !           309:       fl1.l <<= 1;
        !           310:       exp--;
        !           311:     }
        !           312: 
        !           313:   /* now we perform repeated subtraction of fl2.l from fl1.l */
        !           314:   mask = 0x1000000;
        !           315:   result = 0;
        !           316:   while (mask)
        !           317:     {
        !           318:       if (fl1.l >= fl2.l)
        !           319:        {
        !           320:          result |= mask;
        !           321:          fl1.l -= fl2.l;
        !           322:        }
        !           323:       fl1.l <<= 1;
        !           324:       mask >>= 1;
        !           325:     }
        !           326: 
        !           327:   /* round */
        !           328:   result += 1;
        !           329: 
        !           330:   /* normalize down */
        !           331:   exp++;
        !           332:   result >>= 1;
        !           333: 
        !           334:   result &= ~HIDDEN;
        !           335: 
        !           336:   /* pack up and go home */
        !           337:   fl1.l = PACK (sign, exp, result);
        !           338:   return (fl1.f);
        !           339: }
        !           340: 
        !           341: /* convert int to double */
        !           342: double
        !           343: __floatsidf (register long a1)
        !           344: {
        !           345:   register int sign = 0, exp = 31 + EXCESSD;
        !           346:   union double_long dl;
        !           347: 
        !           348:   if (!a1)
        !           349:     {
        !           350:       dl.l.upper = dl.l.lower = 0;
        !           351:       return (dl.d);
        !           352:     }
        !           353: 
        !           354:   if (a1 < 0)
        !           355:     {
        !           356:       sign = SIGNBIT;
        !           357:       a1 = -a1;
        !           358:     }
        !           359: 
        !           360:   while (a1 < 0x1000000)
        !           361:     {
        !           362:       a1 <<= 4;
        !           363:       exp -= 4;
        !           364:     }
        !           365: 
        !           366:   while (a1 < 0x40000000)
        !           367:     {
        !           368:       a1 <<= 1;
        !           369:       exp--;
        !           370:     }
        !           371: 
        !           372:   /* pack up and go home */
        !           373:   dl.l.upper = sign;
        !           374:   dl.l.upper |= exp << 20;
        !           375:   dl.l.upper |= (a1 >> 10) & ~HIDDEND;
        !           376:   dl.l.lower = a1 << 22;
        !           377: 
        !           378:   return (dl.d);
        !           379: }
        !           380: 
        !           381: /* negate a float */
        !           382: float
        !           383: __negsf2 (float a1)
        !           384: {
        !           385:   register union float_long fl1;
        !           386: 
        !           387:   fl1.f = a1;
        !           388:   if (!fl1.l)
        !           389:     return (0);
        !           390: 
        !           391:   fl1.l ^= SIGNBIT;
        !           392:   return (fl1.f);
        !           393: }
        !           394: 
        !           395: /* negate a double */
        !           396: double
        !           397: __negdf2 (double a1)
        !           398: {
        !           399:   register union double_long dl1;
        !           400: 
        !           401:   dl1.d = a1;
        !           402: 
        !           403:   if (!dl1.l.upper && !dl1.l.lower)
        !           404:       return (dl1.d);
        !           405: 
        !           406:   dl1.l.upper ^= SIGNBIT;
        !           407:   return (dl1.d);
        !           408: }
        !           409: 
        !           410: /* convert float to double */
        !           411: double
        !           412: __extendsfdf2 (float a1)
        !           413: {
        !           414:   register union float_long fl1;
        !           415:   register union double_long dl;
        !           416:   register int exp;
        !           417: 
        !           418:   fl1.f = a1;
        !           419: 
        !           420:   if (!fl1.l)
        !           421:     {
        !           422:       dl.l.upper = dl.l.lower = 0;
        !           423:       return (dl.d);
        !           424:     }
        !           425: 
        !           426:   dl.l.upper = SIGN (fl1.l);
        !           427:   exp = EXP (fl1.l) - EXCESS + EXCESSD;
        !           428:   dl.l.upper |= exp << 20;
        !           429:   dl.l.upper |= (MANT (fl1.l) & ~HIDDEN) >> 3;
        !           430:   dl.l.lower = MANT (fl1.l) << 29;
        !           431: 
        !           432:   return (dl.d);
        !           433: }
        !           434: 
        !           435: /* convert double to float */
        !           436: float
        !           437: __truncdfsf2 (double a1)
        !           438: {
        !           439:   register int exp;
        !           440:   register long mant;
        !           441:   register union float_long fl;
        !           442:   register union double_long dl1;
        !           443: 
        !           444:   dl1.d = a1;
        !           445: 
        !           446:   if (!dl1.l.upper && !dl1.l.lower)
        !           447:     return (0);
        !           448: 
        !           449:   exp = EXPD (dl1) - EXCESSD + EXCESS;
        !           450: 
        !           451:   /* shift double mantissa 6 bits so we can round */
        !           452:   mant = MANTD (dl1) >> 6;
        !           453: 
        !           454:   /* now round and shift down */
        !           455:   mant += 1;
        !           456:   mant >>= 1;
        !           457: 
        !           458:   /* did the round overflow? */
        !           459:   if (mant & 0xFF000000)
        !           460:     {
        !           461:       mant >>= 1;
        !           462:       exp++;
        !           463:     }
        !           464: 
        !           465:   mant &= ~HIDDEN;
        !           466: 
        !           467:   /* pack up and go home */
        !           468:   fl.l = PACK (SIGND (dl1), exp, mant);
        !           469:   return (fl.f);
        !           470: }
        !           471: 
        !           472: /* compare two doubles */
        !           473: long
        !           474: __cmpdf2 (double a1, double a2)
        !           475: {
        !           476:   register union double_long dl1, dl2;
        !           477: 
        !           478:   dl1.d = a1;
        !           479:   dl2.d = a2;
        !           480: 
        !           481:   if (SIGND (dl1) && SIGND (dl2))
        !           482:     {
        !           483:       dl1.l.upper ^= SIGNBIT;
        !           484:       dl2.l.upper ^= SIGNBIT;
        !           485:     }
        !           486:   if (dl1.l.upper < dl2.l.upper)
        !           487:     return (-1);
        !           488:   if (dl1.l.upper > dl2.l.upper)
        !           489:     return (1);
        !           490:   if (dl1.l.lower < dl2.l.lower)
        !           491:     return (-1);
        !           492:   if (dl1.l.lower > dl2.l.lower)
        !           493:     return (1);
        !           494:   return (0);
        !           495: }
        !           496: 
        !           497: /* convert double to int */
        !           498: long
        !           499: __fixdfsi (double a1)
        !           500: {
        !           501:   register union double_long dl1;
        !           502:   register int exp;
        !           503:   register long l;
        !           504: 
        !           505:   dl1.d = a1;
        !           506: 
        !           507:   if (!dl1.l.upper && !dl1.l.lower)
        !           508:     return (0);
        !           509: 
        !           510:   exp = EXPD (dl1) - EXCESSD - 31;
        !           511:   l = MANTD (dl1);
        !           512: 
        !           513:   if (exp > 0)
        !           514:     return (0x7FFFFFFF | SIGND (dl1)); /* largest integer */
        !           515: 
        !           516:   /* shift down until exp = 0 or l = 0 */
        !           517:   if (exp < 0 && exp > -32 && l)
        !           518:     l >>= -exp;
        !           519:   else
        !           520:     return (0);
        !           521: 
        !           522:   return (SIGND (dl1) ? -l : l);
        !           523: }
        !           524: 
        !           525: /* convert double to unsigned int */
        !           526: unsigned
        !           527: long __fixunsdfsi (double a1)
        !           528: {
        !           529:   register union double_long dl1;
        !           530:   register int exp;
        !           531:   register unsigned long l;
        !           532: 
        !           533:   dl1.d = a1;
        !           534: 
        !           535:   if (!dl1.l.upper && !dl1.l.lower)
        !           536:     return (0);
        !           537: 
        !           538:   exp = EXPD (dl1) - EXCESSD - 32;
        !           539:   l = (((((dl1.l.upper) & 0xFFFFF) | HIDDEND) << 11) | (dl1.l.lower >> 21));
        !           540: 
        !           541:   if (exp > 0)
        !           542:     return (0xFFFFFFFF);       /* largest integer */
        !           543: 
        !           544:   /* shift down until exp = 0 or l = 0 */
        !           545:   if (exp < 0 && exp > -32 && l)
        !           546:     l >>= -exp;
        !           547:   else
        !           548:     return (0);
        !           549: 
        !           550:   return (l);
        !           551: }
        !           552: 
        !           553: /* For now, the hard double-precision routines simply
        !           554:    punt and do it in single */
        !           555: /* addtwo doubles */
        !           556: double
        !           557: __adddf3 (double a1, double a2)
        !           558: {
        !           559:   return ((float) a1 + (float) a2);
        !           560: }
        !           561: 
        !           562: /* subtract two doubles */
        !           563: double
        !           564: __subdf3 (double a1, double a2)
        !           565: {
        !           566:   return ((float) a1 - (float) a2);
        !           567: }
        !           568: 
        !           569: /* multiply two doubles */
        !           570: double
        !           571: __muldf3 (double a1, double a2)
        !           572: {
        !           573:   return ((float) a1 * (float) a2);
        !           574: }
        !           575: 
        !           576: /* divide two doubles */
        !           577: double
        !           578: __divdf3 (double a1, double a2)
        !           579: {
        !           580:   return ((float) a1 / (float) a2);
        !           581: }

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.