← Back to branch summary

~ubuntu-branches/ubuntu/natty/eglibc/natty-security

« back to all changes in this revision

Viewing changes to sysdeps/i386/fpu/s_cexp.S

  • Committer: Bazaar Package Importer
  • Author(s): Aurelien Jarno
  • Date: 2009-05-05 09:54:14 UTC
  • Revision ID: james.westby@ubuntu.com-20090505095414-c45qsg9ixjheohru
Import upstream version 2.9

Show diffs side-by-side

added added

removed removed

Lines of Context:
sysdeps/i386/fpu/s_cexp.S
 
1
/* ix87 specific implementation of complex exponential function for double.
 
2
   Copyright (C) 1997, 2005 Free Software Foundation, Inc.
 
3
   This file is part of the GNU C Library.
 
4
   Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
 
5
 
 
6
   The GNU C Library is free software; you can redistribute it and/or
 
7
   modify it under the terms of the GNU Lesser General Public
 
8
   License as published by the Free Software Foundation; either
 
9
   version 2.1 of the License, or (at your option) any later version.
 
10
 
 
11
   The GNU C Library is distributed in the hope that it will be useful,
 
12
   but WITHOUT ANY WARRANTY; without even the implied warranty of
 
13
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 
14
   Lesser General Public License for more details.
 
15
 
 
16
   You should have received a copy of the GNU Lesser General Public
 
17
   License along with the GNU C Library; if not, write to the Free
 
18
   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
 
19
   02111-1307 USA.  */
 
20
 
 
21
#include <sysdep.h>
 
22
 
 
23
#ifdef __ELF__
 
24
        .section .rodata
 
25
#else
 
26
        .text
 
27
#endif
 
28
        .align ALIGNARG(4)
 
29
        ASM_TYPE_DIRECTIVE(huge_nan_null_null,@object)
 
30
huge_nan_null_null:
 
31
        .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f
 
32
        .byte 0, 0, 0, 0, 0, 0, 0xff, 0x7f
 
33
        .double 0.0
 
34
zero:   .double 0.0
 
35
infinity:
 
36
        .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f
 
37
        .byte 0, 0, 0, 0, 0, 0, 0xff, 0x7f
 
38
        .double 0.0
 
39
        .byte 0, 0, 0, 0, 0, 0, 0, 0x80
 
40
        ASM_SIZE_DIRECTIVE(huge_nan_null_null)
 
41
 
 
42
        ASM_TYPE_DIRECTIVE(twopi,@object)
 
43
twopi:
 
44
        .byte 0x35, 0xc2, 0x68, 0x21, 0xa2, 0xda, 0xf, 0xc9, 0x1, 0x40
 
45
        .byte 0, 0, 0, 0, 0, 0
 
46
        ASM_SIZE_DIRECTIVE(twopi)
 
47
 
 
48
        ASM_TYPE_DIRECTIVE(l2e,@object)
 
49
l2e:
 
50
        .byte 0xbc, 0xf0, 0x17, 0x5c, 0x29, 0x3b, 0xaa, 0xb8, 0xff, 0x3f
 
51
        .byte 0, 0, 0, 0, 0, 0
 
52
        ASM_SIZE_DIRECTIVE(l2e)
 
53
 
 
54
        ASM_TYPE_DIRECTIVE(one,@object)
 
55
one:    .double 1.0
 
56
        ASM_SIZE_DIRECTIVE(one)
 
57
 
 
58
 
 
59
#ifdef PIC
 
60
#define MO(op) op##@GOTOFF(%ecx)
 
61
#define MOX(op,x,f) op##@GOTOFF(%ecx,x,f)
 
62
#else
 
63
#define MO(op) op
 
64
#define MOX(op,x,f) op(,x,f)
 
65
#endif
 
66
 
 
67
        .text
 
68
ENTRY(__cexp)
 
69
        fldl    8(%esp)                 /* x */
 
70
        fxam
 
71
        fnstsw
 
72
        fldl    16(%esp)                /* y : x */
 
73
#ifdef  PIC
 
74
        LOAD_PIC_REG (cx)
 
75
#endif
 
76
        movb    %ah, %dh
 
77
        andb    $0x45, %ah
 
78
        cmpb    $0x05, %ah
 
79
        je      1f                      /* Jump if real part is +-Inf */
 
80
        cmpb    $0x01, %ah
 
81
        je      2f                      /* Jump if real part is NaN */
 
82
 
 
83
        fxam                            /* y : x */
 
84
        fnstsw
 
85
        /* If the imaginary part is not finite we return NaN+i NaN, as
 
86
           for the case when the real part is NaN.  A test for +-Inf and
 
87
           NaN would be necessary.  But since we know the stack register
 
88
           we applied `fxam' to is not empty we can simply use one test.
 
89
           Check your FPU manual for more information.  */
 
90
        andb    $0x01, %ah
 
91
        cmpb    $0x01, %ah
 
92
        je      20f
 
93
 
 
94
        /* We have finite numbers in the real and imaginary part.  Do
 
95
           the real work now.  */
 
96
        fxch                    /* x : y */
 
97
        fldt    MO(l2e)         /* log2(e) : x : y */
 
98
        fmulp                   /* x * log2(e) : y */
 
99
        fld     %st             /* x * log2(e) : x * log2(e) : y */
 
100
        frndint                 /* int(x * log2(e)) : x * log2(e) : y */
 
101
        fsubr   %st, %st(1)     /* int(x * log2(e)) : frac(x * log2(e)) : y */
 
102
        fxch                    /* frac(x * log2(e)) : int(x * log2(e)) : y */
 
103
        f2xm1                   /* 2^frac(x * log2(e))-1 : int(x * log2(e)) : y */
 
104
        faddl   MO(one)         /* 2^frac(x * log2(e)) : int(x * log2(e)) : y */
 
105
        fscale                  /* e^x : int(x * log2(e)) : y */
 
106
        fst     %st(1)          /* e^x : e^x : y */
 
107
        fxch    %st(2)          /* y : e^x : e^x */
 
108
        fsincos                 /* cos(y) : sin(y) : e^x : e^x */
 
109
        fnstsw
 
110
        testl   $0x400, %eax
 
111
        jnz     7f
 
112
        fmulp   %st, %st(3)     /* sin(y) : e^x : e^x * cos(y) */
 
113
        fmulp   %st, %st(1)     /* e^x * sin(y) : e^x * cos(y) */
 
114
        movl    4(%esp), %eax           /* Pointer to memory for result.  */
 
115
        fstpl   8(%eax)
 
116
        fstpl   (%eax)
 
117
        ret     $4
 
118
 
 
119
        /* We have to reduce the argument to fsincos.  */
 
120
        .align ALIGNARG(4)
 
121
7:      fldt    MO(twopi)       /* 2*pi : y : e^x : e^x */
 
122
        fxch                    /* y : 2*pi : e^x : e^x */
 
123
8:      fprem1                  /* y%(2*pi) : 2*pi : e^x : e^x */
 
124
        fnstsw
 
125
        testl   $0x400, %eax
 
126
        jnz     8b
 
127
        fstp    %st(1)          /* y%(2*pi) : e^x : e^x */
 
128
        fsincos                 /* cos(y) : sin(y) : e^x : e^x */
 
129
        fmulp   %st, %st(3)
 
130
        fmulp   %st, %st(1)
 
131
        movl    4(%esp), %eax           /* Pointer to memory for result.  */
 
132
        fstpl   8(%eax)
 
133
        fstpl   (%eax)
 
134
        ret     $4
 
135
 
 
136
        /* The real part is +-inf.  We must make further differences.  */
 
137
        .align ALIGNARG(4)
 
138
1:      fxam                    /* y : x */
 
139
        fnstsw
 
140
        movb    %ah, %dl
 
141
        testb   $0x01, %ah      /* See above why 0x01 is usable here.  */
 
142
        jne     3f
 
143
 
 
144
 
 
145
        /* The real part is +-Inf and the imaginary part is finite.  */
 
146
        andl    $0x245, %edx
 
147
        cmpb    $0x40, %dl      /* Imaginary part == 0?  */
 
148
        je      4f              /* Yes ->  */
 
149
 
 
150
        fxch                    /* x : y */
 
151
        shrl    $5, %edx
 
152
        fstp    %st(0)          /* y */ /* Drop the real part.  */
 
153
        andl    $16, %edx       /* This puts the sign bit of the real part
 
154
                                   in bit 4.  So we can use it to index a
 
155
                                   small array to select 0 or Inf.  */
 
156
        fsincos                 /* cos(y) : sin(y) */
 
157
        fnstsw
 
158
        testl   $0x0400, %eax
 
159
        jnz     5f
 
160
        fldl    MOX(huge_nan_null_null,%edx,1)
 
161
        movl    4(%esp), %edx           /* Pointer to memory for result.  */
 
162
        fstl    8(%edx)
 
163
        fstpl   (%edx)
 
164
        ftst
 
165
        fnstsw
 
166
        shll    $23, %eax
 
167
        andl    $0x80000000, %eax
 
168
        orl     %eax, 4(%edx)
 
169
        fstp    %st(0)
 
170
        ftst
 
171
        fnstsw
 
172
        shll    $23, %eax
 
173
        andl    $0x80000000, %eax
 
174
        orl     %eax, 12(%edx)
 
175
        fstp    %st(0)
 
176
        ret     $4
 
177
        /* We must reduce the argument to fsincos.  */
 
178
        .align ALIGNARG(4)
 
179
5:      fldt    MO(twopi)
 
180
        fxch
 
181
6:      fprem1
 
182
        fnstsw
 
183
        testl   $0x400, %eax
 
184
        jnz     6b
 
185
        fstp    %st(1)
 
186
        fsincos
 
187
        fldl    MOX(huge_nan_null_null,%edx,1)
 
188
        movl    4(%esp), %edx           /* Pointer to memory for result.  */
 
189
        fstl    8(%edx)
 
190
        fstpl   (%edx)
 
191
        ftst
 
192
        fnstsw
 
193
        shll    $23, %eax
 
194
        andl    $0x80000000, %eax
 
195
        orl     %eax, 4(%edx)
 
196
        fstp    %st(0)
 
197
        ftst
 
198
        fnstsw
 
199
        shll    $23, %eax
 
200
        andl    $0x80000000, %eax
 
201
        orl     %eax, 12(%edx)
 
202
        fstp    %st(0)
 
203
        ret     $4
 
204
 
 
205
        /* The real part is +-Inf and the imaginary part is +-0.  So return
 
206
           +-Inf+-0i.  */
 
207
        .align ALIGNARG(4)
 
208
4:      movl    4(%esp), %eax           /* Pointer to memory for result.  */
 
209
        fstpl   8(%eax)
 
210
        shrl    $5, %edx
 
211
        fstp    %st(0)
 
212
        andl    $16, %edx
 
213
        fldl    MOX(huge_nan_null_null,%edx,1)
 
214
        fstpl   (%eax)
 
215
        ret     $4
 
216
 
 
217
        /* The real part is +-Inf, the imaginary is also is not finite.  */
 
218
        .align ALIGNARG(4)
 
219
3:      fstp    %st(0)
 
220
        fstp    %st(0)          /* <empty> */
 
221
        andb    $0x45, %ah
 
222
        andb    $0x47, %dh
 
223
        xorb    %dh, %ah
 
224
        jnz     30f
 
225
        fldl    MO(infinity)    /* Raise invalid exception.  */
 
226
        fmull   MO(zero)
 
227
        fstp    %st(0)
 
228
30:     movl    %edx, %eax
 
229
        shrl    $5, %edx
 
230
        shll    $4, %eax
 
231
        andl    $16, %edx
 
232
        andl    $32, %eax
 
233
        orl     %eax, %edx
 
234
        movl    4(%esp), %eax           /* Pointer to memory for result.  */
 
235
 
 
236
        fldl    MOX(huge_nan_null_null,%edx,1)
 
237
        fldl    MOX(huge_nan_null_null+8,%edx,1)
 
238
        fxch
 
239
        fstpl   (%eax)
 
240
        fstpl   8(%eax)
 
241
        ret     $4
 
242
 
 
243
        /* The real part is NaN.  */
 
244
        .align ALIGNARG(4)
 
245
20:     fldl    MO(infinity)            /* Raise invalid exception.  */
 
246
        fmull   MO(zero)
 
247
        fstp    %st(0)
 
248
2:      fstp    %st(0)
 
249
        fstp    %st(0)
 
250
        movl    4(%esp), %eax           /* Pointer to memory for result.  */
 
251
        fldl    MO(huge_nan_null_null+8)
 
252
        fstl    (%eax)
 
253
        fstpl   8(%eax)
 
254
        ret     $4
 
255
 
 
256
END(__cexp)
 
257
weak_alias (__cexp, cexp)