lib/ack/libfp/compact.fc

   1 /*
   2   (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
   3   See the copyright notice in the ACK home directory, in the file "Copyright".
   4 */
   5
   6 /* $Header$ */
   7
   8 /*
   9         COMPACT EXTEND FORMAT INTO FLOAT OF PROPER SIZE
  10 */
  11
  12 # include "FP_bias.h"
  13 # include "FP_shift.h"
  14 # include "FP_trap.h"
  15 # include "FP_types.h"
  16 # include "get_put.h"
  17
  18 void
  19 compact(f,to,size)
  20 EXTEND  *f;
  21 unsigned long   *to;
  22 int     size;
  23 {
  24         int     error = 0;
  25
  26         if (size == sizeof(DOUBLE)) {
  27         /*
  28          * COMPACT EXTENDED INTO DOUBLE
  29          */
  30                 DOUBLE *DBL = (DOUBLE *) (void *) to;
  31
  32                 if ((f->m1|(f->m2 & DBL_ZERO)) == 0L)   {
  33                         zrf8(DBL);
  34                         return;
  35                 }
  36                 f->exp += DBL_BIAS;     /* restore proper bias  */
  37                 if (f->exp > DBL_MAX)   {
  38 dbl_over:                       trap(EFOVFL);
  39                         f->exp = DBL_MAX+1;
  40                         f->m1 = 0;
  41                         f->m2 = 0;
  42                         if (error++)
  43                                 return;
  44                 }
  45                 else if (f->exp < DBL_MIN)      {
  46                         b64_rsft(&(f->mantissa));
  47                         if (f->exp < 0) {
  48                                 b64_sft(&(f->mantissa), -f->exp);
  49                                 f->exp = 0;
  50                         }
  51                         /* underflow ??? */
  52                 }
  53
  54                 /* local CAST conversion                */
  55
  56                 /* because of special format shift only 10 bits */
  57                 /* bit shift mantissa 10 bits           */
  58
  59                 /* first align within words, then do store operation */
  60
  61                 DBL->d[0] = f->m1 >> DBL_RUNPACK;   /* plus 22 == 32 */
  62                 DBL->d[1] = f->m2 >> DBL_RUNPACK;   /* plus 22 == 32 */
  63                 DBL->d[1] |= (f->m1 << DBL_LUNPACK); /* plus 10 == 32 */
  64
  65                 /* if not exact then round to nearest   */
  66                 /* on a tie, round to even */
  67
  68 #ifdef EXCEPTION_INEXACT
  69                 if ((f->m2 & DBL_EXACT) != 0) {
  70                     INEXACT();
  71 #endif
  72                     if (((f->m2 & DBL_EXACT) > DBL_ROUNDUP)
  73                         || ((f->m2 & DBL_EXACT) == DBL_ROUNDUP
  74                             && (f->m2 & (DBL_ROUNDUP << 1)))) {
  75                         DBL->d[1]++;    /* rounding up  */
  76                         if (DBL->d[1] == 0L) { /* carry out     */
  77                             DBL->d[0]++;
  78
  79                             if (f->exp == 0 && (DBL->d[0] & ~DBL_MASK)) {
  80                                         f->exp++;
  81                                 }
  82                             if (DBL->d[0] & DBL_CARRYOUT) { /* carry out */
  83                                 if (DBL->d[0] & 01)
  84                                     DBL->d[1] = CARRYBIT;
  85                                 DBL->d[0] >>= 1;
  86                                 f->exp++;
  87                             }
  88                         }
  89                         /*      check for overflow                      */
  90                         if (f->exp > DBL_MAX)
  91                                 goto dbl_over;
  92                     }
  93 #ifdef EXCEPTION_INEXACT
  94                 }
  95 #endif
  96
  97                 /*
  98                  * STORE EXPONENT AND SIGN:
  99                  *
 100                  * 1) clear leading bits (B4-B15)
 101                  * 2) shift and store exponent
 102                  */
 103
 104                 DBL->d[0] &= DBL_MASK;
 105                 DBL->d[0] |=
 106                         ((long) (f->exp << DBL_EXPSHIFT) << EXP_STORE);
 107                 if (f->sign)
 108                         DBL->d[0] |= CARRYBIT;
 109
 110                 /*
 111                  * STORE MANTISSA
 112                  */
 113
 114 #if FL_MSL_AT_LOW_ADDRESS
 115                 put4(DBL->d[0], (char *) &DBL->d[0]);
 116                 put4(DBL->d[1], (char *) &DBL->d[1]);
 117 #else
 118                 { unsigned long l;
 119                   put4(DBL->d[1], (char *) &l);
 120                   put4(DBL->d[0], (char *) &DBL->d[1]);
 121                   DBL->d[0] = l;
 122                 }
 123 #endif
 124         }
 125         else {
 126                 /*
 127                  * COMPACT EXTENDED INTO FLOAT
 128                  */
 129                 SINGLE  *SGL;
 130
 131                 /* local CAST conversion                */
 132                 SGL = (SINGLE *) (void *) to;
 133                 if ((f->m1 & SGL_ZERO) == 0L)   {
 134                         *SGL = 0L;
 135                         return;
 136                 }
 137                 f->exp += SGL_BIAS;     /* restore bias */
 138                 if (f->exp > SGL_MAX)   {
 139 sgl_over:                       trap(EFOVFL);
 140                         f->exp = SGL_MAX+1;
 141                         f->m1 = 0L;
 142                         f->m2 = 0L;
 143                         if (error++)
 144                                 return;
 145                 }
 146                 else if (f->exp < SGL_MIN)      {
 147                         b64_rsft(&(f->mantissa));
 148                         if (f->exp < 0) {
 149                                 b64_sft(&(f->mantissa), -f->exp);
 150                                 f->exp = 0;
 151                         }
 152                         /* underflow ??? */
 153                 }
 154
 155                 /* shift mantissa and store     */
 156                 *SGL = (f->m1 >> SGL_RUNPACK);
 157
 158                 /* check for rounding to nearest        */
 159                 /* on a tie, round to even              */
 160 #ifdef EXCEPTION_INEXACT
 161                 if (f->m2 != 0 ||
 162                     (f->m1 & SGL_EXACT) != 0L) {
 163                         INEXACT();
 164 #endif
 165                         if (((f->m1 & SGL_EXACT) > SGL_ROUNDUP)
 166                             || ((f->m1 & SGL_EXACT) == SGL_ROUNDUP
 167                                 && (f->m1 & (SGL_ROUNDUP << 1)))) {
 168                                 (*SGL)++;
 169                                 if (f->exp == 0 && (*SGL & ~SGL_MASK)) {
 170                                         f->exp++;
 171                                 }
 172                         /* check normal */
 173                                 if (*SGL & SGL_CARRYOUT)        {
 174                                         *SGL >>= 1;
 175                                         f->exp++;
 176                                 }
 177                                 if (f->exp > SGL_MAX)
 178                                         goto sgl_over;
 179                         }
 180 #ifdef EXCEPTION_INEXACT
 181                 }
 182 #endif
 183
 184                 /*
 185                  * STORE EXPONENT AND SIGN:
 186                  *
 187                  * 1) clear leading bit of fraction
 188                  * 2) shift and store exponent
 189                  */
 190
 191                 *SGL &= SGL_MASK; /* B23-B31 are 0 */
 192                 *SGL |= ((long) (f->exp << SGL_EXPSHIFT) << EXP_STORE);
 193                 if (f->sign)
 194                         *SGL |= CARRYBIT;
 195
 196                 /*
 197                  * STORE MANTISSA
 198                  */
 199
 200                 put4(*SGL, (char *) &SGL);
 201         }
 202 }