110550Sdlw /* 222997Skre * Copyright (c) 1980 Regents of the University of California. 322997Skre * All rights reserved. The Berkeley software License Agreement 422997Skre * specifies the terms and conditions for redistribution. 522997Skre * 6*45967Sbostic * @(#)z_div.c 5.3 01/15/91 710550Sdlw */ 810550Sdlw 910550Sdlw #include "complex" 1020187Slibs #include <stdio.h> 1120187Slibs #include <errno.h> 1229971Smckusick #ifdef tahoe 13*45967Sbostic #include <tahoe/math/FP.h> 14*45967Sbostic #endif 1510550Sdlw 1610550Sdlw z_div(c, a, b) 1710550Sdlw dcomplex *a, *b, *c; 1810550Sdlw { 1910550Sdlw double ratio, den; 2010550Sdlw double abr, abi; 2110550Sdlw 2229971Smckusick #ifndef tahoe 2310550Sdlw if( (abr = b->dreal) < 0.) 2410550Sdlw abr = - abr; 2510550Sdlw if( (abi = b->dimag) < 0.) 2610550Sdlw abi = - abi; 2729971Smckusick #else tahoe 2829971Smckusick if( (abr = b->dreal) < 0.) 2929971Smckusick *((long int *)&abr) ^= SIGN_BIT; 3029971Smckusick if( (abi = b->dimag) < 0.) 3129971Smckusick *((long int *)&abi) ^= SIGN_BIT; 3229971Smckusick #endif tahoe 3310550Sdlw if( abr <= abi ) 3410550Sdlw { 3520187Slibs if(abi == 0) { 3620187Slibs fprintf( stderr, "Double complex division by zero\n" ); 3720187Slibs f77_abort(EDOM); 3820187Slibs } 3910550Sdlw ratio = b->dreal / b->dimag ; 4010550Sdlw den = b->dimag * (1 + ratio*ratio); 4110550Sdlw c->dreal = (a->dreal*ratio + a->dimag) / den; 4210550Sdlw c->dimag = (a->dimag*ratio - a->dreal) / den; 4310550Sdlw } 4410550Sdlw 4510550Sdlw else 4610550Sdlw { 4710550Sdlw ratio = b->dimag / b->dreal ; 4810550Sdlw den = b->dreal * (1 + ratio*ratio); 4910550Sdlw c->dreal = (a->dreal + a->dimag*ratio) / den; 5010550Sdlw c->dimag = (a->dimag - a->dreal*ratio) / den; 5110550Sdlw } 5210550Sdlw 5310550Sdlw } 54