TBCI Numerical high perf. C++ Library  2.8.0
TOMS_707.C
Go to the documentation of this file.
1 
5 /* toms_707.f -- translated by f2c (version 19980516).
6  You must link the resulting object file with the libraries:
7  -lf2c -lm (in that order)
8 */
9 
10 /* $Id: TOMS_707.C,v 1.2.2.11 2019/05/28 11:13:02 garloff Exp $ */
11 
12 #if defined(__GNUC__) && __GNUC__ == 2 && defined(__cplusplus)
13 # undef __cplusplus /* workaround for buggy gcc-2.96 */
14 #endif
15 
16 #include "tbci/specfun/prototypes.h"
17 #include "tbci/specfun/prototypes2.h"
18 
19 #ifdef __cplusplus
20 extern "C" {
21 #endif
22 #include "tbci/lapack/f2c.h"
23 
24 /* Common Block Declarations */
25 
26 struct {
27  doublereal y;
28 } stcom_;
29 
30 #define stcom_1 stcom_
31 
32 /* Table of constant values */
33 
34 static doublereal c_b8 = 2.;
35 static doublereal c_b53 = 1.;
36 static doublereal c_b65 = 10.;
37 static integer c__7 = 7;
38 static integer c__1 = 1;
39 static integer c__9 = 9;
40 static integer c__3 = 3;
41 static integer c__2 = 2;
42 
43 /* **************************************************************** */
44 int aradd_(const doublereal *a, const doublereal *b,
45  doublereal *c__, const integer *l, const doublereal *rmax);
46 int arsub_(const doublereal *a, const doublereal *b,
47  doublereal *c__, const integer *l, const doublereal *rmax);
48 
49 int armult_(const doublereal *a, const doublereal *b, doublereal *c__,
50  const integer *l, const doublereal *rmax);
51 int arydiv_(const doublereal *ar, const doublereal *ai,
52  const doublereal *br, const doublereal *bi,
53  doublecomplex *c__, const integer *l, const integer *lnchf,
54  const doublereal *rmax, const integer *bit);
55 
56 int cmpadd_(const doublereal *ar, const doublereal *ai,
57  const doublereal *br, const doublereal *bi,
58  doublereal *cr, doublereal *ci,
59  const integer *l, const doublereal *rmax);
60 int cmpmul_(const doublereal *ar, const doublereal *ai,
61  const doublereal *br, const doublereal *bi,
62  doublereal *cr, doublereal *ci,
63  const integer *l, const doublereal *rmax);
64 
65  extern integer bits_(void);
66  extern doublereal store_(doublereal *x);
67 
68  extern /* Subroutine */ int conv12_(doublecomplex *cn, doublereal *cae), conv21_(doublereal *cae, doublecomplex *cn);
69 
70  extern /* Subroutine */ int emult_(doublereal *n1, doublereal *e1, doublereal *n2, doublereal *e2, doublereal *nf, doublereal *ef);
71  extern /* Subroutine */ int ecpmul_(doublereal *a, doublereal *b, doublereal *c__), eadd_(doublereal *n1, doublereal *e1, doublereal *n2, doublereal *e2, doublereal *nf, doublereal *ef);
72  extern /* Subroutine */ int ecpdiv_(doublereal *a, doublereal *b, doublereal *c__);
73  extern /* Subroutine */ int ediv_(doublereal *n1, doublereal *e1, doublereal *n2, doublereal *e2, doublereal *nf, doublereal *ef);
74  extern /* Subroutine */ int eadd_(doublereal *n1, doublereal *e1, doublereal *n2, doublereal *e2, doublereal *nf, doublereal *ef);
75  extern /* Subroutine */ int esub_(doublereal *n1, doublereal *e1, doublereal *n2, doublereal *e2, doublereal *nf, doublereal *ef);
76 
77 
78 
79 /* Double Complex */
80 VOID chgf_(doublecomplex *ret_val, const doublecomplex *a, const doublecomplex *b, const doublecomplex *z__, const integer *l, const integer *lnchf);
81 
82 /* ALGORITHM 707, COLLECTED ALGORITHMS FROM ACM. */
83 /* THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, */
84 /* VOL. 18, NO. 3, SEPTEMBER, 1992, PP. 345-349. */
85 /* **************************************************************** */
86 /* * * */
87 /* * SOLUTION TO THE CONFLUENT HYPERGEOMETRIC FUNCTION * */
88 /* * * */
89 /* * by * */
90 /* * * */
91 /* * MARK NARDIN, * */
92 /* * * */
93 /* * W. F. PERGER and ATUL BHALLA * */
94 /* * * */
95 /* * * */
96 /* * Michigan Technological University, Copyright 1989 * */
97 /* * * */
98 /* * * */
99 /* * Description : A numerical evaluator for the confluent * */
100 /* * hypergeometric function for complex arguments with large * */
101 /* * magnitudes using a direct summation of the Kummer series. * */
102 /* * The method used allows an accuracy of up to thirteen * */
103 /* * decimal places through the use of large real arrays * */
104 /* * and a single final division. LNCHF is a variable which * */
105 /* * selects how the result should be represented. A '0' will * */
106 /* * return the value in standard exponential form. A '1' * */
107 /* * will return the LOG of the result. IP is an integer * */
108 /* * variable that specifies how many array positions are * */
109 /* * desired (usually 10 is sufficient). Setting IP=0 causes * */
110 /* * the program to estimate the number of array positions. * */
111 /* * * */
112 /* * The confluent hypergeometric function is the solution to * */
113 /* * the differential equation: * */
114 /* * * */
115 /* * zf"(z) + (a-z)f'(z) - bf(z) = 0 * */
116 /* * * */
117 /* * Subprograms called: BITS, CHGF * */
118 /* * * */
119 /* **************************************************************** */
120 /* Double Complex */ /*VOID*/
121 void conhyp_(doublecomplex *ret_val,
122  const doublecomplex *a, const doublecomplex *b,
123  const doublecomplex *z__, const integer *lnchf, const integer *ip)
124 {
125  /* System generated locals */
126  doublecomplex z__1, z__2, z__3, z__4, z__5, z__6, z__7;
127 
128  /* Builtin functions: prototypes.h */
129 
130  /* Local variables */
131  static doublereal term1, term2;
132  static integer i__;
133  static doublereal nterm, fx, ang, max__;
134 
135  if (z_abs(z__) != 0.) {
136  ang = atan2(d_imag(z__), (doublereal) z__->r);
137  } else {
138  ang = 1.;
139  }
140  if (abs(ang) < 1.570796) {
141  ang = 1.;
142  } else {
143  ang = sin(abs(ang) - 1.570796325) + 1.;
144  }
145  max__ = 0.;
146  nterm = 0.;
147  fx = 0.;
148  term1 = 0.;
149 L10:
150  nterm += 1;
151  z__4.r = a->r + nterm, z__4.i = a->i;
152  z__3.r = z__4.r - 1, z__3.i = z__4.i;
153  z__2.r = z__3.r * z__->r - z__3.i * z__->i, z__2.i = z__3.r * z__->i +
154  z__3.i * z__->r;
155  z__7.r = b->r + nterm, z__7.i = b->i;
156  z__6.r = z__7.r - 1, z__6.i = z__7.i;
157  z__5.r = nterm * z__6.r, z__5.i = nterm * z__6.i;
158  z_div(&z__1, &z__2, &z__5);
159  term2 = z_abs(&z__1);
160  if (term2 == 0.) {
161  goto L20;
162  }
163  if (term2 < 1.) {
164  if (a->r + nterm - 1 > 1.) {
165  if (b->r + nterm - 1 > 1.) {
166  if (term2 - term1 < 0.) {
167  goto L20;
168  }
169  }
170  }
171  }
172  fx += log(term2);
173  if (fx > max__) {
174  max__ = fx;
175  }
176  term1 = term2;
177  goto L10;
178 L20:
179  max__ = max__ * 2 / (bits_() * .69314718056);
180  i__ = (integer) (max__ * ang) + 7;
181  if (i__ < 5) {
182  i__ = 5;
183  }
184  if (*ip > i__) {
185  i__ = *ip;
186  }
187  chgf_(&z__1, a, b, z__, &i__, lnchf);
188  ret_val->r = z__1.r, ret_val->i = z__1.i;
189  return ;
190 } /* conhyp_ */
191 
192 /* **************************************************************** */
193 /* * * */
194 /* * FUNCTION BITS * */
195 /* * * */
196 /* * * */
197 /* * Description : Determines the number of significant figures * */
198 /* * of machine precision to arrive at the size of the array * */
199 /* * the numbers must must be stored in to get the accuracy * */
200 /* * of the solution. * */
201 /* * * */
202 /* * Subprogram called: STORE * */
203 /* * * */
204 /* **************************************************************** */
205 integer bits_(void)
206 {
207  /* System generated locals */
208  integer ret_val;
209  doublereal d__1;
210 
211  /* Local variables */
212  static integer count;
213  static doublereal bit, bit2;
214 
215  bit = (float)1.;
216  count = 0;
217 L10:
218  ++count;
219  d__1 = bit * (float)2.;
220  bit2 = store_(&d__1);
221  d__1 = bit2 + (float)1.;
222  bit = store_(&d__1);
223  if (bit - bit2 != (float)0.) {
224  goto L10;
225  }
226  ret_val = count;
227  return ret_val;
228 } /* bits_ */
229 
230 doublereal store_(doublereal *x)
231 {
232  /* System generated locals */
233  doublereal ret_val;
234 
235 
236 /* *********************************************************** */
237 
238 
239 /* This function forces its argument X to be stored in a */
240 /* memory location, thus providing a means of determining */
241 /* floating point number characteristics (such as the machine */
242 /* precision) when it is necessary to avoid computation in */
243 /* high precision registers. */
244 
245 /* On input: */
246 
247 /* X = Value to be stored. */
248 
249 /* X is not altered by this function. */
250 
251 /* On output: */
252 
253 /* STORE = Value of X after it has been stored and */
254 /* possibly truncated or rounded to the double */
255 /* precision word length. */
256 
257 /* Modules required by STORE: None */
258 
259 /* *********************************************************** */
260 
261  stcom_1.y = *x;
262  ret_val = stcom_1.y;
263  return ret_val;
264 } /* store_ */
265 
266 /* **************************************************************** */
267 /* * * */
268 /* * FUNCTION CHGF * */
269 /* * * */
270 /* * * */
271 /* * Description : Function that sums the Kummer series and * */
272 /* * returns the solution of the confluent hypergeometric * */
273 /* * function. * */
274 /* * * */
275 /* * Subprograms called: ARMULT, ARYDIV, BITS, CMPADD, CMPMUL * */
276 /* * * */
277 /* **************************************************************** */
278 /* Double Complex */ VOID chgf_(doublecomplex *ret_val, const doublecomplex *a, const doublecomplex *b, const doublecomplex *z__, const integer *l, const integer *lnchf)
279 {
280  /* System generated locals */
281  integer i__1;
282  doublereal d__1, d__2;
283  doublecomplex z__1, z__2, z__3, z__4, z__5, z__6;
284 
285  /* Builtin functions */
286 
287  /* Local variables */
288  static doublereal rmax, numi[779], sumi[779], numr[779], sumr[779];
289  static integer i__;
290  static doublecomplex final;
291  static doublereal ai, ci, ar, cr;
292  static doublereal xi, sigfig, xr, denomi[779], denomr[779], ai2;
293  static doublereal ci2;
294  static doublereal ar2, cr2, qi1[779], qi2[779], xi2, qr1[779], qr2[779],
295  mx1, mx2, xr2;
296  static integer bit;
297  static doublereal cnt;
298 
299  bit = bits_();
300  i__1 = bit / 2;
301  rmax = pow_di(&c_b8, &i__1);
302  i__1 = bit / 4;
303  sigfig = pow_di(&c_b8, &i__1);
304  ar2 = a->r * sigfig;
305  ar = d_int(&ar2);
306  d__1 = (ar2 - ar) * rmax;
307  ar2 = d_nint(&d__1);
308  ai2 = d_imag(a) * sigfig;
309  ai = d_int(&ai2);
310  d__1 = (ai2 - ai) * rmax;
311  ai2 = d_nint(&d__1);
312  cr2 = b->r * sigfig;
313  cr = d_int(&cr2);
314  d__1 = (cr2 - cr) * rmax;
315  cr2 = d_nint(&d__1);
316  ci2 = d_imag(b) * sigfig;
317  ci = d_int(&ci2);
318  d__1 = (ci2 - ci) * rmax;
319  ci2 = d_nint(&d__1);
320  xr2 = z__->r * sigfig;
321  xr = d_int(&xr2);
322  d__1 = (xr2 - xr) * rmax;
323  xr2 = d_nint(&d__1);
324  xi2 = d_imag(z__) * sigfig;
325  xi = d_int(&xi2);
326  d__1 = (xi2 - xi) * rmax;
327  xi2 = d_nint(&d__1);
328  sumr[0] = 1.;
329  sumi[0] = 1.;
330  numr[0] = 1.;
331  numi[0] = 1.;
332  denomr[0] = 1.;
333  denomi[0] = 1.;
334  i__1 = *l + 1;
335  for (i__ = 0; i__ <= i__1; ++i__) {
336  sumr[i__ + 1] = 0.;
337  sumi[i__ + 1] = 0.;
338  numr[i__ + 1] = 0.;
339  numi[i__ + 1] = 0.;
340  denomr[i__ + 1] = 0.;
341  denomi[i__ + 1] = 0.;
342 /* L100: */
343  }
344  sumr[2] = 1.;
345  numr[2] = 1.;
346  denomr[2] = 1.;
347  cnt = sigfig;
348 L110:
349  if (sumr[2] < (float).5) {
350  mx1 = sumi[*l + 2];
351  } else if (sumi[2] < (float).5) {
352  mx1 = sumr[*l + 2];
353  } else {
354 /* Computing MAX */
355  d__1 = sumr[*l + 2], d__2 = sumi[*l + 2];
356  mx1 = max(d__1,d__2);
357  }
358  if (numr[2] < (float).5) {
359  mx2 = numi[*l + 2];
360  } else if (numi[2] < (float).5) {
361  mx2 = numr[*l + 2];
362  } else {
363 /* Computing MAX */
364  d__1 = numr[*l + 2], d__2 = numi[*l + 2];
365  mx2 = max(d__1,d__2);
366  }
367  if (mx1 - mx2 > (float)2.) {
368  if (cr > 0.) {
369  z__3.r = ar, z__3.i = ai;
370  z__4.r = xr, z__4.i = xi;
371  z__2.r = z__3.r * z__4.r - z__3.i * z__4.i, z__2.i = z__3.r *
372  z__4.i + z__3.i * z__4.r;
373  z__6.r = cr, z__6.i = ci;
374  z__5.r = cnt * z__6.r, z__5.i = cnt * z__6.i;
375  z_div(&z__1, &z__2, &z__5);
376  if (z_abs(&z__1) <= 1.) {
377  goto L190;
378  }
379  }
380  }
381  cmpmul_(sumr, sumi, &cr, &ci, qr1, qi1, l, &rmax);
382  cmpmul_(sumr, sumi, &cr2, &ci2, qr2, qi2, l, &rmax);
383  qr2[*l + 2] += -1;
384  qi2[*l + 2] += -1;
385  cmpadd_(qr1, qi1, qr2, qi2, sumr, sumi, l, &rmax);
386  armult_(sumr, &cnt, sumr, l, &rmax);
387  armult_(sumi, &cnt, sumi, l, &rmax);
388  cmpmul_(denomr, denomi, &cr, &ci, qr1, qi1, l, &rmax);
389  cmpmul_(denomr, denomi, &cr2, &ci2, qr2, qi2, l, &rmax);
390  qr2[*l + 2] += -1;
391  qi2[*l + 2] += -1;
392  cmpadd_(qr1, qi1, qr2, qi2, denomr, denomi, l, &rmax);
393  armult_(denomr, &cnt, denomr, l, &rmax);
394  armult_(denomi, &cnt, denomi, l, &rmax);
395  cmpmul_(numr, numi, &ar, &ai, qr1, qi1, l, &rmax);
396  cmpmul_(numr, numi, &ar2, &ai2, qr2, qi2, l, &rmax);
397  qr2[*l + 2] += -1;
398  qi2[*l + 2] += -1;
399  cmpadd_(qr1, qi1, qr2, qi2, numr, numi, l, &rmax);
400  cmpmul_(numr, numi, &xr, &xi, qr1, qi1, l, &rmax);
401  cmpmul_(numr, numi, &xr2, &xi2, qr2, qi2, l, &rmax);
402  qr2[*l + 2] += -1;
403  qi2[*l + 2] += -1;
404  cmpadd_(qr1, qi1, qr2, qi2, numr, numi, l, &rmax);
405  cmpadd_(sumr, sumi, numr, numi, sumr, sumi, l, &rmax);
406  cnt += sigfig;
407  ar += sigfig;
408  cr += sigfig;
409  goto L110;
410 L190:
411  arydiv_(sumr, sumi, denomr, denomi, &final, l, lnchf, &rmax, &bit);
412  ret_val->r = final.r, ret_val->i = final.i;
413  return ;
414 } /* chgf_ */
415 
416 /* **************************************************************** */
417 /* * * */
418 /* * SUBROUTINE ARADD * */
419 /* * * */
420 /* * * */
421 /* * Description : Accepts two arrays of numbers and returns * */
422 /* * the sum of the array. Each array is holding the value * */
423 /* * of one number in the series. The parameter L is the * */
424 /* * size of the array representing the number and RMAX is * */
425 /* * the actual number of digits needed to give the numbers * */
426 /* * the desired accuracy. * */
427 /* * * */
428 /* * Subprograms called: none * */
429 /* * * */
430 /* **************************************************************** */
431 /* Subroutine */
432 int aradd_(const doublereal *a, const doublereal *b,
433  doublereal *c__, const integer *l, const doublereal *rmax)
434 {
435  /* System generated locals */
436  integer i__1;
437  doublereal d__1;
438 
439  /* Builtin functions */
440 
441  /* Local variables */
442  static integer i__, j, ediff;
443  static doublereal z__[779];
444 
445  /* Parameter adjustments */
446  ++c__;
447  ++b;
448  ++a;
449 
450  /* Function Body */
451  i__1 = *l + 1;
452  for (i__ = 0; i__ <= i__1; ++i__) {
453  z__[i__ + 1] = 0.;
454 /* L110: */
455  }
456  d__1 = a[*l + 1] - b[*l + 1];
457  ediff = (integer) d_nint(&d__1);
458  if (abs(a[1]) < (float).5 || ediff <= -(*l)) {
459  goto L111;
460  }
461  if (abs(b[1]) < (float).5 || ediff >= *l) {
462  goto L113;
463  }
464  goto L115;
465 L111:
466  i__1 = *l + 1;
467  for (i__ = -1; i__ <= i__1; ++i__) {
468  c__[i__] = b[i__];
469 /* L112: */
470  }
471  goto L311;
472 L113:
473  i__1 = *l + 1;
474  for (i__ = -1; i__ <= i__1; ++i__) {
475  c__[i__] = a[i__];
476 /* L114: */
477  }
478  goto L311;
479 L115:
480  z__[0] = a[-1];
481  if ((d__1 = a[-1] - b[-1], abs(d__1)) < (float).5) {
482  goto L200;
483  }
484  if (ediff > 0) {
485  z__[*l + 2] = a[*l + 1];
486  goto L233;
487  }
488  if (ediff < 0) {
489  z__[*l + 2] = b[*l + 1];
490  z__[0] = b[-1];
491  goto L266;
492  }
493  i__1 = *l;
494  for (i__ = 1; i__ <= i__1; ++i__) {
495  if (a[i__] > b[i__]) {
496  z__[*l + 2] = a[*l + 1];
497  goto L233;
498  }
499  if (a[i__] < b[i__]) {
500  z__[*l + 2] = b[*l + 1];
501  z__[0] = b[-1];
502  goto L266;
503  }
504 /* L120: */
505  }
506  goto L300;
507 L200:
508  if (ediff > 0) {
509  goto L203;
510  }
511  if (ediff < 0) {
512  goto L207;
513  }
514  z__[*l + 2] = a[*l + 1];
515  for (i__ = *l; i__ >= 1; --i__) {
516  z__[i__ + 1] = a[i__] + b[i__] + z__[i__ + 1];
517  if (z__[i__ + 1] >= *rmax) {
518  z__[i__ + 1] -= *rmax;
519  z__[i__] = 1.;
520  }
521 /* L201: */
522  }
523  if (z__[1] > (float).5) {
524  for (i__ = *l; i__ >= 1; --i__) {
525  z__[i__ + 1] = z__[i__];
526 /* L202: */
527  }
528  z__[*l + 2] += 1.;
529  z__[1] = 0.;
530  }
531  goto L300;
532 L203:
533  z__[*l + 2] = a[*l + 1];
534  i__1 = ediff + 1;
535  for (i__ = *l; i__ >= i__1; --i__) {
536  z__[i__ + 1] = a[i__] + b[i__ - ediff] + z__[i__ + 1];
537  if (z__[i__ + 1] >= *rmax) {
538  z__[i__ + 1] -= *rmax;
539  z__[i__] = 1.;
540  }
541 /* L204: */
542  }
543  for (i__ = ediff; i__ >= 1; --i__) {
544  z__[i__ + 1] = a[i__] + z__[i__ + 1];
545  if (z__[i__ + 1] >= *rmax) {
546  z__[i__ + 1] -= *rmax;
547  z__[i__] = 1.;
548  }
549 /* L205: */
550  }
551  if (z__[1] > (float).5) {
552  for (i__ = *l; i__ >= 1; --i__) {
553  z__[i__ + 1] = z__[i__];
554 /* L206: */
555  }
556  z__[*l + 2] += 1;
557  z__[1] = 0.;
558  }
559  goto L300;
560 L207:
561  z__[*l + 2] = b[*l + 1];
562  i__1 = 1 - ediff;
563  for (i__ = *l; i__ >= i__1; --i__) {
564  z__[i__ + 1] = a[i__ + ediff] + b[i__] + z__[i__ + 1];
565  if (z__[i__ + 1] >= *rmax) {
566  z__[i__ + 1] -= *rmax;
567  z__[i__] = 1.;
568  }
569 /* L208: */
570  }
571  for (i__ = -ediff; i__ >= 1; --i__) {
572  z__[i__ + 1] = b[i__] + z__[i__ + 1];
573  if (z__[i__ + 1] >= *rmax) {
574  z__[i__ + 1] -= *rmax;
575  z__[i__] = 1.;
576  }
577 /* L209: */
578  }
579  if (z__[1] > (float).5) {
580  for (i__ = *l; i__ >= 1; --i__) {
581  z__[i__ + 1] = z__[i__];
582 /* L210: */
583  }
584  z__[*l + 2] += 1.;
585  z__[1] = 0.;
586  }
587  goto L300;
588 L233:
589  if (ediff > 0) {
590  goto L243;
591  }
592  for (i__ = *l; i__ >= 1; --i__) {
593  z__[i__ + 1] = a[i__] - b[i__] + z__[i__ + 1];
594  if (z__[i__ + 1] < 0.) {
595  z__[i__ + 1] += *rmax;
596  z__[i__] = -1.;
597  }
598 /* L234: */
599  }
600  goto L290;
601 L243:
602  i__1 = ediff + 1;
603  for (i__ = *l; i__ >= i__1; --i__) {
604  z__[i__ + 1] = a[i__] - b[i__ - ediff] + z__[i__ + 1];
605  if (z__[i__ + 1] < 0.) {
606  z__[i__ + 1] += *rmax;
607  z__[i__] = -1.;
608  }
609 /* L244: */
610  }
611  for (i__ = ediff; i__ >= 1; --i__) {
612  z__[i__ + 1] = a[i__] + z__[i__ + 1];
613  if (z__[i__ + 1] < 0.) {
614  z__[i__ + 1] += *rmax;
615  z__[i__] = -1.;
616  }
617 /* L245: */
618  }
619  goto L290;
620 L266:
621  if (ediff < 0) {
622  goto L276;
623  }
624  for (i__ = *l; i__ >= 1; --i__) {
625  z__[i__ + 1] = b[i__] - a[i__] + z__[i__ + 1];
626  if (z__[i__ + 1] < 0.) {
627  z__[i__ + 1] += *rmax;
628  z__[i__] = -1.;
629  }
630 /* L267: */
631  }
632  goto L290;
633 L276:
634  i__1 = 1 - ediff;
635  for (i__ = *l; i__ >= i__1; --i__) {
636  z__[i__ + 1] = b[i__] - a[i__ + ediff] + z__[i__ + 1];
637  if (z__[i__ + 1] < 0.) {
638  z__[i__ + 1] += *rmax;
639  z__[i__] = -1.;
640  }
641 /* L277: */
642  }
643  for (i__ = -ediff; i__ >= 1; --i__) {
644  z__[i__ + 1] = b[i__] + z__[i__ + 1];
645  if (z__[i__ + 1] < 0.) {
646  z__[i__ + 1] += *rmax;
647  z__[i__] = -1.;
648  }
649 /* L278: */
650  }
651 L290:
652  if (z__[2] > (float).5) {
653  goto L300;
654  }
655  i__ = 1;
656 L291:
657  ++i__;
658  if (z__[i__ + 1] < (float).5 && i__ < *l + 1) {
659  goto L291;
660  }
661  if (i__ == *l + 1) {
662  z__[0] = 1.;
663  z__[*l + 2] = 0.;
664  goto L300;
665  }
666 /* L292: */
667  i__1 = *l + 1 - i__;
668  for (j = 1; j <= i__1; ++j) {
669  z__[j + 1] = z__[j + i__];
670 /* L293: */
671  }
672  i__1 = *l;
673  for (j = *l + 2 - i__; j <= i__1; ++j) {
674  z__[j + 1] = 0.;
675 /* L294: */
676  }
677  z__[*l + 2] = z__[*l + 2] - i__ + 1;
678 L300:
679  i__1 = *l + 1;
680  for (i__ = -1; i__ <= i__1; ++i__) {
681  c__[i__] = z__[i__ + 1];
682 /* L310: */
683  }
684 L311:
685  if (c__[1] < (float).5) {
686  c__[-1] = 1.;
687  c__[*l + 1] = 0.;
688  }
689  return 0;
690 } /* aradd_ */
691 
692 /* **************************************************************** */
693 /* * * */
694 /* * SUBROUTINE ARSUB * */
695 /* * * */
696 /* * * */
697 /* * Description : Accepts two arrays and subtracts each element * */
698 /* * in the second array from the element in the first array * */
699 /* * and returns the solution. The parameters L and RMAX are * */
700 /* * the size of the array and the number of digits needed for * */
701 /* * the accuracy, respectively. * */
702 /* * * */
703 /* * Subprograms called: ARADD * */
704 /* * * */
705 /* **************************************************************** */
706 /* Subroutine */
707 int arsub_(const doublereal *a, const doublereal *b,
708  doublereal *c__, const integer *l, const doublereal *rmax)
709 {
710  /* System generated locals */
711  integer i__1;
712 
713  /* Local variables */
714  static integer i__;
715  static doublereal b2[779];
716 
717  /* Parameter adjustments */
718  ++c__;
719  ++b;
720  ++a;
721 
722  /* Function Body */
723  i__1 = *l + 1;
724  for (i__ = -1; i__ <= i__1; ++i__) {
725  b2[i__ + 1] = b[i__];
726 /* L100: */
727  }
728  b2[0] *= -1.;
729  aradd_(&a[-1], b2, &c__[-1], l, rmax);
730  return 0;
731 } /* arsub_ */
732 
733 /* **************************************************************** */
734 /* * * */
735 /* * SUBROUTINE ARMULT * */
736 /* * * */
737 /* * * */
738 /* * Description : Accepts two arrays and returns the product. * */
739 /* * L and RMAX are the size of the arrays and the number of * */
740 /* * digits needed to represent the numbers with the required * */
741 /* * accuracy. * */
742 /* * * */
743 /* * Subprograms called: none * */
744 /* * * */
745 /* **************************************************************** */
746 /* Subroutine */ int armult_(const doublereal *a, const doublereal *b, doublereal *c__, const integer *l, const doublereal *rmax)
747 {
748  /* System generated locals */
749  integer i__1;
750  doublereal d__1;
751 
752  /* Builtin functions */
753 
754  /* Local variables */
755  static doublereal rmax2;
756  static integer i__;
757  static doublereal z__[779], carry, b2;
758 
759  /* Parameter adjustments */
760  ++c__;
761  ++a;
762 
763  /* Function Body */
764  rmax2 = 1. / *rmax;
765  z__[0] = d_sign(&c_b53, b) * a[-1];
766  b2 = abs(*b);
767  z__[*l + 2] = a[*l + 1];
768  i__1 = *l;
769  for (i__ = 0; i__ <= i__1; ++i__) {
770  z__[i__ + 1] = 0.;
771 /* L100: */
772  }
773  if (b2 <= 1e-10 || a[1] <= 1e-10) {
774  z__[0] = 1.;
775  z__[*l + 2] = 0.;
776  goto L198;
777  }
778  for (i__ = *l; i__ >= 1; --i__) {
779  z__[i__ + 1] = a[i__] * b2 + z__[i__ + 1];
780  if (z__[i__ + 1] >= *rmax) {
781  d__1 = z__[i__ + 1] / *rmax;
782  carry = d_int(&d__1);
783  z__[i__ + 1] -= carry * *rmax;
784  z__[i__] = carry;
785  }
786 /* L110: */
787  }
788  if (z__[1] < (float).5) {
789  goto L150;
790  }
791  for (i__ = *l; i__ >= 1; --i__) {
792  z__[i__ + 1] = z__[i__];
793 /* L120: */
794  }
795  z__[*l + 2] += 1.;
796  z__[1] = 0.;
797 L150:
798 L198:
799  i__1 = *l + 1;
800  for (i__ = -1; i__ <= i__1; ++i__) {
801  c__[i__] = z__[i__ + 1];
802 /* L199: */
803  }
804  if (c__[1] < (float).5) {
805  c__[-1] = 1.;
806  c__[*l + 1] = 0.;
807  }
808  return 0;
809 } /* armult_ */
810 
811 /* **************************************************************** */
812 /* * * */
813 /* * SUBROUTINE CMPADD * */
814 /* * * */
815 /* * * */
816 /* * Description : Takes two arrays representing one real and * */
817 /* * one imaginary part, and adds two arrays representing * */
818 /* * another complex number and returns two array holding the * */
819 /* * complex sum. * */
820 /* * (CR,CI) = (AR+BR, AI+BI) * */
821 /* * * */
822 /* * Subprograms called: ARADD * */
823 /* * * */
824 /* **************************************************************** */
825 /* Subroutine */ int cmpadd_(const doublereal *ar, const doublereal *ai, const doublereal *br, const doublereal *bi, doublereal *cr, doublereal *ci, const integer *l, const doublereal *rmax)
826 {
827 
828  /* Parameter adjustments */
829  ++ci;
830  ++cr;
831  ++bi;
832  ++br;
833  ++ai;
834  ++ar;
835 
836  /* Function Body */
837  aradd_(&ar[-1], &br[-1], &cr[-1], l, rmax);
838  aradd_(&ai[-1], &bi[-1], &ci[-1], l, rmax);
839  return 0;
840 } /* cmpadd_ */
841 
842 /* **************************************************************** */
843 /* * * */
844 /* * SUBROUTINE CMPSUB * */
845 /* * * */
846 /* * * */
847 /* * Description : Takes two arrays representing one real and * */
848 /* * one imaginary part, and subtracts two arrays representing * */
849 /* * another complex number and returns two array holding the * */
850 /* * complex sum. * */
851 /* * (CR,CI) = (AR+BR, AI+BI) * */
852 /* * * */
853 /* * Subprograms called: ARADD * */
854 /* * * */
855 /* **************************************************************** */
856 /* Subroutine */ int cmpsub_(doublereal *ar, doublereal *ai, doublereal *br, doublereal *bi, doublereal *cr, doublereal *ci, integer *l, doublereal *rmax)
857 {
858  /* Parameter adjustments */
859  ++ci;
860  ++cr;
861  ++bi;
862  ++br;
863  ++ai;
864  ++ar;
865 
866  /* Function Body */
867  arsub_(&ar[-1], &br[-1], &cr[-1], l, rmax);
868  arsub_(&ai[-1], &bi[-1], &ci[-1], l, rmax);
869  return 0;
870 } /* cmpsub_ */
871 
872 /* **************************************************************** */
873 /* * * */
874 /* * SUBROUTINE CMPMUL * */
875 /* * * */
876 /* * * */
877 /* * Description : Takes two arrays representing one real and * */
878 /* * one imaginary part, and multiplies it with two arrays * */
879 /* * representing another complex number and returns the * */
880 /* * complex product. * */
881 /* * * */
882 /* * Subprograms called: ARMULT, ARSUB, ARADD * */
883 /* * * */
884 /* **************************************************************** */
885 /* Subroutine */ int cmpmul_(const doublereal *ar, const doublereal *ai, const doublereal *br, const doublereal *bi, doublereal *cr, doublereal *ci, const integer *l, const doublereal *rmax)
886 {
887  static doublereal d1[779], d2[779];
888 
889  /* Parameter adjustments */
890  ++ci;
891  ++cr;
892  ++ai;
893  ++ar;
894 
895  /* Function Body */
896  armult_(&ar[-1], br, d1, l, rmax);
897  armult_(&ai[-1], bi, d2, l, rmax);
898  arsub_(d1, d2, &cr[-1], l, rmax);
899  armult_(&ar[-1], bi, d1, l, rmax);
900  armult_(&ai[-1], br, d2, l, rmax);
901  aradd_(d1, d2, &ci[-1], l, rmax);
902  return 0;
903 } /* cmpmul_ */
904 
905 /* **************************************************************** */
906 /* * * */
907 /* * SUBROUTINE ARYDIV * */
908 /* * * */
909 /* * * */
910 /* * Description : Returns the double precision complex number * */
911 /* * resulting from the division of four arrays, representing * */
912 /* * two complex numbers. The number returned will be in one * */
913 /* * two different forms. Either standard scientific or as * */
914 /* * the log of the number. * */
915 /* * * */
916 /* * Subprograms called: CONV21, CONV12, EADD, ECPDIV, EMULT * */
917 /* * * */
918 /* **************************************************************** */
919 /* Subroutine */ int arydiv_(const doublereal *ar, const doublereal *ai, const doublereal *br, const doublereal *bi, doublecomplex *c__, const integer *l, const integer *lnchf, const doublereal *rmax, const integer *bit)
920 {
921  /* System generated locals */
922  doublereal d__1;
923  doublecomplex z__1;
924 
925  /* Builtin functions */
926 
927  /* Local variables */
928  static integer rexp;
929  static doublereal x;
930  static doublereal e1, e2, e3;
931  static doublereal n1, n2, n3, x1, x2, ae[4] /* was [2][2] */, be[4] /*
932  was [2][2] */, ce[4] /* was [2][2] */;
933  static integer ii10, ir10;
934  static doublereal ri10, phi, rr10, dum1, dum2;
935  /* Parameter adjustments */
936  ++bi;
937  ++br;
938  ++ai;
939  ++ar;
940 
941  /* Function Body */
942  rexp = *bit / 2;
943  x = rexp * (ar[*l + 1] - 2);
944  rr10 = x * d_lg10(&c_b8) / d_lg10(&c_b65);
945  ir10 = (integer) rr10;
946  rr10 -= ir10;
947  x = rexp * (ai[*l + 1] - 2);
948  ri10 = x * d_lg10(&c_b8) / d_lg10(&c_b65);
949  ii10 = (integer) ri10;
950  ri10 -= ii10;
951  d__1 = ar[1] * *rmax * *rmax + ar[2] * *rmax + ar[3];
952  dum1 = d_sign(&d__1, &ar[-1]);
953  d__1 = ai[1] * *rmax * *rmax + ai[2] * *rmax + ai[3];
954  dum2 = d_sign(&d__1, &ai[-1]);
955  dum1 *= pow_dd(&c_b65, &rr10);
956  dum2 *= pow_dd(&c_b65, &ri10);
957  z__1.r = dum1, z__1.i = dum2;
958  conv12_(&z__1, ae);
959  ae[2] += ir10;
960  ae[3] += ii10;
961  x = rexp * (br[*l + 1] - 2);
962  rr10 = x * d_lg10(&c_b8) / d_lg10(&c_b65);
963  ir10 = (integer) rr10;
964  rr10 -= ir10;
965  x = rexp * (bi[*l + 1] - 2);
966  ri10 = x * d_lg10(&c_b8) / d_lg10(&c_b65);
967  ii10 = (integer) ri10;
968  ri10 -= ii10;
969  d__1 = br[1] * *rmax * *rmax + br[2] * *rmax + br[3];
970  dum1 = d_sign(&d__1, &br[-1]);
971  d__1 = bi[1] * *rmax * *rmax + bi[2] * *rmax + bi[3];
972  dum2 = d_sign(&d__1, &bi[-1]);
973  dum1 *= pow_dd(&c_b65, &rr10);
974  dum2 *= pow_dd(&c_b65, &ri10);
975  z__1.r = dum1, z__1.i = dum2;
976  conv12_(&z__1, be);
977  be[2] += ir10;
978  be[3] += ii10;
979  ecpdiv_(ae, be, ce);
980  if (*lnchf == 0) {
981  conv21_(ce, c__);
982  } else {
983  emult_(ce, &ce[2], ce, &ce[2], &n1, &e1);
984  emult_(&ce[1], &ce[3], &ce[1], &ce[3], &n2, &e2);
985  eadd_(&n1, &e1, &n2, &e2, &n3, &e3);
986  n1 = ce[0];
987  e1 = ce[2] - ce[3];
988  x2 = ce[1];
989  if (e1 > 74.) {
990  x1 = 1e75;
991  } else if (e1 < -74.) {
992  x1 = 0.;
993  } else {
994  x1 = n1 * pow_dd(&c_b65, &e1);
995  }
996  phi = atan2(x2, x1);
997  d__1 = (log(n3) + e3 * log(10.)) * .5;
998  z__1.r = d__1, z__1.i = phi;
999  c__->r = z__1.r, c__->i = z__1.i;
1000  }
1001  return 0;
1002 } /* arydiv_ */
1003 
1004 /* **************************************************************** */
1005 /* * * */
1006 /* * SUBROUTINE EMULT * */
1007 /* * * */
1008 /* * * */
1009 /* * Description : Takes one base and exponent and multiplies it * */
1010 /* * by another numbers base and exponent to give the product * */
1011 /* * in the form of base and exponent. * */
1012 /* * * */
1013 /* * Subprograms called: none * */
1014 /* * * */
1015 /* **************************************************************** */
1016 /* Subroutine */ int emult_(doublereal *n1, doublereal *e1, doublereal *n2, doublereal *e2, doublereal *nf, doublereal *ef)
1017 {
1018  *nf = *n1 * *n2;
1019  *ef = *e1 + *e2;
1020  if (abs(*nf) >= 10.) {
1021  *nf /= 10.;
1022  *ef += 1.;
1023  }
1024  return 0;
1025 } /* emult_ */
1026 
1027 /* **************************************************************** */
1028 /* * * */
1029 /* * SUBROUTINE EDIV * */
1030 /* * * */
1031 /* * * */
1032 /* * Description : returns the solution in the form of base and * */
1033 /* * exponent of the division of two exponential numbers. * */
1034 /* * * */
1035 /* * Subprograms called: none * */
1036 /* * * */
1037 /* **************************************************************** */
1038 /* Subroutine */ int ediv_(doublereal *n1, doublereal *e1, doublereal *n2, doublereal *e2, doublereal *nf, doublereal *ef)
1039 {
1040  *nf = *n1 / *n2;
1041  *ef = *e1 - *e2;
1042  if (abs(*nf) < 1. && *nf != 0.) {
1043  *nf *= 10.;
1044  *ef += -1.;
1045  }
1046  return 0;
1047 } /* ediv_ */
1048 
1049 /* **************************************************************** */
1050 /* * * */
1051 /* * SUBROUTINE EADD * */
1052 /* * * */
1053 /* * * */
1054 /* * Description : Returns the sum of two numbers in the form * */
1055 /* * of a base and an exponent. * */
1056 /* * * */
1057 /* * Subprograms called: none * */
1058 /* * * */
1059 /* **************************************************************** */
1060 /* Subroutine */ int eadd_(doublereal *n1, doublereal *e1, doublereal *n2, doublereal *e2, doublereal *nf, doublereal *ef)
1061 {
1062  /* Builtin functions */
1063 
1064  /* Local variables */
1065  static doublereal ediff;
1066 
1067  ediff = *e1 - *e2;
1068  if (ediff > 36.) {
1069  *nf = *n1;
1070  *ef = *e1;
1071  } else if (ediff < -36.) {
1072  *nf = *n2;
1073  *ef = *e2;
1074  } else {
1075  *nf = *n1 * pow_dd(&c_b65, &ediff) + *n2;
1076  *ef = *e2;
1077 L400:
1078  if (abs(*nf) < 10.) {
1079  goto L410;
1080  }
1081  *nf /= 10.;
1082  *ef += 1.;
1083  goto L400;
1084 L410:
1085  if (abs(*nf) >= 1. || *nf == 0.) {
1086  goto L420;
1087  }
1088  *nf *= 10.;
1089  *ef += -1.;
1090  goto L410;
1091  }
1092 L420:
1093  return 0;
1094 } /* eadd_ */
1095 
1096 /* **************************************************************** */
1097 /* * * */
1098 /* * SUBROUTINE ESUB * */
1099 /* * * */
1100 /* * * */
1101 /* * Description : Returns the solution to the subtraction of * */
1102 /* * two numbers in the form of base and exponent. * */
1103 /* * * */
1104 /* * Subprograms called: EADD * */
1105 /* * * */
1106 /* **************************************************************** */
1107 /* Subroutine */ int esub_(doublereal *n1, doublereal *e1, doublereal *n2, doublereal *e2, doublereal *nf, doublereal *ef)
1108 {
1109  /* System generated locals */
1110  doublereal d__1;
1111 
1112  /* Local variables */
1113 
1114  d__1 = *n2 * -1.;
1115  eadd_(n1, e1, &d__1, e2, nf, ef);
1116  return 0;
1117 } /* esub_ */
1118 
1119 /* **************************************************************** */
1120 /* * * */
1121 /* * SUBROUTINE CONV12 * */
1122 /* * * */
1123 /* * * */
1124 /* * Description : Converts a number from complex notation to a * */
1125 /* * form of a 2x2 real array. * */
1126 /* * * */
1127 /* * Subprograms called: none * */
1128 /* * * */
1129 /* **************************************************************** */
1130 /* Subroutine */ int conv12_(doublecomplex *cn, doublereal *cae)
1131 {
1132  /* Builtin functions */
1133 
1134  /* Parameter adjustments */
1135  cae -= 3;
1136 
1137  /* Function Body */
1138  cae[3] = cn->r;
1139  cae[5] = 0.;
1140 L300:
1141  if (abs(cae[3]) < 10.) {
1142  goto L310;
1143  }
1144  cae[3] /= 10.;
1145  cae[5] += 1.;
1146  goto L300;
1147 L310:
1148  if (abs(cae[3]) >= 1. || cae[3] == 0.) {
1149  goto L320;
1150  }
1151  cae[3] *= 10.;
1152  cae[5] += -1.;
1153  goto L310;
1154 L320:
1155  cae[4] = d_imag(cn);
1156  cae[6] = 0.;
1157 L330:
1158  if (abs(cae[4]) < 10.) {
1159  goto L340;
1160  }
1161  cae[4] /= 10.;
1162  cae[6] += 1.;
1163  goto L330;
1164 L340:
1165  if (abs(cae[4]) >= 1. || cae[4] == 0.) {
1166  goto L350;
1167  }
1168  cae[4] *= 10.;
1169  cae[6] += -1.;
1170  goto L340;
1171 L350:
1172  return 0;
1173 } /* conv12_ */
1174 
1175 /* **************************************************************** */
1176 /* * * */
1177 /* * SUBROUTINE CONV21 * */
1178 /* * * */
1179 /* * * */
1180 /* * Description : Converts a number represented in a 2x2 real * */
1181 /* * array to the form of a complex number. * */
1182 /* * * */
1183 /* * Subprograms called: none * */
1184 /* * * */
1185 /* **************************************************************** */
1186 /* Subroutine */ int conv21_(doublereal *cae, doublecomplex *cn)
1187 {
1188  /* System generated locals */
1189  doublereal d__1, d__2;
1190  doublecomplex z__1;
1191 
1192  /* Builtin functions */
1193 
1194  /* Parameter adjustments */
1195  cae -= 3;
1196 
1197  /* Function Body */
1198  if (cae[5] > 75. || cae[6] > 75.) {
1199  cn->r = 1e75, cn->i = 1e75;
1200  } else if (cae[6] < -75.) {
1201  d__1 = cae[3] * pow_dd(&c_b65, &cae[5]);
1202  z__1.r = d__1, z__1.i = 0.;
1203  cn->r = z__1.r, cn->i = z__1.i;
1204  } else {
1205  d__1 = cae[3] * pow_dd(&c_b65, &cae[5]);
1206  d__2 = cae[4] * pow_dd(&c_b65, &cae[6]);
1207  z__1.r = d__1, z__1.i = d__2;
1208  cn->r = z__1.r, cn->i = z__1.i;
1209  }
1210  return 0;
1211 } /* conv21_ */
1212 
1213 /* **************************************************************** */
1214 /* * * */
1215 /* * SUBROUTINE ECPMUL * */
1216 /* * * */
1217 /* * * */
1218 /* * Description : Multiplies two numbers which are each * */
1219 /* * represented in the form of a two by two array and returns * */
1220 /* * the solution in the same form. * */
1221 /* * * */
1222 /* * Subprograms called: EMULT, ESUB, EADD * */
1223 /* * * */
1224 /* **************************************************************** */
1225 /* Subroutine */ int ecpmul_(doublereal *a, doublereal *b, doublereal *c__)
1226 {
1227  static doublereal c2[4] /* was [2][2] */, e1, e2;
1228  static doublereal n1, n2;
1229  /* Parameter adjustments */
1230  c__ -= 3;
1231  b -= 3;
1232  a -= 3;
1233 
1234  /* Function Body */
1235  emult_(&a[3], &a[5], &b[3], &b[5], &n1, &e1);
1236  emult_(&a[4], &a[6], &b[4], &b[6], &n2, &e2);
1237  esub_(&n1, &e1, &n2, &e2, c2, &c2[2]);
1238  emult_(&a[3], &a[5], &b[4], &b[6], &n1, &e1);
1239  emult_(&a[4], &a[6], &b[3], &b[5], &n2, &e2);
1240  eadd_(&n1, &e1, &n2, &e2, &c__[4], &c__[6]);
1241  c__[3] = c2[0];
1242  c__[5] = c2[2];
1243  return 0;
1244 } /* ecpmul_ */
1245 
1246 /* **************************************************************** */
1247 /* * * */
1248 /* * SUBROUTINE ECPDIV * */
1249 /* * * */
1250 /* * * */
1251 /* * Description : Divides two numbers and returns the solution. * */
1252 /* * All numbers are represented by a 2x2 array. * */
1253 /* * * */
1254 /* * Subprograms called: EADD, ECPMUL, EDIV, EMULT * */
1255 /* * * */
1256 /* **************************************************************** */
1257 /* Subroutine */ int ecpdiv_(doublereal *a, doublereal *b, doublereal *c__)
1258 {
1259  static doublereal b2[4] /* was [2][2] */, c2[4] /* was [2][2] */, e1,
1260  e2, e3;
1261  static doublereal n1, n2, n3;
1262  /* Parameter adjustments */
1263  c__ -= 3;
1264  b -= 3;
1265  a -= 3;
1266 
1267  /* Function Body */
1268  b2[0] = b[3];
1269  b2[2] = b[5];
1270  b2[1] = b[4] * -1.;
1271  b2[3] = b[6];
1272  ecpmul_(&a[3], b2, c2);
1273  emult_(&b[3], &b[5], &b[3], &b[5], &n1, &e1);
1274  emult_(&b[4], &b[6], &b[4], &b[6], &n2, &e2);
1275  eadd_(&n1, &e1, &n2, &e2, &n3, &e3);
1276  ediv_(c2, &c2[2], &n3, &e3, &c__[3], &c__[5]);
1277  ediv_(&c2[1], &c2[3], &n3, &e3, &c__[4], &c__[6]);
1278  return 0;
1279 } /* ecpdiv_ */
1280 
1281 /* Main program */ int MAIN__(void)
1282 {
1283  /* Format strings */
1284  static char fmt_300[] = "(1x,\002 RESULT FROM CONHYP=\002,1p,2d25.12)";
1285  static char fmt_301[] = "(1x,\002 EXPECTED RESULT=\002,3x,1p,2d25.12)";
1286 
1287  /* System generated locals */
1288  doublecomplex z__1;
1289  olist o__1;
1290 
1291  /* Builtin functions */
1292 
1293  /* Local variables */
1294  static doublecomplex a, b, z__;
1295  static integer lnchf;
1296  static doublereal ai, ar;
1297  static integer ip;
1298  static doublecomplex chf;
1299 
1300  /* Fortran I/O blocks */
1301  static cilist io___88 = { 0, 5, 0, 0, 0 };
1302  static cilist io___90 = { 0, 6, 0, 0, 0 };
1303  static cilist io___91 = { 0, 5, 0, 0, 0 };
1304  static cilist io___93 = { 0, 6, 0, 0, 0 };
1305  static cilist io___94 = { 0, 5, 0, 0, 0 };
1306  static cilist io___96 = { 0, 6, 0, 0, 0 };
1307  static cilist io___97 = { 0, 5, 0, 0, 0 };
1308  static cilist io___99 = { 0, 6, 0, 0, 0 };
1309  static cilist io___100 = { 0, 5, 0, 0, 0 };
1310  static cilist io___102 = { 0, 6, 0, 0, 0 };
1311  static cilist io___104 = { 0, 6, 0, fmt_300, 0 };
1312  static cilist io___107 = { 0, 6, 0, fmt_301, 0 };
1313 
1314 
1315 
1316  o__1.oerr = 0;
1317  o__1.ounit = 5;
1318  o__1.ofnmlen = 4;
1319  o__1.ofnm = "DATA";
1320  o__1.orl = 0;
1321  o__1.osta = "OLD";
1322  o__1.oacc = 0;
1323  o__1.ofm = 0;
1324  o__1.oblnk = 0;
1325  f_open(&o__1);
1326  s_rsle(&io___88);
1327  do_lio(&c__7, &c__1, (char *)&a, (ftnlen)sizeof(doublecomplex));
1328  e_rsle();
1329  s_wsle(&io___90);
1330  do_lio(&c__9, &c__1, " A= ", (ftnlen)4);
1331  do_lio(&c__7, &c__1, (char *)&a, (ftnlen)sizeof(doublecomplex));
1332  e_wsle();
1333  s_rsle(&io___91);
1334  do_lio(&c__7, &c__1, (char *)&b, (ftnlen)sizeof(doublecomplex));
1335  e_rsle();
1336  s_wsle(&io___93);
1337  do_lio(&c__9, &c__1, " B= ", (ftnlen)4);
1338  do_lio(&c__7, &c__1, (char *)&b, (ftnlen)sizeof(doublecomplex));
1339  e_wsle();
1340  s_rsle(&io___94);
1341  do_lio(&c__7, &c__1, (char *)&z__, (ftnlen)sizeof(doublecomplex));
1342  e_rsle();
1343  s_wsle(&io___96);
1344  do_lio(&c__9, &c__1, " Z= ", (ftnlen)4);
1345  do_lio(&c__7, &c__1, (char *)&z__, (ftnlen)sizeof(doublecomplex));
1346  e_wsle();
1347  s_rsle(&io___97);
1348  do_lio(&c__3, &c__1, (char *)&lnchf, (ftnlen)sizeof(integer));
1349  e_rsle();
1350  s_wsle(&io___99);
1351  do_lio(&c__9, &c__1, " LNCHF= ", (ftnlen)8);
1352  do_lio(&c__3, &c__1, (char *)&lnchf, (ftnlen)sizeof(integer));
1353  e_wsle();
1354  s_rsle(&io___100);
1355  do_lio(&c__3, &c__1, (char *)&ip, (ftnlen)sizeof(integer));
1356  e_rsle();
1357  s_wsle(&io___102);
1358  do_lio(&c__9, &c__1, " IP= ", (ftnlen)5);
1359  do_lio(&c__3, &c__1, (char *)&ip, (ftnlen)sizeof(integer));
1360  e_wsle();
1361  conhyp_(&z__1, &a, &b, &z__, &lnchf, &ip);
1362  chf.r = z__1.r, chf.i = z__1.i;
1363  s_wsfe(&io___104);
1364  do_fio(&c__2, (char *)&chf, (ftnlen)sizeof(doublereal));
1365  e_wsfe();
1366  ar = 2.31145634403113e-12;
1367  ai = -1.96169649634905e-11;
1368  s_wsfe(&io___107);
1369  do_fio(&c__1, (char *)&ar, (ftnlen)sizeof(doublereal));
1370  do_fio(&c__1, (char *)&ai, (ftnlen)sizeof(doublereal));
1371  e_wsfe();
1372  s_stop("", (ftnlen)0);
1373  return 0;
1374 } /* MAIN__ */
1375 
1376 /* Main program alias */ int sample_ (void) { return MAIN__ (); }
1377 #ifdef __cplusplus
1378 }
1379 #endif
1380 
ecpdiv_
int ecpdiv_(doublereal *a, doublereal *b, doublereal *c__)
Definition: TOMS_707.C:1257
d_lg10
double d_lg10(const doublereal *)
s_wsfe
integer s_wsfe(cilist *)
d_sign
double d_sign(const doublereal *, const doublereal *)
x
const Vector< T > const Vector< T > & x
Definition: LM_fit.h:97
olist::oblnk
char * oblnk
Definition: f2c.h:100
z_abs
double z_abs(const doublecomplex *)
c__2
static integer c__2
Definition: TOMS_707.C:41
do_lio
integer do_lio(integer *, integer *, char *, ftnlen)
log
TBCI__ cplx< T > log(const TBCI__ cplx< T > &z)
Definition: cplx.h:771
e_wsfe
integer e_wsfe()
conv21_
int conv21_(doublereal *cae, doublecomplex *cn)
Definition: TOMS_707.C:1186
atan2
double atan2(const double, const double)
olist::ounit
ftnint ounit
Definition: f2c.h:93
armult_
int armult_(const doublereal *a, const doublereal *b, doublereal *c__, const integer *l, const doublereal *rmax)
Definition: TOMS_707.C:746
s_wsle
integer s_wsle(cilist *)
olist::ofnmlen
ftnlen ofnmlen
Definition: f2c.h:95
ediv_
int ediv_(doublereal *n1, doublereal *e1, doublereal *n2, doublereal *e2, doublereal *nf, doublereal *ef)
Definition: TOMS_707.C:1038
eadd_
int eadd_(doublereal *n1, doublereal *e1, doublereal *n2, doublereal *e2, doublereal *nf, doublereal *ef)
Definition: TOMS_707.C:1060
c__7
static integer c__7
Definition: TOMS_707.C:37
MAIN__
int MAIN__(void)
Definition: TOMS_707.C:1281
aradd_
int aradd_(const doublereal *a, const doublereal *b, doublereal *c__, const integer *l, const doublereal *rmax)
Definition: TOMS_707.C:432
max
T max(const FS_Vector< dims, T > &fv)
Definition: fs_vector.h:594
olist::oerr
flag oerr
Definition: f2c.h:92
olist::ofm
char * ofm
Definition: f2c.h:98
ecpmul_
int ecpmul_(doublereal *a, doublereal *b, doublereal *c__)
Definition: TOMS_707.C:1225
doublecomplex::r
doublereal r
Definition: f2c.h:34
do_fio
integer do_fio(integer *, char *, ftnlen)
olist
Definition: f2c.h:91
olist::oacc
char * oacc
Definition: f2c.h:97
cmpsub_
int cmpsub_(doublereal *ar, doublereal *ai, doublereal *br, doublereal *bi, doublereal *cr, doublereal *ci, integer *l, doublereal *rmax)
Definition: TOMS_707.C:856
arydiv_
int arydiv_(const doublereal *ar, const doublereal *ai, const doublereal *br, const doublereal *bi, doublecomplex *c__, const integer *l, const integer *lnchf, const doublereal *rmax, const integer *bit)
Definition: TOMS_707.C:919
d_int
double d_int(const doublereal *)
c_b65
static doublereal c_b65
Definition: TOMS_707.C:36
c_b8
static doublereal c_b8
Definition: TOMS_707.C:34
f_open
integer f_open(olist *)
c_b53
static doublereal c_b53
Definition: TOMS_707.C:35
return
return
Definition: LM_fit.h:82
VOID
#define VOID
Definition: f2c.h:146
cilist
Definition: f2c.h:72
conv12_
int conv12_(doublecomplex *cn, doublereal *cae)
Definition: TOMS_707.C:1130
sin
TBCI__ cplx< T > sin(const TBCI__ cplx< T > &z)
Definition: cplx.h:776
b
F_TMatrix< T > b
Definition: f_matrix.h:736
olist::ofnm
char * ofnm
Definition: f2c.h:94
doublereal
double doublereal
Definition: f2c.h:32
c__1
static integer c__1
Definition: TOMS_707.C:38
e_wsle
integer e_wsle()
s_stop
int s_stop(const char *, ftnlen)
pow_di
double pow_di(const doublereal *, const integer *)
doublecomplex::i
doublereal i
Definition: f2c.h:34
integer
int integer
barf [ba:rf] 2.
Definition: f2c.h:27
esub_
int esub_(doublereal *n1, doublereal *e1, doublereal *n2, doublereal *e2, doublereal *nf, doublereal *ef)
Definition: TOMS_707.C:1107
d_imag
double d_imag(const doublecomplex *)
d_nint
double d_nint(const doublereal *)
conhyp_
void conhyp_(doublecomplex *ret_val, const doublecomplex *a, const doublecomplex *b, const doublecomplex *z, const int *lnchf, const int *ip)
Definition: TOMS_707.C:121
pow_dd
double pow_dd(const doublereal *, const doublereal *)
e_rsle
integer e_rsle()
store_
doublereal store_(doublereal *x)
Definition: TOMS_707.C:230
bits_
integer bits_(void)
Definition: TOMS_707.C:205
olist::osta
char * osta
Definition: f2c.h:96
abs
#define abs(x)
Definition: f2c.h:178
c__3
static integer c__3
Definition: TOMS_707.C:40
arsub_
int arsub_(const doublereal *a, const doublereal *b, doublereal *c__, const integer *l, const doublereal *rmax)
Definition: TOMS_707.C:707
y
const Vector< T > Vector< T > Vector< T > Vector< T > & y
Definition: LM_fit.h:172
c__9
static integer c__9
Definition: TOMS_707.C:39
cmpmul_
int cmpmul_(const doublereal *ar, const doublereal *ai, const doublereal *br, const doublereal *bi, doublereal *cr, doublereal *ci, const integer *l, const doublereal *rmax)
Definition: TOMS_707.C:885
s_rsle
integer s_rsle(cilist *)
emult_
int emult_(doublereal *n1, doublereal *e1, doublereal *n2, doublereal *e2, doublereal *nf, doublereal *ef)
Definition: TOMS_707.C:1016
sample_
int sample_(void)
Definition: TOMS_707.C:1376
a
const unsigned TMatrix< T > const Matrix< T > * a
Definition: matrix_kernels.h:25
stcom_
struct @4 stcom_
ftnlen
long int ftnlen
Definition: f2c.h:67
z_div
void z_div(doublecomplex *, const doublecomplex *, const doublecomplex *)
chgf_
VOID chgf_(doublecomplex *ret_val, const doublecomplex *a, const doublecomplex *b, const doublecomplex *z__, const integer *l, const integer *lnchf)
Definition: TOMS_707.C:278
olist::orl
ftnint orl
Definition: f2c.h:99
stcom_1
#define stcom_1
Definition: TOMS_707.C:30
cmpadd_
int cmpadd_(const doublereal *ar, const doublereal *ai, const doublereal *br, const doublereal *bi, doublereal *cr, doublereal *ci, const integer *l, const doublereal *rmax)
Definition: TOMS_707.C:825
Generated by   doxygen 1.8.5