git.gag.com Git - debian/cc1111/blob - device/lib/sincosf.c

   1 /*  sincosf.c: Computes sin or cos 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 r1      -0.1666665668E+0
  29 #define r2       0.8333025139E-2
  30 #define r3      -0.1980741872E-3
  31 #define r4       0.2601903036E-5
  32
  33 /* PI=C1+C2 */
  34 #define C1       3.140625
  35 #define C2       9.676535897E-4
  36
  37 /*A reasonable value for YMAX is the int part of PI*B**(t/2)=3.1416*2**(12)*/
  38 #define YMAX     12867.0
  39
  40 float sincosf(const float x, const int iscos)
  41 {
  42     float y, f, r, g, XN;
  43     int N;
  44         BOOL sign;
  45
  46     if(iscos)
  47     {
  48         y=fabsf(x)+HALF_PI;
  49         sign=0;
  50     }
  51     else
  52     {
  53         if(x<0.0)
  54             { y=-x; sign=1; }
  55         else
  56             { y=x; sign=0; }
  57     }
  58
  59     if(y>YMAX)
  60     {
  61         errno=ERANGE;
  62         return 0.0;
  63     }
  64
  65     /*Round y/PI to the nearest integer*/
  66     N=((y*iPI)+0.5); /*y is positive*/
  67
  68     /*If N is odd change sign*/
  69     if(N&1) sign=!sign;
  70
  71     XN=N;
  72     /*Cosine required? (is done here to keep accuracy)*/
  73     if(iscos) XN-=0.5;
  74
  75     y=fabsf(x);
  76     r=(int)y;
  77     g=y-r;
  78     f=((r-XN*C1)+g)-XN*C2;
  79
  80     g=f*f;
  81     if(g>EPS2) //Used to be if(fabsf(f)>EPS)
  82     {
  83         r=(((r4*g+r3)*g+r2)*g+r1)*g;
  84         f+=f*r;
  85     }
  86     return (sign?-f:f);
  87 }