workbench/libs/mathieeedoubbas/ieeedpadd.c

   1 /*
   2     Copyright © 1995-2017, 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         y - first addend.
  26         z - first addend.
  27
  28     RESULT
  29         x - sum of y and z.
  30
  31         Flags:
  32           zero     : result is zero
  33           negative : result is negative
  34           overflow : result is too large or too small for IEEE DP format
  35
  36     NOTES
  37
  38     EXAMPLE
  39
  40     BUGS
  41
  42     SEE ALSO
  43
  44     INTERNALS
  45
  46 *****************************************************************************/
  47 {
  48     AROS_LIBFUNC_INIT
  49
  50     QUAD Res;
  51     QUAD Mant1, Mant2;
  52     QUAD Qtmp1, Qtmp2;
  53     LONG Shift;
  54     QUAD Exponent;
  55     LONG Ltmp1, Ltmp2;
  56     QUAD * Qy = (QUAD *)&y;
  57     QUAD * Qz = (QUAD *)&z;
  58     double * DRes = (double *)&Res;
  59
  60     SetSR(0, Zero_Bit | Overflow_Bit | Negative_Bit );
  61
  62     AND64C(Qtmp1, *Qy, IEEEDPExponent_Mask_Hi, IEEEDPExponent_Mask_Lo);
  63     AND64C(Qtmp2, *Qz, IEEEDPExponent_Mask_Hi, IEEEDPExponent_Mask_Lo);
  64     SHRU32(Ltmp1, Qtmp1, 52);
  65     SHRU32(Ltmp2, Qtmp2, 52);
  66     Shift = Ltmp1 - Ltmp2;
  67
  68     if
  69     (
  70            is_neqC(*Qy, 0x0, 0x0)
  71         && is_neqC(*Qy, IEEEDPSign_Mask_Hi, IEEEDPSign_Mask_Lo)
  72     )
  73     {
  74         AND64C(Mant1, *Qy, IEEEDPMantisse_Mask_Hi, IEEEDPMantisse_Mask_Lo);
  75         OR64QC(Mant1, 0x00100000, 0x00000000);
  76     }
  77     else
  78     {
  79         Set_Value64C(Mant1, 0x0, 0x0);
  80     }
  81
  82     if
  83     (
  84            is_neqC(*Qz, 0x0, 0x0)
  85         && is_neqC(*Qz, IEEEDPSign_Mask_Hi, IEEEDPSign_Mask_Lo)
  86     )
  87     {
  88         /* Mant2 = (*Qz & IEEESPMantisse_Mask) | 0x00800000; */
  89         AND64C(Mant2, *Qz, IEEEDPMantisse_Mask_Hi, IEEEDPMantisse_Mask_Lo);
  90         OR64QC(Mant2, 0x00100000, 0x00000000);
  91     }
  92     else
  93     {
  94         Set_Value64C(Mant2, 0x0, 0x0);
  95     }
  96
  97     if (Shift > 0)
  98     {
  99         SHRU64(Mant2, Mant2, Shift);
 100         AND64C(Exponent, *Qy, IEEEDPExponent_Mask_Hi, IEEEDPExponent_Mask_Lo);
 101     }
 102     else
 103     {
 104         SHRU64(Mant1, Mant1, (-Shift));
 105         AND64C(Exponent, *Qz, IEEEDPExponent_Mask_Hi, IEEEDPExponent_Mask_Lo);
 106     }
 107
 108     // sign(fnum1) == sign(fnum2)
 109     // simple addition
 110     // 0.25 <= res < 1
 111     AND64C(Qtmp1, *Qz, IEEEDPSign_Mask_Hi, IEEEDPSign_Mask_Lo);
 112     AND64C(Qtmp2, *Qy, IEEEDPSign_Mask_Hi, IEEEDPSign_Mask_Lo);
 113     if ( is_eq(Qtmp1, Qtmp2))
 114     {
 115         AND64C(Res, *Qy, IEEEDPSign_Mask_Hi, IEEEDPSign_Mask_Lo);
 116         ADD64Q(Mant1, Mant2);
 117         AND64C(Qtmp1, Mant1, 0x00200000, 0x0);
 118         if ( is_neqC(Qtmp1, 0x0, 0x0))
 119         {
 120             ADD64QC(Exponent, 0x00100000, 0x0);
 121             SHRU64(Mant1, Mant1, 1);
 122         }
 123         AND64QC(Mant1, IEEEDPMantisse_Mask_Hi, IEEEDPMantisse_Mask_Lo);
 124     }
 125     // second case: sign(fnum1) != sign(fnum2)
 126     // -1 <= res < 1
 127     else
 128     {
 129         //printf("Exponent: %x\n",Exponent);
 130         if ( is_lessSC(*Qy, 0x0, 0x0))
 131         {
 132             SUB64(Mant1, Mant2, Mant1);
 133         }
 134         else // fnum2 < 0
 135         {
 136             SUB64(Mant1, Mant1, Mant2);
 137         }
 138
 139         //if the result is below zero
 140         if (  is_lessSC(Mant1, 0x0, 0x0))
 141         {
 142             Set_Value64C(Res, IEEEDPSign_Mask_Hi, IEEEDPSign_Mask_Lo);
 143             NEG64(Mant1);
 144             SetSR(Negative_Bit, Zero_Bit | Negative_Bit | Overflow_Bit);
 145         }
 146         else
 147         {
 148             Set_Value64C(Res, 0x0, 0x0);
 149         }
 150
 151         if (is_eqC(Mant1, 0x0, 0x0))
 152         {
 153             union { QUAD i; double d; } tmp;
 154             SetSR(Zero_Bit, Zero_Bit | Negative_Bit | Overflow_Bit);
 155             Set_Value64C(Res, 0x0, 0x0);
 156             tmp.i = Res;
 157             return tmp.d;
 158         }
 159         else
 160         {
 161             /* normalize the mantisse */
 162             AND64C(Qtmp1, Mant1, 0x00100000, 0x0);
 163             while ( is_eqC(Qtmp1, 0x0, 0x0))
 164             {
 165                 SHL64(Mant1, Mant1, 1);  //one bit to the left.
 166                 SUB64QC(Exponent, 0x00100000, 0x0);
 167                 AND64C(Qtmp1, Mant1, 0x00100000, 0x0);
 168             }
 169             AND64QC(Mant1, IEEEDPMantisse_Mask_Hi, IEEEDPMantisse_Mask_Lo);
 170         } /* else */
 171
 172     } /* else */
 173
 174     if ( is_lessSC(Exponent, 0x0, 0x0))
 175     {
 176         SetSR(Overflow_Bit, Zero_Bit | Overflow_Bit); //do not change Negative_Bit!
 177         OR64QC
 178         (
 179             Res,
 180             (IEEEDPMantisse_Mask_Hi | IEEEDPExponent_Mask_Hi),
 181             (IEEEDPMantisse_Mask_Lo | IEEEDPExponent_Mask_Lo)
 182         );
 183     }
 184     OR64Q(Res, Mant1);
 185     OR64Q(Res, Exponent);
 186
 187     return *DRes;
 188
 189     AROS_LIBFUNC_EXIT
 190 }