1  /*
       2    Program to test complex divide for correct results on selected values.
       3    Checking known failure points.
       4  */
       5  
       6  #include <float.h>
       7  
       8  extern void abort (void);
       9  extern void exit (int);
      10  
      11  extern int ilogb (double);
      12  int match (double _Complex, double _Complex);
      13  
      14  #define SMALL DBL_MIN
      15  #define MAXBIT DBL_MANT_DIG
      16  #define ERRLIM 6
      17  
      18  /*
      19    Compare c (computed value) with z (expected value).
      20    Return 0 if within allowed range.  Return 1 if not.
      21  */
      22  int match (double _Complex c, double _Complex z)
      23  {
      24    double rz, iz, rc, ic;
      25    double rerr, ierr, rmax;
      26    int biterr;
      27    rz = __real__ z;
      28    iz = __imag__ z;
      29    rc = __real__ c;
      30    ic = __imag__ c;
      31  
      32    if (__builtin_fabs (rz) > SMALL)
      33      {
      34        rerr = __builtin_fabs (rz - rc) / __builtin_fabs (rz);
      35      }
      36    else if (__builtin_fabs (rz) == 0.0)
      37      {
      38        rerr = __builtin_fabs (rc);
      39      }
      40    else
      41      {
      42        rerr = __builtin_fabs (rz - rc) / SMALL;
      43      }
      44  
      45    if (__builtin_fabs (iz) > SMALL)
      46      {
      47        ierr = __builtin_fabs (iz - ic) / __builtin_fabs (iz);
      48      }
      49    else if (__builtin_fabs (iz) == 0.0)
      50      {
      51        ierr = __builtin_fabs (ic);
      52      }
      53    else
      54      {
      55        ierr = __builtin_fabs (iz - ic) / SMALL;
      56      }
      57    rmax = __builtin_fmax(rerr, ierr);
      58    biterr = 0;
      59    if ( rmax != 0.0)      
      60      {
      61        biterr = ilogb (rmax) + MAXBIT + 1;
      62      }
      63  
      64    if (biterr >= ERRLIM)
      65      return 0;
      66    else
      67      return 1;
      68  }
      69  
      70  
      71  int main (int argc, char** argv)
      72  {
      73    double _Complex a,b,c,z;
      74    double xr[4], xi[4], yr[4], yi[4], zr[4], zi[4];
      75    double cr, ci;
      76    int i;
      77    int ok = 1;
      78    xr[0] = -0x1.16e7fad79e45ep+651;
      79    xi[0] = -0x1.f7f75b94c6c6ap-860;
      80    yr[0] = -0x1.2f40d8ff7e55ep+245;
      81    yi[0] = -0x0.0000000004ebcp-1022;
      82    zr[0] = 0x1.d6e4b0e282869p+405;
      83    zi[0] = -0x1.e9095e311e706p-900;
      84  
      85    xr[1] = -0x1.21ff587f953d3p-310;
      86    xi[1] = -0x1.5a526dcc59960p+837;
      87    yr[1] = 0x1.b88b8b552eaadp+735;
      88    yi[1] = -0x1.873e2d6544d92p-327;
      89    zr[1] = 0x1.65734a88b2de0p-961;
      90    zi[1] =  -0x1.927e85b8b5770p+101;
      91  
      92    xr[2] = 0x1.4612e41aa8080p-846;
      93    xi[2] = -0x0.0000000613e07p-1022;
      94    yr[2] = 0x1.df9cd0d58caafp-820;
      95    yi[2] = -0x1.e47051a9036dbp-584;
      96    zr[2] = 0x1.9b194f3fffa32p-469;
      97    zi[2] = 0x1.58a00ab740a6bp-263;
      98  
      99    xr[3] = 0x1.cb27eece7c585p-355;
     100    xi[3] = 0x0.000000223b8a8p-1022;
     101    yr[3] = -0x1.74e7ed2b9189fp-22;
     102    yi[3] = 0x1.3d80439e9a119p-731;
     103    zr[3] = -0x1.3b35ed806ae5ap-333;
     104    zi[3] = -0x0.05e01bcbfd9f6p-1022;
     105  
     106  
     107    for (i = 0; i < 4; i++)
     108      {
     109        __real__ a = xr[i];
     110        __imag__ a = xi[i];
     111        __real__ b = yr[i];
     112        __imag__ b = yi[i];
     113        __real__ z = zr[i];
     114        __imag__ z = zi[i];
     115        c = a / b;
     116        cr = __real__ c;
     117        ci = __imag__ c;
     118  
     119        if (!match (c,z)){
     120  	ok = 0;
     121        }
     122      }
     123    if (!ok)
     124      abort ();
     125    exit (0);
     126  }