workbench/libs/mathieeedoubbas/ieeedpmul.c

   1 /*
   2     Copyright © 1995-2017, The AROS Development Team. All rights reserved.
   3     $Id$
   4 */
   5 #include <aros/debug.h>
   6 #include <proto/utility.h>
   7 #include "mathieeedoubbas_intern.h"
   8
   9 static void add128(ULONG *s, ULONG *d)
  10 {
  11     ULONG ovl = 0;
  12     WORD i;
  13
  14     for (i = 3; i >= 0; i--) {
  15        ULONG o = d[i];
  16        d[i] += s[i] + ovl;
  17        ovl = 0;
  18        if (o > d[i])
  19            ovl = 1;
  20     }
  21 }
  22
  23 /*****************************************************************************
  24
  25     NAME */
  26
  27         AROS_LHQUAD2(double, IEEEDPMul,
  28
  29 /*  SYNOPSIS */
  30         AROS_LHAQUAD(double, y, D0, D1),
  31         AROS_LHAQUAD(double, z, D2, D3),
  32
  33 /*  LOCATION */
  34         struct MathIeeeDoubBasBase *, MathIeeeDoubBasBase, 13, MathIeeeDoubBas)
  35
  36 /*  FUNCTION
  37         Multiplies two IEEE double precision numbers.
  38
  39     INPUTS
  40         y - first multiplicand.
  41         z - second multiplicand.
  42
  43     RESULT
  44         x - product.
  45
  46     NOTES
  47
  48     EXAMPLE
  49
  50     BUGS
  51
  52     SEE ALSO
  53
  54     INTERNALS
  55
  56 *****************************************************************************/
  57 {
  58     AROS_LIBFUNC_INIT
  59
  60     BOOL sign;
  61     QUAD res;
  62     QUAD Qtmp1, Qtmp2;
  63     WORD exp1, exp2, exp;
  64     double *Dres = (double *)&res;
  65
  66     QUAD * Qy = (QUAD *)&y;
  67     QUAD * Qz = (QUAD *)&z;
  68
  69     D(bug("%08x %08x * %08x %08x\n",
  70         (ULONG)Get_High32of64(*Qy), (ULONG)Get_Low32of64(*Qy),
  71         (ULONG)Get_High32of64(*Qz), (ULONG)Get_Low32of64(*Qz)));
  72
  73     sign = ((Get_High32of64(*Qy) ^ Get_High32of64(*Qz)) & IEEEDPSign_Mask_Hi) != 0;
  74
  75     AND64C(Qtmp1, *Qy, IEEEDPExponent_Mask_Hi, IEEEDPExponent_Mask_Lo);
  76     AND64C(Qtmp2, *Qz, IEEEDPExponent_Mask_Hi, IEEEDPExponent_Mask_Lo);
  77
  78     SHRU32(exp1, Qtmp1, 52);
  79     SHRU32(exp2, Qtmp2, 52);
  80
  81     exp1 = (exp1 & 0x7ff) - 1023;
  82     exp2 = (exp2 & 0x7ff) - 1023;
  83     exp = exp1 + exp2;
  84
  85     D(bug("EXP %d + EXP %d = %d (%08x)\n", exp1, exp2, exp, (exp + 1023) << 20));
  86
  87     AND64C(Qtmp1, *Qy, IEEEDPMantisse_Mask_Hi, IEEEDPMantisse_Mask_Lo);
  88     AND64C(Qtmp2, *Qz, IEEEDPMantisse_Mask_Hi, IEEEDPMantisse_Mask_Lo);
  89
  90     /* (53 bit) * (53 bit) multiplication */
  91
  92     ULONG x1, x2, y1, y2;
  93     ULONG t1[4], t2[4];
  94
  95     x1 = Get_High32of64(Qtmp1);
  96     x1 |= 0x00100000; /* Set "hidden" 53rd mantissa bit */
  97     x2 = Get_Low32of64(Qtmp1);
  98     y1 = Get_High32of64(Qtmp2);
  99     y1 |= 0x00100000; /* Set "hidden" 53rd mantissa bit */
 100     y2 = Get_Low32of64(Qtmp2);
 101
 102     D(bug("%08x %08x %08x %08x\n", x1, x2, y1, y2));
 103
 104     Qtmp1 = UMult64(x1, y1);
 105     t1[0] = Get_High32of64(Qtmp1);
 106     t1[1] = Get_Low32of64(Qtmp1);
 107
 108     Qtmp1 = UMult64(x2, y2);
 109     t1[2] = Get_High32of64(Qtmp1);
 110     t1[3] = Get_Low32of64(Qtmp1);
 111
 112     Qtmp1 = UMult64(x1, y2);
 113     t2[3] = 0;
 114     t2[2] = Get_High32of64(Qtmp1);
 115     t2[1] = Get_Low32of64(Qtmp1);
 116     t2[0] = 0;
 117
 118     add128(t2, t1);
 119
 120     Qtmp1 = UMult64(x2, y1);
 121     t2[3] = 0;
 122     t2[2] = Get_High32of64(Qtmp1);
 123     t2[1] = Get_Low32of64(Qtmp1);
 124     t2[0] = 0;
 125
 126     add128(t2, t1);
 127
 128     D(bug("%08x %08x %08x %08x\n", t1[0], t1[1], t1[2], t1[3]));
 129
 130     /* Normalize. Probably could be more optimal */
 131     if (t1[0] || t1[1] || t1[2] || t1[3]) {
 132         while (t1[0] & 0xfffffe00) {
 133             t1[3] >>= 1;
 134             t1[3] |= (t1[2] & 1) ? 0x80000000 : 0;
 135             t1[2] >>= 1;
 136             t1[2] |= (t1[1] & 1) ? 0x80000000 : 0;
 137             t1[1] >>= 1;
 138             t1[1] |= (t1[0] & 1) ? 0x80000000 : 0;
 139             t1[0] >>= 1;
 140             exp++;
 141         }
 142         while (!(t1[0] & 0x00000100)) {
 143             t1[0] <<= 1;
 144             t1[0] |= (t1[1] & 0x80000000) ? 1 : 0;
 145             t1[1] <<= 1;
 146             t1[1] |= (t1[2] & 0x80000000) ? 1 : 0;
 147             t1[2] <<= 1;
 148             t1[2] |= (t1[3] & 0x80000000) ? 1 : 0;
 149             t1[3] <<= 1;
 150             exp--;
 151         }
 152     }
 153
 154     D(bug("%08x %08x %08x %08x\n", t1[0], t1[1], t1[2], t1[3]));
 155
 156     if (exp <= -1023) {
 157         if (sign)
 158             Set_Value64C(res, 0x0, 0x0);
 159         else
 160             Set_Value64C(res, IEEEDPSign_Mask_Hi, IEEEDPSign_Mask_Lo);
 161         return *Dres;
 162     }
 163
 164     if (exp >= 1023) {
 165         if (sign)
 166             Set_Value64C(res, IEEEDPPInfty_Hi, IEEEDPPInfty_Lo);
 167         else
 168             Set_Value64C(res, IEEEDPPInfty_Hi | IEEEDPSign_Mask_Hi, IEEEDPPInfty_Lo | IEEEDPSign_Mask_Lo);
 169         return *Dres;
 170     }
 171
 172     Set_Value64C(res, (t1[0] << 12) | (t1[1] >> 20), (t1[1] << 12) | (t1[2] >> 20));
 173     AND64QC(res, IEEEDPMantisse_Mask_Hi, IEEEDPMantisse_Mask_Lo);
 174     SHL64(Qtmp1, (QUAD)(exp + 1023), 52);
 175     OR64Q(res, Qtmp1);
 176
 177     if (sign)
 178         OR64(res, IEEEDPSign_Mask_Hi, IEEEDPSign_Mask_Lo);
 179
 180     return *Dres;
 181
 182     AROS_LIBFUNC_EXIT
 183 }