git.gag.com Git - debian/cc1111/blob - device/lib/pic/libm/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 /*
  25 ** $Id: sincosf.c 4776 2007-04-29 13:15:51Z borutr $
  26 */
  27
  28 #include <math.h>
  29 #include <errno.h>
  30
  31 #define r1      (-0.1666665668E+0)
  32 #define r2      (0.8333025139E-2)
  33 #define r3      (-0.1980741872E-3)
  34 #define r4       (0.2601903036E-5)
  35
  36 /* PI=C1+C2 */
  37 #define C1       3.140625
  38 #define C2       9.676535897E-4
  39
  40 /*A reasonable value for YMAX is the int part of PI*B**(t/2)=3.1416*2**(12)*/
  41 #define YMAX     12867.0
  42
  43 float sincosf(float x, int iscos)
  44 {
  45     float y, f, r, g, XN;
  46     int N;
  47     char sign;
  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 }