/* hclib.c ****************************************************************************** DIRECTORY accept get variable values from STDIN, as if fortran ACCEPT aloadb loads buffer with binary array (as leading low-byte strings) alpho compare two strings for alphabetical order amod returns remainder from division of two doubles aunload unload buffer to list of variables including arrays bclose | =close file bcreate| =create a new file for buffered block output bopen |bio =open file for buffered block output bread | =read buffer from file in blocks bwrite | =write buffer to file in blocks btof convert buffer to float cacept get character variable values from STDIN (accept for char) *callc =macro routine calls c routine from fortran (not in hclib.obj) call_emt call emt375 center_text center a string cget get character from STDIN without waiting for \n date =return the date datest return system date as a "dd-mmm-yy" string define_file| =open or create binary direct access file dread |dio =read a record from a binary direct access file dwrite | =write a record to a binary direct access file defname insert default filename and/or filename extension delay =delay given number of seconds dtols convert double to leading low-byte string errext type error message and exit ftols convert float to leading low-byte string ftor50 convert ascii string filename to RADIX-50 getln get line of input from a file getword get one word from a file greet type greeting message at beginning of program *istsx =returns 1 if TSX+, 0 if rt11 or TSX version < 5 *limits =returns limits of available system memory loadbf loads buffer with binary data (as leading low-byte strings) lstod convert leading low-byte string to double lstof convert leading low-byte string to float month string of month name nint nearest integer option was option specified on call to rt11 csi password =get password, exit if incorrect pause type message and pause until "RETURN" pflush empty printer output buffer ("REWIND 6") putln put line of input to a file putword put word to a file rename =rename a file *rng =random number generator rtcsi =call the rt11 csi *maccsi =issue .csispc emt and pop options off the stack,used by rtcsi r50tof convert RADIX-50 filename to ascii string shell shell sort *swap_byte =swap the bytes of an integer stopcc =disable (and restore) control-C substi substitute one substring for another substring tics =elapsed time in tics since the last call time =return time of day timest return system time as a "hh-mm-ss" string *ttflush =flush rt11 input ring buffer ungetc unget a character unload unload buffer to list of variables unrad convert RADIX-50 integer to three ascii characters transcendental functions dabs absolute value of double darctan(x) arctangent x dcos(x) cosine x subcos used by dsin and dcos dexp(x) e to the x power (antiloge) d10exp(x) 10 to the x power (antilog10) dloge(x) natural logarithm (base e) x dlog10(x) logarithm (base 10) x dpower(y,x) y to the x power droot(x,y) x root y dsin(x) sine x subsin used by dsin and dcos dsqrt(x) square root x dtento(x,i) x times 10 to the i power (whitesmiths) fabs absolute value of float iabs absolute value of integer polynm calculate polynomial, used by transcendental functions power(x, i) x to the i power, used by transcendental functions * these routines are macros at least in part and require editing of *.mac code. they are contained in hclib2.C and hclib2.mac. = these routines are rt11 dependent and thus non-portable installation instructions (for starting from scratch): hclib sources are gathered in 2 files: hclib1.c which contains this directory and source code for all routines and hclib2.c which contains source code for istsx, limits, maccsi, rng, and swap_byte. these routines in hclib2.c are also in hclib1.c (inactivated as comments) just to have everything in one place. hclib2.mac must be edited prior to assembly. run c hclib2 !contains limits, maccsi, rng, istsx, and swap_byte, !stop before assembly edit hclib2.mac !makes changes as per comments with source code macro hclib2 hclib1.c must then be copied into individual *.c files names with the first 6 characters of the routine name. three files contain more than one subroutine. DIO has define_file, dread, dwrite, and the external definitions they use. BIO has bcreate, bclose, bopen, bread, and bwrite. DSINCO has dsin, dcos, sinsub, and cossub. use hclib.bat to compile all of these modules. use @hclib.bld to build them into library hclib.obj. installation instructions for a new routine (not macro): copy *.c source code module to this disk. add include statements to module update directory in hclib1.c. add source code to hclib1.c (positioned in alphabetical order). add module name to hclib.bld. run c * (thru macro, not link) library/insert hclib * *****************************************************************************/ #include #include /****************************************************************************/ int accept(s, t, add1) /* accepts data from STDIN in manner similar to FORTRAN IV ACCEPT statement. s is string typed as a prompt. t is a string indicating how to convert the data. it only uses i, f, and d. any number and combination of variables can be used. returns the number of variables accepted. if too many are entered they are ignored, if too few 0 is entered. exponents are ok with f or d. */ char *s; /* prompt string */ char *t; /* "format" string */ unsigned add1; /* first argument, an address */ { #define IBUFSIZE 90 register int nc, *p; /* p is pointer to next address */ int *pi; float *pf; double *pd; /* pointers to output variables */ char buff[IBUFSIZE], *pbuff; /* input buffer and pointer to it */ int i, nchars, n_variables = 0; BYTES btod(), btoi(), inbuf(), notbuf(), getlin(); double tempd; /* temporary double for float input */ errfmt(s); /* type prompt */ nchars = getlin(buff, IBUFSIZE); /* get input line */ while((i = inbuf(buff, nchars, ",")) < nchars) buff[i] = ' '; /* change ',' to ' ' */ while((i = notbuf(buff, nchars, " 0123456789.eE+-@")) < nchars) buff[i] = '@'; /* change bad chars to '@' */ nchars = squeeze(buff, nchars, '@'); /* remove all '@' */ pbuff = buff; /* pbuff points at bottom of buff */ p = &add1; /* point p at top of argument stack */ while(*t) { /* read the 'format' string */ switch(*t++) { case 'i' : /* integer */ if(nchars <= 0) { /* if no more text return 0 */ pi = *p++; *pi = 0; break; } nc = btoi(pbuff, nchars, *p++, 10); pbuff += nc; nchars -= nc; n_variables++; break; case 'f' : /* float */ pf = *p++; if(nchars <= 0) { *pf = 0.; break; } nc = btod(pbuff, nchars, &tempd); *pf = (float) tempd; pbuff += nc; nchars -= nc; n_variables++; break; case 'd' : /* double */ if(nchars <= 0) { pd = *p++; *pd = 0.; break; } nc = btod(pbuff, nchars, *p++); pbuff += nc; nchars -= nc; n_variables++; break; default : break; } } return(n_variables); } /****************************************************************************/ int *aloadb(s, topadd, wcnt, starta) /* loads wcnt 'integer' values from an array into the buffer starting at s. the buffer goes to a (top address + 1) of topadd. starta is the address of the first word of the array to load. char arrays can be done but must end on an even boundary. it returns the address for the next addition to buffer. */ register int *s; int *topadd; int wcnt; unsigned starta; { register int *p; /* pointer to next word to load */ char *itols(); p = starta; /* point p to start of array */ for( ; wcnt > 0; wcnt--) { if(s < topadd) itols(s++, *p++); else { errfmt("aloadb - attempt to over fill buffer\n"); return; } } return(s); } /****************************************************************************/ int alpho(s, t) /* compares two strings s and t as to alphabetic order. returns < 0 if s < t, 0 if they are equal, or > 0 if s > t. */ register char *s, *t; { for( ; *s == *t; s++, t++) { if(*s == '\0') return(0); } return(*s - *t); } /****************************************************************************/ double amod(x, y) /* returns the remainder of x divided by y */ double x, y; { return(x - (y * (int) (x / y))); } /****************************************************************************/ int *aunload(buffer, topadd, s, add1) /* unloads the buffer starting at buffer with a (top address + 1) of topadd into any combination and number of int, float, or double. char can not be done. add1 is the address of the first variable. it returns the next address to use (in buffer) unless there is an error in which case 0 is returned. s is a string indicating how to do the unload. i indicates int, f float, d double. all other non-numeric characters are ignored. a number preceding the i, f, or d indicates an array of this many variables */ register int *buffer; int *topadd; char *s; unsigned add1; { register int *p; /* pointer to next address */ char t[6], *message; int i, narray; int *pi, lstoi(); unsigned btoi(); float *pf, lstof(); double *pd, lstod(); message = "aunload - attempt to unload past end of buffer\n"; p = &add1; /* point p at top of argument stack */ while(*s) { /* read the 'format' string */ for(i = 0; *s >= '0' && *s <= '9'; i++) { /* # in array */ t[i] = *s++; if(!*s) return(buffer); /* if eostr */ } if(i) btoi(t, i, &narray, 10); else narray = 1; if(*s == 'i') { /* integer */ pi = *p++; for( ; narray > 0; narray--) { if(buffer >= topadd) { errfmt(message); return(0); } *pi++ = lstoi(buffer++); } } else if(*s == 'f') { /* float */ pf = *p++; for( ; narray > 0; narray--) { if(buffer + 1 >= topadd) { errfmt(message); return(0); } *pf++ = lstof(buffer++); buffer += 1; } } else if(*s == 'd') { /* double */ pd = *p++; for( ; narray > 0; narray--) { if(buffer + 3 >= topadd) { errfmt(message); return(0); } *pd++ = lstod(buffer++); buffer += 3; } } s++; } return(buffer); } /****************************************************************************/ /* there are five block structured I/O subroutines: bclose, bcreate, bopen, bread, bwrite, and also a routine ftor50 to convert ascii filename string to rad50 chan is the channel number. EMTERR is byte 052. */ int bclose(chan) /* close file number chan. returns: 3 created duplicate files. emt call is to .close */ int chan; { int emt(); if(emt(0374, 03000 + chan) < 0) { errfmt("bclose - duplicate files\n"); return(EMTERR); /* 3 */ } return(-1); /* success! */ } /****************************************************************************/ int bcreate(chan, name, length) /* create a new file chan. length is blocks. returns the error: 0 channel already open, 1 no space if length given or full if not, 3 already exists. emt call is to .enter */ int chan, length; char *name; /* ASCII */ { int emt375(), name_rad50[4]; ftor50(name, name_rad50); if(emt375(01000 + chan, name_rad50, length) < 0) { if(!EMTERR)errfmt("bcreate - channel already open\n"); if(EMTERR == 1)errfmt("bcreate - not enough space\n"); if(EMTERR == 3) errfmt("bcreate - protected file already exists\n"); return(EMTERR); /* 0, 1, or 3 */ } return(-1); /* success! */ } /****************************************************************************/ int bopen(chan, name) /* opens file chan for buffered block output. returns the error: 0 channel already open, 1 channel not found on device. emt call is to .lookup */ int chan; char *name; /* ASCII */ { int emt375(), name_rad50[4]; ftor50(name, name_rad50); if(emt375(0400 + chan, name_rad50) < 0) { if(!EMTERR) errfmt("bopen - channel already open\n"); if(EMTERR == 1) errfmt("bopen - file not found\n"); return(EMTERR); /* 0 or 1 */ } return(-1); /* success! */ } /****************************************************************************/ int bread(chan, p, n, b) /* read #n 16-bit words to buffer p from block number #b of file chan. returns: 0 attempted read past end of file, 1 hardware error, 2 channel not open. emt call is to .readw */ int chan, *p, n, b; { int emt375(); if(emt375(04000 + chan, b, p, n, 0) < 0) { if(!EMTERR)errfmt("bread - attempted read past EOF\n"); if(EMTERR == 1) errfmt("bread - hardware error\n"); if(EMTERR == 2) errfmt("bread - channel not open\n"); return(EMTERR); /* 0, 1, or 2 */ } return(-1); /* success! */ } /****************************************************************************/ int bwrite(chan, p, n, b) /* write #n 16-bit words from buffer p to block number #b of file chan. returns: 0 attempted write past end of file, 1 hardware error, 2 channel not open. emt call is to .writw */ int chan, *p, n, b; { int emt375(); if(emt375(04400 + chan, b, p, n, 0) < 0) { if(!EMTERR) errfmt("bwrite - attempted write past EOF\n"); if(EMTERR == 1) errfmt("bwrite - hardware error\n"); if(EMTERR == 2)errfmt("bwrite - channel not open\n"); return(EMTERR); /* 0, 1, or 2 */ } return(-1); /* success! */ } /****************************************************************************/ unsigned btof(s, n, pfnum) /* converts the n char buffer starting at s into a float stored at address pfnum. the string is a text representation of the number. an exponent is allowed. leading whitespace is skipped and an optional sign and fraction are allowed. conversion stops at the end of n char or at the first unrecognizable char. it returns the number of char consumed. if it gets no digits, *pfnum is zero. see "the C programming language" page 69. */ char *s; /* starting address of buffer */ unsigned n; /* number of characters */ float *pfnum; /* address at which to store float */ { double dtento(), val, power; int power10, i, powsign = 1, sign = 1; for(i = 0; s[i] == ' ' || s[i] == '\n' || s[i] == '\t'; i++) if(i >= n) { /* skip whitespace */ *pfnum = 0.; return(0); } if(s[i] == '+' || s[i] == '-') /* deal with sign */ sign = (s[i++] == '+') ? 1 : -1; for(val = 0; s[i] >= '0' && s[i] <= '9'; i++) { if(i > n) { /* get digits before '.' */ *pfnum = sign * val; return(i); } val = 10 * val + s[i] - '0'; } if(s[i] == '.') i++; for(power = 1; s[i] >= '0' && s[i] <= '9'; i++) { if(i >= n) { /* get digits after '.' */ *pfnum = sign * val / power; return(i); } val = 10 * val + s[i] - '0'; power *= 10; } if(s[i] == 'e' || s[i] == 'E') { i++; if(s[i] == '+' || s[i] == '-') /* deal with exponent sign */ powsign = (s[i++] == '+') ? 1 : -1; } for(power10 = 0; s[i] >= '0' && s[i] <= '9'; i++) { if(i >= n) { /* get digits after 'e' */ *pfnum = (float) dtento((sign * val / power), powsign * power10); return(i); } power10 = 10 * power10 + s[i] - '0'; } *pfnum = (float) dtento((sign * val / power), powsign * power10); return(i); } /****************************************************************************/ int cacept(s, n, add1) /* accepts char variables in a manner similar to FORTRAN IV ACCEPT statement. s is a string to be typed as a prompt. n is the number of character variables. it returns the number of variables accepted. if too many are entered, they are ignored, if too few - 0 is entered. */ char *s; /* prompt string */ int n; /* number of variables */ unsigned add1; /* first argument, an address */ { #define IBUFSIZE 82 register *p; /* p is pointer to next address */ char *pc; /* pointer to output variables */ char buff[IBUFSIZE], *pbuff; /* input buffer and pointer to it */ int i, nchars, n_variables = 0; BYTES squeeze(), inbuf(), getlin(); errfmt(s); /* type prompt */ nchars = getlin(buff, IBUFSIZE); /* get input line */ while((i = inbuf(buff, nchars, " ,")) < nchars) buff[i] = '@'; /* change ' ,' to '@' */ nchars = squeeze(buff, nchars, '@'); /* remove all '@' */ pbuff = buff; /* pbuff points at bottom of buff */ p = &add1; /* point p at top of argument stack */ for(i = 0; i < n; i++) { /* read the 'format' string */ pc = *p++; if(nchars <= 0) { /* if no more text return 0 */ *pc = 0; break; } *pc = *pbuff++; nchars--; n_variables++; } return(n_variables); } /****************************************************************************/ /* callc this macro subroutine is not in hclib.obj. it must be edited and reassembled for each use. it can be used to allow a fortran program to call a C subroutine. fortran calls callc as "call callc(...)". any number of arguments can be used. since addresses are passed from fortran, double can be passed as well as int. the name of the c routine to call must be edited into callc before assembly. the c routine looks like "int csub(x,y,z)\nint *x;\ndouble *y;\n,int *z;\n" etc. return (ie fortran functions) doesnt work. the c routine can do no output to the terminal and very complicated c routines (calling others) will crash. fortran common can easily be had by passing the address of the first member and adding offsets. .dsabl reg .enabl gbl .psect c$text .globl callc callc: mov (%5),%4 ;number of arguments l1: mov %4,%3 ;2 * nargs is add %4,%3 ;the offset add %5,%3 ;added to %5 -> address of argument mov (%3),-(%6) ;push argument sob %4,l1 ;repeat if any more arguments jsr %7,csub ;edit in the name of the C routine to call mov (%5),%4 ;number of arguments l2: tst (%6)+ ;restore stack pointer sob %4,l2 ;repeat if more arguments rts %7 .end */ /****************************************************************************/ int call_emt(code, channel, arg1, arg2, arg3, arg4) /* issue a emt375 call */ unsigned code, channel, arg1, arg2, arg3, arg4; { return(emt375((code<<8) + channel, arg1, arg2, arg3, arg4)); } /****************************************************************************/ center_text(s, t, n, c) /* centers the string t within the first n characters of string s using fill character c before and after the centered text. n does not include the '\0'. if t is too long nothing is done. */ char *s, *t, c; int n; { unsigned lenstr(); int front, back; register int i, diff; if((diff = n - lenstr(t)) < 0) return; back = diff / 2; front = (diff % 2) + back; for(i = 0; i < front; i++) *s++ = c; for( ; *t != '\0'; ) *s++ = *t++; for(i = 0; i < back; i++) *s++ = c; *s = '\0'; } /****************************************************************************/ char cget() /* get character from STDIN without waiting for \n!! */ { char c; while(!read(STDIN, &c, 1)) ; return(c & BYTMASK); } /****************************************************************************/ int date(n) /* returns the date into a three int array, the first has the month, second the date, third the year. the argument is the address of int a[3]; the emt call is .date */ register int *n; { register int k, i; int emt(); for(k = 0; k < 3; k++) *n++ = 0; if((i = emt(0374, 05000)) <= 0) { if(!i) errfmt("date - no date entered on system\n"); if(i < 0) errfmt("date - failed\n"); return; } *--n = 1972 + (i & 037); *--n = (i >>= 5) & 037; *--n = (i >>= 5); } /****************************************************************************/ char *datest(s) /* return system date as a string "dd-mmm-yy" */ char *s; { char *p, *month(); int i, k = 0, n[3], date(); unsigned decode(); date(n); k = decode(s, 2, "%i", n[1]); s[k++] = '-'; p = month(n[0]); for(i = 0; i < 3; i++) s[k++] = p[i]; s[k++] = '-'; k += decode(&s[k], 2, "%i", n[2] - 1900); s[k] = '\0'; return(s); } /***************************************************************************/ /* the following extern declarations are required by define_file, dread, and dwrite, which follow */ extern int *inputb = NULL; /* block read buffer */ extern int outputb[BLOCKSIZE] = {0} /* block write buffer */ struct { int nrecs; /* number of records per file */ int words; /* number of 16-bit wrds per record */ } dfile[10] = {0}; /****************************************************************************/ int define_file(dfn, dfname, dnrecs, dwords, fnew) /* sets up a direct access binary disk file in the manner of the FORTRAN IV "DEFINE FILE" statement. dfn is file number (used as channel number), dfname is the file name, dnrecs is the number of records, dwords is the number of 16-bit words per record. fnew carries instructions regarding new files: if fnew is 0, no new file is created (ie return if unsuccessful open); if fnew is < 0, a new file must be created (ie delete an existing file); and if fnew is > 0 the file is created if it can not be opened. the default filename is dk:cfile*.dat where * is dfn. up to 10 files can be open at a time. allocation of space for the input buffer is based on the largest record size ever used in the program. note that a new file is only as large as the highest number block actually written to. returns -1 if successful, EMTERR if open or create fail, 2 if file number too high, or 0 if fnew is 0 and the file doesn't exist. use dread and dwrite with this routine */ int dfn; char *dfname; int dnrecs, dwords; int fnew; { int defname(), bcreate(), bopen(); int size, nblocks; static int bufsize = 0; char *alloc(), *free(); if(dfn > 9) { errfmt("define_file - file number > 9\n"); return(2); } defname(dfname, ".dat", dfn); if(fnew < 0) goto label1; if(bopen(dfn, dfname) != -1) { /* attempt to open */ if(EMTERR != 1) { if(!EMTERR) errfmt("define_file - channel in use\n"); else errfmt("define_file - unable to open file\n"); return(EMTERR); } if(fnew == 0) { /* if fnew=0, do not create */ errfmt("define_file - fnew = 0, unable to open file\n"); return(0); } label1: nblocks = (dnrecs * dwords) / BLOCKSIZE + 1; /*not,create it*/ if(bcreate(dfn, dfname, nblocks) != -1) { errfmt("define_file - unable to create file\n"); return(EMTERR); } } size = dwords + BLOCKSIZE; if(size > bufsize) { free(inputb, 0); inputb = alloc(size * 2, 0); bufsize = size; } dfile[dfn].nrecs = dnrecs; dfile[dfn].words = dwords; return(-1); } /****************************************************************************/ int dread(dfn, recnum, rwords, buffer) /* direct access read similar to a FORTRAN IV "READ" statement. dfn is the file number(channel number), recnum is the number of the record to read(starting at 0), rewords is the number of words of this record to read, and buffer is the integer buffer to fill. returns -1 if successful, the .readw error code if this fails, or a local error code (3 or 4). file must have been opened with define_file. */ int dfn, recnum, rwords; int *buffer; { int bread(); int bufwrd, wrdpos, block, i, j; if(dfn < 1 || dfn > 9) { errfmt("dread - file number out of range\n"); return(5); } if(recnum >= dfile[dfn].nrecs) { errfmt("dread - record index to high\n"); return(3); } if(rwords > dfile[dfn].words) { errfmt("dread - attempt to read more words than record\n"); return(4); } wrdpos = recnum * dfile[dfn].words; block = (wrdpos + 1) / BLOCKSIZE; /* starting block number */ wrdpos = wrdpos % BLOCKSIZE; bufwrd = wrdpos + rwords; /* how many to read */ if(bread(dfn, inputb, bufwrd, block) != -1) { errfmt("dread - unable to read\n"); return(EMTERR); } for(j = wrdpos, i = 0; i < rwords; j++, i++) buffer[i] = inputb[j]; return(-1); } /****************************************************************************/ int dwrite(dfn, recnum, rwords, buffer) /* direct access write similar to a FORTRAN IV "WRITE" statement. dfn is the file number(channel number), recnum is the number of the record to write(starting at 0), rewords is the number of words of this record to write, and buffer is the integer buffer with the values to write. returns -1 if successful, the .readw or .writew error code that was the failure, or a local error code (3 or 4). file must have been opened with define_file. */ int dfn, recnum, rwords; int *buffer; { int bread(), bwrite(); int nblocks, wrdpos, block, i, j, k; if(dfn < 1 || dfn > 9) { errfmt("dwrite - file number out of range\n"); return(5); } if(recnum >= dfile[dfn].nrecs) { errfmt("dwrite - record index too high\n"); return(3); } if(rwords > dfile[dfn].words) { errfmt("dwrite - attempt to write more words than record\n"); return(4); } wrdpos = recnum * dfile[dfn].words; /* word position to start */ block = (wrdpos + 1) / BLOCKSIZE; /* starting block number */ wrdpos = wrdpos % BLOCKSIZE; if(rwords % BLOCKSIZE == 0) /* number of blocks to write*/ nblocks = (wrdpos + rwords - 1) / BLOCKSIZE + 1; else nblocks = (wrdpos + rwords) / BLOCKSIZE + 1; for(k = 0; k < nblocks; k++) { if(k == 0 || k == nblocks - 1) { if(bread(dfn, outputb, BLOCKSIZE, block + k) != -1) { errfmt("dwrite - unable to read block\n"); return(EMTERR); } } if(k == 0) { j = wrdpos; i = 0; } else j = 0; for( ; j < BLOCKSIZE && i < rwords; j++, i++) outputb[j] = buffer[i]; if(bwrite(dfn, outputb, BLOCKSIZE, block + k) != -1) { errfmt("dwrite - unable to write\n"); return(EMTERR); } } return(-1); } /****************************************************************************/ int defname(file_name, file_ext, n) /* if a filename string is NULL add default file name based on n (ie cfile*.dat where * is n). if no file extension is present, it is added. if no device name is present "dk:" is added. returns 0 if no change, 1 if default name added, +2 if file extension added, +3 if device name added */ char file_name[], file_ext[]; int n; { char s[15], *cpystr(); int ret = 0, t; register int i, k; unsigned decode(); /* * change first \n to \0, convert lower to upper case */ for(i = 0; file_name[i] != '\0'; i++) { if(file_name[i] == '\n') { file_name[i] = '\0'; break; } } /* * default file name */ if(file_name[0] == '\0' || file_name[0] == '\n') { /* default name */ n %= 100; /* 0 <= n <= 99 */ cpystr(file_name, "DK:CFILE*.DAT", NULL); t = n > 9; decode(file_name + 8 - t, t + 1, "%i", n); return(1); } /* * is device name too long, if so shorten */ for(i = 0; file_name[i] != 0; i++) { if(file_name[i] == ':') { if(i > 3) { for(k = 3; file_name[i] != '\0'; i++, k++) file_name[k] = file_name[i]; file_name[k] = '\0'; } break; } } /* * is name too long, if so shorten */ for(i = 0; file_name[i] != 0; i++) { if(file_name[i] == ':') break; } if(file_name[i] != ':') i = 0; else i++; for(k = 0; file_name[k] != 0; k++) { if(file_name[k] == '.') break; } if(k - i > 6) { for(i += 6; file_name[k] != '\0'; i++, k++) file_name[i] = file_name[k]; file_name[i] = '\0'; } /* * is file extension too long, if so shorten */ for(i = 0; file_name[i] != 0; i++) { if(file_name[i] == '.') { if((lenstr(file_name) - i) > 4) { file_name[i + 4] = '\0'; } break; } } /* * add file extension */ for(i = 0; file_name[i] != '\0'; i++) { if(file_name[i] == '.') /* if '.' found, no add ext */ goto label1; } ret = 2; file_ext[4] = '\0'; /* extension only 4 char */ cpystr(file_name, file_name, file_ext, NULL); /* add file ext */ /* * add device name */ label1: /* for(i = 0; file_name[i] != '\0'; i++) { if(file_name[i] <= 'z' && file_name[i] >= 'a') file_name[i] -= 32; } */ /* * add device name */ for(i = 0; file_name[i] != '\0'; i++) { /* device name? */ if(file_name[i] == ':') return(ret); } cpystr(s, "DK:", file_name, NULL); /* add device name */ cpystr(file_name, s, NULL); return(ret + 3); } /****************************************************************************/ int delay(seconds) /* delay for seconds seconds */ unsigned seconds; { int time_1, time_2; int emt375(); LONG i; if(emt375(010400, &i) < 0) return; time_1 = i / 60; time_2 = time_1; while(seconds > time_2 - time_1) { if(emt375(010400, &i) < 0) return; time_2 = i / 60; } } /****************************************************************************/ char *dtols(s, d) /* the double version of itols. converts the value of a double into the 8 byte string starting at s. returns s (ie the address). this is faster than using union */ register int *s; double d; { register int *p; char *itols(); p = &d; itols(s++, *p++); itols(s++, *p++); itols(s++, *p++); itols(s, *p); return(s - 3); } /****************************************************************************/ int errext(s) /* prints error message and exits. */ char *s; { errfmt(s); exit(); } /****************************************************************************/ char *ftols(s, a) /* the float version of itols. converts the value of a float into the 4 byte string starting at s. returns s (ie the address). this is faster than using union */ register int *s; float a; { register int *p; char *itols(); p = &a; itols(s++, *p++); itols(s, *p); return(--s); } /****************************************************************************/ int ftor50(ascii, radix) /* converts an ascii string filename to RADIX-50 format in 4 integers pointed to by 'radix'. taken from c manual on rad50 routine */ char *ascii; int *radix; { int rad50(), n; unsigned lenstr(), scnstr(); n = scnstr(ascii, ':'); if(ascii[n]) { radix[0] = rad50(ascii, n); ascii += n + 1; } else radix[0] = 0; n = scnstr(ascii, '.'); radix[1] = rad50(ascii, n); radix[2] = rad50(ascii + 3, n - 3); radix[3] = (ascii[n] == '.') ? rad50(ascii + n + 1,\ lenstr(ascii + n + 1)) : 0; } /****************************************************************************/ int getln(fp, line, limit) /* get line (up to limit number of characters including \n\0) of input from the file with address fp. returns the number of characters (excluding the terminal '\0') or EOF. line is terminated with '\n' '\0'. '\0' is treated like any other char. */ FIO *fp; register char *line; int limit; { register int c, nc = 1; int getc(); while(nc < limit - 1 && (c = getc(fp)) != '\n') { if(c == EOF) return(EOF); *line++ = c; nc++; } *line++ = '\n'; *line = '\0'; return(nc); } /****************************************************************************/ int getword(fp, word, limit) /* gets one word (up to limit number of characters including \0) of input from the file with FIO address fp and places the word at address word. assumes that words are delimited by ascii <= ' '. returns the number of characters (excluding the terminal '\0') or EOF. the word is terminated with '\0' and thus always contains at least one character. */ FIO *fp; register char *word; register int limit; { int c; register int nc = 0; int getc(); while(nc < limit - 1 && (c = getc(fp)) > ' ') { *word++ = c; nc++; } if(c == EOF) return(EOF); *word = '\0'; return(nc); } /****************************************************************************/ int greet(message) /* types a greeting message for the begining of programs. */ char *message; { char *str1 = "\n=============================================\ ===================================\n\n"; errfmt(str1); errfmt(message); errfmt(str1); } /****************************************************************************/ /* int istsx() /* returns 1 if running under TSX+ v5 or later. returns 0 if running under RT11 or an older version of TSX+. .dsabl reg .enabl gbl .globl istsx .psect c$text,i,ro ; rmon = 54 sysgen = 372 ; istsx: ; mov %5,-(%6) mov %6,%5 ; clr %0 mov rmon,%1 tst sysgen(%1) bpl out ;rt11 ; mov #1,%0 ; out: jmp c$rets ; .title istsx .psect c$text,i,ro .even .psect c$data,d,rw .even .end */ /****************************************************************************/ /* int limits(offset, plow, ptop, top) return bottom address of program, top address of program, and top address available for use. offset is how much space to leave at the top for the heap. the *.mac file must be edited as indicated below. unsigned offset, *plow, *ptop, *top; { register x; x = offset; *plow = 1; *ptop = 2; *top = 3; } */ /* .mcall .settop ;hzb added limits: jsr %5,c$sav .settop #-2 ;hzb added get it all tst 4(%5) ;hzb added leave heap space? beq label1 ;hzb added if not sub 4(%5),%0 ;hzb added .settop %0 ;hzb added reset top for heap label1: mov limmt,@6(%5) ;hzb changed return addresses mov limmt+2,@10(%5) ;hzb changed mov %0,@12(%5) ;hzb changed jmp c$ret .psect c$data,d,rw .even limmt: .limit ;hzb added .even */ /****************************************************************************/ int *loadbf(s, topadd, wcnt, arg1) /* loads wcnt 'integer' values into the buffer starting at s. the buffer goes to a (top address + 1) of topadd. arg1 is the first value to load. any number and combination of integer and double can be loaded as long as wcnt is correct and there is room in the buffer. float can not be done because it is converted to double when passed as an argument. also char can not be done. it returns the address for the next addition to buffer. */ register int *s; char *topadd; int wcnt; int arg1; { register int *p; /* pointer to next word to load */ char *itols(); p = &arg1; /* point p at top of stack */ for( ; wcnt > 0; wcnt--) { if(s < topadd) itols(s++, *p++); else { errfmt("loadbf - attempt to over fill buffer\n"); return; } } return(s); } /****************************************************************************/ double lstod(s) /* the double version of lstoi. converts the 8 byte string starting at s to a double which is then returned. this is faster than using union */ register int *s; { double a; register int *p; int lstoi(); p = &a; *p++ = lstoi(s++); *p++ = lstoi(s++); *p++ = lstoi(s++); *p = lstoi(s); return(a); } /****************************************************************************/ float lstof(s) /* the float version of lstoi. converts the 4 byte string starting at s to a float which is then returned. this is faster than using union */ register int *s; { register int *p; int lstoi(); float a; p = &a; *p++ = lstoi(s++); *p = lstoi(s); return(a); } /****************************************************************************/ /* int *maccsi(defr50, filspec) this is a macro-11 routine to be called from a C program to invoke the rt11 csi and thus return the filenames and options to the calling program. it allows for 10 (decimal) words of option results which are popped off the stack after the EMT is issued. the basic routine was generated by the C compiler based on code at the end of this section; the *.mac file was then edited. this routine is used by rtcsi. int defr50[4]; int filspec[39]; { static int *x, *y, optbuf[10]; x = defr50; y = filspec; return(optbuf); } */ /* ; .psect c$data,d,rw .even $1: ; this is x .blkb 2 .even $3: ; this is y .blkb 2 .even $5: ; this is optbuf .blkb 24 .psect c$text,i,ro .mcall .csispc ; hzb added ; maccsi: mov %5,-(%6) mov %6,%5 mov 4(%5),$1 ; put address of defr50 in $1 mov 6(%5),$3 ; put address of filspec in $1 ; ; hzb added (must have reg enabled!!!) ; .csispc $3,$1 ; call csi in special mode mov #12,%4 ; set counter to move 10. words mov #$5,%3 ; point at optbuf l1: mov (sp)+,(%3)+ ; pop off stack, move to optbuf sob %4,l1 ; if < 4 done, do another ; ; end hzb added ; mov #$5,%0 ; return address of option block jmp c$ret .psect c$text .even .psect c$data .even */ /****************************************************************************/ char *month(n) /* returns a pointer to a string containing the text representation of the month name. n is the number (1-12) of the month. taken from "the C programming language" */ int n; { static char *name[] = { "illegal month", "January", "February", "March", "April", "May", "June", "July", "August", "September", "October", "November", "December" }; return((n < 1 || n > 12) ? name[0] : name[n]); } /****************************************************************************/ int nint(x) /* return the nearest integer */ double x; { if(x >= 0.) return((int) (x + 0.5)); else return((int) (x - 0.5)); } /****************************************************************************/ int option(os, c) /* was an option selected in a call to the rt11 csi? rtcsi() converts option information to a null terminated string. option scans a string to see if a character is in it (returns 1) or not (returns 0). */ register char *os, c; { while(*os != '\0') { if(*os++ == c) return(YES); } return(NO); } /****************************************************************************/ int password(s, n) /* prompt for and accept password from operator. s is the string containing the password which must be matched. no more than 10 chars of the string will be used. n is the number of times the operator is allowed to try it; afterward the program exit()'s. */ char *s; int n; { char sinput[11]; int i, cmpstr(); unsigned ijob, getlin(); /* set terminal to special mode */ ijob = JSW; /* remember job status word */ JSW = ijob | 010000; /* set bit 12 */ /* loop n times */ for( ; n > 0; n--) { /* get input */ errfmt("\nenter password "); getlin(sinput, 10); /* terminate input string */ for(i = 0; i < 10; i++) { if(sinput[i] == '\n') break; } sinput[i] = '\0'; /* exit if password incorrect */ if(cmpstr(sinput, s) == YES) { JSW = ijob; /* terminal to previous mode */ errfmt("\n"); return; } } JSW = ijob; /* terminal to previous mode */ exit(); } /****************************************************************************/ int pause(s) /* type string s to terminal, then await return. */ char *s; { errfmt("pause - %p", s); getfmt("%x"); } /****************************************************************************/ int pflush(plp) /* flush LP: buffer by flooding it with nulls. plp is the address of the FIO of the printer. */ FIO *plp; { int i; for(i = 0; i < plp->_nleft; i++) putc(plp, '\0'); } /****************************************************************************/ int putln(fp, line, limit) /* put line of output (up to limit number of characters) to the file with address fp. returns the number of characters (including the terminal '\n'). line is terminated with '\n'. treats '\0' like any other char. */ FIO *fp; register char *line; register int limit; { register int nc = 1; int putc(); while(nc < limit && *line != '\n') { putc(fp, *line++); nc++; } putc(fp, '\n'); return(nc); } /****************************************************************************/ int putword(fp, word, limit) /* put word of output (up to limit number of characters) to the file with address fp. assumes that words are delimited by ascii <= ' '. returns the number of characters (including the terminal '\n'). the output word is terminated with '\n' thus always containing at least one character. */ FIO *fp; register char *word; register int limit; { register int nc = 1; int putc(); while(nc < limit && *word > ' ') { putc(fp, *word++); nc++; } putc(fp, '\n'); return(nc); } /****************************************************************************/ int rename(old_name, new_name) /* rename a file. returns the error: -1 if it worked, 1 channel already open, 2 file not found, 3 protected file already exists. emt call is to .rename. */ char *old_name, *new_name; /* ASCII */ { int chan = 15; int bclose(), emt375(), name_rad50[8]; ftor50(old_name, name_rad50); ftor50(new_name, &name_rad50[4]); label1: if(emt375(02000 + chan, name_rad50) < 0) { if(EMTERR == 1) { errfmt("rename - channel already open\n"); chan += 1; goto label1; } if(EMTERR == 2) errfmt("rename - file not found\n"); if(EMTERR == 3) errfmt("rename - protected file already exists\n"); return(EMTERR); } if(bclose(chan) != -1) errfmt("rename - unable to close channel\n"); return(-1); /* success */ } /****************************************************************************/ /* float rng(addsd1, addsd2) pseudo-random number generator using the DEC RAN algorithm. see the mini-tasker, november 1979, Martin Ackroyd. the program below is converted to macro, %4 and %2 are changed to %2 and %3 respectively, and then the large block of text that follows must be added. sd1 and sd1 are int. contained in hclib2. register *addsd1; register *addsd2; { long ll; return; } */ /* the following code is present at start rng: jsr %5,c$sav add #-4,%6 mov 4(%5),%4 mov 6(%5),%2 the latter two lines are changed to mov 4(%5),%2 mov 6(%5),%3 the following is then added mov (%2),%0 mov (%3),%1 beq init ; asl %1 rol %0 ; add (%2),%0 add (%3),%1 adc %0 ; add (%3),%0 ; bpl next add #100000,%0 ; next: mov %0,(%2) mov %1,(%3) ; mov #201,%2 ; again: asl %1 rol %0 bcs done dec %2 br again ; done: clrb %1 bisb %0,%1 swab %1 clrb %0 bisb %2,%0 swab %0 ; ror %0 ror %1 ; mov %0,-12(%5) mov %1,-10(%5) ldf -12(%5),%0 ; br return init: mov #3,%1 inc %0 br next ; return: the last line was already there jmp c$ret */ /****************************************************************************/ int rtcsi(out_files, in_files, deftype, opts) /* call the rt11 csi in special mode (no lookup, fetch, or enter) to accept and process a command string from the terminal. returns -1 if successful, the error code if not. produces the names (as strings) of three output files and six input files. the names are blank if not entered. deftype is a 12 char string (4 times 3 char) with the default file types (input, output1, output2, output3). opts is a string which will contain the letter of each option selected (terminated be '\0' eg if options A and X were selected, opts would be "AX\0". although the rt11 csi returns more option information (ie which options go with which files and values of options) use of this information is not yet supported. */ char out_files[3][15], in_files[6][15], deftype[12], opts[10]; { int *poption, i; int filspec[39], defr50[4]; int r50tof(), rad50(), emt(), *maccsi(); /* * convert default type string to radix 50 */ for(i = 0; i < 4; i++) defr50[i] = rad50(&deftype[i * 3], 3); /* * maccsi is a macro11 subroutine which issues * an emt call to the rt11 csi in special mode */ poption = maccsi(defr50, filspec); /* * convert file names from radix 50 to strings */ for(i = 0; i < 3; i++) r50tof(&filspec[i * 5], &out_files[i][0]); for(i = 0; i < 6; i++) r50tof(&filspec[i * 4 + 15], &in_files[i][0]); /* * convert option information into a string */ for(i = 0; *(++poption) != NULL; ) { opts[i++] = *poption & 0177; /* use lower 7 bits */ if((*poption & ~077777) != 0) /* if bit 15 set */ poption++; /* burn the next word */ if(i == 9) break; } opts[i] = '\0'; for(i = 0; opts[i] != '\0'; i++) { if(opts[i] >= 'A' && opts[i] <= 'Z') opts[i] += 32; } /* * done */ return(-1); } /****************************************************************************/ int r50tof(radix, ascii) /* convert a 4 integer array containing a rt11 filename in radix 50 to an ascii string. */ int *radix; char *ascii; { int i; int unrad(); unsigned lenstr(), squeeze(); unrad(*radix++, &ascii[0]); ascii[3] = ':'; unrad(*radix++, &ascii[4]); unrad(*radix++, &ascii[7]); ascii[10] = '.'; unrad(*radix, &ascii[11]); ascii[14] = '\0'; /* * remove spaces */ squeeze(ascii, 15, ' '); /* * remove trailing '.', convert ":." to null string */ i = lenstr(ascii); if(ascii[--i] == '.') ascii[i] = '\0'; if(ascii[0] == ':' && ascii[1] == '\0') ascii[0] = '\0'; } /****************************************************************************/ int shell(n, pord, pexch, offset) /* shell sort. based on Software Tools and Knuth (algorithm for picking gap). pord is a routine with compares two items and pexch exchanges them. each is called with two arguments which are the indices (0 - n-1) of the items. offset is the index offset of where to start. n is the number of items. */ int n, offset; int (*pord)(), (*pexch)(); { int k; register i, j, gap; static int h[14] = { 0,1,4,13,40,121,364,1093,3280,9841,29524 88573,265720,797161 }; /* h[s+1] = 3*h[s]+1 */ for(k = 0; h[k] < n && k <= 14; k++) /* h[t+2] >= n, start h[t] */ ; if((k -= 2) < 1) k = 1; for( ; (gap = h[k]) > 0; k--) { for(i = gap + offset; i < n + offset; i++) { for(j = i - gap; j >= offset; j -= gap) { if((*pord) (j, j + gap) <= 0) break; (*pexch)(j, j + gap); } } } } /****************************************************************************/ /* int swap_byte(n) swap the bytes of n. the library (*.obj) version has had "swab %4" added at the appropriate place by changing "mov 4(%5),%4\nmov %4,%0\n" to "swab 4(%5)\nmov 4(%5),%0\n" in the .mac version before assembly. contained in hclib2. register int n; { return(n); } */ /****************************************************************************/ int *stopcc(n) /* calls .scca to inhibit control-c if n = 1, restore it if n = 0. it returns the address of termstat. if bit 15 of termstat is set a double cc has come in. it may put a control-c into terminal input buffer. emt call is to .scca */ int n; { INTERN termstat; int emt375(); if(n == 1) { if(emt375(016400, &termstat) < 0) /* turn off */ errfmt("stopcc - unable to turn off\n"); } else if(n == 0) { if(emt375(016400, 0) < 0) /* restore */ errfmt("stopcc - unable to turn on\n"); } else errfmt("stopcc - incorrect argument\n"); return(&termstat); } /****************************************************************************/ int substi(s, t, u, temp) /* if string s contains string t as a substring, substitute string u in place of string t. uses string temp as a temporary buffer. temp must thus be at least as long as (s - t) + u. returns YES if a substitution made, NO if not (to tell if all done if more than one t. */ register char *u; char *s, *t, *temp; { register i, k; unsigned substr(), lenstr(); i = substr(s, t); /* is it there? */ if(!s[i]) return(NO); /* if not, return unchanged */ for(k = 0; k < i; k++) temp[k] = s[k]; /* copy s prior to t */ while(*u != '\0') temp[k++] = *u++; /* copy in u */ for(i += lenstr(t); s[i] != '\0'; i++) temp[k++] = s[i]; /* copy in rest of s */ for(i = 0; i < k; i++) s[i] = temp[i]; /* copy temp back to s */ s[i] = '\0'; /* terminate it */ return(YES); } /****************************************************************************/ long tics() /* returns the elapsed time in tics since the last call. there are no arguments. the first call should be ignored, it is just for setup. emt call is to .gtim */ { long newtics, difftics; static long oldtics; int emt375(); if(emt375(010400, &newtics) < 0) { /* # of tics from midn */ errfmt("tics - failed\n"); return; } difftics = newtics - oldtics; /* minus tics at last call */ oldtics = newtics; return(difftics); } /****************************************************************************/ int time(n) /* returns the time into a four int array, first the hours, second the minutes, third the seconds, fourth the tics. the argument is the address of int a[4]; emt call is to .gtim */ FAST *n; { FAST k; LONG i; int emt375(); for(k = 0; k < 4; k++) /* zero the array */ *(n + k) = 0; if(emt375(010400, &i) < 0) { errfmt("time - failed\n"); return; } *n++ = i / 216000; /* hours */ *n++ = (i %= 216000) / 3600; /* minutes */ *n++ = (i %= 3600) / 60; /* seconds */ *n++ = i % 60; /* tics */ } /****************************************************************************/ char *timest(s) /* return system time as a string "hh-mm-ss" */ char *s; { int k, n[4], time(); unsigned decode(); time(n); k = decode(s, 2, "%i", n[0]); s[k++] = ':'; k += decode(&s[k], 2, "%i", n[1]); s[k++] = ':'; k += decode(&s[k], 2, "%i", n[2]); s[k] = '\0'; return(s); } /****************************************************************************/ int ttflush() /* flush the rt11 input ring buffer */ { } /* to build: compile the code above, then add the indented lines below .enabl gbl .globl ttflus .psect c$text,i,ro .mcall .ttinr ttflus: mov %5,-(%6) mov %6,%5 mov 44, r2 ; set bit 6 of JSW mov 44, r1 bis #100, r1 mov r1, 44 ; label: .ttinr bcc label ; go until no more characters ; mov r2, 44 ; reset JSW jmp c$rets .title ttflus .psect c$text,i,ro .even .psect c$data,d,rw .even .end */ /****************************************************************************/ int ungetc(fp, c) /* return character c to the file with FIO of address fp. I am not certain how this functions when getc() has just read across a block boundary. if getc() does not read the disk until the next call after exhausting the buffer, it should be ok. */ FIO *fp; int c; { if(fp->_fmode == READ) { if(fp->_nleft < BUFSIZE) { /* space to put it */ *--fp->_pnext = c; /* push char on buffer */ fp->_nleft++; /* bump the number left */ } } } /****************************************************************************/ int *unload(buffer, topadd, s, add1) /* unloads the buffer starting at buffer with a (top address + 1) of topadd into any combination and number of int, float, or double. add1 is the address of the first variable. char can not be done. it returns the next address to use (in buffer) unless there is an error in which case 0 is returned. s is a string indicating how to do the unload. i indicates int, f float, d double. all other characters are ignored. */ register int *buffer; int *topadd; char *s; unsigned add1; { register int *p; /* pointer to next address */ char *message; int *pi, lstoi(); float *pf, lstof(); double *pd, lstod(); message = "unload - attempt to unload past end of buffer\n"; p = &add1; /* point p at top of argument stack */ while(*s) { /* read the 'format' string */ switch(*s++) { case 'i' : if(buffer >= topadd) { errfmt(message); return(0); } pi = *p++; *pi = lstoi(buffer++); break; case 'f' : if(buffer + 1 >= topadd) { errfmt(message); return(0); } pf = *p++; *pf = lstof(buffer++); buffer += 1; break; case 'd' : if(buffer + 3 >= topadd) { errfmt(message); return(0); } pd = *p++; *pd = lstod(buffer++); buffer += 3; break; default : break; } } return(buffer); } /****************************************************************************/ int unrad(irad50, s) /* place in the 3 character string starting at s, the radix 50 contents of irad50. */ unsigned irad50; char *s; { int i; s[2] = irad50 % 050; /* r50 = (c0*050)+c1)*050+c2 */ irad50 = (irad50 - s[2]) / 050; s[1] = irad50 % 050; s[0] = (irad50 - s[1]) / 050; /* * convert radix 50 code to ascii */ for(i = 0; i < 3; i++) { if(s[i] <= 0 || s[i] > 047 || s[i] == 035) /* all illegals - 040 */ s[i] = ' '; else if(s[i] <= 032) /* A - Z, 1 - 32, 101 - 132 */ s[i] += 0100; else if(s[i] == 033) /* $, 033, 044 in ascii */ s[i] += 011; else /* s[i] >= 34 */ s[i] += 022; /* . or 0 - 9, 34&36-47 */ } } /****************************************************************************/ double dabs(x) /* returns absolute value of double x */ double x; { if(x >= 0) return(x); else return(-x); } /****************************************************************************/ double darctan(x) /* returns the arctangent of x (tan-1 x). the algorithm is from a DG fortran V manual */ double x; { int rflag = NO, zflag = NO, sign; double polynm(), r; static double a[] = { 197.20309568493503, 298.92803806939622, 125.57916643798066, 12.595802263029547}, b[] = { 197.20309568493503, 364.66240329770776, 207.69268173360463, 37.066086322091024, 1.0000}; sign = (x >= 0)? 1 : -1; if((x = abs(x)) > 1.) { rflag = YES; x = 1. / x; } if(x > 0.57735026918962577) { /* sqrt(3)/3 */ zflag = YES; x = (0.7320508075688773 * x - 1 + x) / (x + 1.7320508075688773); } /* sqrt(3), sqrt(3)-1 */ r = x * polynm(3, x, a) / polynm(4, x, b); if(zflag == YES) r += 0.52359877559829887; /* PI/6 */ if(rflag == YES) r = 1.5707963267948966 - r; /* PI/2 */ return(sign * r); } /****************************************************************************/ double dcos(x) /* returns the cosine of x. the algorithm is from an old DG fortran V manual. calls sinsub or cossub. */ double x; { int q, sign = 1; double r, sinsub(), cossub(); x = abs(x); x /= 0.7853981633974483; /* PI/4 */ q = (int) x; r = x - q; q %= 8; if(q >= 2 && q <= 5) sign = -1; if(q % 2) r = 1 - r; if(q == 0 || q == 3 || q == 4|| q == 7) return(sign * cossub(r)); else return(sign * sinsub(r)); } double cossub(x) /* used by dsin and dcos */ double x; { double polynm(); static double b[] = { 0.9999999999999999443, -0.30842513753403723, 0.15854344243734568E-01, -0.32599188645404001E-03, 0.359085912336036E-05, -0.2460945716614E-07, 0.11363812697E-09}; return(polynm(6, x, b)); } /****************************************************************************/ double dexp(x) /* returns the exponential (e to the x power) of double x. if error returns 0. the algorithm is from an old DG fortran V manual. */ double x; { int m, n, q, signfl = NO; double z, t, twoz1, twoz2, result; double power(), polynm(); static double a[] = { 1513.9067990543389, 20.202065651286927, 0.23093347753750234E-01}, b[] = { 4368.2116627275585, 233.18421142748162, 1.0}; if(!x) return(1.); if(x < 0.) { signfl = YES; x = -x; } if(x > 88.029691) { /* DG says 174.67308950110618 */ errfmt("dexp - argument too large\n"); return(0); } t = x * 1.4426950408889634; /* log2 e */ q = (int) t; z = t - q - 0.5; n = q % 4; m = (q - n) / 4.; twoz1 = z * polynm(2, z, a) / polynm(2, z, b); twoz2 = (.5 + .5 + twoz1) / (1 - twoz1); result = power(16., m) * power(2., n) * twoz2 * 1.4142135623730950; if(signfl == YES) return(1. / result); return(result); } /****************************************************************************/ double d10exp(x) /* returns 10 to the x power */ double x; { double dexp(); return(dexp(x * 2.302585092994046)); /* loge 10 */ } /****************************************************************************/ double dloge(x) /* returns the natural (base e) log of double x. if error returns 0. the algorithm is from an old DG fortran V manual. */ double x; { int n, m; double y, aa, logf, f = 0., r = 0., recp16 = 0.0625; double power(), polynm(); static double a[] = { -90.174691662040536, 93.463900642858538, -18.327870372215932}, b[] = { -45.087345831020306, 61.761065598471302, -20.733487895513939, 1.0}; if(x <= 0) { errfmt("dloge - argument is less than of equal to 0\n"); return(0); } if(x == 1) return(0.); if(x >= recp16) { /* next 10 lines */ for(n = 0; r < recp16 || r >= 1.; n++) /* x = 16to n *r */ r = x / power(16., n); n--; } else { for(n = 0; r < recp16 || r >= 1.; n--) r = x / power(16., n); n++; } for(m = 0; f < .5 || f >= 1.; m++) f = r * power(2., m); if(f <= 0.70710678118654752) { /* sqrt(.5) */ aa = 0.5; } /* --exponent ie(4n - m) by not --m */ else { aa = 1.0; --m; } y = (f - aa) / (f + aa); logf = y * polynm(2, y, a) / polynm(3, y, b); return(0.69314718055994531 * (4 * n - m) + logf); /* loge 2 */ } /****************************************************************************/ double dlog10(x) /* returns the log base 10 of double x. the algorithm is from an old DG fortran V manual. */ double x; { double dloge(); return(dloge(x) * 0.43429448190325183); /* log10 e */ } /****************************************************************************/ double dpower(y, x) /* returns the x-th power of y (ie y to the x). if error returns 0. the algorithm is from an old DG fortran V manual. */ double y; double x; { double dloge(), dexp(); if(y < 0 && x != 0) { errfmt("dpowr - argument is less than zero\n"); return(0); } if(y == 0 && x<= 0) { errfmt("dpowr - argument is zero\n"); return(0); } x *= dloge(y); if(abs(x) >= 88.029691) { /* DG says 174.67308950110618 */ errfmt("dpowr - arguments (x * ln(y)) too large\n"); return(0); } return(dexp(x)); } /****************************************************************************/ double droot(x, y) /* returns the x-th root of y */ double x; double y; { double dpower(); return(dpower(y, 1. / x)); } /****************************************************************************/ double dsin(x) /* returns the sine of x. the algorithm is from an old DG fortran V manual. calls sinsub or cossub. */ double x; { int q, sign = 1; double r, sinsub(), cossub(); if(x < 0) x = -x + PI; x /= 0.7853981633974483; /* PI/4 */ q = (int) x; r = x - q; q %= 8; if(q >= 4) sign = -1; if(q % 2) r = 1 - r; if(q == 0 || q == 3 || q == 4|| q == 7) return(sign * sinsub(r)); else return(sign * cossub(r)); } double sinsub(x) /* used by dsin and dcos */ double x; { double polynm(); static double a[] = { 0.78539816339744831, -0.80745512188280530E-01, 0.24903945701887361E-02, -0.36576204158455695E-04, 0.31336162166190400E-06, -0.1757149292755E-08, 0.6877100349E-11}; return(x * polynm(6, x, a)); } /****************************************************************************/ double dsqrt(x) /* return the square root of double x. the algorithm is from an old DG fortran V manual. if error returns -1. */ double x; { int i, p, m, n; double y, power(), r = 0, recp16 = 0.0625; double a = 0.22222222, b = 0.88888889; if(!x) return(0.); if(x < 0) { errfmt("sqrt - argument less than zero\n"); return(-1); } if(x >= recp16) { /* next 10 lines */ for(p = 0; r < recp16 || r >= 1.; p++) /* x = 16to p *r */ r = x / power(16., p); p--; } else { for(p = 0; r < recp16 || r >= 1.; p--) r = x / power(16., p); p++; } y = a + b * r; for(i = 0; i < 4; i++) /* iterate 4 times */ y = (r / y + y) / 2.; n = p % 2; /* convert p to 2m + n */ m = (p - n) /2.; return(power(16., m) * power(4., n) * y); } /****************************************************************************/ float fabs(x) /* returns absolute value of double x */ float x; { if(x >= 0) return(x); else return(-x); } /****************************************************************************/ int iabs(x) /* returns absolute value of double x */ int x; { if(x >= 0) return(x); else return(-x); } /****************************************************************************/ double polynm(order, x, cf_add) /* calculates a polynomial of degree order in x2 using a double array of coefficients with a starting address of cf_add. it returns the result. */ register int order; double x; double *cf_add; { register int i; static double lim_low[] = { 1e-19, 3.16228e-10, 4.64159e-7, 1.7783e-5, 1.58489e-4, 6.81292e-4 }; static double lim_high[] = { 1e19, 3.16228e9, 2.15443e6, 5.6234e4, 6.3096e3, 1.4678e3 }; double xabs, poly; double power(); xabs = abs(x); poly = *cf_add++; /* 0th order */ for(i = 0; i < order;) { if(xabs < lim_low[i] || xabs > lim_high[i]) break; /* if power would overflow */ poly += *cf_add++ * power(x, ++i * 2); } return(poly); } /****************************************************************************/ double power(x, pow) /* raises a double x to the int power pow */ double x; int pow; { double power = 1; if(pow >= 0) { for( ; pow > 0; pow--) power *= x; return(power); } else { for( ; pow < 0; pow++) power /= x; return(power); } } /****************************************************************************/