/*  Basic complex arithmetic.  Everything assumed to be in Cartesian form	and made to stay that way.*/#include <stdio.h>#include "bigfloat.h"/*  set a float to a constant 1  */void one( FLOAT *x){	null( x);	x->expnt = 1;	x->mntsa.e[MS_MNTSA] = 0x40000000;}/*  zero out a complex storage location  */void null_cmplx( COMPLEX *x){	null( & x->real);	null( & x->imag);}/*  add two complex numbers.	c = a + b */void add_cmplx( COMPLEX *a, COMPLEX *b, COMPLEX *c){	COMPLEX mya, myb;	copy_cmplx( a, &mya);	copy_cmplx( b, &myb);		add( &mya.real, &myb.real, &c->real);	add( &mya.imag, &myb.imag, &c->imag);}/*  subtract two complex numbers	c = a - b*/void subtract_cmplx( COMPLEX *a, COMPLEX *b, COMPLEX *c){	COMPLEX myb;	copy_cmplx( b, &myb);	negate( &myb.real);	negate( &myb.imag);	add_cmplx( a, &myb, c);}/*  multiply two complex numbers.	c = (a.real * b.real - a.imag * b.imag) + i(a.imag * b.real + a.real * b.imag)*/void multiply_cmplx(  COMPLEX *a, COMPLEX *b, COMPLEX *c){	COMPLEX mya, myb;	FLOAT	temp1, temp2;	copy_cmplx( a, &mya);	copy_cmplx( b, &myb);	multiply( &mya.real, &myb.real, &temp1);	multiply( &mya.imag, &myb.imag, &temp2);	subtract( &temp1, &temp2, &c->real);	multiply( &mya.real, &myb.imag, &temp1);	multiply( &mya.imag, &myb.real, &temp2);	add( &temp1, &temp2, &c->imag);}/*  divide two complex numbers.	To keep things in Cartesian form multiply top by	conjugate of bottom.  Scale result by magnitude of	bottom.			output is c = ( a * b^*) / |b|		returns 1 if b != 0, 0 if |b| = 0*/int divide_cmplx( COMPLEX *a, COMPLEX *b, COMPLEX *c){	FLOAT	mag1, mag2;	COMPLEX	myb;		copy_cmplx( b, &myb);	multiply( &myb.real, &myb.real, &mag1);	multiply( &myb.imag, &myb.imag, &mag2);	add( &mag1, &mag2, &mag1);	negate( &myb.imag);	multiply_cmplx( a, &myb, c);	if( ! divide( &c->real, &mag1, &c->real)) return 0;	divide( &c->imag, &mag1, &c->imag);	return 1;}/*  compute y = x^k	where k is a signed integer in range +/-2^31 and x, y are complex.	returns 1 if ok, 0 if x = 0 and k < 0*/int intpwr_cmplx( COMPLEX *x, int k, COMPLEX  *y){	int signflag, n;	COMPLEX z, t;	/*	FLOAT seven, temp;		null(&seven);	seven.expnt = 3;	seven.mntsa.e[MS_MNTSA] = 0x70000000;	square_root( &seven, &seven);	/*  initialize Knuth's algorithm A pg 442 semi-numerical algorithms  */	copy_cmplx( x, &z);	if ( k < 0 )	{		signflag = 1;		n = -k;	}	else	{		signflag = 0;		n = k;	}	null_cmplx( &t);	one( &t.real);	while (n)	{		if ( n & 1 ) multiply_cmplx( &t, &z, &t);		multiply_cmplx( &z, &z, &z);		n >>= 1;	}	if ( signflag)	{		null_cmplx( &z);		one( &z.real);		return divide_cmplx( &z, &t, y);	}	copy_cmplx( &t, y);	return 1;}/*  compute magnitude of a complex number.  Returns	FLOAT result.*/void magnitude_cmplx( COMPLEX *x, FLOAT *m){	FLOAT x2, y2;		multiply( &x->real, &x->real, &x2);	multiply( &x->imag, &x->imag, &y2);	add( &x2, &y2, m);	square_root( m, m);}/*  complex exponential.  Change limit depending on how accurate	you want to be.  This is crude brute force based on definition	of exponential.  Can be done in place.  	Note that if it takes more than 2^31 terms to converge, this	WILL FAIL!!  So magnitude should be < 1 if you want this to work.*/void	exp_cmplx( COMPLEX *z, COMPLEX *xp){	COMPLEX	zn, factor, myz;	int	i, done;	FLOAT	limit, check;		copy_cmplx( z, &myz);	copy_cmplx( z, &zn);	null_cmplx( xp);	int_to_float( 1, &xp->real);	int_to_float( 1, &limit);	limit.expnt = - 64;		done = 0;	i = 1;	while( !done)	{/*  add in next term of z^n/n!  */		add_cmplx( &zn, xp, xp);/*  compute z^n+1  */		multiply_cmplx( &myz, &zn, &zn);/*  divide by n+1 to give (n+1)! on bottom  */		i++;		int_to_float(i, &factor.real);		null( &factor.imag);		divide_cmplx( &zn, &factor, &zn);/*  check if this next term is smaller than limit  */		magnitude_cmplx( &zn, &check);		if( compare( &check, &limit) < 0) done = 1;		printfloat("check =", &check);	}	add_cmplx( &zn, xp, xp);}