fp.c

   1 #ifndef pdp11
   2 #include <errno.h>
   3 #include <stdlib.h>
   4 #include "c1.h"
   5
   6 #ifndef TRUE
   7 #define TRUE 1
   8 #endif
   9 #ifndef FALSE
  10 #define FALSE 1
  11 #endif
  12
  13 /* Floating point accumulators */
  14
  15 #if 0
  16 /* PSW */
  17
  18 #define PSW_V_C         0                               /* condition codes */
  19 #define PSW_V_V         1
  20 #define PSW_V_Z         2
  21 #define PSW_V_N         3
  22 #endif
  23
  24 /* FPS */
  25
  26 #define FPS_V_C         0                               /* condition codes */
  27 #define FPS_V_V         1
  28 #define FPS_V_Z         2
  29 #define FPS_V_N         3
  30 #define FPS_V_T         5                               /* truncate */
  31 #define FPS_V_L         6                               /* long */
  32 #define FPS_V_D         7                               /* double */
  33 #define FPS_V_IC        8                               /* ic err int */
  34 #define FPS_V_IV        9                               /* overflo err int */
  35 #define FPS_V_IU        10                              /* underflo err int */
  36 #define FPS_V_IUV       11                              /* undef var err int */
  37 #define FPS_V_ID        14                              /* int disable */
  38 #define FPS_V_ER        15                              /* error */
  39
  40 /* Floating point status register */
  41
  42 #define FPS_ER          (1u << FPS_V_ER)                /* error */
  43 #define FPS_ID          (1u << FPS_V_ID)                /* interrupt disable */
  44 #define FPS_IUV         (1u << FPS_V_IUV)               /* int on undef var */
  45 #define FPS_IU          (1u << FPS_V_IU)                /* int on underflow */
  46 #define FPS_IV          (1u << FPS_V_IV)                /* int on overflow */
  47 #define FPS_IC          (1u << FPS_V_IC)                /* int on conv error */
  48 #define FPS_D           (1u << FPS_V_D)                 /* single/double */
  49 #define FPS_L           (1u << FPS_V_L)                 /* word/long */
  50 #define FPS_T           (1u << FPS_V_T)                 /* round/truncate */
  51 #define FPS_N           (1u << FPS_V_N)
  52 #define FPS_Z           (1u << FPS_V_Z)
  53 #define FPS_V           (1u << FPS_V_V)
  54 #define FPS_C           (1u << FPS_V_C)
  55 #define FPS_CC          (FPS_N + FPS_Z + FPS_V + FPS_C)
  56 #define FPS_RW          (FPS_ER + FPS_ID + FPS_IUV + FPS_IU + FPS_IV + \
  57                         FPS_IC + FPS_D + FPS_L + FPS_T + FPS_CC)
  58
  59 /* Floating point exception codes */
  60
  61 #define FEC_OP          2                               /* illegal op/mode */
  62 #define FEC_DZRO        4                               /* divide by zero */
  63 #define FEC_ICVT        6                               /* conversion error */
  64 #define FEC_OVFLO       8                               /* overflow */
  65 #define FEC_UNFLO       10                              /* underflow */
  66 #define FEC_UNDFV       12                              /* undef variable */
  67
  68 /* Floating point format, all assignments 32b relative */
  69
  70 #define FP_V_SIGN       (63 - 32)                       /* high lw: sign */
  71 #define FP_V_EXP        (55 - 32)                       /* exponent */
  72 #define FP_V_HB         FP_V_EXP                        /* hidden bit */
  73 #define FP_V_F0         (48 - 32)                       /* fraction 0 */
  74 #define FP_V_F1         (32 - 32)                       /* fraction 1 */
  75 #define FP_V_FROUND     (31 - 32)                       /* f round point */
  76 #define FP_V_F2         16                              /* low lw: fraction 2 */
  77 #define FP_V_F3         0                               /* fraction 3 */
  78 #define FP_V_DROUND     (-1)                            /* d round point */
  79 #define FP_M_EXP        0377
  80 #define FP_SIGN         (1u << FP_V_SIGN)
  81 #define FP_EXP          (FP_M_EXP << FP_V_EXP)
  82 #define FP_HB           (1u << FP_V_HB)
  83 #define FP_FRACH        ((1u << FP_V_HB) - 1)
  84 #define FP_FRACL        0xFFFFFFFF
  85 #define FP_BIAS         0200                            /* exponent bias */
  86 #define FP_GUARD        3                               /* guard bits */
  87
  88 #if 0
  89 /* Data lengths */
  90
  91 #define _WORD   2
  92 #define _LONG   4
  93 #define _QUAD   8
  94 #endif
  95
  96 /* Double precision operations on 64b quantities */
  97
  98 #define F_LOAD(qd,ac,ds) ds.h = ac.h; ds.l = (qd)? ac.l: 0
  99 #define F_LOAD_P(qd,ac,ds) ds->h = ac.h; ds->l = (qd)? ac.l: 0
 100 #define F_LOAD_FRAC(qd,ac,ds) ds.h = (ac.h & FP_FRACH) | FP_HB; \
 101         ds.l = (qd)? ac.l: 0
 102 #define F_STORE(qd,sr,ac) ac.h = sr.h; if ((qd)) ac.l = sr.l
 103 #define F_STORE_P(qd,sr,ac) ac.h = sr->h; if ((qd)) ac.l = sr->l
 104 #define F_GET_FRAC_P(sr,ds) ds.l = sr->l; \
 105         ds.h = (sr->h & FP_FRACH) | FP_HB
 106 #define F_ADD(s2,s1,ds) ds.l = (s1.l + s2.l) & 0xFFFFFFFF; \
 107         ds.h = (s1.h + s2.h + (ds.l < s2.l)) & 0xFFFFFFFF
 108 #define F_SUB(s2,s1,ds) ds.h = (s1.h - s2.h - (s1.l < s2.l)) & 0xFFFFFFFF; \
 109         ds.l = (s1.l - s2.l) & 0xFFFFFFFF
 110 #define F_LT(x,y) ((x.h < y.h) || ((x.h == y.h) && (x.l < y.l)))
 111 #define F_LT_AP(x,y) (((x->h & ~FP_SIGN) < (y->h & ~FP_SIGN)) || \
 112         (((x->h & ~FP_SIGN) == (y->h & ~FP_SIGN)) && (x->l < y->l)))
 113 #define F_LSH_V(sr,n,ds) \
 114         ds.h = (((n) >= 32)? (sr.l << ((n) - 32)): \
 115                 (sr.h << (n)) | ((sr.l >> (32 - (n))) & and_mask[n])) \
 116                 & 0xFFFFFFFF; \
 117         ds.l = ((n) >= 32)? 0: (sr.l << (n)) & 0xFFFFFFFF
 118 #define F_RSH_V(sr,n,ds) \
 119         ds.l = (((n) >= 32)? (sr.h >> ((n) - 32)) & and_mask[64 - (n)]: \
 120                 ((sr.l >> (n)) & and_mask[32 - (n)]) | \
 121                 (sr.h << (32 - (n)))) & 0xFFFFFFFF; \
 122         ds.h = ((n) >= 32)? 0: \
 123                 ((sr.h >> (n)) & and_mask[32 - (n)]) & 0xFFFFFFFF
 124
 125 /* For the constant shift macro, arguments must in the range [2,31] */
 126
 127 #define F_LSH_1(ds) ds.h = ((ds.h << 1) | ((ds.l >> 31) & 1)) & 0xFFFFFFFF; \
 128         ds.l = (ds.l << 1) & 0xFFFFFFFF
 129 #define F_RSH_1(ds) ds.l = ((ds.l >> 1) & 0x7FFFFFFF) | ((ds.h & 1) << 31); \
 130         ds.h = ((ds.h >> 1) & 0x7FFFFFFF)
 131 #define F_LSH_K(sr,n,ds) \
 132         ds.h =  ((sr.h << (n)) | ((sr.l >> (32 - (n))) & and_mask[n])) \
 133                 & 0xFFFFFFFF; \
 134         ds.l = (sr.l << (n)) & 0xFFFFFFFF
 135 #define F_RSH_K(sr,n,ds) \
 136         ds.l =  (((sr.l >> (n)) & and_mask[32 - (n)]) | \
 137                 (sr.h << (32 - (n)))) & 0xFFFFFFFF; \
 138         ds.h =  ((sr.h >> (n)) & and_mask[32 - (n)]) & 0xFFFFFFFF
 139 #define F_LSH_GUARD(ds) F_LSH_K(ds,FP_GUARD,ds)
 140 #define F_RSH_GUARD(ds) F_RSH_K(ds,FP_GUARD,ds)
 141
 142 #define GET_BIT(ir,n)   (((ir) >> (n)) & 1)
 143 #define GET_SIGN(ir)    GET_BIT((ir), FP_V_SIGN)
 144 #define GET_EXP(ir)     (((ir) >> FP_V_EXP) & FP_M_EXP)
 145 #define GET_SIGN_L(ir)  GET_BIT((ir), 31)
 146 #define GET_SIGN_W(ir)  GET_BIT((ir), 15)
 147
 148 int32_t FEC = 0;
 149 int32_t FPS = FPS_D; /* default to double precision */
 150
 151 _DOUBLE zero_fac = { 0, 0 };
 152 _DOUBLE one_fac = { 1, 0 };
 153 _DOUBLE fround_fac = { (1u << (FP_V_FROUND + 32)), 0 };
 154 _DOUBLE fround_guard_fac = { 0, (1u << (FP_V_FROUND + FP_GUARD)) };
 155 _DOUBLE dround_guard_fac = { (1u << (FP_V_DROUND + FP_GUARD)), 0 };
 156 _DOUBLE fmask_fac = { 0xFFFFFFFF, (1u << (FP_V_HB + FP_GUARD + 1)) - 1 };
 157 static const uint32_t and_mask[33] = { 0,
 158         0x1, 0x3, 0x7, 0xF,
 159         0x1F, 0x3F, 0x7F, 0xFF,
 160         0x1FF, 0x3FF, 0x7FF, 0xFFF,
 161         0x1FFF, 0x3FFF, 0x7FFF, 0xFFFF,
 162         0x1FFFF, 0x3FFFF, 0x7FFFF, 0xFFFFF,
 163         0x1FFFFF, 0x3FFFFF, 0x7FFFFF, 0xFFFFFF,
 164         0x1FFFFFF, 0x3FFFFFF, 0x7FFFFFF, 0xFFFFFFF,
 165         0x1FFFFFFF, 0x3FFFFFFF, 0x7FFFFFFF, 0xFFFFFFFF };
 166
 167 #ifndef __P
 168 #ifdef __STDC__
 169 #define __P(params) params
 170 #else
 171 #define __P(params) ()
 172 #endif
 173 #endif
 174
 175 static void tstfp11 __P((_DOUBLE *fsrc));
 176 static void absfp11 __P((_DOUBLE *fsrc));
 177 static void negfp11 __P((_DOUBLE *fsrc));
 178 static void cmpfp11 __P((_DOUBLE *fac, _DOUBLE *fsrc));
 179 static int32_t moviefp11 __P((int32_t val, _DOUBLE *fac));
 180 static int32_t moveifp11 __P((_DOUBLE *val));
 181 static int32_t moviffp11 __P((int32_t val, _DOUBLE *fac));
 182 static int32_t movfifp11 __P((_DOUBLE *val, int32_t *dst));
 183 static int32_t addfp11 __P((_DOUBLE *facp, _DOUBLE *fsrcp));
 184 static int32_t mulfp11 __P((_DOUBLE *facp, _DOUBLE *fsrcp));
 185 static void frac_mulfp11 __P((_DOUBLE *f1p, _DOUBLE *f2p));
 186 static int32_t divfp11 __P((_DOUBLE *facp, _DOUBLE *fsrcp));
 187 static int32_t roundfp11 __P((_DOUBLE *fptr));
 188 static int32_t round_and_pack __P((_DOUBLE *facp, int32_t exp, _DOUBLE *fracp, int r));
 189 static int32_t fpnotrap __P((int32_t code));
 190
 191 int fp_tst(val) _DOUBLE val; {
 192         tstfp11(&val);
 193         return (FPS & FPS_Z) != 0;
 194 }
 195
 196 _DOUBLE fp_abs(val) _DOUBLE val; {
 197         absfp11(&val);
 198         return val;
 199 }
 200
 201 _DOUBLE fp_neg(val) _DOUBLE val; {
 202         negfp11(&val);
 203         return val;
 204 }
 205
 206 int fp_le(val0, val1) _DOUBLE val0; _DOUBLE val1; {
 207         cmpfp11(&val0, &val1);
 208         return (FPS & FPS_N) == 0;
 209 }
 210
 211 int fp_ge(val0, val1) _DOUBLE val0; _DOUBLE val1; {
 212         cmpfp11(&val1, &val0);
 213         return (FPS & FPS_N) == 0;
 214 }
 215
 216 int fp_gt(val0, val1) _DOUBLE val0; _DOUBLE val1; {
 217         cmpfp11(&val0, &val1);
 218         return (FPS & FPS_N) != 0;
 219 }
 220
 221 int fp_lt(val0, val1) _DOUBLE val0; _DOUBLE val1; {
 222         cmpfp11(&val1, &val0);
 223         return (FPS & FPS_N) != 0;
 224 }
 225
 226 _INT fp_double_to_int(val) _DOUBLE val; {
 227         int32_t res;
 228
 229         if (movfifp11(&val, &res))
 230                 abort();
 231         return (_INT)res;
 232 }
 233
 234 _LONG fp_double_to_long(val) _DOUBLE val; {
 235         int32_t res;
 236
 237         FPS |= FPS_L;
 238         if (movfifp11(&val, &res))
 239                 abort();
 240         FPS &= ~FPS_L;
 241         return (_LONG)res;
 242 }
 243
 244 _FLOAT fp_double_to_float(val) _DOUBLE val; {
 245         _FLOAT res;
 246
 247         if (roundfp11(&val))
 248                 abort();
 249         res.h = val.h;
 250         return res;
 251 }
 252
 253 _DOUBLE fp_int_to_double(val) _INT val; {
 254         _DOUBLE res;
 255
 256         if (moviffp11((int32_t)val << 16, &res))
 257                 abort();
 258         return res;
 259 }
 260
 261 _DOUBLE fp_long_to_double(val) _LONG val; {
 262         _DOUBLE res;
 263
 264         FPS |= FPS_L;
 265         if (moviffp11((int32_t)val, &res))
 266                 abort();
 267         FPS &= ~FPS_L;
 268         return res;
 269 }
 270
 271 _DOUBLE fp_float_to_double(val) _FLOAT val; {
 272         _DOUBLE res;
 273
 274         res.h = val.h;
 275         res.l = 0;
 276         return res;
 277 }
 278
 279 _DOUBLE fp_add(val0, val1) _DOUBLE val0; _DOUBLE val1; {
 280         if (addfp11(&val0, &val1))
 281                 abort();
 282         return val0;
 283 }
 284
 285 _DOUBLE fp_sub(val0, val1) _DOUBLE val0; _DOUBLE val1; {
 286         negfp11(&val1);
 287         if (addfp11(&val0, &val1))
 288                 abort();
 289         return val0;
 290 }
 291
 292 _DOUBLE fp_mul(val0, val1) _DOUBLE val0; _DOUBLE val1; {
 293         if (mulfp11(&val0, &val1))
 294                 abort();
 295         return val0;
 296 }
 297
 298 _DOUBLE fp_div(val0, val1) _DOUBLE val0; _DOUBLE val1; {
 299         if (divfp11(&val0, &val1))
 300                 abort();
 301         return val0;
 302 }
 303
 304 _DOUBLE fp_ldexp(val, exp) _DOUBLE val; _INT exp; {
 305         tstfp11(&val);
 306         if ((FPS & FPS_Z) == 0 && moviefp11(moveifp11(&val) + (int32_t)exp, &val)) {
 307                 val.l = 0xFFFFFFFF;
 308                 val.h = ((FPS << (31 - FPS_V_N)) | 0x7FFFFFFF) & 0xFFFFFFFF;
 309                 errno = ERANGE;
 310         }
 311         return val;
 312 }
 313
 314 /* cut-down routines by Nick, previously inline to instruction decode */
 315 static void tstfp11(fsrc) _DOUBLE *fsrc; {
 316         FPS = FPS & ~FPS_CC;
 317         if (GET_SIGN (fsrc->h)) FPS = FPS | FPS_N;
 318         if (GET_EXP (fsrc->h) == 0) FPS = FPS | FPS_Z;
 319 }
 320
 321 static void absfp11(fsrc) _DOUBLE *fsrc; {
 322         if (GET_EXP (fsrc->h) == 0) *fsrc = zero_fac;
 323         else fsrc->h = fsrc->h & ~FP_SIGN;
 324 }
 325
 326 static void negfp11(fsrc) _DOUBLE *fsrc; {
 327         if (GET_EXP (fsrc->h) == 0) *fsrc = zero_fac;
 328         else fsrc->h = fsrc->h ^ FP_SIGN;
 329 }
 330
 331 /* from PDP-11/70 processor manual:
 332  *   FC <- 0.
 333  *   FV <- 0.
 334  *   FZ <- 1 If (FSRC) - (AC) = 0, else FZ <- 0.
 335  *   FN <- 1 If (FSRC) - (AC) < 0, else FN <- 0.
 336  * (backwards compared to what I expected)
 337  */
 338 static void cmpfp11(fac, fsrc) _DOUBLE *fac; _DOUBLE *fsrc; {
 339         if (GET_EXP (fsrc->h) == 0) *fsrc = zero_fac;
 340         if (GET_EXP (fac->h) == 0) *fac = zero_fac;
 341         if ((fsrc->h == fac->h) && (fsrc->l == fac->l)) {       /* equal? */
 342             FPS = (FPS & ~FPS_CC) | FPS_Z;
 343             return;  }
 344         FPS = (FPS & ~FPS_CC) | ((fsrc->h >> (FP_V_SIGN - FPS_V_N)) & FPS_N);
 345         if ((GET_SIGN (fsrc->h ^ fac->h) == 0) && (fac->h != 0) &&
 346             F_LT ((*fsrc), (*fac))) FPS = FPS ^ FPS_N;
 347 }
 348
 349 static int32_t moviefp11(val, fac) int32_t val; _DOUBLE *fac; {
 350         fac->h = (fac->h & ~FP_EXP) | (((val + FP_BIAS) & FP_M_EXP) << FP_V_EXP);
 351         if ((val > 0177) && (val <= 0177600)) {
 352             if (val < 0100000) {
 353                 if (fpnotrap (FEC_OVFLO)) *fac = zero_fac;
 354                 return FPS_V;  }
 355             if (fpnotrap (FEC_UNFLO)) *fac = zero_fac;  }
 356         return 0;
 357 }
 358
 359 static int32_t moveifp11(val) _DOUBLE *val; {
 360         return (GET_EXP (val->h) - FP_BIAS) & 0177777;
 361 }
 362
 363 static int32_t moviffp11(val, fac) int32_t val; _DOUBLE *fac; {
 364         int32_t i;
 365         int32_t exp, sign;
 366
 367         fac->l = val;
 368         fac->h = 0;
 369         if (fac->l) {
 370             if (sign = GET_SIGN_L (fac->l)) fac->l = (fac->l ^ 0xFFFFFFFF) + 1;
 371             for (i = 0; GET_SIGN_L (fac->l) == 0; i++) fac->l = fac->l << 1;
 372             exp = ((FPS & FPS_L)? FP_BIAS + 32: FP_BIAS + 16) - i;
 373             fac->h = (sign << FP_V_SIGN) | (exp << FP_V_EXP) |
 374                 ((fac->l >> (31 - FP_V_HB)) & FP_FRACH);
 375             fac->l = (fac->l << (FP_V_HB + 1)) & FP_FRACL;
 376             if ((FPS & (FPS_D + FPS_T)) == 0) return roundfp11 (fac);  }
 377         return 0;
 378 }
 379
 380 static int32_t movfifp11(val, dst) _DOUBLE *val; int32_t *dst; {
 381         int32_t i, qdouble, tolong;
 382         int32_t exp, sign;
 383         _DOUBLE fac, fsrc;
 384         static const uint32_t i_limit[2][2] =
 385                 { { 0x80000000, 0x80010000 }, { 0x80000000, 0x80000001 } };
 386
 387         qdouble = FPS & FPS_D;
 388         sign = GET_SIGN (val->h);                       /* get sign, */
 389         exp = GET_EXP (val->h);                         /* exponent, */
 390         F_LOAD_FRAC (qdouble, (*val), fac);             /* fraction */
 391         if (FPS & FPS_L) {
 392             tolong = 1;
 393             i = FP_BIAS + 32;  }
 394         else {
 395             tolong = 0;
 396             i = FP_BIAS + 16;  }
 397         if (exp <= FP_BIAS) *dst = 0;
 398         else if (exp > i) {
 399             *dst = 0;
 400             fpnotrap (FEC_ICVT);
 401             return FPS_V;  }
 402         F_RSH_V (fac, FP_V_HB + 1 + i - exp, fsrc);
 403         if (!tolong) fsrc.l = fsrc.l & ~0177777;
 404         if (fsrc.l >= i_limit[tolong][sign]) {
 405             *dst = 0;
 406             fpnotrap (FEC_ICVT);
 407             return FPS_V;  }
 408         *dst = fsrc.l;
 409         if (sign) *dst = -*dst;
 410         return 0;
 411 }
 412
 413 /* Floating point add
 414
 415    Inputs:
 416         facp    =       pointer to src1 (output)
 417         fsrcp   =       pointer to src2
 418    Outputs:
 419         ovflo   =       overflow variable
 420 */
 421
 422 static int32_t addfp11(facp, fsrcp) _DOUBLE *facp; _DOUBLE *fsrcp; {
 423 int32_t facexp, fsrcexp, ediff;
 424 _DOUBLE facfrac, fsrcfrac;
 425
 426 if (F_LT_AP (facp, fsrcp)) {                            /* if !fac! < !fsrc! */
 427         facfrac = *facp;
 428         *facp = *fsrcp;                                 /* swap operands */
 429         *fsrcp = facfrac;  }
 430 facexp = GET_EXP (facp->h);                             /* get exponents */
 431 fsrcexp = GET_EXP (fsrcp->h);
 432 if (facexp == 0) {                                      /* fac = 0? */
 433         *facp = fsrcexp? *fsrcp: zero_fac;              /* result fsrc or 0 */
 434         return 0;  }
 435 if (fsrcexp == 0) return 0;                             /* fsrc = 0? no op */
 436 ediff = facexp - fsrcexp;                               /* exponent diff */
 437 if (ediff >= 60) return 0;                              /* too big? no op */
 438 F_GET_FRAC_P (facp, facfrac);                           /* get fractions */
 439 F_GET_FRAC_P (fsrcp, fsrcfrac);
 440 F_LSH_GUARD (facfrac);                                  /* guard fractions */
 441 F_LSH_GUARD (fsrcfrac);
 442 if (GET_SIGN (facp->h) != GET_SIGN (fsrcp->h)) {        /* signs different? */
 443         if (ediff) { F_RSH_V (fsrcfrac, ediff, fsrcfrac);  } /* sub, shf fsrc */
 444         F_SUB (fsrcfrac, facfrac, facfrac);             /* sub fsrc from fac */
 445         if ((facfrac.h | facfrac.l) == 0)  {            /* result zero? */
 446             *facp = zero_fac;                           /* no overflow */
 447             return 0;  }
 448         if (ediff <= 1) {                               /* big normalize? */
 449             if ((facfrac.h & (0x00FFFFFF << FP_GUARD)) == 0) {
 450                 F_LSH_K (facfrac, 24, facfrac);
 451                 facexp = facexp - 24;  }
 452             if ((facfrac.h & (0x00FFF000 << FP_GUARD)) == 0) {
 453                 F_LSH_K (facfrac, 12, facfrac);
 454                 facexp = facexp - 12;  }
 455             if ((facfrac.h & (0x00FC0000 << FP_GUARD)) == 0) {
 456                 F_LSH_K (facfrac, 6, facfrac);
 457                 facexp = facexp - 6;  }  }
 458         while (GET_BIT (facfrac.h, FP_V_HB + FP_GUARD) == 0) {
 459             F_LSH_1 (facfrac);
 460             facexp = facexp - 1;  }  }
 461 else {  if (ediff) { F_RSH_V (fsrcfrac, ediff, fsrcfrac);  } /* add, shf fsrc */
 462         F_ADD (fsrcfrac, facfrac, facfrac);             /* add fsrc to fac */
 463         if (GET_BIT (facfrac.h, FP_V_HB + FP_GUARD + 1)) {
 464             F_RSH_1 (facfrac);                          /* carry out, shift */
 465             facexp = facexp + 1;  }  }
 466 return round_and_pack (facp, facexp, &facfrac, 1);
 467 }
 468
 469 /* Floating point multiply
 470
 471    Inputs:
 472         facp    =       pointer to src1 (output)
 473         fsrcp   =       pointer to src2
 474    Outputs:
 475         ovflo   =       overflow indicator
 476 */
 477
 478 static int32_t mulfp11(facp, fsrcp) _DOUBLE *facp; _DOUBLE *fsrcp; {
 479 int32_t facexp, fsrcexp;
 480 _DOUBLE facfrac, fsrcfrac;
 481
 482 facexp = GET_EXP (facp->h);                             /* get exponents */
 483 fsrcexp = GET_EXP (fsrcp->h);
 484 if ((facexp == 0) || (fsrcexp == 0)) {                  /* test for zero */
 485         *facp = zero_fac;
 486         return 0;  }
 487 F_GET_FRAC_P (facp, facfrac);                           /* get fractions */
 488 F_GET_FRAC_P (fsrcp, fsrcfrac);
 489 facexp = facexp + fsrcexp - FP_BIAS;                    /* calculate exp */
 490 facp->h = facp->h  ^ fsrcp->h;                          /* calculate sign */
 491 frac_mulfp11 (&facfrac, &fsrcfrac);                     /* multiply fracs */
 492
 493 /* Multiplying two numbers in the range [.5,1) produces a result in the
 494    range [.25,1).  Therefore, at most one bit of normalization is required
 495    to bring the result back to the range [.5,1).
 496 */
 497
 498 if (GET_BIT (facfrac.h, FP_V_HB + FP_GUARD) == 0) {
 499         F_LSH_1 (facfrac);
 500         facexp = facexp - 1;  }
 501 return round_and_pack (facp, facexp, &facfrac, 1);
 502 }
 503
 504 /* Fraction multiply
 505
 506    Inputs:
 507         f1p     =       pointer to multiplier (output)
 508         f2p     =       pointer to multiplicand fraction
 509
 510    Note: the inputs are unguarded; the output is guarded.
 511
 512    This routine performs a classic shift-and-add multiply.  The low
 513    order bit of the multiplier is tested; if 1, the multiplicand is
 514    added into the high part of the double precision result.  The
 515    result and the multiplier are both shifted right 1.
 516
 517    For the 24b x 24b case, this routine develops 48b of result.
 518    For the 56b x 56b case, this routine only develops the top 64b
 519    of the the result.  Because the inputs are normalized fractions,
 520    the interesting part of the result is the high 56+guard bits.
 521    Everything shifted off to the right, beyond 64b, plays no part
 522    in rounding or the result.
 523
 524    There are many possible optimizations in this routine: scanning
 525    for groups of zeroes, particularly in the 56b x 56b case; using
 526    "extended multiply" capability if available in the hardware.
 527 */
 528
 529 static void frac_mulfp11(f1p, f2p) _DOUBLE *f1p; _DOUBLE *f2p; {
 530 _DOUBLE result, mpy, mpc;
 531 int32_t i;
 532
 533 result = zero_fac;                                      /* clear result */
 534 mpy = *f1p;                                             /* get operands */
 535 mpc = *f2p;
 536 F_LSH_GUARD (mpc);                                      /* guard multipicand */
 537 if ((mpy.l | mpc.l) == 0) {                             /* 24b x 24b? */
 538         for (i = 0; i < 24; i++) {
 539             if (mpy.h & 1) result.h = result.h + mpc.h;
 540             F_RSH_1 (result);
 541             mpy.h = mpy.h >> 1;  }  }
 542 else {  if (mpy.l != 0) {                               /* 24b x 56b? */
 543             for (i = 0; i < 32; i++) {
 544                 if (mpy.l & 1) { F_ADD (mpc, result, result);  }
 545                 F_RSH_1 (result);
 546                 mpy.l = mpy.l >> 1;  }  }
 547         for (i = 0; i < 24; i++) {
 548             if (mpy.h & 1) { F_ADD (mpc, result, result);  }
 549             F_RSH_1 (result);
 550             mpy.h = mpy.h >> 1;  }  }
 551 *f1p = result;
 552 return;
 553 }
 554
 555 /* Floating point divide
 556
 557    Inputs:
 558         facp    =       pointer to dividend (output)
 559         fsrcp   =       pointer to divisor
 560    Outputs:
 561         ovflo   =       overflow indicator
 562
 563    Source operand must be checked for zero by caller!
 564 */
 565
 566 static int32_t divfp11(facp, fsrcp) _DOUBLE *facp; _DOUBLE *fsrcp; {
 567 int32_t facexp, fsrcexp, i, count, qd;
 568 _DOUBLE facfrac, fsrcfrac, quo;
 569
 570 fsrcexp = GET_EXP (fsrcp->h);                           /* get divisor exp */
 571 facexp = GET_EXP (facp->h);                             /* get dividend exp */
 572 if (facexp == 0) {                                      /* test for zero */
 573         *facp = zero_fac;                               /* result zero */
 574         return 0;  }
 575 F_GET_FRAC_P (facp, facfrac);                           /* get fractions */
 576 F_GET_FRAC_P (fsrcp, fsrcfrac);
 577 F_LSH_GUARD (facfrac);                                  /* guard fractions */
 578 F_LSH_GUARD (fsrcfrac);
 579 facexp = facexp - fsrcexp + FP_BIAS + 1;                /* calculate exp */
 580 facp->h = facp->h ^ fsrcp->h;                           /* calculate sign */
 581 qd = FPS & FPS_D;
 582 count = FP_V_HB + FP_GUARD + (qd? 33: 1);               /* count = 56b/24b */
 583
 584 quo = zero_fac;
 585 for (i = count; (i > 0) && ((facfrac.h | facfrac.l) != 0); i--) {
 586         F_LSH_1 (quo);                                  /* shift quotient */
 587         if (!F_LT (facfrac, fsrcfrac)) {                /* divd >= divr? */
 588             F_SUB (fsrcfrac, facfrac, facfrac);         /* divd - divr */
 589             if (qd) quo.l = quo.l | 1;                  /* double or single? */
 590             else quo.h = quo.h | 1;  }
 591         F_LSH_1 (facfrac);  }                           /* shift divd */
 592 if (i > 0) { F_LSH_V (quo, i, quo);  }                  /* early exit? */
 593
 594 /* Dividing two numbers in the range [.5,1) produces a result in the
 595    range [.5,2).  Therefore, at most one bit of normalization is required
 596    to bring the result back to the range [.5,1).  The choice of counts
 597    and quotient bit positions makes this work correctly.
 598 */
 599
 600 if (GET_BIT (quo.h, FP_V_HB + FP_GUARD) == 0) {
 601         F_LSH_1 (quo);
 602         facexp = facexp - 1;  }
 603 return round_and_pack (facp, facexp, &quo, 1);
 604 }
 605
 606 /* Round (in place) floating point number to f_floating
 607
 608    Inputs:
 609         fptr    =       pointer to floating number
 610    Outputs:
 611         ovflow  =       overflow
 612 */
 613
 614 static int32_t roundfp11(fptr) _DOUBLE *fptr; {
 615 _DOUBLE outf;
 616
 617 outf = *fptr;                                           /* get argument */
 618 F_ADD (fround_fac, outf, outf);                         /* round */
 619 if (GET_SIGN (outf.h ^ fptr->h)) {                      /* flipped sign? */
 620         outf.h = (outf.h ^ FP_SIGN) & 0xFFFFFFFF;       /* restore sign */
 621         if (fpnotrap (FEC_OVFLO)) *fptr = zero_fac;     /* if no int, clear */
 622         else *fptr = outf;                              /* return rounded */
 623         return FPS_V;  }                                /* overflow */
 624 else {  *fptr = outf;                                   /* round was ok */
 625         return 0;  }                                    /* no overflow */
 626 }
 627
 628 /* Round result of calculation, test overflow, pack
 629
 630    Input:
 631         facp    =       pointer to result, sign in place
 632         exp     =       result exponent, right justified
 633         fracp   =       pointer to result fraction, right justified with
 634                         guard bits
 635         r       =       round (1) or truncate (0)
 636    Outputs:
 637         ovflo   =       overflow indicator
 638 */
 639
 640 static int32_t round_and_pack(facp, exp, fracp, r) _DOUBLE *facp; int32_t exp; _DOUBLE *fracp; int r; {
 641 _DOUBLE frac;
 642
 643 frac = *fracp;                                          /* get fraction */
 644 if (r && ((FPS & FPS_T) == 0)) {
 645         if (FPS & FPS_D) { F_ADD (dround_guard_fac, frac, frac);  }
 646         else { F_ADD (fround_guard_fac, frac, frac);  }
 647         if (GET_BIT (frac.h, FP_V_HB + FP_GUARD + 1)) {
 648             F_RSH_1 (frac);
 649             exp = exp + 1;  }  }
 650 F_RSH_GUARD (frac);
 651 facp->l = frac.l & FP_FRACL;
 652 facp->h = (facp->h & FP_SIGN) | ((exp & FP_M_EXP) << FP_V_EXP) |
 653         (frac.h & FP_FRACH);
 654 if (exp > 0377) {
 655         if (fpnotrap (FEC_OVFLO)) *facp = zero_fac;
 656         return FPS_V;  }
 657 if ((exp <= 0) && (fpnotrap (FEC_UNFLO))) *facp = zero_fac;
 658 return 0;
 659 }
 660
 661 /* Process floating point exception
 662
 663    Inputs:
 664         code    =       exception code
 665    Outputs:
 666         int     =       FALSE if interrupt enabled, TRUE if disabled
 667 */
 668
 669 static int32_t fpnotrap(code) int32_t code; {
 670 static const int32_t test_code[] = { 0, 0, 0, FPS_IC, FPS_IV, FPS_IU, FPS_IUV };
 671
 672 if ((code >= FEC_ICVT) && (code <= FEC_UNDFV) &&
 673     ((FPS & test_code[code >> 1]) == 0)) return TRUE;
 674 FPS = FPS | FPS_ER;
 675 FEC = code;
 676 /*FEA = (backup_PC - 2) & 0177777;*/
 677 /*if ((FPS & FPS_ID) == 0) setTRAP (TRAP_FPE);*/
 678 return FALSE;
 679 }
 680 #endif