phobos/internal/cmath2.d

// Copyright (C) 2001-2003 by Digital Mars
// All Rights Reserved
// www.digitalmars.com

// Runtime support for complex arithmetic code generation
// (for linux)

import std.math;

extern (C):

/****************************
 * Multiply two complex floating point numbers, x and y.
 * Input:
 *	x.re	ST3
 *	x.im	ST2
 *	y.re	ST1
 *	y.im	ST0
 * Output:
 *	ST1	real part
 *	ST0	imaginary part
 */

void _Cmul()
{
    // p.re = x.re * y.re - x.im * y.im;
    // p.im = x.im * y.re + x.re * y.im;
    asm
    {	naked			;
	fld	ST(1)		; // x.re
	fmul	ST,ST(4)	; // ST0 = x.re * y.re

	fld	ST(1)		; // y.im
	fmul	ST,ST(4)	; // ST0 = x.im * y.im

	fsubp	ST(1),ST	; // ST0 = x.re * y.re - x.im * y.im

	fld	ST(3)		; // x.im
	fmul	ST,ST(3)	; // ST0 = x.im * y.re

	fld	ST(5)		; // x.re
	fmul	ST,ST(3)	; // ST0 = x.re * y.im

	faddp	ST(1),ST	; // ST0 = x.im * y.re + x.re * y.im

	fxch	ST(4),ST	;
	fstp	ST(0)		;
	fxch	ST(4),ST	;
	fstp	ST(0)		;
	fstp	ST(0)		;
	fstp	ST(0)		;

	ret			;
    }
/+
    if (isnan(x) && isnan(y))
    {
	// Recover infinities that computed as NaN+ iNaN ...
	int recalc = 0;
	if ( isinf( a) || isinf( b) )
	{   // z is infinite
	    // "Box" the infinity and change NaNs in the other factor to 0
	    a = copysignl( isinf( a) ? 1.0 : 0.0, a);
	    b = copysignl( isinf( b) ? 1.0 : 0.0, b);
	    if (isnan( c)) c = copysignl( 0.0, c);
	    if (isnan( d)) d = copysignl( 0.0, d);
	    recalc = 1;
	}
	if (isinf(c) || isinf(d))
	{   // w is infinite
	    // "Box" the infinity and change NaNs in the other factor to 0
	    c = copysignl( isinf( c) ? 1.0 : 0.0, c);
	    d = copysignl( isinf( d) ? 1.0 : 0.0, d);
	    if (isnan( a)) a = copysignl( 0.0, a);
	    if (isnan( b)) b = copysignl( 0.0, b);
	    recalc = 1;
	}
	if (!recalc && (isinf(ac) || isinf(bd) ||
	    isinf(ad) || isinf(bc)))
	{
	    // Recover infinities from overflow by changing NaNs to 0 ...
	    if (isnan( a)) a = copysignl( 0.0, a);
	    if (isnan( b)) b = copysignl( 0.0, b);
	    if (isnan( c)) c = copysignl( 0.0, c);
	    if (isnan( d)) d = copysignl( 0.0, d);
	    recalc = 1;
	}
	if (recalc)
	{
	    x = INFINITY * (a * c - b * d);
	    y = INFINITY * (a * d + b * c);
	}
    }
+/
}

/****************************
 * Divide two complex floating point numbers, x / y.
 * Input:
 *	x.re	ST3
 *	x.im	ST2
 *	y.re	ST1
 *	y.im	ST0
 * Output:
 *	ST1	real part
 *	ST0	imaginary part
 */

void _Cdiv()
{
    real x_re, x_im;
    real y_re, y_im;
    real q_re, q_im;
    real r;
    real den;

    asm
    {
	fstp	y_im	;
	fstp	y_re	;
	fstp	x_im	;
	fstp	x_re	;
    }

    if (fabs(y_re) < fabs(y_im))
    {
	r = y_re / y_im;
	den = y_im + r * y_re;
	q_re = (x_re * r + x_im) / den;
	q_im = (x_im * r - x_re) / den;
    }
    else
    {
	r = y_im / y_re;
	den = y_re + r * y_im;
	q_re = (x_re + r * x_im) / den;
	q_im = (x_im - r * x_re) / den;
    }
//printf("q.re = %g, q.im = %g\n", (double)q_re, (double)q_im);
/+
    if (isnan(q_re) && isnan(q_im))
    {
	real denom = y_re * y_re + y_im * y_im;

	// non-zero / zero
	if (denom == 0.0 && (!isnan(x_re) || !isnan(x_im)))
	{
	    q_re = copysignl(INFINITY, y_re) * x_re;
	    q_im = copysignl(INFINITY, y_re) * x_im;
	}
	// infinite / finite
	else if ((isinf(x_re) || isinf(x_im)) && isfinite(y_re) && isfinite(y_im))
	{
	    x_re = copysignl(isinf(x_re) ? 1.0 : 0.0, x_re); 
	    x_im = copysignl(isinf(x_im) ? 1.0 : 0.0, x_im); 
	    q_re = INFINITY * (x_re * y_re + x_im * y_im); 
	    q_im = INFINITY * (x_im * y_re - x_re * y_im); 
	}
	// finite / infinite
	else if (isinf(logbw) && isfinite(x_re) && isfinite(x_im))
	{
	    y_re = copysignl(isinf(y_re) ? 1.0 : 0.0, y_re); 
	    y_im = copysignl(isinf(y_im) ? 1.0 : 0.0, y_im); 
	    q_re = 0.0 * (x_re * y_re + x_im * y_im); 
	    q_im = 0.0 * (x_im * y_re - x_re * y_im); 
	}
    }
    return q_re + q_im * 1.0i;
+/
    asm
    {
	fld	q_re;
	fld	q_im;
    }
}

/****************************
 * Compare two complex floating point numbers, x and y.
 * Input:
 *	x.re	ST3
 *	x.im	ST2
 *	y.re	ST1
 *	y.im	ST0
 * Output:
 *	8087 stack is cleared
 *	flags set
 */

void _Ccmp()
{
    asm
    {	naked			;
	fucomp	ST(2)		; // compare x.im and y.im
	fstsw	AX		;
	sahf			;
	jne	L1		;
	jp	L1		; // jmp if NAN
	fucomp	ST(2)		; // compare x.re and y.re
	fstsw	AX		;
	sahf			;
	fstp	ST(0)		; // pop
	fstp	ST(0)		; // pop
	ret			;

      L1:
	fstp	ST(0)		; // pop
	fstp	ST(0)		; // pop
	fstp	ST(0)		; // pop
	ret			;
    }
}