OpenCores
URL https://opencores.org/ocsvn/or1k/or1k/trunk

Subversion Repositories or1k

[/] [or1k/] [trunk/] [linux/] [linux-2.4/] [arch/] [mips/] [math-emu/] [sp_sqrt.c] - Blame information for rev 1765

Details | Compare with Previous | View Log

Line No. Rev Author Line
1 1275 phoenix 
/* IEEE754 floating point arithmetic
2 
 * single precision square root
3 
 */
4 
/*
5 
 * MIPS floating point support
6 
 * Copyright (C) 1994-2000 Algorithmics Ltd.  All rights reserved.
7 
 * http://www.algor.co.uk
8 
 *
9 
 * ########################################################################
10 
 *
11 
 *  This program is free software; you can distribute it and/or modify it
12 
 *  under the terms of the GNU General Public License (Version 2) as
13 
 *  published by the Free Software Foundation.
14 
 *
15 
 *  This program is distributed in the hope it will be useful, but WITHOUT
16 
 *  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 
 *  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
18 
 *  for more details.
19 
 *
20 
 *  You should have received a copy of the GNU General Public License along
21 
 *  with this program; if not, write to the Free Software Foundation, Inc.,
22 
 *  59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 
 *
24 
 * ########################################################################
25 
 */
26 
 
27 
 
28 
#include "ieee754sp.h"
29 
 
30 
ieee754sp ieee754sp_sqrt(ieee754sp x)
31 
{
32 
        int ix, s, q, m, t, i;
33 
        unsigned int r;
34 
        COMPXSP;
35 
 
36 
        /* take care of Inf and NaN */
37 
 
38 
        EXPLODEXSP;
39 
        CLEARCX;
40 
        FLUSHXSP;
41 
 
42 
        /* x == INF or NAN? */
43 
        switch (xc) {
44 
        case IEEE754_CLASS_QNAN:
45 
                /* sqrt(Nan) = Nan */
46 
                return ieee754sp_nanxcpt(x, "sqrt");
47 
        case IEEE754_CLASS_SNAN:
48 
                SETCX(IEEE754_INVALID_OPERATION);
49 
                return ieee754sp_nanxcpt(ieee754sp_indef(), "sqrt");
50 
        case IEEE754_CLASS_ZERO:
51 
                /* sqrt(0) = 0 */
52 
                return x;
53 
        case IEEE754_CLASS_INF:
54 
                if (xs) {
55 
                        /* sqrt(-Inf) = Nan */
56 
                        SETCX(IEEE754_INVALID_OPERATION);
57 
                        return ieee754sp_nanxcpt(ieee754sp_indef(), "sqrt");
58 
                }
59 
                /* sqrt(+Inf) = Inf */
60 
                return x;
61 
        case IEEE754_CLASS_DNORM:
62 
        case IEEE754_CLASS_NORM:
63 
                if (xs) {
64 
                        /* sqrt(-x) = Nan */
65 
                        SETCX(IEEE754_INVALID_OPERATION);
66 
                        return ieee754sp_nanxcpt(ieee754sp_indef(), "sqrt");
67 
                }
68 
                break;
69 
        }
70 
 
71 
        ix = x.bits;
72 
 
73 
        /* normalize x */
74 
        m = (ix >> 23);
75 
        if (m == 0) {            /* subnormal x */
76 
                for (i = 0; (ix & 0x00800000) == 0; i++)
77 
                        ix <<= 1;
78 
                m -= i - 1;
79 
        }
80 
        m -= 127;               /* unbias exponent */
81 
        ix = (ix & 0x007fffff) | 0x00800000;
82 
        if (m & 1)              /* odd m, double x to make it even */
83 
                ix += ix;
84 
        m >>= 1;                /* m = [m/2] */
85 
 
86 
        /* generate sqrt(x) bit by bit */
87 
        ix += ix;
88 
        q = s = 0;               /* q = sqrt(x) */
89 
        r = 0x01000000;         /* r = moving bit from right to left */
90 
 
91 
        while (r != 0) {
92 
                t = s + r;
93 
                if (t <= ix) {
94 
                        s = t + r;
95 
                        ix -= t;
96 
                        q += r;
97 
                }
98 
                ix += ix;
99 
                r >>= 1;
100 
        }
101 
 
102 
        if (ix != 0) {
103 
                SETCX(IEEE754_INEXACT);
104 
                switch (ieee754_csr.rm) {
105 
                case IEEE754_RP:
106 
                        q += 2;
107 
                        break;
108 
                case IEEE754_RN:
109 
                        q += (q & 1);
110 
                        break;
111 
                }
112 
        }
113 
        ix = (q >> 1) + 0x3f000000;
114 
        ix += (m << 23);
115 
        x.bits = ix;
116 
        return x;
117 
}

powered by: WebSVN 2.1.0

© copyright 1999-2024 OpenCores.org, equivalent to Oliscience, all rights reserved. OpenCores®, registered trademark.