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