Open64: osprey/libm/mips/fmodf.c Source File

00001 /*
00002  * Copyright 2003, 2004 PathScale, Inc.  All Rights Reserved.
00003  */
00004 
00005 
00006 /*
00007 
00008   Copyright (C) 2000, 2001 Silicon Graphics, Inc.  All Rights Reserved.
00009 
00010   This program is free software; you can redistribute it and/or modify it
00011   under the terms of version 2.1 of the GNU Lesser General Public License 
00012   as published by the Free Software Foundation.
00013 
00014   This program is distributed in the hope that it would be useful, but
00015   WITHOUT ANY WARRANTY; without even the implied warranty of
00016   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
00017 
00018   Further, this software is distributed without any warranty that it is
00019   free of the rightful claim of any third person regarding infringement 
00020   or the like.  Any license provided herein, whether implied or 
00021   otherwise, applies only to this software file.  Patent licenses, if
00022   any, provided herein do not apply to combinations of this program with 
00023   other software, or any other product whatsoever.  
00024 
00025   You should have received a copy of the GNU Lesser General Public 
00026   License along with this program; if not, write the Free Software 
00027   Foundation, Inc., 59 Temple Place - Suite 330, Boston MA 02111-1307, 
00028   USA.
00029 
00030   Contact information:  Silicon Graphics, Inc., 1600 Amphitheatre Pky,
00031   Mountain View, CA 94043, or:
00032 
00033   http://www.sgi.com
00034 
00035   For further information regarding this notice, see:
00036 
00037   http://oss.sgi.com/projects/GenInfo/NoticeExplan
00038 
00039 */
00040 
00041 /* ====================================================================
00042  * ====================================================================
00043  *
00044  * Module: fmodf.c
00045  * $Revision: 1.1.1.1 $
00046  * $Date: 2005/10/21 19:00:00 $
00047  * $Author: marcel $
00048  * $Source: /proj/osprey/CVS/open64/osprey1.0/libm/mips/fmodf.c,v $
00049  *
00050  * Revision history:
00051  *  09-Jun-93 - Original Version
00052  *
00053  * Description: source code for fmodf function
00054  *
00055  * ====================================================================
00056  * ====================================================================
00057  */
00058 
00059 static char *rcs_id = "$Source: /proj/osprey/CVS/open64/osprey1.0/libm/mips/fmodf.c,v $ $Revision: 1.1.1.1 $";
00060 
00061 #ifdef _CALL_MATHERR
00062 #include <stdio.h>
00063 #include <math.h>
00064 #include <errno.h>
00065 #endif
00066 
00067 #include "libm.h"
00068 
00069 #if defined(mips) && !defined(__GNUC__)
00070 
00071 #include <ieeefp.h>
00072 
00073 extern  fp_rnd  swapRM(fp_rnd);
00074 #endif
00075 
00076 #ifdef __GNUC__
00077 
00078 #include <fenv.h>
00079 
00080 extern  int swapRM(int);
00081 #endif
00082 
00083 #if defined(mips) && !defined(__GNUC__)
00084 extern  float fmodf(float, float);
00085 
00086 #pragma weak fmodf = __fmodf
00087 #endif
00088 
00089 #if defined(BUILD_OS_DARWIN) /* Mach-O doesn't support aliases */
00090 extern float __fmodf(float, float);
00091 #pragma weak fmodf
00092 float fmodf( float x, float y ) {
00093   return __fmodf( x, y );
00094 }
00095 #elif defined(__GNUC__)
00096 extern  float  __fmodf(float, float);
00097 float   fmodf(float, float) __attribute__ ((weak, alias ("__fmodf")));
00098 #endif
00099 
00100 static const  fu  Qnan = {QNANF};
00101 
00102 
00103 /* ====================================================================
00104  *
00105  * FunctionName   fmodf
00106  *
00107  * Description    computes fmodf function of args
00108  *
00109  * ====================================================================
00110  */
00111 
00112 float
00113 __fmodf( float x, float y )
00114 {
00115 #ifdef _32BIT_MACHINE
00116 
00117 int ix, iy, xptx, xpty, m;
00118 
00119 #else
00120 
00121 long long ix, iy, xptx, xpty, m;
00122 
00123 #endif
00124 
00125 #if defined(mips) && !defined(__GNUC__)
00126 fp_rnd  rm;
00127 #endif
00128 
00129 #ifdef __GNUC__
00130 int rm;
00131 #endif
00132 
00133 double  dx, dy;
00134 double  nd;
00135 double  dy1;
00136 #ifdef _CALL_MATHERR
00137 struct exception  exstruct;
00138 #endif
00139 
00140   /* filter out Nans */
00141 
00142   if ( (x != x) || (y != y) )
00143   {
00144     /* x or y is a NaN; return a quiet NaN */
00145 
00146 #ifdef _CALL_MATHERR
00147 
00148                 exstruct.type = DOMAIN;
00149                 exstruct.name = "fmodf";
00150                 exstruct.arg1 = x;
00151                 exstruct.arg2 = y;
00152                 exstruct.retval = Qnan.f;
00153 
00154                 if ( matherr( &exstruct ) == 0 )
00155                 {
00156                         fprintf(stderr, "domain error in fmodf\n");
00157                         SETERRNO(EDOM);
00158                 }
00159 
00160                 return ( exstruct.retval );
00161 #else
00162     NAN_SETERRNO(EDOM);
00163 
00164     return ( Qnan.f );
00165 #endif
00166   }
00167 
00168   dx = x;
00169   dy = y;
00170 
00171   /* extract exponent of dx for some quick screening */
00172 
00173 #ifdef _32BIT_MACHINE
00174 
00175   DBLHI2INT(dx, ix);  /* copy MSW of dx to ix */
00176   DBLHI2INT(dy, iy);  /* copy MSW of dy to iy */
00177 #else
00178   DBL2LL(dx, ix);   /* copy dx to ix */
00179   DBL2LL(dy, iy);   /* copy dy to iy */
00180 #endif
00181   xptx = (ix >> DMANTWIDTH);
00182   xptx &= 0x7ff;
00183 
00184   xpty = (iy >> DMANTWIDTH);
00185   xpty &= 0x7ff;
00186 
00187   if ( (xptx == 0x7ff) || ((xpty == 0) && (y == 0.0f)) )
00188   {
00189     /* x == +/- infinity or y == +/- 0.0; return a quiet NaN */
00190 
00191 #ifdef _CALL_MATHERR
00192 
00193     exstruct.type = DOMAIN;
00194     exstruct.name = "fmodf";
00195     exstruct.arg1 = x;
00196     exstruct.arg2 = y;
00197     exstruct.retval = Qnan.f;
00198 
00199     if ( matherr( &exstruct ) == 0 )
00200     {
00201       fprintf(stderr, "domain error in fmodf\n");
00202       SETERRNO(EDOM);
00203     }
00204 
00205     return ( exstruct.retval );
00206 #else
00207     SETERRNO(EDOM);
00208 
00209     return ( Qnan.f );
00210 #endif
00211   }
00212 
00213   dy = fabs(dy);
00214 
00215   if ( fabs(dx) < dy )
00216     return ( x );
00217 
00218   /* set rounding mode to round to zero */
00219   
00220 #if defined(mips) && !defined(__GNUC__)
00221   rm = swapRM( FP_RZ );
00222 #endif
00223 
00224 #ifdef __GNUC__
00225   rm = swapRM( FE_TOWARDZERO );
00226 #endif
00227 
00228 
00229   if ( xptx < xpty + 24 )
00230   {
00231     /* compute dx = dx - floor(dx/dy)*dy
00232     */
00233 
00234     nd = dx/dy;
00235 
00236     nd = (int)nd;
00237 
00238     dx = dx - nd*dy;
00239 
00240     rm = swapRM( rm );
00241 
00242     goto last;
00243   }
00244   else do
00245   {
00246     /* scale dy up and compute dx = dx - floor(dx/(k*dy))*k*dy
00247     */
00248 
00249     dy1 = dy;
00250 
00251 #ifdef _32BIT_MACHINE
00252 
00253     DBLHI2INT(dy1, m);
00254   
00255 #else
00256     DBL2LL(dy1, m);
00257 #endif
00258     m &= DEXPMASK;  /* clear exponent of dy1  */
00259     m |= ((xptx - 23) << DMANTWIDTH);
00260 
00261     /* set exponent of dy1 to xptx - 23 */
00262 
00263 #ifdef _32BIT_MACHINE
00264 
00265     INT2DBLHI(m, dy1);
00266 #else
00267     LL2DBL(m, dy1);
00268 #endif
00269     nd = dx/dy1;
00270 
00271     nd = (int)nd;
00272 
00273     dx = dx - nd*dy1;
00274 
00275 #ifdef _32BIT_MACHINE
00276 
00277     DBLHI2INT(dx, ix);  /* copy MSW of dx to ix */
00278 #else
00279     DBL2LL(dx, ix);   /* copy dx to ix */
00280 #endif
00281     xptx = (ix >> DMANTWIDTH);
00282     xptx &= 0x7ff;
00283 
00284   } 
00285   while ( xptx >= xpty + 24 );
00286 
00287   if ( fabs(dx) < dy )
00288   {
00289     rm = swapRM( rm );
00290 
00291     goto last;
00292   }
00293 
00294   /* compute dx = dx - floor(dx/dy)*dy
00295   */
00296 
00297   nd = dx/dy;
00298 
00299   nd = (int)nd;
00300 
00301   dx = dx - nd*dy;
00302 
00303   rm = swapRM( rm );
00304 
00305 last: if ( (x < 0.0f) && (dx == 0.0) )
00306     dx = -0.0;
00307 
00308   return ( (float)dx );
00309 }
00310