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  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  23  */
  24 /*
  25  * Copyright 2005 Sun Microsystems, Inc.  All rights reserved.
  26  * Use is subject to license terms.
  27  */
  28 
  29 #pragma weak fmodf = __fmodf
  30 
  31 #include "libm.h"
  32 
  33 /* INDENT OFF */
  34 static const int
  35         is = (int)0x80000000,
  36         im = 0x007fffff,
  37         ii = 0x7f800000,
  38         iu = 0x00800000;
  39 /* INDENT ON */
  40 
  41 static const float zero = 0.0;
  42 
  43 float
  44 fmodf(float x, float y) {
  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                 /* INDENT OFF */
  63                 /*
  64                  * scale x,y to "normal" with
  65                  *      ny = exponent of y
  66                  *      nd = exponent of x minus exponent of y
  67                  */
  68                 /* INDENT ON */
  69                 ny = iy >> 23;
  70                 k = ix >> 23;
  71 
  72                 /* special case for subnormal y or x */
  73                 if (ny == 0) {
  74                         ny = 1;
  75                         while (iy < iu) {
  76                                 ny -= 1;
  77                                 iy += iy;
  78                         }
  79                         nd = k - ny;
  80                         if (k == 0) {
  81                                 nd += 1;
  82                                 while (ix < iu) {
  83                                         nd -= 1;
  84                                         ix += ix;
  85                                 }
  86                         } else {
  87                                 ix = iu | (ix & im);
  88                         }
  89                 } else {
  90                         nd = k - ny;
  91                         ix = iu | (ix & im);
  92                         iy = iu | (iy & im);
  93                 }
  94 
  95                 /* fix point fmod for normalized ix and iy */
  96                 /* INDENT OFF */
  97                 /*
  98                  * while (nd--) {
  99                  *      iz = ix - iy;
 100                  * if (iz < 0)
 101                  *      ix = ix + ix;
 102                  * else if (iz == 0) {
 103                  *      *(int *) &w = is & hx;
 104                  *      return w;
 105                  * }
 106                  * else
 107                  *      ix = iz + iz;
 108                  * }
 109                  */
 110                 /* INDENT ON */
 111                 /* unroll the above loop 4 times to gain performance */
 112                 k = nd >> 2;
 113                 nd -= k << 2;
 114                 while (k--) {
 115                         iz = ix - iy;
 116                         if (iz >= 0)
 117                                 ix = iz + iz;
 118                         else
 119                                 ix += ix;
 120                         iz = ix - iy;
 121                         if (iz >= 0)
 122                                 ix = iz + iz;
 123                         else
 124                                 ix += ix;
 125                         iz = ix - iy;
 126                         if (iz >= 0)
 127                                 ix = iz + iz;
 128                         else
 129                                 ix += ix;
 130                         iz = ix - iy;
 131                         if (iz >= 0)
 132                                 ix = iz + iz;
 133                         else
 134                                 ix += ix;
 135                         if (iz == 0) {
 136                                 *(int *)&w = is & hx;
 137                                 return (w);
 138                         }
 139                 }
 140                 while (nd--) {
 141                         iz = ix - iy;
 142                         if (iz >= 0)
 143                                 ix = iz + iz;
 144                         else
 145                                 ix += ix;
 146                 }
 147                 /* end of unrolling */
 148 
 149                 iz = ix - iy;
 150                 if (iz >= 0)
 151                         ix = iz;
 152 
 153                 /* convert back to floating value and restore the sign */
 154                 if (ix == 0) {
 155                         *(int *)&w = is & hx;
 156                         return (w);
 157                 }
 158                 while (ix < iu) {
 159                         ix += ix;
 160                         ny -= 1;
 161                 }
 162                 while (ix > (iu + iu)) {
 163                         ny += 1;
 164                         ix >>= 1;
 165                 }
 166                 if (ny > 0) {
 167                         *(int *)&w = (is & hx) | (ix & im) | (ny << 23);
 168                 } else {
 169                         /* subnormal output */
 170                         k = -ny + 1;
 171                         ix >>= k;
 172                         *(int *)&w = (is & hx) | ix;
 173                 }
 174         }
 175         return (w);
 176 }