1 /*
   2  * CDDL HEADER START
   3  *
   4  * The contents of this file are subject to the terms of the
   5  * Common Development and Distribution License (the "License").
   6  * You may not use this file except in compliance with the License.
   7  *
   8  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
   9  * or http://www.opensolaris.org/os/licensing.
  10  * See the License for the specific language governing permissions
  11  * and limitations under the License.
  12  *
  13  * When distributing Covered Code, include this CDDL HEADER in each
  14  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
  15  * If applicable, add the following below this CDDL HEADER, with the
  16  * fields enclosed by brackets "[]" replaced with your own identifying
  17  * information: Portions Copyright [yyyy] [name of copyright owner]
  18  *
  19  * CDDL HEADER END
  20  */
  21 
  22 /*
  23  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  24  */
  25 
  26 /*
  27  * Copyright 2005 Sun Microsystems, Inc.  All rights reserved.
  28  * Use is subject to license terms.
  29  */
  30 
  31 #pragma weak __fmodf = fmodf
  32 
  33 #include "libm.h"
  34 
  35 static const int is = (int)0x80000000,
  36         im = 0x007fffff,
  37         ii = 0x7f800000,
  38         iu = 0x00800000;
  39 
  40 static const float zero = 0.0;
  41 
  42 float
  43 fmodf(float x, float y)
  44 {
  45         float w;
  46         int hx, ix, iy, iz, k, ny, nd;
  47 
  48         hx = *(int *)&x;
  49         ix = hx & 0x7fffffff;
  50         iy = *(int *)&y & 0x7fffffff;
  51 
  52         /* purge off exception values */
  53         if (ix >= ii || iy > ii || iy == 0) {
  54                 w = x * y;
  55                 w = w / w;
  56         } else if (ix <= iy) {
  57                 if (ix < iy)
  58                         w = x;          /* return x if |x|<|y| */
  59                 else
  60                         w = zero * x;   /* return sign(x)*0.0 */
  61         } else {
  62 
  63                 /*
  64                  * scale x,y to "normal" with
  65                  *      ny = exponent of y
  66                  *      nd = exponent of x minus exponent of y
  67                  */
  68                 ny = iy >> 23;
  69                 k = ix >> 23;
  70 
  71                 /* special case for subnormal y or x */
  72                 if (ny == 0) {
  73                         ny = 1;
  74 
  75                         while (iy < iu) {
  76                                 ny -= 1;
  77                                 iy += iy;
  78                         }
  79 
  80                         nd = k - ny;
  81 
  82                         if (k == 0) {
  83                                 nd += 1;
  84 
  85                                 while (ix < iu) {
  86                                         nd -= 1;
  87                                         ix += ix;
  88                                 }
  89                         } else {
  90                                 ix = iu | (ix & im);
  91                         }
  92                 } else {
  93                         nd = k - ny;
  94                         ix = iu | (ix & im);
  95                         iy = iu | (iy & im);
  96                 }
  97 
  98                 /*
  99                  * fix point fmod for normalized ix and iy
 100                  */
 101 
 102                 /* BEGIN CSTYLED */
 103                 /*
 104                  * while (nd--) {
 105                  *      iz = ix - iy;
 106                  * if (iz < 0)
 107                  *      ix = ix + ix;
 108                  * else if (iz == 0) {
 109                  *      *(int *) &w = is & hx;
 110                  *      return w;
 111                  * }
 112                  * else
 113                  *      ix = iz + iz;
 114                  * }
 115                  */
 116                 /* END CSTYLED */
 117 
 118                 /*
 119                  * unroll the above loop 4 times to gain performance
 120                  */
 121                 k = nd >> 2;
 122                 nd -= k << 2;
 123 
 124                 while (k--) {
 125                         iz = ix - iy;
 126 
 127                         if (iz >= 0)
 128                                 ix = iz + iz;
 129                         else
 130                                 ix += ix;
 131 
 132                         iz = ix - iy;
 133 
 134                         if (iz >= 0)
 135                                 ix = iz + iz;
 136                         else
 137                                 ix += ix;
 138 
 139                         iz = ix - iy;
 140 
 141                         if (iz >= 0)
 142                                 ix = iz + iz;
 143                         else
 144                                 ix += ix;
 145 
 146                         iz = ix - iy;
 147 
 148                         if (iz >= 0)
 149                                 ix = iz + iz;
 150                         else
 151                                 ix += ix;
 152 
 153                         if (iz == 0) {
 154                                 *(int *)&w = is & hx;
 155                                 return (w);
 156                         }
 157                 }
 158 
 159                 while (nd--) {
 160                         iz = ix - iy;
 161 
 162                         if (iz >= 0)
 163                                 ix = iz + iz;
 164                         else
 165                                 ix += ix;
 166                 }
 167 
 168                 /* end of unrolling */
 169 
 170                 iz = ix - iy;
 171 
 172                 if (iz >= 0)
 173                         ix = iz;
 174 
 175                 /* convert back to floating value and restore the sign */
 176                 if (ix == 0) {
 177                         *(int *)&w = is & hx;
 178                         return (w);
 179                 }
 180 
 181                 while (ix < iu) {
 182                         ix += ix;
 183                         ny -= 1;
 184                 }
 185 
 186                 while (ix > (iu + iu)) {
 187                         ny += 1;
 188                         ix >>= 1;
 189                 }
 190 
 191                 if (ny > 0) {
 192                         *(int *)&w = (is & hx) | (ix & im) | (ny << 23);
 193                 } else {
 194                         /* subnormal output */
 195                         k = -ny + 1;
 196                         ix >>= k;
 197                         *(int *)&w = (is & hx) | ix;
 198                 }
 199         }
 200 
 201         return (w);
 202 }