Return to fyl2x.c CVS log

Up to [Previous NeXT emulator] / previous / src / softfloat

Annotation of previous/src/softfloat/fyl2x.c, revision 1.1

1.1     ! root        1: /*============================================================================
        !             2: This source file is an extension to the SoftFloat IEC/IEEE Floating-point
        !             3: Arithmetic Package, Release 2b, written for Bochs (x86 achitecture simulator)
        !             4: floating point emulation.
        !             5: float_raise(float_flag_invalid)
        !             6: THIS SOFTWARE IS DISTRIBUTED AS IS, FOR FREE.  Although reasonable effort has
        !             7: been made to avoid it, THIS SOFTWARE MAY CONTAIN FAULTS THAT WILL AT TIMES
        !             8: RESULT IN INCORRECT BEHAVIOR.  USE OF THIS SOFTWARE IS RESTRICTED TO PERSONS
        !             9: AND ORGANIZATIONS WHO CAN AND WILL TAKE FULL RESPONSIBILITY FOR ALL LOSSES,
        !            10: COSTS, OR OTHER PROBLEMS THEY INCUR DUE TO THE SOFTWARE, AND WHO FURTHERMORE
        !            11: EFFECTIVELY INDEMNIFY JOHN HAUSER AND THE INTERNATIONAL COMPUTER SCIENCE
        !            12: INSTITUTE (possibly via similar legal warning) AGAINST ALL LOSSES, COSTS, OR
        !            13: OTHER PROBLEMS INCURRED BY THEIR CUSTOMERS AND CLIENTS DUE TO THE SOFTWARE.
        !            14: 
        !            15: Derivative works are acceptable, even for commercial purposes, so long as
        !            16: (1) the source code for the derivative work includes prominent notice that
        !            17: the work is derivative, and (2) the source code includes prominent notice with
        !            18: these four paragraphs for those parts of this code that are retained.
        !            19: =============================================================================*/
        !            20: 
        !            21: /*============================================================================
        !            22:  * Written for Bochs (x86 achitecture simulator) by
        !            23:  *            Stanislav Shwartsman [sshwarts at sourceforge net]
        !            24:  * Adapted for lib/softfloat in MESS by Hans Ostermeyer (03/2012)
        !            25:  * ==========================================================================*/
        !            26: 
        !            27: #define FLOAT128
        !            28: 
        !            29: #define USE_estimateDiv128To64
        !            30: #include "mamesf.h"
        !            31: #include "softfloat.h"
        !            32: #include "fpu_constant.h"
        !            33: 
        !            34: static const floatx80 floatx80_log10_2 = { 0x3ffd, 0x9a209a84fbcff798U };
        !            35: static const floatx80 floatx80_ln_2 = { 0x3ffe, 0xb17217f7d1cf79acU };
        !            36: static const floatx80 floatx80_one = { 0x3fff, 0x8000000000000000U };
        !            37: static const floatx80 floatx80_default_nan = { floatx80_default_nan_high, floatx80_default_nan_low };
        !            38: 
        !            39: #define packFloat_128(zHi, zLo) {(zHi), (zLo)}
        !            40: #define PACK_FLOAT_128(hi,lo) packFloat_128(LIT64(hi),LIT64(lo))
        !            41: 
        !            42: #define EXP_BIAS 0x3FFF
        !            43: 
        !            44: 
        !            45: static const float128 float128_one =
        !            46:        packFloat_128(0x3fff000000000000U, 0x0000000000000000U);
        !            47: static const float128 float128_two =
        !            48:        packFloat_128(0x4000000000000000U, 0x0000000000000000U);
        !            49: 
        !            50: static const float128 float128_ln2inv2 =
        !            51:        packFloat_128(0x400071547652b82fU, 0xe1777d0ffda0d23aU);
        !            52: 
        !            53: #define SQRT2_HALF_SIG  0xb504f333f9de6484U
        !            54: 
        !            55: extern float128 OddPoly(float128 x, float128 *arr, unsigned n);
        !            56: 
        !            57: #define L2_ARR_SIZE 9
        !            58: 
        !            59: static float128 ln_arr[L2_ARR_SIZE] =
        !            60: {
        !            61:        PACK_FLOAT_128(0x3fff000000000000, 0x0000000000000000), /*  1 */
        !            62:        PACK_FLOAT_128(0x3ffd555555555555, 0x5555555555555555), /*  3 */
        !            63:        PACK_FLOAT_128(0x3ffc999999999999, 0x999999999999999a), /*  5 */
        !            64:        PACK_FLOAT_128(0x3ffc249249249249, 0x2492492492492492), /*  7 */
        !            65:        PACK_FLOAT_128(0x3ffbc71c71c71c71, 0xc71c71c71c71c71c), /*  9 */
        !            66:        PACK_FLOAT_128(0x3ffb745d1745d174, 0x5d1745d1745d1746), /* 11 */
        !            67:        PACK_FLOAT_128(0x3ffb3b13b13b13b1, 0x3b13b13b13b13b14), /* 13 */
        !            68:        PACK_FLOAT_128(0x3ffb111111111111, 0x1111111111111111), /* 15 */
        !            69:        PACK_FLOAT_128(0x3ffae1e1e1e1e1e1, 0xe1e1e1e1e1e1e1e2)  /* 17 */
        !            70: };
        !            71: 
        !            72: static float128 poly_ln(float128 x1)
        !            73: {
        !            74: /*
        !            75:     //
        !            76:     //                     3     5     7     9     11     13     15
        !            77:     //        1+u         u     u     u     u     u      u      u
        !            78:     // 1/2 ln ---  ~ u + --- + --- + --- + --- + ---- + ---- + ---- =
        !            79:     //        1-u         3     5     7     9     11     13     15
        !            80:     //
        !            81:     //                     2     4     6     8     10     12     14
        !            82:     //                    u     u     u     u     u      u      u
        !            83:     //       = u * [ 1 + --- + --- + --- + --- + ---- + ---- + ---- ] =
        !            84:     //                    3     5     7     9     11     13     15
        !            85:     //
        !            86:     //           3                          3
        !            87:     //          --       4k                --        4k+2
        !            88:     //   p(u) = >  C  * u           q(u) = >  C   * u
        !            89:     //          --  2k                     --  2k+1
        !            90:     //          k=0                        k=0
        !            91:     //
        !            92:     //          1+u                 2
        !            93:     //   1/2 ln --- ~ u * [ p(u) + u * q(u) ]
        !            94:     //          1-u
        !            95:     //
        !            96: */
        !            97:        return OddPoly(x1, ln_arr, L2_ARR_SIZE);
        !            98: }
        !            99: 
        !           100: /* required sqrt(2)/2 < x < sqrt(2) */
        !           101: static float128 poly_l2(float128 x)
        !           102: {
        !           103:        /* using float128 for approximation */
        !           104:        float128 x_p1 = float128_add(x, float128_one);
        !           105:        float128 x_m1 = float128_sub(x, float128_one);
        !           106:        x = float128_div(x_m1, x_p1);
        !           107:        x = poly_ln(x);
        !           108:        x = float128_mul(x, float128_ln2inv2);
        !           109:        return x;
        !           110: }
        !           111: 
        !           112: static float128 poly_l2p1(float128 x)
        !           113: {
        !           114:        /* using float128 for approximation */
        !           115:        float128 x_p2 = float128_add(x, float128_two);
        !           116:        x = float128_div(x, x_p2);
        !           117:        x = poly_ln(x);
        !           118:        x = float128_mul(x, float128_ln2inv2);
        !           119:        return x;
        !           120: }
        !           121: 
        !           122: // =================================================
        !           123: // FYL2X                   Compute y * log (x)
        !           124: //                                        2
        !           125: // =================================================
        !           126: 
        !           127: //
        !           128: // Uses the following identities:
        !           129: //
        !           130: // 1. ----------------------------------------------------------
        !           131: //              ln(x)
        !           132: //   log (x) = -------,  ln (x*y) = ln(x) + ln(y)
        !           133: //      2       ln(2)
        !           134: //
        !           135: // 2. ----------------------------------------------------------
        !           136: //                1+u             x-1
        !           137: //   ln (x) = ln -----, when u = -----
        !           138: //                1-u             x+1
        !           139: //
        !           140: // 3. ----------------------------------------------------------
        !           141: //                        3     5     7           2n+1
        !           142: //       1+u             u     u     u           u
        !           143: //   ln ----- = 2 [ u + --- + --- + --- + ... + ------ + ... ]
        !           144: //       1-u             3     5     7           2n+1
        !           145: //
        !           146: 
        !           147: static floatx80 fyl2x(floatx80 a, floatx80 b)
        !           148: {
        !           149:        uint64_t aSig = extractFloatx80Frac(a);
        !           150:        int32_t aExp = extractFloatx80Exp(a);
        !           151:        int aSign = extractFloatx80Sign(a);
        !           152:        uint64_t bSig = extractFloatx80Frac(b);
        !           153:        int32_t bExp = extractFloatx80Exp(b);
        !           154:        int bSign = extractFloatx80Sign(b);
        !           155: 
        !           156:        int zSign = bSign ^ 1;
        !           157: 
        !           158:        if (aExp == 0x7FFF) {
        !           159:                if ((uint64_t) (aSig<<1)
        !           160:                                || ((bExp == 0x7FFF) && (uint64_t) (bSig<<1)))
        !           161:                {
        !           162:                        return propagateFloatx80NaN(a, b);
        !           163:                }
        !           164:                if (aSign)
        !           165:                {
        !           166: invalid:
        !           167:                        float_raise(float_flag_invalid);
        !           168:                        return floatx80_default_nan;
        !           169:                }
        !           170:                else {
        !           171:                        if (bExp == 0) {
        !           172:                                if (bSig == 0) goto invalid;
        !           173:                                float_raise(float_flag_denormal);
        !           174:                        }
        !           175:                        return packFloatx80(bSign, 0x7FFF, floatx80_default_infinity_low);
        !           176:                }
        !           177:        }
        !           178:        if (bExp == 0x7FFF)
        !           179:        {
        !           180:                if ((uint64_t) (bSig<<1)) return propagateFloatx80NaN(a, b);
        !           181:                if (aSign && (uint64_t)(aExp | aSig)) goto invalid;
        !           182:                if (aSig && (aExp == 0))
        !           183:                        float_raise(float_flag_denormal);
        !           184:                if (aExp < 0x3FFF) {
        !           185:                        return packFloatx80(zSign, 0x7FFF, floatx80_default_infinity_low);
        !           186:                }
        !           187:                if (aExp == 0x3FFF && ((uint64_t) (aSig<<1) == 0)) goto invalid;
        !           188:                return packFloatx80(bSign, 0x7FFF, floatx80_default_infinity_low);
        !           189:        }
        !           190:        if (aExp == 0) {
        !           191:                if (aSig == 0) {
        !           192:                        if ((bExp | bSig) == 0) goto invalid;
        !           193:                        float_raise(float_flag_divbyzero);
        !           194:                        return packFloatx80(zSign, 0x7FFF, floatx80_default_infinity_low);
        !           195:                }
        !           196:                if (aSign) goto invalid;
        !           197:                float_raise(float_flag_denormal);
        !           198:                normalizeFloatx80Subnormal(aSig, &aExp, &aSig);
        !           199:        }
        !           200:        if (aSign) goto invalid;
        !           201:        if (bExp == 0) {
        !           202:                if (bSig == 0) {
        !           203:                        if (aExp < 0x3FFF) return packFloatx80(zSign, 0, 0);
        !           204:                        return packFloatx80(bSign, 0, 0);
        !           205:                }
        !           206:                float_raise(float_flag_denormal);
        !           207:                normalizeFloatx80Subnormal(bSig, &bExp, &bSig);
        !           208:        }
        !           209:        if (aExp == 0x3FFF && ((uint64_t) (aSig<<1) == 0))
        !           210:                return packFloatx80(bSign, 0, 0);
        !           211: 
        !           212:        float_raise(float_flag_inexact);
        !           213: 
        !           214:        int ExpDiff = aExp - 0x3FFF;
        !           215:        aExp = 0;
        !           216:        if (aSig >= SQRT2_HALF_SIG) {
        !           217:                ExpDiff++;
        !           218:                aExp--;
        !           219:        }
        !           220: 
        !           221:        /* ******************************** */
        !           222:        /* using float128 for approximation */
        !           223:        /* ******************************** */
        !           224: 
        !           225:        uint64_t zSig0, zSig1;
        !           226:        shift128Right(aSig<<1, 0, 16, &zSig0, &zSig1);
        !           227:        float128 x = packFloat128(0, aExp+0x3FFF, zSig0, zSig1);
        !           228:        x = poly_l2(x);
        !           229:        x = float128_add(x, int64_to_float128((int64_t) ExpDiff));
        !           230:        return floatx80_mul(b, float128_to_floatx80(x));
        !           231: }
        !           232: 
        !           233: // =================================================
        !           234: // FYL2XP1                 Compute y * log (x + 1)
        !           235: //                                        2
        !           236: // =================================================
        !           237: 
        !           238: //
        !           239: // Uses the following identities:
        !           240: //
        !           241: // 1. ----------------------------------------------------------
        !           242: //              ln(x)
        !           243: //   log (x) = -------
        !           244: //      2       ln(2)
        !           245: //
        !           246: // 2. ----------------------------------------------------------
        !           247: //                  1+u              x
        !           248: //   ln (x+1) = ln -----, when u = -----
        !           249: //                  1-u             x+2
        !           250: //
        !           251: // 3. ----------------------------------------------------------
        !           252: //                        3     5     7           2n+1
        !           253: //       1+u             u     u     u           u
        !           254: //   ln ----- = 2 [ u + --- + --- + --- + ... + ------ + ... ]
        !           255: //       1-u             3     5     7           2n+1
        !           256: //
        !           257: 
        !           258: floatx80 fyl2xp1(floatx80 a, floatx80 b)
        !           259: {
        !           260:        int32_t aExp, bExp;
        !           261:        uint64_t aSig, bSig, zSig0, zSig1, zSig2;
        !           262:        int aSign, bSign;
        !           263: 
        !           264:        aSig = extractFloatx80Frac(a);
        !           265:        aExp = extractFloatx80Exp(a);
        !           266:        aSign = extractFloatx80Sign(a);
        !           267:        bSig = extractFloatx80Frac(b);
        !           268:        bExp = extractFloatx80Exp(b);
        !           269:        bSign = extractFloatx80Sign(b);
        !           270:        int zSign = aSign ^ bSign;
        !           271: 
        !           272:        if (aExp == 0x7FFF) {
        !           273:                if ((uint64_t) (aSig<<1)
        !           274:                                || ((bExp == 0x7FFF) && (uint64_t) (bSig<<1)))
        !           275:                {
        !           276:                        return propagateFloatx80NaN(a, b);
        !           277:                }
        !           278:                if (aSign)
        !           279:                {
        !           280: invalid:
        !           281:                        float_raise(float_flag_invalid);
        !           282:                        return floatx80_default_nan;
        !           283:                }
        !           284:                        else {
        !           285:                        if (bExp == 0) {
        !           286:                                if (bSig == 0) goto invalid;
        !           287:                                float_raise(float_flag_denormal);
        !           288:                        }
        !           289:                        return packFloatx80(bSign, 0x7FFF, floatx80_default_infinity_low);
        !           290:                }
        !           291:        }
        !           292:        if (bExp == 0x7FFF)
        !           293:        {
        !           294:                if ((uint64_t) (bSig<<1))
        !           295:                        return propagateFloatx80NaN(a, b);
        !           296: 
        !           297:                if (aExp == 0) {
        !           298:                        if (aSig == 0) goto invalid;
        !           299:                        float_raise(float_flag_denormal);
        !           300:                }
        !           301: 
        !           302:                return packFloatx80(zSign, 0x7FFF, floatx80_default_infinity_low);
        !           303:        }
        !           304:        if (aExp == 0) {
        !           305:                if (aSig == 0) {
        !           306:                        if (bSig && (bExp == 0)) float_raise(float_flag_denormal);
        !           307:                        return packFloatx80(zSign, 0, 0);
        !           308:                }
        !           309:                float_raise(float_flag_denormal);
        !           310:                normalizeFloatx80Subnormal(aSig, &aExp, &aSig);
        !           311:        }
        !           312:        if (bExp == 0) {
        !           313:                if (bSig == 0) return packFloatx80(zSign, 0, 0);
        !           314:                float_raise(float_flag_denormal);
        !           315:                normalizeFloatx80Subnormal(bSig, &bExp, &bSig);
        !           316:        }
        !           317: 
        !           318:        float_raise(float_flag_inexact);
        !           319: 
        !           320:        if (aSign && aExp >= 0x3FFF)
        !           321:                return a;
        !           322: 
        !           323:        if (aExp >= 0x3FFC) // big argument
        !           324:        {
        !           325:                return fyl2x(floatx80_add(a, floatx80_one), b);
        !           326:        }
        !           327: 
        !           328:        // handle tiny argument
        !           329:        if (aExp < EXP_BIAS-70)
        !           330:        {
        !           331:                // first order approximation, return (a*b)/ln(2)
        !           332:                int32_t zExp = aExp + FLOAT_LN2INV_EXP - 0x3FFE;
        !           333: 
        !           334:        mul128By64To192(FLOAT_LN2INV_HI, FLOAT_LN2INV_LO, aSig, &zSig0, &zSig1, &zSig2);
        !           335:                if (0 < (int64_t) zSig0) {
        !           336:                        shortShift128Left(zSig0, zSig1, 1, &zSig0, &zSig1);
        !           337:                        --zExp;
        !           338:                }
        !           339: 
        !           340:                zExp = zExp + bExp - 0x3FFE;
        !           341:        mul128By64To192(zSig0, zSig1, bSig, &zSig0, &zSig1, &zSig2);
        !           342:                if (0 < (int64_t) zSig0) {
        !           343:                        shortShift128Left(zSig0, zSig1, 1, &zSig0, &zSig1);
        !           344:                        --zExp;
        !           345:                }
        !           346: 
        !           347:                return
        !           348:                        roundAndPackFloatx80(80, aSign ^ bSign, zExp, zSig0, zSig1);
        !           349:        }
        !           350: 
        !           351:        /* ******************************** */
        !           352:        /* using float128 for approximation */
        !           353:        /* ******************************** */
        !           354: 
        !           355:        shift128Right(aSig<<1, 0, 16, &zSig0, &zSig1);
        !           356:        float128 x = packFloat128(aSign, aExp, zSig0, zSig1);
        !           357:        x = poly_l2p1(x);
        !           358:        return floatx80_mul(b, float128_to_floatx80(x));
        !           359: }
        !           360: 
        !           361: floatx80 floatx80_flognp1(floatx80 a)
        !           362: {
        !           363:        return fyl2xp1(a, floatx80_ln_2);
        !           364: }
        !           365: 
        !           366: floatx80 floatx80_flogn(floatx80 a)
        !           367: {
        !           368:        return fyl2x(a, floatx80_ln_2);
        !           369: }
        !           370: 
        !           371: floatx80 floatx80_flog2(floatx80 a)
        !           372: {
        !           373:        return fyl2x(a, floatx80_one);
        !           374: }
        !           375: 
        !           376: floatx80 floatx80_flog10(floatx80 a)
        !           377: {
        !           378:        return fyl2x(a, floatx80_log10_2);
        !           379: }

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