                    



                       PDP-1 COMPUTER
             ELECTRICAL ENGINEERING DEPARTMENT
           MASSACHUSETTS INSTITUTE OF TECHNOLOGY
                  CAMBRIDGE, MASSACHUSETTS
                           02139

















                           PDP-50

               RELOCATABLE SUBROUTINE LIBRARY














                      October 17, 1971






















                
>>76<<                                

                          Preface

   This memo describes the programs available in the library
of relocatable subroutines.  These subroutines reside on the
FORTRAN  tape or  on public  DECtape 3  under the  file name
"fortlib" (See Appendix II for  the current location of  all
routines  described in this library).   Those on the FORTRAN
tape will be available to the loader automatically. Those on
DECtape  3 are not  available to the  loader unless the user
specifies to the loader that a private library tape is to be
used.

   The  programs described  herein are  divided into several
groups according to their function.  This is done to make it
easier  to  find  the  description  of  a  desired  program.

   The most frequent use of the library will be in  conjunc-
tion  with FORTRAN programs, however, it  may be used by any
program which is in the  relocatable format required by  the
loader.   In all cases, the programs  in the library use the
calling conventions required by FORTRAN (See Appendix I  for
a  description  of  the FORTRAN  subroutine  linkage conven-
tions), except  for those  in Section  4, which  may not  be
explicitly called from FORTRAN programs.

   New  subroutines will from  time to time  be added to the
library.  Anyone who writes a subroutine of general interest
is  invited  to  make it  available  in the  library  and to
document it in this memo.


                  Use of Library Routines

   FORTRAN does not  do automatic  argument mode  conversion
for calls to library routines, and does not know the mode of
the result  of library  functions.  Hence  all arguments  to
library routines must be in the mode required by the library
routine.  Library functions which  do not correspond to  the
naming  conventions of FORTRAN must be declared using any of
the type declaration statements  available in FORTRAN.   All
logical,  doubleprecision, and complex functions must there-
fore be declared.  Some real and integer functions may  also
need to be decclared.


             Conventions used in this Document

   The  modes of arguments and functions which correspond to
the usual FORTRAN convention (names beginning with I, J,  K,
L,  M, and N are integers_, all  others are real) will not be
expicitly mentioned  in this  document.   The modes  of  any
arguments  or function names which do not correspond to this
convention will be given explicitly.











                
                                 
                   Input-Output Routines


1.1       Typewriter

   All of the  typewriter output routines  type nothing  but
the  information  specified  in  the  description,  i.e., no
carriage returns, tabs, or  unmentioned spaces are  printed.
Input  to all input routines, with  the exception of tyi, is
terminated by either a tab or a carriage return.

subroutine fip(x)
          Read a  real  number and  store  it in  "x".   The
          number  may  include  a minus  (-)  sign  and/or a
          decimal point.   Exponential notation  is not  al-
          lowed.

subroutine fop(x)
          Print the real number stored in "x".

subroutine cip(c)
          Read  a complex number  and store it  in "c".  The
          format is  that  of  fip,  except  that  two  real
          numbers  separated by a  comma are read.  Example_.
          3.14159,-2.71828

subroutine cop(c)
          Print the complex number "c".

subroutine tin(j)
          Read a  decimal integer  into "j".   A minus  sign
          indicates a negative number.

subroutine dpt(j)
          Print the integer "j" in decimal.

subroutine opt(j)
          Print the integer "j" in octal.

subroutine ret
          Print a carriage return.

subroutine tab
          Print a tab.

subroutine spc
          Print a space.

subroutine txx(a)
          Print  the text  string stored  in the  array "a".
          The text  must be  packed  3 characters  per  word
          (left justified). The end of the text is indicated
          by the character centerdot (_ ).  This may be  used
          to  print  a  message  using  the  following call_.
          call txx('text  string  ended  with  centerdot_  ')










                
y                                
subroutine tyi(j)
          Read  one character  and store  it right justified
          in "j".

subroutine tyo(j)
          Type the character  stored in the  right third  of
          the integer variable "j".

subroutine printi(w,i)
          Type  the integer stored in "i" right justified in
          a field  of  "w" positions.   Leading  zeroes  are
          suppressed  and  the  field  is  filled  out  with
          spaces.  Parameter "w" must be a positive integer.
          If  the number cannot fit into the field, fill out
          the field with 'x''s.

subroutine printf(w,d,x)
          Type the  real number  "x"  right justified  in  a
          field  of  "w" positions  with  "d" digits  to the
          right of the decimal point.  Exponential (e or  E)
          notation  is  not used.   Leading zeroes  are sup-
          pressed and the field  is filled out with  spaces.
          Parameters  "w" and "d" must be positive integers.
          If the number cannot fit into the field, fill  out
          the field with 'x''s.

subroutine printe(w,d,p,x)
          Type  the  real number  "x"  right justified  in a
          field of  "w" positions  with  "d" digits  to  the
          right of the decimal point.  The number is printed
          using exponential  notation  (e or  E),  with  "p"
          significant  digits  to   he left  of  the decimal
          point. Leading zeroes are suppressed and the field
          is filled out with spaces.  Parameters "w" and "d"
          are positive integers.  Parameter "p" is an  inte-
          ger.  If "p" is negative, an appropriate number of
          leading zeroes will be printed to the right of the
          decimal  point.  If the number cannot fit into the
          field, fill out the field with 'x''s.


























                
p                                
1.2       Calcomp Plotter

   The conventions  for  naming the  x  and y  axis  on  the
plotter are as follows_.

          x  increases as  the paper moves  forward over the
          plotter drum.
          y  increases  as  the  pen  moves  the  the  left.

subroutine getccm
          Assign the calcomp plotter and use the present pen
          position as coordinate (0,0).

subroutine relccm
          Dismiss the calcomp plotter.  The plotter must  be
          dismissed  when  a program  is finished  using it.

common block pxp
          Common block pxp is 3 words long and contains  the
          current pen coordinates (jx,jy,jz).  "jx" and "jy"
          are  the  x  and  y  position  of  the  pen,  "jz"
          indicates  whether the pen is up  or down (jz = 32
          if pen is up, 16 if  pen is down).  "jx" and  "jy"
          are set to (0,0) by getccm. (jx,jy,jz) are updated
          whenever a  command  to  move the  pen  is  given,
          however,  when the pen is in the "up" position, it
          will not be moved until it is lowered or relccm is
          called.

subroutine penup
          Raise  the plotter pen.  If the pen is already up,
          do nothing.

subroutine pendwn
          Lower the  plotter pen.   If  the pen  is  already
          down, do nothing.

subroutine plot(ix,iy)
          Move  the  pen from  the  current pen  position to
          (ix,iy).  The  coordinates  (ix,iy)  are  integers
          corresponding  to 0.005  inch steps.   Pen motions
          made when  the pen  is  raised are  not  performed
          until the pen is lowered.

subroutine drawto(x,y)
          Move  the  pen from  the  current pen  position to
          (x,y).  The  coordinates  (x,y) are  real  numbers
          corresponding  to 0.005 inch  steps.  In all other
          respects, drawto is similar to plot.
















                

                                subroutine cpl(ih,io,k)
          Draw one character on the plotter.  Parameter "ih"
          is  the character  height in units  of 0.005 inch.
          Parameter "io"  is  the character  orientation  in
          units  of pi/4 radians from the x axis ("io" is 0,
          2, 4, or 6). "k" is an integer variable containing
          the  character to  be plotted  in its  rightmost 6
          bits.

subroutine pxx(ih,io,a)
          Plot a character string.  Parameters "ih" and "io"
          are  as in  subroutine cpl.   Parameter "a"  is an
          array containing  the text  to be  plotted  packed
          three  characters per  word (right  justified) and
          terminated with a centerdot (_ ).  This  subroutine
          may be used as follows to plot a message_.
          call pxx(i,j,'text ended with centerdot_ ')

subroutine ploti(ih,io,w,i)
          Plot the integer "i" right justified in a field of
          width "w".  Leading zeroes are suppressed and  the
          field  is filled out with spaces.  Parameters "ih"
          and "io" are the same as the corresponding parame-
          ters  in subroutine  cpl.  "w" must  be a positive
          integer.  If  the number  will  not fit  into  the
          specified  field,  the  field is  filled  out with
          'x''s.

subroutine plotf(ih,io,w,d,x)
          Plot the  real number  "x"  right justified  in  a
          field of width "w" with "d" digits to the right of
          the decimal point.  Exponential (e or E)  notation
          is not used. Leading zeroes are suppressed and the
          field is filled out  with spaces.  Parameters  "w"
          and "d" must be positive integers. Parameters "ih"
          and "io" are the same as the corresponding parame-
          ters in subroutine cpl. If the number will not fit
          into the specified field, the field is filled  out
          with 'x''s.

subroutine plote(ih,io,w,d,p,x)
          Plot  the  real number  "x"  right justified  in a
          field  of  width  "w".   "d"  decimal  places  are
          printed.   "p" is the number of significant digits
          to be plotted  to the left  of the decimal  point.
          Leading  zeroes  are suppressed  and the  field is
          filled out with  spaces.  Parameters  "w" and  "d"
          must  be  positive  integers.  "p"  may  be either
          positive or  negative.   If  it  is  negative,  an
          appropriate  number  of  leading  zeroes  will  be
          plotted to the right of the decimal point. Parame-
          ters "ih" and "io" are the same as the correspond-
          ing parameters in subroutine  cpl.  If the  number
          will  not fit into the  specified field, the field
          is filled out with 'x''s.











                
e                                1.3       Display

   All of the following  display generating subroutines  use
the  same scaling conventions.   (0,0) is the  center of the
display screen.  The height and width of the screen are both
1023 units. Hence the edges of the screen are at coordinates
-511 and 511.

subroutine lds(ix,iy,jx,jy)
          Draw a line from (ix,iy) to (jx,jy).

subroutine dpy(ix,iy)
          Display a point at (ix,iy).

subroutine dlist(ix,iy,n)
          Display "n" points.  Parameters "ix" and "iy"  are
          arrays whose length is given by parameter "n". For
          each element j of these arrays, display the  point
          (ix(j),iy(j)).

subroutine dlist1(iy,n)
          Display  a list of ordinates.  The elements of the
          array "iy"  are  displayed  as a  set  of  equally
          spaced  ordinates.   The parameter  "n"  gives the
          length of the array "iy".

subroutine dlist2(ixy,n)
          Display a compact  list of  cooordinates.  "n"  is
          the length of the list.  "ixy" is an integer array
          containing coordinates packed 2  per word.  The  x
          coordinate divided by 2 is stored in the left half
          of a word, the y coordinate divided by 2 is stored
          in  the right half of a word.  This format is very
          compact and  quite useful  when space  is  scarce.
































                
>>52<<                                
1.4       Paper Tape Reader and Punch

subroutine getptr(lab)
          Assign  the paper  tape reader.   If it  cannot be
          assigned, transfer to statement label "lab".

subroutine relptr
          Deassign the paper  tape reader.  This  subroutine
          should always be called when a program is finished
          using the reader.

subroutine rpa(n,eot)
          Read one character from paper tape into the  least
          significant  8 bits of the parameter "n".  When an
          attempt is made to read  past the physical end  of
          the  tape, the program will transfer to the state-
          ment label "eot".

subroutine getptp(lab)
          Assign the  paper tape  punch.   If it  cannot  be
          assigned, transfer to statement label "lab".

subrouting relptp
          Deassign  the paper  tape punch.   This subroutine
          should always be called when a program is finished
          using the punch.

subroutine ppa(n)
          Punch the character contained in the least signifi-
          cant of parameter "n".



































                
b                                1.5       Miscellaneous

subroutine lat(j)
          Read the test word switches into the "j".

function switch(j)
          Returns  the  logical value  true is  sense switch
          "j" is  up, false  if sense  switch "j"  is  down.

function sltest(j)
          Returns  the logical value of sense light (program
          flag) "j" (1.<j.<6).  The  logical value is true  if
          the sense light is on.  Sense light 3 is turned on
          if the  light pen  "sees" light  during a  display
          instruction.

subroutine slight(j)
          If  1.<j.<6, turn  off sense light  "j".  If 9.<j.<14,
          turn on sense light "j"-8. Sense light 3 is turned
          on  if the light pen "sees" light during a display
          instruction.

subroutine rbt(j,k)
          Read the buttons  and switches  of button  console
          "j"  into "k".  The left half of "k" is set by the
          buttons and the right  half of "k"  is set by  the
          switches of the button console. A depressed button
          or switch  places a  "1"  in the  appropriate  bit
          position.   Currently,  button  consoles  0  and 1
          exist.

subroutine ckn(j,k)
          Read knob "j" into "k".   The value read will  lie
          between  0  and  255.  Knobs  0  through  7 exist.
































                
>>60<<                                
1.6       DECtape File System

   These subroutines make  it possible for  programs to  use
the DECtape file system (See PDP-42, MICROTAPE FILE SYSTEM).
Since these subroutines work through the normal file  system
program,  there is little danger  of destoying files through
program errors.

subroutine dtopen(filename,drive,fileno,rw)
          Open a DECtape file.  "filename" is an array  con-
          taining  the name of the file packed three charac-
          ters per word (left justified) and terminated with
          a centerdot (_ ).  "drive" is an integer containing
          the number  of  the  tape  drive  upon  which  the
          DECtape  containing the file is mounted.  "fileno"
          is an integer which will be used in all references
          to  the file after it has been opened.  "rw" is an
          integer which indicates whether the file is  being
          opened  for  input  or  for  output.   If  "rw" is
          positive, the file  may only be  read_, if "rw"  is
          negative,  the file may only be written.  One file
          may  simultaneously   be  open   under   different
          "fileno"'s.   Since all of these subroutines cause
          physical actions  to take  place immediately,  the
          effects  of writing and  reading several "fileno"s
          which all correspond  to the same  file on a  tape
          are  calculable and  deterministic.  If  a file is
          opened under a "fileno"  which already is in  use,
          the previous use is pushed, that is, only the most
          recently opened file is available.  When that file
          is  closed,  the previous  use  is popped  and the
          original file  opened  as "fileno"  may  again  be
          manipulated.

subroutine dtclose(fileno)
          Close a DECtape file.  The file opened as "fileno"
          (an integer) closed.  If this is the last file  on
          the  tape that was  open, the tape  is rewound and
          the directory rewritten if necessary. After a file
          has  been closed,  no more use  may be  made of it
          until it is reopened.  A file must be closed  when
          it  is  no  longer  required  to  ensure  that the
          directory is rewritten properly.   A file must  be
          open before it may be closed.

subroutine rtape(fileno,buffer,eof)
          Read  from a DECtape file.  One block of 256 words
          is read from file  "fileno" (an integer) into  the
          array "buffer". Blocks of a file are read one at a
          time starting with the first block. Reading starts
          from  the first block when  a file is opened.  Any
          attempt to read past the end of a file will result
          in an error return to the statement labeled "eof".












                
/                                subroutine wtape(fileno,buffer)
          Write on a DECtape file.   One block of 256  words
          is  written from the array  "buffer" onto the file
          "fileno" (an  integer).   This  subroutine  always
          appends to the end of a file. If the file does not
          exist on the tape (no  name in the tape  directory
          corresponds to the "filename" under which the file
          was opened), the file is created with length 0 and
          then  this  file  is  appended  to.   Thus,  it is
          impossible to write  over any file  existing on  a
          tape using only the wtape subroutine.

subroutine dltape(fileno)
          Delete  a  DECtape  file.   The  file  "fileno" is
          deleted.  A file must be open for output in  order
          for it to be deleted. If one desires to rewrite an
          existing file, it must first be deleted before  it
          can  be rewritten.  If it is not first deleted, it
          will be appended to, not rewritten.















































                
b                                                  Mathematical Subroutines


2.1       Trigonometric Functions

   All 72 combinations  of real,  doubleprecision, and  com-
plex_,  normal and inverse_, sine, cosine, tangent, cotangent,
secant, and cosecant_, circular and hyperbolic are  available
in this library.

            FUNCTION                    MODE

                         REAL      DOUBLEPRECISION   COMPLEX

            CIRCULAR     sin       dsin              csin
                         cos       dcos              ccos
                         tan       dtan              ctan
                         cot       dcot              ccot
                         sec       dsec              csec
                         csc       dcsc              ccsc

            HYPERBOLIC   sinh      dsinh             csinh
                         cosh      dcosh             ccosh
                         tanh      dtanh             ctanh
                         coth      dcoth             ccoth
                         sech      dsech             csech
                         csch      dcsch             ccsch

            INVERSE      asin      dasin             casin
            CIRCULAR     acos      dacos             cacos
                         atan      datan             catan
                         acot      dacot             cacot
                         asec      dasec             casec
                         acsc      dacsc             cacsc

            INVERSE      asinh     dasinh            casinh
            HYPERBOLIC   acosh     dacosh            cacosh
                         atanh     datanh            catanh
                         acoth     dacoth            cacoth
                         asech     dasech            casech
                         acsch     dacsch            cacsch

          Each  of these functions  takes one argument.  The
          argument mode  is  the same  as  the mode  of  the
          result in each case.  All inverse functions return
          the principal value.

          In addition, the functions atan2 (real) and datan2
          (doubleprecision)  are  available.  Each  of these
          functions takes two  arguments.  The  mode of  the
          arguments is the same as the mode of the result in
          each case.  atan2 is real and datan2 is doublepre-
          cision.  atan2(y,x) returns the value of atan(y/x)
          with the  sign corresponding  to the  quadrant  of
          (x,y).











                
n                                
2.2       Additional Elementary Functions

   The  library contains functions  to compute square roots,
cube roots,  logs, and  exponentials of  real, complex,  and
doubleprecision  quantities.  Each of  these functions takes
one argument. The mode of the result is the same as the mode
of the argument in all cases.

          FUNCTION   VALUE MODE         DESCRIPTION

           sqrt      real               square root
           dsqrt     doubleprecision    square root
           csqrt     complex            square root

           cbrt      real               cube root
           dcbrt     doubleprecision    cube root

           alog      real               natural (base e) log
           alog10    real               common (base 10) log
           dlog      doubleprecision    natural (base e) log
           dlog10    doubleprecision    common (base 10) log
           clog      complex            natural (base e) log

           exp       real               exponential funcon
           dexp      doubleprecision    exponential function
           cexp      complex            exponential function


2.3       Miscellaneous Elementary Functions

function carg(c)
          The function carg takes a complex argument "c" and
          returns a real value which  is the phase (arg)  of
          the complex number "c".

function cabs(c)
          The function cabs takes a complex argument "c" and
          returns a real value which  is the modulus of  the
          complex number "c".

function dmod(dx,dy)
          The  function dmod  takes 2  doubleprecision argu-
          ments and  returns a  doubleprecision value.   The
          value returned is (dx-integerpart(dx/dy)xdy).





















                
u                                
2.4       Random Number Generators

subroutine randu(ix,iy,yfl)
          Subroutine  randu generates a pseudo-random number
          using a multiplicative congruential method.   "ix"
          is  the seed.  "iy" is set to the new value of the
          seed by the subroutine.  "yfl" is a random  number
          uniformly  distributed on  the interval (0.0,1.0].
          "ix" and "iy" are  integers, "yfl" is real.   "ix"
          and "iy" may be the same variable.

subroutine ran(x)
          Subroutine ran generates a pseudo-random mumber on
          the interval (0.0,1.0].  The random number  gener-
          ated is stored in the real variable "x".  There is
          no provision for resetting the generator to regen-
          erate  a  particular sequence  of  random numbers.
          Every sequence of calls to ran generates the  same
          sequence of random numbers.

subroutine gauss(ix,s,am,y)
          This subroutine generates a gaussian pseudo-random
          number through calls to randu.  The integer  vari-
          able  "ix" is the seed and is updated by each call
          to gauss. The real variable "s" is the variance of
          the  random numbers  produced.  The  real variable
          "am" is  the  mean  of the  random  numbers.   The
          gaussian  random  number  is  stored  in  the real
          variable "y".

subroutine sysrng(n)
          This subroutine generates a somewhat random number
          and stores it in the variable "n". This routine is
          useful for providing random seeds for other random
          number generators.  It is fairly random so long as
          it is not called very frequently or regularly.  It
          generates a value in the range (0,4095).




























                
>>53<<                                
2.5       Matrix Operations

subroutine invrt(n,a,b,s,r)
subroutine dinvrt(n,a,b,s,r)
subroutine cinvrt(n,a,b,s,r)
          These subroutines may be used to invert real, dou-
          bleprecision, and complex matrices,  respectively.

                 Parameter   Description
                 n           Matrix dimension (integer)
                 a           Matrix  to be  inverted  (real,
                               doubleprecision,  or  complex
                               nxn array)
                 b           To become inverse matrix (real,
                               doubleprecision,  or  complex
                               nxn array)
                 s           Scratch space (real, doublepre-
                               cision, or complex nxn array)
                 r           Scratch  space  (logical or in-
                               teger n array)


          The  mode of the parameters "a", "b", and "s" cor-
          respond to the mode of the matrix to be  inverted.
          These  subroutines  compute  the  inverse  of  the
          matrix "a" and store it in the matrix "b".  Arrays
          "s" and "r" are used in the computation. Parameter
          "n" contains  the size  of the  matrix to  be  in-
          verted.

subroutine detrm(n,a,s,r)
subroutine ddetrm(n,a,s,r)
subroutine cdetrm(n,a,s,r)
          These  functions may be used to compute the deter-
          minant of real,  doubleprecision, and complex  ma-
          trices,  respectively.   The  mode  of  the  value
          returned by the function is  the same as the  mode
          of  the matrix  whose determinant  is being evalu-
          ated.  Parameters "n", "a", "s", and "r" have  the
          same  meaning as  the corresponding  parameters in
          invrt, dinvrt, and cinvrt.


2.6       Sorting Programs

subroutine sortac(x,n)
          Subroutine sortac sorts the real array "x" of  "n"
          elements into ascending order.

















                
i                                
2.7       Elliptic Integrals

function elpk(m)
          Real elliptic integral K(m)
function elpkc(m)
          Real complementary elliptic integral K'(m)
function elpe(m)
          Real elliptic integral E(m)
function elpec(m)
          Real complementary elliptic integral E'(m)

          K'(m) = K(1.0-m)
          E'(m) = E(1.0-m)

          m  is parameter, or square  of modulus (m is real)
                    0.0 .< m < 1.0 for elpk
                    0.0 < m .< 1.0 for elpkc
                    0.0  .<  m  .<  1.0  for  elpe  and  elpec

function delpk(m)
          Doubleprecision elliptic integral K(m)
function delpkc(m)
          Doubleprecision  complementary  elliptic  integral
          K'(m)

          K'(m) = K(1.0-m)

          m is parameter, or square of modulus (m is double-
          precision)
                    m < 1.0 for delpk
                    m > 0.0 for delpkc

function celpk(m)
          Complex elliptic integral K(m)
function celpkc(m)
          Complex complementary elliptic integral K'(m)

          K'(m) = K(1.0-m)

          m  is  parameter, or  square of  modulus (complex)

                    For celpkc, | m|  > 0.0, and the  function
                    is  analytic  for  -pi  <  arg(m)  < pi.

                    For celpk, | 1.0-m|  > 0.0, and the  func-
                    tion  is analytic for -pi < arg(1.0-m) <
                    pi.


















                
(                                
2.8       Jacobi Elliptic Functions

function sn(x,m)
function cn(x,m)
function dn(x,m)
function sd(x,m)
function cd(x,m)
function nd(x,m)
function sc(x,m)
function dc(x,m)
function nc(x,m)
function cs(x,m)
function ds(x,m)
function ns(x,m)
          Each of these subroutines  computes one of the  12
          Jacobi elliptic functions.  All of these functions
          return real values. Parameter "x" is the argument.
          Real  parameter "m" is the parameter, or square of
          the modulus.  0.0 .< m .< 1.0


2.9       Elliptic Modular Functions

function emodf(tau)
function emodg(tau)
function emodh(tau)
          These three  functions  are used  to  compute  the
          elliptic modular functions.  Parameter "tau" (com-
          plex) is the period ratio. Each of these functions
          returns a complex value.

                       Im(tau) >~0.03

          In terms of theta functions_.

          emodf(tau) = (theta2(0,tau)/theta3(0,tau))^4
          emodg(tau) = (theta4(0,tau)/theta3(0,tau))^4
          emodh(tau) = -(theta2(0,tau)/theta4(0,tau))^4


2.10 Gamma Function

gamma(x)
          The  function gamma returns a  real value which is
          the gamma function of the parameter "x".

                 gamma(x) = factorial(x-1)


















                
p                                
2.11      Bessel Functions

function dj0(d)
          Doubleprecision Bessel function of the first  kind
          and  0th order.  Parameter "d" is doubleprecision.

function dy0(d)
          Doubleprecision Bessel function of the second kind
          and  0th order.  Parameter "d" is doubleprecision.

function dj1(d)
          Doubleprecision Bessel function of the first  kind
          and  1st order.  Parameter "d" is doubleprecision.

function dy1(d)
          Doubleprecision Bessel function of the second kind
          and  1st order.  Parameter "d" is doubleprecision.

function djn(n,d,da)
          Doubleprecision Bessel function of the first  kind
          and "n"th order. "n" must be at least 2. Parameter
          "n" is  an  integer, "d"  is  the  doubleprecision
          argument  of the  Bessel function,  and "da"  is a
          doubleprecision array of  "n"+1 words  for use  by
          djn.


2.12      Error Function

function erf(x)
          Real error function.
function derf(d)
          Doubleprecision error function.
function erfc(x)
          Real complementary error function.
function derfc(d)
          Doubleprecision complementary error function.

          These  functions  compute the  error  function and
          the complemented error function. In all cases, the
          mode  of the argument  is the same  as the mode of
          the result.























                
k                                2.13      Hermite Polynomials

function herm(x,n)
          Real Hermite polynomial.
function dherm(x,n)
          Doubleprecision Hermite polynomial.

          These functions compute the hermite polynomial  of
          order  "n" of the  argument "x".  The  mode of the
          argument "x" must be the  same as the mode of  the
          result.  Integer parameter "n" must be at least 0.


2.14      Fast Fourier Transform

subroutine fft(real,imag,exp,size)
subroutine ffti(real,imag,exp,size)
          These subroutines  compute  the  discrete  fourier
          transform of a set of complex data. Subroutine fft
          computes the forward transform_, ffti computes  the
          inverse  transform.  These subroutines work on the
          data in place, that is,  the result of the  trans-
          form replaces the data being transformed.

          These  subroutines use a  block-floating point im-
          plementation of the Cooley-Tukey algorithm.  As  a
          result,  neither the  input nor the  output of the
          transform is in any form that can be directly used
          numerically  by  FORTRAN.  However,  conversion is
          quite simple.

          The parameter "size" is  a positive integer  which
          must be a power of 2 and no greater than 1024. The
          parameters "real" and "imag" are integer arrays of
          length "size". Parameter "exp" is an integer which
          is an  exponent (base  2) which  applies to  every
          element  in  "real"  and "imag".   This  is block-
          floating point format.  The  input to fft or  ffti
          need  not be  normalized.  However,  the output is
          always normalized by making the elements in "real"
          and  "imag"  as  large  as  possible  and suitably
          modifying "exp".  The  reslt of  the transform  is
          returned in "real", "imag", and "exp".























                
1                                
2.15      Numerical Integration

subroutine ktmbeg(array,deriv,n,err,lim,init,max,scratch)
subroutine ktm
subroutine ktmret(lab)
          These  subroutines are different entry points to a
          package that will integrate systems of first order
          differential equations using a Kutta-Merson scheme
          with automatic step size adjustment.

          The package is  initialized by a  call to  ktmbeg.
          The  significance of the parameters is as follows_.

          array -- real array  of length "n" containing  the
                    variables to be integrated and the inde-
                    pendent  variable.   "array(1)"  is  the
                    independent variable.

          deriv  -- real array of  length "n" containing the
                    derivatives of the dependent  variables.
                    "deriv(1)"  must  be set  to  1.0 before
                    ktmbeg is called.

          n -- integer containing the number of variables in
                    "array" and "deriv". "n" is one plus the
                    number of  equations  being  integrated.

          err -- real variable containing the maximum allow-
                    able error in any  step of the  integra-
                    tion.   If it  is exceeded,  the program
                    will do  the  step  over  again  with  a
                    smaller step size.

          lim  -- real variable containing  the limit of in-
                    tegration for the independent  variable.
                    Integration  starts at the initial value
                    of the independent variable and proceeds
                    until  the independent  variable has in-
                    creased by "lim".

          init -- real variable containing the initial  step
                    size.   The program will adjust the step
                    size as necessary.

          max -- real variable  containing the maximum  step
                    size.  The program will never adjust the
                    step size in such a way that it  exceeds
                    "max".

          scratch  -- real  array of  length "n"x4  used for
                    temporary storage by the program.














                
9                                          A call  to ktm  will begin  the integration.   The
          integrator  will periodically  (about 4  times per
          step) ask for the derivatives of the variables  in
          "array". The first time it will do so by returning
          from ktm.  The calling  program must then  compute
          the  derivatives of the  current values of "array"
          and store then in "deriv".  Then it must return to
          the  integrator  through  a call  to  ktmret.  The
          argument "lab" of ktmret  is a statement label  to
          which  the integrator will  transfer when it wants
          more values of the derivatives.  When it does  so,
          the  program must  then recompute  the derivatives
          and store  them in  "deriv" and  then call  ktmret
          again.  When ktm makes a normal exit, the integra-
          tion is complete, and the results are in  "array".
          It  is  then  possible  to  call  ktm  again,  and
          integrate further without another call to  ktmbeg.

                          EXAMPLE

                dimension y(2),d(2),scr(8)
                call ktmbeg(y,d,2,1.e-7,.1,.01,1.,scr)
          c initial step size = 0.01, stop every 0.1
                d(1)=1.0
                y(1)=0.0
          c to begin the integration at 0.0
                y(2)=0.0
          c to give it the constant of integration
          10    call ktm
          20    d(2)=1.0/(1.0+y(1)xx2)
                call ktmret(20s)
          c dy/dx = 1/(1+x^2)
                write(0,30) y(1),y(2)
          30    format(2f10.5)
          c this makes a table of arctangents in steps of 0.1
                goto 10
                end





























                
                                              Pointer and Free Storage Package


3.1       Free Storage Management

subroutine getfre(ptr,size,core)
          This  subroutine will  obtain a  block of storage.
          All arguments are integers.  "size" is the  length
          of the block being requested.  If "core" is not 0,
          the block  will lie  entirely  in 1  core  module.
          Otherwise, the block may cross core module bounda-
          ries. The address of the block is stored in "ptr".

subroutine relfre(ptr,size)
          Calling  this  subroutine will  free the  block of
          memory whose address and size are contained in the
          integer parameters "ptr" and "size".


3.2       Pointer manipulation

function iaddr(n)
          This function returns the address of variable "n".
          "n" may be any  data type.  If it  is a data  type
          that  requires several  words of  storage, the ad-
          dress of  the  first (smallest  address)  word  is
          returned.

function iget(ptr)
function get(ptr)
function dget(ptr)
function cget(ptr)
          These  functions return the contents of the memory
          location whose address is in the integer parameter
          "ptr".  Functions iget, get, dget, and cget return
          integer, real, doubleprecision,  and complex  val-
          ues,  respectively.   For  those  data  types that
          occupy several words of storage, "ptr" is taken as
          the  address of the first (smallest address) word.

subroutine iput(x,ptr)
subroutine put(x,ptr)
subroutine dput(x,ptr)
subroutine cput(x,ptr)
          The subroutines store the  contents of the  varia-
          ble "x" in the memory location whose address is in
          the integer  parameter "ptr".   Subroutines  iput,
          put,  dput, and cput  store integer, real, double-
          precision, and complex values, respectively.   For
          those  data  types  that occupy  several  words of
          storage, "ptr"  is taken  as  the address  of  the
          first  (smallest address) word  to be stored into.














                
d                                            Arithmetic and Input-Output Routines


   This section  describes the  arithmetic routines  in  the
library.   These are  not called  with the  standard calling
sequence.  They  may not  be  explicitly called  by  FORTRAN
programs,  however, the compiler generates  calls to them to
do real, doubleprecision, and complex arithmetic.

   The calling sequence for these subroutines is as follows.
One argument may be in an accumulator.  This is the hardware
AC for integer and  logical data, the  AC-IO for real  data,
the doubleprecision accumulator (symbol .dac) for doublepre-
cision data, and the  complex accumulator (symbol .cac)  for
complex data.  Another argument (if it exists) is pointed to
by the word  following the jdp  instruction which calls  the
subroutine.  Results are returned in the appropriate accumu-
lator.

   In this  section  of  this document,  a  subroutine  name
written  in the  form 'function  foo(n)' indicates  that one
argument  to  subroutine  foo  is  indicated  by  a  pointer
following the 'jdp foo' instruction. Another argument may be
in an accumulator.


4.1       Real Arithmetic

   A real number is stored in 2 consecutive words of memory.
A  real quantity is  always addressed by  the first (smaller
address) word.  Bit 0 of the  first word is the sign of  the
number. When this bit is 0, the number is positive. Bits 1-9
are the exponent, which is stored in excess 256 form.   Bits
10-17  of the first word and the entire second word are a 26
bit fraction.   Real  numbers  are normalized  so  that  the
fraction  is less than one and greater than or equal to one-
half.  To  negate a  real number,  complement all  36  bits.

function .fad(x)
          Add  "x" to the real accumulator and store the sum
          in the real accumulator.

function .fmp(x)
          Multiply the real accumulator by "x" and store the
          product in the real accumulator.

function .fdv(x)
          Divide  the number in the  real accumulator by "x"
          and place the  quotient in  the real  accumulator.

function .fnm
          Normalize the real accumulator.














                
7                                4.2 Doubleprecision Arithmetic

   A doubleprecision number is stored in 5 consecutive words
of memory.  The last word (highest address) is the  exponent
(base  2), stored as an integer.   The first 4 words are the
fraction  stored  in  sign-magnitude  form,  with  the  most
significant  word last.   A doubleprecision  quantity is ad-
dressed using the address of  the first (lowest address)  of
the five words.

function .dld(d)
          Load  the  doubleprecision  accumulator  from "d".

function .dst(d)
          Store the  doubleprecision accumulator  into  "d".

function .dad(d)
          Add  "d"  to the  doubleprecision  accumulator and
          store the sum in the doubleprecision  accumulator.

function .dmp(d)
          Multiply  the  doubleprecision accumulator  by "d"
          and  store  the  product  in  the  doubleprecision
          accumulator.

function .ddv(d)
          Divide  the doubleprecision accumulator by "d" and
          store the quotient in the doubleprecision  accumu-
          lator.

function .dnm
          Normalize the doubleprecision accumulator.

function .dne
          Negate the contents of the doubleprecision accumu-
          lator.

common block .dac
          Doubleprecision accumulator.

function .dsp
          Skip iff the doubleprecision accumulator  contains
          a positive number.

function .dsn
          Skip  iff the doubleprecision accumulator contains
          a negative number.

function .dsi(d)
          The number in  the doubleprecision accumulator  is
          negated  if its  sign is not  the same  as that of
          "d".














                
_                                
4.3       Complex Arithmetic

   A complex number is  a pair of  real numbers.  The  first
real  number  is the  real part  of  the complex  number.  A
complex number is addressed by the first (lowest address) of
the four words.

function .cld(c)
          Load the complex accumulator from "c".

function .cst(c)
          Store the complex accumulator into "c".

function .cad(c)
          Add  "c" to the complex  accumulator and store the
          sum in the complex accumulator.

function .cmp(c)
          Multiply the complex accumulator by "c" and  store
          the product in the complex accumulator.

function .cdv(c)
          Divide  the complex  accumulator by  "c" and store
          the quotient in the complex accumulator.

function .cne
          Negate the complex accumulator.

common block .cac
          Complex accumulator.


4.4       Exponentiation Functions

function .iti(n)
          The integer  in  the AC  is  raised to  the  "n"th
          power.

function .rti(n)
          The  real  number in  the AC-IO  is raised  to the
          "n"th power.

function .dti(n)
          The number in  the doubleprecision accumulator  is
          raised to the "n"th power.

function .cti(n)
          The  number in  the complex  accumulator is raised
          to the "n"th power.
















                
7                                
4.5       Conversion Routines

   All conversion  routines  take their  argument  from  the
appropriate  accumulator  and  leave  their  result  in  the
appropriate accumulator.

function .flo
          Convert from integer to real.

function .fix
          Convert from real to integer.

function .sng
          Convert from doubleprecision to real.

function .dbl
          Convert from real to doubleprecision.

function .dnt
          Convert from doubleprecision to integer.

function .dlo
          Convert from integer to doubleprecision.










































                
                                
                         Appendix I


                 FORTRAN Subroutine Linkage

      The FORTRAN calling conventions are as follows_.

1.  FORTRAN  programs  run  in  extended  address  mode  and
      deferred address mode (eem and dam).

2.   All  subroutines  preserve the  contents  of  the index
      register.

3.  All subroutines must be  capable of running in any  core
      module  that the loader may put them into.  This means
      that 15-bit addresses must  be used in all  references
      to other subroutines and to arguments.  Also, it makes
      it necessary  to  be  careful  when  using  the  index
      register.

4. All subroutines and functions are called in the following
      manner_.

              jdp i (15-bit. subroutine. address
              15. bit. argument. address         
              15. bit. argument. address         
                              
              _ _ _                              
                              
              15. bit. argument. address         

      A subroutine has no way of knowing the number of argu-
      ments with which it was called.

5.  Functions return results in the appropriate accumulator.

       Data Type          Accumulator

       Integer            AC
       Logical            AC
       Real               AC-IO
       Doubleprecision    doubleprecision accumulator
       Complex            complex accumulator

      The AC  and AC-IO  are  the hardware  accumulator  and
      input-output registers of the PDP-1.  The doublepreci-
      sion accumulator  is a  set  of locations  inside  the
      doubleprecision  arithmetic package  referenced by the
      symbolic name ".dac". The complex accumulator is a set
      of  locations  inside the  complex  arithmetic package
      referenced by the symbolic name "cac".














                
u                                                        Appendix II


Location of subroutines in the library

   All of  the subroutines  in  the subroutine  library  are
stored  on  the  FORTRAN  tape  with  the  exception  of the
following which are available on public microtape 3.

Section     Subroutine

1.1         printi
            printf
            printe
1.2         ploti
            plotf
            plote
2.5         detrm
            ddetrm
            cdetrm
2.7         delpk
            delpkc
            celpk
            celpkc
2.8         sn
            cn
            dn
            sd
            cd
            nd
            sc
            dc
            nc
            cs
            ds
            ns
2.11        dj0
            dy0
            dj1
            dy1
            djn
2.12        erf
            derf
            erfc
            derfc
2.13        herm
            dherm
2.14        fft
            ffti

















                
y                                                                                                                                      h
