Annotation of nono/fpe/fpu_log.c, revision 1.1

1.1     ! root        1: /*     $NetBSD: fpu_log.c,v 1.18 2014/01/04 13:23:22 isaki Exp $       */
        !             2: 
        !             3: /*
        !             4:  * Copyright (c) 1995  Ken Nakata
        !             5:  *     All rights reserved.
        !             6:  *
        !             7:  * Redistribution and use in source and binary forms, with or without
        !             8:  * modification, are permitted provided that the following conditions
        !             9:  * are met:
        !            10:  * 1. Redistributions of source code must retain the above copyright
        !            11:  *    notice, this list of conditions and the following disclaimer.
        !            12:  * 2. Redistributions in binary form must reproduce the above copyright
        !            13:  *    notice, this list of conditions and the following disclaimer in the
        !            14:  *    documentation and/or other materials provided with the distribution.
        !            15:  * 3. Neither the name of the author nor the names of its contributors
        !            16:  *    may be used to endorse or promote products derived from this software
        !            17:  *    without specific prior written permission.
        !            18:  *
        !            19:  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
        !            20:  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
        !            21:  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
        !            22:  * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
        !            23:  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
        !            24:  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
        !            25:  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
        !            26:  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
        !            27:  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
        !            28:  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
        !            29:  * SUCH DAMAGE.
        !            30:  *
        !            31:  *     @(#)fpu_log.c   10/8/95
        !            32:  */
        !            33: 
        !            34: #include "fpu_emulate.h"
        !            35: 
        !            36: static uint32_t logA6[] = { 0x3FC2499A, 0xB5E4040B };
        !            37: static uint32_t logA5[] = { 0xBFC555B5, 0x848CB7DB };
        !            38: static uint32_t logA4[] = { 0x3FC99999, 0x987D8730 };
        !            39: static uint32_t logA3[] = { 0xBFCFFFFF, 0xFF6F7E97 };
        !            40: static uint32_t logA2[] = { 0x3FD55555, 0x555555A4 };
        !            41: static uint32_t logA1[] = { 0xBFE00000, 0x00000008 };
        !            42: 
        !            43: static uint32_t logB5[] = { 0x3F175496, 0xADD7DAD6 };
        !            44: static uint32_t logB4[] = { 0x3F3C71C2, 0xFE80C7E0 };
        !            45: static uint32_t logB3[] = { 0x3F624924, 0x928BCCFF };
        !            46: static uint32_t logB2[] = { 0x3F899999, 0x999995EC };
        !            47: static uint32_t logB1[] = { 0x3FB55555, 0x55555555 };
        !            48: 
        !            49: /* sfpn = shortened fp number; can represent only positive numbers */
        !            50: static struct sfpn {
        !            51:        int             sp_exp;
        !            52:        uint32_t        sp_m0, sp_m1;
        !            53: } logtbl[] = {
        !            54:        { 0x3FFE - 0x3fff, 0xFE03F80FU, 0xE03F80FEU },
        !            55:        { 0x3FF7 - 0x3fff, 0xFF015358U, 0x833C47E2U },
        !            56:        { 0x3FFE - 0x3fff, 0xFA232CF2U, 0x52138AC0U },
        !            57:        { 0x3FF9 - 0x3fff, 0xBDC8D83EU, 0xAD88D549U },
        !            58:        { 0x3FFE - 0x3fff, 0xF6603D98U, 0x0F6603DAU },
        !            59:        { 0x3FFA - 0x3fff, 0x9CF43DCFU, 0xF5EAFD48U },
        !            60:        { 0x3FFE - 0x3fff, 0xF2B9D648U, 0x0F2B9D65U },
        !            61:        { 0x3FFA - 0x3fff, 0xDA16EB88U, 0xCB8DF614U },
        !            62:        { 0x3FFE - 0x3fff, 0xEF2EB71FU, 0xC4345238U },
        !            63:        { 0x3FFB - 0x3fff, 0x8B29B775U, 0x1BD70743U },
        !            64:        { 0x3FFE - 0x3fff, 0xEBBDB2A5U, 0xC1619C8CU },
        !            65:        { 0x3FFB - 0x3fff, 0xA8D839F8U, 0x30C1FB49U },
        !            66:        { 0x3FFE - 0x3fff, 0xE865AC7BU, 0x7603A197U },
        !            67:        { 0x3FFB - 0x3fff, 0xC61A2EB1U, 0x8CD907ADU },
        !            68:        { 0x3FFE - 0x3fff, 0xE525982AU, 0xF70C880EU },
        !            69:        { 0x3FFB - 0x3fff, 0xE2F2A47AU, 0xDE3A18AFU },
        !            70:        { 0x3FFE - 0x3fff, 0xE1FC780EU, 0x1FC780E2U },
        !            71:        { 0x3FFB - 0x3fff, 0xFF64898EU, 0xDF55D551U },
        !            72:        { 0x3FFE - 0x3fff, 0xDEE95C4CU, 0xA037BA57U },
        !            73:        { 0x3FFC - 0x3fff, 0x8DB956A9U, 0x7B3D0148U },
        !            74:        { 0x3FFE - 0x3fff, 0xDBEB61EEU, 0xD19C5958U },
        !            75:        { 0x3FFC - 0x3fff, 0x9B8FE100U, 0xF47BA1DEU },
        !            76:        { 0x3FFE - 0x3fff, 0xD901B203U, 0x6406C80EU },
        !            77:        { 0x3FFC - 0x3fff, 0xA9372F1DU, 0x0DA1BD17U },
        !            78:        { 0x3FFE - 0x3fff, 0xD62B80D6U, 0x2B80D62CU },
        !            79:        { 0x3FFC - 0x3fff, 0xB6B07F38U, 0xCE90E46BU },
        !            80:        { 0x3FFE - 0x3fff, 0xD3680D36U, 0x80D3680DU },
        !            81:        { 0x3FFC - 0x3fff, 0xC3FD0329U, 0x06488481U },
        !            82:        { 0x3FFE - 0x3fff, 0xD0B69FCBU, 0xD2580D0BU },
        !            83:        { 0x3FFC - 0x3fff, 0xD11DE0FFU, 0x15AB18CAU },
        !            84:        { 0x3FFE - 0x3fff, 0xCE168A77U, 0x25080CE1U },
        !            85:        { 0x3FFC - 0x3fff, 0xDE1433A1U, 0x6C66B150U },
        !            86:        { 0x3FFE - 0x3fff, 0xCB8727C0U, 0x65C393E0U },
        !            87:        { 0x3FFC - 0x3fff, 0xEAE10B5AU, 0x7DDC8ADDU },
        !            88:        { 0x3FFE - 0x3fff, 0xC907DA4EU, 0x871146ADU },
        !            89:        { 0x3FFC - 0x3fff, 0xF7856E5EU, 0xE2C9B291U },
        !            90:        { 0x3FFE - 0x3fff, 0xC6980C69U, 0x80C6980CU },
        !            91:        { 0x3FFD - 0x3fff, 0x82012CA5U, 0xA68206D7U },
        !            92:        { 0x3FFE - 0x3fff, 0xC4372F85U, 0x5D824CA6U },
        !            93:        { 0x3FFD - 0x3fff, 0x882C5FCDU, 0x7256A8C5U },
        !            94:        { 0x3FFE - 0x3fff, 0xC1E4BBD5U, 0x95F6E947U },
        !            95:        { 0x3FFD - 0x3fff, 0x8E44C60BU, 0x4CCFD7DEU },
        !            96:        { 0x3FFE - 0x3fff, 0xBFA02FE8U, 0x0BFA02FFU },
        !            97:        { 0x3FFD - 0x3fff, 0x944AD09EU, 0xF4351AF6U },
        !            98:        { 0x3FFE - 0x3fff, 0xBD691047U, 0x07661AA3U },
        !            99:        { 0x3FFD - 0x3fff, 0x9A3EECD4U, 0xC3EAA6B2U },
        !           100:        { 0x3FFE - 0x3fff, 0xBB3EE721U, 0xA54D880CU },
        !           101:        { 0x3FFD - 0x3fff, 0xA0218434U, 0x353F1DE8U },
        !           102:        { 0x3FFE - 0x3fff, 0xB92143FAU, 0x36F5E02EU },
        !           103:        { 0x3FFD - 0x3fff, 0xA5F2FCABU, 0xBBC506DAU },
        !           104:        { 0x3FFE - 0x3fff, 0xB70FBB5AU, 0x19BE3659U },
        !           105:        { 0x3FFD - 0x3fff, 0xABB3B8BAU, 0x2AD362A5U },
        !           106:        { 0x3FFE - 0x3fff, 0xB509E68AU, 0x9B94821FU },
        !           107:        { 0x3FFD - 0x3fff, 0xB1641795U, 0xCE3CA97BU },
        !           108:        { 0x3FFE - 0x3fff, 0xB30F6352U, 0x8917C80BU },
        !           109:        { 0x3FFD - 0x3fff, 0xB7047551U, 0x5D0F1C61U },
        !           110:        { 0x3FFE - 0x3fff, 0xB11FD3B8U, 0x0B11FD3CU },
        !           111:        { 0x3FFD - 0x3fff, 0xBC952AFEU, 0xEA3D13E1U },
        !           112:        { 0x3FFE - 0x3fff, 0xAF3ADDC6U, 0x80AF3ADEU },
        !           113:        { 0x3FFD - 0x3fff, 0xC2168ED0U, 0xF458BA4AU },
        !           114:        { 0x3FFE - 0x3fff, 0xAD602B58U, 0x0AD602B6U },
        !           115:        { 0x3FFD - 0x3fff, 0xC788F439U, 0xB3163BF1U },
        !           116:        { 0x3FFE - 0x3fff, 0xAB8F69E2U, 0x8359CD11U },
        !           117:        { 0x3FFD - 0x3fff, 0xCCECAC08U, 0xBF04565DU },
        !           118:        { 0x3FFE - 0x3fff, 0xA9C84A47U, 0xA07F5638U },
        !           119:        { 0x3FFD - 0x3fff, 0xD2420487U, 0x2DD85160U },
        !           120:        { 0x3FFE - 0x3fff, 0xA80A80A8U, 0x0A80A80BU },
        !           121:        { 0x3FFD - 0x3fff, 0xD7894992U, 0x3BC3588AU },
        !           122:        { 0x3FFE - 0x3fff, 0xA655C439U, 0x2D7B73A8U },
        !           123:        { 0x3FFD - 0x3fff, 0xDCC2C4B4U, 0x9887DACCU },
        !           124:        { 0x3FFE - 0x3fff, 0xA4A9CF1DU, 0x96833751U },
        !           125:        { 0x3FFD - 0x3fff, 0xE1EEBD3EU, 0x6D6A6B9EU },
        !           126:        { 0x3FFE - 0x3fff, 0xA3065E3FU, 0xAE7CD0E0U },
        !           127:        { 0x3FFD - 0x3fff, 0xE70D785CU, 0x2F9F5BDCU },
        !           128:        { 0x3FFE - 0x3fff, 0xA16B312EU, 0xA8FC377DU },
        !           129:        { 0x3FFD - 0x3fff, 0xEC1F392CU, 0x5179F283U },
        !           130:        { 0x3FFE - 0x3fff, 0x9FD809FDU, 0x809FD80AU },
        !           131:        { 0x3FFD - 0x3fff, 0xF12440D3U, 0xE36130E6U },
        !           132:        { 0x3FFE - 0x3fff, 0x9E4CAD23U, 0xDD5F3A20U },
        !           133:        { 0x3FFD - 0x3fff, 0xF61CCE92U, 0x346600BBU },
        !           134:        { 0x3FFE - 0x3fff, 0x9CC8E160U, 0xC3FB19B9U },
        !           135:        { 0x3FFD - 0x3fff, 0xFB091FD3U, 0x8145630AU },
        !           136:        { 0x3FFE - 0x3fff, 0x9B4C6F9EU, 0xF03A3CAAU },
        !           137:        { 0x3FFD - 0x3fff, 0xFFE97042U, 0xBFA4C2ADU },
        !           138:        { 0x3FFE - 0x3fff, 0x99D722DAU, 0xBDE58F06U },
        !           139:        { 0x3FFE - 0x3fff, 0x825EFCEDU, 0x49369330U },
        !           140:        { 0x3FFE - 0x3fff, 0x9868C809U, 0x868C8098U },
        !           141:        { 0x3FFE - 0x3fff, 0x84C37A7AU, 0xB9A905C9U },
        !           142:        { 0x3FFE - 0x3fff, 0x97012E02U, 0x5C04B809U },
        !           143:        { 0x3FFE - 0x3fff, 0x87224C2EU, 0x8E645FB7U },
        !           144:        { 0x3FFE - 0x3fff, 0x95A02568U, 0x095A0257U },
        !           145:        { 0x3FFE - 0x3fff, 0x897B8CACU, 0x9F7DE298U },
        !           146:        { 0x3FFE - 0x3fff, 0x94458094U, 0x45809446U },
        !           147:        { 0x3FFE - 0x3fff, 0x8BCF55DEU, 0xC4CD05FEU },
        !           148:        { 0x3FFE - 0x3fff, 0x92F11384U, 0x0497889CU },
        !           149:        { 0x3FFE - 0x3fff, 0x8E1DC0FBU, 0x89E125E5U },
        !           150:        { 0x3FFE - 0x3fff, 0x91A2B3C4U, 0xD5E6F809U },
        !           151:        { 0x3FFE - 0x3fff, 0x9066E68CU, 0x955B6C9BU },
        !           152:        { 0x3FFE - 0x3fff, 0x905A3863U, 0x3E06C43BU },
        !           153:        { 0x3FFE - 0x3fff, 0x92AADE74U, 0xC7BE59E0U },
        !           154:        { 0x3FFE - 0x3fff, 0x8F1779D9U, 0xFDC3A219U },
        !           155:        { 0x3FFE - 0x3fff, 0x94E9BFF6U, 0x15845643U },
        !           156:        { 0x3FFE - 0x3fff, 0x8DDA5202U, 0x37694809U },
        !           157:        { 0x3FFE - 0x3fff, 0x9723A1B7U, 0x20134203U },
        !           158:        { 0x3FFE - 0x3fff, 0x8CA29C04U, 0x6514E023U },
        !           159:        { 0x3FFE - 0x3fff, 0x995899C8U, 0x90EB8990U },
        !           160:        { 0x3FFE - 0x3fff, 0x8B70344AU, 0x139BC75AU },
        !           161:        { 0x3FFE - 0x3fff, 0x9B88BDAAU, 0x3A3DAE2FU },
        !           162:        { 0x3FFE - 0x3fff, 0x8A42F870U, 0x5669DB46U },
        !           163:        { 0x3FFE - 0x3fff, 0x9DB4224FU, 0xFFE1157CU },
        !           164:        { 0x3FFE - 0x3fff, 0x891AC73AU, 0xE9819B50U },
        !           165:        { 0x3FFE - 0x3fff, 0x9FDADC26U, 0x8B7A12DAU },
        !           166:        { 0x3FFE - 0x3fff, 0x87F78087U, 0xF78087F8U },
        !           167:        { 0x3FFE - 0x3fff, 0xA1FCFF17U, 0xCE733BD4U },
        !           168:        { 0x3FFE - 0x3fff, 0x86D90544U, 0x7A34ACC6U },
        !           169:        { 0x3FFE - 0x3fff, 0xA41A9E8FU, 0x5446FB9FU },
        !           170:        { 0x3FFE - 0x3fff, 0x85BF3761U, 0x2CEE3C9BU },
        !           171:        { 0x3FFE - 0x3fff, 0xA633CD7EU, 0x6771CD8BU },
        !           172:        { 0x3FFE - 0x3fff, 0x84A9F9C8U, 0x084A9F9DU },
        !           173:        { 0x3FFE - 0x3fff, 0xA8489E60U, 0x0B435A5EU },
        !           174:        { 0x3FFE - 0x3fff, 0x83993052U, 0x3FBE3368U },
        !           175:        { 0x3FFE - 0x3fff, 0xAA59233CU, 0xCCA4BD49U },
        !           176:        { 0x3FFE - 0x3fff, 0x828CBFBEU, 0xB9A020A3U },
        !           177:        { 0x3FFE - 0x3fff, 0xAC656DAEU, 0x6BCC4985U },
        !           178:        { 0x3FFE - 0x3fff, 0x81848DA8U, 0xFAF0D277U },
        !           179:        { 0x3FFE - 0x3fff, 0xAE6D8EE3U, 0x60BB2468U },
        !           180:        { 0x3FFE - 0x3fff, 0x80808080U, 0x80808081U },
        !           181:        { 0x3FFE - 0x3fff, 0xB07197A2U, 0x3C46C654U },
        !           182: };
        !           183: 
        !           184: static struct fpn *__fpu_logn(struct fpemu *fe);
        !           185: 
        !           186: /*
        !           187:  * natural log - algorithm taken from Motorola FPSP,
        !           188:  * except this doesn't bother to check for invalid input.
        !           189:  */
        !           190: static struct fpn *
        !           191: __fpu_logn(struct fpemu *fe)
        !           192: {
        !           193:        static struct fpn X, F, U, V, W, KLOG2;
        !           194:        struct fpn *d;
        !           195:        int i, k;
        !           196: 
        !           197:        CPYFPN(&X, &fe->fe_f2);
        !           198: 
        !           199:        /* see if |X-1| < 1/16 approx. */
        !           200:        if ((-1 == X.fp_exp && (0xf07d0000U >> (31 - FP_LG)) <= X.fp_mant[0]) ||
        !           201:            (0 == X.fp_exp && X.fp_mant[0] <= (0x88410000U >> (31 - FP_LG)))) {
        !           202:                /* log near 1 */
        !           203: #if FPE_DEBUG
        !           204:                printf("__fpu_logn: log near 1\n");
        !           205: #endif
        !           206: 
        !           207:                fpu_const(&fe->fe_f1, FPU_CONST_1);
        !           208:                /* X+1 */
        !           209:                d = fpu_add(fe);
        !           210:                CPYFPN(&V, d);
        !           211: 
        !           212:                CPYFPN(&fe->fe_f1, &X);
        !           213:                fpu_const(&fe->fe_f2, FPU_CONST_1);
        !           214:                fe->fe_f2.fp_sign = 1; /* -1.0 */
        !           215:                /* X-1 */
        !           216:                d = fpu_add(fe);
        !           217:                CPYFPN(&fe->fe_f1, d);
        !           218:                /* 2(X-1) */
        !           219:                fe->fe_f1.fp_exp++; /* *= 2 */
        !           220:                CPYFPN(&fe->fe_f2, &V);
        !           221:                /* U=2(X-1)/(X+1) */
        !           222:                d = fpu_div(fe);
        !           223:                CPYFPN(&U, d);
        !           224:                CPYFPN(&fe->fe_f1, d);
        !           225:                CPYFPN(&fe->fe_f2, d);
        !           226:                /* V=U*U */
        !           227:                d = fpu_mul(fe);
        !           228:                CPYFPN(&V, d);
        !           229:                CPYFPN(&fe->fe_f1, d);
        !           230:                CPYFPN(&fe->fe_f2, d);
        !           231:                /* W=V*V */
        !           232:                d = fpu_mul(fe);
        !           233:                CPYFPN(&W, d);
        !           234: 
        !           235:                /* calculate U+U*V*([B1+W*(B3+W*B5)]+[V*(B2+W*B4)]) */
        !           236: 
        !           237:                /* B1+W*(B3+W*B5) part */
        !           238:                CPYFPN(&fe->fe_f1, d);
        !           239:                fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logB5);
        !           240:                /* W*B5 */
        !           241:                d = fpu_mul(fe);
        !           242:                CPYFPN(&fe->fe_f1, d);
        !           243:                fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logB3);
        !           244:                /* B3+W*B5 */
        !           245:                d = fpu_add(fe);
        !           246:                CPYFPN(&fe->fe_f1, d);
        !           247:                CPYFPN(&fe->fe_f2, &W);
        !           248:                /* W*(B3+W*B5) */
        !           249:                d = fpu_mul(fe);
        !           250:                CPYFPN(&fe->fe_f1, d);
        !           251:                fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logB1);
        !           252:                /* B1+W*(B3+W*B5) */
        !           253:                d = fpu_add(fe);
        !           254:                CPYFPN(&X, d);
        !           255: 
        !           256:                /* [V*(B2+W*B4)] part */
        !           257:                CPYFPN(&fe->fe_f1, &W);
        !           258:                fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logB4);
        !           259:                /* W*B4 */
        !           260:                d = fpu_mul(fe);
        !           261:                CPYFPN(&fe->fe_f1, d);
        !           262:                fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logB2);
        !           263:                /* B2+W*B4 */
        !           264:                d = fpu_add(fe);
        !           265:                CPYFPN(&fe->fe_f1, d);
        !           266:                CPYFPN(&fe->fe_f2, &V);
        !           267:                /* V*(B2+W*B4) */
        !           268:                d = fpu_mul(fe);
        !           269:                CPYFPN(&fe->fe_f1, d);
        !           270:                CPYFPN(&fe->fe_f2, &X);
        !           271:                /* B1+W*(B3+W*B5)+V*(B2+W*B4) */
        !           272:                d = fpu_add(fe);
        !           273:                CPYFPN(&fe->fe_f1, d);
        !           274:                CPYFPN(&fe->fe_f2, &V);
        !           275:                /* V*(B1+W*(B3+W*B5)+V*(B2+W*B4)) */
        !           276:                d = fpu_mul(fe);
        !           277:                CPYFPN(&fe->fe_f1, d);
        !           278:                CPYFPN(&fe->fe_f2, &U);
        !           279:                /* U*V*(B1+W*(B3+W*B5)+V*(B2+W*B4)) */
        !           280:                d = fpu_mul(fe);
        !           281:                CPYFPN(&fe->fe_f1, d);
        !           282:                CPYFPN(&fe->fe_f2, &U);
        !           283:                /* U+U*V*(B1+W*(B3+W*B5)+V*(B2+W*B4)) */
        !           284:                d = fpu_add(fe);
        !           285:        } else /* the usual case */ {
        !           286: #if FPE_DEBUG
        !           287:                printf("__fpu_logn: the usual case. X=(%d,%08x,%08x...)\n",
        !           288:                    X.fp_exp, X.fp_mant[0], X.fp_mant[1]);
        !           289: #endif
        !           290: 
        !           291:                k = X.fp_exp;
        !           292:                /* X <- Y */
        !           293:                X.fp_exp = fe->fe_f2.fp_exp = 0;
        !           294: 
        !           295:                /* get the most significant 7 bits of X */
        !           296:                F.fp_class = FPC_NUM;
        !           297:                F.fp_sign = 0;
        !           298:                F.fp_exp = X.fp_exp;
        !           299:                F.fp_mant[0] = X.fp_mant[0] & (0xfe000000U >> (31 - FP_LG));
        !           300:                F.fp_mant[0] |= (0x01000000U >> (31 - FP_LG));
        !           301:                F.fp_mant[1] = F.fp_mant[2] = 0;
        !           302:                F.fp_sticky = 0;
        !           303: 
        !           304: #if FPE_DEBUG
        !           305:                printf("__fpu_logn: X=Y*2^k=(%d,%08x,%08x...)*2^%d\n",
        !           306:                    fe->fe_f2.fp_exp, fe->fe_f2.fp_mant[0],
        !           307:                    fe->fe_f2.fp_mant[1], k);
        !           308:                printf("__fpu_logn: F=(%d,%08x,%08x...)\n",
        !           309:                    F.fp_exp, F.fp_mant[0], F.fp_mant[1]);
        !           310: #endif
        !           311: 
        !           312:                /* index to the table */
        !           313:                i = (F.fp_mant[0] >> (FP_LG - 7)) & 0x7e;
        !           314: 
        !           315: #if FPE_DEBUG
        !           316:                printf("__fpu_logn: index to logtbl i=%d(%x)\n", i, i);
        !           317: #endif
        !           318: 
        !           319:                CPYFPN(&fe->fe_f1, &F);
        !           320:                /* -F */
        !           321:                fe->fe_f1.fp_sign = 1;
        !           322:                /* Y-F */
        !           323:                d = fpu_add(fe);
        !           324:                CPYFPN(&fe->fe_f1, d);
        !           325: 
        !           326:                /* fe_f2 = 1/F */
        !           327:                fe->fe_f2.fp_class = FPC_NUM;
        !           328:                fe->fe_f2.fp_sign = fe->fe_f2.fp_sticky = fe->fe_f2.fp_mant[2]
        !           329:                    = 0;
        !           330:                fe->fe_f2.fp_exp = logtbl[i].sp_exp;
        !           331:                fe->fe_f2.fp_mant[0] = (logtbl[i].sp_m0 >> (31 - FP_LG));
        !           332:                fe->fe_f2.fp_mant[1] = (logtbl[i].sp_m0 << (FP_LG + 1)) |
        !           333:                    (logtbl[i].sp_m1 >> (31 - FP_LG));
        !           334:                fe->fe_f2.fp_mant[2] =
        !           335:                        (uint32_t)(logtbl[i].sp_m1 << (FP_LG + 1));
        !           336: 
        !           337: #if FPE_DEBUG
        !           338:                printf("__fpu_logn: 1/F=(%d,%08x,%08x...)\n", fe->fe_f2.fp_exp,
        !           339:                    fe->fe_f2.fp_mant[0], fe->fe_f2.fp_mant[1]);
        !           340: #endif
        !           341: 
        !           342:                /* U = (Y-F) * (1/F) */
        !           343:                d = fpu_mul(fe);
        !           344:                CPYFPN(&U, d);
        !           345: 
        !           346:                /* KLOG2 = K * ln(2) */
        !           347:                /* fe_f1 == (fpn)k */
        !           348:                fpu_explode(fe, &fe->fe_f1, FTYPE_LNG, &k);
        !           349:                (void)fpu_const(&fe->fe_f2, FPU_CONST_LN_2);
        !           350: #if FPE_DEBUG
        !           351:                printf("__fpu_logn: fp(k)=(%d,%08x,%08x...)\n",
        !           352:                    fe->fe_f1.fp_exp,
        !           353:                    fe->fe_f1.fp_mant[0], fe->fe_f1.fp_mant[1]);
        !           354:                printf("__fpu_logn: ln(2)=(%d,%08x,%08x...)\n",
        !           355:                    fe->fe_f2.fp_exp,
        !           356:                    fe->fe_f2.fp_mant[0], fe->fe_f2.fp_mant[1]);
        !           357: #endif
        !           358:                /* K * LOGOF2 */
        !           359:                d = fpu_mul(fe);
        !           360:                CPYFPN(&KLOG2, d);
        !           361: 
        !           362:                /* V=U*U */
        !           363:                CPYFPN(&fe->fe_f1, &U);
        !           364:                CPYFPN(&fe->fe_f2, &U);
        !           365:                d = fpu_mul(fe);
        !           366:                CPYFPN(&V, d);
        !           367: 
        !           368:                /*
        !           369:                 * approximation of LOG(1+U) by
        !           370:                 * (U+V*(A1+V*(A3+V*A5)))+(U*V*(A2+V*(A4+V*A6)))
        !           371:                 */
        !           372: 
        !           373:                /* (U+V*(A1+V*(A3+V*A5))) part */
        !           374:                CPYFPN(&fe->fe_f1, d);
        !           375:                fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logA5);
        !           376:                /* V*A5 */
        !           377:                d = fpu_mul(fe);
        !           378: 
        !           379:                CPYFPN(&fe->fe_f1, d);
        !           380:                fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logA3);
        !           381:                /* A3+V*A5 */
        !           382:                d = fpu_add(fe);
        !           383: 
        !           384:                CPYFPN(&fe->fe_f1, d);
        !           385:                CPYFPN(&fe->fe_f2, &V);
        !           386:                /* V*(A3+V*A5) */
        !           387:                d = fpu_mul(fe);
        !           388: 
        !           389:                CPYFPN(&fe->fe_f1, d);
        !           390:                fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logA1);
        !           391:                /* A1+V*(A3+V*A5) */
        !           392:                d = fpu_add(fe);
        !           393: 
        !           394:                CPYFPN(&fe->fe_f1, d);
        !           395:                CPYFPN(&fe->fe_f2, &V);
        !           396:                /* V*(A1+V*(A3+V*A5)) */
        !           397:                d = fpu_mul(fe);
        !           398: 
        !           399:                CPYFPN(&fe->fe_f1, d);
        !           400:                CPYFPN(&fe->fe_f2, &U);
        !           401:                /* U+V*(A1+V*(A3+V*A5)) */
        !           402:                d = fpu_add(fe);
        !           403: 
        !           404:                CPYFPN(&X, d);
        !           405: 
        !           406:                /* (U*V*(A2+V*(A4+V*A6))) part */
        !           407:                CPYFPN(&fe->fe_f1, &V);
        !           408:                fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logA6);
        !           409:                /* V*A6 */
        !           410:                d = fpu_mul(fe);
        !           411:                CPYFPN(&fe->fe_f1, d);
        !           412:                fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logA4);
        !           413:                /* A4+V*A6 */
        !           414:                d = fpu_add(fe);
        !           415:                CPYFPN(&fe->fe_f1, d);
        !           416:                CPYFPN(&fe->fe_f2, &V);
        !           417:                /* V*(A4+V*A6) */
        !           418:                d = fpu_mul(fe);
        !           419:                CPYFPN(&fe->fe_f1, d);
        !           420:                fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logA2);
        !           421:                /* A2+V*(A4+V*A6) */
        !           422:                d = fpu_add(fe);
        !           423:                CPYFPN(&fe->fe_f1, d);
        !           424:                CPYFPN(&fe->fe_f2, &V);
        !           425:                /* V*(A2+V*(A4+V*A6)) */
        !           426:                d = fpu_mul(fe);
        !           427:                CPYFPN(&fe->fe_f1, d);
        !           428:                CPYFPN(&fe->fe_f2, &U);
        !           429:                /* U*V*(A2+V*(A4+V*A6)) */
        !           430:                d = fpu_mul(fe);
        !           431:                CPYFPN(&fe->fe_f1, d);
        !           432:                i++;
        !           433:                /* fe_f2 = logtbl[i+1] (== LOG(F)) */
        !           434:                fe->fe_f2.fp_class = FPC_NUM;
        !           435:                fe->fe_f2.fp_sign = fe->fe_f2.fp_sticky = fe->fe_f2.fp_mant[2]
        !           436:                    = 0;
        !           437:                fe->fe_f2.fp_exp = logtbl[i].sp_exp;
        !           438:                fe->fe_f2.fp_mant[0] = (logtbl[i].sp_m0 >> (31 - FP_LG));
        !           439:                fe->fe_f2.fp_mant[1] = (logtbl[i].sp_m0 << (FP_LG + 1)) |
        !           440:                    (logtbl[i].sp_m1 >> (31 - FP_LG));
        !           441:                fe->fe_f2.fp_mant[2] = (logtbl[i].sp_m1 << (FP_LG + 1));
        !           442: 
        !           443: #if FPE_DEBUG
        !           444:                printf("__fpu_logn: ln(F)=(%d,%08x,%08x,...)\n",
        !           445:                    fe->fe_f2.fp_exp,
        !           446:                    fe->fe_f2.fp_mant[0], fe->fe_f2.fp_mant[1]);
        !           447: #endif
        !           448: 
        !           449:                /* LOG(F)+U*V*(A2+V*(A4+V*A6)) */
        !           450:                d = fpu_add(fe);
        !           451:                CPYFPN(&fe->fe_f1, d);
        !           452:                CPYFPN(&fe->fe_f2, &X);
        !           453:                /* LOG(F)+U+V*(A1+V*(A3+V*A5))+U*V*(A2+V*(A4+V*A6)) */
        !           454:                d = fpu_add(fe);
        !           455: 
        !           456: #if FPE_DEBUG
        !           457:                printf("__fpu_logn: ln(Y)=(%c,%d,%08x,%08x,%08x)\n",
        !           458:                    d->fp_sign ? '-' : '+', d->fp_exp,
        !           459:                    d->fp_mant[0], d->fp_mant[1], d->fp_mant[2]);
        !           460: #endif
        !           461: 
        !           462:                CPYFPN(&fe->fe_f1, d);
        !           463:                CPYFPN(&fe->fe_f2, &KLOG2);
        !           464:                /* K*LOGOF2+LOG(F)+U+V*(A1+V*(A3+V*A5))+U*V*(A2+V*(A4+V*A6)) */
        !           465:                d = fpu_add(fe);
        !           466:        }
        !           467: 
        !           468:        return d;
        !           469: }
        !           470: 
        !           471: struct fpn *
        !           472: fpu_log10(struct fpemu *fe)
        !           473: {
        !           474:        struct fpn *fp = &fe->fe_f2;
        !           475:        uint32_t fpsr;
        !           476: 
        !           477:        fpsr = fe->fe_fpsr & ~FPSR_EXCP;        /* clear all exceptions */
        !           478: 
        !           479:        if (fp->fp_class >= FPC_NUM) {
        !           480:                if (fp->fp_sign) {      /* negative number or Inf */
        !           481:                        fp = fpu_newnan(fe);
        !           482:                        fpsr |= FPSR_OPERR;
        !           483:                } else if (fp->fp_class == FPC_NUM) {
        !           484:                        /* the real work here */
        !           485:                        fp = __fpu_logn(fe);
        !           486:                        if (fp != &fe->fe_f1)
        !           487:                                CPYFPN(&fe->fe_f1, fp);
        !           488:                        (void)fpu_const(&fe->fe_f2, FPU_CONST_LN_10);
        !           489:                        fp = fpu_div(fe);
        !           490:                } /* else if fp == +Inf, return +Inf */
        !           491:        } else if (fp->fp_class == FPC_ZERO) {
        !           492:                /* return -Inf */
        !           493:                fp->fp_class = FPC_INF;
        !           494:                fp->fp_sign = 1;
        !           495:                fpsr |= FPSR_DZ;
        !           496:        } else if (fp->fp_class == FPC_SNAN) {
        !           497:                fpsr |= FPSR_SNAN;
        !           498:                fp = fpu_newnan(fe);
        !           499:        } else {
        !           500:                fp = fpu_newnan(fe);
        !           501:        }
        !           502: 
        !           503:        fe->fe_fpsr = fpsr;
        !           504: 
        !           505:        return fp;
        !           506: }
        !           507: 
        !           508: struct fpn *
        !           509: fpu_log2(struct fpemu *fe)
        !           510: {
        !           511:        struct fpn *fp = &fe->fe_f2;
        !           512:        uint32_t fpsr;
        !           513: 
        !           514:        fpsr = fe->fe_fpsr & ~FPSR_EXCP;        /* clear all exceptions */
        !           515: 
        !           516:        if (fp->fp_class >= FPC_NUM) {
        !           517:                if (fp->fp_sign) {      /* negative number or Inf */
        !           518:                        fp = fpu_newnan(fe);
        !           519:                        fpsr |= FPSR_OPERR;
        !           520:                } else if (fp->fp_class == FPC_NUM) {
        !           521:                        /* the real work here */
        !           522:                        if (fp->fp_mant[0] == FP_1 && fp->fp_mant[1] == 0 &&
        !           523:                            fp->fp_mant[2] == 0) {
        !           524:                                /* fp == 2.0 ^ exp <--> log2(fp) == exp */
        !           525:                                fpu_explode(fe, &fe->fe_f3, FTYPE_LNG,
        !           526:                                    &fp->fp_exp);
        !           527:                                fp = &fe->fe_f3;
        !           528:                        } else {
        !           529:                                fp = __fpu_logn(fe);
        !           530:                                if (fp != &fe->fe_f1)
        !           531:                                        CPYFPN(&fe->fe_f1, fp);
        !           532:                                (void)fpu_const(&fe->fe_f2, FPU_CONST_LN_2);
        !           533:                                fp = fpu_div(fe);
        !           534:                        }
        !           535:                } /* else if fp == +Inf, return +Inf */
        !           536:        } else if (fp->fp_class == FPC_ZERO) {
        !           537:                /* return -Inf */
        !           538:                fp->fp_class = FPC_INF;
        !           539:                fp->fp_sign = 1;
        !           540:                fpsr |= FPSR_DZ;
        !           541:        } else if (fp->fp_class == FPC_SNAN) {
        !           542:                fpsr |= FPSR_SNAN;
        !           543:                fp = fpu_newnan(fe);
        !           544:        } else {
        !           545:                fp = fpu_newnan(fe);
        !           546:        }
        !           547: 
        !           548:        fe->fe_fpsr = fpsr;
        !           549:        return fp;
        !           550: }
        !           551: 
        !           552: struct fpn *
        !           553: fpu_logn(struct fpemu *fe)
        !           554: {
        !           555:        struct fpn *fp = &fe->fe_f2;
        !           556:        uint32_t fpsr;
        !           557: 
        !           558:        fpsr = fe->fe_fpsr & ~FPSR_EXCP;        /* clear all exceptions */
        !           559: 
        !           560:        if (fp->fp_class >= FPC_NUM) {
        !           561:                if (fp->fp_sign) {      /* negative number or Inf */
        !           562:                        fp = fpu_newnan(fe);
        !           563:                        fpsr |= FPSR_OPERR;
        !           564:                } else if (fp->fp_class == FPC_NUM) {
        !           565:                        /* the real work here */
        !           566:                        fp = __fpu_logn(fe);
        !           567:                } /* else if fp == +Inf, return +Inf */
        !           568:        } else if (fp->fp_class == FPC_ZERO) {
        !           569:                /* return -Inf */
        !           570:                fp->fp_class = FPC_INF;
        !           571:                fp->fp_sign = 1;
        !           572:                fpsr |= FPSR_DZ;
        !           573:        } else if (fp->fp_class == FPC_SNAN) {
        !           574:                fpsr |= FPSR_SNAN;
        !           575:                fp = fpu_newnan(fe);
        !           576:        } else {
        !           577:                fp = fpu_newnan(fe);
        !           578:        }
        !           579: 
        !           580:        fe->fe_fpsr = fpsr;
        !           581: 
        !           582:        return fp;
        !           583: }
        !           584: 
        !           585: struct fpn *
        !           586: fpu_lognp1(struct fpemu *fe)
        !           587: {
        !           588:        struct fpn *fp;
        !           589: 
        !           590:        /* if src is +0/-0, return +0/-0 */
        !           591:        if (ISZERO(&fe->fe_f2))
        !           592:                return &fe->fe_f2;
        !           593: 
        !           594:        /* build a 1.0 */
        !           595:        fp = fpu_const(&fe->fe_f1, FPU_CONST_1);
        !           596:        /* fp = 1.0 + f2 */
        !           597:        fp = fpu_add(fe);
        !           598: 
        !           599:        /* copy the result to the src opr */
        !           600:        CPYFPN(&fe->fe_f2, fp);
        !           601: 
        !           602:        return fpu_logn(fe);
        !           603: }

unix.superglobalmegacorp.com

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