new from Kahan

SCCS-vsn: old/libm/libm/Makefile 4.9
SCCS-vsn: old/libm/libm/floor.c 4.2
SCCS-vsn: old/libm/libm/README 1.5
SCCS-vsn: old/libm/libm/lgamma.c 4.4
SCCS-vsn: include/math.h 4.6
SCCS-vsn: lib/libm/common_source/infnan.3 6.1
SCCS-vsn: lib/libm/common_source/exp.3 6.7
SCCS-vsn: lib/libm/common_source/floor.3 6.3
SCCS-vsn: lib/libm/common_source/lgamma.3 6.1
SCCS-vsn: lib/libm/common_source/j0.3 6.5
SCCS-vsn: lib/libm/common_source/hypot.3 6.4
SCCS-vsn: lib/libm/common_source/sin.3 6.5
SCCS-vsn: lib/libm/common_source/sinh.3 6.4
SCCS-vsn: lib/libm/common_source/math.3 6.5
SCCS-vsn: lib/libm/common_source/asinh.3 6.1
SCCS-vsn: lib/libm/common_source/erf.3 6.1
SCCS-vsn: lib/libm/common_source/ieee.3 6.1
SCCS-vsn: lib/libm/common_source/sqrt.3 6.1

diff --git a/usr/src/include/math.h b/usr/src/include/math.h
index 7b0ec9e..f483a37 100644 (file)
--- a/usr/src/include/math.h
+++ b/usr/src/include/math.h
@@ -1,10 +1,10 @@
-/*     math.h  4.5     %G%     */
+/*     math.h  4.6     %G%     */

 
 extern double asinh(), acosh(), atanh();
 extern double erf(), erfc();
 extern double exp(), expm1(), log(), log10(), log1p(), pow();
 
 extern double asinh(), acosh(), atanh();
 extern double erf(), erfc();
 extern double exp(), expm1(), log(), log10(), log1p(), pow();

-extern double fabs(), floor(), ceil();
-extern double gamma();
+extern double fabs(), floor(), ceil(), rint();
+extern double lgamma();

 extern double hypot(), cabs();
 extern double copysign(), drem(), logb(), scalb();
 extern int finite();
 extern double hypot(), cabs();
 extern double copysign(), drem(), logb(), scalb();
 extern int finite();

diff --git a/usr/src/lib/libm/common_source/asinh.3 b/usr/src/lib/libm/common_source/asinh.3
index 4894108..d5963d7 100644 (file)
--- a/usr/src/lib/libm/common_source/asinh.3
+++ b/usr/src/lib/libm/common_source/asinh.3
@@ -1,6 +1,12 @@
 From Prof. Kahan at UC at Berkeley
 From Prof. Kahan at UC at Berkeley

-.TH ASINH 3M "4 August 1985"
-.UC 4
+.\" Copyright (c) 1985 Regents of the University of California.
+.\" All rights reserved.  The Berkeley software License Agreement
+.\" specifies the terms and conditions for redistribution.
+.\"
+.\"    @(#)asinh.3     6.1 (Berkeley) %G%
+.\"
+.TH ASINH 3M  ""
+.UC 6

 .ds up \fIulp\fR
 .SH NAME
 asinh, acosh, atanh \- inverse hyperbolic functions
 .ds up \fIulp\fR
 .SH NAME
 asinh, acosh, atanh \- inverse hyperbolic functions

diff --git a/usr/src/lib/libm/common_source/erf.3 b/usr/src/lib/libm/common_source/erf.3
index cb7560a..34b7497 100644 (file)
--- a/usr/src/lib/libm/common_source/erf.3
+++ b/usr/src/lib/libm/common_source/erf.3
@@ -1,5 +1,12 @@
 From Prof. Kahan at UC at Berkeley
 From Prof. Kahan at UC at Berkeley

-.TH ERF 3M  "4 August 1985"
+.\" Copyright (c) 1985 Regents of the University of California.
+.\" All rights reserved.  The Berkeley software License Agreement
+.\" specifies the terms and conditions for redistribution.
+.\"
+.\"    @(#)erf.3       6.1 (Berkeley) %G%
+.\"
+.TH ERF 3M  ""
+.UC 6

 .SH NAME
 erf, erfc \- error functions
 .SH SYNOPSIS
 .SH NAME
 erf, erfc \- error functions
 .SH SYNOPSIS

diff --git a/usr/src/lib/libm/common_source/exp.3 b/usr/src/lib/libm/common_source/exp.3
index 3b480e3..bdb0eba 100644 (file)
--- a/usr/src/lib/libm/common_source/exp.3
+++ b/usr/src/lib/libm/common_source/exp.3
@@ -1,6 +1,13 @@
-.TH EXP 3M  "19 August 1985"
+.\" Copyright (c) 1985 Regents of the University of California.
+.\" All rights reserved.  The Berkeley software License Agreement
+.\" specifies the terms and conditions for redistribution.
+.\"
+.\"    @(#)exp.3       6.7 (Berkeley) %G%
+.\"
+.TH EXP 3M  ""

 .UC 4
 .ds nn \fINaN\fR
 .UC 4
 .ds nn \fINaN\fR

+.ds up \fIulp\fR

 .SH NAME
 exp, expm1, log, log10, log1p, pow \- exponential, logarithm, power
 .SH SYNOPSIS
 .SH NAME
 exp, expm1, log, log10, log1p, pow \- exponential, logarithm, power
 .SH SYNOPSIS


@@ -46,15 +53,18 @@ x\u\s8y\s10\d.
 intro(3M), infnan(3M)
 .SH ERROR (due to Roundoff etc.)
 exp(x), log(x), expm1(x) and log1p(x) are accurate to within 
 intro(3M), infnan(3M)
 .SH ERROR (due to Roundoff etc.)
 exp(x), log(x), expm1(x) and log1p(x) are accurate to within 

-an \fIulp\fR, and log10(x) to within about 2 \fIulp\fRs;
-an \fIulp\fR is one \fIU\fRnit in the \fIL\fRast \fIP\fRlace.
-The error in pow(x,y) is below about 2 \fIulp\fRs when its
+an \*(up, and log10(x) to within about 2 \*(ups;
+an \*(up is one \fIU\fRnit in the \fIL\fRast \fIP\fRlace.
+The error in pow(x,y) is below about 2 \*(ups when its

 magnitude is moderate, but increases as pow(x,y) approaches
 magnitude is moderate, but increases as pow(x,y) approaches

-the over/underflow thresholds until almost as many bits are
+the over/underflow thresholds until almost as many bits could be

 lost as are occupied by the floating\-point format's exponent
 field; that is 8 bits for VAX D and 11 bits for IEEE 754 Double.
 lost as are occupied by the floating\-point format's exponent
 field; that is 8 bits for VAX D and 11 bits for IEEE 754 Double.

+No such drastic loss has been exposed by testing; the worst
+errors observed have been below 20 \*(ups for VAX D,
+300 \*(ups for IEEE 754 Double.

 Moderate values of pow are accurate enough that pow(integer,integer)
 Moderate values of pow are accurate enough that pow(integer,integer)

-is exact until it is bigger than 2**56 on a VAX, 2**53 for IEEE.
+is exact until it is bigger than 2**56 on a VAX, 2**53 for IEEE 754.

 .SH DIAGNOSTICS
 Exp, expm1 and pow return the reserved operand on a VAX when the correct
 value would overflow, and they set \fIerrno\fR to ERANGE.
 .SH DIAGNOSTICS
 Exp, expm1 and pow return the reserved operand on a VAX when the correct
 value would overflow, and they set \fIerrno\fR to ERANGE.

diff --git a/usr/src/lib/libm/common_source/floor.3 b/usr/src/lib/libm/common_source/floor.3
index c66c3a5..1f20704 100644 (file)
--- a/usr/src/lib/libm/common_source/floor.3
+++ b/usr/src/lib/libm/common_source/floor.3
@@ -1,6 +1,14 @@
-.TH FLOOR 3M  "19 January 1983"
+.\" Copyright (c) 1985 Regents of the University of California.
+.\" All rights reserved.  The Berkeley software License Agreement
+.\" specifies the terms and conditions for redistribution.
+.\"
+.\"    @(#)floor.3     6.3 (Berkeley) %G%
+.\"
+.TH FLOOR 3M  ""
+.UC 4

 .SH NAME
 .SH NAME

-fabs, floor, ceil \- absolute value, floor, ceiling functions
+fabs, floor, ceil, rint \- absolute value, floor, ceiling, and
+round-to-nearest functions

 .SH SYNOPSIS
 .nf
 .B #include <math.h>
 .SH SYNOPSIS
 .nf
 .B #include <math.h>


@@ -13,19 +21,41 @@ fabs, floor, ceil \- absolute value, floor, ceiling functions
 .PP
 .B double fabs(x)
 .B double x;
 .PP
 .B double fabs(x)
 .B double x;

-.nf
+.PP
+.B double rint(x)
+.B double x;
+.fi

 .SH DESCRIPTION
 .SH DESCRIPTION

-.I Fabs
-returns the absolute value
-.RI | \|x\| |.
+Fabs returns the absolute value |\|x\||.

 .PP
 .PP

-.I Floor
-returns the largest integer not greater than
-.IR x .
+Floor returns the largest integer no greater than x.

 .PP
 .PP

-.I Ceil
-returns the smallest integer not less than
-.IR x .
+Ceil returns the smallest integer no less than x.
+.PP
+Rint returns the integer (represented as a double precision number)
+nearest x in the direction of the prevailing rounding mode.

 .SH SEE ALSO
 abs(3),
 .SH SEE ALSO
 abs(3),

+ieee(3M),

 intro(3M)
 intro(3M)

+.SH NOTES
+On a VAX, rint(x) is equivalent to adding half to the magnitude
+and then rounding towards zero.
+.PP
+In the default rounding mode, to nearest,
+on a machine that conforms to IEEE 754,
+rint(x) is the integer nearest x with the additional stipulation
+that if |rint(x)\-x|=1/2 then rint(x) is even.
+Other rounding modes can make rint act like floor, or like ceil,
+or round towards zero.
+.PP
+Another way to obtain an integer near x is to declare (in C)
+.RS
+double x;\0\0\0\0 int k;\0\0\0\0k\0=\0x;
+.RE
+Most C compilers round x towards 0 to get the integer k, but
+some do otherwise.
+If in doubt, use floor, ceil, or rint first, whichever you intend.
+Also note that, if x is larger than k can accommodate, the value of
+k and the presence or absence of an integer overflow are hard to
+predict.

diff --git a/usr/src/lib/libm/common_source/hypot.3 b/usr/src/lib/libm/common_source/hypot.3
index df7e6c7..85a3d7c 100644 (file)
--- a/usr/src/lib/libm/common_source/hypot.3
+++ b/usr/src/lib/libm/common_source/hypot.3
@@ -1,4 +1,10 @@
-.TH HYPOT 3M  "7 August 1985"
+.\" Copyright (c) 1985 Regents of the University of California.
+.\" All rights reserved.  The Berkeley software License Agreement
+.\" specifies the terms and conditions for redistribution.
+.\"
+.\"    @(#)hypot.3     6.4 (Berkeley) %G%
+.\"
+.TH HYPOT 3M  ""

 .UC 4
 .ds up \fIulp\fR
 .ds nn \fINaN\fR
 .UC 4
 .ds up \fIulp\fR
 .ds nn \fINaN\fR

diff --git a/usr/src/lib/libm/common_source/ieee.3 b/usr/src/lib/libm/common_source/ieee.3
index ff384cc..4db2f55 100644 (file)
--- a/usr/src/lib/libm/common_source/ieee.3
+++ b/usr/src/lib/libm/common_source/ieee.3
@@ -1,6 +1,12 @@
 From Prof. Kahan at UC at Berkeley
 From Prof. Kahan at UC at Berkeley

-.TH IEEE 3M  "7 August 1985"
-.UC 4
+.\" Copyright (c) 1985 Regents of the University of California.
+.\" All rights reserved.  The Berkeley software License Agreement
+.\" specifies the terms and conditions for redistribution.
+.\"
+.\"    @(#)ieee.3      6.1 (Berkeley) %G%
+.\"
+.TH IEEE 3M  ""
+.UC 6

 .ds nn \fINaN\fR
 .SH NAME
 copysign, drem, finite, logb, scalb \- copysign, remainder,
 .ds nn \fINaN\fR
 .SH NAME
 copysign, drem, finite, logb, scalb \- copysign, remainder,


@@ -65,7 +71,7 @@ or x lies between 0 and the Underflow Threshold; see below under "BUGS".
 Scalb(x,n) = x\(**(2**n) computed, for integer n, without first computing
 2**n.
 .SH SEE ALSO
 Scalb(x,n) = x\(**(2**n) computed, for integer n, without first computing
 2**n.
 .SH SEE ALSO

-intro(3M), infnan(3M)
+floor(3M), intro(3M), infnan(3M)

 .SH DIAGNOSTICS
 IEEE 754 defines drem(x,0) and
 .if n \
 .SH DIAGNOSTICS
 IEEE 754 defines drem(x,0) and
 .if n \

diff --git a/usr/src/lib/libm/common_source/infnan.3 b/usr/src/lib/libm/common_source/infnan.3
index 13fe93e..e93236d 100644 (file)
--- a/usr/src/lib/libm/common_source/infnan.3
+++ b/usr/src/lib/libm/common_source/infnan.3
@@ -1,13 +1,14 @@
-.TH INFNAN 3M  "20 August 1985"
-.UC 4
+.\" Copyright (c) 1985 Regents of the University of California.
+.\" All rights reserved.  The Berkeley software License Agreement
+.\" specifies the terms and conditions for redistribution.
+.\"
+.\"    @(#)infnan.3    6.1 (Berkeley) %G%
+.\"
+.TH INFNAN 3M  ""
+.UC 6

 .ds nn \fINaN\fR
 .SH NAME
 .ds nn \fINaN\fR
 .SH NAME

-.nf
-.ta \w'infnan \- 'u+1n
-infnan \-      signals invalid floating\-point operations on a VAX
-       TEMPORARY! ... Subject to Change.
-.ta
-.fi
+infnan \- signals invalid floating-point operations on a VAX ... TEMPORARY

 .SH SYNOPSIS
 .nf
 .B #include <math.h>
 .SH SYNOPSIS
 .nf
 .B #include <math.h>


@@ -54,8 +55,13 @@ implement that suggestion follows.
 .ta \w'Div\-by\-0'u+2n +\w'+Infinity'u+1n +\w'+ERANGE'u+1n +\w'ERANGE or EDOM'u+4n +\w'+HUGE'u+1n
 IEEE   IEEE
 Signal Default \fIiarg\fR      \fIerrno\fR     \fIinfnan\fR
 .ta \w'Div\-by\-0'u+2n +\w'+Infinity'u+1n +\w'+ERANGE'u+1n +\w'ERANGE or EDOM'u+4n +\w'+HUGE'u+1n
 IEEE   IEEE
 Signal Default \fIiarg\fR      \fIerrno\fR     \fIinfnan\fR

-.if t \
-\l'4i'
+.if t \{\
+.sp -0.5
+.ta \w'Div\-by\-0'u+2n+\w'+Infinity'u+1n+\w'+ERANGE'u+1n+\w'ERANGE or EDOM'u+4n+\w'+HUGE'u+1n
+.tc \(ru
+       
+.ta \w'Div\-by\-0'u+2n +\w'+Infinity'u+1n +\w'+ERANGE'u+1n +\w'ERANGE or EDOM'u+4n +\w'+HUGE'u+1n
+.tc \}

 .if n \
 \l'5i'
 Invalid        \*(nn   EDOM    EDOM    0
 .if n \
 \l'5i'
 Invalid        \*(nn   EDOM    EDOM    0

diff --git a/usr/src/lib/libm/common_source/j0.3 b/usr/src/lib/libm/common_source/j0.3
index f9f99a1..28173c3 100644 (file)
--- a/usr/src/lib/libm/common_source/j0.3
+++ b/usr/src/lib/libm/common_source/j0.3
@@ -1,4 +1,11 @@
-.TH J0 3M  "19 January 1983"
+.\" Copyright (c) 1985 Regents of the University of California.
+.\" All rights reserved.  The Berkeley software License Agreement
+.\" specifies the terms and conditions for redistribution.
+.\"
+.\"    @(#)j0.3        6.5 (Berkeley) %G%
+.\"
+.TH J0 3M  ""
+.UC 4

 .SH NAME
 j0, j1, jn, y0, y1, yn \- bessel functions
 .SH SYNOPSIS
 .SH NAME
 j0, j1, jn, y0, y1, yn \- bessel functions
 .SH SYNOPSIS

diff --git a/usr/src/lib/libm/common_source/lgamma.3 b/usr/src/lib/libm/common_source/lgamma.3
index a3c4ca9..b87e10d 100644 (file)
--- a/usr/src/lib/libm/common_source/lgamma.3
+++ b/usr/src/lib/libm/common_source/lgamma.3
@@ -1,30 +1,87 @@
-.TH GAMMA 3M  "19 January 1983"
+.\" Copyright (c) 1985 Regents of the University of California.
+.\" All rights reserved.  The Berkeley software License Agreement
+.\" specifies the terms and conditions for redistribution.
+.\"
+.\"    @(#)lgamma.3    6.1 (Berkeley) %G%
+.\"
+.TH LGAMMA 3M  ""
+.UC 6

 .SH NAME
 .SH NAME

-gamma \- log gamma function
+lgamma \- log gamma function

 .SH SYNOPSIS
 .nf
 .B #include <math.h>
 .PP
 .SH SYNOPSIS
 .nf
 .B #include <math.h>
 .PP

-.B double gamma(x)
+.B double lgamma(x)

 .B double x;
 .fi
 .SH DESCRIPTION
 .B double x;
 .fi
 .SH DESCRIPTION

-.I Gamma
-returns ln \||\|\(*G(\||\|\fIx\fR\||\|)\||\|.
-The sign of \(*G(\||\|\fIx\fR\||\|) is returned in the external integer
-.IR signgam .
-The following C program might be used to calculate \(*G:
+.nf
+.ta \w'Lgamma returns ln\||\(*G(x)| where'u+1n +1.7i
+.if t \{\
+Lgamma returns ln\||\(*G(x)| where     \(*G(x) = \(is\d\s8\z0\s10\u\u\s8\(if\s10\d t\u\s8x\-1\s10\d e\u\s8\-t\s10\d dt for x > 0 and
+.br
+       \(*G(x) = \(*p/(\(*G(1\-x)\|sin(\(*px)) for x < 1.  \}
+.if n \
+Lgamma returns ln\||\(*G(x)|.
+.ta
+.fi
+.PP
+The external integer signgam returns the sign of
+\(*G(x) .
+.SH IDIOSYNCRASIES
+Do \fBnot\fR use the expression signgam\(**exp(lgamma(x))
+to compute g := \(*G(x).  Instead use a program like this (in C):
+.RS
+lg = lgamma(x); g = signgam\(**exp(lg);
+.RE
+.PP
+Only after lgamma has returned can signgam be correct.
+Note too that \(*G(x) must overflow when x is large enough,
+underflow when \-x is large enough, and spawn a division by zero
+when x is a nonpositive integer.
+.PP
+Only in the UNIX math library for C was the name gamma ever attached
+to ln\(*G.  Elsewhere, for instance in IBM's FORTRAN library, the name
+GAMMA belongs to \(*G and the name ALGAMA to ln\(*G in single precision;
+in double the names are DGAMMA and DLGAMA.  Why should C be different?
+.PP
+Archaeological records suggest that C's gamma originally delivered
+ln(\(*G(|x|)).  Later, the program gamma was changed to
+cope with negative arguments x in a more conventional way, but
+the documentation did not reflect that change correctly.  The most
+recent change corrects inaccurate values when x is almost a
+negative integer, and lets \(*G(x) be computed without
+conditional expressions.  Programmers should not assume that
+lgamma has settled down.
+.PP
+At some time in the future, the name \fIgamma\fR will be rehabilitated
+and used for the gamma function, just as is done in FORTRAN.
+The reason for this is not so much compatibility with FORTRAN as a
+desire to achieve greater speed for smaller values of |x| and greater
+accuracy for larger values.

 .PP
 .PP

+Meanwhile, programmers who have to use the name \fIgamma\fR in its former
+sense, for what is now \fIlgamma\fR, have two choices:
+.IP 1) \w'1)\0'u
+Use the old math library, \fIlibom\fR.
+.IP 2) \w'1)\0'u
+Add the following program to your others:
+.RS

 .nf
 .nf

-       y = gamma(x);
-       if (y > 88.0)
-               error();
-       y = exp(y);
-       if(signgam)
-               y = \-y;
+\fB#include <math.h>
+double gamma(x)
+double x;
+{
+.RS
+\fBreturn (lgamma(x));
+.RE
+}\fR
+.RE

 .fi
 .fi

+.SH SEE ALSO
+intro(3M), infnan(3M)

 .SH DIAGNOSTICS
 .SH DIAGNOSTICS

-The reserved operand is returned on a VAX for negative integer
-arguments, \fIerrno\fR is set to ERANGE.
-.SH BUGS
-There should be a positive indication of error.
+The reserved operand is returned on a VAX for negative integer arguments,
+\fIerrno\fR is set to ERANGE; for very large arguments over/underflows will
+occur inside the program lgamma.

diff --git a/usr/src/lib/libm/common_source/math.3 b/usr/src/lib/libm/common_source/math.3
index 504429b..d538608 100644 (file)
--- a/usr/src/lib/libm/common_source/math.3
+++ b/usr/src/lib/libm/common_source/math.3
@@ -1,4 +1,10 @@
-.TH INTRO 3M "19 August 1985"
+.\" Copyright (c) 1985 Regents of the University of California.
+.\" All rights reserved.  The Berkeley software License Agreement
+.\" specifies the terms and conditions for redistribution.
+.\"
+.\"    @(#)math.3      6.5 (Berkeley) %G%
+.\"
+.TH INTRO 3M ""

 .UC 4
 .ds up \fIulp\fR
 .ds nn \fINaN\fR
 .UC 4
 .ds up \fIulp\fR
 .ds nn \fINaN\fR


@@ -45,17 +51,18 @@ exp exp.3m  exponential     1
 expm1  exp.3m  exp(x)\-1       1
 fabs   floor.3m        absolute value  0
 floor  floor.3m        integer no greater than 0
 expm1  exp.3m  exp(x)\-1       1
 fabs   floor.3m        absolute value  0
 floor  floor.3m        integer no greater than 0

-gamma  gamma.3m        log gamma function

 hypot  hypot.3m        Euclidean distance      1
 infnan infnan.3m       signals exceptions
 j0     j0.3m   bessel function ???
 j1     j0.3m   bessel function ???
 jn     j0.3m   bessel function ???
 hypot  hypot.3m        Euclidean distance      1
 infnan infnan.3m       signals exceptions
 j0     j0.3m   bessel function ???
 j1     j0.3m   bessel function ???
 jn     j0.3m   bessel function ???

+lgamma lgamma.3m       log gamma function; (formerly gamma.3m)

 log    exp.3m  natural logarithm       1
 logb   ieee.3m exponent extraction     0
 log10  exp.3m  logarithm to base 10    3
 log1p  exp.3m  log(1+x)        1
 pow    exp.3m  exponential x**y        60\-500
 log    exp.3m  natural logarithm       1
 logb   ieee.3m exponent extraction     0
 log10  exp.3m  logarithm to base 10    3
 log1p  exp.3m  log(1+x)        1
 pow    exp.3m  exponential x**y        60\-500

+rint   floor.3m        round to nearest integer        0

 scalb  ieee.3m exponent adjustment     0
 sin    sin.3m  trigonometric function  1
 sinh   sinh.3m hyperbolic function     3
 scalb  ieee.3m exponent adjustment     0
 sin    sin.3m  trigonometric function  1
 sinh   sinh.3m hyperbolic function     3


@@ -114,7 +121,17 @@ Properties of D_floating\-point:
 .RS
 Wordsize: 64 bits, 8 bytes.  Radix: Binary.
 .br
 .RS
 Wordsize: 64 bits, 8 bytes.  Radix: Binary.
 .br

-Precision: 56 bits; equivalent to about 17 significant decimals.
+Precision: 56
+.if n \
+sig.
+.if t \
+significant
+bits, roughly like 17
+.if n \
+sig.
+.if t \
+significant
+decimals.

 .RS
 If x and x' are consecutive positive D_floating\-point
 numbers (they differ by 1 \*(up), then
 .RS
 If x and x' are consecutive positive D_floating\-point
 numbers (they differ by 1 \*(up), then


@@ -230,13 +247,13 @@ produced by a host of manufacturers, among them ...
 .nf
 .ta 0.5i +\w'Intel i8070, i80287'u+6n
        Intel i8087, i80287     National Semiconductor  32081
 .nf
 .ta 0.5i +\w'Intel i8070, i80287'u+6n
        Intel i8087, i80287     National Semiconductor  32081

-       Motorola 68881  Weitek WTL-1032, ... , -1065
-       Zilog Z8070     Western Electric (AT&T) 32106.
+       Motorola 68881  Weitek WTL-1032, ... , -1165
+       Zilog Z8070     Western Electric (AT&T) WE32106.

 .ta
 .fi
 Other implementations range from software, done thoroughly
 in the Apple Macintosh, through VLSI in the Hewlett\-Packard
 .ta
 .fi
 Other implementations range from software, done thoroughly
 in the Apple Macintosh, through VLSI in the Hewlett\-Packard

-9000 series, to the ELXSI 6400 running at 3 Megaflops.
+9000 series, to the ELXSI 6400 running ECL at 3 Megaflops.

 Several other companies have adopted the formats
 of IEEE 754 without, alas, adhering to the standard's way
 of handling rounding and exceptions like over/underflow.
 Several other companies have adopted the formats
 of IEEE 754 without, alas, adhering to the standard's way
 of handling rounding and exceptions like over/underflow.


@@ -248,13 +265,23 @@ nobody has volunteered to do that yet.
 .PP
 The codes in 4.3 BSD's \fIlibm\fR for machines that conform to
 IEEE 754 are intended primarily for the National Semi. 32081
 .PP
 The codes in 4.3 BSD's \fIlibm\fR for machines that conform to
 IEEE 754 are intended primarily for the National Semi. 32081

-and Weitek 1065.  To use these codes with the Intel or Zilog
+and WTL 1164/65.  To use these codes with the Intel or Zilog

 chips, or with the Apple Macintosh or ELXSI 6400, is to
 chips, or with the Apple Macintosh or ELXSI 6400, is to

-forego the use of better codes available (for a price) from
+forego the use of better codes provided (perhaps freely) by

 those companies and designed by some of the authors of the
 those companies and designed by some of the authors of the

-codes above.  The Motorola 68881 has all the \fIelementary\fR
-functions above except \fIpow\fR implemented on chip already,
-and faster and more accurate.  The main virtue of 4.3 BSD's
+codes above.
+Except for \fIatan\fR, \fIcabs\fR, \fIcbrt\fR, \fIerf\fR,
+\fIerfc\fR, \fIhypot\fR, \fIj0\-jn\fR, \fIlgamma\fR, \fIpow\fR
+and \fIy0\-yn\fR,
+the Motorola 68881 has all the functions in \fIlibm\fR on chip,
+and faster and more accurate;
+it, Apple, the i8087, Z8070 and WE32106 all use 64
+.if n \
+sig.
+.if t \
+significant
+bits.
+The main virtue of 4.3 BSD's

 \fIlibm\fR codes is that they are intended for the public domain;
 they may be copied freely provided their provenance is always
 acknowledged, and provided users assist the authors in their
 \fIlibm\fR codes is that they are intended for the public domain;
 they may be copied freely provided their provenance is always
 acknowledged, and provided users assist the authors in their


@@ -266,7 +293,17 @@ Properties of IEEE 754 Double\-Precision:
 .RS
 Wordsize: 64 bits, 8 bytes.  Radix: Binary.
 .br
 .RS
 Wordsize: 64 bits, 8 bytes.  Radix: Binary.
 .br

-Precision: 53 bits; equivalent to about 16 significant decimals.
+Precision: 53
+.if n \
+sig.
+.if t \
+significant
+bits, roughly like 16
+.if n \
+sig.
+.if t \
+significant
+decimals.

 .RS
 If x and x' are consecutive positive Double\-Precision
 numbers (they differ by 1 \*(up), then
 .RS
 If x and x' are consecutive positive Double\-Precision
 numbers (they differ by 1 \*(up), then


@@ -308,8 +345,8 @@ finite x = y then
 .If
 is signed.
 .RS
 .If
 is signed.
 .RS

-it persists after addition of itself
-or of any finite number.  Its sign transforms
+it persists when added to itself
+or to any finite number.  Its sign transforms

 correctly through multiplication and division, and
 .If (finite)/\(+- \0=\0\(+-0
 (nonzero)/0 =
 correctly through multiplication and division, and
 .If (finite)/\(+- \0=\0\(+-0
 (nonzero)/0 =


@@ -498,9 +535,9 @@ the data supplied to that function.
 Any exception signaled should be identified with that
 function rather than with one of its subroutines.
 .IP iii) \w'iii)'u+2n
 Any exception signaled should be identified with that
 function rather than with one of its subroutines.
 .IP iii) \w'iii)'u+2n

-The internal behavior of an atomic function cannot be
-disrupted when a calling program changes from one
-to another of the five or so ways of handling
+The internal behavior of an atomic function should not
+be disrupted when a calling program changes from
+one to another of the five or so ways of handling

 exceptions listed above, although the definition
 of the function may be correlated intentionally
 with exception handling.
 exceptions listed above, although the definition
 of the function may be correlated intentionally
 with exception handling.


@@ -520,7 +557,7 @@ Over/Underflow
 .RS
 when a result, if properly computed, might have lain barely within range, and
 .RE
 .RS
 when a result, if properly computed, might have lain barely within range, and
 .RE

-Inexact in pow(x,y)
+Inexact in \fIcabs\fR, \fIcbrt\fR, \fIhypot\fR, \fIlog10\fR and \fIpow\fR

 .RS
 when it happens to be exact thanks to fortuitous cancellation of errors.
 .RE
 .RS
 when it happens to be exact thanks to fortuitous cancellation of errors.
 .RE


@@ -552,7 +589,7 @@ greater range or precision would be needed to represent the exact result.
 When signals are appropriate, they are emitted by certain
 operations within the codes, so a subroutine\-trace may be
 needed to identify the function with its signal in case
 When signals are appropriate, they are emitted by certain
 operations within the codes, so a subroutine\-trace may be
 needed to identify the function with its signal in case

-method 6) above is in use.  And the codes all take the
+method 5) above is in use.  And the codes all take the

 IEEE 754 defaults for granted; this means that a decision to
 trap all divisions by zero could disrupt a code that would
 otherwise get correct results despite division by zero.
 IEEE 754 defaults for granted; this means that a decision to
 trap all divisions by zero could disrupt a code that would
 otherwise get correct results despite division by zero.


@@ -568,5 +605,5 @@ Articles in the IEEE magazine COMPUTER vol. 14 no. 3 (Mar.
 Oct. 1979, may be helpful although they pertain to
 superseded drafts of the standard.
 .SH AUTHOR
 Oct. 1979, may be helpful although they pertain to
 superseded drafts of the standard.
 .SH AUTHOR

-W. Kahan, with the help of Alex Z\-S. Liu, Stuart I. McDonald,
+W. Kahan, with the help of Z\-S. Alex Liu, Stuart I. McDonald,

 Dr. Kwok\-Choi Ng, Peter Tang.
 Dr. Kwok\-Choi Ng, Peter Tang.

diff --git a/usr/src/lib/libm/common_source/sin.3 b/usr/src/lib/libm/common_source/sin.3
index f35e0aa..5232c8b 100644 (file)
--- a/usr/src/lib/libm/common_source/sin.3
+++ b/usr/src/lib/libm/common_source/sin.3
@@ -1,4 +1,10 @@
-.TH SIN 3M  "6 August 1985"
+.\" Copyright (c) 1985 Regents of the University of California.
+.\" All rights reserved.  The Berkeley software License Agreement
+.\" specifies the terms and conditions for redistribution.
+.\"
+.\"    @(#)sin.3       6.5 (Berkeley) %G%
+.\"
+.TH SIN 3M  ""

 .UC 4
 .de Pi         \" PI stuff sign
 .if n \\
 .UC 4
 .de Pi         \" PI stuff sign
 .if n \\


@@ -183,8 +189,9 @@ sin(x)**2+cos(x)**2\0=\01
 .if t \
 sin\u\s62\s10\d(x)+cos\u\s62\s10\d(x)\0=\01
 and sin(2x)\0=\02\|sin(x)cos(x) to within a few \*(ups no matter how big
 .if t \
 sin\u\s62\s10\d(x)+cos\u\s62\s10\d(x)\0=\01
 and sin(2x)\0=\02\|sin(x)cos(x) to within a few \*(ups no matter how big

-x may be.  Therefore scientific and engineering calculations
-are most unlikely to be affected by the difference between P and
-.Pi .
+x may be.  Therefore the difference between P and
+.Pi
+is most unlikely to affect scientific and engineering computations.

 .SH AUTHOR
 .SH AUTHOR

-Robert P. Corbett, W. Kahan, Stuart I. McDonald, Kwok\-Choi\0Ng
+Robert P. Corbett, W. Kahan, Stuart\0I.\0McDonald, Peter\0Tang and,
+for the codes for IEEE 754, Dr. Kwok\-Choi\0Ng.

diff --git a/usr/src/lib/libm/common_source/sinh.3 b/usr/src/lib/libm/common_source/sinh.3
index af4fe77..db5c171 100644 (file)
--- a/usr/src/lib/libm/common_source/sinh.3
+++ b/usr/src/lib/libm/common_source/sinh.3
@@ -1,4 +1,10 @@
-.TH SINH 3M "7 August 1985"
+.\" Copyright (c) 1985 Regents of the University of California.
+.\" All rights reserved.  The Berkeley software License Agreement
+.\" specifies the terms and conditions for redistribution.
+.\"
+.\"    @(#)sinh.3      6.4 (Berkeley) %G%
+.\"
+.TH SINH 3M  ""

 .UC 4
 .SH NAME
 sinh, cosh, tanh \- hyperbolic functions
 .UC 4
 .SH NAME
 sinh, cosh, tanh \- hyperbolic functions

diff --git a/usr/src/lib/libm/common_source/sqrt.3 b/usr/src/lib/libm/common_source/sqrt.3
index 153e87f..2e8b223 100644 (file)
--- a/usr/src/lib/libm/common_source/sqrt.3
+++ b/usr/src/lib/libm/common_source/sqrt.3
@@ -1,5 +1,11 @@
-.TH SQRT 3M  "7 August 1985"
-.UC 4
+.\" Copyright (c) 1985 Regents of the University of California.
+.\" All rights reserved.  The Berkeley software License Agreement
+.\" specifies the terms and conditions for redistribution.
+.\"
+.\"    @(#)sqrt.3      6.1 (Berkeley) %G%
+.\"
+.TH SQRT 3M  ""
+.UC 6

 .ds up \fIulp\fR
 .SH NAME
 cbrt, sqrt \- cube root, square root
 .ds up \fIulp\fR
 .SH NAME
 cbrt, sqrt \- cube root, square root

diff --git a/usr/src/old/libm/libm/Makefile b/usr/src/old/libm/libm/Makefile
index 0e44a6c..619d78b 100644 (file)
--- a/usr/src/old/libm/libm/Makefile
+++ b/usr/src/old/libm/libm/Makefile
@@ -1,7 +1,7 @@
 #
 #

-#      @(#)Makefile    4.8     %G%
+#      @(#)Makefile    4.9     %G%

 #
 #

-SCCSID = "@(#)Makefile 4.8 %G%"
+SCCSID = "@(#)Makefile 4.9 %G%"

 #
 # This high quality math library is intended to run on either a VAX in
 # D_floating format or a machine that conforms to the IEEE standard 754
 #
 # This high quality math library is intended to run on either a VAX in
 # D_floating format or a machine that conforms to the IEEE standard 754


@@ -35,11 +35,11 @@ DESTDIR=
 # actually there are more under ${MACH}/ subdirectory.
 #
 SRCS = acosh.c asincos.c asinh.c atan.c atanh.c cosh.c erf.c \
 # actually there are more under ${MACH}/ subdirectory.
 #
 SRCS = acosh.c asincos.c asinh.c atan.c atanh.c cosh.c erf.c \

-       exp.c exp__E.c expm1.c floor.c gamma.c j0.c j1.c jn.c \
+       exp.c exp__E.c expm1.c floor.c lgamma.c j0.c j1.c jn.c \

        log.c log10.c log1p.c log__L.c pow.c sinh.c tanh.c
 
 FILES =        acosh.o asincos.o asinh.o atan.o atanh.o cosh.o erf.o \
        log.c log10.c log1p.c log__L.c pow.c sinh.c tanh.c
 
 FILES =        acosh.o asincos.o asinh.o atan.o atanh.o cosh.o erf.o \

-       exp.o exp__E.o expm1.o floor.o gamma.o j0.o j1.o jn.o \
+       exp.o exp__E.o expm1.o floor.o lgamma.o j0.o j1.o jn.o \

        log.o log10.o log1p.o log__L.o pow.o sinh.o tanh.o 
 
 TAGSFILE=tags
        log.o log10.o log1p.o log__L.o pow.o sinh.o tanh.o 
 
 TAGSFILE=tags

diff --git a/usr/src/old/libm/libm/README b/usr/src/old/libm/libm/README
index 4026378..4b74d00 100644 (file)
--- a/usr/src/old/libm/libm/README
+++ b/usr/src/old/libm/libm/README
@@ -1,99 +1,113 @@
+***************************************************************************
+*                                                                         * 
+* Copyright (c) 1985 Regents of the University of California.             *
+*                                                                         * 
+* Use and reproduction of this software are granted  in  accordance  with *
+* the terms and conditions specified in  the  Berkeley  Software  License *
+* Agreement (in particular, this entails acknowledgement of the programs' *
+* source, and inclusion of this notice) with the additional understanding *
+* that  all  recipients  should regard themselves as participants  in  an *
+* ongoing  research  project and hence should  feel  obligated  to report *
+* their  experiences (good or bad) with these elementary function  codes, *
+* using "sendbug 4bsd-bugs@BERKELEY", to the authors.                     *
+*                                                                         *
+* K.C. Ng, with Z-S. Alex Liu, S. McDonald, P. Tang, W. Kahan.            *
+* Revised on 5/10/85, 5/13/85, 6/14/85, 8/20/85, 8/27/85.                 *
+*                                                                         *
+***************************************************************************

 
 

--1.  The machine-independent Version 7 math library found in 4.2BSD
-     is now /usr/lib/libom.a.  To compile with these routines use -lom.
+/*     @(#)README      1.5 (Berkeley) %G%      */

 
 

-K.C. Ng, March 7, 1985, with Z-S. Liu,  S. McDonald,  P. Tang,  W. Kahan.
-Revised on 5/10/85, 5/13/85, 6/14/85, 8/20/85.
+NB.  The machine-independent Version 7 math library found in 4.2BSD
+     is now /usr/lib/libom.a.  To compile with those routines use -lom.

 
 ******************************************************************************
 *  This is a description of the upgraded elementary functions (listed in 1). *
 *  Bessel functions (j0, j1, jn, y0, y1, yn), floor, and fabs passed over    *
 
 ******************************************************************************
 *  This is a description of the upgraded elementary functions (listed in 1). *
 *  Bessel functions (j0, j1, jn, y0, y1, yn), floor, and fabs passed over    *

-*  from 4.2BSD without change. Three lines that contain "errno in erf.c      *
-*  (error function erf, erfc) have been deleted to prevent overwriting the   *
-*  system "errno".                                                          *
+*  from 4.2BSD without change except perhaps for the way floating point      *
+*  exception is signaled on a VAX.  Three lines that contain "errno in erf.c *
+*  (error function erf, erfc) have been deleted to prevent overriding the    *
+*  system "errno".                                                           *

 ******************************************************************************
 
 0. Total number of files: 40
 
 ******************************************************************************
 
 0. Total number of files: 40
 

-  ( unchanges files from 4.2BSD, total 4 )
-    floor.c    gamma.c
-
-  ( new files for 4.3BSD, total 36 )
-    IEEE/Makefile      VAX/cabs.s      asinh.c         jn.c
-    IEEE/atan2.c       VAX/cbrt.s      atan.c          log.c
-    IEEE/cabs.c                VAX/infnan.s    cosh.c          log__L.c
-    IEEE/cbrt.c                VAX/support.s   erf.c           log1p.c
-    IEEE/trig.c                VAX/sincos.s    exp.c           pow.c
-    Makefile           VAX/tan.s       exp__E.c        sinh.c
-    README             VAX/argred.s    expm1.c         tanh.c
-    VAX/Makefile       acosh.c         j0.c
-    VAX/atan2.s                asincos.c       j1.c
+        IEEE/Makefile   VAX/Makefile    VAX/support.s   erf.c       lgamma.c
+        IEEE/atan2.c    VAX/argred.s    VAX/tan.s       exp.c       log.c
+        IEEE/cabs.c     VAX/atan2.s     acosh.c         exp__E.c    log10.c
+        IEEE/cbrt.c     VAX/cabs.s      asincos.c       expm1.c     log1p.c
+        IEEE/support.c  VAX/cbrt.s      asinh.c         floor.c     log__L.c
+        IEEE/trig.c     VAX/infnan.s    atan.c          j0.c        pow.c
+        Makefile        VAX/sincos.s    atanh.c         j1.c        sinh.c
+        README          VAX/sqrt.s      cosh.c          jn.c        tanh.c

 
 1. Functions implemented :
     (A). Standard elementary functions (total 22) :
 
 1. Functions implemented :
     (A). Standard elementary functions (total 22) :

-       acos(x)                 ...in file  asincos.c 
-       asin(x)                 ...in file  asincos.c
-       atan(x)                 ...in file  atan.c
-       atan2(x,y)              ...in files IEEE/atan2.c, VAX/atan2.s
-       sin(x)                  ...in files IEEE/trig.c , VAX/sincos.s
-       cos(x)                  ...in files IEEE/trig.c , VAX/sincos.s
-       tan(x)                  ...in files IEEE/trig.c , VAX/tan.s
-       cabs(x,y)               ...in files IEEE/cabs.c , VAX/cabs.s
-       hypot(x,y)              ...in files IEEE/cabs.c , VAX/cabs.s
-       cbrt(x)                 ...in files IEEE/cbrt.c , VAX/cbrt.s
-       exp(x)                  ...in file  exp.c
-       expm1(x):=exp(x)-1      ...in file  expm1.c
-       log(x)                  ...in file  log.c
-       log10(x)                ...in file  log10.c
-       log1p(x):=log(1+x)      ...in file  log1p.c
-       pow(x,y)                ...in file  pow.c
-       sinh(x)                 ...in file  sinh.c
-       cosh(x)                 ...in file  cosh.c
-       tanh(x)                 ...in file  cosh.c
-       asinh(x)                ...in file  asinh.c
-       acosh(x)                ...in file  acosh.c
-       atanh(x)                ...in file  atanh.c
-                       
+        acos(x)                 ...in file  asincos.c 
+        asin(x)                 ...in file  asincos.c
+        atan(x)                 ...in file  atan.c
+        atan2(x,y)              ...in files IEEE/atan2.c, VAX/atan2.s
+        sin(x)                  ...in files IEEE/trig.c , VAX/sincos.s
+        cos(x)                  ...in files IEEE/trig.c , VAX/sincos.s
+        tan(x)                  ...in files IEEE/trig.c , VAX/tan.s
+        cabs(x,y)               ...in files IEEE/cabs.c , VAX/cabs.s
+        hypot(x,y)              ...in files IEEE/cabs.c , VAX/cabs.s
+        cbrt(x)                 ...in files IEEE/cbrt.c , VAX/cbrt.s
+        exp(x)                  ...in file  exp.c
+        expm1(x):=exp(x)-1      ...in file  expm1.c
+        log(x)                  ...in file  log.c
+        log10(x)                ...in file  log10.c
+        log1p(x):=log(1+x)      ...in file  log1p.c
+        pow(x,y)                ...in file  pow.c
+        sinh(x)                 ...in file  sinh.c
+        cosh(x)                 ...in file  cosh.c
+        tanh(x)                 ...in file  cosh.c
+        asinh(x)                ...in file  asinh.c
+        acosh(x)                ...in file  acosh.c
+        atanh(x)                ...in file  atanh.c
+                        

     (B). Kernel functions :
     (B). Kernel functions :

-       exp__E(x,c) ...in file exp__E.c, used by expm1/exp/pow/cosh
-       log__L(s)   ...in file log__L.c, used by log1p/log/pow
-       libm$argred ...in file VAX/argred.s, used by VAX/tan.s and VAX/sincos.s
+        exp__E(x,c) ...in file exp__E.c, used by expm1/exp/pow/cosh
+        log__L(s)   ...in file log__L.c, used by log1p/log/pow
+        libm$argred ...in file VAX/argred.s, used by VAX/tan.s and VAX/sincos.s

 
     (C). System supported functions :
 
     (C). System supported functions :

-       sqrt()      ...in files IEEE/support.c , VAX/sqrt.s
-       drem()      ...in files IEEE/support.c , VAX/support.s
-       finite()    ...in files IEEE/support.c , VAX/support.s
-       logb()      ...in files IEEE/support.c , VAX/support.s
-       scalb()     ...in files IEEE/support.c , VAX/support.s
-       copysign()  ...in files IEEE/support.c , VAX/support.s
+        sqrt()      ...in files IEEE/support.c , VAX/sqrt.s
+        drem()      ...in files IEEE/support.c , VAX/support.s
+        finite()    ...in files IEEE/support.c , VAX/support.s
+        logb()      ...in files IEEE/support.c , VAX/support.s
+        scalb()     ...in files IEEE/support.c , VAX/support.s
+        copysign()  ...in files IEEE/support.c , VAX/support.s
+        rint()      ...in file  floor.c

 
 
    Notes: 
        i. The codes in files ending with .s are written in VAX assembly 
 
 
    Notes: 
        i. The codes in files ending with .s are written in VAX assembly 

-         language. They are intended for VAX computers.
+          language. They are intended for VAX computers.

 
 

-         Files that end with .c are written in C. They are intended
-         for either a VAX or a machine that conforms to the IEEE 
-         standard 754 for (double-precision) floating-point arithmetic.
+          Files that end with .c are written in C. They are intended
+          for either a VAX or a machine that conforms to the IEEE 
+          standard 754 for (double-precision) floating-point arithmetic.

 
       ii. On other than VAX or IEEE machines, run the original math 
 
       ii. On other than VAX or IEEE machines, run the original math 

-         library, formerly libm.a, now libom.a, if nothing better
-         is available.
+          library, formerly libm.a, now libom.a, if nothing better
+          is available.

 
      iii. The trigonometric functions sin/cos/tan/atan2 in files "VAX/sincos.s",
 
      iii. The trigonometric functions sin/cos/tan/atan2 in files "VAX/sincos.s",

-         "VAX/tan.s" and "VAX/atan2.s" are different from those in
-         "IEEE/trig.c" and "IEEE/atan2.c".  The VAX assembler code uses the
-         true value of pi to perform argument reduction, while the C code uses
-         a machine value of PI (see "IEEE/trig.c").
+          "VAX/tan.s" and "VAX/atan2.s" are different from those in
+          "IEEE/trig.c" and "IEEE/atan2.c".  The VAX assembler code uses the
+          true value of pi to perform argument reduction, while the C code uses
+          a machine value of PI (see "IEEE/trig.c").

 
 
 2. A computer system that conforms to IEEE standard 754 should provide 
 
 
 2. A computer system that conforms to IEEE standard 754 should provide 

-               sqrt(x),
-               drem(x,p), (double precision remainder function)
-               copysign(x,y),
-               finite(x),
-               scalb(x,N),
-               logb(x) .
+                sqrt(x),
+                drem(x,p), (double precision remainder function)
+                copysign(x,y),
+                finite(x),
+                scalb(x,N),
+                logb(x) and
+                rint(x).

    These functions are required or recommended by the standard.
    For convenience, a (slow) C implementation of these functions is 
    provided in the file IEEE/support.c.
    These functions are required or recommended by the standard.
    For convenience, a (slow) C implementation of these functions is 
    provided in the file IEEE/support.c.


@@ -104,21 +118,21 @@ Revised on 5/10/85, 5/13/85, 6/14/85, 8/20/85.
    National 32000, a Zilog 8000, a VAX, and a SUN (cf. the "Makefile" in
    this directory). Invoke the C compiler thus:
 
    National 32000, a Zilog 8000, a VAX, and a SUN (cf. the "Makefile" in
    this directory). Invoke the C compiler thus:
 

-       cc -c -DVAX IEEE/support.c              ... on a VAX, D-format
-       cc -c -DNATIONAL IEEE/support.c         ... on a National 32081
-       cc -c  IEEE/support.c                   ... on other IEEE machines,
-                                                   we hope.
+        cc -c -DVAX IEEE/support.c              ... on a VAX, D-format
+        cc -c -DNATIONAL IEEE/support.c         ... on a National 32081
+        cc -c  IEEE/support.c                   ... on other IEEE machines,
+                                                    we hope.

 
    Notes: 
       1. Faster versions of "drem" and "sqrt" for IEEE double precision
 
    Notes: 
       1. Faster versions of "drem" and "sqrt" for IEEE double precision

-        (coded in C but intended for assembly language) are given at the
-        end of support.c but commented out since they require certain
-        machine-dependent functions.
+         (coded in C but intended for assembly language) are given at the
+         end of support.c but commented out since they require certain
+         machine-dependent functions.

 
       2. A fast VAX assembler version of the system supported functions
 
       2. A fast VAX assembler version of the system supported functions

-        copysign(), logb(), scalb(), finite(), and drem() appears in file 
-        VAX/support.s.  A fast VAX assembler version of sqrt() is in
-        file sqrt.s .
+         copysign(), logb(), scalb(), finite(), and drem() appears in file 
+         VAX/support.s.  A fast VAX assembler version of sqrt() is in
+         file sqrt.s .

 
 3. Two formats are supported by all the standard elementary functions: 
    the VAX D-format (56 bits' precision), and the IEEE double format 
 
 3. Two formats are supported by all the standard elementary functions: 
    the VAX D-format (56 bits' precision), and the IEEE double format 


@@ -130,140 +144,121 @@ Revised on 5/10/85, 5/13/85, 6/14/85, 8/20/85.
    Makefile in this directory). 
 
     Example:
    Makefile in this directory). 
 
     Example:

-       cc -c -DVAX sin.c               ... for VAX D-format
+        cc -c -DVAX sin.c               ... for VAX D-format

 
        Warning: The values of floating-point constants used in the code are
 
        Warning: The values of floating-point constants used in the code are

-               given in both hexadecimal and decimal.  The hexadecimal values
-               are the intended ones. The decimal values may be use provided 
-               that the compiler converts from decimal to binary accurately
-               enough to produce the hexadecimal values shown. If the
-               conversion is inaccurate, then one must know the exact machine 
-               representation of the constants and alter the assembly-
-               language output from the compiler, or apply tricks like 
-               the following in a C program.
+                given in both hexadecimal and decimal.  The hexadecimal values
+                are the intended ones. The decimal values may be use provided 
+                that the compiler converts from decimal to binary accurately
+                enough to produce the hexadecimal values shown. If the
+                conversion is inaccurate, then one must know the exact machine 
+                representation of the constants and alter the assembly-
+                language output from the compiler, or apply tricks like 
+                the following in a C program.

 
 

-                       Example: to store the floating-point constant 
+                        Example: to store the floating-point constant 

 
 

-                            p1= 2^-6 * .F83ABE67E1066A (Hexadecimal)
+                             p1= 2^-6 * .F83ABE67E1066A (Hexadecimal)

 
 

-                       on a VAX in C, we use two long word to store its 
-                       machine value and define p1 to be the double constant 
-                       at the location of these two long words:
+                        on a VAX in C, we use two long word to store its 
+                        machine value and define p1 to be the double constant 
+                        at the location of these two long words:

 
 

-                       static long  p1x[] = { 0x3abe3d78, 0x066a67e1};
-                       #define      p1      (*(double*)p1x)
+                        static long  p1x[] = { 0x3abe3d78, 0x066a67e1};
+                        #define      p1      (*(double*)p1x)

 
     Note:  On a VAX, some functions have two codes. For example, cabs() 
 
     Note:  On a VAX, some functions have two codes. For example, cabs() 

-          has one implementation in cabs.c, and another in VAX/cabs.s. 
-          In this case, the assembly version is preferred.
+           has one implementation in cabs.c, and another in VAX/cabs.s. 
+           In this case, the assembly version is preferred.

 
 
 4. Accuracy. 
 
             The errors in expm1(), log1p(), exp(), log(), cabs(), hypot()
 
 
 4. Accuracy. 
 
             The errors in expm1(), log1p(), exp(), log(), cabs(), hypot()

-           and cbrt() are below 1 ULP (Unit in the Last Place).
+            and cbrt() are below 1 ULP (Unit in the Last Place).

 
 

-           The error in pow(x,y) grows with the size of y. Nevertheless,
-           for integers x and y, pow(x,y) returns the correct integer value 
-           on all tested machines (VAX, SUN, NATIONAL, ZILOG), provided that 
-           x to the power of y is representable exactly.
+            The error in pow(x,y) grows with the size of y. Nevertheless,
+            for integers x and y, pow(x,y) returns the correct integer value 
+            on all tested machines (VAX, SUN, NATIONAL, ZILOG), provided that 
+            x to the power of y is representable exactly.

 
 

-           cosh, sinh, acosh, asinh, tanh, atanh and log10 have errors below
-           about 3 ULPs. 
+            cosh, sinh, acosh, asinh, tanh, atanh and log10 have errors below
+            about 3 ULPs. 

 
 

-           For trigonometric and inverse trigonometric functions: 
+            For trigonometric and inverse trigonometric functions: 

 
 

-               Let [trig(x)] denote the value actually computed for trig(x),
+                Let [trig(x)] denote the value actually computed for trig(x),

 
 

-               1) Those codes using the machine's value PI (true pi rounded):
-                  (source codes: IEEE/{trig.c,atan2.c}, asincos.c and atan.c)
+                1) Those codes using the machine's value PI (true pi rounded):
+                   (source codes: IEEE/{trig.c,atan2.c}, asincos.c and atan.c)

 
 

-                  The errors in [sin(x)], [cos(x)], and [atan(x)] are below 
-                  1 ULP compared with sin(x*pi/PI), cos(x*pi/PI), and 
-                  atan(x)*PI/pi respectively, where PI is the machine's
-                  value of pi rounded. [Tan(x)] returns tan(x*pi/PI) within
-                  about 2 ULPs; [acos(x)], [asin(x)], and [atan2(y,x)] 
-                  return acos(x)*PI/pi, asin(x)*PI/pi, and atan2(y,x)*PI/pi
-                  respectively to similar accuracy.
+                   The errors in [sin(x)], [cos(x)], and [atan(x)] are below 
+                   1 ULP compared with sin(x*pi/PI), cos(x*pi/PI), and 
+                   atan(x)*PI/pi respectively, where PI is the machine's
+                   value of pi rounded. [Tan(x)] returns tan(x*pi/PI) within
+                   about 2 ULPs; [acos(x)], [asin(x)], and [atan2(y,x)] 
+                   return acos(x)*PI/pi, asin(x)*PI/pi, and atan2(y,x)*PI/pi
+                   respectively to similar accuracy.

 
 
 
 

-               2) Those using true pi (for VAX D-format only)
-                  (source codes: VAX/{sincos.s,tan.s,atan2.s}, asincos.c and
-                  atan.c)
+                2) Those using true pi (for VAX D-format only)
+                   (source codes: VAX/{sincos.s,tan.s,atan2.s}, asincos.c and
+                   atan.c)

 
 

-                  The errors in [sin(x)], [cos(x)], and [atan(x)] are below
-                  1 ULP. [Tan(x)], [atan2(y,x)], [acos(x)], and [asin(x)] 
-                  have errors below about 2 ULPs. 
+                   The errors in [sin(x)], [cos(x)], and [atan(x)] are below
+                   1 ULP. [Tan(x)], [atan2(y,x)], [acos(x)], and [asin(x)] 
+                   have errors below about 2 ULPs. 

 
 
 
 

-           Here are the results of some test runs to find worst errors on 
-           the VAX :
+            Here are the results of some test runs to find worst errors on 
+            the VAX :

 
 

-    tan   :  2.09 ULPs         ...1,024,000 random arguments (machine PI)
-    sin   :  .861 ULPs         ...1,024,000 random arguments (machine PI)
-    cos   :  .857 ULPs         ...1,024,000 random arguments (machine PI)
+    tan   :  2.09 ULPs          ...1,024,000 random arguments (machine PI)
+    sin   :  .861 ULPs          ...1,024,000 random arguments (machine PI)
+    cos   :  .857 ULPs          ...1,024,000 random arguments (machine PI)

     (compared with tan, sin, cos of (x*pi/PI))
 
     (compared with tan, sin, cos of (x*pi/PI))
 

-    acos  :  2.07 ULPs         .....200,000 random arguments (machine PI)
-    asin  :  2.06 ULPs         .....200,000 random arguments (machine PI)
-    atan2 :  1.41 ULPs         .....356,000 random arguments (machine PI)
-    atan  :  0.86 ULPs         ...1,536,000 random arguments (machine PI)
+    acos  :  2.07 ULPs          .....200,000 random arguments (machine PI)
+    asin  :  2.06 ULPs          .....200,000 random arguments (machine PI)
+    atan2 :  1.41 ULPs          .....356,000 random arguments (machine PI)
+    atan  :  0.86 ULPs          ...1,536,000 random arguments (machine PI)

     (compared with (PI/pi)*(atan, asin, acos, atan2 of x))
 
     (compared with (PI/pi)*(atan, asin, acos, atan2 of x))
 

-    tan   :  2.15 ULPs         ...1,024,000 random arguments (true pi)
-    sin   :  .814 ULPs         ...1,024,000 random arguments (true pi)
-    cos   :  .792 ULPs         ...1,024,000 random arguments (true pi)
-    acos  :  2.15 ULPs         ...1,024,000 random arguments (true pi)
-    asin  :  1.99 ULPs         ...1,024,000 random arguments (true pi)
-    atan2 :  1.48 ULPs         ...1,024,000 random arguments (true pi)
-    atan  :  .850 ULPs         ...1,024,000 random arguments (true pi)
-
-    acosh :  3.30 ULPs         .....512,000 random arguments
-    asinh :  1.58 ULPs         .....512,000 random arguments
-    atanh :  1.71 ULPs         .....512,000 random arguments  
-    cosh  :  1.23 ULPs         .....768,000 random arguments
-    sinh  :  1.93 ULPs         ...1,024,000 random arguments
-    tanh  :  2.22 ULPs         ...1,024,000 random arguments
-    log10 :  1.74 ULPs         ...1,536,000 random arguments
-    pow   :  1.79 ULPs         .....100,000 random arguments, 0 < x, y < 20.
-       
-    exp   :  .768 ULPs         ...1,156,000 random arguments
-    expm1 :  .844 ULPs         ...1,166,000 random arguments
-    log1p :  .846 ULPs         ...1,536,000 random arguments
-    log   :  .826 ULPs         ...1,536,000 random arguments
-    cabs  :  .959 ULPs         .....500,000 random arguments
+    tan   :  2.15 ULPs          ...1,024,000 random arguments (true pi)
+    sin   :  .814 ULPs          ...1,024,000 random arguments (true pi)
+    cos   :  .792 ULPs          ...1,024,000 random arguments (true pi)
+    acos  :  2.15 ULPs          ...1,024,000 random arguments (true pi)
+    asin  :  1.99 ULPs          ...1,024,000 random arguments (true pi)
+    atan2 :  1.48 ULPs          ...1,024,000 random arguments (true pi)
+    atan  :  .850 ULPs          ...1,024,000 random arguments (true pi)
+
+    acosh :  3.30 ULPs          .....512,000 random arguments
+    asinh :  1.58 ULPs          .....512,000 random arguments
+    atanh :  1.71 ULPs          .....512,000 random arguments  
+    cosh  :  1.23 ULPs          .....768,000 random arguments
+    sinh  :  1.93 ULPs          ...1,024,000 random arguments
+    tanh  :  2.22 ULPs          ...1,024,000 random arguments
+    log10 :  1.74 ULPs          ...1,536,000 random arguments
+    pow   :  1.79 ULPs          .....100,000 random arguments, 0 < x, y < 20.
+        
+    exp   :  .768 ULPs          ...1,156,000 random arguments
+    expm1 :  .844 ULPs          ...1,166,000 random arguments
+    log1p :  .846 ULPs          ...1,536,000 random arguments
+    log   :  .826 ULPs          ...1,536,000 random arguments
+    cabs  :  .959 ULPs          .....500,000 random arguments

     cbrt  :  .666 ULPs          ...5,120,000 random arguments
 
 
 5. Speed.
 
     cbrt  :  .666 ULPs          ...5,120,000 random arguments
 
 
 5. Speed.
 

-       Some functions coded in VAX assembly language (cabs, hypot, sqrt)
-       are significantly faster than the corresponding ones in 4.2BSD.
-       In general, to improve performance, all functions in IEEE/support.c
-       should be written in assembler and, whenever possible, should be
-       called via short subroutine calls.
+        Some functions coded in VAX assembly language (cabs, hypot, sqrt)
+        are significantly faster than the corresponding ones in 4.2BSD.
+        In general, to improve performance, all functions in IEEE/support.c
+        should be written in assembler and, whenever possible, should be
+        called via short subroutine calls.

 
 
 6. j0,j1,jn.
 
 
 
 6. j0,j1,jn.
 

-       The modications to these routines were only in how an invalide
-       floating point operations is signaled.
-
-
-7. Copyright notice, and Disclaimer:
-
-***************************************************************************
-*                                                                        * 
-* Copyright (c) 1985 Regents of the University of California.            *
-*                                                                        * 
-* Use and reproduction of this software are granted  in  accordance  with *
-* the terms and conditions specified in  the  Berkeley  Software  License *
-* Agreement (in particular, this entails acknowledgement of the programs' *
-* source, and inclusion of this notice) with the additional understanding *
-* that  all  recipients  should regard themselves as participants  in  an *
-* ongoing  research  project and hence should  feel  obligated  to report *
-* their  experiences (good or bad) with these elementary function  codes, *
-* using "sendbug 4bsd-bugs@BERKELEY", to the authors.                    *
-*                                                                        *
-***************************************************************************
-
+        The modications to these routines were only in how an invalide
+        floating point operations is signaled.

diff --git a/usr/src/old/libm/libm/floor.c b/usr/src/old/libm/libm/floor.c
index 7fce6b8..56916da 100644 (file)
--- a/usr/src/old/libm/libm/floor.c
+++ b/usr/src/old/libm/libm/floor.c
@@ -1,4 +1,4 @@
-/*     @(#)floor.c     4.1     %G%     */
+/*     @(#)floor.c     4.2     %G% */

 
 /*
  * floor and ceil-- greatest integer <= arg
 
 /*
  * floor and ceil-- greatest integer <= arg


@@ -30,3 +30,46 @@ double d;
 {
        return(-floor(-d));
 }
 {
        return(-floor(-d));
 }

+
+/*
+ * algorithm for rint(x) in pseudo-pascal form ...
+ *
+ * real rint(x): real x;
+ *     ... delivers integer nearest x in direction of prevailing rounding
+ *     ... mode
+ * const       L = (last consecutive integer)/2
+ *       = 2**55; for VAX D
+ *       = 2**52; for IEEE 754 Double
+ * real        s,t;
+ * begin
+ *     if x != x then return x;                ... NaN
+ *     if |x| >= L then return x;              ... already an integer
+ *     s := copysign(L,x);
+ *     t := x + s;                             ... = (x+s) rounded to integer
+ *     return t - s
+ * end;
+ *
+ * Note: Inexact will be signaled if x is not an integer, as is
+ *     customary for IEEE 754.  No other signal can be emitted.
+ */
+#ifdef VAX
+static long Lx[] = {0x5c00,0x0};               /* 2**55 */
+#define L *(double *) Lx
+#else  /* IEEE double */
+static double L = 4503599627370496.0E0;                /* 2**52 */
+#endif
+double
+rint(x)
+double x;
+{
+       double s,t,one = 1.0,copysign();
+#ifndef VAX
+       if (x != x)                             /* NaN */
+               return (x);
+#endif
+       if (copysign(x,one) >= L)               /* already an integer */
+           return (x);
+       s = copysign(L,x);
+       t = x + s;                              /* x+s rounded to integer */
+       return (t - s);
+}

diff --git a/usr/src/old/libm/libm/lgamma.c b/usr/src/old/libm/libm/lgamma.c
index 309d056..5cdd53e 100644 (file)
--- a/usr/src/old/libm/libm/lgamma.c
+++ b/usr/src/old/libm/libm/lgamma.c
@@ -1,11 +1,13 @@
-/*     @(#)lgamma.c    4.4     %G% */
+#ifndef lint
+static char sccsid[] = "@(#)lgamma.c   4.4 (Berkeley) %G%";
+#endif not lint

 
 /*
 
 /*

-       C program for floating point log gamma function
+       C program for floating point log Gamma function

 
 

-       gamma(x) computes the log of the absolute
-       value of the gamma function.
-       The sign of the gamma function is returned in the
+       lgamma(x) computes the log of the absolute
+       value of the Gamma function.
+       The sign of the Gamma function is returned in the

        external quantity signgam.
 
        The coefficients for expansion around zero
        external quantity signgam.
 
        The coefficients for expansion around zero


@@ -15,15 +17,10 @@
        Calls log, floor and sin.
 */
 
        Calls log, floor and sin.
 */
 

-#include <errno.h>

 #include <math.h>
 #ifdef VAX
 #include <math.h>
 #ifdef VAX

-static long    NaN_[] = {0x8000, 0x0};
-#define NaN    (*(double *) NaN_)
+#include <errno.h>

 #endif
 #endif

-
-
-int    errno;

 int    signgam = 0;
 static double goobie   = 0.9189385332046727417803297;  /* log(2*pi)/2 */
 static double pi       = 3.1415926535897932384626434;
 int    signgam = 0;
 static double goobie   = 0.9189385332046727417803297;  /* log(2*pi)/2 */
 static double pi       = 3.1415926535897932384626434;


@@ -60,7 +57,7 @@ static double q2[] = {
 };
 
 double
 };
 
 double

-gamma(arg)
+lgamma(arg)

 double arg;
 {
        double log(), pos(), neg(), asym();
 double arg;
 {
        double log(), pos(), neg(), asym();


@@ -106,18 +103,13 @@ double arg;
       */
        t = floor(arg);
        if (arg - t  > 0.5e0)
       */
        t = floor(arg);
        if (arg - t  > 0.5e0)

-           t += 1.e0;                  /* t := integer nearest arg */
-#if !(IEEE|NATIONAL)
-       if(arg == t) {
-               errno = EDOM;
+           t += 1.e0;                          /* t := integer nearest arg */

 #ifdef VAX
 #ifdef VAX

-               return (NaN);
-#else
-               return(HUGE);
-#endif
+       if (arg == t) {
+           extern double infnan();
+           return(infnan(ERANGE));             /* +INF */

        }
 #endif
        }
 #endif

-

        signgam = (int) (t - 2*floor(t/2));     /* signgam =  1 if t was odd, */
                                                /*            0 if t was even */
        signgam = signgam - 1 + signgam;        /* signgam =  1 if t was odd, */
        signgam = (int) (t - 2*floor(t/2));     /* signgam =  1 if t was odd, */
                                                /*            0 if t was even */
        signgam = signgam - 1 + signgam;        /* signgam =  1 if t was odd, */