



		    asin_.pl1                       10/03/83  1719.9rew 10/03/83  1004.8       12888



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1974 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

asin_: procedure (value) returns (float binary (27));

/*     compute the arcsine, arccosine, or arctangent of a single-precision floating-point number     */

declare	value float binary (27),
	code_ entry(fixed bin),
	(abs, atan, atand, sqrt) builtin;

	if abs(value) > 1.e0 then go to out_of_range;
	return(atan(value, sqrt(-value*value+1.e0)));

acos_:	entry(value) returns(float bin(27));
	if abs(value) > 1.e0
	then do;

out_of_range:
	     call code_(58);
	     return (0e0);
	     end;
	return(atan(sqrt(-value*value+1.e0), value));

asind_:	entry(value) returns(float bin(27));
	if abs(value) > 1.e0 then go to out_of_range;
	return(atand(value, sqrt(-value*value+1.e0)));

acosd_:	entry(value) returns(float bin(27));
	if abs(value) > 1.e0 then go to out_of_range;
	return(atand(sqrt(-value*value+1.e0), value));

atan_:	entry(value) returns(float bin(27));
	return(atan(value));

atand_:	entry(value) returns(float bin(27));
	return(atand(value));
	end;




		    asinh_.pl1                      10/03/83  1719.9rew 10/03/83  1004.8       11349



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

/* Modified by A. Downing to use round builtin function 07/16/73 */

asinh_: procedure (number) returns (float binary (27));

/*      compute the hyperbolic arcsine or arccosine of a single-precision floating-point number      */
declare (number, r) float binary (27),
        (f, n, p) float binary (63);
dcl  code_ ext entry (fixed bin),
     (abs, log, round, sqrt) builtin;
	p = 1.e0;
	n = number;
	f = abs (n);
	r = f;
	if f < 5.e-5 then go to negate;
asinhs:	if f >= 11586.e0 then p = f;
	else p = sqrt (f*f + p);
	p = p + f - 1.e0;
	r = round (p, 28);
	r = log (r+1);
negate:	if n < 0.0e0 then r = -r;
	return (r);
acosh_: entry (number) returns (float binary (27));
	p = -1.e0;
	n, f = abs (number);
	if f >= 1.e0 then go to asinhs;
	call code_ (34);
	return (0.0e0);
     end asinh_;




		    atan2_.pl1                      10/03/83  1719.9rew 10/03/83  1004.8        7272



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

atan2_: procedure (value1, value2) returns (float binary (27));

/*      compute the arctangent of two single-precision floating-point numbers	*/
declare	(value1, value2) float binary (27),
	(atan, atand) builtin;

	return(atan(value1, value2));

atand2_:	entry(value1, value2) returns(float bin(27));
	return(atand(value1, value2));
	end;




		    cabs_.pl1                       10/03/83  1719.9rew 10/03/83  1004.8        5859



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

/* calls to round_ removed by A. downing 07/16/73 */
cabs_: proc (number) returns (float bin (27));

declare	number complex float bin (27),
	abs builtin;

	return(abs(number));
	end;




		    casin_.pl1                      10/03/83  1719.9rew 10/03/83  1004.8        9819



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

casin_: proc (number) returns (complex float bin (27));

dcl (number, a, b, c) complex float bin (27);
dcl	(imag, log, real, sqrt) builtin;

	real (a) = -imag (number) + 0.0e0;
	imag (a) = real (number) + 0.0e0;

	b = 1.0e0;

trig:	
	c = -1.0e0i;

ret:	
	return (log (sqrt (a*a+b)+a)*c);

cacos_: entry (number) returns (complex float bin (27));

	a = number +0.0e0;
	b = -1.0e0;

	goto trig;

casinh_: entry (number) returns (complex float bin (27));

	b = 1.0e0;

hyper:	
	a = number +0.0e0;
	c = 1.0e0;

	goto ret;

cacosh_: entry (number) returns (complex float bin (27));

	b = -1.0e0;

	goto hyper;

     end casin_;




		    catan_.pl1                      10/03/83  1719.9rew 10/03/83  1004.5        8928



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

catan_: proc (number) returns (complex float bin (27));

dcl (number, a, b, c) complex float bin (27);
dcl	(imag, log, real) builtin;

dcl  code_ ext entry (fixed bin (17));

	b = 1.0e0i;
	c = 0.5e0i;

atans:	
	a = number + 0.0e0;

	if a = b
	     then do;
err:	     
	     call code_ (59);
	     return ((real (a)+imag (a))*170141182.0e30);
	end;

	if a = -b then goto err;

	return (log ((b+a)/ (b-a))*c);

catanh_: entry (number) returns (complex float bin (27));

	b = 1.0e0;
	c = 0.5e0;

	goto atans;

     end catan_;




		    cexp_.pl1                       10/03/83  1719.9rew 10/03/83  1004.8        9162



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

cexp_: proc (number) returns (complex float bin (27));

dcl (number, a, r) complex float bin (27);
dcl	(abs, cos, exp, imag, real, sin) builtin;
dcl  code_ ext entry (fixed bin (17));

	a = number + 0.0e0;

	if real (a)>88.028e0
	     then do;
	     call code_ (26);

	     real (a) = 170141182.0e30;
	end;
	else real (a) = exp (real (a));

	if abs (imag (a)) >= 134217728.0e0
	     then do;
	     call code_ (27);

	     return (0.0e0);
	end;

	imag (r) = sin (imag (a))*real (a);
	real (r) = cos (imag (a))*real (a);

	return (r);

     end cexp_;




		    cfdp_.pl1                       10/03/83  1719.9rew 10/03/83  1004.8        5319



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

cfdp_: proc (a, b) returns (complex float bin (27));

dcl	(a, b) complex float bin(27);

	return (a / b);

     end cfdp_;




		    cfmp_.pl1                       10/03/83  1719.9rew 10/03/83  1004.8        5328



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

cfmp_: proc (a, b) returns (complex float bin (27));

dcl (a, b) complex float bin (27);

	return (a * b);

     end cfmp_;




		    clog_.pl1                       10/03/83  1719.9rew 10/03/83  1004.8        7632



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

clog_: proc (number) returns (complex float bin (27));

dcl (number, a, b) complex float bin (27);

dcl  code_ ext entry (fixed bin (17)),
	(abs, atan, imag, log, real) builtin;

	a = number;

	if a = 0.0e0
	     then do;
	     call code_ (28);
	     return (-170141182.0e30);
	end;

	imag (b) = atan(imag (a), real (a));
	real (b) = log (abs (a));

	return (b);

     end clog_;




		    code_.pl1                       10/03/83  1719.9rew 10/03/83  1004.8       12375



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

code_: proc (n);

dcl  n fixed bin;					/* error message number */

dcl (p, q, called_pt, caller_pt) ptr,
    (called_size, caller_size) fixed bin,
     bn bit (18);

dcl  cu_$stack_frame_ptr entry returns (ptr),
     pl1_frame_$name entry (ptr, ptr, fixed bin),
     math_error_ entry (fixed bin, char (*) aligned, char (*) aligned, ptr);

dcl 1 frame based,
    2 skip (8) ptr,
    2 back ptr,
    2 forward ptr,
    2 return ptr;

dcl  called char (called_size) aligned based (called_pt),
     caller char (caller_size) aligned based (caller_pt);

	p = cu_$stack_frame_ptr() -> frame.back;
	bn = baseno (p -> frame.return);

	do while (baseno (p -> frame.return) = bn);
	     p = p -> frame.back;
	end;

	q = p -> frame.forward;

	call pl1_frame_$name (p, caller_pt, caller_size);
	call pl1_frame_$name (q, called_pt, called_size);

	call math_error_ (n, (caller), called, p -> frame.return);

     end;




		    csin_.pl1                       10/03/83  1719.9rew 10/03/83  1004.5       18234



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

csin_: proc (number) returns (complex float bin (27));

dcl (number, a, b) complex float bin (27);
dcl (sinx, cosx, sinhy, coshy) float bin (27),
	(abs, cos, cosh, imag, real, sin, sinh) builtin,
     i fixed bin (17);

dcl  code_ ext entry (fixed bin (17));

	i = 1;

csins:	
	a = number;

test:	
	if abs (imag (a))>88.028e0
	     then do;
	     call code_ (29);

	     sinhy,
	     coshy = 170141182.0e30;
	end;
	else do;
	     sinhy = sinh (imag (a));
	     coshy = cosh (imag (a));
	end;

	if abs (real (a)) >= 170141182.0e30
	     then do;
	     call code_ (30);

	     return (0.0e0);
	end;

	sinx = sin (real (a));
	cosx = cos (real (a));

	if i>0
	     then do;
	     real (a) = sinx*coshy;
	     imag (a) = cosx*sinhy;
	end;
	else if i<0
	     then do;
	     real (a) = cosx*sinhy;
	     imag (a) = -sinx*coshy;

	     i = -i;
	end;

	if i = 1 then return (a);

	real (b) = cosx*coshy;
	imag (b) = -sinx*sinhy;

	if i = 0 then goto ret;
	if b = 0.0e0
	     then do;
	     call code_ (36);

	     return (170141182.0e30*sinx);
	end;

	b = a/b;

ret:	
	return (b);

ccos_: entry (number) returns (complex float bin (27));

	i = 0;

	goto csins;

ctan_: entry (number) returns (complex float bin (27));

	i = 2;

	goto csins;

csinh_: entry (number) returns (complex float bin (27));

	i = -1;

csinhs:	
	real (a) = -imag (number);
	imag (a) = real (number);

	goto test;

ccosh_: entry (number) returns (complex float bin (27));

	i = 0;
	goto csinhs;

ctanh_: entry (number) returns (complex float bin (27));

	i = -2;
	goto csinhs;

     end csin_;




		    csqrt_.pl1                      10/03/83  1719.9rew 10/03/83  1004.8        9135



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

csqrt_: proc (number) returns (complex float bin (27));

dcl (number, a, b) complex float bin (27);

	a = number;
	real (b) = abs (real (a));

	if imag (a) = 0.0e0
	then do;
	     real (b) = sqrt (real (b));
	     imag (b) = 0.0e0;
	end;
	else do;
	     real (b) = sqrt ((abs (a) + real (b))*0.5e0);

	     if real (a)<0.0e0
		& imag (a)<0.0e0
		then real (b) = -real (b);

	     imag (b) = imag (a) * 0.5e0/real (b);
	end;

	if real (a) >= 0.0e0
	     then return (b);

	real (a) = imag (b);
	imag (a) = real (b);

	return (a);

     end csqrt_;




		    cxp12_.pl1                      10/03/83  1719.9rew 10/03/83  1004.8        6129



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

cxp12_: proc (base, exponent) returns (complex float bin (27));

dcl base complex float bin (27),
     exponent fixed bin (71),
     f complex float bin(63);

	f = base;
	return (f ** exponent);

     end cxp12_;




		    cxp1_.pl1                       10/03/83  1719.9rew 10/03/83  1004.9        6129



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

cxp1_: proc (base, exponent) returns (complex float bin (27));

dcl base complex float bin (27),
     f complex float bin (63),
     exponent fixed bin (17);


	f = base;
	return (f ** exponent);

     end cxp1_;




		    cxp2_.pl1                       10/03/83  1719.9rew 10/03/83  1004.5        6066



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

cxp2_: proc (base, exponent) returns (complex float bin (27));

dcl (base, exponent) complex float bin (27),
    (a, b) complex float bin (63);

	a = base;
	b = exponent;
	return (a ** b);

     end cxp2_;




		    dasin_.pl1                      10/03/83  1719.9rew 10/03/83  1004.9       14004



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1974 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

dasin_: procedure (number) returns (float binary (63));

/*     compute the arcsine, arccosine, or arctangent of a double-precision floating-point number     */

declare	(number, value) float binary (63),
	code_ entry(fixed bin),
	(abs, atan, atand, sqrt) builtin;

	value = number;
	if abs(value) > 1.e0 then go to out_of_range;
	return(atan(value, sqrt(-value*value+1.e0)));

dacos_:	entry(number) returns(float bin(63));
	value = number;
	if abs(value) > 1.e0
	then do;

out_of_range:
	     call code_(58);
	     return (0e0);
	     end;
	return(atan(sqrt(-value*value+1.e0), value));

dasind_:	entry(number) returns(float bin(63));
	value = number;
	if abs(value) > 1.e0 then go to out_of_range;
	return(atand(value, sqrt(-value*value+1.e0)));

dacosd_:	entry(number) returns(float bin(63));
	value = number;
	if abs(value) > 1.e0 then go to out_of_range;
	return(atand(sqrt(-value*value+1.e0), value));

datan_:	entry(number) returns(float bin(63));
	value = number;
	return(atan(value));

datand_:	entry(number) returns(float bin(63));
	value = number;
	return(atand(value));
	end;




		    dasinh_.pl1                     10/03/83  1719.9rew 10/03/83  1004.9       11331



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

dasinh_: procedure (number) returns (float binary (63));

/*       compute the hyperbolic arcsine or arccosine of a double-precision floating-point number       */
declare (number, f, n, p, q, r) float binary (63),
	(abs, log, sqrt) builtin;
dcl  code_ ext entry (fixed bin);
	n = number;
	f = abs (n);
	if f < 1.e-10 then go to negate;
	if f >= 3037000500.e0 then go to setup;
	p = f*f;
	q = sqrt (p + 1.e0);
	r = q - 1.e0;
	p = (r*r + p) * 0.5e0 / q;
loner:	f = log (f + p + 1);
negate:	if n < 0.0e0 then f = -f;
	return (f);
dacosh_: entry (number) returns (float binary (63));
	f, n = abs (number);
if f < 1.e0 then err: do;
	     call code_ (35);
	     return (0.0e0);
	end err;
setup:	p = f - 1.e0;
	if f < 3037000500.e0 then f = sqrt ((f + 1.e0) * p);
	go to loner;
     end dasinh_;




		    datan2_.pl1                     10/03/83  1719.9rew 10/03/83  1004.9        7290



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

datan2_: procedure (value1, value2) returns (float binary (63));

/*      compute the arctangent of two double-precision floating-point numbers	*/
declare	(value1, value2) float binary (63),
	(atan, atand) builtin;

	return(atan(value1, value2));

datand2_:	entry(value1, value2) returns(float bin(63));
	return(atand(value1, value2));
	end;




		    dcabs_.pl1                      10/03/83  1719.9rew 10/03/83  1004.5        7596



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

dcabs_: proc (number) returns (float bin (63));

dcl  number complex float bin (63),
	(abs, imag, real, sqrt) builtin,
    (r, x, y) float bin (63);

	r,
	x = abs (real (number));
	y = abs (imag (number));

	if y<x
	     then do;
	     x = y;
	     y = r;
	end;

	if y ^= 0.0e0
	     then do;
	     r = x/y;
	     y = sqrt (r*r+1.0e0)*y;
	end;

	return (y);

     end dcabs_;




		    dcasin_.pl1                     10/03/83  1719.9rew 10/03/83  1004.5        9594



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

dcasin_: proc (number) returns (complex float bin (63));

dcl (number, a, b, c) complex float bin (63);
dcl	(imag, log, real, sqrt) builtin;

	real (a) = -imag (number);
	imag (a) = real (number);

	b = 1.0e0;

trig:	
	c = -1.0e0i;

ret:	
	return (log (sqrt (a*a+b)+a)*c);

dcacos_: entry (number) returns (complex float bin (63));

	a = number;
	b = -1.0e0;

	goto trig;

dcasinh_: entry (number) returns (complex float bin (63));

	b = 1.0e0;

hyper:	
	a = number;

	c = 1.0e0;

	goto ret;

dcacosh_: entry (number) returns (complex float bin (63));

	b = -1.0e0;
	goto hyper;

     end dcasin_;




		    dcatan_.pl1                     10/03/83  1719.9rew 10/03/83  1004.9        8883



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

dcatan_: proc (number) returns (complex float bin (63));

dcl (number, a, b, c) complex float bin (63);
dcl	(imag, log, real) builtin;

dcl  code_ ext entry (fixed bin (17));

	b = 1.0e0i;
	c = 0.5e0i;

atans:	
	a = number;

	if a = b
	     then do;
err:	     
	     call code_ (32);
	     return ((real (a)+imag (a))*170141182.0e30);
	end;

	if a = -b then goto err;

	return (log ((b+a)/ (b-a))*c);

dcatanh_: entry (number) returns (complex float bin (63));

	b = 1.0e0;
	c = 0.5e0;

	goto atans;

     end dcatan_;




		    dcexp_.pl1                      10/03/83  1719.9rew 10/03/83  1004.9        9198



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

dcexp_: proc (number) returns (complex float bin (63));

dcl (number, a, r) complex float bin (63);

dcl	(abs, cos, exp, imag, real, sin) builtin;

dcl  code_ ext entry (fixed bin (17));

	a = number;

	if real (a)>88.028e0
	     then do;
	     call code_ (68);

	     real (a) = 170141182.0e30;
	end;
	else real (a) = exp (real (a));

	if abs (imag (a)) >= 18104398509481984.0e0
	     then do;
	     call code_ (69);

	     return (0.0e0);
	end;

	imag (r) = sin (imag (a))*real (a);
	real (r) = cos (imag (a))*real (a);

	return (r);

     end dcexp_;




		    dcfdp_.pl1                      10/03/83  1719.9rew 10/03/83  1004.9        5346



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

dcfdp_: proc (a, b) returns (complex float bin (63));

dcl (a, b) complex float bin (63);

	return (a / b);

     end dcfdp_;




		    dcfmp_.pl1                      10/03/83  1719.9rew 10/03/83  1004.9        5346



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

dcfmp_: proc (a, b) returns (complex float bin (63));

dcl (a, b) complex float bin (63);

	return (a * b);

     end dcfmp_;




		    dclog_.pl1                      10/03/83  1719.9rew 10/03/83  1004.9        7632



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

dclog_: proc (number) returns (complex float bin (63));

dcl (number, a, b) complex float bin (63);

dcl  code_ ext entry (fixed bin (17)),
	(abs, atan, imag, log, real) builtin;

	a = number;
	if a = 0.0e0
	     then do;
	     call code_ (31);
	     return (170141182.0e30);
	end;

	imag (b) = atan(imag (a), real (a));
	real (b) = log (abs (a));

	return (b);

     end dclog_;




		    dcsin_.pl1                      10/03/83  1719.9rew 10/03/83  1004.9       18342



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

dcsin_: proc (number) returns (complex float bin (63));

dcl (number, a, b) complex float bin (63),
    (sinx, cosx, sinhy, coshy) float bin (63),
	(abs, cos, cosh, imag, real, sin, sinh) builtin,
     i fixed bin (17);

dcl  code_ ext entry (fixed bin (17));

	i = 1;

csins:	
	a = number;

test:	
	if abs (imag (a))>88.018e0
	     then do;
	     call code_ (61);

	     sinhy,
	     coshy = 170141182.0e30;
	end;
	else do;
	     sinhy = sinh (imag (a));
	     coshy = cosh (imag (a));
	end;

	if abs (real (a)) >= 18104398509481984.0e0
	     then do;
	     call code_ (62);

	     return (0.0e0);
	end;

	sinx = sin (real (a));
	cosx = cos (real (a));

	if i>0
	     then do;
	     real (a) = sinx*coshy;
	     imag (a) = cosx*sinhy;
	end;
	else if i<0
	     then do;
	     real (a) = cosx*sinhy;
	     imag (a) = -sinx*coshy;

	     i = -i;
	end;

	if i = 1 then return (a);

	real (b) = cosx*coshy;
	imag (b) = -sinx*sinhy;

	if i = 0 then goto ret;

	if b = 0.0e0
	     then do;
	     call code_ (64);
	     return (170141182.0e30*sinx);
	end;

	b = a/b;

ret:	
	return (b);

dccos_: entry (number) returns (complex float bin (63));

	i = 0;
	goto csins;

dctan_: entry (number) returns (complex float bin (63));

	i = 2;
	goto csins;

dcsinh_: entry (number) returns (complex float bin (63));

	i = -1;

csinhs:	
	real (a) = -imag (number);
	imag (a) = real (number);

	goto test;

dccosh_: entry (number) returns (complex float bin (63));

	i = 0;
	goto csinhs;

dctanh_: entry (number) returns (complex float bin (63));

	i = -2;
	goto csinhs;

     end dcsin_;




		    dcsqrt_.pl1                     10/03/83  1719.9rew 10/03/83  1004.9        9081



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

dcsqrt_: proc (number) returns (complex float bin (63));

dcl (number, a, b) complex float bin (63);

	a = number;
	real (b) = abs (real (a));

	if imag (a) = 0.0e0
	then do;
	     real (b) = sqrt (real (b));
	     imag (b) = 0.0e0;
	end;
	else do;
	     real (b) = sqrt ((abs (a)+real (b))*0.5e0);

	     if real (a)<0.0e0
		& imag (a)<0.0e0
		then real (b) = -real (b);

	     imag (b) = imag (a) * 0.5e0/real (b);
	end;

	if real (a) >= 0.0e0 then return (b);

	real (a) = imag (b);
	imag (a) = real (b);

	return (a);

     end dcsqrt_;




		    dcxp12_.pl1                     10/03/83  1719.9rew 10/03/83  1004.9       10098



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

dcxp12_: proc (base, exponent) returns (complex float bin (63));

dcl (base, a, f) complex float bin (63),
    (exponent, h, k, m) fixed bin (71);

dcl	(abs, divide, sign) builtin;

dcl  code_ ext entry (fixed bin (17));

	a = base;
	k = exponent;

	f = 1.0e0;
	if a = 0.0e0
	     then do;
	     if k>0 then goto clear;

	     call code_ (14-sign (k));
	     goto clear;
	end;

	if k = 0 then goto ret;
	m = abs (k);

loop:	
	h = divide(m,2,71,0);

	if h+h ^= m
	     then f = f*a;

	if h ^= 0
	     then do;
	     m = h;
	     a = a*a;
	     goto loop;
	end;

	if k<0
	     then f = 1.0e0/f;

ret:	
	return (f);

clear:	
	return (a);

     end dcxp12_;




		    dcxp1_.pl1                      10/03/83  1719.9rew 10/03/83  1004.9        9963



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

dcxp1_: proc (base, exponent) returns (complex float bin (63));

dcl (base, a, f) complex float bin (63),
    (exponent, h, k, m) fixed bin (17);

dcl	divide builtin;

dcl  code_ ext entry (fixed bin (17));

	a = base;
	k = exponent;

	f = 1.0e0;
	if a = 0.0e0
	     then do;
	     if k>0 then goto clear;

	     call code_ (14-sign (k));
	     goto clear;
	end;

	if k = 0 then goto ret;
	m = abs (k);

loop:	
	h = divide(m,2,17,0);

	if h+h ^= m
	     then f = f*a;

	if h ^= 0
	     then do;
	     m = h;
	     a = a*a;
	     goto loop;
	end;

	if k<0
	     then f = 1.0e0/f;

ret:	
	return (f);

clear:	
	return (a);

     end dcxp1_;




		    dcxp2_.pl1                      10/03/83  1719.9rew 10/03/83  1004.9        7974



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

dcxp2_: proc (base, exponent) returns (complex float bin (63));

dcl (base, exponent, a, b) complex float bin (63);

dcl  code_ ext entry (fixed bin (17));

	a = base;
	b = exponent;

	if a = 0.0e0
	     then do;
	     if real (b)>0.0e0
		& imag (b) = 0.0e0
		then goto ret;

	     call code_ (57);
	     goto ret;
	end;

	if b = 0.0e0 then return (1.0e0);

	a = exp (log (a)*b);

ret:	
	return (a);

     end dcxp2_;




		    derf_.pl1                       10/03/83  1719.9rew 10/03/83  1004.9       34542



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

derf_: procedure (number) returns (float binary (63));

/*     compute the error (or complementary error) function of a double-precision floating-point number     */
declare (number, r) float binary (63),
        (f, n, p, q) float binary (63),
	(abs, exp) builtin,
        dexerfc_ entry (float binary (63)) returns (float binary (63));
	r = 0.0e0;
erfs:	n = number + 0.0e0;
	f = abs (n);
if f < 0.47693628e0 then small: do;
	     q = 1.128379167095512574e0 * n;
	     if f < 1.e-10 then go to comp;
	     p = f*f;
	     n = (((( 0.8350702795147239592e-10 *p - 0.1089222103714857338e-8)*p + 0.1312253296380280507e-7)*p
	     - 0.1450385222315046876e-6)*p + 0.1458916900093370682e-5)*p;
	     q = ((((((((n - 0.1322751322751322751e-4)*p + 0.1068376068376068376e-3)*p - 0.7575757575757575758e-3)*p
	     + 0.4629629629629629630e-2)*p - 0.2380952380952380952e-1)*p + 0.1e0)*p
	     - 0.3333333333333333333e0)*p + 1.e0)*q;
comp:	     if r ^= 0.0e0 then q = 1.e0 - q;
	     go to finis;
	end small;
if f >= 2.5e0 then large: do;
	     if f >= 9.30630096e0 then q = 0.0e0;
	     else q = exp (-f*f) * dexerfc_ (f);
end large; else middle: do;
if f < 1.5e0 then lower: do;
		p = f - 1.e0;
		q = ((((( - 0.8445354974538876839e-12 *p + 0.2035190183702655211e-11)*p + 0.1016478011836914822e-10)*p
		- 0.3914544150228919332e-10)*p - 0.9687156253803097249e-10)*p + 0.6116760848521522320e-9)*p;
		q = ((((((( q + 0.5758182413135021396e-9)*p - 0.7972478608404360084e-8)*p + 0.1720401851653628375e-8)*p
		+ 0.8630591951220226579e-7)*p - 0.1092604901414462467e-6)*p - 0.7524484361200326476e-6)*p
		+ 0.1844279900355114423e-5)*p;
		q = (((((((((q + 0.4854967260442840621e-5)*p - 0.2100986406780402429e-4)*p - 0.1658718964594168655e-4)*p
		+ 0.1772942579639506779e-3)*p - 0.7579366392581423493e-4)*p - 0.1086769072243678816e-2)*p
		+ 0.1670704693150730120e-2)*p + 0.4233533251576121955e-2)*p - 0.1343051928086217999e-1)*p;
		q = ((((((( q - 0.4087549346349359129e-2)*p + 0.6131324019524038693e-1)*p - 0.6131324019524038693e-1)*p
		- 0.1226264803904807739e0)*p + 0.3678794411714423216e0)*p - 0.3678794411714423216e0)*p
		+ 0.1394027926403309882e0) * 1.128379167095512574e0;
end lower; else upper: do;
		p = f - 2.e0;
		q = (( - 0.1856500126366457284e-12 *p - 0.1944964938417663150e-12)*p + 0.3103020133929132565e-11)*p;
		q = (((((((((q - 0.3723550879099818846e-11)*p - 0.3426266649193350739e-10)*p + 0.1200727503356961138e-9)*p
		+ 0.1787795198978365218e-9)*p - 0.1821176414744498095e-8)*p + 0.1759567016164738904e-8)*p
		+ 0.1642444033947200479e-7)*p - 0.5324702588728536810e-7)*p - 0.5245213918278798165e-7)*p;
		q = (((((((((q + 0.6206354808105919169e-6)*p - 0.8484562971386512658e-6)*p - 0.3501612055936534974e-5)*p
		+ 0.1439484990506010484e-4)*p - 0.4634950036778170257e-5)*p - 0.9195995379200109923e-4)*p
		+ 0.2329025685851383420e-3)*p + 0.4441623187303262417e-4)*p - 0.1410013470005726578e-2)*p;
		q = ((((((((q + 0.2994461596094635826e-2)*p - 0.4070141975274262287e-3)*p - 0.1159990462953164752e-1)*p
		+ 0.3052606481455696716e-1)*p - 0.4273649074037975402e-1)*p + 0.3663127777746836059e-1)*p
		- 0.1831563888873418029e-1)*p + 0.4145534690336333682e-2) * 1.128379167095512574e0;
	end upper; end middle;
	if r = 0.0e0 then q = 1.e0 - q;
	if n < 0.0e0 then q = r - q;
finis:	return (q);
derfc_: entry (number) returns (float binary (63));
	r = 2.e0;
	go to erfs;
     end derf_;




		    dexerfc_.pl1                    10/03/83  1719.9rew 10/03/83  1004.9       16830



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

dexerfc_: procedure (number) returns (float binary (63));

/*        compute the special error function of a double-precision floating-point number        */
declare (number, f, n, hn, dh, dl, dm, ph, rh, rl, rm, sh, sl, sm, th, tl, tm) float binary (63),
	(abs, exp, erfc) builtin;
dcl  code_ ext entry (fixed bin);
	n = number + 0.0e0;
if n < -9.30630096e0 then err: do;
	     call code_ (66);
	     return (170141182.e30);
	end err;
	f = abs (n);
	if f < 2.5e0 then th = erfc (n) * exp (f*f);
else large: do;
	     rm, th = 0.5e0 / f;
	     if f >= 1.e11 then go to done;
	     ph = f;
	     rl = 0.0e0;
	     hn = 0.5e0;
	     th = -1.e10;
loop:	     dm = hn / ph;
	     ph = dm + f;
	     rh = (rl * dm + rm * f) / ph;
	     dl = rl - rm;
	     dh = rh - rm;
	     dm = dh + dl;
	     if dm = 0.0e0 then go to dvc1;
	     sh = (dh/dm) * dl + rm;
	     if hn < 1.25e0 then go to step;
	     dl = sl - sm;
	     dh = sh - sm;
	     dm = dh + dl;
	     if dm = 0.0e0 then go to dvc2;
	     th = (dh/dm) * dl + sm;
	     if th = tm then go to done;
	     if th = tl then go to done;
step:	     hn = hn + 0.5e0;
	     rl = rm;
	     sl = sm;
	     tl = tm;
	     rm = rh;
	     sm = sh;
	     tm = th;
	     go to loop;
dvc1:	     sh = rh;
dvc2:	     th = sh;
done:	     th = 1.128379167095512574e0 * th;
	     if n < 0.0e0 then th = 2.e0 * exp (f*f) - th;
	end large;
	return (th);
     end dexerfc_;




		    dexp_.pl1                       10/03/83  1719.9rew 10/03/83  1004.9        6660



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

/* calls to round_ removed on 07/16/73 by A. Downing */
dexp_: procedure (number) returns (float binary (63));

/*     compute the exponential of a double-precision floating-point number    */
declare (number) float binary (63),
	exp builtin;

	return(exp(number));
	end;




		    diexp_.pl1                      10/03/83  1719.9rew 10/03/83  1004.9        9990



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

diexp_: procedure (base, exponent) returns (fixed binary (71));
						/*     compute integer base ** integer exponent     */
declare (base, exponent, h, i, j) fixed binary (17),
         f fixed binary (71),
         code_ entry (fixed binary);
	i = base;
	j = exponent;
	f = 1;
if i = 0 then test: do;
	     if j > 0 then go to clear;
	     call code_ (5 - sign (j));
	     go to clear;
	end test;
	if j = 0 then go to finis;
	if abs (i) = 1 then j = mod (j, 2);
else if j < 0 then clear: return (0);
loop:	h = divide (j, 2, 17, 0);
	if h+h ^= j then f = f*i;
if h = 0 then finis: return (f);
	j = h;
	i = i*i;
	go to loop;
     end diexp_;




		    dlog_.pl1                       10/03/83  1719.9rew 10/03/83  1005.0       21078



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

/* calls to round_ removed 07/16/73 by A. Downing */
dlog_: procedure (number) returns (float binary (63));

/*     compute the logarithm or hyperbolic arctangent of a double-precision floating-point number     */
declare (number, a, f, n, p, r) float binary (63),
         i fixed binary (17);
dcl  code_ ext entry (fixed bin),
     1 word aligned based,
     2 exponent fixed bin (7) unal;
declare	(abs, addr, log, log10, log2, sign) builtin;

	return(log(number));		/* Natural log of value. */

long:	i = addr (f) -> word.exponent;
	addr (f) -> word.exponent = 0;
if f < 0.7071067811865475244e0 then lower: do;
	     a = 0.5946035575013605334e0;
	     n = 0.75e0;
end lower; else upper: do;
	     a = 0.840896415253714543e0;
	     n = 0.25e0;
	end upper;
	r = f + a;
	f = (f - a) / r;
	n = i - n;
short:	if abs (f) < 1.e-19 then go to finis;
	a = f*f;
	f = ((((((((0.1176470588235294118e0 * a + 0.1333333333333333333e0) * a + 0.1538461538461538462e0) * a
	+ 0.1818181818181818182e0)* a + 0.2222222222222222222e0) * a + 0.2857142857142857143e0) * a
	+ 0.4e0) * a + 0.6666666666666666667e0) * a + 2.e0) * f;
finis:	return ((0.6931471805599453094e0 * n + f) * p);

dlog2_: entry (number) returns (float binary (63));
	return(log2(number));			/* Log (2) of value. */

dlog10_: entry (number) returns (float binary (63));
	return(log10(number));			/* Log (10) of value. */

dlone_: entry (number) returns (float binary (63));
	return(log(number+1.0e0));			/* Natural log of x+1. */

datanh_: entry (number) returns (float binary (63));
	p = 0.5e0;
	f = number;
	a = abs (f);
	if a < 0.1e0
	then do;
	     n = 0.0e0;
	     go to short;
	     end;
	a = a - 1.e0;
	if a >= 0.0e0 then err2: do;
	     call code_ (sign (a) + 44);
	     if a = 0.0e0 then f = 170141182.e30 * f; else f = 0.0e0;
	     return (f);
	end err2;
	f = (1.e0 + f) / (1.e0 - f);
	go to long;
     end dlog_;




		    dmod_.fortran                   10/03/83  1719.9rew 10/03/83  1004.5        6858



c  ******************************************************
c  *                                                    *
c  *                                                    *
c  * Copyright (c) 1972 by Massachusetts Institute of   *
c  * Technology and Honeywell Information Systems, Inc. *
c  *                                                    *
c  *                                                    *
c  ******************************************************
c
   double precision function dmod_(x,y)
c
c     double precision modulus function
c
c     Modified: 30 August 1979 by C R Davis to just call the dmod intrinsic,
c          since it is now properly implemented in pl1_operators_.
c

   double precision x,y

   dmod_ = dmod (x, y)
   return
   end




		    dsin_.pl1                       10/03/83  1719.9rew 10/03/83  1005.0        8217



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

dsin_: procedure (number) returns (float binary (63));

/*     compute the sine or cosine of a double-precision floating-point number     */
declare number float binary (63),
	(cos, cosd, sin, sind) builtin;

	return(sin(number));

dcos_:	entry(number) returns(float bin(63));
	return(cos(number));

dsind_:	entry(number) returns(float bin(63));
	return(sind(number));

dcosd_:	entry(number) returns(float bin(63));
	return(cosd(number));
	end;




		    dsinh_.pl1                      10/03/83  1719.9rew 10/03/83  1005.0       13365



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

dsinh_: procedure (number) returns (float binary (63));

/*      compute the hyperbolic sine or cosine of a double-precision floating-point number      */
declare (number, b, n, p, r) float binary (63),
	(abs, exp) builtin;
dcl  code_ ext entry (fixed bin);
	n = number;
	r = abs (n);
if r >= 0.7135236e0 then large: do;
	     p = -1.e0;
	     go to sinhs;
	end large;
	if r < 1.e-10 then go to negate;
	b = r*r;
	r = ((((((((0.2811457254345520763e-14 * b + 0.7647163731819816476e-12) * b + 0.1605904383682161460e-9) * b
	+ 0.2505210838544171878e-7) * b + 0.2755731922398589065e-5) * b + 0.1984126984126984127e-3) * b
	+ 0.8333333333333333333e-2) * b + 0.1666666666666666667e0) * b + 1.e0) * n;
	go to finis;
dcosh_: entry (number) returns (float binary (63));
	n, r = abs (number);
	p = 1.e0;
sinhs: if r > 88.028e0 then err: do;
	     call code_ (50);
	     r = 170141182.e30;
	     go to negate;
	end err;
	r = exp (r);
	r = (p/r + r) * 0.5e0;
negate:	if n < 0.0e0 then r = -r;
finis:	return (r);
     end dsinh_;




		    dsqrt_.pl1                      10/03/83  1719.9rew 10/03/83  1005.0        6669



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

/* calls to round_ removed 07/16/73 by A. Downing */
dsqrt_: procedure (number) returns (float binary (63));

/*      compute the square root of a double-precision floating-point number      */
declare  number float binary (63),
	sqrt builtin;
	return(sqrt(number));
	end;




		    dtan_.pl1                       10/03/83  1719.9rew 10/03/83  1005.0        6813



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

dtan_: procedure (number) returns (float binary (63));

/*     compute the tangent of a double-precision floating-point number     */
declare number float binary (63),
	(tan, tand) builtin;

	return(tan(number));

dtand_:	entry(number) returns(float bin(63));
	return(tand(number));
	end;




		    dtanh_.pl1                      10/03/83  1719.9rew 10/03/83  1005.0       14094



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

dtanh_: procedure (number) returns (float binary (63));

/*      compute the hyperbolic tangent of a double-precision floating-point number      */
declare (number, f, n, p, q) float binary (63),
	(abs, exp) builtin;
	n = number;
	f = abs (n);
if f >= 0.55e0 then if f >= 24.5e0 then f = 1.e0; else large: do;
	     f = f + f;
	     f = exp (f);
	     p = f + 1.e0;
	     f = (f - 1.e0) / p;
end large; else if f >= 1.e-10 then small: do;
	     p = f*f;
	     q = ((((((( 0.4779477332387385297e-13 * p + 0.1147074559772972471e-10) * p
	     + 0.2087675698786809898e-8) * p + 0.2755731922398589065e-6) * p + 0.2480158730158730159e-4) * p
	     + 0.1388888888888888889e-2) * p + 0.4166666666666666667e-1) * p + 0.5e0) * p + 1.e0;
	     f = ((((((((0.2811457254345520763e-14 * p + 0.7647163731819816476e-12) * p + 0.1605904383682161460e-9) * p
	     + 0.2505210838544171878e-7) * p + 0.2755731922398589065e-5) * p + 0.1984126984126984127e-3) * p
	     + 0.8333333333333333333e-2) * p + 0.1666666666666666667e0) * p + 1.e0) * f / q;
	end small;
	if n < 0.0e0 then f = -f;
	return (f);
     end dtanh_;




		    dxp12_.pl1                      10/03/83  1719.9rew 10/03/83  1005.0        9729



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

dxp12_: procedure (base, exponent) returns (float binary (63));

declare (base, a, f) float binary (63),
        (exponent, h, k, m) fixed binary (71);
dcl  code_ ext entry (fixed bin);
dcl	(abs, divide, sign) builtin;
	a = base;
	k = exponent;
	f = 1.e0;
if a = 0.0e0 then test: do;
if k > 0 then clear: return (a);
	     call code_ (3 - sign (k));
	     go to clear;
	end test;
	if k = 0 then go to finis;
	m = abs (k);
loop:	h = divide (m, 2, 71, 0);
	if h+h ^= m then f = f*a;
	if h = 0 then go to invert;
	m = h;
	a = a*a;
	go to loop;
invert:	if k < 0 then f = 1.e0 / f;
finis:	return (f);
     end dxp12_;




		    dxp1_.pl1                       10/03/83  1719.9rew 10/03/83  1005.0        6417



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

dxp1_: procedure (base, exponent) returns (float binary (63));

/*     compute double-precision floating-point base ** integer exponent     */
declare base float binary (63),
        exponent fixed binary (17);

	return(base ** exponent);
	end;




		    dxp2_.pl1                       10/03/83  1719.9rew 10/03/83  1005.0        6237



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

dxp2_: procedure (base, exponent) returns (float binary (63));

/*     compute base ** exponent for double-precision floating-point numbers     */
declare (base, exponent) float binary (63);

	return(base ** exponent);
	end;




		    erf_.pl1                        10/03/83  1719.9rew 10/03/83  1005.0       26136



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

/* calls to round_ removed 07/16/73 by A. Downing */
erf_: procedure (number) returns (float binary (27));

/*    compute the error (or complementary error) function of a single-precision floating-point number    */
declare (number, s) float binary (27),
        (f, n, p, q, r) float binary (63),
         exerfc_ entry (float binary (27)) returns (float binary (27)),
	(abs, exp, round) builtin;
	r = 0.0e0;
erfs:	n = number;
	f = abs (n);
if f < 0.5e0 then small: do;
	     q = 1.128379167095512574e0 * n;
	     if f < 5.e-5 then go to comp;
	     p = f*f;
	     q = ((((((( - 0.1322751322751322751e-4 *p + 0.1068376068376068376e-3)*p - 0.7575757575757575758e-3)*p
	     + 0.4629629629629629630e-2)*p - 0.2380952380952380952e-1)*p + 0.1e0)*p
	     - 0.3333333333333333333e0)*p + 1.e0)*q;
comp:	     if r ^= 0.0e0 then q = 1.e0 - q;
	     go to finis;
	end small;
if f >= 2.5e0 then large: do;
	     s = f;
	     if f >= 9.30630096e0 then q = 0.0e0;
	     else q = exp (-f*f) * exerfc_ (s);
end large; else middle: do;
if f < 1.5e0 then lower: do;
		p = f - 1.e0;
		q = (((((((( 0.4854967260442840621e-5 *p - 0.2100986406780402429e-4)*p - 0.1658718964594168655e-4)*p
		+ 0.1772942579639506779e-3)*p - 0.7579366392581423493e-4)*p - 0.1086769072243678816e-2)*p
		+ 0.1670704693150730120e-2)*p + 0.4233533251576121955e-2)*p - 0.1343051928086217999e-1)*p;
		q = ((((((( q - 0.4087549346349359129e-2)*p + 0.6131324019524038693e-1)*p - 0.6131324019524038693e-1)*p
		- 0.1226264803904807739e0)*p + 0.3678794411714423216e0)*p - 0.3678794411714423216e0)*p
		+ 0.1394027926403309882e0) * 1.128379167095512574e0;
end lower; else upper: do;
		p = f - 2.e0;
		q = (((((((( 0.6206354808105919169e-6 *p - 0.8484562971386512658e-6)*p - 0.3501612055936534974e-5)*p
		+ 0.1439484990506010484e-4)*p - 0.4634950036778170257e-5)*p - 0.9195995379200109923e-4)*p
		+ 0.2329025685851383420e-3)*p + 0.4441623187303262417e-4)*p - 0.1410013470005726578e-2)*p;
		q = ((((((((q + 0.2994461596094635826e-2)*p - 0.4070141975274262287e-3)*p - 0.1159990462953164752e-1)*p
		+ 0.3052606481455696716e-1)*p - 0.4273649074037975402e-1)*p + 0.3663127777746836059e-1)*p
		- 0.1831563888873418029e-1)*p + 0.4145534690336333682e-2) * 1.128379167095512574e0;
	end upper; end middle;
	if r = 0.0e0 then q = 1.e0 - q;
	if n < 0.0e0 then q = r - q;
finis:		s = round (q, 28);
	return (s);
erfc_: entry (number) returns (float binary (27));
	r = 2.e0;
	go to erfs;
     end erf_;




		    exerfc_.pl1                     10/03/83  1719.9rew 10/03/83  1005.0       17964



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

exerfc_: procedure (number) returns (float binary (27));

/*       compute the special error function of a single-precision floating-point number       */
declare (number, r, uh, ul, um) float binary (27),
	(abs, exp, erfc) builtin,
        (f, n, hn, dh, dl, dm, ph, rh, rl, rm, sh, sl, sm, th, tl, tm) float binary (63);
dcl  code_ ext entry (fixed bin),
     round builtin;
	n, r = number + 0.0e0;
if n < -9.30630096e0 then err: do;
	     call code_ (41);
	     return (170141182.e30);
	end err;
	f = abs (n);
	if f < 2.5e0 then th = erfc (r) * exp (f*f);
else large: do;
	     rm, th = 0.5e0 / f;
	     if f >= 1.e5 then go to done;
	     ph = f;
	     rl = 0.0e0;
	     hn = 0.5e0;
	     uh = -1.e10;
loop:	     dm = hn / ph;
	     ph = dm + f;
	     rh = (rl * dm + rm * f) / ph;
	     dl = rl - rm;
	     dh = rh - rm;
	     dm = dh + dl;
	     if dm = 0.0e0 then go to dvc1;
	     sh = (dh/dm) * dl + rm;
	     if hn < 1.25e0 then go to step;
	     dl = sl - sm;
	     dh = sh - sm;
	     dm = dh + dl;
	     if dm = 0.0e0 then go to dvc2;
	     th = (dh/dm) * dl + sm;
	     uh = th;
	     if uh = um then go to done;
	     if uh = ul then go to done;
step:	     hn = hn + 0.5e0;
	     rl = rm;
	     sl = sm;
	     tl = tm;
	     ul = um;
	     rm = rh;
	     sm = sh;
	     tm = th;
	     um = uh;
	     go to loop;
dvc1:	     sh = rh;
dvc2:	     th = sh;
done:	     th = 1.128379167095512574e0 * th;
	     if n < 0.0e0 then th = 2.e0 * exp (f*f) - th;
	end large;
	r = round (th, 28);
	return (r);
     end exerfc_;




		    exp_.pl1                        10/03/83  1719.9rew 10/03/83  1005.0        6597



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

/* calls to round_ remove by A. Downing 07/16/73 */

exp_: procedure (number) returns (float binary (27));

/*    compute the exponential of a single-precision floating-point number    */
declare number float binary (27),
	exp builtin;

	return(exp(number));
	end;




		    iexp_.pl1                       10/03/83  1719.9rew 10/03/83  1005.1        5967



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

iexp_: procedure (base, exponent) returns (fixed binary (35));

/*     compute integer base ** integer exponent     */
declare (base, exponent) fixed binary (35);

	return(base**exponent);
	end;




		    log_.pl1                        10/03/83  1719.9rew 10/03/83  1005.1       20043



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

/* modified by A. Downing to remove calls to round_ */

log_: procedure (number) returns (float binary (27));

/*    compute the logarithm or hyperbolic arctangent of a single-precision floating-point number     */
declare (number, r) float binary (27),
        (a, f, n, p) float binary (63),
         i fixed binary (17);
dcl  code_ ext entry (fixed bin),
     (abs, addr, log, log10, log2, round, sign) builtin,
     1 word aligned based,
     2 exponent fixed bin (7) unal;
	return(log(number));			/* Natural log of value. */

long:	i = addr (f) -> exponent;
	addr (f) -> exponent = 0;
if f < 0.7071067811865475244e0 then lower: do;
	     a = 0.5946035575013605334e0;
	     n = 0.75e0;
end lower; else upper: do;
	     a = 0.840896415253714543e0;
	     n = 0.25e0;
	end upper;
	f = (f - a) / (f + a);
	n = i - n;
short:	if abs (f) < 0.7450580597e-8 then go to finis;
	a = f*f;
	f = (((0.2857142857142857143e0 * a + 0.4e0) * a + 0.6666666666666666667e0) * a + 2.e0) * f;
finis:	f = (0.6931471805599453094e0 * n + f) * p;
	r = round (f, 28);
	return (r);

log2_: entry (number) returns (float binary (27));
	return(log2(number));			/* log(2) of value. */

log10_: entry (number) returns (float binary (27));
	return(log10(number));			/* log(10) of value. */

lone_: entry (number) returns (float binary (27));
	return(log(number+1.0e0));			/* Natural log of x+1. */

atanh_: entry (number) returns (float binary (27));
	p = 0.5e0;
	f = number;
	a = abs (f);
	if a < 0.1e0	then do;
	     n = 0.0e0;
	     go to short;
	     end;
	a = a - 1.e0;
	if a >= 0.0e0 then err2: do;
	     call code_ (sign (a) + 42);
	     if a = 0.0e0 then f = 170141182.e30 * f; else f = 0.0e0;
	     return (f);
	end err2;
	f = (1.e0 + f) / (1.e0 - f);
	go to long;
     end log_;




		    random_.alm                     10/03/83  1719.9rew 10/03/83  1005.1       83610



"  ******************************************************
"  *                                                    *
"  *                                                    *
"  * Copyright (c) 1972 by Massachusetts Institute of   *
"  * Technology and Honeywell Information Systems, Inc. *
"  *                                                    *
"  *                                                    *
"  ******************************************************

"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
" This procedure generates pseudo-random numbers using the 
" Tausworth method.  36 bits are used in the generation.
"
" There are multiple entry points.  For all entry points:
"	The first argument is a fixed binary input argument,
"     which is a non-zero integer.  This is an optional argument--
"     if not provided by caller, a value maintained in internal
"     static is used.  This value, from either source,
"     is the seed for the random number generator.  Its value is
"     modified so that upon return it has the value that should
"     be used as the seed for the next call.
"
" There are a set of entry points with two arguments which
" are used to generate a single random number.  For these:
"	The second argument is a floating point return argument
"     that returns the value of the random number generated.
"
" There are a set of entry points with three arguments which
" are used to generate a sequence of random numbers.  For these:
"	The second argument is an array of single precision
"     floating point numbers.  This array returns a sequence of
"     of random numbers, beginning at the base of the array.
"	The third argument is a fixed binary(17) input
"     argument the specifies the size of the array.
"
"	Coded 1 January 1970 by Roger R. Schell
""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""

	name	random_

" Table of contents
	entry	set_seed
	entry	get_seed
	entry	random_
	entry	uniform
	entry	uniform_ant
	entry	uniform_seq
	entry	uniform_ant_seq
	entry	normal
	entry	normal_ant
	entry	normal_seq
	entry	normal_ant_seq
	entry	exponential
	entry	exponential_seq


	equ	shift,11
	equ	size,36
	equ	expon,0		exponent to convert integer to floating point


" " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " "
"
" CODING CONVENTIONS
"
"     XR0	used for return address for specific distribution subroutines
"     XR1	used for return address for generator primitive subroutine
"     XR2	general purpose register for distribution subroutines
"     XR3	usedto indicate: 1=> antithetic variable, 0=> usual
"     XR4	contains the address of the distribution subroutine for this call
"     XR5	index into return array for the next random number
"     XR6	count of the number of values to be generated after current one
"     XR7	general purpose register
"
"     A-reg distribution routine uses to return floating point value
"     Q-reg always has the seed used by primitive generator
"
"     BP	pointer to base of return arguments
"     AP	pointer to the seed
"
" " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " "




"
"	call random_$set_seed(seed);


set_seed:
	ldq	ap|2,*	qet new seed into Q-reg
	stq	lp|internal_seed	save as new value of seed
	tra	return	return to caller



"
"	call random_$get_seed(seed);


get_seed:
	ldq	lp|internal_seed	get current value of seed
	stq	ap|2,*	return value to caller
	tra	return	return to caller



"
"	call random_$uniform(seed,random_no);


random_:
uniform:
	eax4	uniform_	set XR4 to the address of uniform distribution
	tra	single	this entry generates a single random number



"
"	call random_$uniform_ant(seed,random_no);
"		This entry gives negatively correlated value.



uniform_ant:
	eax4	uniform_ant_	set up the proper distribution
	tra	single	this entry generates a single random number



"
"	call random_$uniform_seq(seed,array,array_size);
"		This entry gives an array of return values



uniform_seq:
	eax4	uniform_	we generate sequence from uniform distribution
	tra	sequence	this entry gives a sequence of numbers



"
"	call random_$uniform_ant_seq(seed,array,array_size);


uniform_ant_seq:
	eax4	uniform_ant_	negatively correlated generator
	tra	sequence	this entry gives a sequence of numbers



"
"	call random_$normal(seed,random_no);


normal:
	eax4	normal_	normal distribution
	eax3	0	not negatively correlated value
	tra	single	this entry gives a single number



"
"	call random_$normal_ant(seed,random_no);


normal_ant:
	eax4	normal_	normal distribution
	eax3	1	negatively correlated
	tra	single	this entry gives single number



"
"	call random_$normal_seq(seed,array,array_size);



normal_seq:
	eax4	normal_	normal distribution
	eax3	0	not negatively correlated value
	tra	sequence	this entry gives a sequence of numbers



"
"	call random_$normal_ant_seq(seed,array,array_size);


normal_ant_seq:
	eax4	normal_	normal distribution
	eax3	1	negatively correlated
	tra	sequence	this entry gives a sequence of numbers



"
"	call random_$exponential(seed,random_no);


exponential:
	eax4	exponential_	exponential distribution
	tra	single	this entry gives a single value



"
"	call random_$exponential_seq(seed,array,array_size);


exponential_seq:
	eax4	exponential_	exponential distribution
	tra	sequence	this entry gives a sequence of numbers



"!!!!!!!!!!--set up the number of values to be generated--!!!!!!!!!!


sequence:
	ldx7	ap|0	twice number of arguments in XR7
	lxl6	ap|0,7*	length of sequence to XR6
	eax7	-2,7	subtract two from XR7
	eax6	-1,6	decrement by one
	tpl	common	if positive value, use it
	tra	return	if zero or negative, return

single:
	ldx7	ap|0	twice number of arguments in XR7
	eax6	0	use sequence of length one

common:
	eppbp	ap|0,7*	set bp to point to first return value
	eax5	0	index into array is in XR5
	eaa	-2,7	upper A-reg has offset of seed in arglist
	ars	19	should be one or zero in A-reg
	xec	set_ap,al	set ap to point to the seed
	szn	ap|0	test for a seed of zero
	tnz	loop	if non-zero continue
	eax4	zero_arg	if zero, generate zero return values

loop:
	tsx0	0,4	go to appropriate generator
	fst	bp|0,5	store value returned by generator
	eax5	1,5	increment index into array
	eax6	-1,6	decrement count of remaining
	tpl	loop	if not done, loop again

return:
	short_return



set_ap:		"get pointer to seed--from caller or default
	eppap	lp|internal_seed	use internal value if not provided in call
	eppap	ap|2,*	seed is the first argument if provided




"$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
"$	This is the primitive that actually generates the
"$	random number in integer form from the seed.
"$
"$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$



generate:
	ldq	ap|0	load seed into Q-reg
	qrl	shift	shift right the seed
	ersq	ap|0	exclusive or to the seed
	ldq	ap|0	put same value in Q-reg
	qls	size-shift	shift the result left
	ersq	ap|0	exclusive or to previous result
	tra	0,1	return to the caller of primitive




"!!!!!!!!!!--zero argument generator--!!!!!!!!!!

zero_arg:		"used if input seed is zero
	fld	=0.,du	load a floating point zero
	tra	0,0	return




"!!!!!!!!!!--uniform generator--!!!!!!!!!!

uniform_:
	tsx1	generate	generate one random number
	lda	ap|0	load A-reg with integer value
	arl	1	make it a positive number
	lde	expon,du	convert to floating point
	fad	=0.,du	normalize
	tra	0,0	return


uniform_ant_:
	tsx1	generate	generate one random number
	lda	ap|0	load integer value into A-reg
	arl	1	make it a positive number
	lde	expon,du	convert to floating point
	fneg		"take negative value
	fad	=1.,du	normalize
	tra	0,0	return




"!!!!!!!!!!--exponential generator--!!!!!!!!!!

exponential_:
	eax7	-1	count number of 'runs' with XR7
outer:
	eax7	1,7	add one to count of 'runs'
	eax2	1	use as counter of 'run' length
			"initialize XR2 with a count of one
	tsx1	generate	go to primitive generator
	lda	ap|0	get seed in A-reg
	arl	1	make it a positive number
	lde	expon,du	convert to floating point
	fst	bp|0,5	store it temporarily in return value
inner:
	lda	ap|0	keep value in A-reg for comparison
	tsx1	generate	generate another value
	eax2	1,2	add one to count of 'run'length
	cmpa	ap|0	compare last number with new one
	trc	inner	if still a run down,loop again
	anx2	=1,du	check if 'run' has even length
	tnz	outer	if not even, get another run
	eaa	0,7	no of runs before even length
	lde	=17b25,du	convert to floating point
	fad	bp|0,5	add first random number to number of 'runs'
	tra 	0,0	return




"!!!!!!!!!!--normal distribution generator--!!!!!!!!!!

normal_:
	fld	=0.,du	load a zero
	eax2	12	use XR2 to count 12 times thru loop
n_loop:
	fst	bp|0,5	store the new sum
	tsx1	generate	generate the next random number
	lda	ap|0	load seed into A-reg
	arl	1	make it a positive number
	lde	expon,du	convert to floating point
	fad	bp|0,5	add random number to sum
	eax2	-1,2	decrement counter by one
	tnz	n_loop	accumulate twelve numbers
	fsb	=6.,du	give a mean of zero
	xec	n_norm,3	antithetic if appropriate
	tra	0,0	return

n_norm:
	nop	"o.k. as is if not antithetic
	fneg	"take negative for antithetic





"
"	INTERNAL STATIC DATA
"

	use	.lkstat.
	join	/link/.lkstat.

internal_seed:
	dec	4084114320	"initial internal seed for a new process

	use	main


	end




		    round_.alm                      10/03/83  1719.9rew 10/03/83  1005.1       11403



"  ******************************************************
"  *                                                    *
"  *                                                    *
"  * Copyright (c) 1972 by Massachusetts Institute of   *
"  * Technology and Honeywell Information Systems, Inc. *
"  *                                                    *
"  *                                                    *
"  ******************************************************

	name	round_
	entry	round_
round_:	
	fld	ap|2,*	get high-order fl.pt. number.
	tze	store	special case if zero.
	lrs	72	mantissa = 0 if +, -2**(-71) if -.
	adla	128,dl	add 2**(-28) to mantissa.
	dfad	ap|2,*	round d.p. fl.pt. number.
store:	fst	ap|4,*	store rounded s.p. number.
	short_return
	entry	expon_
expon_:	
	lda	ap|2,*	get high-order fl.pt. number.
	lrs	28	make exponent an integer.
	qrl	8	then isolate high order mantissa.
	stq	ap|2,*	store high-order mantissa.
	sta	ap|4,*	store integer exponent.
	short_return		number = mantissa * 2**exponent.
	entry	adexp_
adexp_:	
	lda	ap|2,*	get high-order fl.pt. number.
	lrs	28	make exponent an integer.
	adla	ap|4,*	add power to exponent.
	lrs	8	attach high-order mantissa.
	stq	ap|2,*	store high-order result:
	short_return		number * 2**power.
	end




		    sin_.pl1                        10/03/83  1719.9rew 10/03/83  1005.1        8955



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

/* modified 07/16/73 by A. Downing to use round builtin function */

sin_: procedure (number) returns (float binary (27));

/*    compute the sine or cosine of a single-precision floating-point number     */
declare number float binary (27),
	(cos, cosd, sin, sind) builtin;

	return(sin(number));

cos_: entry (number) returns (float binary (27));
	return(cos(number));

sind_: entry (number) returns (float binary (27));
	return(sind(number));

cosd_: entry (number) returns (float binary (27));
	return(cosd(number));
	end;




		    sinh_.pl1                       10/03/83  1719.9rew 10/03/83  1005.1       13356



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

/* modified by A. Downing 07/16/73 to use round builtin function */

sinh_: procedure (number) returns (float binary (27));

/*     compute the hyperbolic sine or cosine of a single-precision floating-point number     */
declare (number, n, p, r) float binary (27),
        (a, b) float binary (63);
dcl  code_ ext entry (fixed bin),
     (abs, exp, round) builtin;
	n = number;
	r = abs (n);
if r >= 0.67943378e0 then large: do;
	     p = -1.e0;
	     go to sinhs;
	end large;
	if r < 5.e-5 then go to negate;
	a = r;
	b = a*a;
	a = ((((0.2755731922398589065e-5 * b + 0.1984126984126984127e-3) * b + 0.8333333333333333333e-2) * b
	+ 0.1666666666666666667e0) * b + 1.e0) * a;
	go to finis;
cosh_: entry (number) returns (float binary (27));
	n, r = abs (number);
	p = 1.e0;
sinhs: if r > 88.028e0 then err: do;
	     call code_ (39);
	     r = 170141182.e30;
	     go to negate;
	end err;
	a = exp (r);
	a = (p/a + a) * 0.5e0;
finis:	r = round (a, 28);
negate:	if n < 0.0e0 then r = -r;
	return (r);
     end sinh_;




		    sqrt_.pl1                       10/03/83  1719.9rew 10/03/83  1005.1        6777



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

/* modified by A. Downing on 07/16/73 to remove calls to round_ */

sqrt_: procedure (number) returns (float binary (27));

/*     compute the square root of a single-precision floating-point number     */
declare number float binary (27),
	sqrt builtin;

	return(sqrt(number));
	end;




		    tan_.pl1                        10/03/83  1719.9rew 10/03/83  1005.1        6777



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

tan_: procedure (number) returns (float binary (27));

/*    compute the tangent of a single-precision floating-point number    */
declare number float binary (27),
	(tan, tand) builtin;

	return(tan(number));

tand_:	entry(number) returns(float bin(27));
	return(tand(number));
	end;




		    tanh_.pl1                       10/03/83  1719.9rew 10/03/83  1005.1       12537



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

/* modified by A. Downing on 07/16/73 to use round builtin */

tanh_: procedure (number) returns (float binary (27));

/*     compute the hyperbolic tangent of a single-precision floating-point number     */
declare (number, n, r) float binary (27),
        (f, p, q) float binary (63);
dcl  (abs, exp, round) builtin;
	n = number;
	f = abs (n);
if f >= 0.55e0 then if f >= 10.5e0 then f = 1.e0; else large: do;
	     r = f + f;
	     f = exp (r);
	     f = (f - 1.e0) / (f + 1.e0);
end large; else if f >= 5.e-5 then small: do;
	     p = f*f;
	     q = (((0.2480158730158730159e-4 * p + 0.1388888888888888889e-2) * p + 0.4166666666666666667e-1) * p
	     + 0.5e0) * p + 1.e0;
	     f = ((((0.2755731922398589065e-5 * p + 0.1984126984126984127e-3) * p
	     + 0.8333333333333333333e-2)* p + 0.1666666666666666667e0) * p + 1.e0) * f / q;
	end small;
	r = round (f, 28);
	if n < 0.0e0 then r = -r;
	return (r);
     end tanh_;




		    xp22_.pl1                       10/03/83  1719.9rew 10/03/83  1005.1        6444



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

xp22_:	procedure (base, exponent) returns (float binary (27));

declare base float binary (27),
         f float binary (63),
         exponent fixed binary (71);

dcl  round builtin;

	f = base;
	return (round(f ** exponent, 28));
     end xp22_;




		    xp2_.pl1                        10/03/83  1719.9rew 10/03/83  1005.2        6948



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

/* modified 07/16/73 by A. Downing to use round builtin */

xp2_: procedure (base, exponent) returns (float binary (27));

/*    compute single-precision floating-point base ** integer exponent     */
declare base float binary (27),
         exponent fixed binary (17);

	return(base ** exponent);
	end;




		    xp3_.pl1                        10/03/83  1719.9rew 10/03/83  1005.2        7056



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

/* modified by A. Downing to use round builtin function 07/16/73 */

xp3_: procedure (base, exponent) returns (float binary (27));

/*    compute base ** exponent for single-precision floating-point numbers    */
declare base float binary (27),
         exponent float binary (27);

	return(base ** exponent);
	end;




		    fort_index_.pl1                 10/03/83  1719.9rew 10/03/83  1005.2        6669



/* ******************************************************
   *                                                    *
   *                                                    *
   * Copyright (c) 1972 by Massachusetts Institute of   *
   * Technology and Honeywell Information Systems, Inc. *
   *                                                    *
   *                                                    *
   ****************************************************** */

fort_index_:
	procedure (str1, str2) returns (fixed binary (35));

	/* This procedure is the external equivalent of the Fortran
	   intrinsic function INDEX. */

dcl	(str1, str2) character (*) parameter;
dcl	index builtin;

	return (index (str1, str2));

	end fort_index_;




		    fort_math_builtins_.fortran     10/03/83  1719.9rew 10/03/83  1005.3       22248



c ***********************************************************
c *                                                         *
c * Copyright, (C) Honeywell Information Systems Inc., 1982 *
c *                                                         *
c ***********************************************************
      %global card, ansi77;
      function abs_ (x)
c     perform external equivalent of in-line abs ()
      real x
      abs_ = abs (x)
      return
      end
      function iabs_ (i)
      integer i
      iabs_ = iabs (i)
      return
      end
      double precision function dabs_ (d)
      double precision d
      dabs_ = dabs (d)
      return
      end
      function dim_ (a, b)
      real a, b
      dim_ = dim (a, b)
      return
      end
      function idim_ (i, j)
      integer i, j
      idim_ = idim (i, j)
      return
      end
      double precision function ddim_ (x, y)
      double precision x, y
      ddim_ = ddim (x, y)
      return
      end
      function sign_ (x, y)
      real x, y
      sign_ = sign (x, y)
      return
      end
      function isign_ (i, j)
      integer i, j
      isign_ = isign (i, j)
      return
      end
      double precision function dsign_ (x, y)
      double precision x, y
      dsign_ = dsign (x, y)
      return
      end
      function aint_ (x)
      real x
      aint_ = aint (x)
      return
      end
      function aimag_ (c)
      complex c
      aimag_ = aimag (c)
      return
      end
      complex function conjg_ (c)
      complex c
      conjg_ = conjg (c)
      return
      end
      function len_ (chars)
      character* (*) chars
      len_ = len (chars)
      return
      end
      double precision function dint_ (d)
      double precision d
      dint_ = dint (d)
      return
      end
      function anint_ (x)
      real x
      anint_ = anint (x)
      return
      end
      double precision function dnint_ (d)
      double precision d
      dnint_ = dnint (d)
      return
      end
      function nint_ (x)
      real x
      nint_ = nint (x)
      return
      end
      function idnint_ (d)
      double precision d
      idnint_ = idnint (d)
      return
      end
      double precision function dprod_ (x, y)
      real x, y
      dprod_ = dprod (x, y)
      return
      end
      function mod_ (i, j)
      integer i, j
      mod_ = mod (i, j)
      return
      end
      function amod_ (x, y)
      real x, y
      amod_ = amod (x, y)
      return
      end



		    bull_copyright_notice.txt       08/30/05  1008.4r   08/30/05  1007.3    00020025

                                          -----------------------------------------------------------


Historical Background

This edition of the Multics software materials and documentation is provided and donated
to Massachusetts Institute of Technology by Group Bull including Bull HN Information Systems Inc. 
as a contribution to computer science knowledge.  
This donation is made also to give evidence of the common contributions of Massachusetts Institute of Technology,
Bell Laboratories, General Electric, Honeywell Information Systems Inc., Honeywell Bull Inc., Groupe Bull
and Bull HN Information Systems Inc. to the development of this operating system. 
Multics development was initiated by Massachusetts Institute of Technology Project MAC (1963-1970),
renamed the MIT Laboratory for Computer Science and Artificial Intelligence in the mid 1970s, under the leadership
of Professor Fernando Jose Corbato.Users consider that Multics provided the best software architecture for 
managing computer hardware properly and for executing programs. Many subsequent operating systems
incorporated Multics principles.
Multics was distributed in 1975 to 2000 by Group Bull in Europe , and in the U.S. by Bull HN Information Systems Inc., 
as successor in interest by change in name only to Honeywell Bull Inc. and Honeywell Information Systems Inc. .

                                          -----------------------------------------------------------

Permission to use, copy, modify, and distribute these programs and their documentation for any purpose and without
fee is hereby granted,provided that the below copyright notice and historical background appear in all copies
and that both the copyright notice and historical background and this permission notice appear in supporting
documentation, and that the names of MIT, HIS, Bull or Bull HN not be used in advertising or publicity pertaining
to distribution of the programs without specific prior written permission.
    Copyright 1972 by Massachusetts Institute of Technology and Honeywell Information Systems Inc.
    Copyright 2006 by Bull HN Information Systems Inc.
    Copyright 2006 by Bull SAS
    All Rights Reserved
