interface/zrotg.c

   1 #include <math.h>
   2 #include "common.h"
   3 #ifdef FUNCTION_PROFILE
   4 #include "functable.h"
   5 #endif
   6
   7 void NAME(FLOAT *DA, FLOAT *DB, FLOAT *C, FLOAT *S){
   8
   9 #if defined(__i386__) || defined(__x86_64__) || defined(__ia64__) || defined(_M_X64) || defined(_M_IX86)
  10
  11   long double da_r = *(DA + 0);
  12   long double da_i = *(DA + 1);
  13   long double db_r = *(DB + 0);
  14   long double db_i = *(DB + 1);
  15   long double r;
  16
  17   long double ada = fabsl(da_r) + fabsl(da_i);
  18
  19   PRINT_DEBUG_NAME;
  20
  21   IDEBUG_START;
  22
  23   FUNCTION_PROFILE_START();
  24
  25   if (ada == ZERO) {
  26     *C        = ZERO;
  27     *(S  + 0) = ONE;
  28     *(S  + 1) = ZERO;
  29     *(DA + 0) = db_r;
  30     *(DA + 1) = db_i;
  31   } else {
  32     long double alpha_r, alpha_i;
  33
  34     ada = sqrt(da_r * da_r + da_i * da_i);
  35
  36     r = sqrt(da_r * da_r + da_i * da_i + db_r * db_r + db_i * db_i);
  37
  38     alpha_r = da_r / ada;
  39     alpha_i = da_i / ada;
  40
  41     *(C + 0)  = ada / r;
  42     *(S + 0)  = (alpha_r * db_r + alpha_i *db_i) / r;
  43     *(S + 1)  = (alpha_i * db_r - alpha_r *db_i) / r;
  44     *(DA + 0) = alpha_r * r;
  45     *(DA + 1) = alpha_i * r;
  46   }
  47 #else
  48   FLOAT da_r = *(DA + 0);
  49   FLOAT da_i = *(DA + 1);
  50   FLOAT db_r = *(DB + 0);
  51   FLOAT db_i = *(DB + 1);
  52   FLOAT r;
  53
  54   FLOAT ada = fabs(da_r) + fabs(da_i);
  55   FLOAT adb;
  56
  57   PRINT_DEBUG_NAME;
  58
  59   IDEBUG_START;
  60
  61   FUNCTION_PROFILE_START();
  62
  63   if (ada == ZERO) {
  64     *C        = ZERO;
  65     *(S  + 0) = ONE;
  66     *(S  + 1) = ZERO;
  67     *(DA + 0) = db_r;
  68     *(DA + 1) = db_i;
  69   } else {
  70     FLOAT scale;
  71     FLOAT aa_r, aa_i, bb_r, bb_i;
  72     FLOAT alpha_r, alpha_i;
  73
  74     aa_r = fabs(da_r);
  75     aa_i = fabs(da_i);
  76
  77     if (aa_i > aa_r) {
  78       aa_r = fabs(da_i);
  79       aa_i = fabs(da_r);
  80     }
  81
  82     scale = (aa_i / aa_r);
  83     ada = aa_r * sqrt(ONE + scale * scale);
  84
  85     bb_r = fabs(db_r);
  86     bb_i = fabs(db_i);
  87
  88     if (bb_i > bb_r) {
  89       bb_r = fabs(bb_i);
  90       bb_i = fabs(bb_r);
  91     }
  92
  93     scale = (bb_i / bb_r);
  94     adb = bb_r * sqrt(ONE + scale * scale);
  95
  96     scale = ada + adb;
  97
  98     aa_r    = da_r / scale;
  99     aa_i    = da_i / scale;
 100     bb_r    = db_r / scale;
 101     bb_i    = db_i / scale;
 102
 103     r = scale * sqrt(aa_r * aa_r + aa_i * aa_i + bb_r * bb_r + bb_i * bb_i);
 104
 105     alpha_r = da_r / ada;
 106     alpha_i = da_i / ada;
 107
 108     *(C + 0)  = ada / r;
 109     *(S + 0)  = (alpha_r * db_r + alpha_i *db_i) / r;
 110     *(S + 1)  = (alpha_i * db_r - alpha_r *db_i) / r;
 111     *(DA + 0) = alpha_r * r;
 112     *(DA + 1) = alpha_i * r;
 113   }
 114 #endif
 115
 116   FUNCTION_PROFILE_END(4, 4, 4);
 117
 118   IDEBUG_END;
 119
 120   return;
 121 }