Subversion Repositories HelenOS-historic

Rev

Rev 534 | Go to most recent revision | 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;
625 palkovsky 65 
    __u64 mask;
66 
    // __u64 mantisa;
230 cejka 67 
    int i;
68 
 
69 
    exp=fmath_get_binary_exponent(num);
70 
 
71 
    if (exp<0) {
72 
        *intp = 0.0;
73 
        return num;
74 
        }
75 
 
76 
 
77 
    if (exp>51) {
78 
        *intp=num;
79 
        num=0.0;
80 
        return num;
81 
    }
82 
 
83 
    fmath_ld_union_num.bf = num;
84 
 
85 
    mask = FMATH_MANTISA_MASK>>exp;
86 
    //mantisa = (fmath_get-binary_mantisa(num))&(~mask);
87 
 
88 
    for (i=0;i<7;i++) {
89 
        /* Ugly construction for obtain sign, exponent and integer part from num */
90 
        fmath_ld_union_int.ldd[i]=fmath_ld_union_num.ldd[i]&(((~mask)>>(i*8))&0xff);
91 
    }
92 
 
93 
    fmath_ld_union_int.ldd[6]|=((fmath_ld_union_num.ldd[6])&(0xf0));
94 
    fmath_ld_union_int.ldd[7]=fmath_ld_union_num.ldd[7];
95 
 
96 
    *intp=fmath_ld_union_int.bf;
97 
    return fmath_ld_union_num.bf-fmath_ld_union_int.bf;
98 
};
99 
 
100 
 
101 
double fmath_dpow(double base, double exponent)
102 
{
103 
    double value=1.0;
104 
    if (base<=0.0) return base;
105 
 
106 
    //2^(x*log2(10)) = 2^y = 10^x
107 
 
108 
    __asm__ __volatile__ (      \
109 
        "fyl2x # ST(1):=ST(1)*log2(ST(0)), pop st(0) \n\t "     \
110 
        "fld    %%st(0) \n\t"   \
111 
        "frndint \n\t"      \
112 
        "fxch %%st(1) \n\t"     \
113 
        "fsub %%st(1),%%st(0) \n\t" \
114 
        "f2xm1  # ST := 2^ST -1\n\t"            \
115 
        "fld1 \n\t"         \
116 
        "faddp %%st(0),%%st(1) \n\t"    \
117 
        "fscale #ST:=ST*2^(ST(1))\n\t"      \
118 
        "fstp %%st(1) \n\t"     \
119 
    "" : "=t" (value) :  "0" (base), "u" (exponent) );
120 
    return value;
121 
}
122 
 
266 cejka 123 
int fmath_is_nan(double num)
124 
{
125 
    __u16 exp;
126 
    fmath_ld_union_t fmath_ld_union;
127 
    fmath_ld_union.bf = num;
128 
    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 */
129 
 
130 
    if (exp!=0x07ff) return 0;
131 
    if (fmath_get_binary_mantisa(num)>=FMATH_NAN) return 1;
132 
 
133 
 
134 
    return 0;
135 
}
136 
 
137 
int fmath_is_infinity(double num)
138 
{
139 
    __u16 exp;
140 
    fmath_ld_union_t fmath_ld_union;
141 
    fmath_ld_union.bf = num;
142 
    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 */
143 
 
144 
    if (exp!=0x07ff) return 0;
145 
    if (fmath_get_binary_mantisa(num)==0x0) return 1;
146 
    return 0;
147 
}