Subversion Repositories HelenOS

Rev

Rev 1962 | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
230 cejka 1 
/*
2 
 * Copyright (C) 2005 Josef Cejka
3 
 * All rights reserved.
4 
 *
5 
 * Redistribution and use in source and binary forms, with or without
6 
 * modification, are permitted provided that the following conditions
7 
 * are met:
8 
 *
9 
 * - Redistributions of source code must retain the above copyright
10 
 *   notice, this list of conditions and the following disclaimer.
11 
 * - Redistributions in binary form must reproduce the above copyright
12 
 *   notice, this list of conditions and the following disclaimer in the
13 
 *   documentation and/or other materials provided with the distribution.
14 
 * - The name of the author may not be used to endorse or promote products
15 
 *   derived from this software without specific prior written permission.
16 
 *
17 
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
18 
 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
19 
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
20 
 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
21 
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
22 
 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
23 
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
24 
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
25 
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
26 
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27 
 */
28 
 
29 
#include <arch/fmath.h>
30 
#include <print.h>
31 
 
32 
#define FMATH_MANTISA_MASK ( 0x000fffffffffffffLL )
266 cejka 33 
#define FMATH_NAN ( 0x0001000000000001LL )
230 cejka 34 
signed short fmath_get_binary_exponent(double num)
35 
{
36 
    fmath_ld_union_t fmath_ld_union;
37 
    fmath_ld_union.bf = num;
38 
    return (signed short)((((fmath_ld_union.ldd[7])&0x7f)<<4) + (((fmath_ld_union.ldd[6])&0xf0)>>4)) -FMATH_EXPONENT_BIAS; /* exponent is 11 bits lenght, so sevent bits is in 8th byte and 4 bits in 7th */
39 
}
40 
 
41 
double fmath_get_decimal_exponent(double num)
42 
{
43 
    double value;
44 
    /* log10(2)*log2(x) => log10(x) */
45 
    __asm__ __volatile__ ( \
46 
    "fldlg2     #load log10(2)  \n\t"   \
47 
    "fxch %%st(1)       \n\t" \
48 
    "fyl2x      #count st(0)*log2(st(1))->st(1); pop st(0)  \n\t" \
49 
    : "=t" (value) : "0"(num) );
50 
    return value;
51 
}
52 
 
53 
__u64 fmath_get_binary_mantisa(double num)
54 
{
55 
    union { __u64 _u; double _d;} un = { _d : num };
56 
    un._u=un._u &(FMATH_MANTISA_MASK); /* mask 52 bits of mantisa*/
57 
    return un._u;
58 
}
59 
 
60 
double fmath_fint(double num, double *intp)
61 
{
62 
    fmath_ld_union_t fmath_ld_union_num;
63 
    fmath_ld_union_t fmath_ld_union_int;
64 
    signed short exp;
65 
    __u64 mask,mantisa;
66 
    int i;
67 
 
68 
    exp=fmath_get_binary_exponent(num);
69 
 
70 
    if (exp<0) {
71 
        *intp = 0.0;
72 
        return num;
73 
        }
74 
 
75 
 
76 
    if (exp>51) {
77 
        *intp=num;
78 
        num=0.0;
79 
        return num;
80 
    }
81 
 
82 
    fmath_ld_union_num.bf = num;
83 
 
84 
    mask = FMATH_MANTISA_MASK>>exp;
85 
    //mantisa = (fmath_get-binary_mantisa(num))&(~mask);
86 
 
87 
    for (i=0;i<7;i++) {
88 
        /* Ugly construction for obtain sign, exponent and integer part from num */
89 
        fmath_ld_union_int.ldd[i]=fmath_ld_union_num.ldd[i]&(((~mask)>>(i*8))&0xff);
90 
    }
91 
 
92 
    fmath_ld_union_int.ldd[6]|=((fmath_ld_union_num.ldd[6])&(0xf0));
93 
    fmath_ld_union_int.ldd[7]=fmath_ld_union_num.ldd[7];
94 
 
95 
    *intp=fmath_ld_union_int.bf;
96 
    return fmath_ld_union_num.bf-fmath_ld_union_int.bf;
97 
};
98 
 
99 
 
100 
double fmath_dpow(double base, double exponent)
101 
{
102 
    double value=1.0;
103 
    if (base<=0.0) return base;
104 
 
105 
    //2^(x*log2(10)) = 2^y = 10^x
106 
 
107 
    __asm__ __volatile__ (      \
108 
        "fyl2x # ST(1):=ST(1)*log2(ST(0)), pop st(0) \n\t "     \
109 
        "fld    %%st(0) \n\t"   \
110 
        "frndint \n\t"      \
111 
        "fxch %%st(1) \n\t"     \
112 
        "fsub %%st(1),%%st(0) \n\t" \
113 
        "f2xm1  # ST := 2^ST -1\n\t"            \
114 
        "fld1 \n\t"         \
115 
        "faddp %%st(0),%%st(1) \n\t"    \
116 
        "fscale #ST:=ST*2^(ST(1))\n\t"      \
117 
        "fstp %%st(1) \n\t"     \
118 
    "" : "=t" (value) :  "0" (base), "u" (exponent) );
119 
    return value;
120 
}
121 
 
266 cejka 122 
int fmath_is_nan(double num)
123 
{
124 
    __u16 exp;
125 
    fmath_ld_union_t fmath_ld_union;
126 
    fmath_ld_union.bf = num;
127 
    exp=(((fmath_ld_union.ldd[7])&0x7f)<<4) + (((fmath_ld_union.ldd[6])&0xf0)>>4); /* exponent is 11 bits lenght, so sevent bits is in 8th byte and 4 bits in 7th */
128 
 
129 
    if (exp!=0x07ff) return 0;
130 
    if (fmath_get_binary_mantisa(num)>=FMATH_NAN) return 1;
131 
 
132 
 
133 
    return 0;
134 
}
135 
 
136 
int fmath_is_infinity(double num)
137 
{
138 
    __u16 exp;
139 
    fmath_ld_union_t fmath_ld_union;
140 
    fmath_ld_union.bf = num;
141 
    exp=(((fmath_ld_union.ldd[7])&0x7f)<<4) + (((fmath_ld_union.ldd[6])&0xf0)>>4); /* exponent is 11 bits lenght, so sevent bits is in 8th byte and 4 bits in 7th */
142 
 
143 
    if (exp!=0x07ff) return 0;
144 
    if (fmath_get_binary_mantisa(num)==0x0) return 1;
145 
    return 0;
146 
}