DECUS C LANGUAGE SYSTEM Cmath Library Reference Manual by Hamish Ross Department of Physiology Medical School Birmingham University BIRMINGHAM B15 2TJ This document describes an implementation of the elementary maths functions for the DECUS C language. Version of 15-Mar-85 Copyright (C) 1985, DECUS. General permission to copy or modify, but not for profit, is hereby granted, provided that the above copyright notice is included and reference made to the fact that reproduction privileges were granted by DECUS. The information in this document is subject to change without notice and should not be construed as a commitment by Digital Equipment Corporation or by DECUS. Neither Digital Equipment Corporation, DECUS, nor the author assumes any responsibility for the use or reliability of this document or the described software. This software is made available without any support whatsoever. UNIX is a trademark of Bell Telephone Laboratories. RSX, RSTS/E, RT11 and VMS are trademarks of Digital Equipment Corporation. Cmath Library Reference Manual Page 4 1.0 Introduction and Scope ____________ ___ _____ This collection of programs provides a library of the elementary maths functions for the DECUS C system. It differs from the PML library, supplied on the 11-SP-18 tape, in the major design objectives which it attempts to meet. These are, to provide a set of functions which have the same names and functionality as their counterparts in UNIX, and therefore in a number of other C compilers, and which have a high level of accuracy throughout the entire range of their arguments. A collection of algorithms which seems to satisfy the latter objective is published in the book "Software Manual for the Elementary Functions" by William J. Cody, Jr. and William Waite, Prentice-Hall Series in Computational Mathematics (1980). This kit is an implementation of all the algorithms in that book together with some of the closely associated functions in the UNIX math library. The functions provided are sqrt, log, log10, exp, pow, sin, cos, tan, cotan, asin, acos, atan, atan2, sinh, cosh and tanh, based on the Cody and Waite algorithms, together with the additional functions fabs, floor, ceil, frexp, ldexp and modf. Since there were some problems with the i/o functions atof and $$dtoa in the 11-SP-18 set, particularly in dealing with double precision numbers, new versions of those, written in C, have been included also, as has an error handling package. There were found to be a number of bugs in the supplied 11-SP-18 kit, in the handling of the /d switch by the compiler and in the handling of some i/o buffers in the run-time library. These bugs have been fixed and the patches are included in this kit. More details about that are included in section 4 and in the text files (.TXT) associated with each of the patches. Error handling is described more fully in the next section and portability and maintainability in section 3. Instructions for building the kit are provided in section 4, and further information about accuracy and testing in section 5. A full description of each of the functions, including useage instructions, the meanings of all diagnostic messages and notes about any differences between the actual implementation and the one described in Cody and Waite, appears in section 6. An appendix contains a KWIK index. 2.0 Error handling _____ ________ In dealing with errors, a number of somewhat incompatible aims had to be reconciled. On the one hand the UNIX convention seems to be to do relatively little about errors except to set the EDOM and ERANGE flags for domain and range values which lie outside the permissible limits, and to generate a floating exception signal for floating point exceptions. Having become used to something rather more friendly in the Fortran system on the PDP-11, it seemed desirable to do more than the UNIX minimum, but not in such a way as to cause problems in moving Cmath Library Reference Manual Page 5 programs between different systems. Error handling is therefore done by calling a common function cmemsg(), with two arguments; one is an error code number and the other is a pointer to some other piece of useful data, such as an argument value. A full list of the error codes and their associated messages is given in the documentation for the function cmathe(), which initialises and reinitialises the error handling mechanism. Cmathe() must be called at least once to start up the error handling system. To avoid having to remember to do this at the beginning of each main program, and to eliminate the need for a program edit when moving between operating systems a call to cmathe() has been included in an initialiser program called cminit(). Thus, as long as cminit is loaded with each program, everything will be taken care of automatically. Cminit plays a role similar to SUPORT in the LINK argument for RT11 and it could easily be included within it, but it seemed tidier to keep the two separate. Cminit also contains a call to the CLIB function fps(), which sets the floating point exception flags, and code, which is conditionally compiled for RT11 only, to set up a vector to catch the floating point exception interrupts. That is done in $$traps for other operating systems but not for RT11, so these errors are caught and handled by the cmemsg() error handling function also. Further details about error handling, including instructions on how to add a new error to the system, are to be found in the documentation for cmathe and cminit in section 6. 3.0 Portability and Maintainability ___________ ___ _______________ 3.1 Portability of Application Programs ___________ __ ___________ ________ As mentioned in section 1, one of the main objectives of this package was to facilitate the movement of application programs between other C math packages and this one. Therefore, all functions have been given the same names and argument lists as their counterparts in UNIX. In all cases the functions use doubles rather than floats to pass arguments and return values. Since DECUS C does not always use the UNIX convention of promoting floats to doubles, for argument passing and return of results, there is an source of potential problems there. The simple solution for new programs, is to ensure that all functions which return a double result are declared double rather than float and likewise for any float/double arguments. This method must be used if the DECUS compiler /f switch is used. In the case of existing programs which declare their maths functions and their arguments to be of type float everything should still work properly if the /d switch of the DECUS compiler is used. However, in practice, it does not, Cmath Library Reference Manual Page 6 unless the compiler is patched. In the supplied version all the type conversions are not applied automatically. A patch kit for the CC103 module of the compiler which fixes that problem is included with this kit (see the Build instructions in section 4 for details). One further consequence of the DECUS compiler's acceptance of arguments which are floats rather than doubles is that printf and friends need to be able to distinguish those two types. That is done by using lower case format symbols for float conversions, e.g., "%f" and upper case for doubles, e.g., "%F". The latter is not a normal usage for other C compilers and may thus cause portability problems. Two solutions are offered. One is to abandon the DECUS /f facility and to stick to the more widely used convention. In that case the /d switch must be made the default (and the patched version of CC103 used), then a patched version of the support library function doprnt which does not distinguish between upper and lower case format must be used also. A patch for that purpose is included in this kit (DOPRNT.TX2). The other solution to the upper/lower case format problem is to ensure that the correct case is provided for whatever maths/io package is being used. To facilitate conversion between the two format conventions the 2 utility programs FORMUC.C and FORMLC.C are supplied. These are C programs built using LEX. (The 11-SP-18 version of LEX seems to have the outermost layer of braces missing, this causes it to fail while decoding the command line). FORMUC will convert all formats to upper case, FORMLC to lower case. Mixed case programs will require hand editing. 3.2 Portability of the cmath package ___________ __ ___ _____ _______ The portability of this package to other computers was not one of the primary aims of this exercise. The material in Cody and Waite shows that to obtain the best levels of accuracy it is necessary to know how the target machine does its arithmetic and thus portability and accuracy are to some extent incompatible. In this case accuracy has been given precedence - this is a package designed to run on PDP-11s with FPU floating point hardware. On the other hand portability has not been ignored. The programs have been kept reasonably clean and tidy so it is easy to describe where the machine dependencies lie. The main set of functions, that is to say the ones described as the Cody and Waite set in section 1, together with ceil, floor and fabs do not contain any explicit assumptions about the form of the floating point representation of a number within the machine. Where they need to use such information, as, for example in extracting the mantissa and exponent of a number or the integer and fractional parts they use the machine dependent functions ldexp(), frexp() or modf() which are written in MACRO. In the functions sqrt, sin, cos, tan and cotan it is necessary to determine whether or not an integer, n, was odd. That has been done with a test on n & 1 which makes an assumption about the Cmath Library Reference Manual Page 7 representation of integers. A package of this type inevitably contains a great many magic numbers, as the constants for the polynomials used to calculate the functions in the basic range and as other factors such as pi, e or the square root of three. Constants such as the latter, which are used in a number of places, reside in the constants file math.h. That file also contains machine dependent constants such the largest and smallest numbers representable in the machine. (see the documentation on math.h in section 6 for further details). Other constants, which are specific to individual functions, are initialised as static variables or in #define statements at the head of the function. 3.3 Note to maintainers/improvers ____ __ _____________________ A casual examination of the code will reveal some inefficient constructs like x = (y - 0.5) - 0.5. Do not 'improve' them unless you really understand the consequences of your action. In many cases the code is written in what might appear to be an odd way because Cody and Waite say that is the way it should be done to get the best accuracy from this type of hardware. Therefore, if you are attempting to use other hardware with a different type of arithmetic processor (decimal, octal or hexadecimal) you will need to make many changes to the code to retain the same level of accuracy in the algorithms. It would be prudent not to make any changes without a careful reading of Cody and Waite first. 4.0 Building ________ Before compiling the functions and installing them in a library it is advisable to apply some patches to the original 11-SP-18 kit. These patches are included in this package in the files CC103.TXT, DOPRNT.TXT, DOSCAN.TXT, FCLOSE.TXT and IOABUF.TXT each of which includes an explanation of the purpose of the patch and the SLP text to implement it. Very briefly the reasons for each of the patches is as follows:- CC103 - floats passed to and from functions not promoted to doubles on /d. DOPRNT - increase buffer size for doubles and pass size to $$dtoa. DOPRNT - optionally treat e, f, g formats as E, F, G. DOSCAN - bad handling of EOF means malloc data overwritten. FCLOSE - incomplete terminal line not flushed on closedown (RT-11). IOABUF - tty buffer can overfill on lines > 80 chars (RT-11). All the C files of the math package contain the /*)LIBRARY command for BUILD. In addition, the MACRO functions frexp, ldexp and modf (and in the case of RT11, fptrap also) must be Cmath Library Reference Manual Page 8 assembled. The main decision which must be made is where to put the functions. If almost all applications need to use the package then it is simplest to add them to CLIB, in which case modules ddtoa and atofd should be removed. If it is desired to leave CLIB unaltered put the new functions in a separate library, and ensure that it is called before CLIB. In that case copies of atof.obj and $$dtoa.obj should be left on device C: and should be protected. In either case a copy of cminit should be left on that device also, remember it plays a role similar to SUPORT in RT11. A command file which does all these things for RT11 is included in the kit as cmath.com. As it stands, it will build a system with CMATH separate from CLIB. Removing the comment marks on the indicated sections will incorporate CMATH in CLIB. The penalty for including CMATH in CLIB is not only longer search times but inclusion of larger versions of $$dtoa and atof in programs which do not use floating point. The benefit is simpler link instructions and fewer libraries to search for programs which use floating point. 5.0 Accuracy and testing ________ ___ _______ All the test programs from the Cody and Waite book have been included in the kit. They are in the archive file cwtest.arc. They were transcribed straight into C from the Fortran listings in the book, without any effort at making them portable. There is a program, machar, which obtains values of machine dependent constants which the other programs need. It seems to work ok, but it is worth checking its output against the FPU manual. The other programs do a series of consistency checks on each of the functions using sets of random arguments which exercise different parts of the algorithm so that inaccuracies may tracked to a subcomponent of the algorithm. The test programs in the kit all have BUILD instructions for a CMATH library incorporated into CLIB in which case commands of the type RUN BUILD TSQRT will run the test. If CMATH is separate the BUILD sequence should be modified to /*)BUILD $(LIBS) = {cmath cminit} $(DTOA) = 1 $(ATOF) = 1 */ and again RUN BUILD TSQRT will run the test. In each case no args are needed in the command line, but if argv[1] is supplied with a string 'xxxx' then the message 'Performance test on xxxx' will appear at the head of the output. The root mean square relative error for all the random number tests was less than one binary digit, and the maximum relative error less than three bits, in most cases it was less than 2 bits. The package has been developed and tested on a PDP11/34 with FP11-A FPU. The operating system was RT-11 v4. It has not been tested on any other operating system. The only part which Cmath Library Reference Manual Page 9 should be operating system dependent is the code for handling FPU interrupts. It is conditionally included if the operating system is RT-11 (see documentation for cminit() in section 6). It is assumed that $$traps in CLIB will deal with floating point exceptions in other operating systems. 6.0 Function descriptions ________ ____________ This section gives useage, description, diagnostic and some internal information for each of the supplied functions. It also gives some further background information about the error handling package cmathe(), the initialisation function cminit() and the maths constants file math.h. The fullest description of the meanings of the error messages is found here too, in the diagnostics section of the function concerned. For examples of the use of the functions see the listings of the test programs in the file cwtest.arc. Cmath Library Reference Manual Page 10 $$dtoa double to ascii conversion 6.1 double to ascii conversion ______ __ _____ __________ ********** * $$dtoa * ********** File name: dtoa.c NAME: $$dtoa -- double to ascii conversion USAGE: char buff[]; /* space for result string */ char conv; /* conversion type 'e', 'f', or 'g' others default to e */ int bsize; /* size of buffer */ int dplace; /* places after the point */ double value; /* value to converted */ $$dtoa(buff, conv, bsize, dplace, value); DESCRIPTION: $$dtoa is a (temporary?) replacement for the program of the same name on the 11-SP-18 tape and described on page 5-9 of the Wizard document. The rest of this description is a copy of that. $$dtoa() is called by printf to convert a double precision number to ascii. The text is written to buff[]. Bsize is the maximum size of the result string; it is determined by the amount of space allocated to buff by the calling program (normally doprnt). Dplace is the number of decimal places after the decimal point. If the result is too large to fit into the space allocated to buff, digits are lost from the low end. With the buffer size of 30 bytes allocated by the modified version of doprnt() (edit b1) only rubbish digits can be lost in this way. Conv is 'e', 'f' or 'g' to define the format. DAIGNOSTICS: None INTERNAL: The program uses the following algorithm:- set default value of dplace Cmath Library Reference Manual Page 11 $$dtoa double to ascii conversion strip off sign scale and compute no of leading digits decide if 'g' goes to 'e' or 'f' compute no of digits to print change 'f' to 'e' if insufficient space for result scale to range 1.0 - 10.0 and round up build digit string if not 'f' print exponent part print final null AUTHOR: Hamish Ross. DATE: 9-Mar-85 Cmath Library Reference Manual Page 12 arsico arc sin and arc cosine functions 6.2 arc sin and arc cosine functions ___ ___ ___ ___ ______ _________ ********** * arsico * ********** File name: arsico.c NAME: arsico -- arc sin and arc cosine functions USAGE: double x, f, asin(); double g, acos(); f = asin(x); g = acos(x); DESCRIPTION: asin returns arcsin of x. acos returns arccosine of x. DIAGNOSTICS: If the absolute value of the argument is greater than 1 the error message 'asin, acos arg bad', followed by the value of the argument, is written to stderr. The algorithm returns a result of HUGE. INTERNAL: Algorithm from Cody and Waite pp. 174-193. AUTHOR: Hamish Ross. DATE: 8-Jan-85 Cmath Library Reference Manual Page 13 atan arctan function 6.3 arctan function ______ ________ ******** * atan * ******** File name: atan.c NAME: atan -- arctan function USAGE: double f, x, atan(); f = atan(x); double g, u, v, atan2(); g = atan2(v, u); DESCRIPTION: atan(x) returns the arctangent of x. atan2(v, u) returns the arctangent of v / u. DIAGNOSTICS: If both arguments of atan2 are zero the result is indeterminate. The message 'atan2 args both zero' followed by the value of the denominator, is written to stderr and a value of zero is returned. INTERNAL: Algorithm based on Cody and Waite pp. 194-216. The algorithm here has been modified slightly to avoid the test for overflow in v/u. By checking which of u or v is greater in magnitude all overflows become underflows which are set to zero without a message. This has meant that the division of the work between atan(), atan2() and the internal function _atan() is slightly different from that suggested in Cody and Waite. AUTHOR: Hamish Ross. DATE: Cmath Library Reference Manual Page 14 atan arctan function 1-Jan-85 Cmath Library Reference Manual Page 15 atof ascii to double conversion 6.4 ascii to double conversion _____ __ ______ __________ ******** * atof * ******** File name: atof.c NAME: atof -- ascii to double conversion USAGE: double d, atof(); char *buffer; d = atof(buffer); DESCRIPTION: This program is a (temporary?) replacement for the one of the same name described on page 5-52 of the Wizard document on the 11-SP-18 tape. The type of conversion from ascii string to double floating point number given there is reproduced here. The BNF for numbers which can be decoded by that routine is :- := := + | - | := | | := | . | . | . := | := := E | e | D | d := 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 This routine differs slightly from the above, in that a number may not start at the exponent part without a leading digit. That means that the description of a real nuber is now := | The reason for doing it this way is that it stops scanf("%f", &x) eating its way into a following string if that string starts with a 'd' or an 'e'. It will still do that if the following string begins with a '.'. K & R does not define what atof's idea of a floating number is, but the above definition is closer than the DECUS C version to their definition of a floating constant. The conversion stops on the first non-legal character, but unlike the original Cmath Library Reference Manual Page 16 atof ascii to double conversion version this one prints an error message and the offending string through the cmathe error reporting system. DIAGNOSTICS: If a non-legal character is encountered the message 'illegal character in floating input string', followed by the entire string in the buffer argument, is written to stderr. The value of the number built up to that point is returned. If the input produces a number which would result in floating overflow the message 'atof input too large', followed by the entire string in the buffer argument, is written to stderr. The value HUGE is returned. INTERNAL: Algorithm is straightforward, see comments in source code for details. AUTHOR: Hamish Ross. DATE: 22-Aug-84 original version 25-Mar-85 insist numbers start with a digit not 'e' or 'd'. Cmath Library Reference Manual Page 17 ceil smallest integer not less than argument 6.5 smallest integer not less than argument ________ _______ ___ ____ ____ ________ ******** * ceil * ******** File name: ceil.c NAME: ceil -- smallest integer not less than argument USAGE: double x, y, ceil(); y = ceil(x); DESCRIPTION: Returns the smallest integer not less than the argument. DIAGNOSTICS: If the argument is greater than the integer part of HUGE a floating point overflow exception will occur. INTERNAL: AUTHOR: Hamish Ross. DATE: 31-Jan-85 Cmath Library Reference Manual Page 18 cmathe cmath error package 6.6 cmath error package _____ _____ _______ ********** * cmathe * ********** File name: cmathe.c NAME: cmathe -- cmath error package USAGE: There are three function calls associated with this package. int ercode; cmathe(ercode); int erno; char *datap; cmemsg(erno, datap); cmelog(); DESCRIPTION: This set of functions provides an error reporting package for the floating point math package. For RT-11 it catches the FP error interrupts also. It may be connected to any other function which needs to report errors (such as square root of negative number). It provides the option, present in Fortran, of reporting n errors and then either aborting or ignoring subsequent errors. In that case they are counted but no messages are printed. Unless there is an abort, control is returned to the calling function, which must be able to continue or else abort itself. The cmathe system is started or restarted by the call cmathe(n); where is n is usually a literal constant. If n is 0 all errors are reported and there is no aborting. If n > 0, n errors are reported then there is an abort, so cmathe(1) will abort after reporting the first error. If n < 0, n errors are reported then the reporting Cmath Library Reference Manual Page 19 cmathe cmath error package stops, though the errors are still counted. The errors are reported on stderr, and on an abort a calltrace is done to stderr. This produces a trace of function calls leading to the error if profiling has been enabled, otherwise it produces a sequence of base addresses of the functions called. The program then exits with exit code 4. Cmathe sets all the error counts to zero. The cmath initialisation function cminit() calls cmathe with argument zero, so that if that argument is appropriate there is no need to call cmathe explicitly. Cminit also calls fps, and under RT-11 takes over the FPU interrupt vector to deal with the FPU exceptions. See the documentation on cminit for further details. The function which invokes the error reporting mechanism is cmemsg(erno, datap); where erno is a symbol which defines the error number. These definitions are in math.h. Datap is a pointer to an object whose value might help to throw further light on the cause of the error. Further information on both those items may be found in the internal section below. The cumulative error list may be printed on stdout via the call cmelog(); This call does not clear the error counts, to do that use cmathe. DIAGNOSTICS: If cmemsg() is called with an error number which is undefined or out of range the message 'error code out of range' is written to stderr. INTERNAL: The error table currently looks like :- Error Caller Sev Eno Message Argument FP_OPER FPU F E floating op code error PC FP_ZDIV FPU F E floating divide by zero PC FP_FTOI FPU F E floating integer PC conversion error FP_OFLO FPU F E floating overflow PC Cmath Library Reference Manual Page 20 cmathe cmath error package FP_UFLO FPU W E floating underflow PC FP_UDEF FPU W E floating undefined PC variable FP_BIGI atof F R atof input too large string FP_BADC atof F D illegal character in string floating input string FP_NESQ sqrt F D sqrt arg negative value FP_LEXP exp F R exp overflow value FP_SEXP exp W R exp underflow value FP_NLOG log F D log arg negative or value zero FP_TANE tancot F D tan, cot arg too large value FP_TRIG sincos F D sin, cos arg too large value FP_ATAN atan2 F D atan2 args both zero value FP_COTE cotan F D cotan arg too small value FP_ARSC arsico F D asin, acos arg bad value FP_SINH sinh F D sinh arg too large value FP_COSH cosh F D cosh arg too large value FP_POWN pow F D pow arg negative value FP_POWO pow F R pow overflow value FP_POWU pow W R pow underflow value Error is the error number as defined in the math.h file. Caller is the name of the function in which the error was encountered, or it is FPU for floating point exceptions. Sev is the severity of the error. Errors of severity W(arning) do not print their messages and their error counts do not count towards the abort limit. Their individual counts are updated and may be checked using cmelog(). Eno is the error number for the perror() function. Only three have been used E (SIGFPE) is the floating point exception error, D (EDOM) is a domain error (inadmissable argument) and R (ERANGE) is a range error (result out of range). Message is the error message printed and argument is an argument of the function which encountered the error, or the PC in the Cmath Library Reference Manual Page 21 cmathe cmath error package case of FPU errors. To add a new error enter the error message in the error message table in the cmathe.c file, define the error code in math.h, then place a cmemsg(erno, datap); statement in the function which needs to call the error system. Erno is the index number of the error in the error table (the first col in the table above) and datap is a pointer to the argument passed back (the last col in the table above). For FPU generated errors (error nos up to FP_NFPU) it points to an int (value of PC), for bad string errors (error nos up to FP_NFPU + FP_NSTR) it points to a string and for all higher values it points to a double precision argument. That organisation allows the addition of new errors passing any of those argument types simply by adding to the error tables in this file and in math.h. AUTHOR: Hamish Ross. DATE: 10-Feb-85 Cmath Library Reference Manual Page 22 cminit cmath initialisation procedure 6.7 cmath initialisation procedure _____ ______________ _________ ********** * cminit * ********** File name: cminit.c NAME: cminit -- cmath initialisation procedure USAGE: Do not call this procedure explicitly, but be sure that it is loaded by including its name in the linkers command string. The program will then be called during the normal initialisation procedure (via SUPORT in RT-11). DESCRIPTION: Used by the floating point math package to initialise the error handling system and to set the standard default parameters for the floating point status. If the operating system is rt11, it arranges to catch the interrupts from the FPU since that is not done in traps for rt11 (c.f. Wizard and source code). Though cmathe() and fps() are called here they can also be called at any stage by the users program to invoke different values of the parameters set by these functions. DIAGNOSTICS: None INTERNAL: AUTHOR: Hamish Ross DATE: 31-Jan-85 Cmath Library Reference Manual Page 23 cosh hyperbolic cosine function 6.8 hyperbolic cosine function __________ ______ ________ ******** * cosh * ******** File name: cosh.c NAME: cosh -- hyperbolic cosine function USAGE: double x, f, cosh(); f = cosh(x); DESCRIPTION: Returns hyperbolic cosine of argument. DIAGNOSTICS: If the argument is large enough to cause overflow in the result the message 'cosh arg too large', followed by the argument, is written to stderr. A value of HUGE is returned. INTERNAL: Algorithm from pp. 217-238 of Cody and Waite. The constant v which appears in 2 define statements is chosen so that log(v) is an exact octal constant slightly larger than log(2) (See C and W pp. 218 and 227. AUTHOR: Hamish Ross. DATE: 12-Jan-85 Cmath Library Reference Manual Page 24 exp exponential function 6.9 exponential function ___________ ________ ******* * exp * ******* File name: exp.c NAME: exp -- exponential function USAGE: double x, f, exp(); f = exp(x); DESCRIPTION: Returns exponential of x. DIAGNOSTICS: If the argument is large enough to cause overflow in the result the message 'exp overflow', followed by the value of the argument, is written to stderr. A value of HUGE is returned. If the argument would produce a result exponent too small to accomodate, a value of zero is returned. A warning message 'exp underflow', followed by the value of the argument, is written to stderr, if warnings are enabled (normally not). INTERNAL: Algorithm from Cody and Waite pp. 60-83. AUTHOR: Hamish Ross. DATE: 28-Dec-84 Cmath Library Reference Manual Page 25 fabs floating absolute value of floating argument 6.10 floating absolute value of floating argument ________ ________ _____ __ ________ ________ ******** * fabs * ******** File name: fabs.c NAME: fabs -- floating absolute value of floating argument USAGE: double x, y, fabs(); y = fabs(x); DESCRIPTION: Returns the absolute value of its floating argument. DIAGNOSTICS: None INTERNAL: AUTHOR: Hamish Ross. DATE: 31-Jan-85 Cmath Library Reference Manual Page 26 floor largest integer not greater than argument 6.11 largest integer not greater than argument _______ _______ ___ _______ ____ ________ ********* * floor * ********* File name: floor.c NAME: floor -- largest integer not greater than argument USAGE: double x, y, floor(); y = floor(x); DESCRIPTION: Returns the largest integer not greater than the argument. DIAGNOSTICS: If the argument is less than the integer part of -HUGE a floating point exception (overflow) will occur. INTERNAL: AUTHOR: Hamish Ross. DATE: 31-Jan-85 Cmath Library Reference Manual Page 27 frexp get mantissa and exponent of a number to base 2 6.12 get mantissa and exponent of a number to base 2 ___ ________ ___ ________ __ _ ______ __ ____ _ ********* * frexp * ********* File name: frexp.c NAME: frexp -- get mantissa and exponent of a number to base 2 USAGE: double x, mant, frexp(); int xi; mant = frexp(x, &xi); DESCRIPTION: This is the UNIX function frexp. It returns the mantissa of x as its value and the exponent of x as an integer via xi. DIAGNOSTICS: None INTERNAL: This does not exist as a C program. The program is written in macro since AS does not seem to handle the STEXP instruction properly. This file is included only for the documentation. AUTHOR: Hamish Ross. DATE: 23-Aug-84 Cmath Library Reference Manual Page 28 ldexp load mantissa and exponent of a number to base 2 6.13 load mantissa and exponent of a number to base 2 ____ ________ ___ ________ __ _ ______ __ ____ _ ********* * ldexp * ********* File name: ldexp.c NAME: ldexp -- load mantissa and exponent of a number to base 2 USAGE: double x, mant, ldexp(); int xi; x = ldexp(mant, xi); DESCRIPTION: This is the UNIX function ldexp. It returns the floating point number built from the exponent xi and mantissa mant. DIAGNOSTICS: None INTERNAL: This does not exist as a C program. The program is written in macro, as ldexp.mac, since AS does not seem to handle the LDEXP instruction properly. This file is included only for the documentation. AUTHOR: Hamish Ross. DATE: 23-Aug-84 Cmath Library Reference Manual Page 29 log logarithm function 6.14 logarithm function _________ ________ ******* * log * ******* File name: log.c NAME: log -- logarithm function USAGE: double x, f, log(); f = log(x); DESCRIPTION: Returns logarithm of x. DIAGNOSTICS: If the argument is negative or zero the message 'log arg negative or zero', followed by the value of the argument, is written to stderr. If the argument was zero a value of -HUGE is returned, otherwise the value of log(fabs(x)) is returned. INTERNAL: Algorithm from pp. 35-59 of Cody and Waite. AUTHOR: Hamish Ross DATE: 26-Dec-84 Cmath Library Reference Manual Page 30 log10 logarithm function base 10 6.15 logarithm function base 10 _________ ________ ____ __ ********* * log10 * ********* File name: log10.c NAME: log10 -- logarithm function base 10 USAGE: double x, f, log10(); f = log10(x); DESCRIPTION: Returns logarithm of x to base 10. DIAGNOSTICS: If the argument is negative or zero the message 'log arg negative or zero', followed by the value of the argument, is written to stderr. If the argument was zero a value of -HUGE is returned, otherwise the value of log(fabs(x)) is returned. INTERNAL: returns value of log(x) * log10(e) AUTHOR: Hamish Ross. DATE: 26-Dec-84 Cmath Library Reference Manual Page 31 modf modulus function 6.16 modulus function _______ ________ ******** * modf * ******** File name: modf.c NAME: modf -- modulus function USAGE: double x, xi, fr, modf(); fr = modf(x, &xi); DESCRIPTION: This is the UNIX function modf. It returns the fractional part of x as its value and the integer part of x as a double value in xi. The result has the same sign as x. DIAGNOSTICS: None INTERNAL: The program is written in macro as modf.mac. It uses the FPU MODD instruction. This C file is included only for the documentation. AUTHOR: Hamish Ross. DATE: 22-Aug-84 Cmath Library Reference Manual Page 32 pow x raised to power y 6.17 x raised to power y _ ______ __ _____ _ ******* * pow * ******* File name: pow.c NAME: pow -- x raised to power y USAGE: double x, y, f, pow(); f = pow(x, y); DESCRIPTION: Returns value of x raised to power y DIAGNOSTICS: There are three error possible error messages from this function. If the x argument is negative the message 'pow arg negative', followed by the value of x, is written to stderr. The value of pow for |x| is returned. If x = 0.0 and y <= 0.0 or if result overflows the message 'pow overflow', followed by the value of y, is written to stderr. The value of HUGE is returned. If the result underflows and if warnings are enabled (normally not), the message 'pow underflow', followed by the value of y, is written to stderr. The value of 0 is returned. The suggestion of Cody and Waite, that the domain be reduced to simplify the overflow test, has been adopted, consequently overflow is reported if the result would exceed HUGE * 2**(-1/16). 2**(-1/16) is approximately 0.9576. INTERNAL: Algorithm from Cody and Waite pp. 84-124. This algorithm required two auxiliary programs POWGA1 and Cmath Library Reference Manual Page 33 pow x raised to power y POWGA2 to calculate, respectively, the arrays a1[] and a2[] used to represent the powers of 2**(-1/16) to more than machine precision. The source code for these programs are in the files POWGA1.AUX and POWGA2.AUX. The octal table on page 98 of Cody and Waite is in the file POWOCT.DAT which is required on stdin by POWGA2. AUTHOR: Hamish Ross. DATE: 27-Jan-85 Cmath Library Reference Manual Page 34 powten table of powers of 10 6.18 table of powers of 10 _____ __ ______ __ __ ********** * powten * ********** File name: powten.c NAME: powten -- table of powers of 10 USAGE: extern double powten[]; DESCRIPTION: This is an extern array of the powers of ten from 1e-38 to 1e+38. It is used by dtoa and atof. DIAGNOSTICS: None INTERNAL: powten[0] is 1e-38 powten[38] is 1e0 powten[76] is 1e38 AUTHOR: Hamish Ross. DATE: 25-Mar-85 Cmath Library Reference Manual Page 35 sincos sin and cosine functions 6.19 sin and cosine functions ___ ___ ______ _________ ********** * sincos * ********** File name: sincos.c NAME: sincos -- sin and cosine functions USAGE: double x, f, sin(); double g, cos(); f = sin(x); g = cos(x); DESCRIPTION: sin returns sin of x. cos returns cosine of x. DIAGNOSTICS: If the absolute value of the argument is greater than TRIG_MAX the error message 'sin, cos arg too large', followed by the value of the argument, is written to stderr. The algorithm returns a result of lower accuracy. The number of reliable decimal digits will be roughly 20 - log10(|x / pi|). INTERNAL: Algorithm from Cody and Waite pp. 125-149. A small modification to the published algorithm has meant that the value of the argument is not restricted to values whose integer part can be represented by a long integer. The accuracy of the result in those cases will, however, be very low (c.f. above). AUTHOR: Hamish Ross. DATE: 8-Jan-85 Cmath Library Reference Manual Page 36 sincos sin and cosine functions 6.20 hyperbolic sine function __________ ____ ________ ******** * sinh * ******** File name: sinh.c NAME: sinh -- hyperbolic sine function USAGE: double x, f, sinh(); f = sinh(x); DESCRIPTION: Returns hyperbolic sine of argument. DIAGNOSTICS: If the argument is large enough to cause overflow in the result the message 'sinh arg too large', followed by the value of the argument, is written to stderr. A value of HUGE, with the same sign as the argument, is returned. INTERNAL: Algorithm from pp. 217-238 of Cody and Waite. The constant v which appears in 2 define statements is chosen so that log(v) is an exact octal constant slightly larger than log(2) (See C and W pp. 218 and 227. AUTHOR: Hamish Ross. DATE: 12-Jan-85 Cmath Library Reference Manual Page 37 sqrt square root function 6.21 square root function ______ ____ ________ ******** * sqrt * ******** File name: sqrt.c NAME: sqrt -- square root function USAGE: double x, f, sqrt(); f = sqrt(x); DESCRIPTION: Returns square root of argument. DIAGNOSTICS: If the argument is negative the message 'sqrt arg negative', followed by the value of the argument, is written to stderr. The square root of the absolute value of the argument is returned. INTERNAL: Algorithm from pp. 17-34 of Cody and Waite. Uses 3 Newton iterations. AUTHOR: Hamish Ross. DATE: 26-Dec-84 Cmath Library Reference Manual Page 38 tancot tangent and cotangent functions 6.22 tangent and cotangent functions _______ ___ _________ _________ ********** * tancot * ********** File name: tancot.c NAME: tancot -- tangent and cotangent functions USAGE: double x, f, tan(); double g, cotan(); f = tan(x); g = cotan(x); DESCRIPTION: tan returns tan of x. cotan returns cotangent of x. DIAGNOSTICS: If the absolute value of the argument is greater than TRIG_MAX the error message 'tan, cotan arg too large', followed by the value of the argument, is written to stderr. The algorithm returns a result of lower accuracy. The number of reliable decimal digits will be roughly 20 - log10(|x|). If the argument of cotan is less in magnitude than 1 / HUGE the error message 'cotan arg too small', followed by the value of the argument, is written to stderr. A value of HUGE with the same sign as the argument is returned. INTERNAL: Algorithm from Cody and Waite pp. 150-173. A small modification to the published algorithm has meant that the value of the argument is not restricted to values whose integer part can be represented by a long integer. The accuracy of the result in those cases will, however, be very low (c.f. above). AUTHOR: Hamish Ross. DATE: Cmath Library Reference Manual Page 39 tancot tangent and cotangent functions 8-Jan-85 Cmath Library Reference Manual Page 40 tanh hyperbolic tanh function 6.23 hyperbolic tanh function __________ ____ ________ ******** * tanh * ******** File name: tanh.c NAME: tanh -- hyperbolic tanh function USAGE: double x, f, tanh(); f = tanh(x); DESCRIPTION: Returns hyperbolic tanh of argument. DIAGNOSTICS: None INTERNAL: Algorithm from pp. 239-255 of Cody and Waite. AUTHOR: Hamish Ross. DATE: 30-Jan-85 Cmath Library Reference Manual Page 41 math.h cmath constants and definitions 6.24 cmath constants and definitions _____ _________ ___ ___________ ********** * math.h * ********** File name: math.h NAME: math.h -- cmath constants and definitions USAGE: #include DESCRIPTION: is a file of constants and definitions used by the cmath library routines. As such it retains the commonly occurring magic numbers from that package in one convenient location. If any of those definitions have to be modified the cmath functions will need to be recompiled also. In some cases the accuracy of these constants is critical to the accuracy of some component of the cmath package, so it would be wise not to touch them unless you understand the significance for that package of your actions. It will usually not be necessary to include this file in a program which uses the cmath functions. However, since it contains useful constants like PI, E, PI_BY_2, ROOT_2 and so on it may be included to provide those. Full details of its contents appear in the 'internal' section of this documentation. DIAGNOSTICS: Not applicable INTERNAL: contains 4 types of constants :- 1) Absolute constants such as pi and e. 2) Checked octal constants; these are absolute constants for which it has been verified that the compiler converts the decimal representation given here to a fully accurate binary number inside the machine. Cmath Library Reference Manual Page 42 math.h cmath constants and definitions This level of accuracy is essential for at least one of the algorithms in the package. 3) System dependent constants, which includes such things as the largest and smallest numbers, number of bits in the exponent and mantissa and so on. 4) The symbolic values of the error codes for the cmathe error package. The symbolic names, values and descriptions of each of those groups of constants appears below. /* * absolute constants */ PI 3.14159265358979324 pi HALF_PI 1.57079632679489662 pi/2 REC_PI 0.318309886183790672 1/pi RPIBY2 0.636619772367581343 reciprocal of pi/2 E 2.718281828459045235 e LOGBE2 0.69314718055994530942 log of 2 to base e LOGB2E 1.44269504088896341 log of e to base 2 ROOT_2 1.4142135623730950488 square root of 2 ROOT_3 1.732050807568877293 square root of 3 /* * checked octal constants */ PI_BY_2 1.57079632679489662 pi/2 PI_BY_4 0.785398163397448310 pi/4 LOGB10E 0.434294481903251828 log of e to base 10 ROOT_05 0.70710678118654752440 square root of 0.5 /* * system dependent constants */ TINY 0.293873587705571877E-38 smallest no = 2**-128 HUGE 0.170141183460469230E+39 largest no = 2**+127 LOG_HUGE 0.880296919311130543E+02 log of HUGE LOG_TINY -0.887228391116729997E+02 log of TINY MIN_EXP -128 minimum base 2 exponent MAX_EXP 127 maximum base 2 exponent MAXLONG 017777777777L largest long integer SIGFIGS 18 max no useful digits for dtoa TRIG_MAX 3.1415926535897932385e12 trig fn limit ROOT_EPS 0.372529029846191406e-8 2**-28 REC_HUGE 0.587747175411143754E-38 1 / HUGE Cmath Library Reference Manual Page 43 math.h cmath constants and definitions /* * error codes to communicate with cmathe */ FP_OPER 1 /* FPU op code error */ FP_ZDIV 2 /* FPU divide by zero */ FP_FTOI 3 /* FPU float to integer conv error */ FP_OFLO 4 /* FPU overflow */ FP_UFLO 5 /* FPU underflow */ FP_UDEF 6 /* FPU undefined variable (-0) */ FP_BIGI 7 /* Atof input too large */ FP_BADC 8 /* Bad character in atof input string */ FP_NESQ 9 /* Square root of negative number */ FP_LEXP 10 /* Exp argument too large */ FP_SEXP 11 /* Exp argument too small */ FP_NLOG 12 /* Log argument zero or negative */ FP_TANE 13 /* Argument of tan too large */ FP_TRIG 14 /* Argument of sin/cos too large */ FP_ATAN 15 /* Atan2 arguments both zero */ FP_COTE 16 /* Argument of cotan too small */ FP_ARSC 17 /* Bad argument for asin/acos */ FP_SINH 18 /* Argument of sinh too large */ FP_COSH 19 /* Argument of cosh too large */ FP_POWN 20 /* Negative argument in pow */ FP_POWO 21 /* Result of pow overflows */ FP_POWU 22 /* Result of pow underflows */ /* numbers of each type of error - to determine argument type */ FP_NFPU 6 /* No of FPU generated errors */ FP_NSTR 2 /* No of string argument errors */ FP_NMAR 14 /* No of math routine double arg errors */ AUTHOR: Hamish Ross DATE: 24-Feb-85 APPENDIX A CMATH LIBRARY INDEX The following is a keyword in context index to the cmath library. The entry in the left-hand column is the title of the routine in the cmath library documentation. fabs argument floating absolute value of floating arsico arc cosine function arsico arc sin function atan arctan function atan atan arctan function atan2 $$dtoa double to ascii conversion atof ascii to double conversion atan arctan function atan atan arctan function atan2 log10 logarithm function base 10 frexp exponent of a number to base 2 get mantissa and ldexp exponent of a number to base 2 load mantissa and math.h cmath constants and definitions atof ascii to double conversion $$dtoa double to ascii conversion sincos cosine function arsico arc cosine function cosh hyperbolic cosine function tancot cotangent function math.h cmath constants and definitions atof ascii to double conversion $$dtoa double to ascii conversion cmathe cmath error package frexp 2 get mantissa and exponent of a number to base ldexp 2 load mantissa and exponent of a number to base exp exponential function fabs floating argument floating absolute value of fabs absolute value of floating argument floating floor largest integer not greater than argument cosh hyperbolic cosine function sinh hyperbolic sine function tanh hyperbolic tanh function cminit cmath initialisation floor argument largest integer not greater than ceil argument smallest integer not less than Cmath Library Reference Manual Page A-2 Cmath Library Index floor than argument largest integer not greater ceil smallest integer not less than argument log logarithm function log10 logarithm function base 10 frexp number to base 2 get mantissa and exponent of a ldexp number to base 2 load mantissa and exponent of a modf modulus function pow x raised to power y powten table of powers of 10 sqrt square root function sincos sin function arsico arc sin function sinh hyperbolic sine function ceil than argument smallest integer not less sqrt square root function powten table of powers of 10 tancot tangent function tanh hyperbolic tanh function floor integer not greater than argument largest ceil smallest integer not less than argument