src/compiler/dmd/cmath2.d

   1
   2 // Copyright (C) 2001-2003 by Digital Mars
   3 // All Rights Reserved
   4 // www.digitalmars.com
   5
   6 // Runtime support for complex arithmetic code generation
   7 // (for linux)
   8
   9 /*
  10  *  Modified by Sean Kelly for use with the D Runtime Project
  11  */
  12
  13 module rt.cmath2;
  14
  15 private import stdc.math;
  16
  17 extern (C):
  18
  19 /****************************
  20  * Multiply two complex floating point numbers, x and y.
  21  * Input:
  22  *      x.re    ST3
  23  *      x.im    ST2
  24  *      y.re    ST1
  25  *      y.im    ST0
  26  * Output:
  27  *      ST1     real part
  28  *      ST0     imaginary part
  29  */
  30
  31 void _Cmul()
  32 {
  33     // p.re = x.re * y.re - x.im * y.im;
  34     // p.im = x.im * y.re + x.re * y.im;
  35     asm
  36     {   naked                   ;
  37         fld     ST(1)           ; // x.re
  38         fmul    ST,ST(4)        ; // ST0 = x.re * y.re
  39
  40         fld     ST(1)           ; // y.im
  41         fmul    ST,ST(4)        ; // ST0 = x.im * y.im
  42
  43         fsubp   ST(1),ST        ; // ST0 = x.re * y.re - x.im * y.im
  44
  45         fld     ST(3)           ; // x.im
  46         fmul    ST,ST(3)        ; // ST0 = x.im * y.re
  47
  48         fld     ST(5)           ; // x.re
  49         fmul    ST,ST(3)        ; // ST0 = x.re * y.im
  50
  51         faddp   ST(1),ST        ; // ST0 = x.im * y.re + x.re * y.im
  52
  53         fxch    ST(4),ST        ;
  54         fstp    ST(0)           ;
  55         fxch    ST(4),ST        ;
  56         fstp    ST(0)           ;
  57         fstp    ST(0)           ;
  58         fstp    ST(0)           ;
  59
  60         ret                     ;
  61     }
  62 /+
  63     if (isnan(x) && isnan(y))
  64     {
  65         // Recover infinities that computed as NaN+ iNaN ...
  66         int recalc = 0;
  67         if ( isinf( a) || isinf( b) )
  68         {   // z is infinite
  69             // "Box" the infinity and change NaNs in the other factor to 0
  70             a = copysignl( isinf( a) ? 1.0 : 0.0, a);
  71             b = copysignl( isinf( b) ? 1.0 : 0.0, b);
  72             if (isnan( c)) c = copysignl( 0.0, c);
  73             if (isnan( d)) d = copysignl( 0.0, d);
  74             recalc = 1;
  75         }
  76         if (isinf(c) || isinf(d))
  77         {   // w is infinite
  78             // "Box" the infinity and change NaNs in the other factor to 0
  79             c = copysignl( isinf( c) ? 1.0 : 0.0, c);
  80             d = copysignl( isinf( d) ? 1.0 : 0.0, d);
  81             if (isnan( a)) a = copysignl( 0.0, a);
  82             if (isnan( b)) b = copysignl( 0.0, b);
  83             recalc = 1;
  84         }
  85         if (!recalc && (isinf(ac) || isinf(bd) ||
  86             isinf(ad) || isinf(bc)))
  87         {
  88             // Recover infinities from overflow by changing NaNs to 0 ...
  89             if (isnan( a)) a = copysignl( 0.0, a);
  90             if (isnan( b)) b = copysignl( 0.0, b);
  91             if (isnan( c)) c = copysignl( 0.0, c);
  92             if (isnan( d)) d = copysignl( 0.0, d);
  93             recalc = 1;
  94         }
  95         if (recalc)
  96         {
  97             x = INFINITY * (a * c - b * d);
  98             y = INFINITY * (a * d + b * c);
  99         }
 100     }
 101 +/
 102 }
 103
 104 /****************************
 105  * Divide two complex floating point numbers, x / y.
 106  * Input:
 107  *      x.re    ST3
 108  *      x.im    ST2
 109  *      y.re    ST1
 110  *      y.im    ST0
 111  * Output:
 112  *      ST1     real part
 113  *      ST0     imaginary part
 114  */
 115
 116 void _Cdiv()
 117 {
 118     real x_re, x_im;
 119     real y_re, y_im;
 120     real q_re, q_im;
 121     real r;
 122     real den;
 123
 124     asm
 125     {
 126         fstp    y_im    ;
 127         fstp    y_re    ;
 128         fstp    x_im    ;
 129         fstp    x_re    ;
 130     }
 131
 132     if (fabs(y_re) < fabs(y_im))
 133     {
 134         r = y_re / y_im;
 135         den = y_im + r * y_re;
 136         q_re = (x_re * r + x_im) / den;
 137         q_im = (x_im * r - x_re) / den;
 138     }
 139     else
 140     {
 141         r = y_im / y_re;
 142         den = y_re + r * y_im;
 143         q_re = (x_re + r * x_im) / den;
 144         q_im = (x_im - r * x_re) / den;
 145     }
 146 //printf("q.re = %g, q.im = %g\n", (double)q_re, (double)q_im);
 147 /+
 148     if (isnan(q_re) && isnan(q_im))
 149     {
 150         real denom = y_re * y_re + y_im * y_im;
 151
 152         // non-zero / zero
 153         if (denom == 0.0 && (!isnan(x_re) || !isnan(x_im)))
 154         {
 155             q_re = copysignl(INFINITY, y_re) * x_re;
 156             q_im = copysignl(INFINITY, y_re) * x_im;
 157         }
 158         // infinite / finite
 159         else if ((isinf(x_re) || isinf(x_im)) && isfinite(y_re) && isfinite(y_im))
 160         {
 161             x_re = copysignl(isinf(x_re) ? 1.0 : 0.0, x_re);
 162             x_im = copysignl(isinf(x_im) ? 1.0 : 0.0, x_im);
 163             q_re = INFINITY * (x_re * y_re + x_im * y_im);
 164             q_im = INFINITY * (x_im * y_re - x_re * y_im);
 165         }
 166         // finite / infinite
 167         else if (isinf(logbw) && isfinite(x_re) && isfinite(x_im))
 168         {
 169             y_re = copysignl(isinf(y_re) ? 1.0 : 0.0, y_re);
 170             y_im = copysignl(isinf(y_im) ? 1.0 : 0.0, y_im);
 171             q_re = 0.0 * (x_re * y_re + x_im * y_im);
 172             q_im = 0.0 * (x_im * y_re - x_re * y_im);
 173         }
 174     }
 175     return q_re + q_im * 1.0i;
 176 +/
 177     asm
 178     {
 179         fld     q_re;
 180         fld     q_im;
 181     }
 182 }
 183
 184 /****************************
 185  * Compare two complex floating point numbers, x and y.
 186  * Input:
 187  *      x.re    ST3
 188  *      x.im    ST2
 189  *      y.re    ST1
 190  *      y.im    ST0
 191  * Output:
 192  *      8087 stack is cleared
 193  *      flags set
 194  */
 195
 196 void _Ccmp()
 197 {
 198     asm
 199     {   naked                   ;
 200         fucomp  ST(2)           ; // compare x.im and y.im
 201         fstsw   AX              ;
 202         sahf                    ;
 203         jne     L1              ;
 204         jp      L1              ; // jmp if NAN
 205         fucomp  ST(2)           ; // compare x.re and y.re
 206         fstsw   AX              ;
 207         sahf                    ;
 208         fstp    ST(0)           ; // pop
 209         fstp    ST(0)           ; // pop
 210         ret                     ;
 211
 212       L1:
 213         fstp    ST(0)           ; // pop
 214         fstp    ST(0)           ; // pop
 215         fstp    ST(0)           ; // pop
 216         ret                     ;
 217     }
 218 }