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