|
|
1.1 root 1: /*
2: * @(#)z_div.c 2.1 6/7/85
3: */
4:
5: #include "complex"
6: #include <stdio.h>
7: #include <errno.h>
8:
9: /* THIS IS BASED ON TAHOE FP REPRESENTATION! */
10: #include <FP.h>
11:
12: z_div(c, a, b)
13: dcomplex *a, *b, *c;
14: {
15: double ratio, den;
16: double abr, abi;
17:
18: if( (abr = b->dreal) < 0.)
19: *((long int *)&abr) ^= SIGN_BIT;
20: if( (abi = b->dimag) < 0.)
21: *((long int *)&abi) ^= SIGN_BIT;
22: if( abr <= abi )
23: {
24: if(abi == 0) {
25: fprintf( stderr, "Double complex division by zero\n" );
26: /* f77_abort(EDOM); */
27: }
28: ratio = b->dreal / b->dimag ;
29: den = b->dimag * (1 + ratio*ratio);
30: c->dreal = (a->dreal*ratio + a->dimag) / den;
31: c->dimag = (a->dimag*ratio - a->dreal) / den;
32: }
33:
34: else
35: {
36: ratio = b->dimag / b->dreal ;
37: den = b->dreal * (1 + ratio*ratio);
38: c->dreal = (a->dreal + a->dimag*ratio) / den;
39: c->dimag = (a->dimag - a->dreal*ratio) / den;
40: }
41:
42: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.