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 y88.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