Complex Arithmetic (complex.c) 949 c.1=-z.1; return c; fcomplex Cdiv(fcomplex a,fcomplex b) fcomplex c; float r,den; if (fabs(b.r)>=fabs(b.i)){ r=b.i/b.r; den=b.r+r*b.i; c.r=(a.r+r*a.i)/den; Permission is Sample page c.i=(a.i-r*a.r)/den; else r=b.r/b.1 den=b.i+r*b.r; c.r=(a.r*r+a.i)/den; c.i=(a.i*r-a.r)/den; return c; float Cabs(fcomplex z) http://ww.nr.com or call 1-800-872-7423(North America float x,y,ans,temp; x=fabs(z.r); y=fabs(z.i); 1f(x=0.0) ans=yi e1se1f(y==0.0) readable files (including this one)to any server computer,is strictly prohibited. ans=x; else if (x>y){ temp=y/x; ans=x*sgrt (1.0+temp*temp); else temp=x/yi ans=y*sqrt (1.0+temp*temp); return ans; Copyright (C)1988-1992 by Cambridge University Press.Programs Copyright(C)1988-1992 by Numerical Recipes from NUMERICAL RECIPES IN C:THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) fcomplex Csqrt(fcomplex z) fcomplex c; float x,y,W,ri 1f(z.r=0.0)&(z.1=0.0){ personal use.Further reproduction,or c.r=0.0 c.1=0.0; return c; else only),orsend email to directcustserv@cambridge.org(outside North America). x=fabs(z.r): Software. y=fabs(z.i）; 1f(x>=y) r-y/x; w=sqrt(x)*gqrt(0.5*(1.0+sqrt(1.0+r*r))); else r=x/y; w=sqrt(y)*sqrt(0.5*(r+sgrt(1.0+r*r))); 1f(z.r>0.0)[ c.r=w; c.1=z.1/(2.0*w); else c.1=(z.i>=0)？:-wComplex Arithmetic (complex.c) 949 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). c.i = -z.i; return c; } fcomplex Cdiv(fcomplex a, fcomplex b) { fcomplex c; float r,den; if (fabs(b.r) >= fabs(b.i)) { r=b.i/b.r; den=b.r+r*b.i; c.r=(a.r+r*a.i)/den; c.i=(a.i-r*a.r)/den; } else { r=b.r/b.i; den=b.i+r*b.r; c.r=(a.r*r+a.i)/den; c.i=(a.i*r-a.r)/den; } return c; } float Cabs(fcomplex z) { float x,y,ans,temp; x=fabs(z.r); y=fabs(z.i); if (x == 0.0) ans=y; else if (y == 0.0) ans=x; else if (x > y) { temp=y/x; ans=x*sqrt(1.0+temp*temp); } else { temp=x/y; ans=y*sqrt(1.0+temp*temp); } return ans; } fcomplex Csqrt(fcomplex z) { fcomplex c; float x,y,w,r; if ((z.r == 0.0) && (z.i == 0.0)) { c.r=0.0; c.i=0.0; return c; } else { x=fabs(z.r); y=fabs(z.i); if (x >= y) { r=y/x; w=sqrt(x)*sqrt(0.5*(1.0+sqrt(1.0+r*r))); } else { r=x/y; w=sqrt(y)*sqrt(0.5*(r+sqrt(1.0+r*r))); } if (z.r >= 0.0) { c.r=w; c.i=z.i/(2.0*w); } else { c.i=(z.i >= 0) ? w : -w;