110550Sdlw /* 2*20187Slibs * "@(#)z_div.c 1.2" 310550Sdlw */ 410550Sdlw 510550Sdlw #include "complex" 6*20187Slibs #include <stdio.h> 7*20187Slibs #include <errno.h> 810550Sdlw 910550Sdlw z_div(c, a, b) 1010550Sdlw dcomplex *a, *b, *c; 1110550Sdlw { 1210550Sdlw double ratio, den; 1310550Sdlw double abr, abi; 1410550Sdlw 1510550Sdlw if( (abr = b->dreal) < 0.) 1610550Sdlw abr = - abr; 1710550Sdlw if( (abi = b->dimag) < 0.) 1810550Sdlw abi = - abi; 1910550Sdlw if( abr <= abi ) 2010550Sdlw { 21*20187Slibs if(abi == 0) { 22*20187Slibs fprintf( stderr, "Double complex division by zero\n" ); 23*20187Slibs f77_abort(EDOM); 24*20187Slibs } 2510550Sdlw ratio = b->dreal / b->dimag ; 2610550Sdlw den = b->dimag * (1 + ratio*ratio); 2710550Sdlw c->dreal = (a->dreal*ratio + a->dimag) / den; 2810550Sdlw c->dimag = (a->dimag*ratio - a->dreal) / den; 2910550Sdlw } 3010550Sdlw 3110550Sdlw else 3210550Sdlw { 3310550Sdlw ratio = b->dimag / b->dreal ; 3410550Sdlw den = b->dreal * (1 + ratio*ratio); 3510550Sdlw c->dreal = (a->dreal + a->dimag*ratio) / den; 3610550Sdlw c->dimag = (a->dimag - a->dreal*ratio) / den; 3710550Sdlw } 3810550Sdlw 3910550Sdlw } 40