sys/src/libc/port/frexp.c

   1 #include <u.h>
   2 #include <libc.h>
   3
   4 /*
   5  * this is big/little endian non-portable
   6  * it gets the endian from the FPdbleword
   7  * union in u.h.
   8  */
   9 #define MASK    0x7ffL
  10 #define SHIFT   20
  11 #define BIAS    1022L
  12 #define SIG     52
  13
  14 double
  15 frexp(double d, int *ep)
  16 {
  17         FPdbleword x, a;
  18
  19         *ep = 0;
  20         /* order matters: only isNaN can operate on NaN */
  21         if(isNaN(d) || isInf(d, 0) || d == 0)
  22                 return d;
  23         x.x = d;
  24         a.x = fabs(d);
  25         if((a.hi >> SHIFT) == 0){       /* normalize subnormal numbers */
  26                 x.x = (double)(1ULL<<SIG) * d;
  27                 *ep = -SIG;
  28         }
  29         *ep += ((x.hi >> SHIFT) & MASK) - BIAS;
  30         x.hi &= ~(MASK << SHIFT);
  31         x.hi |= BIAS << SHIFT;
  32         return x.x;
  33 }
  34
  35 double
  36 ldexp(double d, int deltae)
  37 {
  38         int e, bits;
  39         FPdbleword x;
  40         ulong z;
  41
  42         if(d == 0)
  43                 return 0;
  44         x.x = d;
  45         e = (x.hi >> SHIFT) & MASK;
  46         if(deltae >= 0 || e+deltae >= 1){       /* no underflow */
  47                 e += deltae;
  48                 if(e >= MASK){          /* overflow */
  49                         if(d < 0)
  50                                 return Inf(-1);
  51                         return Inf(1);
  52                 }
  53         }else{  /* underflow gracefully */
  54                 deltae = -deltae;
  55                 /* need to shift d right deltae */
  56                 if(e > 1){              /* shift e-1 by exponent manipulation */
  57                         deltae -= e-1;
  58                         e = 1;
  59                 }
  60                 if(deltae > 0 && e==1){ /* shift 1 by switch from 1.fff to 0.1ff */
  61                         deltae--;
  62                         e = 0;
  63                         x.lo >>= 1;
  64                         x.lo |= (x.hi&1)<<31;
  65                         z = x.hi & ((1<<SHIFT)-1);
  66                         x.hi &= ~((1<<SHIFT)-1);
  67                         x.hi |= (1<<(SHIFT-1)) | (z>>1);
  68                 }
  69                 while(deltae > 0){              /* shift bits down */
  70                         bits = deltae;
  71                         if(bits > SHIFT)
  72                                 bits = SHIFT;
  73                         x.lo >>= bits;
  74                         x.lo |= (x.hi&((1<<bits)-1)) << (32-bits);
  75                         z = x.hi & ((1<<SHIFT)-1);
  76                         x.hi &= ~((1<<SHIFT)-1);
  77                         x.hi |= z>>bits;
  78                         deltae -= bits;
  79                 }
  80         }
  81         x.hi &= ~(MASK << SHIFT);
  82         x.hi |= (long)e << SHIFT;
  83         return x.x;
  84 }
  85
  86 double
  87 modf(double d, double *ip)
  88 {
  89         FPdbleword x;
  90         int e;
  91
  92         x.x = d;
  93         e = (x.hi >> SHIFT) & MASK;
  94         if(e == MASK){
  95                 *ip = d;
  96                 if(x.lo != 0 || (x.hi & 0xfffffL != 0)) /* NaN */
  97                         return d;
  98                 /* ±Inf */
  99                 x.hi &= 0x80000000L;
 100                 return x.x;
 101         }
 102         if(d < 1) {
 103                 if(d < 0) {
 104                         x.x = modf(-d, ip);
 105                         *ip = -*ip;
 106                         return -x.x;
 107                 }
 108                 *ip = 0;
 109                 return d;
 110         }
 111         e -= BIAS;
 112         if(e <= SHIFT+1) {
 113                 x.hi &= ~(0x1fffffL >> e);
 114                 x.lo = 0;
 115         } else
 116         if(e <= SHIFT+33)
 117                 x.lo &= ~(0x7fffffffL >> (e-SHIFT-2));
 118         *ip = x.x;
 119         return d - x.x;
 120 }