workbench/libs/mathieeedoubbas/ieeedpadd.c

   1 /*
   2     Copyright © 1995-2003, The AROS Development Team. All rights reserved.
   3     $Id$
   4 */
   5
   6 #include "mathieeedoubbas_intern.h"
   7
   8 /*****************************************************************************
   9
  10     NAME */
  11
  12         AROS_LHQUAD2(double, IEEEDPAdd,
  13
  14 /*  SYNOPSIS */
  15         AROS_LHAQUAD(double, y, D0, D1),
  16         AROS_LHAQUAD(double, z, D2, D3),
  17
  18 /*  LOCATION */
  19         struct MathIeeeDoubBasBase *, MathIeeeDoubBasBase, 11, MathIeeeDoubBas)
  20
  21 /*  FUNCTION
  22         Calculate the sum of two IEEE double precision numbers
  23
  24     INPUTS
  25
  26     RESULT
  27         sum of y and z
  28
  29         Flags:
  30           zero     : result is zero
  31           negative : result is negative
  32           overflow : result is too large or too small for ffp format
  33
  34     BUGS
  35
  36     INTERNALS
  37
  38 *****************************************************************************/
  39 {
  40     AROS_LIBFUNC_INIT
  41
  42     QUAD Res;
  43     QUAD Mant1, Mant2;
  44     QUAD Qtmp1, Qtmp2;
  45     LONG Shift;
  46     QUAD Exponent;
  47     LONG Ltmp1, Ltmp2;
  48     QUAD * Qy = (QUAD *)&y;
  49     QUAD * Qz = (QUAD *)&z;
  50     double * DRes = (double *)&Res;
  51
  52     SetSR(0, Zero_Bit | Overflow_Bit | Negative_Bit );
  53
  54     AND64C(Qtmp1, *Qy, IEEEDPExponent_Mask_Hi, IEEEDPExponent_Mask_Lo);
  55     AND64C(Qtmp2, *Qz, IEEEDPExponent_Mask_Hi, IEEEDPExponent_Mask_Lo);
  56     SHRU32(Ltmp1, Qtmp1, 52);
  57     SHRU32(Ltmp2, Qtmp2, 52);
  58     Shift = Ltmp1 - Ltmp2;
  59
  60     if
  61     (
  62            is_neqC(*Qy, 0x0, 0x0)
  63         && is_neqC(*Qy, IEEEDPSign_Mask_Hi, IEEEDPSign_Mask_Lo)
  64     )
  65     {
  66         AND64C(Mant1, *Qy, IEEEDPMantisse_Mask_Hi, IEEEDPMantisse_Mask_Lo);
  67         OR64QC(Mant1, 0x00100000, 0x00000000);
  68     }
  69     else
  70     {
  71         Set_Value64C(Mant1, 0x0, 0x0);
  72     }
  73
  74     if
  75     (
  76            is_neqC(*Qz, 0x0, 0x0)
  77         && is_neqC(*Qz, IEEEDPSign_Mask_Hi, IEEEDPSign_Mask_Lo)
  78     )
  79     {
  80         /* Mant2 = (*Qz & IEEESPMantisse_Mask) | 0x00800000; */
  81         AND64C(Mant2, *Qz, IEEEDPMantisse_Mask_Hi, IEEEDPMantisse_Mask_Lo);
  82         OR64QC(Mant2, 0x00100000, 0x00000000);
  83     }
  84     else
  85     {
  86         Set_Value64C(Mant2, 0x0, 0x0);
  87     }
  88
  89     if (Shift > 0)
  90     {
  91         SHRU64(Mant2, Mant2, Shift);
  92         AND64C(Exponent, *Qy, IEEEDPExponent_Mask_Hi, IEEEDPExponent_Mask_Lo);
  93     }
  94     else
  95     {
  96         SHRU64(Mant1, Mant1, (-Shift));
  97         AND64C(Exponent, *Qz, IEEEDPExponent_Mask_Hi, IEEEDPExponent_Mask_Lo);
  98     }
  99
 100     // sign(fnum1) == sign(fnum2)
 101     // simple addition
 102     // 0.25 <= res < 1
 103     AND64C(Qtmp1, *Qz, IEEEDPSign_Mask_Hi, IEEEDPSign_Mask_Lo);
 104     AND64C(Qtmp2, *Qy, IEEEDPSign_Mask_Hi, IEEEDPSign_Mask_Lo);
 105     if ( is_eq(Qtmp1, Qtmp2))
 106     {
 107         AND64C(Res, *Qy, IEEEDPSign_Mask_Hi, IEEEDPSign_Mask_Lo);
 108         ADD64Q(Mant1, Mant2);
 109         AND64C(Qtmp1, Mant1, 0x00200000, 0x0);
 110         if ( is_neqC(Qtmp1, 0x0, 0x0))
 111         {
 112             ADD64QC(Exponent, 0x00100000, 0x0);
 113             SHRU64(Mant1, Mant1, 1);
 114         }
 115         AND64QC(Mant1, IEEEDPMantisse_Mask_Hi, IEEEDPMantisse_Mask_Lo);
 116     }
 117     // second case: sign(fnum1) != sign(fnum2)
 118     // -1 <= res < 1
 119     else
 120     {
 121         //printf("Exponent: %x\n",Exponent);
 122         if ( is_lessSC(*Qy, 0x0, 0x0))
 123         {
 124             SUB64(Mant1, Mant2, Mant1);
 125         }
 126         else // fnum2 < 0
 127         {
 128             SUB64(Mant1, Mant1, Mant2);
 129         }
 130
 131         //if the result is below zero
 132         if (  is_lessSC(Mant1, 0x0, 0x0))
 133         {
 134             Set_Value64C(Res, IEEEDPSign_Mask_Hi, IEEEDPSign_Mask_Lo);
 135             NEG64(Mant1);
 136             SetSR(Negative_Bit, Zero_Bit | Negative_Bit | Overflow_Bit);
 137         }
 138         else
 139         {
 140             Set_Value64C(Res, 0x0, 0x0);
 141         }
 142
 143         if (is_eqC(Mant1, 0x0, 0x0))
 144         {
 145             union { QUAD i; double d; } tmp;
 146             SetSR(Zero_Bit, Zero_Bit | Negative_Bit | Overflow_Bit);
 147             Set_Value64C(Res, 0x0, 0x0);
 148             tmp.i = Res;
 149             return tmp.d;
 150         }
 151         else
 152         {
 153             /* normalize the mantisse */
 154             AND64C(Qtmp1, Mant1, 0x00100000, 0x0);
 155             while ( is_eqC(Qtmp1, 0x0, 0x0))
 156             {
 157                 SHL64(Mant1, Mant1, 1);  //one bit to the left.
 158                 SUB64QC(Exponent, 0x00100000, 0x0);
 159                 AND64C(Qtmp1, Mant1, 0x00100000, 0x0);
 160             }
 161             AND64QC(Mant1, IEEEDPMantisse_Mask_Hi, IEEEDPMantisse_Mask_Lo);
 162         } /* else */
 163
 164     } /* else */
 165
 166     if ( is_lessSC(Exponent, 0x0, 0x0))
 167     {
 168         SetSR(Overflow_Bit, Zero_Bit | Overflow_Bit); //do not change Negative_Bit!
 169         OR64QC
 170         (
 171             Res,
 172             (IEEEDPMantisse_Mask_Hi | IEEEDPExponent_Mask_Hi),
 173             (IEEEDPMantisse_Mask_Lo | IEEEDPExponent_Mask_Lo)
 174         );
 175     }
 176     OR64Q(Res, Mant1);
 177     OR64Q(Res, Exponent);
 178
 179     return *DRes;
 180
 181     AROS_LIBFUNC_EXIT
 182 }