Imported Upstream version 2.9.0
[debian/cc1111] / device / lib / sincoshf.c
1 /*  sincoshf.c: Computes sinh or cosh of a 32-bit float as outlined in [1]
2
3     Copyright (C) 2001, 2002  Jesus Calvino-Fraga, jesusc@ieee.org
4
5     This library is free software; you can redistribute it and/or
6     modify it under the terms of the GNU Lesser General Public
7     License as published by the Free Software Foundation; either
8     version 2.1 of the License, or (at your option) any later version.
9
10     This library is distributed in the hope that it will be useful,
11     but WITHOUT ANY WARRANTY; without even the implied warranty of
12     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
13     Lesser General Public License for more details.
14
15     You should have received a copy of the GNU Lesser General Public
16     License along with this library; if not, write to the Free Software
17     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA */
18
19 /* [1] William James Cody and W.  M.  Waite.  _Software manual for the
20    elementary functions_, Englewood Cliffs, N.J.:Prentice-Hall, 1980. */
21
22 /* Version 1.0 - Initial release */
23
24 #include <math.h>
25 #include <errno.h>
26 #include <stdbool.h>
27
28 #define P0 -0.713793159E+1
29 #define P1 -0.190333999E+0
30 #define Q0 -0.428277109E+2
31 #define Q1  0.100000000E+1
32
33 #define P(z) (P1*z+P0)
34 #define Q(z) (Q1*z+Q0)
35
36 #define K1 0.69316101074218750000E+0 /* ln(v)   */
37 #define K2 0.24999308500451499336E+0 /* v**(-2) */
38 #define K3 0.13830277879601902638E-4 /* v/2-1   */
39
40 //WMAX is defined as ln(HUGE_VALF)-ln(v)+0.69
41 #define WMAX 44.93535952E+0
42 //WBAR 0.35*(b+1)
43 #define WBAR 1.05
44 #define YBAR 9.0 /*Works for me*/
45
46 float sincoshf(const float x, const int iscosh)
47 {
48     float y, w, z;
49         BOOL sign;
50
51     if (x<0.0) { y=-x; sign=1; }
52           else { y=x;  sign=0; }
53
54     if ((y>1.0) || iscosh)
55     {
56         if(y>YBAR)
57         {
58             w=y-K1;
59             if (w>WMAX)
60             {
61                 errno=ERANGE;
62                 z=HUGE_VALF;
63             }
64             else
65             {
66                 z=expf(w);
67                 z+=K3*z;
68             }
69         }
70         else
71         {
72             z=expf(y);
73             w=1.0/z;
74             if(!iscosh) w=-w;
75             z=(z+w)*0.5;
76         }
77         if(sign) z=-z;
78     }
79     else
80     {
81         if (y<EPS)
82             z=x;
83         else
84         {
85             z=x*x;
86             z=x+x*z*P(z)/Q(z);
87         }
88     }
89     return z;
90 }