double tanh(double);
double exp(double);
+double exp2(double);
double frexp(double, int *);
double ldexp(double, int);
double log(double);
/* 7.12.6 exp / log */
float expf(float);
+float exp2f(float);
float expm1f(float);
float frexpf(float, int *);
long double frexpl(long double, int *);
float cbrtf(float);
float fabsf(float);
+long double fabsl(long double);
float hypotf(float, float);
float powf(float, float);
float sqrtf(float);
float fmodf(float, float);
float remainderf(float, float);
+/* 7.12.10.3 The remquo functions */
+double remquo(double, double, int *);
+float remquof(float, float, int *);
+
/* 7.12.11 manipulation */
float copysignf(float, float);
-long double copysignl(long double x, long double y);
+long double copysignl(long double, long double);
double nan(const char *);
float nanf(const char *);
long double nanl(const char *);
float nextafterf(float, float);
+long double nextafterl(long double, long double);
+double nexttoward (double, long double);
/* 7.12.12 maximum, minimum, positive difference */
double fdim(double, double);
int __isnanl(long double);
int __signbitl(long double);
#endif
+
+int ilogbl(long double);
+long double logbl(long double);
+long double scalbnl(long double, int);
__END_DECLS
#endif /* _MATH_H_ */
s_log1p.S s_log1pf.S s_logb.S s_logbf.S s_rint.S s_rintf.S \
s_scalbn.S s_scalbnf.S s_significand.S s_significandf.S s_sin.S \
s_sinf.S s_tan.S s_tanf.S
+SRCS+= s_ilogbl.S s_logbl.S s_modf.S s_scalbnl.S
# This file is included by arch/x86_64/Makefile.inc, so pull in the right
# fenv.[ch].
+/* $NetBSD: abi.h,v 1.7 2011/06/18 22:19:52 joerg Exp $ */
+
/*
* Written by Frank van der Linden (fvdl@wasabisystems.com)
- *
- * $NetBSD: abi.h,v 1.3 2004/03/22 13:41:09 wiz Exp $
- * $DragonFly: src/lib/libm/arch/i386/abi.h,v 1.1 2005/07/26 21:15:19 joerg Exp $
*/
/*
- * The x86-64 ABI specifies that float, double and long double
- * arguments are passed in SSE2 (xmm) registers. Unfortunately,
+ * The x86-64 ABI specifies that float and double arguments
+ * are passed in SSE2 (xmm) registers. Unfortunately,
* there is no way to push those on to the FP stack, which is
* where the fancier instructions get their arguments from.
*
#ifdef __x86_64__
+#define ARG_LONG_DOUBLE_ONE 8(%rsp)
+#define ARG_LONG_DOUBLE_TWO 24(%rsp)
#define ARG_DOUBLE_ONE -8(%rsp)
+#define ARG_DOUBLE_ONE_LSW -8(%rsp)
+#define ARG_DOUBLE_ONE_MSW -4(%rsp)
#define ARG_DOUBLE_TWO -16(%rsp)
#define ARG_FLOAT_ONE -4(%rsp)
#define ARG_FLOAT_TWO -8(%rsp)
#else
+#define ARG_LONG_DOUBLE_ONE 4(%esp)
+#define ARG_LONG_DOUBLE_TWO 16(%esp)
#define ARG_DOUBLE_ONE 4(%esp)
+#define ARG_DOUBLE_ONE_LSW 4(%esp)
+#define ARG_DOUBLE_ONE_MSW 8(%esp)
#define ARG_DOUBLE_TWO 12(%esp)
#define ARG_FLOAT_ONE 4(%esp)
#define ARG_FLOAT_TWO 8(%esp)
+/* $NetBSD: e_exp.S,v 1.14 2008/06/23 10:24:13 drochner Exp $ */
+
/*
* Copyright (c) 1993,94 Winning Strategies, Inc.
* All rights reserved.
* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
- *
- * $NetBSD: e_exp.S,v 1.12 2002/02/27 16:32:46 christos Exp $
- * $FreeBSD: src/lib/msun/i387/e_exp.S,v 1.8.2.1 2000/07/10 09:16:28 obrien Exp $
*/
/*
/* e^x = 2^(x * log2(e)) */
ENTRY(exp)
-#ifndef __i386__
- /*
- * XXX: This code is broken and needs to be merged with the i386 case.
- */
- fstcw -12(%rsp)
- movw -12(%rsp),%dx
- orw $0x0180,%dx
- movw %dx,-16(%rsp)
- fldcw -16(%rsp)
- movsd %xmm0,-8(%rsp)
- fldl -8(%rsp)
-
- fldl2e
- fmulp /* x * log2(e) */
- fld %st(0)
- frndint /* int(x * log2(e)) */
- fxch %st(1)
- fsub %st(1),%st /* fract(x * log2(e)) */
- f2xm1 /* 2^(fract(x * log2(e))) - 1 */
- fld1
- faddp /* 2^(fract(x * log2(e))) */
- fscale /* e^x */
- fstp %st(1)
-
- fstpl -8(%rsp)
- movsd -8(%rsp),%xmm0
- fldcw -12(%rsp)
- ret
-#else
+ XMM_ONE_ARG_DOUBLE_PROLOGUE
/*
* If x is +-Inf, then the subtraction would give Inf-Inf = NaN.
* Avoid this. Also avoid it if x is NaN for convenience.
*/
- movl 8(%esp),%eax
- andl $0x7fffffff,%eax
- cmpl $0x7ff00000,%eax
+ movl ARG_DOUBLE_ONE_MSW, %eax
+ andl $0x7fffffff, %eax
+ cmpl $0x7ff00000, %eax
jae x_Inf_or_NaN
- fldl 4(%esp)
+ fldl ARG_DOUBLE_ONE
/*
* Ensure that the rounding mode is to nearest (to give the smallest
* possible fraction) and that the precision is as high as possible.
* We may as well mask interrupts if we switch the mode.
*/
- fstcw 4(%esp)
- movl 4(%esp),%eax
- andl $0x0300,%eax
- cmpl $0x0300,%eax /* RC == 0 && PC == 3? */
+#define CWSTORE_SAV ARG_DOUBLE_ONE_LSW /* XXX overwrites the argument */
+#define CWSTORE_TMP ARG_DOUBLE_ONE_MSW
+ fstcw CWSTORE_SAV
+ movl CWSTORE_SAV, %eax
+ andl $0x0f00, %eax
+ cmpl $0x0300, %eax /* RC == 0 && PC == 3? */
je 1f /* jump if mode is good */
- movl $0x137f,8(%esp)
- fldcw 8(%esp)
+ movl $0x137f, CWSTORE_TMP
+ fldcw CWSTORE_TMP
1:
fldl2e
fmulp /* x * log2(e) */
fscale /* e^x */
fstp %st(1)
je 1f
- fldcw 4(%esp)
+ fldcw CWSTORE_SAV
1:
+ XMM_DOUBLE_EPILOGUE
ret
-
x_Inf_or_NaN:
/*
* Return 0 if x is -Inf. Otherwise just return x, although the
* C version would return (x + x) (Real Indefinite) if x is a NaN.
*/
- cmpl $0xfff00000,8(%esp)
+ cmpl $0xfff00000, ARG_DOUBLE_ONE_MSW
jne x_not_minus_Inf
- cmpl $0,4(%esp)
+ cmpl $0, ARG_DOUBLE_ONE_LSW
jne x_not_minus_Inf
fldz
+ XMM_DOUBLE_EPILOGUE
ret
x_not_minus_Inf:
- fldl 4(%esp)
+ fldl ARG_DOUBLE_ONE
+ XMM_DOUBLE_EPILOGUE
ret
-#endif
END(exp)
* Written by J.T. Conklin <jtc@NetBSD.org>.
* Public domain.
*
- * $NetBSD: e_expf.S,v 1.5 2003/07/26 19:24:58 salo Exp $
+ * $NetBSD: e_expf.S,v 1.6 2008/06/24 17:27:56 drochner Exp $ $
*/
#include <machine/asm.h>
/* e^x = 2^(x * log2(e)) */
ENTRY(expf)
XMM_ONE_ARG_FLOAT_PROLOGUE
+
+ /*
+ * catch +/-Inf and NaN arguments
+ */
+ movl ARG_FLOAT_ONE,%eax
+ andl $0x7fffffff,%eax
+ cmpl $0x7f800000,%eax
+ jae x_Inf_or_NaN
+
flds ARG_FLOAT_ONE
fldl2e
fmulp /* x * log2(e) */
fstp %st(1)
XMM_FLOAT_EPILOGUE
ret
+
+x_Inf_or_NaN:
+ /*
+ * Return 0 if x is -Inf. Otherwise just return x, although the
+ * C version would return (x + x) (Real Indefinite) if x is a NaN.
+ */
+ movl ARG_FLOAT_ONE,%eax
+ cmpl $0xff800000,%eax
+ jne x_not_minus_Inf
+ fldz
+ XMM_FLOAT_EPILOGUE
+ ret
+
+x_not_minus_Inf:
+ flds ARG_FLOAT_ONE
+ XMM_FLOAT_EPILOGUE
+ ret
END(expf)
* Written by J.T. Conklin <jtc@NetBSD.org>.
* Public domain.
*
- * $NetBSD: s_ceil.S,v 1.7 2003/07/26 19:25:00 salo Exp $
+ * $NetBSD: s_ceil.S,v 1.8 2011/06/18 21:24:51 joerg Exp
*/
#include <machine/asm.h>
movw %dx,-8(%ebp)
fldcw -8(%ebp) /* load modfied control word */
- fldl 8(%ebp); /* round */
+ fldl 8(%ebp) /* round */
frndint
fldcw -4(%ebp) /* restore original control word */
* Written by J.T. Conklin <jtc@NetBSD.org>.
* Public domain.
*
- * $NetBSD: s_ceilf.S,v 1.8 2004/07/16 18:40:24 drochner Exp $
+ * $NetBSD: s_ceilf.S,v 1.9 2011/06/18 21:24:51 joerg Exp $
*/
#include <machine/asm.h>
movw %dx,-8(%ebp)
fldcw -8(%ebp) /* load modfied control word */
- flds 8(%ebp); /* round */
+ flds 8(%ebp) /* round */
frndint
fldcw -4(%ebp) /* restore original control word */
* Written by J.T. Conklin <jtc@NetBSD.org>.
* Public domain.
*
- * $NetBSD: s_copysign.S,v 1.6 2003/07/26 19:25:01 salo Exp $
+ * $NetBSD: s_copysign.S,v 1.7 2011/06/18 20:49:26 joerg Exp
*/
/*
movl %eax,8(%esp)
fldl 4(%esp)
#else
-#if 0
- /*
- * XXXfvdl gas doesn't grok this yet.
- */
movq .Lpos(%rip),%xmm2
movq .Lneg(%rip),%xmm3
pand %xmm2,%xmm1
pand %xmm3,%xmm0
por %xmm1,%xmm0
-#else
- movsd %xmm0,-8(%rsp)
- movsd %xmm1,-16(%rsp)
- movl -12(%rsp),%edx
- andl $0x80000000,%edx
- movl -4(%rsp),%eax
- andl $0x7fffffff,%eax
- orl %edx,%eax
- movl %eax,-4(%rsp)
- movsd -8(%rsp),%xmm0
-#endif
#endif
ret
END(copysign)
* Written by J.T. Conklin <jtc@NetBSD.org>.
* Public domain.
*
- * $NetBSD: s_copysignf.S,v 1.5 2003/07/26 19:25:01 salo Exp $
+ * $NetBSD: s_copysignf.S,v 1.7 2011/06/21 21:52:49 joerg Exp
*/
#include <machine/asm.h>
movl %eax,4(%esp)
flds 4(%esp)
#else
-#if 0
- /*
- * XXXfvdl gas doesn't grok this.
- * but it's legal according to the p4 manual.
- */
- movss .Lpos(%rip),%xmm2
- movss .Lneg(%rip),%xmm3
- pandq %xmm2,%xmm1
- pandq %xmm3,%xmm0
- porq %xmm1,%xmm0
-#else
- movss %xmm0,-4(%rsp)
- movss %xmm1,-8(%rsp)
- movl -8(%rsp),%edx
- andl $0x80000000,%edx
- movl -4(%rsp),%eax
- andl $0x7fffffff,%eax
- orl %edx,%eax
- movl %eax,-4(%rsp)
- movss -4(%rsp),%xmm0
-#endif
+ movss .Lpos(%rip),%xmm2
+ movss .Lneg(%rip),%xmm3
+ pand %xmm2,%xmm1
+ pand %xmm3,%xmm0
+ por %xmm1,%xmm0
#endif
ret
END(copysignf)
* Written by J.T. Conklin <jtc@NetBSD.org>.
* Public domain.
*
- * $NetBSD: s_floor.S,v 1.8 2003/07/26 19:25:02 salo Exp $
+ * $NetBSD: s_floor.S,v 1.9 2011/06/18 21:24:51 joerg Exp $
*/
#include <machine/asm.h>
movw %dx,-8(%ebp)
fldcw -8(%ebp) /* load modfied control word */
- fldl 8(%ebp); /* round */
+ fldl 8(%ebp) /* round */
frndint
fldcw -4(%ebp) /* restore original control word */
* Written by J.T. Conklin <jtc@NetBSD.org>.
* Public domain.
*
- * $NetBSD: s_floorf.S,v 1.7 2004/07/16 18:40:24 drochner Exp $
+ * $NetBSD: s_floorf.S,v 1.8 2011/06/18 21:24:51 joerg Exp $
*/
#include <machine/asm.h>
movw %dx,-8(%ebp)
fldcw -8(%ebp) /* load modfied control word */
- flds 8(%ebp); /* round */
+ flds 8(%ebp) /* round */
frndint
fldcw -4(%ebp) /* restore original control word */
--- /dev/null
+/*
+ * Written by J.T. Conklin <jtc@NetBSD.org>.
+ * Public domain.
+ *
+ * $NetBSD: s_ilogbl.S,v 1.1 2011/07/28 22:32:28 joerg Exp $
+ */
+
+#include <machine/asm.h>
+
+#include "abi.h"
+
+ENTRY(ilogbl)
+ fldt ARG_LONG_DOUBLE_ONE
+ fxtract
+ fstp %st
+#ifdef __i386__
+ pushl %eax
+ fistpl 0(%esp)
+ popl %eax
+#else
+ fistpl -4(%rsp)
+ movl -4(%rsp), %eax
+#endif
+ ret
+END(ilogbl)
\ No newline at end of file
--- /dev/null
+/*
+ * Written by J.T. Conklin <jtc@NetBSD.org>.
+ * Public domain.
+ *
+ * $NetBSD: s_logbl.S,v 1.1 2011/08/03 14:13:07 joerg Exp
+ */
+
+#include <machine/asm.h>
+
+#include "abi.h"
+
+ENTRY(logbl)
+ fldt ARG_LONG_DOUBLE_ONE
+ fxtract
+ fstp %st
+ ret
+END(logbl)
\ No newline at end of file
--- /dev/null
+/* $NetBSD: s_modf.S,v 1.1 2006/03/22 20:45:58 drochner Exp $ */
+
+/*-
+ * Copyright (c) 1990 The Regents of the University of California.
+ * All rights reserved.
+ *
+ * This code is derived from software contributed to Berkeley by
+ * Sean Eric Fagan.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions and the following disclaimer in the
+ * documentation and/or other materials provided with the distribution.
+ * 3. Neither the name of the University nor the names of its contributors
+ * may be used to endorse or promote products derived from this software
+ * without specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ *
+ * from: @(#)modf.s 5.5 (Berkeley) 3/18/91
+ */
+
+#include <machine/asm.h>
+
+/*
+ * modf(value, iptr): return fractional part of value, and stores the
+ * integral part into iptr (a pointer to double).
+ *
+ * Written by Sean Eric Fagan (sef@kithrup.COM)
+ * Sun Mar 11 20:27:30 PST 1990
+ */
+
+/* With CHOP mode on, frndint behaves as TRUNC does. Useful. */
+ENTRY(modf)
+#ifdef __x86_64__
+ pushq %rbp
+ movq %rsp,%rbp
+ subq $24,%rsp
+
+ /* Set chop mode. */
+ fnstcw -12(%rbp)
+ movw -12(%rbp),%dx
+ orw $3072,%dx
+ movw %dx,-16(%rbp)
+ fldcw -16(%rbp)
+
+ /* Get integral part. */
+ movsd %xmm0,-24(%rbp)
+ fldl -24(%rbp)
+ frndint
+ fstpl -8(%rbp)
+
+ /* Restore control word. */
+ fldcw -12(%rbp)
+
+ /* Store integral part. */
+ movsd -8(%rbp),%xmm0
+ movsd %xmm0,(%rdi)
+
+ /* Get fractional part and return it. */
+ fldl -24(%rbp)
+ fsubl -8(%rbp)
+ fstpl -24(%rbp)
+ movsd -24(%rbp),%xmm0
+#else
+ pushl %ebp
+ movl %esp,%ebp
+ subl $16,%esp
+ fnstcw -12(%ebp)
+ movw -12(%ebp),%dx
+ orw $3072,%dx
+ movw %dx,-16(%ebp)
+ fldcw -16(%ebp)
+ fldl 8(%ebp)
+ frndint
+ fstpl -8(%ebp)
+ fldcw -12(%ebp)
+ movl 16(%ebp),%eax
+ movl -8(%ebp),%edx
+ movl -4(%ebp),%ecx
+ movl %edx,(%eax)
+ movl %ecx,4(%eax)
+ fldl 8(%ebp)
+ fsubl -8(%ebp)
+ jmp L1
+L1:
+#endif
+ leave
+ ret
+END(modf)
* Written by J.T. Conklin <jtc@NetBSD.org>.
* Public domain.
*
- * $NetBSD: s_scalbn.S,v 1.8 2006/03/21 11:35:21 drochner Exp $
+ * $NetBSD: s_scalbn.S,v 1.9 2010/04/23 19:17:07 drochner Exp $
*/
#include <machine/asm.h>
* Written by J.T. Conklin <jtc@NetBSD.org>.
* Public domain.
*
- * $NetBSD: s_scalbnf.S,v 1.7 2006/03/21 11:35:21 drochner Exp $
+ * $NetBSD: s_scalbnf.S,v 1.8 2010/04/23 19:17:07 drochner Exp $
*/
#include <machine/asm.h>
--- /dev/null
+/*
+ * Written by J.T. Conklin <jtc@NetBSD.org>.
+ * Public domain.
+ *
+ * $NetBSD: s_scalbnl.S,v 1.1 2011/07/26 17:03:23 joerg Exp $
+ */
+
+#include <machine/asm.h>
+
+ENTRY(scalbnl)
+#ifdef __x86_64__
+ movl %edi,-4(%rsp)
+ fildl -4(%rsp)
+ fldt 8(%rsp)
+ fscale
+ fstp %st(1)
+#else
+ fildl 16(%esp)
+ fldt 4(%esp)
+ fscale
+ fstp %st(1) /* clean up stack */
+#endif
+ ret
+END(scalbnl)
#include <complex.h>
#include <math.h>
-#include "../src/math_private.h"
double
cabs(double complex z)
{
- return hypot(creal(z), cimag(z));
+ return hypot(__real__ z, __imag__ z);
}
#include <complex.h>
#include <math.h>
-#include "../src/math_private.h"
float
cabsf(float complex z)
{
- return hypotf(crealf(z), cimagf(z));
+ return hypotf(__real__ z, __imag__ z);
}
-/* $NetBSD: cacosh.c,v 1.1 2007/08/20 16:01:30 drochner Exp $ */
+/* $NetBSD: cacosh.c,v 1.2 2009/08/03 19:41:32 drochner Exp $ */
/*-
* Copyright (c) 2007 The NetBSD Foundation, Inc.
-/* $NetBSD: cacoshf.c,v 1.1 2007/08/20 16:01:31 drochner Exp $ */
+/* $NetBSD: cacoshf.c,v 1.2 2009/08/03 19:41:32 drochner Exp $ */
/*-
* Copyright (c) 2007 The NetBSD Foundation, Inc.
carg(double complex z)
{
- return atan2(cimag(z), creal(z));
+ return atan2(__imag__ z, __real__ z);
}
cargf(float complex z)
{
- return atan2f(cimag(z), creal(z));
+ return atan2f(__imag__ z, __real__ z);
}
-/* $NetBSD: catan.c,v 1.1 2007/08/20 16:01:32 drochner Exp $ */
+/* $NetBSD: catan.c,v 1.2 2011/07/03 06:45:24 mrg Exp $ */
/*-
* Copyright (c) 2007 The NetBSD Foundation, Inc.
#include <complex.h>
#include <math.h>
+#include <float.h>
#include "cephes_subr.h"
#ifdef __weak_alias
__weak_alias(catan, _catan)
#endif
-#define MAXNUM 1.0e308
+#define MAXNUM DBL_MAX
double complex
catan(double complex z)
-/* $NetBSD: catanf.c,v 1.1 2007/08/20 16:01:32 drochner Exp $ */
+/* $NetBSD: catanf.c,v 1.2 2011/07/03 06:45:24 mrg Exp $ */
/*-
* Copyright (c) 2007 The NetBSD Foundation, Inc.
#include <complex.h>
#include <math.h>
+#include <float.h>
#include "cephes_subrf.h"
#ifdef __weak_alias
__weak_alias(catanf, _catanf)
#endif
-#define MAXNUMF 1.0e38F
+#define MAXNUMF FLT_MAX
float complex
catanf(float complex z)
-.\" $NetBSD: cimag.3,v 1.1 2008/02/20 09:55:38 drochner Exp $
+.\" $NetBSD: cimag.3,v 1.3 2010/09/15 18:40:27 wiz Exp $
.\" Copyright (c) 2001-2003 The Open Group, All Rights Reserved
.TH "CIMAG" 3P 2003 "IEEE/The Open Group" "POSIX Programmer's Manual"
.\" cimag
.br
float cimagf(float complex\fP \fIz\fP\fB);
.br
+long double cimagl(long double complex\fP \fIz\fP\fB);
+.br
\fP
.SH DESCRIPTION
.LP
-/*-
- * Copyright (c) 2004 Stefan Farfeleder
- * All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- * notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- * notice, this list of conditions and the following disclaimer in the
- * documentation and/or other materials provided with the distribution.
- *
- * THIS SOFTWARE IS PROVIDED BY AUTHOR AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED. IN NO EVENT SHALL AUTHOR OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- *
- * $FreeBSD$
+/* $NetBSD: cimag.c,v 1.2 2010/09/15 16:11:29 christos Exp $ */
+
+/*
+ * Written by Matthias Drochner <drochner@NetBSD.org>.
+ * Public domain.
*/
#include <complex.h>
double
cimag(double complex z)
{
- const double_complex z1 = { .f = z };
+ double_complex w = { .z = z };
- return (IMAGPART(z1));
+ return (IMAG_PART(w));
}
-/*-
- * Copyright (c) 2004 Stefan Farfeleder
- * All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- * notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- * notice, this list of conditions and the following disclaimer in the
- * documentation and/or other materials provided with the distribution.
- *
- * THIS SOFTWARE IS PROVIDED BY AUTHOR AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED. IN NO EVENT SHALL AUTHOR OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- *
- * $FreeBSD$
+/* $NetBSD: cimagf.c,v 1.2 2010/09/15 16:11:29 christos Exp $ */
+
+/*
+ * Written by Matthias Drochner <drochner@NetBSD.org>.
+ * Public domain.
*/
#include <complex.h>
float
cimagf(float complex z)
{
- const float_complex z1 = { .f = z };
+ float_complex w = { .z = z };
- return (IMAGPART(z1));
+ return (IMAG_PART(w));
}
+/* $NetBSD: cimagl.c,v 1.3 2010/09/20 16:55:20 christos Exp $ */
+
/*-
- * Copyright (c) 2004 Stefan Farfeleder
+ * Copyright (c) 2010 The NetBSD Foundation, Inc.
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
- * THIS SOFTWARE IS PROVIDED BY AUTHOR AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED. IN NO EVENT SHALL AUTHOR OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- *
- * $FreeBSD$
+ * THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
+ * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+ * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+ * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR CONTRIBUTORS
+ * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
*/
#include <complex.h>
#include "../src/math_private.h"
+/*
+ * cimagl(long double complex z)
+ * This function returns the imaginary part value (as a real) of z.
+ */
long double
cimagl(long double complex z)
{
- const long_double_complex z1 = { .f = z };
+ long_double_complex w = { .z = z };
- return (IMAGPART(z1));
+ return (IMAG_PART(w));
}
-.\" $NetBSD: conj.3,v 1.1 2008/02/20 09:55:38 drochner Exp $
+.\" $NetBSD: conj.3,v 1.3 2010/09/15 18:40:27 wiz Exp $
.\" Copyright (c) 2001-2003 The Open Group, All Rights Reserved
.TH "CONJ" 3P 2003 "IEEE/The Open Group" "POSIX Programmer's Manual"
.\" conj
.SH NAME
-conj, conjf, conjl \- complex conjugate functions
+conj, conjf \- complex conjugate functions
.SH SYNOPSIS
.LP
\fB#include <complex.h>
.br
float complex conjf(float complex\fP \fIz\fP\fB);
.br
-long double complex conjl(long double complex\fP \fIz\fP\fB);
+long double complex conjl(long double complex\fP \fIz\fP\fB, long double complex\fP \fIz\fP\fB);
.br
\fP
.SH DESCRIPTION
-/*-
- * Copyright (c) 2004 Stefan Farfeleder
- * All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- * notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- * notice, this list of conditions and the following disclaimer in the
- * documentation and/or other materials provided with the distribution.
- *
- * THIS SOFTWARE IS PROVIDED BY AUTHOR AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED. IN NO EVENT SHALL AUTHOR OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- *
- * $FreeBSD$
+/* $NetBSD: conj.c,v 1.2 2010/09/15 16:11:29 christos Exp $ */
+
+/*
+ * Written by Matthias Drochner <drochner@NetBSD.org>.
+ * Public domain.
*/
#include <complex.h>
-
#include "../src/math_private.h"
double complex
conj(double complex z)
{
+ double_complex w = { .z = z };
+
+ IMAG_PART(w) = -IMAG_PART(w);
- return (cpack(creal(z), -cimag(z)));
+ return (w.z);
}
-/*-
- * Copyright (c) 2004 Stefan Farfeleder
- * All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- * notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- * notice, this list of conditions and the following disclaimer in the
- * documentation and/or other materials provided with the distribution.
- *
- * THIS SOFTWARE IS PROVIDED BY AUTHOR AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED. IN NO EVENT SHALL AUTHOR OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- *
- * $FreeBSD$
+/* $NetBSD: conjf.c,v 1.2 2010/09/15 16:11:29 christos Exp $ */
+
+/*
+ * Written by Matthias Drochner <drochner@NetBSD.org>.
+ * Public domain.
*/
#include <complex.h>
-
#include "../src/math_private.h"
float complex
conjf(float complex z)
{
+ float_complex w = { .z = z };
+
+ IMAG_PART(w) = -IMAG_PART(w);
- return (cpackf(crealf(z), -cimagf(z)));
+ return (w.z);
}
+/* $NetBSD: conjl.c,v 1.4 2010/09/20 16:55:20 christos Exp $ */
+
/*-
- * Copyright (c) 2004 Stefan Farfeleder
+ * Copyright (c) 2010 The NetBSD Foundation, Inc.
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
- * THIS SOFTWARE IS PROVIDED BY AUTHOR AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED. IN NO EVENT SHALL AUTHOR OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- *
- * $FreeBSD$
+ * THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
+ * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+ * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+ * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR CONTRIBUTORS
+ * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
*/
#include <complex.h>
+#include <math.h>
#include "../src/math_private.h"
+/*
+ * conjl(long double complex z)
+ * This function returns the complex conjugate value of its argument, z.
+ */
long double complex
conjl(long double complex z)
{
+ long_double_complex w = { .z = z };
+
+ IMAG_PART(w) = -IMAG_PART(w);
- return (cpackl(creall(z), -cimagl(z)));
+ return (w.z);
}
+.\" $NetBSD: cproj.3,v 1.3 2011/11/29 13:17:04 drochner Exp $
.\" Copyright (c) 2001-2003 The Open Group, All Rights Reserved
.TH "CPROJ" 3P 2003 "IEEE/The Open Group" "POSIX Programmer's Manual"
.\" cproj
\fB#include <complex.h>
.br
.sp
-double complex cproj(double complex\fP \fIz\fP\fB);
+double cproj(double complex\fP \fIz\fP\fB);
.br
-float complex cprojf(float complex\fP \fIz\fP\fB);
+float cprojf(float complex\fP \fIz\fP\fB);
.br
-long double complex cprojl(long double complex\fP \fIz\fP\fB);
+long double cprojl(long double complex\fP \fIz\fP\fB);
.br
\fP
.SH DESCRIPTION
.LP
-These functions compute a projection of \fIz\fP onto the Riemann
-sphere: \fIz\fP projects to \fIz\fP, except that all
-complex infinities (even those with one infinite part and one NaN
-part) project to positive infinity on the real axis. If \fIz\fP
-has an infinite part, then \fIcproj\fP( \fIz\fP) shall be equivalent
-to:
-.sp
-.RS
-.nf
-
-\fBINFINITY + I * copysign(0.0, cimag(z))
-\fP
-.fi
-.RE
+These functions compute a projection of \fIz\fP onto the Riemann sphere:
+\fIz\fP projects to \fIz\fP , except that all complex infinities (even those
+with one infinite part and one NaN part) project to positive infinity on the
+real axis. If \fIz\fP has an infinite part, then cproj(z) shall be equivalent to:
+INFINITY + I * copysign(0.0, cimag(z))
.SH RETURN VALUE
.LP
-These functions return the value of the projection onto the
-Riemann sphere.
+These functions return the value of the projection onto the Riemann sphere.
.SH ERRORS
.LP
No errors are defined.
None.
.SH RATIONALE
.LP
-Two topologies are commonly used in complex mathematics: the complex
-plane with its continuum of infinities, and the Riemann
-sphere with its single infinity. The complex plane is better suited
-for transcendental functions, the Riemann sphere for algebraic
-functions. The complex types with their multiplicity of infinities
-provide a useful (though imperfect) model for the complex plane.
-The \fIcproj\fP() function helps model the Riemann sphere by mapping
-all infinities to one, and should be used just before any
-operation, especially comparisons, that might give spurious results
-for any of the other infinities. Note that a complex value with
-one infinite part and one NaN part is regarded as an infinity, not
-a NaN, because if one part is infinite, the complex value is
-infinite independent of the value of the other part. For the same
-reason, \fIcabs\fP()
-returns an infinity if its argument has an infinite part and a NaN
-part.
+None.
.SH FUTURE DIRECTIONS
.LP
None.
.SH SEE ALSO
.LP
-\fIcarg\fP(), \fIcimag\fP(), \fIconj\fP(), \fIcreal\fP(), the
+\fIcarg\fP(), \fIcimag\fP(), \fIconj\fP(), \fIcreal\fP() the
Base Definitions volume of IEEE\ Std\ 1003.1-2001, \fI<complex.h>\fP
.SH COPYRIGHT
Portions of this text are reprinted and reproduced in electronic form
+/* $NetBSD: cproj.c,v 1.5 2011/11/02 02:34:56 christos Exp $ */
+
/*-
- * Copyright (c) 2008 David Schultz <das@FreeBSD.ORG>
+ * Copyright (c) 2010 The NetBSD Foundation, Inc.
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
- * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
+ * THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
+ * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+ * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+ * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR CONTRIBUTORS
+ * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
*/
-#include <sys/cdefs.h>
-
#include <complex.h>
#include <math.h>
#include "../src/math_private.h"
+/*
+ * cproj(double complex z)
+ *
+ * These functions return the value of the projection (not stereographic!)
+ * onto the Riemann sphere.
+ *
+ * z projects to z, except that all complex infinities (even those with one
+ * infinite part and one NaN part) project to positive infinity on the real axis.
+ * If z has an infinite part, then cproj(z) shall be equivalent to:
+ *
+ * INFINITY + I * copysign(0.0, cimag(z))
+ */
double complex
cproj(double complex z)
{
+ double_complex w = { .z = z };
- if (!isinf(creal(z)) && !isinf(cimag(z)))
- return (z);
- else
- return (cpack(INFINITY, copysign(0.0, cimag(z))));
-}
-
-#if LDBL_MANT_DIG == 53
-__weak_reference(cproj, cprojl);
+ /*CONSTCOND*/
+ if (isinf(creal(z)) || isinf(cimag(z))) {
+#ifdef __INFINITY
+ REAL_PART(w) = HUGE_VAL;
+#else
+ REAL_PART(w) = INFINITY;
#endif
+ IMAG_PART(w) = copysign(0.0, cimag(z));
+ }
+ return (w.z);
+}
+/* $NetBSD: cprojf.c,v 1.5 2011/11/02 02:34:56 christos Exp $ */
+
/*-
- * Copyright (c) 2008 David Schultz <das@FreeBSD.ORG>
+ * Copyright (c) 2010 The NetBSD Foundation, Inc.
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
- * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
+ * THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
+ * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+ * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+ * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR CONTRIBUTORS
+ * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
*/
-#include <sys/cdefs.h>
-
#include <complex.h>
#include <math.h>
#include "../src/math_private.h"
+/*
+ * cprojf(float complex z)
+ *
+ * These functions return the value of the projection (not stereographic!)
+ * onto the Riemann sphere.
+ *
+ * z projects to z, except that all complex infinities (even those with one
+ * infinite part and one NaN part) project to positive infinity on the real axis.
+ * If z has an infinite part, then cproj(z) shall be equivalent to:
+ *
+ * INFINITY + I * copysign(0.0, cimag(z))
+ */
+
float complex
cprojf(float complex z)
{
+ float_complex w = { .z = z };
+
+ /*CONSTCOND*/
+ if (isinf(crealf(z)) || isinf(cimagf(z))) {
+#ifdef __INFINITY
+ REAL_PART(w) = HUGE_VAL;
+#else
+ REAL_PART(w) = INFINITY;
+#endif
+ IMAG_PART(w) = copysignf(0.0, cimagf(z));
+ }
- if (!isinf(crealf(z)) && !isinf(cimagf(z)))
- return (z);
- else
- return (cpackf(INFINITY, copysignf(0.0, cimagf(z))));
+ return (w.z);
}
+/* $NetBSD: cprojl.c,v 1.6 2011/11/02 02:34:56 christos Exp $ */
+
/*-
- * Copyright (c) 2008 David Schultz <das@FreeBSD.ORG>
+ * Copyright (c) 2010 The NetBSD Foundation, Inc.
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
- * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
+ * THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
+ * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+ * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+ * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR CONTRIBUTORS
+ * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
*/
-#include <sys/cdefs.h>
-
#include <complex.h>
#include <math.h>
#include "../src/math_private.h"
+/*
+ * cprojl(long double complex z)
+ *
+ * These functions return the value of the projection (not stereographic!)
+ * onto the Riemann sphere.
+ *
+ * z projects to z, except that all complex infinities (even those with one
+ * infinite part and one NaN part) project to positive infinity on the real axis.
+ * If z has an infinite part, then cproj(z) shall be equivalent to:
+ *
+ * INFINITY + I * copysign(0.0, cimag(z))
+ */
long double complex
cprojl(long double complex z)
{
+ long_double_complex w = { .z = z };
+
+ /*CONSTCOND*/
+ if (isinf(creall(z)) || isinf(cimagl(z))) {
+#ifdef __INFINITY
+ REAL_PART(w) = HUGE_VAL;
+#else
+ REAL_PART(w) = INFINITY;
+#endif
+ IMAG_PART(w) = copysignl(0.0, cimagl(z));
+ }
- if (!isinf(creall(z)) && !isinf(cimagl(z)))
- return (z);
- else
- return (cpackl(INFINITY, copysignl(0.0, cimagl(z))));
+ return (w.z);
}
-.\" $NetBSD: creal.3,v 1.1 2008/02/20 09:55:38 drochner Exp $
+.\" $NetBSD: creal.3,v 1.3 2010/09/15 18:40:27 wiz Exp $
.\" Copyright (c) 2001-2003 The Open Group, All Rights Reserved
.TH "CREAL" 3P 2003 "IEEE/The Open Group" "POSIX Programmer's Manual"
.\" creal
.br
float crealf(float complex\fP \fIz\fP\fB);
.br
+long double creall(long double complex\fP \fIz\fP\fB);
+.br
\fP
.SH DESCRIPTION
.LP
-/*-
- * Copyright (c) 2004 Stefan Farfeleder
- * All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- * notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- * notice, this list of conditions and the following disclaimer in the
- * documentation and/or other materials provided with the distribution.
- *
- * THIS SOFTWARE IS PROVIDED BY AUTHOR AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED. IN NO EVENT SHALL AUTHOR OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- *
- * $FreeBSD$
+/* $NetBSD: creal.c,v 1.2 2010/09/15 16:11:29 christos Exp $ */
+
+/*
+ * Written by Matthias Drochner <drochner@NetBSD.org>.
+ * Public domain.
*/
#include <complex.h>
+#include "../src/math_private.h"
double
creal(double complex z)
{
- return z;
+ double_complex w = { .z = z };
+
+ return (REAL_PART(w));
}
-/*-
- * Copyright (c) 2004 Stefan Farfeleder
- * All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- * notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- * notice, this list of conditions and the following disclaimer in the
- * documentation and/or other materials provided with the distribution.
- *
- * THIS SOFTWARE IS PROVIDED BY AUTHOR AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED. IN NO EVENT SHALL AUTHOR OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- *
- * $FreeBSD$
+/* $NetBSD: crealf.c,v 1.2 2010/09/15 16:11:29 christos Exp $ */
+
+/*
+ * Written by Matthias Drochner <drochner@NetBSD.org>.
+ * Public domain.
*/
#include <complex.h>
+#include "../src/math_private.h"
float
crealf(float complex z)
{
- return z;
+ float_complex w = { .z = z };
+
+ return (REAL_PART(w));
}
+/* $NetBSD: creall.c,v 1.3 2010/09/20 16:55:20 christos Exp $ */
+
/*-
- * Copyright (c) 2004 Stefan Farfeleder
+ * Copyright (c) 2010 The NetBSD Foundation, Inc.
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
- * THIS SOFTWARE IS PROVIDED BY AUTHOR AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED. IN NO EVENT SHALL AUTHOR OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- *
- * $FreeBSD$
+ * THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
+ * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+ * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+ * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR CONTRIBUTORS
+ * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
*/
#include <complex.h>
+#include "../src/math_private.h"
+/*
+ * creall(long double complex z)
+ * This function returns the real part value of z.
+ */
long double
creall(long double complex z)
{
- return z;
+ long_double_complex w = { .z = z };
+
+ return (REAL_PART(w));
}
-/* $NetBSD: ctan.c,v 1.1 2007/08/20 16:01:37 drochner Exp $ */
+/* $NetBSD: ctan.c,v 1.2 2011/07/03 06:45:24 mrg Exp $ */
/*-
* Copyright (c) 2007 The NetBSD Foundation, Inc.
* POSSIBILITY OF SUCH DAMAGE.
*/
-
#include <complex.h>
#include <math.h>
+#include <float.h>
#include "cephes_subr.h"
-#define MAXNUM 1.0e308
+#define MAXNUM DBL_MAX
double complex
ctan(double complex z)
-/* $NetBSD: ctanf.c,v 1.1 2007/08/20 16:01:38 drochner Exp $ */
+/* $NetBSD: ctanf.c,v 1.2 2011/07/03 06:45:24 mrg Exp $ */
/*-
* Copyright (c) 2007 The NetBSD Foundation, Inc.
* POSSIBILITY OF SUCH DAMAGE.
*/
-
#include <complex.h>
#include <math.h>
+#include <float.h>
#include "cephes_subrf.h"
-#define MAXNUMF 1.0e38f
+#define MAXNUMF FLT_MAX
float complex
ctanf(float complex z)
+/* $NetBSD: nan.c,v 1.2 2008/04/28 20:23:01 martin Exp $ */
+
/*-
* Copyright (c) 2006 The NetBSD Foundation, Inc.
* All rights reserved.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
- * 3. All advertising materials mentioning features or use of this software
- * must display the following acknowledgement:
- * This product includes software developed by the NetBSD
- * Foundation, Inc. and its contributors.
- * 4. Neither the name of The NetBSD Foundation nor the names of its
- * contributors may be used to endorse or promote products derived
- * from this software without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
* ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
* POSSIBILITY OF SUCH DAMAGE.
- *
- * $NetBSD: nan.c,v 1.1 2006/03/15 22:07:09 kleink Exp $
- * $DragonFly: src/lib/libm/gen/nan.c,v 1.1 2007/06/17 17:46:01 pavalos Exp $
*/
#include <assert.h>
+/* $NetBSD: nanf.c,v 1.2 2008/04/28 20:23:01 martin Exp $ */
+
/*-
* Copyright (c) 2006 The NetBSD Foundation, Inc.
* All rights reserved.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
- * 3. All advertising materials mentioning features or use of this software
- * must display the following acknowledgement:
- * This product includes software developed by the NetBSD
- * Foundation, Inc. and its contributors.
- * 4. Neither the name of The NetBSD Foundation nor the names of its
- * contributors may be used to endorse or promote products derived
- * from this software without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
* ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
* POSSIBILITY OF SUCH DAMAGE.
- *
- * $NetBSD: nanf.c,v 1.1 2006/03/15 22:07:09 kleink Exp $
- * $DragonFly: src/lib/libm/gen/nanf.c,v 1.1 2007/06/17 17:46:01 pavalos Exp $
*/
#define NAN_FUNCTION nanf
#define NAN_TYPE float
-#define NAN_STRTOD strtod /* XXX should be strtof */
+#define NAN_STRTOD strtof
#include "nan.c"
+/* $NetBSD: nanl.c,v 1.2 2008/04/28 20:23:01 martin Exp $ */
+
/*-
* Copyright (c) 2006 The NetBSD Foundation, Inc.
* All rights reserved.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
- * 3. All advertising materials mentioning features or use of this software
- * must display the following acknowledgement:
- * This product includes software developed by the NetBSD
- * Foundation, Inc. and its contributors.
- * 4. Neither the name of The NetBSD Foundation nor the names of its
- * contributors may be used to endorse or promote products derived
- * from this software without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
* ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
* POSSIBILITY OF SUCH DAMAGE.
- *
- * $NetBSD: nanl.c,v 1.1 2006/03/15 22:07:09 kleink Exp $
- * $DragonFly: src/lib/libm/gen/nanl.c,v 1.1 2007/06/17 17:46:01 pavalos Exp $
*/
#define NAN_FUNCTION nanl
#define NAN_TYPE long double
-#define NAN_STRTOD strtod /* XXX should be strtold */
+#define NAN_STRTOD strtold
#include "nan.c"
MAN+= acos.3 acosh.3 asin.3 asinh.3 atan.3 atan2.3 atanh.3 ceil.3 \
cos.3 cosh.3 erf.3 exp.3 fabs.3 fdim.3 feclearexcept.3 feenableexcept.3 \
- fegetenv.3 fegetround.3 fenv.3 floor.3 fmax.3 fmod.3 hypot.3 ieee.3 \
+ fegetenv.3 fegetround.3 fenv.3 floor.3 fmax.3 fmod.3 hypot.3 \
ieee_test.3 j0.3 lgamma.3 math.3 rint.3 round.3 sin.3 \
sinh.3 sqrt.3 tan.3 tanh.3 trunc.3
+MAN+= copysign.3 frexp.3 ilogb.3 log.3 scalbn.3 pow.3 nextafter.3 remainder.3 \
+ finite.3
MLINKS+=acos.3 acosf.3
MLINKS+=acosh.3 acoshf.3
MLINKS+=cos.3 cosf.3
MLINKS+=cosh.3 coshf.3
MLINKS+=erf.3 erff.3 erf.3 erfc.3 erf.3 erfcf.3
-MLINKS+=exp.3 expf.3 exp.3 expm1.3 exp.3 expm1f.3 exp.3 log.3 exp.3 logf.3 \
- exp.3 log10.3 exp.3 log10f.3 exp.3 log1p.3 exp.3 log1pf.3 exp.3 log2.3 \
- exp.3 log2f.3 exp.3 pow.3 exp.3 powf.3
-MLINKS+=fabs.3 fabsf.3
+MLINKS+=log.3 logf.3 log.3 log10.3 log.3 log10f.3 log.3 log1p.3 log.3 log1pf.3 \
+ log.3 log2.3 log.3 log2f.3
+MLINKS+=exp.3 expf.3 exp.3 exp2.3 exp.3 exp2f.3 exp.3 exp2l.3
+
+MLINKS+=fabs.3 fabsf.3 fabs.3 fabsl.3
MLINKS+=fdim.3 fdimf.3 fdim.3 fdiml.3
MLINKS+=feclearexcept.3 fegetexceptflag.3 feclearexcept.3 feraiseexcept.3 \
feclearexcept.3 fesetexceptflag.3 feclearexcept.3 fetestexcept.3
fmax.3 fmin.3 fmax.3 fminf.3 fmax.3 fminl.3
MLINKS+=fmod.3 fmodf.3
MLINKS+=hypot.3 hypotf.3 hypot.3 cabs.3 hypot.3 cabsf.3
-MLINKS+=ieee.3 copysign.3 ieee.3 copysignf.3 ieee.3 finite.3 ieee.3 finitef.3 \
- ieee.3 ilogb.3 ieee.3 ilogbf.3 ieee.3 nextafter.3 ieee.3 nextafterf.3 \
- ieee.3 remainder.3 ieee.3 remainderf.3 ieee.3 scalbn.3 \
- ieee.3 scalbnf.3
MLINKS+=ieee_test.3 logb.3 ieee_test.3 logbf.3 ieee_test.3 scalb.3 \
ieee_test.3 scalbf.3 ieee_test.3 significand.3 \
ieee_test.3 significandf.3
MLINKS+=tan.3 tanf.3
MLINKS+=tanh.3 tanhf.3
MLINKS+=trunc.3 truncf.3
+MLINKS+=copysign.3 copysignf.3 copysign.3 copysignl.3
+MLINKS+=scalbn.3 scalbnf.3 scalbn.3 scalbnl.3
+MLINKS+=pow.3 powf.3
+MLINKS+=frexp.3 frexpf.3 frexp.3 frexpl.3
+MLINKS+=ilogb.3 ilogbf.3 ilogb.3 ilogbl.3
+MLINKS+=nextafter.3 nextafterf.3 nextafter.3 nextafterl.3 nextafter.3 nexttoward.3
+MLINKS+=remainder.3 remainderf.3 remainder.3 remquo.3 remainder.3 remquof.3
+MLINKS+=finite.3 finitef.3
.\" SUCH DAMAGE.
.\"
.\" from: @(#)ceil.3 5.1 (Berkeley) 5/2/91
-.\" $NetBSD: ceil.3,v 1.17 2003/08/07 16:44:47 agc Exp $
-.\" $DragonFly: src/lib/libm/man/ceil.3,v 1.1 2005/07/26 21:15:20 joerg Exp $
+.\" $NetBSD: ceil.3,v 1.19 2011/09/18 05:33:13 jruoho Exp $
.\"
-.Dd March 10, 1994
+.Dd September 18, 2011
.Dt CEIL 3
.Os
.Sh NAME
.Nm ceil ,
-.Nm ceilf
-.Nd "round to smallest integral value greater than or equal to x"
+.Nm ceilf ,
+.Nm floor ,
+.Nm floorf
+.Nd ceiling and floor
.Sh LIBRARY
.Lb libm
.Sh SYNOPSIS
.Fn ceil "double x"
.Ft float
.Fn ceilf "float x"
+.Ft double
+.Fn floor "double x"
+.Ft float
+.Fn floorf "float x"
.Sh DESCRIPTION
The
.Fn ceil
functions return the smallest integral value
greater than or equal to
.Fa x .
+Conversely, the
+.Fn floor
+and
+.Fn floorf
+functions return the largest integral value
+less than or equal to
+.Fa x .
.Sh SEE ALSO
.Xr abs 3 ,
.Xr fabs 3 ,
-.Xr floor 3 ,
-.Xr ieee 3 ,
.Xr math 3 ,
+.Xr nextafter 3 ,
.Xr rint 3
.Sh STANDARDS
-The
-.Fn ceil
-function conforms to
-.St -ansiC .
+The described functions conform to
+.St -isoC-99 .
--- /dev/null
+.\" $NetBSD: copysign.3,v 1.1 2011/04/13 04:57:10 jruoho Exp $
+.\"
+.\" Copyright (c) 2011 Jukka Ruohonen <jruohonen@iki.fi>
+.\" All rights reserved.
+.\"
+.\" Redistribution and use in source and binary forms, with or without
+.\" modification, are permitted provided that the following conditions
+.\" are met:
+.\" 1. Redistributions of source code must retain the above copyright
+.\" notice, this list of conditions and the following disclaimer.
+.\" 2. Redistributions in binary form must reproduce the above copyright
+.\" notice, this list of conditions and the following disclaimer in the
+.\" documentation and/or other materials provided with the distribution.
+.\"
+.\" THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
+.\" ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+.\" TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+.\" PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR CONTRIBUTORS
+.\" BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+.\" CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+.\" SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+.\" INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+.\" CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+.\" ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+.\" POSSIBILITY OF SUCH DAMAGE.
+.\"
+.Dd April 13, 2011
+.Dt COPYSIGN 3
+.Os
+.Sh NAME
+.Nm copysign ,
+.Nm copysignf ,
+.Nm copysignl
+.Nd functions to manipulate signs
+.Sh LIBRARY
+.Lb libm
+.Sh SYNOPSIS
+.In math.h
+.Ft double
+.Fn copysign "double x" "double y"
+.Ft float
+.Fn copysignf "float x" "float y"
+.Ft long double
+.Fn copysignl "long double x" "long double y"
+.Sh DESCRIPTION
+The
+.Fn copysign ,
+.Fn copysignf ,
+and
+.Fn copysignl
+functions return a value whose absolute value matches
+.Fa x ,
+but whose sign bit is taken from
+.Fa y .
+.Sh RETURN VALUES
+Upon successful completion,
+all three functions return a value with the magnitude of
+.Fa x
+and the sign of
+.Fa y .
+If
+.Fa x
+is
+\*(Na ,
+the functions return a
+\*(Na
+with the sign of
+.Fa y .
+.Sh SEE ALSO
+.Xr math 3 ,
+.Xr signbit 3
+.Sh STANDARDS
+The described functions conform to
+.St -isoC-99 .
+.\"
+.\" XXX: Verify this.
+.\"
+.\" The functions are also recommended by
+.\" .St -ieee754
+.\"
+.\" .Sh HISTORY
+.\"
+.\" XXX: Fill this.
+.\"
+.\" These functions first appeared in ???.
+.\"
+.Sh CAVEATS
+Note that on implementations that represent a signed zero
+but do not treat negative zero consistently in arithmetic operations,
+these functions may regard the sign of zero as positive.
.\" SUCH DAMAGE.
.\"
.\" from: @(#)exp.3 6.12 (Berkeley) 7/31/91
-.\" $NetBSD: exp.3,v 1.23 2005/07/21 12:58:22 wiz Exp $
-.\" $DragonFly: src/lib/libm/man/exp.3,v 1.2 2007/06/16 22:26:53 pavalos Exp $
+.\" $FreeBSD: src/lib/msun/man/exp.3,v 1.24 2008/01/18 21:43:00 das Exp $
+.\" $NetBSD: exp.3,v 1.28 2011/09/17 10:52:52 jruoho Exp $
.\"
-.Dd July 21, 2005
+.Dd September 13, 2011
.Dt EXP 3
.Os
.Sh NAME
.Nm exp ,
.Nm expf ,
+.\" The sorting error is intentional. exp and expf should be adjacent.
+.Nm exp2 ,
+.Nm exp2f ,
+.\" .Nm exp2l ,
.Nm expm1 ,
.Nm expm1f ,
-.Nm log ,
-.Nm logf ,
-.Nm log2 ,
-.Nm log2f ,
-.Nm log10 ,
-.Nm log10f ,
-.Nm log1p ,
-.Nm log1pf ,
-.Nm pow ,
-.Nm powf
-.Nd exponential, logarithm, power functions
+.Nd exponential functions
.Sh LIBRARY
.Lb libm
.Sh SYNOPSIS
.Ft float
.Fn expf "float x"
.Ft double
-.Fn expm1 "double x"
-.Ft float
-.Fn expm1f "float x"
-.Ft double
-.Fn log "double x"
-.Ft float
-.Fn logf "float x"
-.Ft double
-.Fn log2 "double x"
+.Fn exp2 "double x"
.Ft float
-.Fn log2f "float x"
+.Fn exp2f "float x"
+.\" .Ft long double
+.\" .Fn exp2l "long double x"
.Ft double
-.Fn log10 "double x"
-.Ft float
-.Fn log10f "float x"
-.Ft double
-.Fn log1p "double x"
-.Ft float
-.Fn log1pf "float x"
-.Ft double
-.Fn pow "double x" "double y"
+.Fn expm1 "double x"
.Ft float
-.Fn powf "float x" "float y"
+.Fn expm1f "float x"
.Sh DESCRIPTION
The
.Fn exp
-function computes the exponential value of the given argument
-.Fa x .
-.Pp
-The
-.Fn expm1
-function computes the value exp(x)\-1 accurately even for tiny argument
-.Fa x .
-.Pp
-The
-.Fn log
-function computes the value of the natural logarithm of argument
+and the
+.Fn expf
+functions compute the base
+.Ms e
+exponential value of the given argument
.Fa x .
.Pp
The
-.Fn log10
-function computes the value of the logarithm of argument
-.Fa x
-to base 10.
-.Pp
-The
-.Fn log1p
-function computes
-the value of log(1+x) accurately even for tiny argument
+.Fn exp2 ,
+and
+.Fn exp2f
+.\" .Fn exp2f ,
+.\" and
+.\" .Fn exp2l
+functions compute the base 2 exponential of the given argument
.Fa x .
.Pp
The
-.Fn log2
+.Fn expm1
and the
-.Fn log2f
-functions compute the value of the logarithm of argument
-.Fa x
-to base 2.
-.Pp
-The
-.Fn pow
-computes the value
-of
-.Ar x
-to the exponent
-.Ar y .
+.Fn expm1f
+functions computes the value exp(x)\-1 accurately even for tiny argument
+.Fa x .
.Sh RETURN VALUES
These functions will return the appropriate computation unless an error
occurs or an argument is out of range.
The functions
-.Fn exp ,
-.Fn expm1
+.Fn exp
and
-.Fn pow
+.Fn expm1
detect if the computed value will overflow,
set the global variable
.Va errno
.Er ERANGE
and cause a reserved operand fault on a
.Tn VAX .
-The function
-.Fn pow x y
-checks to see if
-.Fa x
-\*[Lt] 0 and
-.Fa y
-is not an integer, in the event this is true,
-the global variable
-.Va errno
-is set to
-.Er EDOM
-and on the
-.Tn VAX
-generate a reserved operand fault.
-On a
-.Tn VAX ,
-.Va errno
-is set to
-.Er EDOM
-and the reserved operand is returned
-by log unless
-.Fa x
-\*[Gt] 0, by
-.Fn log1p
-unless
-.Fa x
-\*[Gt] \-1.
-.Sh ERRORS
-exp(x), log(x), expm1(x) and log1p(x) are accurate to within
-an
-.Em ulp ,
-and log10(x) to within about 2
-.Em ulps ;
-an
-.Em ulp
-is one
-.Em Unit
-in the
-.Em Last
-.Em Place .
-The error in
-.Fn pow x y
-is below about 2
-.Em ulps
-when its
-magnitude is moderate, but increases as
-.Fn pow x y
-approaches
-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
-.Tn "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
-.Em ulps
-for
-.Tn "VAX D" ,
-300
-.Em ulps
-for
-.Tn IEEE
-754 Double.
-Moderate values of
-.Fn pow
-are accurate enough that
-.Fn pow integer integer
-is exact until it is bigger than 2**56 on a
-.Tn VAX ,
-2**53 for
-.Tn IEEE
-754.
-.Sh NOTES
-The functions exp(x)\-1 and log(1+x) are called
-expm1 and logp1 in
-.Tn BASIC
-on the Hewlett\-Packard
-.Tn HP Ns \-71B
-and
-.Tn APPLE
-Macintosh,
-.Tn EXP1
-and
-.Tn LN1
-in Pascal, exp1 and log1 in C
-on
-.Tn APPLE
-Macintoshes, where they have been provided to make
-sure financial calculations of ((1+x)**n\-1)/x, namely
-expm1(n\(**log1p(x))/x, will be accurate when x is tiny.
-They also provide accurate inverse hyperbolic functions.
-.Pp
-The function
-.Fn pow x 0
-returns x**0 = 1 for all x including x = 0,
-.if n \
-Infinity
-.if t \
-\(if
-(not found on a
-.Tn VAX ) ,
-and
-.Em NaN
-(the reserved
-operand on a
-.Tn VAX ) .
-Previous implementations of pow may
-have defined x**0 to be undefined in some or all of these
-cases.
-Here are reasons for returning x**0 = 1 always:
-.Bl -enum -width indent
-.It
-Any program that already tests whether x is zero (or
-infinite or \*(Na) before computing x**0 cannot care
-whether 0**0 = 1 or not.
-Any program that depends
-upon 0**0 to be invalid is dubious anyway since that
-expression's meaning and, if invalid, its consequences
-vary from one computer system to another.
-.It
-Some Algebra texts (e.g. Sigler's) define x**0 = 1 for
-all x, including x = 0.
-This is compatible with the convention that accepts a[0]
-as the value of polynomial
-.Bd -literal -offset indent
-p(x) = a[0]\(**x**0 + a[1]\(**x**1 + a[2]\(**x**2 +...+ a[n]\(**x**n
-.Ed
-.Pp
-at x = 0 rather than reject a[0]\(**0**0 as invalid.
-.It
-Analysts will accept 0**0 = 1 despite that x**y can
-approach anything or nothing as x and y approach 0
-independently.
-The reason for setting 0**0 = 1 anyway is this:
-.Bd -filled -offset indent
-If x(z) and y(z) are
-.Em any
-functions analytic (expandable
-in power series) in z around z = 0, and if there
-x(0) = y(0) = 0, then x(z)**y(z) \(-\*[Gt] 1 as z \(-\*[Gt] 0.
-.Ed
-.It
-If 0**0 = 1, then
-.if n \
-infinity**0 = 1/0**0 = 1 too; and
-.if t \
-\(if**0 = 1/0**0 = 1 too; and
-then \*(Na**0 = 1 too because x**0 = 1 for all finite
-and infinite x, i.e., independently of x.
-.El
.Sh SEE ALSO
.Xr math 3
.Sh STANDARDS
The
-.Fn exp ,
-.Fn log ,
-.Fn log10
-and
-.Fn pow
+.Fn exp
functions conform to
.St -ansiC .
-.Sh HISTORY
-A
-.Fn exp ,
-.Fn log
+The
+.Fn exp2 ,
+.Fn exp2f ,
+.Fn expf ,
+.Fn expm1 ,
and
-.Fn pow
-functions
-appeared in
+.Fn expm1f
+functions conform to
+.St -isoC-99 .
+.Sh HISTORY
+The
+.Fn exp
+functions appeared in
.At v6 .
-A
-.Fn log10
-function
-appeared in
-.At v7 .
The
-.Fn log1p
-and
.Fn expm1
-functions appeared in
+function appeared in
.Bx 4.3 .
.\" SUCH DAMAGE.
.\"
.\" from: @(#)fabs.3 5.1 (Berkeley) 5/2/91
-.\" $NetBSD: fabs.3,v 1.14 2003/08/07 16:44:47 agc Exp $
-.\" $DragonFly: src/lib/libm/man/fabs.3,v 1.1 2005/07/26 21:15:20 joerg Exp $
+.\" $NetBSD: fabs.3,v 1.15 2011/09/13 07:11:43 njoly Exp $
.\"
.Dd May 2, 1991
.Dt FABS 3
.Xr abs 3 ,
.Xr ceil 3 ,
.Xr floor 3 ,
-.Xr ieee 3 ,
.Xr math 3 ,
.Xr rint 3
.Sh STANDARDS
The macro
.Dv FE_DFL_ENV
expands to a pointer to the default environment.
-.Sh CAVEATS
-The FENV_ACCESS pragma can be enabled with
-.Dl "#pragma STDC FENV_ACCESS ON"
-and disabled with the
-.Dl "#pragma STDC FENV_ACCESS OFF"
-directive.
-This lexically-scoped annotation tells the compiler that the program
-may access the floating-point environment, so optimizations that would
-violate strict IEEE-754 semantics are disabled.
-If execution reaches a block of code for which
-.Dv FENV_ACCESS
-is off, the floating-point environment will become undefined.
.Sh EXAMPLES
The following routine computes the square root function.
It explicitly raises an invalid exception on appropriate inputs using
double x = 1.0;
fenv_t env;
- if (isnan(n) || n < 0.0) {
+ if (isnan(n) || n \*[Lt] 0.0) {
feraiseexcept(FE_INVALID);
return (NAN);
}
if (isinf(n) || n == 0.0)
return (n);
- feholdexcept(&env);
- while (fabs((x * x) - n) > DBL_EPSILON * 2 * x)
+ feholdexcept(\*[Am]env);
+ while (fabs((x * x) - n) \*[Gt] DBL_EPSILON * 2 * x)
x = (x / 2) + (n / (2 * x));
if (x * x == n)
feclearexcept(FE_INEXACT);
- feupdateenv(&env);
+ feupdateenv(\*[Am]env);
return (x);
}
.Ed
.In ieeefp.h
and documented in
.Xr fpgetround 3 .
+.Sh CAVEATS
+The FENV_ACCESS pragma can be enabled with
+.Dl "#pragma STDC FENV_ACCESS ON"
+and disabled with the
+.Dl "#pragma STDC FENV_ACCESS OFF"
+directive.
+This lexically-scoped annotation tells the compiler that the program
+may access the floating-point environment, so optimizations that would
+violate strict IEEE-754 semantics are disabled.
+If execution reaches a block of code for which
+.Dv FENV_ACCESS
+is off, the floating-point environment will become undefined.
.Sh BUGS
The
.Dv FENV_ACCESS
-.\" Copyright (c) 1991 The Regents of the University of California.
+.\" Copyright (c) 1985, 1991 Regents of the University of California.
.\" All rights reserved.
.\"
.\" Redistribution and use in source and binary forms, with or without
.\" OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
.\" SUCH DAMAGE.
.\"
-.\" from: @(#)fabs.3 5.1 (Berkeley) 5/2/91
-.\" $NetBSD: fabs.3,v 1.14 2003/08/07 16:44:47 agc Exp $
-.\" $DragonFly: src/lib/libm/man/fabs.3,v 1.1 2005/07/26 21:15:20 joerg Exp $
+.\" from: @(#)ieee.3 6.4 (Berkeley) 5/6/91
+.\" $NetBSD: finite.3,v 1.2 2011/08/06 11:09:22 wiz Exp $
.\"
-.Dd May 2, 1991
-.Dt FABS 3
+.Dd July 28, 2011
+.Dt FINITE 3
.Os
.Sh NAME
-.Nm fabs ,
-.Nm fabsf
-.Nd floating-point absolute value function
+.Nm finite ,
+.Nm finitef
+.Nd tests for finite values
.Sh LIBRARY
.Lb libm
.Sh SYNOPSIS
.In math.h
-.Ft double
-.Fn fabs "double x"
-.Ft float
-.Fn fabsf "float x"
+.Ft int
+.Fn finite "double x"
+.Ft int
+.Fn finitef "float x"
.Sh DESCRIPTION
The
-.Fn fabs
-and
-.Fn fabsf
-functions compute the absolute value of a floating-point number
-.Fa x .
-.Sh RETURN VALUES
-The
-.Fn fabs
-function returns the absolute value of
-.Fa x .
+.Fn finite
+function returns the value 1 when
+.Bd -ragged -offset indent
+\-\*(If \*(Lt
+.Fa x
+\*(Lt +\*(If.
+.Ed
+.Pp
+Otherwise a zero is returned
+(that is,
+.Pf \*(Ba Ns Fa x Ns \*(Ba
+= \*(If or
+.Fa x
+is \*(Na).
.Sh SEE ALSO
-.Xr abs 3 ,
-.Xr ceil 3 ,
-.Xr floor 3 ,
-.Xr ieee 3 ,
-.Xr math 3 ,
-.Xr rint 3
+.Xr isfinite 3 ,
+.Xr math 3
.Sh STANDARDS
+The described functions conform to
+.St -ieee754 .
+Note that unlike
+.Xr isfinite 3 ,
+neither function is present in the
+.Dv ISO
+C-language standards or in the
+.Dv IEEE
+.Dv POSIX
+standards.
+.Sh HISTORY
The
-.Fn fabs
-function conforms to
-.St -ansiC .
+.Nm finite
+and
+.Fn finitef
+functions first appeared in
+.Bx 4.3 .
-.\" Copyright (c) 1991 The Regents of the University of California.
-.\" All rights reserved.
+.\" $NetBSD: frexp.3,v 1.2 2010/04/29 08:35:03 joerg Exp $
+.\"
+.\" Copyright (c) 1991, 1993
+.\" The Regents of the University of California. All rights reserved.
+.\"
+.\" This code is derived from software contributed to Berkeley by
+.\" the American National Standards Committee X3, on Information
+.\" Processing Systems.
.\"
.\" Redistribution and use in source and binary forms, with or without
.\" modification, are permitted provided that the following conditions
.\" OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
.\" SUCH DAMAGE.
.\"
-.\" from: @(#)tanh.3 5.1 (Berkeley) 5/2/91
-.\" $NetBSD: tanh.3,v 1.12 2003/08/07 16:44:49 agc Exp $
-.\" $DragonFly: src/lib/libm/man/tanh.3,v 1.1 2005/07/26 21:15:20 joerg Exp $
+.\" @(#)frexp.3 8.1 (Berkeley) 6/4/93
.\"
-.Dd May 2, 1991
-.Dt TANH 3
+.Dd March 21, 2006
+.Dt FREXP 3
.Os
.Sh NAME
-.Nm tanh ,
-.Nm tanhf
-.Nd hyperbolic tangent function
+.Nm frexp
+.Nd convert floating-point number to fractional and integral components
.Sh LIBRARY
.Lb libm
.Sh SYNOPSIS
.In math.h
.Ft double
-.Fn tanh "double x"
+.Fn frexp "double value" "int *exp"
.Ft float
-.Fn tanhf "float x"
+.Fn frexpf "float value" "int *exp"
.Sh DESCRIPTION
The
-.Fn tanh
-and
-.Fn tanhf
-functions compute the hyperbolic tangent of
-.Fa x .
-For a discussion of error due to roundoff, see
-.Xr math 3 .
+.Fn frexp
+function breaks a floating-point number into a normalized
+fraction and an integral power of 2.
+It stores the integer in the
+.Em int
+object pointed to by
+.Fa exp .
.Sh RETURN VALUES
The
-.Fn tanh
-function returns the hyperbolic tangent value.
+.Fn frexp
+function returns the value
+.Em x ,
+such that
+.Em x
+is a
+.Em double
+with magnitude in the interval [1/2, 1) or zero, and
+.Fa value
+equals
+.Em x
+times 2 raised to the power
+.Fa *exp .
+If
+.Fa value
+is zero, both parts of the result are zero.
.Sh SEE ALSO
-.Xr acos 3 ,
-.Xr asin 3 ,
-.Xr atan 3 ,
-.Xr atan2 3 ,
-.Xr cos 3 ,
-.Xr cosh 3 ,
+.Xr ldexp 3 ,
.Xr math 3 ,
-.Xr sin 3 ,
-.Xr sinh 3 ,
-.Xr tan 3
+.Xr modf 3
.Sh STANDARDS
The
-.Fn tanh
+.Fn frexp
function conforms to
.St -ansiC .
+++ /dev/null
-.\" Copyright (c) 1985, 1991 Regents of the University of California.
-.\" All rights reserved.
-.\"
-.\" Redistribution and use in source and binary forms, with or without
-.\" modification, are permitted provided that the following conditions
-.\" are met:
-.\" 1. Redistributions of source code must retain the above copyright
-.\" notice, this list of conditions and the following disclaimer.
-.\" 2. Redistributions in binary form must reproduce the above copyright
-.\" notice, this list of conditions and the following disclaimer in the
-.\" documentation and/or other materials provided with the distribution.
-.\" 3. Neither the name of the University nor the names of its contributors
-.\" may be used to endorse or promote products derived from this software
-.\" without specific prior written permission.
-.\"
-.\" THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
-.\" ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
-.\" IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
-.\" ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
-.\" FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
-.\" DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
-.\" OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
-.\" HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
-.\" LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
-.\" OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
-.\" SUCH DAMAGE.
-.\"
-.\" from: @(#)ieee.3 6.4 (Berkeley) 5/6/91
-.\" $NetBSD: ieee.3,v 1.19 2003/08/07 16:44:48 agc Exp $
-.\" $DragonFly: src/lib/libm/man/ieee.3,v 1.2 2006/02/28 02:25:10 swildner Exp $
-.\"
-.Dd February 25, 1994
-.Dt IEEE 3
-.Os
-.Sh NAME
-.Nm copysign ,
-.Nm copysignf ,
-.Nm finite ,
-.Nm finitef ,
-.Nm ilogb ,
-.Nm ilogbf ,
-.Nm nextafter ,
-.Nm nextafterf ,
-.Nm remainder ,
-.Nm remainderf ,
-.Nm scalbn ,
-.Nm scalbnf
-.Nd functions for IEEE arithmetic
-.Sh LIBRARY
-.Lb libm
-.Sh SYNOPSIS
-.In math.h
-.Ft double
-.Fn copysign "double x" "double y"
-.Ft float
-.Fn copysignf "float x" "float y"
-.Ft int
-.Fn finite "double x"
-.Ft int
-.Fn finitef "float x"
-.Ft int
-.Fn ilogb "double x"
-.Ft int
-.Fn ilogbf "float x"
-.Ft double
-.Fn nextafter "double x" "double y"
-.Ft float
-.Fn nextafterf "float x" "float y"
-.Ft double
-.Fn remainder "double x" "double y"
-.Ft float
-.Fn remainderf "float x" "float y"
-.Ft double
-.Fn scalbn "double x" "int n"
-.Ft float
-.Fn scalbnf "float x" "int n"
-.Sh DESCRIPTION
-These functions are required or recommended by
-.St -ieee754 .
-.Pp
-.Fn copysign
-returns
-.Fa x
-with its sign changed to
-.Fa y Ns 's .
-.Pp
-.Fn finite
-returns the value 1 just when
-\-\*(If \*(Lt
-.Fa x
-\*(Lt +\*(If;
-otherwise a
-zero is returned
-(when
-.Pf \\*(Ba Ns Fa x Ns \\*(Ba
-= \*(If or
-.Fa x
-is \*(Na).
-.Pp
-.Fn ilogb
-returns
-.Fa x Ns 's exponent
-.Fa n ,
-in integer format.
-.Fn ilogb \*(Pm\*(If
-returns
-.Dv INT_MAX
-and
-.Fn ilogb 0
-returns
-.Dv INT_MIN .
-.Pp
-.Fn nextafter
-returns the next machine representable number from
-.Fa x
-in direction
-.Fa y .
-.Pp
-.Fn remainder
-returns the remainder
-.Fa r
-:=
-.Fa x
-\-
-.Fa n\(**y
-where
-.Fa n
-is the integer nearest the exact value of
-.Bk -words
-.Fa x Ns / Ns Fa y ;
-.Ek
-moreover if
-.Pf \\*(Ba Fa n
-\-
-.Sm off
-.Fa x No / Fa y No \\*(Ba
-.Sm on
-=
-\(12
-then
-.Fa n
-is even.
-Consequently the remainder is computed exactly and
-.Sm off
-.Pf \\*(Ba Fa r No \\*(Ba
-.Sm on
-\*(Le
-.Sm off
-.Pf \\*(Ba Fa y No \\*(Ba/2 .
-.Sm on
-But
-.Fn remainder x 0
-and
-.Fn remainder \*(If 0
-are invalid operations that produce a \*(Na.
-.Pp
-.Fn scalbn
-returns
-.Fa x Ns \(**(2** Ns Fa n )
-computed by exponent manipulation.
-.Sh SEE ALSO
-.Xr math 3
-.Sh STANDARDS
-.St -ieee754
-.Sh HISTORY
-The
-.Nm ieee
-functions appeared in
-.Bx 4.3 .
--- /dev/null
+.\" $NetBSD: ilogb.3,v 1.3 2011/08/02 10:15:03 wiz Exp $
+.\"
+.\" Copyright (c) 2011 Jukka Ruohonen <jruohonen@iki.fi>
+.\" All rights reserved.
+.\"
+.\" Redistribution and use in source and binary forms, with or without
+.\" modification, are permitted provided that the following conditions
+.\" are met:
+.\" 1. Redistributions of source code must retain the above copyright
+.\" notice, this list of conditions and the following disclaimer.
+.\" 2. Redistributions in binary form must reproduce the above copyright
+.\" notice, this list of conditions and the following disclaimer in the
+.\" documentation and/or other materials provided with the distribution.
+.\"
+.\" THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
+.\" ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+.\" TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+.\" PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR CONTRIBUTORS
+.\" BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+.\" CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+.\" SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+.\" INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+.\" CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+.\" ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+.\" POSSIBILITY OF SUCH DAMAGE.
+.\"
+.Dd July 29, 2011
+.Dt ILOGB 3
+.Os
+.Sh NAME
+.Nm ilogb ,
+.Nm ilogbf ,
+.Nm ilogbl
+.Nd an unbiased exponent
+.Sh LIBRARY
+.Lb libm
+.Sh SYNOPSIS
+.In math.h
+.Ft int
+.Fn ilogb "double x"
+.Ft int
+.Fn ilogbf "float x"
+.Ft int
+.Fn ilogbl "long double x"
+.Sh DESCRIPTION
+The
+.Fn ilogb ,
+.Fn ilogbf ,
+and
+.Fn ilogbl
+functions return the exponent of the non-zero real floating-point number
+.Fa x
+as a signed integer value.
+Formally the return value is the integral part of
+.Bd -ragged -offset indent
+log_r |
+.Va x | ,
+.Ed
+.Pp
+where
+.Fa r
+is the radix of the machine's floating-point arithmetic defined by the
+.Dv FLT_RADIX
+constant in
+.In float.h .
+.Sh RETURN VALUES
+As described above, upon successful completion,
+the functions return the exponent.
+Functionally this is the same as calling the corresponding
+.Xr logb 3
+function and casting the return value to
+.Vt int .
+.Pp
+The following special cases may occur:
+.Bl -enum -offset indent
+.It
+If
+.Fa x
+is zero, the value of
+.Dv FP_ILOGB0
+is returned and a domain error occurs.
+.It
+If
+.Fa x
+is infinite, a domain error occurs and the value of
+.Dv INT_MAX
+is returned.
+.It
+If
+.Fa x
+is \*(Na, a domain error is raised and the value of
+.Dv FP_ILOGBNAN
+is returned.
+.It
+If the correct value is outside the range of the return type,
+a domain error occurs but an unspecified value is returned.
+.El
+.Sh SEE ALSO
+.Xr ilog2 3 ,
+.Xr logb 3 ,
+.Xr math 3
+.Sh STANDARDS
+The described functions conform to
+.St -isoC-99 .
+.Sh BUGS
+Neither
+.Dv FP_ILOGB0
+nor
+.Dv FP_ILOGBNAN
+is defined currently in
+.Nx .
--- /dev/null
+.\" $NetBSD: log.3,v 1.3 2011/09/13 08:51:32 wiz Exp $
+.\"
+.\" Copyright (c) 2011 Jukka Ruohonen <jruohonen@iki.fi>
+.\" All rights reserved.
+.\"
+.\" Redistribution and use in source and binary forms, with or without
+.\" modification, are permitted provided that the following conditions
+.\" are met:
+.\" 1. Redistributions of source code must retain the above copyright
+.\" notice, this list of conditions and the following disclaimer.
+.\" 2. Redistributions in binary form must reproduce the above copyright
+.\" notice, this list of conditions and the following disclaimer in the
+.\" documentation and/or other materials provided with the distribution.
+.\"
+.\" THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
+.\" ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+.\" TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+.\" PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR CONTRIBUTORS
+.\" BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+.\" CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+.\" SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+.\" INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+.\" CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+.\" ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+.\" POSSIBILITY OF SUCH DAMAGE.
+.\"
+.Dd September 13, 2011
+.Dt LOG 3
+.Os
+.Sh NAME
+.Nm log ,
+.Nm logf ,
+.Nm log10 ,
+.Nm log10f ,
+.Nm log1p ,
+.Nm log1pf
+.Nm log2 ,
+.Nm log2f ,
+.Nd logarithm functions
+.Sh LIBRARY
+.Lb libm
+.Sh SYNOPSIS
+.In math.h
+.Ft double
+.Fn log "double x"
+.Ft float
+.Fn logf "float x"
+.Ft double
+.Fn log10 "double x"
+.Ft float
+.Fn log10f "float x"
+.Ft double
+.Fn log1p "double x"
+.Ft float
+.Fn log1pf "float x"
+.Ft double
+.Fn log2 "double x"
+.Ft float
+.Fn log2f "float x"
+.Sh DESCRIPTION
+The following functions compute logarithms:
+.Bl -bullet -offset 2n
+.It
+The
+.Fn log
+and
+.Fn logf
+functions return the natural logarithm.
+.It
+The
+.Fn log10
+and
+.Fn log10f
+functions return the base 10 logarithm.
+.It
+The
+.Fn log1p
+and
+.Fn log1pf
+functions return the natural logarithm of (1.0 +
+.Fa x )
+accurately even for very small values of
+.Fa x .
+.It
+The
+.Fn log2
+and
+.Fn log2f
+functions return the base 2 logarithm.
+.El
+.Sh RETURN VALUES
+Upon successful completion, the functions return the logarithm of
+.Fa x
+as descibed above.
+Otherwise the following may occur:
+.Bl -enum -offset indent
+.It
+If
+.Fa x
+is \*(Na, all functions return \*(Na.
+.It
+If
+.Fa x
+is positive infinity, all functions return
+.Fa x .
+If
+.Fa x
+is negative infinity, all functions return \*(Na.
+.It
+If
+.Fa x
+is +0.0 or -0.0, the
+.Fn log ,
+.Fn log10 ,
+and
+.Fn log2
+families return either
+.Dv -HUGE_VAL ,
+.Dv -HUGE_VALF ,
+or
+.Dv -HUGE_VALL ,
+whereas the
+.Fn log1p
+family returns
+.Fa x .
+.It
+If
+.Fa x
+is +1.0, the
+.Fn log ,
+.Fn log10 ,
+and
+.Fn log2
+families return +0.0.
+If
+.Fa x
+is -1.0, the
+.Fn log1p
+family returns
+.Dv -HUGE_VAL ,
+.Dv -HUGE_VALF ,
+or
+.Dv -HUGE_VALL .
+.El
+.Pp
+In addition, on a
+.Tn VAX ,
+.Va errno
+is set to
+.Er EDOM
+and the reserved operand is returned
+by
+.Fn log
+unless
+.Fa x
+\*[Gt] 0, by
+.Fn log1p
+unless
+.Fa x
+\*[Gt] \-1.
+.Sh SEE ALSO
+.Xr exp 3 ,
+.Xr ilogb 3 ,
+.Xr math 3
+.Sh STANDARDS
+The described functions conform to
+.St -isoC-99 .
+.Sh HISTORY
+The history of the logarithm functions dates back to
+.At v6 .
+.\" $NetBSD: math.3,v 1.25 2011/09/22 18:14:09 njoly Exp $
+.\"
.\" Copyright (c) 1985 Regents of the University of California.
.\" All rights reserved.
.\"
.\" SUCH DAMAGE.
.\"
.\" from: @(#)math.3 6.10 (Berkeley) 5/6/91
-.\" $NetBSD: math.3,v 1.18 2003/12/03 23:31:21 jschauma Exp $
.\"
-.TH MATH 3 "July 12, 2009"
-.UC 4
-.ds up \fIulp\fR
-.ds nn \fINaN\fR
-.de If
-.if n \\
-\\$1Infinity\\$2
-.if t \\
-\\$1\\(if\\$2
-..
-.SH NAME
-math \- introduction to mathematical library functions
-.SH DESCRIPTION
-These functions constitute the C math library,
-.I libm.
-The link editor searches this library under the \*(lq\-lm\*(rq option.
+.Dd February 23, 2007
+.Dt MATH 3
+.Os
+.Sh NAME
+.Nm math
+.Nd introduction to mathematical library functions
+.Sh LIBRARY
+.Lb libm
+.Sh SYNOPSIS
+.In math.h
+.Sh DESCRIPTION
+These functions constitute the C
+.Lb libm .
Declarations for these functions may be obtained from the include file
-.RI \*[Lt] math.h \*[Gt].
+.In math.h .
.\" The Fortran math library is described in ``man 3f intro''.
-.SH "LIST OF FUNCTIONS"
-.sp 2
-.nf
-.ta \w'copysign'u+2n +\w'lgamma.3'u+10n +\w'inverse trigonometric func'u
-\fIName\fP \fIAppears on Page\fP \fIDescription\fP \fIError Bound (ULPs)\fP
-.ta \w'copysign'u+4n +\w'lgamma.3'u+4n +\w'inverse trigonometric function'u+6nC
-.sp 5p
-acos acos.3 inverse trigonometric function 3
-acosh acosh.3 inverse hyperbolic function 3
-asin asin.3 inverse trigonometric function 3
-asinh asinh.3 inverse hyperbolic function 3
-atan atan.3 inverse trigonometric function 1
-atanh atanh.3 inverse hyperbolic function 3
-atan2 atan2.3 inverse trigonometric function 2
-cabs hypot.3 complex absolute value 1
-cbrt sqrt.3 cube root 1
-ceil ceil.3 integer no less than 0
-copysign ieee.3 copy sign bit 0
-cos cos.3 trigonometric function 1
-cosh cosh.3 hyperbolic function 3
-erf erf.3 error function ???
-erfc erf.3 complementary error function ???
-exp exp.3 exponential 1
-expm1 exp.3 exp(x)\-1 1
-fabs fabs.3 absolute value 0
-fdim fdim.3 positive difference ???
-finite ieee.3 test for finity 0
-floor floor.3 integer no greater than 0
-fmax fmax.3 maximum function ???
-fmin fmin.3 minimum function ???
-fmod fmod.3 remainder ???
-hypot hypot.3 Euclidean distance 1
-ilogb ieee.3 exponent extraction 0
-isinf isinf.3 test for infinity 0
-isnan isnan.3 test for not-a-number 0
-j0 j0.3 Bessel function ???
-j1 j0.3 Bessel function ???
-jn j0.3 Bessel function ???
-lgamma lgamma.3 log gamma function ???
-log exp.3 natural logarithm 1
-log10 exp.3 logarithm to base 10 3
-log1p exp.3 log(1+x) 1
-nan nan.3 return quiet \*(nn 0
-nextafter ieee.3 next representable number 0
-pow exp.3 exponential x**y 60\-500
-remainder ieee.3 remainder 0
-rint rint.3 round to nearest integer 0
-scalbn ieee.3 exponent adjustment 0
-sin sin.3 trigonometric function 1
-sinh sinh.3 hyperbolic function 3
-sqrt sqrt.3 square root 1
-tan tan.3 trigonometric function 3
-tanh tanh.3 hyperbolic function 3
-trunc trunc.3 nearest integral value 3
-y0 j0.3 Bessel function ???
-y1 j0.3 Bessel function ???
-yn j0.3 Bessel function ???
-.ta
-.fi
-.SH "LIST OF DEFINED VALUES"
-.sp 2
-.nf
-.ta \w'M_2_SQRTPI'u+2n +\w'1.12837916709551257390'u+4n +\w'2/sqrt(pi)'u+6nC
-\fIName\fP \fIValue\fP \fIDescription\fP
-.ta \w'M_2_SQRTPI'u+2n +\w'1.12837916709551257390'u+4n +\w'2/sqrt(pi)'u+6nC
-.sp 3p
-M_E 2.7182818284590452354 e
-M_LOG2E 1.4426950408889634074 log 2e
-M_LOG10E 0.43429448190325182765 log 10e
-M_LN2 0.69314718055994530942 log e2
-M_LN10 2.30258509299404568402 log e10
-M_PI 3.14159265358979323846 pi
-M_PI_2 1.57079632679489661923 pi/2
-M_PI_4 0.78539816339744830962 pi/4
-M_1_PI 0.31830988618379067154 1/pi
-M_2_PI 0.63661977236758134308 2/pi
-M_2_SQRTPI 1.12837916709551257390 2/sqrt(pi)
-M_SQRT2 1.41421356237309504880 sqrt(2)
-M_SQRT1_2 0.70710678118654752440 1/sqrt(2)
-.ta
-.fi
-.SH NOTES
+.Ss List of Functions
+.Bl -column "copysignX" "gammaX3XX" "inverse trigonometric funcX"
+.It Sy Name Ta Sy Man page Ta Sy Description Ta Sy Error Bound Dv ( ULP Ns No s)
+.It acos Ta Xr acos 3 Ta inverse trigonometric function Ta 3
+.It acosh Ta Xr acosh 3 Ta inverse hyperbolic function Ta 3
+.It asin Ta Xr asin 3 Ta inverse trigonometric function Ta 3
+.It asinh Ta Xr asinh 3 Ta inverse hyperbolic function Ta 3
+.It atan Ta Xr atan 3 Ta inverse trigonometric function Ta 1
+.It atanh Ta Xr atanh 3 Ta inverse hyperbolic function Ta 3
+.It atan2 Ta Xr atan2 3 Ta inverse trigonometric function Ta 2
+.It cbrt Ta Xr sqrt 3 Ta cube root Ta 1
+.It ceil Ta Xr ceil 3 Ta integer no less than Ta 0
+.It copysign Ta Xr copysign 3 Ta copy sign bit Ta 0
+.It cos Ta Xr cos 3 Ta trigonometric function Ta 1
+.It cosh Ta Xr cosh 3 Ta hyperbolic function Ta 3
+.It erf Ta Xr erf 3 Ta error function Ta ???
+.It erfc Ta Xr erf 3 Ta complementary error function Ta ???
+.It exp Ta Xr exp 3 Ta exponential Ta 1
+.It expm1 Ta Xr exp 3 Ta exp(x)\-1 Ta 1
+.It fabs Ta Xr fabs 3 Ta absolute value Ta 0
+.It finite Ta Xr finite 3 Ta test for finity Ta 0
+.It floor Ta Xr floor 3 Ta integer no greater than Ta 0
+.It fmod Ta Xr fmod 3 Ta remainder Ta ???
+.It hypot Ta Xr hypot 3 Ta Euclidean distance Ta 1
+.It ilogb Ta Xr ilogb 3 Ta exponent extraction Ta 0
+.It isinf Ta Xr isinf 3 Ta test for infinity Ta 0
+.It isnan Ta Xr isnan 3 Ta test for not-a-number Ta 0
+.It j0 Ta Xr j0 3 Ta Bessel function Ta ???
+.It j1 Ta Xr j0 3 Ta Bessel function Ta ???
+.It jn Ta Xr j0 3 Ta Bessel function Ta ???
+.It lgamma Ta Xr lgamma 3 Ta log gamma function Ta ???
+.It log Ta Xr log 3 Ta natural logarithm Ta 1
+.It log10 Ta Xr log 3 Ta logarithm to base 10 Ta 3
+.It log1p Ta Xr log 3 Ta log(1+x) Ta 1
+.It nan Ta Xr nan 3 Ta return quiet \*(Na Ta 0
+.It nextafter Ta Xr nextafter 3 Ta next representable number Ta 0
+.It pow Ta Xr pow 3 Ta exponential x**y Ta 60\-500
+.It remainder Ta Xr remainder 3 Ta remainder Ta 0
+.It rint Ta Xr rint 3 Ta round to nearest integer Ta 0
+.It scalbn Ta Xr scalbn 3 Ta exponent adjustment Ta 0
+.It sin Ta Xr sin 3 Ta trigonometric function Ta 1
+.It sinh Ta Xr sinh 3 Ta hyperbolic function Ta 3
+.It sqrt Ta Xr sqrt 3 Ta square root Ta 1
+.It tan Ta Xr tan 3 Ta trigonometric function Ta 3
+.It tanh Ta Xr tanh 3 Ta hyperbolic function Ta 3
+.It trunc Ta Xr trunc 3 Ta nearest integral value Ta 3
+.It y0 Ta Xr j0 3 Ta Bessel function Ta ???
+.It y1 Ta Xr j0 3 Ta Bessel function Ta ???
+.It yn Ta Xr j0 3 Ta Bessel function Ta ???
+.El
+.Ss List of Defined Values
+.Bl -column "M_2_SQRTPIXX" "1.12837916709551257390XX" "2/sqrt(pi)XXX"
+.It Sy Name Ta Sy Value Ta Sy Description
+.It M_E 2.7182818284590452354 e
+.It M_LOG2E 1.4426950408889634074 log 2e
+.It M_LOG10E 0.43429448190325182765 log 10e
+.It M_LN2 0.69314718055994530942 log e2
+.It M_LN10 2.30258509299404568402 log e10
+.It M_PI 3.14159265358979323846 pi
+.It M_PI_2 1.57079632679489661923 pi/2
+.It M_PI_4 0.78539816339744830962 pi/4
+.It M_1_PI 0.31830988618379067154 1/pi
+.It M_2_PI 0.63661977236758134308 2/pi
+.It M_2_SQRTPI 1.12837916709551257390 2/sqrt(pi)
+.It M_SQRT2 1.41421356237309504880 sqrt(2)
+.It M_SQRT1_2 0.70710678118654752440 1/sqrt(2)
+.El
+.Sh NOTES
In 4.3 BSD, distributed from the University of California
in late 1985, most of the foregoing functions come in two
versions, one for the double\-precision "D" format in the
similarly, as should be expected from programs more accurate
and robust than was the norm when UNIX was born.
For instance, the programs are accurate to within the numbers
-of \*(ups tabulated above; an \*(up is one \fIU\fRnit in the \fIL\fRast
-\fIP\fRlace.
+of
+.Dv ULPs
+tabulated above; an
+.Dv ULP
+is one Unit in the Last Place.
And the programs have been cured of anomalies that
-afflicted the older math library \fIlibm\fR in which incidents like
+afflicted the older math library
+in which incidents like
the following had been reported:
-.RS
+.Bd -literal -offset indent
sqrt(\-1.0) = 0.0 and log(\-1.0) = \-1.7e38.
-.br
cos(1.0e\-11) \*[Gt] cos(0.0) \*[Gt] 1.0.
-.br
-pow(x,1.0)
-.if n \
-!=
-.if t \
-\(!=
-x when x = 2.0, 3.0, 4.0, ..., 9.0.
-.br
+pow(x,1.0) \(!= x when x = 2.0, 3.0, 4.0, ..., 9.0.
pow(\-1.0,1.0e10) trapped on Integer Overflow.
-.br
sqrt(1.0e30) and sqrt(1.0e\-30) were very slow.
-.RE
+.Ed
However the two versions do differ in ways that have to be
explained, to which end the following notes are provided.
-.PP
-\fBDEC VAX\-11 D_floating\-point:\fR
-.PP
-This is the format for which the original math library \fIlibm\fR
+.Ss DEC VAX\-11 D_floating\-point
+This is the format for which the original math library
was developed, and to which this manual is still principally dedicated.
-It is \fIthe\fR double\-precision format for the PDP\-11
+It is
+.Em the
+double\-precision format for the PDP\-11
and the earlier VAX\-11 machines; VAX\-11s after 1983 were
provided with an optional "G" format closer to the IEEE
double\-precision format.
The earlier DEC MicroVAXs have no D format, only G double\-precision.
(Why?
Why not?)
-.PP
+.Pp
Properties of D_floating\-point:
-.RS
-Wordsize: 64 bits, 8 bytes.
-Radix: Binary.
-.br
-Precision: 56
-.if n \
-sig.
-.if t \
-significant
-bits, roughly like 17
-.if n \
-sig.
-.if t \
-significant
-decimals.
-.RS
+.Bl -hang -offset indent
+.It Wordsize :
+64 bits, 8 bytes.
+.It Radix :
+Binary.
+.It Precision :
+56 significant bits, roughly like 17 significant decimals.
If x and x' are consecutive positive D_floating\-point
-numbers (they differ by 1 \*(up), then
-.br
-1.3e\-17 \*[Lt] 0.5**56 \*[Lt] (x'\-x)/x \*[Le] 0.5**55 \*[Lt] 2.8e\-17.
-.RE
-.nf
-.ta \w'Range:'u+1n +\w'Underflow threshold'u+1n +\w'= 2.0**127'u+1n
-Range: Overflow threshold = 2.0**127 = 1.7e38.
- Underflow threshold = 0.5**128 = 2.9e\-39.
- NOTE: THIS RANGE IS COMPARATIVELY NARROW.
-.ta
-.fi
-.RS
+numbers (they differ by 1
+.Dv ULP ) ,
+then
+.Dl 1.3e\-17 \*[Lt] 0.5**56 \*[Lt] (x'\-x)/x \*[Le] 0.5**55 \*[Lt] 2.8e\-17.
+.It Range :
+.Bl -column "Underflow thresholdX" "2.0**127X"
+.It Overflow threshold = 2.0**127 = 1.7e38.
+.It Underflow threshold = 0.5**128 = 2.9e\-39.
+.El
+.Em NOTE: THIS RANGE IS COMPARATIVELY NARROW.
+.Pp
Overflow customarily stops computation.
-.br
Underflow is customarily flushed quietly to zero.
-.br
-CAUTION:
-.RS
+.Em CAUTION :
It is possible to have x
-.if n \
-!=
-.if t \
\(!=
-y and yet
-x\-y = 0 because of underflow.
+y and yet x\-y = 0 because of underflow.
Similarly x \*[Gt] y \*[Gt] 0 cannot prevent either x\(**y = 0
-or y/x = 0 from happening without warning.
-.RE
-.RE
-Zero is represented ambiguously.
-.RS
+or y/x = 0 from happening without warning.
+.It Zero is represented ambiguously :
Although 2**55 different representations of zero are accepted by
the hardware, only the obvious representation is ever produced.
There is no \-0 on a VAX.
-.RE
-.If
-is not part of the VAX architecture.
-.br
-Reserved operands:
-.RS
+.It \*(If is not part of the VAX architecture .
+.It Reserved operands :
of the 2**55 that the hardware
recognizes, only one of them is ever produced.
Any floating\-point operation upon a reserved
operand, even a MOVF or MOVD, customarily stops
computation, so they are not much used.
-.RE
-Exceptions:
-.RS
+.It Exceptions :
Divisions by zero and operations that
overflow are invalid operations that customarily
stop computation or, in earlier machines, produce
reserved operands that will stop computation.
-.RE
-Rounding:
-.RS
+.It Rounding :
Every rational operation (+, \-, \(**, /) on a
VAX (but not necessarily on a PDP\-11), if not an
over/underflow nor division by zero, is rounded to
-within half an \*(up, and when the rounding error is
-exactly half an \*(up then rounding is away from 0.
-.RE
-.RE
-.PP
+within half an
+.Dv ULP ,
+and when the rounding error is
+exactly half an
+.Dv ULP
+then rounding is away from 0.
+.El
+.Pp
Except for its narrow range, D_floating\-point is one of the
better computer arithmetics designed in the 1960's.
Its properties are reflected fairly faithfully in the elementary
like 1/0; and sqrt(\-3) and acos(3) behave like 0/0;
they all produce reserved operands and/or stop computation!
The situation is described in more detail in manual pages.
-.RS
-.ll -0.5i
-\fIThis response seems excessively punitive, so it is destined
-to be replaced at some time in the foreseeable future by a
-more flexible but still uniform scheme being developed to
-handle all floating\-point arithmetic exceptions neatly.\fR
-.ll +0.5i
-.RE
-.PP
-How do the functions in 4.3 BSD's new \fIlibm\fR for UNIX
+.Pp
+.Em This response seems excessively punitive, so it is destined
+.Em to be replaced at some time in the foreseeable future by a
+.Em more flexible but still uniform scheme being developed to
+.Em handle all floating\-point arithmetic exceptions neatly.
+.Pp
+How do the functions in 4.3 BSD's new math library for UNIX
compare with their counterparts in DEC's VAX/VMS library?
Some of the VMS functions are a little faster, some are
a little more accurate, some are more puritanical about
exceptions (like pow(0.0,0.0) and atan2(0.0,0.0)),
and most occupy much more memory than their counterparts in
-\fIlibm\fR.
+libm.
The VMS codes interpolate in large table to achieve
-speed and accuracy; the \fIlibm\fR codes use tricky formulas
+speed and accuracy; the libm codes use tricky formulas
compact enough that all of them may some day fit into a ROM.
-.PP
+.Pp
More important, DEC regards the VMS codes as proprietary
and guards them zealously against unauthorized use.
-But the \fIlibm\fR codes in 4.3 BSD are intended for the public domain;
+But the libm codes in 4.3 BSD 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
researches by reporting experience with the codes.
Therefore no user of UNIX on a machine whose arithmetic resembles
-VAX D_floating\-point need use anything worse than the new \fIlibm\fR.
-.PP
-\fBIEEE STANDARD 754 Floating\-Point Arithmetic:\fR
-.PP
+VAX D_floating\-point need use anything worse than the new libm.
+.Ss IEEE STANDARD 754 Floating\-Point Arithmetic
This standard is on its way to becoming more widely adopted
than any other design for computer arithmetic.
VLSI chips that conform to some version of that standard have been
produced by a host of manufacturers, among them ...
-.nf
-.ta 0.5i +\w'Intel i8070, i80287'u+6n
- Intel i8087, i80287 National Semiconductor 32081
- Motorola 68881 Weitek WTL-1032, ... , -1165
- Zilog Z8070 Western Electric (AT\*[Am]T) WE32106.
-.ta
-.fi
+.Bl -column "Intel i8070, i80287XX"
+.It Intel i8087, i80287 National Semiconductor 32081
+.It 68881 Weitek WTL-1032, ... , -1165
+.It Zilog Z8070 Western Electric (AT\*[Am]T) WE32106.
+.El
Other implementations range from software, done thoroughly
in the Apple Macintosh, through VLSI in the Hewlett\-Packard
9000 series, to the ELXSI 6400 running ECL at 3 Megaflops.
IEEE versions of most of the elementary functions listed
above could easily be converted to run on a MicroVAX, though
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 libm for machines that conform to
+IEEE 754 are intended primarily for the National Semiconductor 32081
and WTL 1164/65.
To use these codes with the Intel or Zilog
chips, or with the Apple Macintosh or ELXSI 6400, is to
forego the use of better codes provided (perhaps freely) by
those companies and designed by some of the authors of the
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,
+Except for
+.Fn atan ,
+.Fn cbrt ,
+.Fn erf ,
+.Fn erfc ,
+.Fn hypot ,
+.Fn j0-jn ,
+.Fn lgamma ,
+.Fn pow ,
+and
+.Fn y0\-yn ,
+the Motorola 68881 has all the functions in libm on chip,
and faster and more accurate;
-it, Apple, the i8087, Z8070 and WE32106 all use 64
-.if n \
-sig.
-.if t \
-significant
-bits.
+it, Apple, the i8087, Z8070 and WE32106 all use 64 significant bits.
The main virtue of 4.3 BSD's
-\fIlibm\fR codes is that they are intended for the public domain;
+libm 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
researches by reporting experience with the codes.
Therefore no user of UNIX on a machine that conforms to
-IEEE 754 need use anything worse than the new \fIlibm\fR.
-.PP
+IEEE 754 need use anything worse than the new libm.
+.Pp
Properties of IEEE 754 Double\-Precision:
-.RS
-Wordsize: 64 bits, 8 bytes.
-Radix: Binary.
-.br
-Precision: 53
-.if n \
-sig.
-.if t \
-significant
-bits, roughly like 16
-.if n \
-sig.
-.if t \
-significant
-decimals.
-.RS
+.Bl -hang -offset indent
+.It Wordsize :
+64 bits, 8 bytes.
+.It Radix :
+Binary.
+.It Precision :
+53 significant bits, roughly like 16 significant decimals.
If x and x' are consecutive positive Double\-Precision
-numbers (they differ by 1 \*(up), then
-.br
-1.1e\-16 \*[Lt] 0.5**53 \*[Lt] (x'\-x)/x \*[Le] 0.5**52 \*[Lt] 2.3e\-16.
-.RE
-.nf
-.ta \w'Range:'u+1n +\w'Underflow threshold'u+1n +\w'= 2.0**1024'u+1n
-Range: Overflow threshold = 2.0**1024 = 1.8e308
- Underflow threshold = 0.5**1022 = 2.2e\-308
-.ta
-.fi
-.RS
-Overflow goes by default to a signed
-.If "" .
-.br
-Underflow is \fIGradual,\fR rounding to the nearest
+numbers (they differ by 1
+.Dv ULP ) ,
+then
+.Dl 1.1e\-16 \*[Lt] 0.5**53 \*[Lt] (x'\-x)/x \*[Le] 0.5**52 \*[Lt] 2.3e\-16.
+.It Range :
+.Bl -column "Underflow thresholdX" "2.0**1024X"
+.It Overflow threshold = 2.0**1024 = 1.8e308
+.It Underflow threshold = 0.5**1022 = 2.2e\-308
+.El
+Overflow goes by default to a signed \*(If.
+Underflow is
+.Sy Gradual ,
+rounding to the nearest
integer multiple of 0.5**1074 = 4.9e\-324.
-.RE
-Zero is represented ambiguously as +0 or \-0.
-.RS
+.It Zero is represented ambiguously as +0 or \-0:
Its sign transforms correctly through multiplication or
division, and is preserved by addition of zeros
with like signs; but x\-x yields +0 for every
sign are division by zero and copysign(x,\(+-0).
In particular, comparison (x \*[Gt] y, x \*[Ge] y, etc.)
cannot be affected by the sign of zero; but if
-finite x = y then
-.If
+finite x = y then \*(If
\&= 1/(x\-y)
-.if n \
-!=
-.if t \
\(!=
\-1/(y\-x) =
-.If \- .
-.RE
-.If
-is signed.
-.RS
+\- \*(If .
+.It \*(If is signed :
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
+\*(If (finite)/\(+- \0=\0\(+-0
(nonzero)/0 =
-.If \(+- .
+\(+- \*(If.
But
-.if n \
-Infinity\-Infinity, Infinity\(**0 and Infinity/Infinity
-.if t \
\(if\-\(if, \(if\(**0 and \(if/\(if
are, like 0/0 and sqrt(\-3),
-invalid operations that produce \*(nn. ...
-.RE
-Reserved operands:
-.RS
+invalid operations that produce \*(Na.
+.It Reserved operands :
there are 2**53\-2 of them, all
-called \*(nn (\fIN\fRot \fIa N\fRumber).
-Some, called Signaling \*(nns, trap any floating\-point operation
+called \*(Na (Not A Number).
+Some, called Signaling \*[Na]s, trap any floating\-point operation
performed upon them; they are used to mark missing
or uninitialized values, or nonexistent elements of arrays.
-The rest are Quiet \*(nns; they are
+The rest are Quiet \*[Na]s; they are
the default results of Invalid Operations, and
propagate through subsequent arithmetic operations.
If x
-.if n \
-!=
-.if t \
\(!=
-x then x is \*(nn; every other predicate
-(x \*[Gt] y, x = y, x \*[Lt] y, ...) is FALSE if \*(nn is involved.
-.br
-NOTE: Trichotomy is violated by \*(nn.
-.RS
+x then x is \*(Na; every other predicate
+(x \*[Gt] y, x = y, x \*[Lt] y, ...) is FALSE if \*(Na is involved.
+.Pp
+.Em NOTE :
+Trichotomy is violated by \*(Na.
Besides being FALSE, predicates that entail ordered
comparison, rather than mere (in)equality,
-signal Invalid Operation when \*(nn is involved.
-.RE
-.RE
-Rounding:
-.RS
+signal Invalid Operation when \*(Na is involved.
+.It Rounding :
Every algebraic operation (+, \-, \(**, /,
-.if n \
-sqrt)
-.if t \
\(sr)
-is rounded by default to within half an \*(up, and
-when the rounding error is exactly half an \*(up then
-the rounded value's least significant bit is zero.
+is rounded by default to within half an
+.Dv ULP ,
+and when the rounding error is exactly half an
+.Dv ULP
+then the rounded value's least significant bit is zero.
This kind of rounding is usually the best kind,
sometimes provably so; for instance, for every
x = 1.0, 2.0, 3.0, 4.0, ..., 2.0**52, we find
But no single kind of rounding can be
proved best for every circumstance, so IEEE 754
provides rounding towards zero or towards
-.If +
++\*(If
or towards
-.If \-
+\-\*(If
at the programmer's option.
And the same kinds of rounding are specified for
Binary\-Decimal Conversions, at least for magnitudes
between roughly 1.0e\-10 and 1.0e37.
-.RE
-Exceptions:
-.RS
+.It Exceptions :
IEEE 754 recognizes five kinds of floating\-point exceptions,
listed below in declining order of probable importance.
-.RS
-.nf
-.ta \w'Invalid Operation'u+6n +\w'Gradual Underflow'u+2n
-Exception Default Result
-.tc \(ru
-
-.tc
-Invalid Operation \*(nn, or FALSE
-.if n \{\
-Overflow \(+-Infinity
-Divide by Zero \(+-Infinity \}
-.if t \{\
-Overflow \(+-\(if
-Divide by Zero \(+-\(if \}
-Underflow Gradual Underflow
-Inexact Rounded value
-.ta
-.fi
-.RE
-NOTE: An Exception is not an Error unless handled badly.
+.Bl -column "Invalid OperationX" "Gradual OverflowX"
+.It Sy Exception Ta Sy Default Result
+.It Invalid Operation \*(Na, or FALSE
+.It Overflow \(+-\(if
+.It Divide by Zero \(+-\(if \}
+.It Underflow Gradual Underflow
+.It Inexact Rounded value
+.El
+.Pp
+.Em NOTE :
+An Exception is not an Error unless handled badly.
What makes a class of exceptions exceptional
is that no single default response can be satisfactory
in every instance.
response will serve most instances satisfactorily,
the unsatisfactory instances cannot justify aborting
computation every time the exception occurs.
-.RE
-.PP
+.El
+.Pp
For each kind of floating\-point exception, IEEE 754
provides a Flag that is raised each time its exception
is signaled, and stays raised until the program resets it.
Thus, IEEE 754 provides three ways by which programs
may cope with exceptions for which the default result
might be unsatisfactory:
-.IP 1) \w'\0\0\0\0'u
+.Bl -enum
+.It
Test for a condition that might cause an exception
later, and branch to avoid the exception.
-.IP 2) \w'\0\0\0\0'u
+.It
Test a flag to see whether an exception has occurred
since the program last reset its flag.
-.IP 3) \w'\0\0\0\0'u
+.It
Test a result to see whether it is a value that only
an exception could have produced.
-.RS
-CAUTION: The only reliable ways to discover
+.Em CAUTION :
+The only reliable ways to discover
whether Underflow has occurred are to test whether
products or quotients lie closer to zero than the
underflow threshold, or to test the Underflow flag.
(Sums and differences cannot underflow in
IEEE 754; if x
-.if n \
-!=
-.if t \
\(!=
y then x\-y is correct to
full precision and certainly nonzero regardless of
underflow threshold, as is almost always the case,
digits lost to gradual underflow will not be missed
because they would have been rounded off anyway.
-So gradual underflows are usually \fIprovably\fR ignorable.
+So gradual underflows are usually
+.Em provably
+ignorable.
The same cannot be said of underflows flushed to 0.
-.RE
-.PP
+.Pp
At the option of an implementor conforming to IEEE 754,
other ways to cope with exceptions may be provided:
-.IP 4) \w'\0\0\0\0'u
+.It
ABORT.
This mechanism classifies an exception in
advance as an incident to be handled by means
statements like "ON ERROR GO TO ...".
Different languages offer different forms of this statement,
but most share the following characteristics:
-.IP \(em \w'\0\0\0\0'u
+.Bl -dash
+.It
No means is provided to substitute a value for
the offending operation's result and resume
computation from what may be the middle of an expression.
An exceptional result is abandoned.
-.IP \(em \w'\0\0\0\0'u
+.It
In a subprogram that lacks an error\-handling
statement, an exception causes the subprogram to
abort within whatever program called it, and so
on back up the chain of calling subprograms until
an error\-handling statement is encountered or the
whole task is aborted and memory is dumped.
-.IP 5) \w'\0\0\0\0'u
+.El
+.It
STOP.
This mechanism, requiring an interactive
debugging environment, is more for the programmer
unexceptionable, so the programmer ought ideally
to be able to resume execution after each one as if
execution had not been stopped.
-.IP 6) \w'\0\0\0\0'u
+.It
\&... Other ways lie beyond the scope of this document.
-.RE
-.PP
+.El
+.Pp
The crucial problem for exception handling is the problem of
Scope, and the problem's solution is understood, but not
enough manpower was available to implement it fully in time
-to be distributed in 4.3 BSD's \fIlibm\fR.
+to be distributed in 4.3 BSD's libm.
Ideally, each elementary function should act
as if it were indivisible, or atomic, in the sense that ...
-.IP i) \w'iii)'u+2n
+.Bl -enum
+.It
No exception should be signaled that is not deserved by
the data supplied to that function.
-.IP ii) \w'iii)'u+2n
+.It
Any exception signaled should be identified with that
function rather than with one of its subroutines.
-.IP iii) \w'iii)'u+2n
+.It
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.
-.PP
-Ideally, every programmer should be able \fIconveniently\fR to
-turn a debugged subprogram into one that appears atomic to
+.El
+.Pp
+Ideally, every programmer should be able
+.Em conveniently
+to turn a debugged subprogram into one that appears atomic to
its users.
But simulating all three characteristics of an
atomic function is still a tedious affair, entailing hosts
of tests and saves\-restores; work is under way to ameliorate
the inconvenience.
-.PP
-Meanwhile, the functions in \fIlibm\fR are only approximately atomic.
+.Pp
+Meanwhile, the functions in libm are only approximately atomic.
They signal no inappropriate exception except possibly ...
-.RS
-Over/Underflow
-.RS
+.Bl -ohang -offset indent
+.It Over/Underflow
when a result, if properly computed, might have lain barely within range, and
-.RE
-Inexact in \fIcabs\fR, \fIcbrt\fR, \fIhypot\fR, \fIlog10\fR and \fIpow\fR
-.RS
+.It Inexact in Fn cbrt , Fn hypot , Fn log10 and Fn pow
when it happens to be exact, thanks to fortuitous cancellation of errors.
-.RE
-.RE
+.El
Otherwise, ...
-.RS
-Invalid Operation is signaled only when
-.RS
-any result but \*(nn would probably be misleading.
-.RE
-Overflow is signaled only when
-.RS
+.Bl -ohang -offset indent
+.It Invalid Operation is signaled only when
+any result but \*(Na would probably be misleading.
+.It Overflow is signaled only when
the exact result would be finite but beyond the overflow threshold.
-.RE
-Divide\-by\-Zero is signaled only when
-.RS
+.It Divide\-by\-Zero is signaled only when
a function takes exactly infinite values at finite operands.
-.RE
-Underflow is signaled only when
-.RS
+.It Underflow is signaled only when
the exact result would be nonzero but tinier than the underflow threshold.
-.RE
-Inexact is signaled only when
-.RS
+.It Inexact is signaled only when
greater range or precision would be needed to represent the exact result.
-.RE
-.RE
+.El
.\" .Sh FILES
.\" .Bl -tag -width /usr/lib/profile/libm.a -compact
.\" .It Pa /usr/lib/libm.a
.\" .It Pa /usr/lib/profile/libm.a
.\" the static math library compiled for profiling
.\" .El
-.SH SEE ALSO
+.Sh SEE ALSO
An explanation of IEEE 754 and its proposed extension p854
was published in the IEEE magazine MICRO in August 1984 under
the title "A Proposed Radix\- and Word\-length\-independent
and in the ACM SIGNUM Newsletter Special Issue of
Oct. 1979, may be helpful although they pertain to
superseded drafts of the standard.
-.SH BUGS
+.Sh BUGS
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
--- /dev/null
+.\" $NetBSD: nextafter.3,v 1.4 2011/09/18 05:33:14 jruoho Exp $
+.\"
+.\" Copyright (c) 2011 Jukka Ruohonen <jruohonen@iki.fi>
+.\" All rights reserved.
+.\"
+.\" Redistribution and use in source and binary forms, with or without
+.\" modification, are permitted provided that the following conditions
+.\" are met:
+.\" 1. Redistributions of source code must retain the above copyright
+.\" notice, this list of conditions and the following disclaimer.
+.\" 2. Redistributions in binary form must reproduce the above copyright
+.\" notice, this list of conditions and the following disclaimer in the
+.\" documentation and/or other materials provided with the distribution.
+.\"
+.\" THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
+.\" ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+.\" TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+.\" PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR CONTRIBUTORS
+.\" BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+.\" CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+.\" SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+.\" INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+.\" CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+.\" ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+.\" POSSIBILITY OF SUCH DAMAGE.
+.\"
+.Dd September 18, 2011
+.Dt NEXTAFTER 3
+.Os
+.Sh NAME
+.Nm nextafter ,
+.Nm nextafterf ,
+.Nm nextafterl ,
+.Nm nexttoward
+.\"
+.\" XXX: Not yet implemented.
+.\"
+.\" .Nm nexttowardf ,
+.\" .Nm nexttowardl
+.\"
+.Nd next representable floating-point number
+.Sh LIBRARY
+.Lb libm
+.Sh SYNOPSIS
+.In math.h
+.Ft double
+.Fn nextafter "double x" "double y"
+.Ft float
+.Fn nextafterf "float x" "float y"
+.Ft long double
+.Fn nextafterl "long double x" "long double y"
+.Ft double
+.Fn nexttoward "double x" "long double y"
+.Sh DESCRIPTION
+The
+.Fn nextafter ,
+.Fn nextafterf ,
+and
+.Fn nextafterl
+functions return the next machine representable number from
+.Fa x
+in direction of
+.Fa y .
+In other words, if
+.Fa y
+is less than
+.Fa x ,
+the functions return the largest representable floating-point number less than
+.Fa x .
+When
+.Fa x
+equals
+.Fa y ,
+the value of
+.Fa y
+is returned.
+The three functions differ only in the type of the return value and
+.Fa x .
+.Pp
+The
+.Fn nexttoward
+function is equivalent to the
+.Fn nextafter
+family of functions with two exceptions:
+.Bl -enum -offset indent
+.It
+The second parameter has a type
+.Vt long double .
+.It
+The return value is
+.Fa y
+converted to the type of the function, provided that
+.Fa x
+equals
+.Fa y .
+.El
+.Sh RETURN VALUES
+Upon successful completion, the described functions return
+the next representable floating-point value as described above.
+If
+.Fa x
+is finite but an overflow would occur,
+a range error follows and the functions return
+.Dv \*(Pm\*HHUGE_VAL ,
+.Dv \*(Pm\*HHUGE_VALF ,
+or
+.Dv \*(Pm\*HHUGE_VALL
+with the same sign as
+.Fa x .
+When either
+.Fa x
+or
+.Fa y
+is \*(Na, a \*(Na is returned.
+When
+.Fa x
+is not
+.Fa y
+but the function value is subnormal, zero, or underflows,
+a range error occurs, and either 0.0 or the correct function
+value (if representable) is returned.
+.Sh SEE ALSO
+.Xr math 3
+.Sh STANDARDS
+The described functions conform to
+.St -isoC-99 .
--- /dev/null
+.\" $NetBSD: pow.3,v 1.1 2011/09/17 10:51:53 jruoho Exp $
+.\"
+.\" Copyright (c) 2011 Jukka Ruohonen <jruohonen@iki.fi>
+.\" All rights reserved.
+.\"
+.\" Redistribution and use in source and binary forms, with or without
+.\" modification, are permitted provided that the following conditions
+.\" are met:
+.\" 1. Redistributions of source code must retain the above copyright
+.\" notice, this list of conditions and the following disclaimer.
+.\" 2. Redistributions in binary form must reproduce the above copyright
+.\" notice, this list of conditions and the following disclaimer in the
+.\" documentation and/or other materials provided with the distribution.
+.\"
+.\" THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
+.\" ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+.\" TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+.\" PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR CONTRIBUTORS
+.\" BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+.\" CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+.\" SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+.\" INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+.\" CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+.\" ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+.\" POSSIBILITY OF SUCH DAMAGE.
+.\"
+.Dd September 13, 2011
+.Dt POW 3
+.Os
+.Sh NAME
+.Nm pow ,
+.Nm powf
+.Nd power functions
+.Sh LIBRARY
+.Lb libm
+.Sh SYNOPSIS
+.In math.h
+.Ft double
+.Fn pow "double x" "double y"
+.Ft float
+.Fn powf "float x" "float y"
+.Sh DESCRIPTION
+The
+.Fn pow
+family of functions return
+.Fa x
+raised to the power of
+.Fa y .
+.Pp
+If
+.Fa x
+is negative and
+.Fa y
+is not an integer, the global variable
+.Va errno
+is set to
+.Er EDOM ,
+and on
+.Tn VAX
+a reserved operand fault is generated.
+A portable application should nevertheless ensure that
+.Fa y
+is an integer value whenever
+.Fa x
+is negative.
+.Sh RETURN VALUES
+.\"
+.\" XXX: List also the special return values?
+.\"
+Upon successful completion, the described functions return
+.Fa x^y .
+.Sh SEE ALSO
+.Xr exp 3 ,
+.Xr log 3
+.Sh STANDARDS
+The described functions conform to
+.St -isoC-99 .
+.Sh HISTORY
+The history of the power functions dates back to
+.At v6 .
--- /dev/null
+.\" $NetBSD: remainder.3,v 1.2 2011/09/18 05:33:14 jruoho Exp $
+.\"
+.\" Copyright (c) 2011 Jukka Ruohonen <jruohonen@iki.fi>
+.\" All rights reserved.
+.\"
+.\" Redistribution and use in source and binary forms, with or without
+.\" modification, are permitted provided that the following conditions
+.\" are met:
+.\" 1. Redistributions of source code must retain the above copyright
+.\" notice, this list of conditions and the following disclaimer.
+.\" 2. Redistributions in binary form must reproduce the above copyright
+.\" notice, this list of conditions and the following disclaimer in the
+.\" documentation and/or other materials provided with the distribution.
+.\"
+.\" THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
+.\" ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+.\" TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+.\" PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR CONTRIBUTORS
+.\" BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+.\" CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+.\" SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+.\" INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+.\" CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+.\" ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+.\" POSSIBILITY OF SUCH DAMAGE.
+.\"
+.Dd September 18, 2011
+.Dt REMAINDER 3
+.Os
+.Sh NAME
+.Nm remainder ,
+.Nm remainderf ,
+.Nm remquo ,
+.Nm remquof
+.Nd remainder functions
+.Sh LIBRARY
+.Lb libm
+.Sh SYNOPSIS
+.In math.h
+.Ft double
+.Fn remainder "double x" "double y"
+.Ft float
+.Fn remainderf "float x" "float y"
+.Ft double
+.Fn remquo "double x" "double y" "int *quo"
+.Ft float
+.Fn remquof "float x" "float y" "int *quo"
+.Sh DESCRIPTION
+Provided that
+.Fa y
+\*(Ne 0 ,
+the
+.Fn remainder
+and
+.Fn remainderf
+functions calculate the floating-point remainder
+.Fa r
+of
+.Bd -ragged -offset indent
+.Va r
+=
+.Va x - ny ,
+.Ed
+.Pp
+where
+.Fa n
+is the integral value nearest to the exact value of
+.Fa x
+/
+.Fa y .
+If
+.Bd -ragged -offset indent
+.Va | n
+-
+.Va x / y |
+= 1/2 ,
+.Ed
+.Pp
+the value
+.Fa n
+is chosen to be even.
+Consequently, the remainder is computed exactly and
+.Va | r |
+\*(Le
+.Fa | y |
+/ 2 .
+.Pp
+Also the
+.Fn remquo
+and
+.Fn remquof
+functions calculate the remainder as described above.
+But these additionally use
+.Fa quo
+to store a value whose sign is the sign of
+.Va x / y
+and whose magnitude is congruent modulo
+.Va 2^k
+to the magnitude of the integral quotient of
+.Va x / y ,
+where
+.Fa k
+is an implementation-defined integer greater than or equal to 3.
+.Pp
+The rationale of the
+.Fn remquo
+family of functions relates to situations where
+only few bits of the quotient are required.
+The exact representation of the quotient may not be meaningful when
+.Fa x
+is large in magnitude compared to
+.Fa y .
+.Sh RETURN VALUES
+The functions return the remainder independent of the rounding mode.
+If
+.Fa y
+is zero ,
+\*(Na
+is returned and a domain error occurs.
+A domain error occurs and a
+\*(Na
+is returned also when
+.Fa x
+is infinite but
+.Fa y
+is not a
+\*(Na.
+If either
+.Fa x
+or
+.Fa y
+is
+\*(Na,
+a
+\*(Na
+is always returned.
+.Sh SEE ALSO
+.Xr div 3 ,
+.Xr fast_remainder32 3 ,
+.Xr fmod 3 ,
+.Xr math 3
+.Sh STANDARDS
+The described functions conform to
+.St -isoC-99 .
+.\" $NetBSD: round.3,v 1.6 2011/09/13 07:11:43 njoly Exp $
+.\"
.\" Copyright (c) 2003, Steven G. Kargl
.\" All rights reserved.
.\"
.\" SUCH DAMAGE.
.\"
.\" $FreeBSD: src/lib/msun/man/round.3,v 1.2 2004/06/20 09:27:17 das Exp $
-.\" $NetBSD: round.3,v 1.3 2004/07/15 12:12:39 junyoung Exp $
-.\" $DragonFly: src/lib/libm/man/round.3,v 1.1 2005/07/26 21:15:20 joerg Exp $
.\"
.Dd July 15, 2004
.Dt ROUND 3
.Sh SEE ALSO
.Xr ceil 3 ,
.Xr floor 3 ,
-.Xr ieee 3 ,
.Xr math 3 ,
-.Xr rint 3
-.\" .Xr trunc 3
+.Xr rint 3 ,
+.Xr trunc 3
.Sh STANDARDS
The
.Fn round
-function conforms to
+and
+.Fn roundf
+functions conform to
.St -isoC-99 .
.Sh HISTORY
The
--- /dev/null
+.\" $NetBSD: scalbn.3,v 1.2 2011/09/18 05:33:14 jruoho Exp $
+.\"
+.\" Copyright (c) 2011 Jukka Ruohonen <jruohonen@iki.fi>
+.\" All rights reserved.
+.\"
+.\" Redistribution and use in source and binary forms, with or without
+.\" modification, are permitted provided that the following conditions
+.\" are met:
+.\" 1. Redistributions of source code must retain the above copyright
+.\" notice, this list of conditions and the following disclaimer.
+.\" 2. Redistributions in binary form must reproduce the above copyright
+.\" notice, this list of conditions and the following disclaimer in the
+.\" documentation and/or other materials provided with the distribution.
+.\"
+.\" THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
+.\" ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+.\" TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+.\" PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR CONTRIBUTORS
+.\" BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+.\" CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+.\" SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+.\" INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+.\" CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+.\" ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+.\" POSSIBILITY OF SUCH DAMAGE.
+.\"
+.Dd September 18, 2011
+.Dt SCALBN 3
+.Os
+.Sh NAME
+.Nm scalbn ,
+.Nm scalbnf ,
+.Nm scalbnl
+.Nd exponent using FLT_RADIX
+.Sh LIBRARY
+.Lb libm
+.Sh SYNOPSIS
+.In math.h
+.Ft double
+.Fn scalbn "double x" "int n"
+.Ft float
+.Fn scalbnf "float x" "int n"
+.Ft long double
+.Fn scalbnl "long double x" "int n"
+.Sh DESCRIPTION
+The
+.Fn scalbn ,
+.Fn scalbnf ,
+and
+.Fn scalbnl
+functions compute
+.Fa x
+*
+.Fa r^n ,
+where
+.Fa r
+is the radix of the machine's floating point arithmetic, defined by the
+.Dv FLT_RADIX
+constant in
+.In float.h .
+The rationale is efficiency;
+.Fa r^n
+is not computed explicitly.
+.Sh RETURN VALUES
+As described above, upon successful completion, the described functions return
+the exponent computed using
+.Dv FLT_RADIX .
+Otherwise the following may occur:
+.Pp
+.Bl -enum -offset indent
+.It
+When the result would cause an overflow, a range error occurs and
+.Dv \*(Pm\*HHUGE_VAL ,
+.Dv \*(Pm\*HHUGE_VALF ,
+or
+.Dv \*(Pm\*HHUGE_VALL
+is returned according to the sign of
+.Fa x
+and the return type of the corresponding function.
+.It
+When the correct value would cause an underflow
+and it is not representable, a range error occurs and
+either 0.0 or an implementation-defined value is returned.
+When an underflow occurs but the correct value is representable,
+a range error occurs but the correct value is returned.
+.It
+If
+.Fa x
+is \*(Pm0 or \*(Pm\Inf,
+.Fa x
+is returned.
+Likewise, if
+.Fa n
+is zero,
+.Fa x
+is returned.
+If
+.Fa x
+is \*(Na, \*(Na is returned.
+.El
+.Sh SEE ALSO
+.Xr exp 3 ,
+.Xr frexp 3 ,
+.Xr ldexp 3 ,
+.Xr math 3
+.Sh STANDARDS
+The described functions conform to
+.St -isoC-99 .
.\" SUCH DAMAGE.
.\"
.\" from: @(#)tanh.3 5.1 (Berkeley) 5/2/91
-.\" $NetBSD: tanh.3,v 1.12 2003/08/07 16:44:49 agc Exp $
-.\" $DragonFly: src/lib/libm/man/tanh.3,v 1.1 2005/07/26 21:15:20 joerg Exp $
+.\" $NetBSD: tanh.3,v 1.15 2011/11/17 23:46:32 wiz Exp $
.\"
-.Dd May 2, 1991
+.Dd September 18, 2011
.Dt TANH 3
.Os
.Sh NAME
For a discussion of error due to roundoff, see
.Xr math 3 .
.Sh RETURN VALUES
-The
-.Fn tanh
-function returns the hyperbolic tangent value.
+Upon successful completion,
+these functions return the hyperbolic tangent value.
+The following may also occur:
+.Bl -enum -offset indent
+.It
+If
+.Fa x
+is \*(Pm 0,
+.Fa x
+is returned.
+.It
+If
+.Fa x
+is \*(Na, a \*(Na is returned.
+.It
+If
+.Fa x
+is positive infinity, a value 1 is returned;
+if
+.Fa x
+is negative infinity, -1 is returned.
+.It
+If
+.Fa x
+is subnormal, a range error can occur and
+.Fa x
+is returned.
+.El
.Sh SEE ALSO
.Xr acos 3 ,
.Xr asin 3 ,
.Xr sinh 3 ,
.Xr tan 3
.Sh STANDARDS
-The
-.Fn tanh
-function conforms to
-.St -ansiC .
+The described functions conform to
+.St -isoC-99 .
s_round.c s_roundf.c s_scalbn.c s_scalbnf.c \
s_signgam.c s_significand.c s_significandf.c s_sin.c s_sinf.c \
s_tan.c s_tanf.c s_tanh.c s_tanhf.c
+MI_FUNCS+= s_exp2.c s_exp2f.c s_frexp.c s_ilogbl.c s_logbl.c s_nextafterl.c s_nexttoward.c \
+ s_remquo.c s_remquof.c s_scalbnl.c s_fabsl.c
.for f in ${MI_FUNCS}
.if empty(SRCS:M${f:H:R})
* is preserved.
* ====================================================
*
- * $NetBSD: e_hypot.c,v 1.12 2002/05/26 22:01:50 wiz Exp $
- * $DragonFly: src/lib/libm/src/e_hypot.c,v 1.1 2005/07/26 21:15:20 joerg Exp $
+ * $NetBSD: e_hypot.c,v 1.13 2008/04/25 22:21:53 christos
*/
/* hypot(x,y)
* x1*x1+(y*y+(x2*(x+x1))) for x*x+y*y
* where x1 = x with lower 32 bits cleared, x2 = x-x1; else
* 2. if x <= 2y use
- * t1*y1+((x-y)*(x-y)+(t1*y2+t2*y))
+ * t1*yy1+((x-y)*(x-y)+(t1*y2+t2*y))
* where t1 = 2x with lower 32 bits cleared, t2 = 2x-t1,
- * y1= y with lower 32 bits chopped, y2 = y-y1.
+ * yy1= y with lower 32 bits chopped, y2 = y-yy1.
*
* NOTE: scaling may be necessary if some argument is too
* large or too tiny
double
hypot(double x, double y)
{
- double a,b,t1,t2,y1_,y2,w;
+ double a=x,b=y,t1,t2,yy1,y2,w;
int32_t j,k,ha,hb;
GET_HIGH_WORD(ha,x);
w = sqrt(t1*t1-(b*(-b)-t2*(a+t1)));
} else {
a = a+a;
- y1_ = 0;
- SET_HIGH_WORD(y1_,hb);
- y2 = b - y1_;
+ yy1 = 0;
+ SET_HIGH_WORD(yy1,hb);
+ y2 = b - yy1;
t1 = 0;
SET_HIGH_WORD(t1,ha+0x00100000);
t2 = a - t1;
- w = sqrt(t1*y1_-(w*(-w)-(t1*y2+t2*b)));
+ w = sqrt(t1*yy1-(w*(-w)-(t1*y2+t2*b)));
}
if(k!=0) {
u_int32_t high;
* is preserved.
* ====================================================
*
- * $NetBSD: e_hypotf.c,v 1.8 2002/05/26 22:01:50 wiz Exp $
- * $DragonFly: src/lib/libm/src/e_hypotf.c,v 1.1 2005/07/26 21:15:20 joerg Exp $
+ * $NetBSD: e_hypotf.c,v 1.9 2008/04/25 22:21:53 christos Exp $
*/
#include <math.h>
float
hypotf(float x, float y)
{
- float a,b,t1,t2,y1_,y2,w;
+ float a=x,b=y,t1,t2,yy1,y2,w;
int32_t j,k,ha,hb;
GET_FLOAT_WORD(ha,x);
ha &= 0x7fffffff;
GET_FLOAT_WORD(hb,y);
hb &= 0x7fffffff;
- if(hb > ha) {j=ha;ha=hb;hb=j;}
+ if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;}
SET_FLOAT_WORD(a,ha); /* a <- |a| */
SET_FLOAT_WORD(b,hb); /* b <- |b| */
if((ha-hb)>0xf000000) {return a+b;} /* x/y > 2**30 */
w = sqrtf(t1*t1-(b*(-b)-t2*(a+t1)));
} else {
a = a+a;
- SET_FLOAT_WORD(y1_,hb&0xfffff000);
- y2 = b - y1_;
+ SET_FLOAT_WORD(yy1,hb&0xfffff000);
+ y2 = b - yy1;
SET_FLOAT_WORD(t1,ha+0x00800000);
t2 = a - t1;
- w = sqrtf(t1*y1_-(w*(-w)-(t1*y2+t2*b)));
+ w = sqrtf(t1*yy1-(w*(-w)-(t1*y2+t2*b)));
}
if(k!=0) {
SET_FLOAT_WORD(t1,0x3f800000+(k<<23));
* is preserved.
* ====================================================
*
- * $NetBSD: e_j0.c,v 1.11 2002/05/26 22:01:50 wiz Exp $
- * $DragonFly: src/lib/libm/src/e_j0.c,v 1.1 2005/07/26 21:15:20 joerg Exp $
+ * $NetBSD: e_j0.c,v 1.12 2007/08/20 16:01:38 drochner Exp $
*/
/* j0(x), y0(x)
* is preserved.
* ====================================================
*
- * $NetBSD: e_j0f.c,v 1.9 2006/03/19 20:42:44 christos Exp $
- * $DragonFly: src/lib/libm/src/e_j0f.c,v 1.3 2007/06/17 00:53:52 pavalos Exp $
+ * $NetBSD: e_j0f.c,v 1.10 2007/08/20 16:01:38 drochner Exp $
*/
#include <math.h>
* is preserved.
* ====================================================
*
- * $NetBSD: e_j1.c,v 1.11 2002/05/26 22:01:50 wiz Exp $
- * $DragonFly: src/lib/libm/src/e_j1.c,v 1.1 2005/07/26 21:15:20 joerg Exp $
+ * $NetBSD: e_j1.c,v 1.12 2007/08/20 16:01:38 drochner Exp $
*/
/* j1(x), y1(x)
* is preserved.
* ====================================================
*
- * $NetBSD: e_j1f.c,v 1.10 2006/03/19 20:54:15 christos Exp $
- * $DragonFly: src/lib/libm/src/e_j1f.c,v 1.2 2007/06/17 01:09:00 pavalos Exp $
+ * $NetBSD: e_j1f.c,v 1.11 2007/08/20 16:01:38 drochner Exp $
*/
#include <math.h>
* is preserved.
* ====================================================
*
- * $NetBSD: e_jn.c,v 1.12 2002/05/26 22:01:50 wiz Exp $
- * $DragonFly: src/lib/libm/src/e_jn.c,v 1.1 2005/07/26 21:15:20 joerg Exp $
+ * $NetBSD: e_jn.c,v 1.14 2010/11/29 15:10:06 drochner Exp $
*/
/*
}
}
}
- b = (t*j0(x)/b);
+ z = j0(x);
+ w = j1(x);
+ if (fabs(z) >= fabs(w))
+ b = (t*z/b);
+ else
+ b = (t*w/a);
}
}
if(sgn==1) return -b; else return b;
* is preserved.
* ====================================================
*
- * $NetBSD: e_jnf.c,v 1.9 2002/05/26 22:01:50 wiz Exp $
- * $DragonFly: src/lib/libm/src/e_jnf.c,v 1.1 2005/07/26 21:15:20 joerg Exp $
+ * $NetBSD: e_jnf.c,v 1.11 2010/11/29 15:10:06 drochner Exp $
*/
#include <math.h>
}
}
}
- b = (t*j0f(x)/b);
+ z = j0f(x);
+ w = j1f(x);
+ if (fabsf(z) >= fabsf(w))
+ b = (t*z/b);
+ else
+ b = (t*w/a);
}
}
if(sgn==1) return -b; else return b;
b = y1f(x);
/* quit if b is -inf */
GET_FLOAT_WORD(ib,b);
- for(i=1;i<n&&ib!=0xff800000;i++){
+ for(i=1;i<n&&(uint32_t)ib!=0xff800000;i++){
temp = b;
b = ((float)(i+i)/x)*b - a;
GET_FLOAT_WORD(ib,b);
* is preserved.
* ====================================================
*
- * $NetBSD: e_pow.c,v 1.13 2004/06/30 18:43:15 drochner Exp $
- * $DragonFly: src/lib/libm/src/e_pow.c,v 1.1 2005/07/26 21:15:20 joerg Exp $
+ * $NetBSD: e_pow.c,v 1.16 2010/04/23 19:17:07 drochner Exp $
*/
/* pow(x,y) return x**y
pow(double x, double y)
{
double z,ax,z_h,z_l,p_h,p_l;
- double y1_,t1,t2,r,s,t,u,v,w;
+ double yy1,t1,t2,r,s,t,u,v,w;
int32_t i,j,k,yisint,n;
int32_t hx,hy,ix,iy;
u_int32_t lx,ly;
k = (iy>>20)-0x3ff; /* exponent */
if(k>20) {
j = ly>>(52-k);
- if((j<<(52-k))==ly) yisint = 2-(j&1);
+ if((uint32_t)(j<<(52-k))==ly) yisint = 2-(j&1);
} else if(ly==0) {
j = iy>>(20-k);
if((j<<(20-k))==iy) yisint = 2-(j&1);
t2 = z_l-(((t1-t)-dp_h[k])-z_h);
}
- /* split up y into y1_+y2 and compute (y1_+y2)*(t1+t2) */
- y1_ = y;
- SET_LOW_WORD(y1_,0);
- p_l = (y-y1_)*t1+y*t2;
- p_h = y1_*t1;
+ /* split up y into yy1+y2 and compute (yy1+y2)*(t1+t2) */
+ yy1 = y;
+ SET_LOW_WORD(yy1,0);
+ p_l = (y-yy1)*t1+y*t2;
+ p_h = yy1*t1;
z = p_l+p_h;
EXTRACT_WORDS(j,i,z);
if (j>=0x40900000) { /* z >= 1024 */
* is preserved.
* ====================================================
*
- * $NetBSD: e_powf.c,v 1.12 2006/03/19 20:46:25 christos Exp $
- * $DragonFly: src/lib/libm/src/e_powf.c,v 1.2 2007/06/17 00:57:06 pavalos Exp $
+ * $NetBSD: e_powf.c,v 1.15 2010/04/23 19:17:07 drochner Exp $
*/
#include <math.h>
powf(float x, float y)
{
float z,ax,z_h,z_l,p_h,p_l;
- float y1_,t1,t2,r,s,t,u,v,w;
+ float yy1,t1,t2,r,s,t,u,v,w;
int32_t i,j,k,yisint,n;
int32_t hx,hy,ix,iy,is;
if(((((u_int32_t)hx>>31)-1)|(yisint-1))==0)
s = -one; /* (-ve)**(odd int) */
- /* split up y into y1_+y2 and compute (y1_+y2)*(t1+t2) */
+ /* split up y into yy1+y2 and compute (yy1+y2)*(t1+t2) */
GET_FLOAT_WORD(is,y);
- SET_FLOAT_WORD(y1_,is&0xfffff000);
- p_l = (y-y1_)*t1+y*t2;
- p_h = y1_*t1;
+ SET_FLOAT_WORD(yy1,is&0xfffff000);
+ p_l = (y-yy1)*t1+y*t2;
+ p_h = yy1*t1;
z = p_l+p_h;
GET_FLOAT_WORD(j,z);
if (j>0x43000000) /* if z > 128 */
else if (j==0x43000000) { /* if z == 128 */
if(p_l+ovt>z-p_h) return s*huge*huge; /* overflow */
}
- else if (j==0xc3160000){ /* z == -150 */
+ else if ((uint32_t)j==0xc3160000){ /* z == -150 */
if(p_l<=z-p_h) return s*tiny*tiny; /* underflow */
}
else if ((j&0x7fffffff)>0x43160000) /* z <= -150 */
* is preserved.
* ====================================================
*
- * $NetBSD: e_rem_pio2f.c,v 1.8 2002/05/26 22:01:52 wiz Exp $
- * $DragonFly: src/lib/libm/src/e_rem_pio2f.c,v 1.1 2005/07/26 21:15:20 joerg Exp $
+ * $NetBSD: e_rem_pio2f.c,v 1.9 2009/01/19 06:00:30 lukem Exp $
*/
/* __libm_rem_pio2f(x,y)
fn = (float)n;
r = t-fn*pio2_1;
w = fn*pio2_1t; /* 1st round good to 40 bit */
- if(n<32&&(ix&0xffffff00)!=npio2_hw[n-1]) {
+ if(n<32&&(int32_t)(ix&0xffffff00)!=npio2_hw[n-1]) {
y[0] = r-w; /* quick check no cancellation */
} else {
u_int32_t high;
* is preserved.
* ====================================================
*
- * $NetBSD: e_scalb.c,v 1.9 2002/05/26 22:01:52 wiz Exp $
- * $DragonFly: src/lib/libm/src/e_scalb.c,v 1.1 2005/07/26 21:15:20 joerg Exp $
+ * $NetBSD: e_scalb.c,v 1.10 2010/04/23 19:17:07 drochner Exp $
*/
/*
* is preserved.
* ====================================================
*
- * $NetBSD: e_scalbf.c,v 1.6 2002/05/26 22:01:52 wiz Exp $
- * $DragonFly: src/lib/libm/src/e_scalbf.c,v 1.1 2005/07/26 21:15:20 joerg Exp $
+ * $NetBSD: e_scalbf.c,v 1.7 2010/04/23 19:17:07 drochner Exp $
*/
#include <math.h>
* is preserved.
* ====================================================
*
- * $NetBSD: e_sqrt.c,v 1.12 2002/05/26 22:01:52 wiz Exp $
- * $DragonFly: src/lib/libm/src/e_sqrt.c,v 1.1 2005/07/26 21:15:20 joerg Exp $
+ * $NetBSD: e_sqrt.c,v 1.13 2009/02/16 01:19:34 lukem Exp $
*/
/* sqrt(x)
t = s0;
if((t<ix0)||((t==ix0)&&(t1<=ix1))) {
s1 = t1+r;
- if(((t1&sign)==sign)&&(s1&sign)==0) s0 += 1;
+ if(((t1&sign)==(u_int32_t)sign)&&(s1&sign)==0) s0 += 1;
ix0 -= t;
if (ix1 < t1) ix0 -= 1;
ix1 -= t1;
* is preserved.
* ====================================================
*
- * $NetBSD: k_rem_pio2.c,v 1.11 2003/01/04 23:43:03 wiz Exp $
- * $DragonFly: src/lib/libm/src/k_rem_pio2.c,v 1.1 2005/07/26 21:15:20 joerg Exp $
+ * $NetBSD: k_rem_pio2.c,v 1.12 2010/04/23 19:17:07 drochner Exp $
*/
/*
* is preserved.
* ====================================================
*
- * $NetBSD: k_rem_pio2f.c,v 1.7 2002/05/26 22:01:53 wiz Exp $
- * $DragonFly: src/lib/libm/src/k_rem_pio2f.c,v 1.1 2005/07/26 21:15:20 joerg Exp $
+ * $NetBSD: k_rem_pio2f.c,v 1.8 2010/04/23 19:17:07 drochner Exp $
*/
#include <math.h>
+/* $NetBSD: lrint.c,v 1.4 2008/04/26 23:49:50 christos Exp $ */
+
/*-
* Copyright (c) 2004
* Matthias Drochner. All rights reserved.
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
- *
- * $NetBSD: lrint.c,v 1.3 2004/10/13 15:18:32 drochner Exp $
- * $DragonFly: src/lib/libm/src/lrint.c,v 1.1 2005/07/26 21:15:20 joerg Exp $
*/
#include <math.h>
GET_HIGH_WORD(i0, x);
e = i0 >> 20;
- s = e >> DBL_EXPBITS;
+ s = (uint32_t)e >> DBL_EXPBITS;
e = (e & 0x7ff) - DBL_EXP_BIAS;
/* 1.0 x 2^-1 is the smallest number which can be rounded to 1 */
+/* $NetBSD: lrintf.c,v 1.5 2008/04/26 23:49:50 christos Exp $ */
+
/*-
* Copyright (c) 2004
* Matthias Drochner. All rights reserved.
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
- *
- * $NetBSD: lrintf.c,v 1.4 2006/08/01 20:14:35 drochner Exp $
- * $DragonFly: src/lib/libm/src/lrintf.c,v 1.2 2007/06/17 02:27:53 pavalos Exp $
*/
#include <math.h>
GET_FLOAT_WORD(i0, x);
e = i0 >> SNG_FRACBITS;
- s = e >> SNG_EXPBITS;
+ s = (uint32_t)e >> SNG_EXPBITS;
e = (e & 0xff) - SNG_EXP_BIAS;
/* 1.0 x 2^-1 is the smallest number which can be rounded to 1 */
+/* $NetBSD: lround.c,v 1.3 2008/04/26 23:49:50 christos Exp $ */
+
/*-
* Copyright (c) 2004
* Matthias Drochner. All rights reserved.
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
- *
- * $NetBSD: lround.c,v 1.2 2004/10/13 15:18:32 drochner Exp $
- * $DragonFly: src/lib/libm/src/lround.c,v 1.1 2005/07/26 21:15:20 joerg Exp $
*/
#include <math.h>
GET_HIGH_WORD(i0, x);
e = i0 >> 20;
- s = e >> DBL_EXPBITS;
+ s = (uint32_t)e >> DBL_EXPBITS;
e = (e & 0x7ff) - DBL_EXP_BIAS;
/* 1.0 x 2^-1 is the smallest number which can be rounded to 1 */
+/* $NetBSD: lroundf.c,v 1.3 2008/04/26 23:49:50 christos Exp $ */
+
/*-
* Copyright (c) 2004
* Matthias Drochner. All rights reserved.
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
- *
- * $NetBSD: lroundf.c,v 1.2 2004/10/13 15:18:32 drochner Exp $
- * $DragonFly: src/lib/libm/src/lroundf.c,v 1.1 2005/07/26 21:15:20 joerg Exp $
*/
#include <math.h>
GET_FLOAT_WORD(i0, x);
e = i0 >> SNG_FRACBITS;
- s = e >> SNG_EXPBITS;
+ s = (uint32_t)e >> SNG_EXPBITS;
e = (e & 0xff) - SNG_EXP_BIAS;
/* 1.0 x 2^-1 is the smallest number which can be rounded to 1 */
/*
* from: @(#)fdlibm.h 5.1 93/09/24
- * $NetBSD: math_private.h,v 1.11 2001/02/21 18:09:26 bjh21 Exp $
- * $DragonFly: src/lib/libm/src/math_private.h,v 1.1 2005/07/26 21:15:20 joerg Exp $
+ * $NetBSD: math_private.h,v 1.16 2010/09/16 20:39:50 drochner Exp $
*/
#ifndef _MATH_PRIVATE_H_
ew_u.value = (d); \
(ix0) = ew_u.parts.msw; \
(ix1) = ew_u.parts.lsw; \
-} while (0)
+} while (/*CONSTCOND*/0)
/* Get the more significant 32 bit int from a double. */
ieee_double_shape_type gh_u; \
gh_u.value = (d); \
(i) = gh_u.parts.msw; \
-} while (0)
+} while (/*CONSTCOND*/0)
/* Get the less significant 32 bit int from a double. */
ieee_double_shape_type gl_u; \
gl_u.value = (d); \
(i) = gl_u.parts.lsw; \
-} while (0)
+} while (/*CONSTCOND*/0)
/* Set a double from two 32 bit ints. */
iw_u.parts.msw = (ix0); \
iw_u.parts.lsw = (ix1); \
(d) = iw_u.value; \
-} while (0)
+} while (/*CONSTCOND*/0)
/* Set the more significant 32 bits of a double from an int. */
sh_u.value = (d); \
sh_u.parts.msw = (v); \
(d) = sh_u.value; \
-} while (0)
+} while (/*CONSTCOND*/0)
/* Set the less significant 32 bits of a double from an int. */
sl_u.value = (d); \
sl_u.parts.lsw = (v); \
(d) = sl_u.value; \
-} while (0)
+} while (/*CONSTCOND*/0)
/* A union which permits us to convert between a float and a 32 bit
int. */
ieee_float_shape_type gf_u; \
gf_u.value = (d); \
(i) = gf_u.word; \
-} while (0)
+} while (/*CONSTCOND*/0)
/* Set a float from a 32 bit int. */
ieee_float_shape_type sf_u; \
sf_u.word = (i); \
(d) = sf_u.value; \
-} while (0)
+} while (/*CONSTCOND*/0)
-#ifdef _COMPLEX_H
+/*
+ * Attempt to get strict C99 semantics for assignment with non-C99 compilers.
+ */
+#if FLT_EVAL_METHOD == 0 || __GNUC__ == 0
+#define STRICT_ASSIGN(type, lval, rval) ((lval) = (rval))
+#else
+#define STRICT_ASSIGN(type, lval, rval) do { \
+ volatile type __lval; \
+ \
+ if (sizeof(type) >= sizeof(double)) \
+ (lval) = (rval); \
+ else { \
+ __lval = (rval); \
+ (lval) = __lval; \
+ } \
+} while (/*CONSTCOND*/0)
+#endif
+
+#ifdef _COMPLEX_H
/*
- * C99 specifies that complex numbers have the same representation as
- * an array of two elements, where the first element is the real part
- * and the second element is the imaginary part.
+ * Quoting from ISO/IEC 9899:TC2:
+ *
+ * 6.2.5.13 Types
+ * Each complex type has the same representation and alignment requirements as
+ * an array type containing exactly two elements of the corresponding real type;
+ * the first element is equal to the real part, and the second element to the
+ * imaginary part, of the complex number.
*/
typedef union {
- float complex f;
- float a[2];
+ float complex z;
+ float parts[2];
} float_complex;
+
typedef union {
- double complex f;
- double a[2];
+ double complex z;
+ double parts[2];
} double_complex;
+
typedef union {
- long double complex f;
- long double a[2];
+ long double complex z;
+ long double parts[2];
} long_double_complex;
-#define REALPART(z) ((z).a[0])
-#define IMAGPART(z) ((z).a[1])
+
+#define REAL_PART(z) ((z).parts[0])
+#define IMAG_PART(z) ((z).parts[1])
/*
* Inline functions that can be used to construct complex values.
static __inline float complex
cpackf(float x, float y)
{
- float_complex z;
+ float_complex fc;
- REALPART(z) = x;
- IMAGPART(z) = y;
- return (z.f);
+ REAL_PART(fc) = x;
+ IMAG_PART(fc) = y;
+ return (fc.z);
}
static __inline double complex
cpack(double x, double y)
{
- double_complex z;
+ double_complex dc;
- REALPART(z) = x;
- IMAGPART(z) = y;
- return (z.f);
+ REAL_PART(dc) = x;
+ IMAG_PART(dc) = y;
+ return (dc.z);
}
static __inline long double complex
cpackl(long double x, long double y)
{
- long_double_complex z;
+ long_double_complex ldc;
- REALPART(z) = x;
- IMAGPART(z) = y;
- return (z.f);
+ REAL_PART(ldc) = x;
+ IMAG_PART(ldc) = y;
+ return (ldc.z);
}
-#endif /* _COMPLEX_H */
+#endif /* _COMPLEX_H */
__BEGIN_DECLS
#pragma GCC visibility push(hidden)
* is preserved.
* ====================================================
*
- * $NetBSD: s_ceil.c,v 1.11 2002/05/26 22:01:54 wiz Exp $
- * $DragonFly: src/lib/libm/src/s_ceil.c,v 1.1 2005/07/26 21:15:20 joerg Exp $
+ * $NetBSD: s_ceil.c,v 1.13 2009/02/16 01:22:18 lukem Exp $
*/
/*
double
ceil(double x)
{
- int32_t i0,i1,j0;
+ int32_t i0,i1,jj0;
u_int32_t i,j;
EXTRACT_WORDS(i0,i1,x);
- j0 = ((i0>>20)&0x7ff)-0x3ff;
- if(j0<20) {
- if(j0<0) { /* raise inexact if x != 0 */
+ jj0 = ((i0>>20)&0x7ff)-0x3ff;
+ if(jj0<20) {
+ if(jj0<0) { /* raise inexact if x != 0 */
if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */
if(i0<0) {i0=0x80000000;i1=0;}
else if((i0|i1)!=0) { i0=0x3ff00000;i1=0;}
}
} else {
- i = (0x000fffff)>>j0;
+ i = (0x000fffff)>>jj0;
if(((i0&i)|i1)==0) return x; /* x is integral */
if(huge+x>0.0) { /* raise inexact flag */
- if(i0>0) i0 += (0x00100000)>>j0;
+ if(i0>0) i0 += (0x00100000)>>jj0;
i0 &= (~i); i1=0;
}
}
- } else if (j0>51) {
- if(j0==0x400) return x+x; /* inf or NaN */
+ } else if (jj0>51) {
+ if(jj0==0x400) return x+x; /* inf or NaN */
else return x; /* x is integral */
} else {
- i = ((u_int32_t)(0xffffffff))>>(j0-20);
+ i = ((u_int32_t)(0xffffffff))>>(jj0-20);
if((i1&i)==0) return x; /* x is integral */
if(huge+x>0.0) { /* raise inexact flag */
if(i0>0) {
- if(j0==20) i0+=1;
+ if(jj0==20) i0+=1;
else {
- j = i1 + (1<<(52-j0));
- if(j<i1) i0+=1; /* got a carry */
+ j = i1 + (1<<(52-jj0));
+ if(j<(u_int32_t)i1) i0+=1; /* got a carry */
i1 = j;
}
}
* is preserved.
* ====================================================
*
- * $NetBSD: s_ceilf.c,v 1.7 2002/05/26 22:01:54 wiz Exp $
- * $DragonFly: src/lib/libm/src/s_ceilf.c,v 1.1 2005/07/26 21:15:20 joerg Exp $
+ * $NetBSD: s_ceilf.c,v 1.8 2008/04/25 22:21:53 christos Exp $
*/
#include <math.h>
float
ceilf(float x)
{
- int32_t i0,j0;
+ int32_t i0,jj0;
u_int32_t i;
GET_FLOAT_WORD(i0,x);
- j0 = ((i0>>23)&0xff)-0x7f;
- if(j0<23) {
- if(j0<0) { /* raise inexact if x != 0 */
+ jj0 = ((i0>>23)&0xff)-0x7f;
+ if(jj0<23) {
+ if(jj0<0) { /* raise inexact if x != 0 */
if(huge+x>(float)0.0) {/* return 0*sign(x) if |x|<1 */
if(i0<0) {i0=0x80000000;}
else if(i0!=0) { i0=0x3f800000;}
}
} else {
- i = (0x007fffff)>>j0;
+ i = (0x007fffff)>>jj0;
if((i0&i)==0) return x; /* x is integral */
if(huge+x>(float)0.0) { /* raise inexact flag */
- if(i0>0) i0 += (0x00800000)>>j0;
+ if(i0>0) i0 += (0x00800000)>>jj0;
i0 &= (~i);
}
}
} else {
- if(j0==0x80) return x+x; /* inf or NaN */
+ if(jj0==0x80) return x+x; /* inf or NaN */
else return x; /* x is integral */
}
SET_FLOAT_WORD(x,i0);
* is preserved.
* ====================================================
*
- * $NetBSD: s_cos.c,v 1.10 2002/05/26 22:01:54 wiz Exp $
- * $DragonFly: src/lib/libm/src/s_cos.c,v 1.1 2005/07/26 21:15:20 joerg Exp $
+ * $NetBSD: s_cos.c,v 1.11 2007/08/20 16:01:39 drochner Exp $
*/
/* cos(x)
* is preserved.
* ====================================================
*
- * $NetBSD: s_cosf.c,v 1.8 2002/05/26 22:01:55 wiz Exp $
- * $DragonFly: src/lib/libm/src/s_cosf.c,v 1.1 2005/07/26 21:15:20 joerg Exp $
+ * $NetBSD: s_cosf.c,v 1.9 2007/08/20 16:01:39 drochner Exp $
*/
#include <math.h>
--- /dev/null
+/*-
+ * Copyright (c) 2005 David Schultz <das@FreeBSD.ORG>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions and the following disclaimer in the
+ * documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ *
+ * $NetBSD: s_exp2.c,v 1.2 2010/01/11 23:38:24 christos Exp $
+ */
+
+#include <float.h>
+
+#include <math.h>
+#include "math_private.h"
+
+#define TBLBITS 8
+#define TBLSIZE (1 << TBLBITS)
+
+static const double
+ huge = 0x1p1000,
+ redux = 0x1.8p52 / TBLSIZE,
+ P1 = 0x1.62e42fefa39efp-1,
+ P2 = 0x1.ebfbdff82c575p-3,
+ P3 = 0x1.c6b08d704a0a6p-5,
+ P4 = 0x1.3b2ab88f70400p-7,
+ P5 = 0x1.5d88003875c74p-10;
+
+static volatile double twom1000 = 0x1p-1000;
+
+static const double tbl[TBLSIZE * 2] = {
+/* exp2(z + eps) eps */
+ 0x1.6a09e667f3d5dp-1, 0x1.9880p-44,
+ 0x1.6b052fa751744p-1, 0x1.8000p-50,
+ 0x1.6c012750bd9fep-1, -0x1.8780p-45,
+ 0x1.6cfdcddd476bfp-1, 0x1.ec00p-46,
+ 0x1.6dfb23c651a29p-1, -0x1.8000p-50,
+ 0x1.6ef9298593ae3p-1, -0x1.c000p-52,
+ 0x1.6ff7df9519386p-1, -0x1.fd80p-45,
+ 0x1.70f7466f42da3p-1, -0x1.c880p-45,
+ 0x1.71f75e8ec5fc3p-1, 0x1.3c00p-46,
+ 0x1.72f8286eacf05p-1, -0x1.8300p-44,
+ 0x1.73f9a48a58152p-1, -0x1.0c00p-47,
+ 0x1.74fbd35d7ccfcp-1, 0x1.f880p-45,
+ 0x1.75feb564267f1p-1, 0x1.3e00p-47,
+ 0x1.77024b1ab6d48p-1, -0x1.7d00p-45,
+ 0x1.780694fde5d38p-1, -0x1.d000p-50,
+ 0x1.790b938ac1d00p-1, 0x1.3000p-49,
+ 0x1.7a11473eb0178p-1, -0x1.d000p-49,
+ 0x1.7b17b0976d060p-1, 0x1.0400p-45,
+ 0x1.7c1ed0130c133p-1, 0x1.0000p-53,
+ 0x1.7d26a62ff8636p-1, -0x1.6900p-45,
+ 0x1.7e2f336cf4e3bp-1, -0x1.2e00p-47,
+ 0x1.7f3878491c3e8p-1, -0x1.4580p-45,
+ 0x1.80427543e1b4ep-1, 0x1.3000p-44,
+ 0x1.814d2add1071ap-1, 0x1.f000p-47,
+ 0x1.82589994ccd7ep-1, -0x1.1c00p-45,
+ 0x1.8364c1eb942d0p-1, 0x1.9d00p-45,
+ 0x1.8471a4623cab5p-1, 0x1.7100p-43,
+ 0x1.857f4179f5bbcp-1, 0x1.2600p-45,
+ 0x1.868d99b4491afp-1, -0x1.2c40p-44,
+ 0x1.879cad931a395p-1, -0x1.3000p-45,
+ 0x1.88ac7d98a65b8p-1, -0x1.a800p-45,
+ 0x1.89bd0a4785800p-1, -0x1.d000p-49,
+ 0x1.8ace5422aa223p-1, 0x1.3280p-44,
+ 0x1.8be05bad619fap-1, 0x1.2b40p-43,
+ 0x1.8cf3216b54383p-1, -0x1.ed00p-45,
+ 0x1.8e06a5e08664cp-1, -0x1.0500p-45,
+ 0x1.8f1ae99157807p-1, 0x1.8280p-45,
+ 0x1.902fed0282c0ep-1, -0x1.cb00p-46,
+ 0x1.9145b0b91ff96p-1, -0x1.5e00p-47,
+ 0x1.925c353aa2ff9p-1, 0x1.5400p-48,
+ 0x1.93737b0cdc64ap-1, 0x1.7200p-46,
+ 0x1.948b82b5f98aep-1, -0x1.9000p-47,
+ 0x1.95a44cbc852cbp-1, 0x1.5680p-45,
+ 0x1.96bdd9a766f21p-1, -0x1.6d00p-44,
+ 0x1.97d829fde4e2ap-1, -0x1.1000p-47,
+ 0x1.98f33e47a23a3p-1, 0x1.d000p-45,
+ 0x1.9a0f170ca0604p-1, -0x1.8a40p-44,
+ 0x1.9b2bb4d53ff89p-1, 0x1.55c0p-44,
+ 0x1.9c49182a3f15bp-1, 0x1.6b80p-45,
+ 0x1.9d674194bb8c5p-1, -0x1.c000p-49,
+ 0x1.9e86319e3238ep-1, 0x1.7d00p-46,
+ 0x1.9fa5e8d07f302p-1, 0x1.6400p-46,
+ 0x1.a0c667b5de54dp-1, -0x1.5000p-48,
+ 0x1.a1e7aed8eb8f6p-1, 0x1.9e00p-47,
+ 0x1.a309bec4a2e27p-1, 0x1.ad80p-45,
+ 0x1.a42c980460a5dp-1, -0x1.af00p-46,
+ 0x1.a5503b23e259bp-1, 0x1.b600p-47,
+ 0x1.a674a8af46213p-1, 0x1.8880p-44,
+ 0x1.a799e1330b3a7p-1, 0x1.1200p-46,
+ 0x1.a8bfe53c12e8dp-1, 0x1.6c00p-47,
+ 0x1.a9e6b5579fcd2p-1, -0x1.9b80p-45,
+ 0x1.ab0e521356fb8p-1, 0x1.b700p-45,
+ 0x1.ac36bbfd3f381p-1, 0x1.9000p-50,
+ 0x1.ad5ff3a3c2780p-1, 0x1.4000p-49,
+ 0x1.ae89f995ad2a3p-1, -0x1.c900p-45,
+ 0x1.afb4ce622f367p-1, 0x1.6500p-46,
+ 0x1.b0e07298db790p-1, 0x1.fd40p-45,
+ 0x1.b20ce6c9a89a9p-1, 0x1.2700p-46,
+ 0x1.b33a2b84f1a4bp-1, 0x1.d470p-43,
+ 0x1.b468415b747e7p-1, -0x1.8380p-44,
+ 0x1.b59728de5593ap-1, 0x1.8000p-54,
+ 0x1.b6c6e29f1c56ap-1, 0x1.ad00p-47,
+ 0x1.b7f76f2fb5e50p-1, 0x1.e800p-50,
+ 0x1.b928cf22749b2p-1, -0x1.4c00p-47,
+ 0x1.ba5b030a10603p-1, -0x1.d700p-47,
+ 0x1.bb8e0b79a6f66p-1, 0x1.d900p-47,
+ 0x1.bcc1e904bc1ffp-1, 0x1.2a00p-47,
+ 0x1.bdf69c3f3a16fp-1, -0x1.f780p-46,
+ 0x1.bf2c25bd71db8p-1, -0x1.0a00p-46,
+ 0x1.c06286141b2e9p-1, -0x1.1400p-46,
+ 0x1.c199bdd8552e0p-1, 0x1.be00p-47,
+ 0x1.c2d1cd9fa64eep-1, -0x1.9400p-47,
+ 0x1.c40ab5fffd02fp-1, -0x1.ed00p-47,
+ 0x1.c544778fafd15p-1, 0x1.9660p-44,
+ 0x1.c67f12e57d0cbp-1, -0x1.a100p-46,
+ 0x1.c7ba88988c1b6p-1, -0x1.8458p-42,
+ 0x1.c8f6d9406e733p-1, -0x1.a480p-46,
+ 0x1.ca3405751c4dfp-1, 0x1.b000p-51,
+ 0x1.cb720dcef9094p-1, 0x1.1400p-47,
+ 0x1.ccb0f2e6d1689p-1, 0x1.0200p-48,
+ 0x1.cdf0b555dc412p-1, 0x1.3600p-48,
+ 0x1.cf3155b5bab3bp-1, -0x1.6900p-47,
+ 0x1.d072d4a0789bcp-1, 0x1.9a00p-47,
+ 0x1.d1b532b08c8fap-1, -0x1.5e00p-46,
+ 0x1.d2f87080d8a85p-1, 0x1.d280p-46,
+ 0x1.d43c8eacaa203p-1, 0x1.1a00p-47,
+ 0x1.d5818dcfba491p-1, 0x1.f000p-50,
+ 0x1.d6c76e862e6a1p-1, -0x1.3a00p-47,
+ 0x1.d80e316c9834ep-1, -0x1.cd80p-47,
+ 0x1.d955d71ff6090p-1, 0x1.4c00p-48,
+ 0x1.da9e603db32aep-1, 0x1.f900p-48,
+ 0x1.dbe7cd63a8325p-1, 0x1.9800p-49,
+ 0x1.dd321f301b445p-1, -0x1.5200p-48,
+ 0x1.de7d5641c05bfp-1, -0x1.d700p-46,
+ 0x1.dfc97337b9aecp-1, -0x1.6140p-46,
+ 0x1.e11676b197d5ep-1, 0x1.b480p-47,
+ 0x1.e264614f5a3e7p-1, 0x1.0ce0p-43,
+ 0x1.e3b333b16ee5cp-1, 0x1.c680p-47,
+ 0x1.e502ee78b3fb4p-1, -0x1.9300p-47,
+ 0x1.e653924676d68p-1, -0x1.5000p-49,
+ 0x1.e7a51fbc74c44p-1, -0x1.7f80p-47,
+ 0x1.e8f7977cdb726p-1, -0x1.3700p-48,
+ 0x1.ea4afa2a490e8p-1, 0x1.5d00p-49,
+ 0x1.eb9f4867ccae4p-1, 0x1.61a0p-46,
+ 0x1.ecf482d8e680dp-1, 0x1.5500p-48,
+ 0x1.ee4aaa2188514p-1, 0x1.6400p-51,
+ 0x1.efa1bee615a13p-1, -0x1.e800p-49,
+ 0x1.f0f9c1cb64106p-1, -0x1.a880p-48,
+ 0x1.f252b376bb963p-1, -0x1.c900p-45,
+ 0x1.f3ac948dd7275p-1, 0x1.a000p-53,
+ 0x1.f50765b6e4524p-1, -0x1.4f00p-48,
+ 0x1.f6632798844fdp-1, 0x1.a800p-51,
+ 0x1.f7bfdad9cbe38p-1, 0x1.abc0p-48,
+ 0x1.f91d802243c82p-1, -0x1.4600p-50,
+ 0x1.fa7c1819e908ep-1, -0x1.b0c0p-47,
+ 0x1.fbdba3692d511p-1, -0x1.0e00p-51,
+ 0x1.fd3c22b8f7194p-1, -0x1.0de8p-46,
+ 0x1.fe9d96b2a23eep-1, 0x1.e430p-49,
+ 0x1.0000000000000p+0, 0x0.0000p+0,
+ 0x1.00b1afa5abcbep+0, -0x1.3400p-52,
+ 0x1.0163da9fb3303p+0, -0x1.2170p-46,
+ 0x1.02168143b0282p+0, 0x1.a400p-52,
+ 0x1.02c9a3e77806cp+0, 0x1.f980p-49,
+ 0x1.037d42e11bbcap+0, -0x1.7400p-51,
+ 0x1.04315e86e7f89p+0, 0x1.8300p-50,
+ 0x1.04e5f72f65467p+0, -0x1.a3f0p-46,
+ 0x1.059b0d315855ap+0, -0x1.2840p-47,
+ 0x1.0650a0e3c1f95p+0, 0x1.1600p-48,
+ 0x1.0706b29ddf71ap+0, 0x1.5240p-46,
+ 0x1.07bd42b72a82dp+0, -0x1.9a00p-49,
+ 0x1.0874518759bd0p+0, 0x1.6400p-49,
+ 0x1.092bdf66607c8p+0, -0x1.0780p-47,
+ 0x1.09e3ecac6f383p+0, -0x1.8000p-54,
+ 0x1.0a9c79b1f3930p+0, 0x1.fa00p-48,
+ 0x1.0b5586cf988fcp+0, -0x1.ac80p-48,
+ 0x1.0c0f145e46c8ap+0, 0x1.9c00p-50,
+ 0x1.0cc922b724816p+0, 0x1.5200p-47,
+ 0x1.0d83b23395dd8p+0, -0x1.ad00p-48,
+ 0x1.0e3ec32d3d1f3p+0, 0x1.bac0p-46,
+ 0x1.0efa55fdfa9a6p+0, -0x1.4e80p-47,
+ 0x1.0fb66affed2f0p+0, -0x1.d300p-47,
+ 0x1.1073028d7234bp+0, 0x1.1500p-48,
+ 0x1.11301d0125b5bp+0, 0x1.c000p-49,
+ 0x1.11edbab5e2af9p+0, 0x1.6bc0p-46,
+ 0x1.12abdc06c31d5p+0, 0x1.8400p-49,
+ 0x1.136a814f2047dp+0, -0x1.ed00p-47,
+ 0x1.1429aaea92de9p+0, 0x1.8e00p-49,
+ 0x1.14e95934f3138p+0, 0x1.b400p-49,
+ 0x1.15a98c8a58e71p+0, 0x1.5300p-47,
+ 0x1.166a45471c3dfp+0, 0x1.3380p-47,
+ 0x1.172b83c7d5211p+0, 0x1.8d40p-45,
+ 0x1.17ed48695bb9fp+0, -0x1.5d00p-47,
+ 0x1.18af9388c8d93p+0, -0x1.c880p-46,
+ 0x1.1972658375d66p+0, 0x1.1f00p-46,
+ 0x1.1a35beb6fcba7p+0, 0x1.0480p-46,
+ 0x1.1af99f81387e3p+0, -0x1.7390p-43,
+ 0x1.1bbe084045d54p+0, 0x1.4e40p-45,
+ 0x1.1c82f95281c43p+0, -0x1.a200p-47,
+ 0x1.1d4873168b9b2p+0, 0x1.3800p-49,
+ 0x1.1e0e75eb44031p+0, 0x1.ac00p-49,
+ 0x1.1ed5022fcd938p+0, 0x1.1900p-47,
+ 0x1.1f9c18438cdf7p+0, -0x1.b780p-46,
+ 0x1.2063b88628d8fp+0, 0x1.d940p-45,
+ 0x1.212be3578a81ep+0, 0x1.8000p-50,
+ 0x1.21f49917ddd41p+0, 0x1.b340p-45,
+ 0x1.22bdda2791323p+0, 0x1.9f80p-46,
+ 0x1.2387a6e7561e7p+0, -0x1.9c80p-46,
+ 0x1.2451ffb821427p+0, 0x1.2300p-47,
+ 0x1.251ce4fb2a602p+0, -0x1.3480p-46,
+ 0x1.25e85711eceb0p+0, 0x1.2700p-46,
+ 0x1.26b4565e27d16p+0, 0x1.1d00p-46,
+ 0x1.2780e341de00fp+0, 0x1.1ee0p-44,
+ 0x1.284dfe1f5633ep+0, -0x1.4c00p-46,
+ 0x1.291ba7591bb30p+0, -0x1.3d80p-46,
+ 0x1.29e9df51fdf09p+0, 0x1.8b00p-47,
+ 0x1.2ab8a66d10e9bp+0, -0x1.27c0p-45,
+ 0x1.2b87fd0dada3ap+0, 0x1.a340p-45,
+ 0x1.2c57e39771af9p+0, -0x1.0800p-46,
+ 0x1.2d285a6e402d9p+0, -0x1.ed00p-47,
+ 0x1.2df961f641579p+0, -0x1.4200p-48,
+ 0x1.2ecafa93e2ecfp+0, -0x1.4980p-45,
+ 0x1.2f9d24abd8822p+0, -0x1.6300p-46,
+ 0x1.306fe0a31b625p+0, -0x1.2360p-44,
+ 0x1.31432edeea50bp+0, -0x1.0df8p-40,
+ 0x1.32170fc4cd7b8p+0, -0x1.2480p-45,
+ 0x1.32eb83ba8e9a2p+0, -0x1.5980p-45,
+ 0x1.33c08b2641766p+0, 0x1.ed00p-46,
+ 0x1.3496266e3fa27p+0, -0x1.c000p-50,
+ 0x1.356c55f929f0fp+0, -0x1.0d80p-44,
+ 0x1.36431a2de88b9p+0, 0x1.2c80p-45,
+ 0x1.371a7373aaa39p+0, 0x1.0600p-45,
+ 0x1.37f26231e74fep+0, -0x1.6600p-46,
+ 0x1.38cae6d05d838p+0, -0x1.ae00p-47,
+ 0x1.39a401b713ec3p+0, -0x1.4720p-43,
+ 0x1.3a7db34e5a020p+0, 0x1.8200p-47,
+ 0x1.3b57fbfec6e95p+0, 0x1.e800p-44,
+ 0x1.3c32dc313a8f2p+0, 0x1.f800p-49,
+ 0x1.3d0e544ede122p+0, -0x1.7a00p-46,
+ 0x1.3dea64c1234bbp+0, 0x1.6300p-45,
+ 0x1.3ec70df1c4eccp+0, -0x1.8a60p-43,
+ 0x1.3fa4504ac7e8cp+0, -0x1.cdc0p-44,
+ 0x1.40822c367a0bbp+0, 0x1.5b80p-45,
+ 0x1.4160a21f72e95p+0, 0x1.ec00p-46,
+ 0x1.423fb27094646p+0, -0x1.3600p-46,
+ 0x1.431f5d950a920p+0, 0x1.3980p-45,
+ 0x1.43ffa3f84b9ebp+0, 0x1.a000p-48,
+ 0x1.44e0860618919p+0, -0x1.6c00p-48,
+ 0x1.45c2042a7d201p+0, -0x1.bc00p-47,
+ 0x1.46a41ed1d0016p+0, -0x1.2800p-46,
+ 0x1.4786d668b3326p+0, 0x1.0e00p-44,
+ 0x1.486a2b5c13c00p+0, -0x1.d400p-45,
+ 0x1.494e1e192af04p+0, 0x1.c200p-47,
+ 0x1.4a32af0d7d372p+0, -0x1.e500p-46,
+ 0x1.4b17dea6db801p+0, 0x1.7800p-47,
+ 0x1.4bfdad53629e1p+0, -0x1.3800p-46,
+ 0x1.4ce41b817c132p+0, 0x1.0800p-47,
+ 0x1.4dcb299fddddbp+0, 0x1.c700p-45,
+ 0x1.4eb2d81d8ab96p+0, -0x1.ce00p-46,
+ 0x1.4f9b2769d2d02p+0, 0x1.9200p-46,
+ 0x1.508417f4531c1p+0, -0x1.8c00p-47,
+ 0x1.516daa2cf662ap+0, -0x1.a000p-48,
+ 0x1.5257de83f51eap+0, 0x1.a080p-43,
+ 0x1.5342b569d4edap+0, -0x1.6d80p-45,
+ 0x1.542e2f4f6ac1ap+0, -0x1.2440p-44,
+ 0x1.551a4ca5d94dbp+0, 0x1.83c0p-43,
+ 0x1.56070dde9116bp+0, 0x1.4b00p-45,
+ 0x1.56f4736b529dep+0, 0x1.15a0p-43,
+ 0x1.57e27dbe2c40ep+0, -0x1.9e00p-45,
+ 0x1.58d12d497c76fp+0, -0x1.3080p-45,
+ 0x1.59c0827ff0b4cp+0, 0x1.dec0p-43,
+ 0x1.5ab07dd485427p+0, -0x1.4000p-51,
+ 0x1.5ba11fba87af4p+0, 0x1.0080p-44,
+ 0x1.5c9268a59460bp+0, -0x1.6c80p-45,
+ 0x1.5d84590998e3fp+0, 0x1.69a0p-43,
+ 0x1.5e76f15ad20e1p+0, -0x1.b400p-46,
+ 0x1.5f6a320dcebcap+0, 0x1.7700p-46,
+ 0x1.605e1b976dcb8p+0, 0x1.6f80p-45,
+ 0x1.6152ae6cdf715p+0, 0x1.1000p-47,
+ 0x1.6247eb03a5531p+0, -0x1.5d00p-46,
+ 0x1.633dd1d1929b5p+0, -0x1.2d00p-46,
+ 0x1.6434634ccc313p+0, -0x1.a800p-49,
+ 0x1.652b9febc8efap+0, -0x1.8600p-45,
+ 0x1.6623882553397p+0, 0x1.1fe0p-40,
+ 0x1.671c1c708328ep+0, -0x1.7200p-44,
+ 0x1.68155d44ca97ep+0, 0x1.6800p-49,
+ 0x1.690f4b19e9471p+0, -0x1.9780p-45,
+};
+
+/*
+ * exp2(x): compute the base 2 exponential of x
+ *
+ * Accuracy: Peak error < 0.503 ulp for normalized results.
+ *
+ * Method: (accurate tables)
+ *
+ * Reduce x:
+ * x = 2**k + y, for integer k and |y| <= 1/2.
+ * Thus we have exp2(x) = 2**k * exp2(y).
+ *
+ * Reduce y:
+ * y = i/TBLSIZE + z - eps[i] for integer i near y * TBLSIZE.
+ * Thus we have exp2(y) = exp2(i/TBLSIZE) * exp2(z - eps[i]),
+ * with |z - eps[i]| <= 2**-9 + 2**-39 for the table used.
+ *
+ * We compute exp2(i/TBLSIZE) via table lookup and exp2(z - eps[i]) via
+ * a degree-5 minimax polynomial with maximum error under 1.3 * 2**-61.
+ * The values in exp2t[] and eps[] are chosen such that
+ * exp2t[i] = exp2(i/TBLSIZE + eps[i]), and eps[i] is a small offset such
+ * that exp2t[i] is accurate to 2**-64.
+ *
+ * Note that the range of i is +-TBLSIZE/2, so we actually index the tables
+ * by i0 = i + TBLSIZE/2. For cache efficiency, exp2t[] and eps[] are
+ * virtual tables, interleaved in the real table tbl[].
+ *
+ * This method is due to Gal, with many details due to Gal and Bachelis:
+ *
+ * Gal, S. and Bachelis, B. An Accurate Elementary Mathematical Library
+ * for the IEEE Floating Point Standard. TOMS 17(1), 26-46 (1991).
+ */
+double
+exp2(double x)
+{
+ double r, t, twopk, twopkp1000, z;
+ uint32_t hx, ix, lx, i0;
+ int k;
+
+ /* Filter out exceptional cases. */
+ GET_HIGH_WORD(hx,x);
+ ix = hx & 0x7fffffff; /* high word of |x| */
+ if(ix >= 0x40900000) { /* |x| >= 1024 */
+ if(ix >= 0x7ff00000) {
+ GET_LOW_WORD(lx,x);
+ if(((ix & 0xfffff) | lx) != 0 || (hx & 0x80000000) == 0)
+ return (x + x); /* x is NaN or +Inf */
+ else
+ return (0.0); /* x is -Inf */
+ }
+ if(x >= 0x1.0p10)
+ return (huge * huge); /* overflow */
+ if(x <= -0x1.0ccp10)
+ return (twom1000 * twom1000); /* underflow */
+ } else if (ix < 0x3c900000) { /* |x| < 0x1p-54 */
+ return (1.0 + x);
+ }
+
+ /* Reduce x, computing z, i0, and k. */
+ STRICT_ASSIGN(double, t, x + redux);
+ GET_LOW_WORD(i0, t);
+ i0 += TBLSIZE / 2;
+ k = (i0 >> TBLBITS) << 20;
+ i0 = (i0 & (TBLSIZE - 1)) << 1;
+ t -= redux;
+ z = x - t;
+
+ /* Compute r = exp2(y) = exp2t[i0] * p(z - eps[i]). */
+ t = tbl[i0]; /* exp2t[i0] */
+ z -= tbl[i0 + 1]; /* eps[i0] */
+ if (k >= -1021 << 20)
+ INSERT_WORDS(twopk, 0x3ff00000 + k, 0);
+ else
+ INSERT_WORDS(twopkp1000, 0x3ff00000 + k + (1000 << 20), 0);
+ r = t + t * z * (P1 + z * (P2 + z * (P3 + z * (P4 + z * P5))));
+
+ /* Scale by 2**(k>>20). */
+ if(k >= -1021 << 20) {
+ if (k == 1024 << 20)
+ return (r * 2.0 * 0x1p1023);
+ return (r * twopk);
+ } else {
+ return (r * twopkp1000 * twom1000);
+ }
+}
+
+#ifdef notyet
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(exp2, exp2l);
+#endif
+#endif
--- /dev/null
+/*-
+ * Copyright (c) 2005 David Schultz <das@FreeBSD.ORG>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions and the following disclaimer in the
+ * documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ *
+ * $NetBSD: s_exp2f.c,v 1.1 2010/01/11 16:28:39 christos Exp $
+ */
+
+#include <float.h>
+
+#include <math.h>
+#include "math_private.h"
+
+#define TBLBITS 4
+#define TBLSIZE (1 << TBLBITS)
+
+static const float
+ huge = 0x1p100f,
+ redux = 0x1.8p23f / TBLSIZE,
+ P1 = 0x1.62e430p-1f,
+ P2 = 0x1.ebfbe0p-3f,
+ P3 = 0x1.c6b348p-5f,
+ P4 = 0x1.3b2c9cp-7f;
+
+static volatile float twom100 = 0x1p-100f;
+
+static const double exp2ft[TBLSIZE] = {
+ 0x1.6a09e667f3bcdp-1,
+ 0x1.7a11473eb0187p-1,
+ 0x1.8ace5422aa0dbp-1,
+ 0x1.9c49182a3f090p-1,
+ 0x1.ae89f995ad3adp-1,
+ 0x1.c199bdd85529cp-1,
+ 0x1.d5818dcfba487p-1,
+ 0x1.ea4afa2a490dap-1,
+ 0x1.0000000000000p+0,
+ 0x1.0b5586cf9890fp+0,
+ 0x1.172b83c7d517bp+0,
+ 0x1.2387a6e756238p+0,
+ 0x1.306fe0a31b715p+0,
+ 0x1.3dea64c123422p+0,
+ 0x1.4bfdad5362a27p+0,
+ 0x1.5ab07dd485429p+0,
+};
+
+/*
+ * exp2f(x): compute the base 2 exponential of x
+ *
+ * Accuracy: Peak error < 0.501 ulp; location of peak: -0.030110927.
+ *
+ * Method: (equally-spaced tables)
+ *
+ * Reduce x:
+ * x = 2**k + y, for integer k and |y| <= 1/2.
+ * Thus we have exp2f(x) = 2**k * exp2(y).
+ *
+ * Reduce y:
+ * y = i/TBLSIZE + z for integer i near y * TBLSIZE.
+ * Thus we have exp2(y) = exp2(i/TBLSIZE) * exp2(z),
+ * with |z| <= 2**-(TBLSIZE+1).
+ *
+ * We compute exp2(i/TBLSIZE) via table lookup and exp2(z) via a
+ * degree-4 minimax polynomial with maximum error under 1.4 * 2**-33.
+ * Using double precision for everything except the reduction makes
+ * roundoff error insignificant and simplifies the scaling step.
+ *
+ * This method is due to Tang, but I do not use his suggested parameters:
+ *
+ * Tang, P. Table-driven Implementation of the Exponential Function
+ * in IEEE Floating-Point Arithmetic. TOMS 15(2), 144-157 (1989).
+ */
+float
+exp2f(float x)
+{
+ double tv, twopk, u, z;
+ float t;
+ uint32_t hx, ix, i0;
+ int32_t k;
+
+ /* Filter out exceptional cases. */
+ GET_FLOAT_WORD(hx, x);
+ ix = hx & 0x7fffffff; /* high word of |x| */
+ if(ix >= 0x43000000) { /* |x| >= 128 */
+ if(ix >= 0x7f800000) {
+ if ((ix & 0x7fffff) != 0 || (hx & 0x80000000) == 0)
+ return (x + x); /* x is NaN or +Inf */
+ else
+ return (0.0); /* x is -Inf */
+ }
+ if(x >= 0x1.0p7f)
+ return (huge * huge); /* overflow */
+ if(x <= -0x1.2cp7f)
+ return (twom100 * twom100); /* underflow */
+ } else if (ix <= 0x33000000) { /* |x| <= 0x1p-25 */
+ return (1.0f + x);
+ }
+
+ /* Reduce x, computing z, i0, and k. */
+ STRICT_ASSIGN(float, t, x + redux);
+ GET_FLOAT_WORD(i0, t);
+ i0 += TBLSIZE / 2;
+ k = (i0 >> TBLBITS) << 20;
+ i0 &= TBLSIZE - 1;
+ t -= redux;
+ z = x - t;
+ INSERT_WORDS(twopk, 0x3ff00000 + k, 0);
+
+ /* Compute r = exp2(y) = exp2ft[i0] * p(z). */
+ tv = exp2ft[i0];
+ u = tv * z;
+ tv = tv + u * (P1 + z * P2) + u * (z * z) * (P3 + z * P4);
+
+ /* Scale by 2**(k>>20). */
+ return (tv * twopk);
+}
-/* $NetBSD: cacosh.c,v 1.1 2007/08/20 16:01:30 drochner Exp $ */
+/* $NetBSD: s_fabsl.c,v 1.2 2010/09/17 20:39:39 christos Exp $ */
/*-
- * Copyright (c) 2007 The NetBSD Foundation, Inc.
+ * Copyright (c) 2010 The NetBSD Foundation, Inc.
* All rights reserved.
*
- * This code is derived from software written by Stephen L. Moshier.
- * It is redistributed by the NetBSD Foundation by permission of the author.
- *
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
* POSSIBILITY OF SUCH DAMAGE.
+ *
+ * $NetBSD: s_fabsl.c,v 1.2 2010/09/17 20:39:39 christos Exp $
*/
-#include <complex.h>
+#include <math.h>
+#include <machine/ieee.h>
-double complex
-cacosh(double complex z)
+/*
+ * fabsl(long double x)
+ * This function returns the absolute value of its argumetn x, |x|.
+ */
+#ifdef EXT_EXP_INFNAN
+long double
+fabsl(long double x)
{
- double complex w;
+ union ieee_ext_u ux;
-#if 0 /* does not give the principal value */
- w = I * cacos(z);
-#else
- w = clog(z + csqrt(z + 1) * csqrt(z - 1));
-#endif
- return w;
+ ux.extu_ld = x;
+ ux.extu_ext.ext_sign = 0;
+
+ return (ux.extu_ld);
}
+#endif
* is preserved.
* ====================================================
*
- * $NetBSD: s_floor.c,v 1.11 2002/05/26 22:01:56 wiz Exp $
- * $DragonFly: src/lib/libm/src/s_floor.c,v 1.1 2005/07/26 21:15:20 joerg Exp $
+ * $NetBSD: s_floor.c,v 1.13 2009/02/16 01:27:36 lukem Exp $
*/
/*
double
floor(double x)
{
- int32_t i0,i1,j0;
+ int32_t i0,i1,jj0;
u_int32_t i,j;
EXTRACT_WORDS(i0,i1,x);
- j0 = ((i0>>20)&0x7ff)-0x3ff;
- if(j0<20) {
- if(j0<0) { /* raise inexact if x != 0 */
+ jj0 = ((i0>>20)&0x7ff)-0x3ff;
+ if(jj0<20) {
+ if(jj0<0) { /* raise inexact if x != 0 */
if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */
if(i0>=0) {i0=i1=0;}
else if(((i0&0x7fffffff)|i1)!=0)
{ i0=0xbff00000;i1=0;}
}
} else {
- i = (0x000fffff)>>j0;
+ i = (0x000fffff)>>jj0;
if(((i0&i)|i1)==0) return x; /* x is integral */
if(huge+x>0.0) { /* raise inexact flag */
- if(i0<0) i0 += (0x00100000)>>j0;
+ if(i0<0) i0 += (0x00100000)>>jj0;
i0 &= (~i); i1=0;
}
}
- } else if (j0>51) {
- if(j0==0x400) return x+x; /* inf or NaN */
+ } else if (jj0>51) {
+ if(jj0==0x400) return x+x; /* inf or NaN */
else return x; /* x is integral */
} else {
- i = ((u_int32_t)(0xffffffff))>>(j0-20);
+ i = ((u_int32_t)(0xffffffff))>>(jj0-20);
if((i1&i)==0) return x; /* x is integral */
if(huge+x>0.0) { /* raise inexact flag */
if(i0<0) {
- if(j0==20) i0+=1;
+ if(jj0==20) i0+=1;
else {
- j = i1+(1<<(52-j0));
- if(j<i1) i0 +=1 ; /* got a carry */
+ j = i1+(1<<(52-jj0));
+ if(j<(u_int32_t)i1) i0 +=1 ; /* got a carry */
i1=j;
}
}
* is preserved.
* ====================================================
*
- * $NetBSD: s_floorf.c,v 1.7 2002/05/26 22:01:56 wiz Exp $
- * $DragonFly: src/lib/libm/src/s_floorf.c,v 1.1 2005/07/26 21:15:20 joerg Exp $
+ * $NetBSD: s_floorf.c,v 1.8 2008/04/25 22:21:53 christos Exp $
*/
/*
float
floorf(float x)
{
- int32_t i0,j0;
+ int32_t i0,jj0;
u_int32_t i;
GET_FLOAT_WORD(i0,x);
- j0 = ((i0>>23)&0xff)-0x7f;
- if(j0<23) {
- if(j0<0) { /* raise inexact if x != 0 */
+ jj0 = ((i0>>23)&0xff)-0x7f;
+ if(jj0<23) {
+ if(jj0<0) { /* raise inexact if x != 0 */
if(huge+x>(float)0.0) {/* return 0*sign(x) if |x|<1 */
if(i0>=0) {i0=0;}
else if((i0&0x7fffffff)!=0)
{ i0=0xbf800000;}
}
} else {
- i = (0x007fffff)>>j0;
+ i = (0x007fffff)>>jj0;
if((i0&i)==0) return x; /* x is integral */
if(huge+x>(float)0.0) { /* raise inexact flag */
- if(i0<0) i0 += (0x00800000)>>j0;
+ if(i0<0) i0 += (0x00800000)>>jj0;
i0 &= (~i);
}
}
} else {
- if(j0==0x80) return x+x; /* inf or NaN */
+ if(jj0==0x80) return x+x; /* inf or NaN */
else return x; /* x is integral */
}
SET_FLOAT_WORD(x,i0);
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
+ *
+ * $NetBSD: s_fmax.c,v 1.2 2010/03/08 01:05:20 snj Exp $
*/
-#include <sys/cdefs.h>
-
#include <math.h>
-#include "fpmath.h"
+#include <machine/ieee.h>
double
fmax(double x, double y)
{
- union IEEEd2bits u[2];
+ union ieee_double_u u[2];
- u[0].d = x;
- u[1].d = y;
+ u[0].dblu_d = x;
+ u[1].dblu_d = y;
/* Check for NaNs to avoid raising spurious exceptions. */
- if (u[0].bits.exp == 2047 && (u[0].bits.manh | u[0].bits.manl) != 0)
+ if (u[0].dblu_dbl.dbl_exp == DBL_EXP_INFNAN &&
+ (u[0].dblu_dbl.dbl_frach | u[0].dblu_dbl.dbl_fracl) != 0)
return (y);
- if (u[1].bits.exp == 2047 && (u[1].bits.manh | u[1].bits.manl) != 0)
+ if (u[1].dblu_dbl.dbl_exp == DBL_EXP_INFNAN &&
+ (u[1].dblu_dbl.dbl_frach | u[1].dblu_dbl.dbl_fracl) != 0)
return (x);
/* Handle comparisons of signed zeroes. */
- if (u[0].bits.sign != u[1].bits.sign)
- return (u[u[0].bits.sign].d);
+ if (u[0].dblu_dbl.dbl_sign != u[1].dblu_dbl.dbl_sign)
+ return (u[u[0].dblu_dbl.dbl_sign].dblu_d);
return (x > y ? x : y);
}
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
+ *
+ * $NetBSD: s_fmaxf.c,v 1.2 2010/03/08 01:05:20 snj Exp $
*/
-#include <sys/cdefs.h>
-
#include <math.h>
-#include "fpmath.h"
+#include <machine/ieee.h>
float
fmaxf(float x, float y)
{
- union IEEEf2bits u[2];
+ union ieee_single_u u[2];
- u[0].f = x;
- u[1].f = y;
+ u[0].sngu_f = x;
+ u[1].sngu_f = y;
/* Check for NaNs to avoid raising spurious exceptions. */
- if (u[0].bits.exp == 255 && u[0].bits.man != 0)
+ if (u[0].sngu_sng.sng_exp == SNG_EXP_INFNAN &&
+ u[0].sngu_sng.sng_frac != 0)
return (y);
- if (u[1].bits.exp == 255 && u[1].bits.man != 0)
+ if (u[1].sngu_sng.sng_exp == SNG_EXP_INFNAN &&
+ u[1].sngu_sng.sng_frac != 0)
return (x);
- /* Handle comparisons of signed zeroes. */
- if (u[0].bits.sign != u[1].bits.sign)
- return (u[u[0].bits.sign].f);
+ /* Handle comparisons of sng_signed zeroes. */
+ if (u[0].sngu_sng.sng_sign != u[1].sngu_sng.sng_sign)
+ return (u[u[0].sngu_sng.sng_sign].sngu_f);
return (x > y ? x : y);
}
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
+ *
+ * $NetBSD: s_fmaxl.c,v 1.3 2011/07/04 11:46:41 mrg Exp $
*/
-#include <sys/cdefs.h>
-
+#include <string.h>
#include <math.h>
-#include "fpmath.h"
-
+#include <machine/ieee.h>
+#ifdef EXT_EXP_INFNAN
long double
fmaxl(long double x, long double y)
{
- union IEEEl2bits u[2];
+ union ieee_ext_u u[2];
- u[0].e = x;
- mask_nbit_l(u[0]);
- u[1].e = y;
- mask_nbit_l(u[1]);
+ memset(&u, 0, sizeof u);
+ u[0].extu_ld = x;
+ u[0].extu_ext.ext_frach &= ~0x80000000;
+ u[1].extu_ld = y;
+ u[1].extu_ext.ext_frach &= ~0x80000000;
/* Check for NaNs to avoid raising spurious exceptions. */
- if (u[0].bits.exp == 32767 && (u[0].bits.manh | u[0].bits.manl) != 0)
+ if (u[0].extu_ext.ext_exp == EXT_EXP_INFNAN &&
+ (u[0].extu_ext.ext_frach | u[0].extu_ext.ext_fracl) != 0)
return (y);
- if (u[1].bits.exp == 32767 && (u[1].bits.manh | u[1].bits.manl) != 0)
+ if (u[1].extu_ext.ext_exp == EXT_EXP_INFNAN &&
+ (u[1].extu_ext.ext_frach | u[1].extu_ext.ext_fracl) != 0)
return (x);
- /* Handle comparisons of signed zeroes. */
- if (u[0].bits.sign != u[1].bits.sign)
- return (u[0].bits.sign ? y : x);
+ /* Handle comparisons of ext_signed zeroes. */
+ if (u[0].extu_ext.ext_sign != u[1].extu_ext.ext_sign)
+ return (u[0].extu_ext.ext_sign ? y : x);
return (x > y ? x : y);
}
+#endif
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
+ *
+ * $NetBSD: s_fmin.c,v 1.1 2009/10/04 22:04:30 christos Exp $
*/
-#include <sys/cdefs.h>
-
#include <math.h>
-#include "fpmath.h"
+#include <machine/ieee.h>
double
fmin(double x, double y)
{
- union IEEEd2bits u[2];
+ union ieee_double_u u[2];
- u[0].d = x;
- u[1].d = y;
+ u[0].dblu_d = x;
+ u[1].dblu_d = y;
/* Check for NaNs to avoid raising spurious exceptions. */
- if (u[0].bits.exp == 2047 && (u[0].bits.manh | u[0].bits.manl) != 0)
+ if (u[0].dblu_dbl.dbl_exp == DBL_EXP_INFNAN &&
+ (u[0].dblu_dbl.dbl_frach | u[0].dblu_dbl.dbl_fracl) != 0)
return (y);
- if (u[1].bits.exp == 2047 && (u[1].bits.manh | u[1].bits.manl) != 0)
+ if (u[1].dblu_dbl.dbl_exp == DBL_EXP_INFNAN &&
+ (u[1].dblu_dbl.dbl_frach | u[1].dblu_dbl.dbl_fracl) != 0)
return (x);
/* Handle comparisons of signed zeroes. */
- if (u[0].bits.sign != u[1].bits.sign)
- return (u[u[1].bits.sign].d);
+ if (u[0].dblu_dbl.dbl_sign != u[1].dblu_dbl.dbl_sign)
+ return (u[u[1].dblu_dbl.dbl_sign].dblu_d);
return (x < y ? x : y);
}
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
+ *
+ * $NetBSD: s_fminf.c,v 1.2 2010/03/08 01:05:20 snj Exp $
*/
-#include <sys/cdefs.h>
-
#include <math.h>
-#include "fpmath.h"
+#include <machine/ieee.h>
float
fminf(float x, float y)
{
- union IEEEf2bits u[2];
+ union ieee_single_u u[2];
- u[0].f = x;
- u[1].f = y;
+ u[0].sngu_f = x;
+ u[1].sngu_f = y;
/* Check for NaNs to avoid raising spurious exceptions. */
- if (u[0].bits.exp == 255 && u[0].bits.man != 0)
+ if (u[0].sngu_sng.sng_exp == SNG_EXP_INFNAN &&
+ u[0].sngu_sng.sng_frac != 0)
return (y);
- if (u[1].bits.exp == 255 && u[1].bits.man != 0)
+ if (u[1].sngu_sng.sng_exp == SNG_EXP_INFNAN &&
+ u[1].sngu_sng.sng_frac != 0)
return (x);
- /* Handle comparisons of signed zeroes. */
- if (u[0].bits.sign != u[1].bits.sign)
- return (u[u[1].bits.sign].f);
+ /* Handle comparisons of sng_singed zeroes. */
+ if (u[0].sngu_sng.sng_sign != u[1].sngu_sng.sng_sign)
+ return (u[u[1].sngu_sng.sng_sign].sngu_f);
return (x < y ? x : y);
}
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
+ *
+ * $NetBSD: s_fminl.c,v 1.3 2011/07/04 11:46:41 mrg Exp $
*/
-#include <sys/cdefs.h>
-
#include <math.h>
+#include <string.h>
-#include "fpmath.h"
-
+#include <machine/ieee.h>
+#ifdef EXT_EXP_INFNAN
long double
fminl(long double x, long double y)
{
- union IEEEl2bits u[2];
+ union ieee_ext_u u[2];
- u[0].e = x;
- mask_nbit_l(u[0]);
- u[1].e = y;
- mask_nbit_l(u[1]);
+ memset(&u, 0, sizeof u);
+ u[0].extu_ld = x;
+ u[0].extu_ext.ext_frach &= ~0x80000000;
+ u[1].extu_ld = y;
+ u[1].extu_ext.ext_frach &= ~0x80000000;
/* Check for NaNs to avoid raising spurious exceptions. */
- if (u[0].bits.exp == 32767 && (u[0].bits.manh | u[0].bits.manl) != 0)
+ if (u[0].extu_ext.ext_exp == EXT_EXP_INFNAN &&
+ (u[0].extu_ext.ext_frach | u[0].extu_ext.ext_fracl) != 0)
return (y);
- if (u[1].bits.exp == 32767 && (u[1].bits.manh | u[1].bits.manl) != 0)
+ if (u[1].extu_ext.ext_exp == EXT_EXP_INFNAN &&
+ (u[1].extu_ext.ext_frach | u[1].extu_ext.ext_fracl) != 0)
return (x);
- /* Handle comparisons of signed zeroes. */
- if (u[0].bits.sign != u[1].bits.sign)
- return (u[1].bits.sign ? y : x);
+ /* Handle comparisons of ext_signed zeroes. */
+ if (u[0].extu_ext.ext_sign != u[1].extu_ext.ext_sign)
+ return (u[1].extu_ext.ext_sign ? y : x);
return (x < y ? x : y);
}
+#endif
--- /dev/null
+/* @(#)s_frexp.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ *
+ * $NetBSD: s_frexp.c,v 1.13 2008/09/28 18:54:55 christos Exp $
+ */
+
+/*
+ * for non-zero x
+ * x = frexp(arg,&exp);
+ * return a double fp quantity x such that 0.5 <= |x| <1.0
+ * and the corresponding binary exponent "exp". That is
+ * arg = x*2^exp.
+ * If arg is inf, 0.0, or NaN, then frexp(arg,&exp) returns arg
+ * with *exp=0.
+ */
+
+#include <math.h>
+#include "math_private.h"
+
+static const double
+two54 = 1.80143985094819840000e+16; /* 0x43500000, 0x00000000 */
+
+double
+frexp(double x, int *eptr)
+{
+ int32_t hx, ix, lx;
+ EXTRACT_WORDS(hx,lx,x);
+ ix = 0x7fffffff&hx;
+ *eptr = 0;
+ if(ix>=0x7ff00000||((ix|lx)==0)) return x; /* 0,inf,nan */
+ if (ix<0x00100000) { /* subnormal */
+ x *= two54;
+ GET_HIGH_WORD(hx,x);
+ ix = hx&0x7fffffff;
+ *eptr = -54;
+ }
+ *eptr += ((uint32_t)ix>>20)-1022;
+ hx = (hx&0x800fffff)|0x3fe00000;
+ SET_HIGH_WORD(x,hx);
+ return x;
+}
* is preserved.
* ====================================================
*
- * $NetBSD: s_frexpf.c,v 1.9 2002/12/05 16:03:42 scw Exp $
- * $DragonFly: src/lib/libm/src/s_frexpf.c,v 1.1 2005/07/26 21:15:20 joerg Exp $
+ * $NetBSD: s_frexpf.c,v 1.10 2007/08/21 20:12:27 drochner Exp $
*/
#include <math.h>
}
*eptr += (ix>>23)-126;
hx = (hx&0x807fffff)|0x3f000000;
- *(int*)(void*)&x = hx;
+ SET_FLOAT_WORD(x,hx);
return x;
}
* is preserved.
* ====================================================
*
- * $NetBSD: s_ilogb.c,v 1.12 2002/05/26 22:01:56 wiz Exp $
- * $DragonFly: src/lib/libm/src/s_ilogb.c,v 1.1 2005/07/26 21:15:20 joerg Exp $
+ * $NetBSD: s_ilogb.c,v 1.13 2011/07/28 22:32:29 joerg Exp $
*/
/* ilogb(double x)
+/* $NetBSD: s_ilogbl.c,v 1.1 2011/07/28 22:32:29 joerg Exp $ */
+
/*-
- * Copyright (c) 2006 The NetBSD Foundation, Inc.
+ * Copyright (c) 2011 The NetBSD Foundation, Inc.
* All rights reserved.
*
* This code is derived from software contributed to The NetBSD Foundation
- * by Klaus Klein.
+ * by Joerg Sonnenberger.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
- * 3. All advertising materials mentioning features or use of this software
- * must display the following acknowledgement:
- * This product includes software developed by the NetBSD
- * Foundation, Inc. and its contributors.
- * 4. Neither the name of The NetBSD Foundation nor the names of its
- * contributors may be used to endorse or promote products derived
- * from this software without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
* ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
* POSSIBILITY OF SUCH DAMAGE.
- *
- * $NetBSD: nanl.c,v 1.1 2006/03/15 22:07:09 kleink Exp $
- * $DragonFly: src/lib/libm/gen/nanl.c,v 1.1 2007/06/17 17:46:01 pavalos Exp $
*/
-#define NAN_FUNCTION nanl
-#define NAN_TYPE long double
-#define NAN_STRTOD strtod /* XXX should be strtold */
+#include <float.h>
+#include <math.h>
+#include <machine/ieee.h>
+
+#ifdef __HAVE_LONG_DOUBLE
+
+#if LDBL_MANT_DIG == 64
+#define FROM_UNDERFLOW 0x1p65L
+#elif LDBL_MANT_DIG == 113
+#define FROM_UNDERFLOW 0x1p114L
+#else
+#error Unsupported long double format
+#endif
+
+int
+ilogbl(long double x)
+{
+ union ieee_ext_u u;
+
+ if (x == 0.0L)
+ return 0x80000001; /* ilogbl(0) = 0x80000001 */
+
+ u.extu_ld = x;
+
+ if (u.extu_ext.ext_exp == EXT_EXP_INFNAN)
+ return 0x7fffffff;
+
+ if (u.extu_ext.ext_exp == 0) {
+ /*
+ * Scale denormalized numbers slightly,
+ * so that they are normal.
+ */
+ u.extu_ld *= FROM_UNDERFLOW;
+ return u.extu_ext.ext_exp - EXT_EXP_BIAS - LDBL_MANT_DIG - 1;
+ }
+ return u.extu_ext.ext_exp - EXT_EXP_BIAS;
+
+}
-#include "nan.c"
+#endif
-/* s_ldexpf.c -- float version of s_ldexp.c.
+/* s_ldexp0f.c -- float version of s_ldexp0.c.
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
*/
* is preserved.
* ====================================================
*
- * $NetBSD: s_ldexpf.c,v 1.6 2002/05/26 22:01:57 wiz Exp $
- * $DragonFly: src/lib/libm/src/s_ldexpf.c,v 1.1 2005/07/26 21:15:20 joerg Exp $
+ * $NetBSD: s_ldexpf.c,v 1.8 2010/04/23 19:17:07 drochner Exp $
*/
#include <math.h>
#include <errno.h>
float
-ldexpf(float value, int expo)
+ldexpf(float value, int exp0)
{
if(!finitef(value)||value==(float)0.0) return value;
- value = scalbnf(value,expo);
+ value = scalbnf(value,exp0);
if(!finitef(value)||value==(float)0.0) errno = ERANGE;
return value;
}
* is preserved.
* ====================================================
*
- * $NetBSD: s_logb.c,v 1.11 2002/05/26 22:01:57 wiz Exp $
- * $DragonFly: src/lib/libm/src/s_logb.c,v 1.1 2005/07/26 21:15:20 joerg Exp $
+ * $NetBSD: s_logb.c,v 1.12 2011/08/03 14:13:07 joerg Exp $
*/
/*
+/* $NetBSD: s_logbl.c,v 1.1 2011/08/03 14:13:07 joerg Exp $ */
+
/*-
- * Copyright (c) 2006 The NetBSD Foundation, Inc.
+ * Copyright (c) 2011 The NetBSD Foundation, Inc.
* All rights reserved.
*
* This code is derived from software contributed to The NetBSD Foundation
- * by Klaus Klein.
+ * by Joerg Sonnenberger.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
- * 3. All advertising materials mentioning features or use of this software
- * must display the following acknowledgement:
- * This product includes software developed by the NetBSD
- * Foundation, Inc. and its contributors.
- * 4. Neither the name of The NetBSD Foundation nor the names of its
- * contributors may be used to endorse or promote products derived
- * from this software without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
* ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
* POSSIBILITY OF SUCH DAMAGE.
- *
- * $NetBSD: nanf.c,v 1.1 2006/03/15 22:07:09 kleink Exp $
- * $DragonFly: src/lib/libm/gen/nanf.c,v 1.1 2007/06/17 17:46:01 pavalos Exp $
*/
-#define NAN_FUNCTION nanf
-#define NAN_TYPE float
-#define NAN_STRTOD strtod /* XXX should be strtof */
+#include <float.h>
+#include <math.h>
+#include <machine/ieee.h>
+
+#ifdef __HAVE_LONG_DOUBLE
+
+#if LDBL_MANT_DIG == 64
+#define FROM_UNDERFLOW 0x1p65L
+#elif LDBL_MANT_DIG == 113
+#define FROM_UNDERFLOW 0x1p114L
+#else
+#error Unsupported long double format
+#endif
+
+long double
+logbl(long double x)
+{
+ union ieee_ext_u u;
+
+ if (x == 0.0L)
+ return -1.0L / fabsl(x); /* -HUGE_VALL + exception */
+
+ u.extu_ld = x;
+
+ if (u.extu_ext.ext_exp == EXT_EXP_INFNAN)
+ return fabsl(x); /* NaN or +Inf */
+
+ if (u.extu_ext.ext_exp == 0) {
+ /*
+ * Scale denormalized numbers slightly,
+ * so that they are normal.
+ */
+ u.extu_ld *= FROM_UNDERFLOW;
+ return u.extu_ext.ext_exp - EXT_EXP_BIAS - LDBL_MANT_DIG - 1;
+ }
+ return u.extu_ext.ext_exp - EXT_EXP_BIAS;
+
+}
-#include "nan.c"
+#endif
* is preserved.
* ====================================================
*
- * $NetBSD: s_modff.c,v 1.7 2002/05/26 22:01:57 wiz Exp $
- * $DragonFly: src/lib/libm/src/s_modff.c,v 1.1 2005/07/26 21:15:20 joerg Exp $
+ * $NetBSD: s_modff.c,v 1.9 2010/01/27 14:07:41 drochner Exp $ $
*/
#include <math.h>
float
modff(float x, float *iptr)
{
- int32_t i0,j0_;
+ int32_t i0,jj0;
u_int32_t i;
GET_FLOAT_WORD(i0,x);
- j0_ = ((i0>>23)&0xff)-0x7f; /* exponent of x */
- if(j0_<23) { /* integer part in x */
- if(j0_<0) { /* |x|<1 */
+ jj0 = ((i0>>23)&0xff)-0x7f; /* exponent of x */
+ if(jj0<23) { /* integer part in x */
+ if(jj0<0) { /* |x|<1 */
SET_FLOAT_WORD(*iptr,i0&0x80000000); /* *iptr = +-0 */
return x;
} else {
- i = (0x007fffff)>>j0_;
+ i = (0x007fffff)>>jj0;
if((i0&i)==0) { /* x is integral */
u_int32_t ix;
*iptr = x;
} else { /* no fraction part */
u_int32_t ix;
*iptr = x*one;
+ if (jj0 == 0x80) /* +-inf or NaN */
+ return 0.0 / x; /* +-0 or NaN */
GET_FLOAT_WORD(ix,x);
SET_FLOAT_WORD(x,ix&0x80000000); /* return +-0 */
return x;
* is preserved.
* ====================================================
*
- * $NetBSD: s_nextafter.c,v 1.11 2002/05/26 22:01:57 wiz Exp $
- * $DragonFly: src/lib/libm/src/s_nextafter.c,v 1.1 2005/07/26 21:15:20 joerg Exp $
+ * $NetBSD: s_nextafter.c,v 1.12 2011/04/18 15:59:09 drochner Exp $
*/
/* IEEE functions
if(((ix>=0x7ff00000)&&((ix-0x7ff00000)|lx)!=0) || /* x is nan */
((iy>=0x7ff00000)&&((iy-0x7ff00000)|ly)!=0)) /* y is nan */
return x+y;
- if(x==y) return x; /* x=y, return x */
+ if(x==y) return y; /* x=y, return y */
if((ix|lx)==0) { /* x == 0 */
INSERT_WORDS(x,hy&0x80000000,1); /* return +-minsubnormal */
y = x*x;
* is preserved.
* ====================================================
*
- * $NetBSD: s_nextafterf.c,v 1.7 2002/05/26 22:01:58 wiz Exp $
- * $DragonFly: src/lib/libm/src/s_nextafterf.c,v 1.1 2005/07/26 21:15:20 joerg Exp $
+ * $NetBSD: s_nextafterf.c,v 1.8 2011/04/18 15:59:09 drochner Exp $
*/
#include <math.h>
if((ix>0x7f800000) || /* x is nan */
(iy>0x7f800000)) /* y is nan */
return x+y;
- if(x==y) return x; /* x=y, return x */
+ if(x==y) return y; /* x=y, return y */
if(ix==0) { /* x == 0 */
SET_FLOAT_WORD(x,(hy&0x80000000)|1);/* return +-minsubnormal */
y = x*x;
--- /dev/null
+/* $NetBSD: s_nextafterl.c,v 1.2 2010/09/17 20:39:39 christos Exp $ */
+
+/* @(#)s_nextafter.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include <float.h>
+#include <math.h>
+#include <machine/ieee.h>
+
+#ifdef EXT_EXP_INFNAN
+#if LDBL_MAX_EXP != 0x4000
+#error "Unsupported long double format"
+#endif
+
+/*
+ * IEEE functions
+ * nextafterl(x,y)
+ * return the next machine floating-point number of x in the
+ * direction toward y.
+ * Special cases:
+ * If x == y, y shall be returned
+ * If x or y is NaN, a NaN shall be returned
+ */
+long double
+nextafterl(long double x, long double y)
+{
+ volatile long double t;
+ union ieee_ext_u ux, uy;
+
+ ux.extu_ld = x;
+ uy.extu_ld = y;
+
+ if ((ux.extu_exp == EXT_EXP_NAN &&
+ ((ux.extu_frach &~ LDBL_NBIT)|ux.extu_fracl) != 0) ||
+ (uy.extu_exp == EXT_EXP_NAN &&
+ ((uy.extu_frach &~ LDBL_NBIT)|uy.extu_fracl) != 0))
+ return x+y; /* x or y is nan */
+
+ if (x == y) return y; /* x=y, return y */
+
+ if (x == 0.0) {
+ ux.extu_frach = 0; /* return +-minsubnormal */
+ ux.extu_fracl = 1;
+ ux.extu_sign = uy.extu_sign;
+ t = ux.extu_ld * ux.extu_ld;
+ if (t == ux.extu_ld)
+ return t;
+ else
+ return ux.extu_ld; /* raise underflow flag */
+ }
+
+ if ((x>0.0) ^ (x<y)) { /* x -= ulp */
+ if (ux.extu_fracl == 0) {
+ if ((ux.extu_frach & ~LDBL_NBIT) == 0)
+ ux.extu_exp -= 1;
+ ux.extu_frach = (ux.extu_frach - 1) |
+ (ux.extu_frach & LDBL_NBIT);
+ }
+ ux.extu_fracl -= 1;
+ } else { /* x += ulp */
+ ux.extu_fracl += 1;
+ if (ux.extu_fracl == 0) {
+ ux.extu_frach = (ux.extu_frach + 1) |
+ (ux.extu_frach & LDBL_NBIT);
+ if ((ux.extu_frach & ~LDBL_NBIT) == 0)
+ ux.extu_exp += 1;
+ }
+ }
+
+ if (ux.extu_exp == EXT_EXP_INF)
+ return x+x; /* overflow */
+
+ if (ux.extu_exp == 0) { /* underflow */
+ mask_nbit_l(ux);
+ t = ux.extu_ld * ux.extu_ld;
+ if (t != ux.extu_ld) /* raise underflow flag */
+ return ux.extu_ld;
+ }
+
+ return ux.extu_ld;
+}
+#endif
--- /dev/null
+/* $NetBSD: s_nexttoward.c,v 1.1 2010/09/15 16:12:05 christos Exp $ */
+
+/* @(#)s_nextafter.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+/*
+ * We assume that a long double has a 15-bit exponent. On systems
+ * where long double is the same as double, nexttoward() is an alias
+ * for nextafter(), so we don't use this routine.
+ */
+#include <float.h>
+
+#include <machine/ieee.h>
+#include <math.h>
+#include "math_private.h"
+
+#if LDBL_MAX_EXP != 0x4000
+#error "Unsupported long double format"
+#endif
+
+/*
+ * The nexttoward() function is equivalent to nextafter() function,
+ * except that the second parameter shall have type long double and
+ * the functions shall return y converted to the type of the function
+ * if x equals y.
+ *
+ * Special cases: XXX
+ */
+double
+nexttoward(double x, long double y)
+{
+ union ieee_ext_u uy;
+ volatile double t;
+ int32_t hx, ix;
+ uint32_t lx;
+
+ EXTRACT_WORDS(hx, lx, x);
+ ix = hx & 0x7fffffff; /* |x| */
+ uy.extu_ld = y;
+
+ if (((ix >= 0x7ff00000) && ((ix - 0x7ff00000) | lx) != 0) ||
+ (uy.extu_exp == 0x7fff &&
+ ((uy.extu_frach & ~LDBL_NBIT) | uy.extu_fracl) != 0))
+ return x+y; /* x or y is nan */
+
+ if (x == y)
+ return (double)y; /* x=y, return y */
+
+ if (x == 0.0) {
+ INSERT_WORDS(x, uy.extu_sign<<31, 1); /* return +-minsubnormal */
+ t = x*x;
+ if (t == x)
+ return t;
+ else
+ return x; /* raise underflow flag */
+ }
+
+ if ((hx > 0.0) ^ (x < y)) { /* x -= ulp */
+ if (lx == 0) hx -= 1;
+ lx -= 1;
+ } else { /* x += ulp */
+ lx += 1;
+ if (lx == 0) hx += 1;
+ }
+ ix = hx & 0x7ff00000;
+ if (ix >= 0x7ff00000) return x+x; /* overflow */
+ if (ix < 0x00100000) { /* underflow */
+ t = x*x;
+ if (t != x) { /* raise underflow flag */
+ INSERT_WORDS(y, hx, lx);
+ return y;
+ }
+ }
+ INSERT_WORDS(x, hx, lx);
+
+ return x;
+}
--- /dev/null
+/* @(#)e_fmod.c 1.3 95/01/18 */
+/*-
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+
+#include <float.h>
+
+#include <math.h>
+#include "math_private.h"
+
+static const double Zero[] = {0.0, -0.0,};
+
+/*
+ * Return the IEEE remainder and set *quo to the last n bits of the
+ * quotient, rounded to the nearest integer. We choose n=31 because
+ * we wind up computing all the integer bits of the quotient anyway as
+ * a side-effect of computing the remainder by the shift and subtract
+ * method. In practice, this is far more bits than are needed to use
+ * remquo in reduction algorithms.
+ */
+double
+remquo(double x, double y, int *quo)
+{
+ int32_t n,hx,hy,hz,ix,iy,sx,i;
+ u_int32_t lx,ly,lz,q,sxy;
+
+ EXTRACT_WORDS(hx,lx,x);
+ EXTRACT_WORDS(hy,ly,y);
+ sxy = (hx ^ hy) & 0x80000000;
+ sx = hx&0x80000000; /* sign of x */
+ hx ^=sx; /* |x| */
+ hy &= 0x7fffffff; /* |y| */
+
+ /* purge off exception values */
+ if((hy|ly)==0||(hx>=0x7ff00000)|| /* y=0,or x not finite */
+ ((hy|((ly|-ly)>>31))>0x7ff00000)) /* or y is NaN */
+ return (x*y)/(x*y);
+ if(hx<=hy) {
+ if((hx<hy)||(lx<ly)) {
+ q = 0;
+ goto fixup; /* |x|<|y| return x or x-y */
+ }
+ if(lx==ly) {
+ *quo = 1;
+ return Zero[(u_int32_t)sx>>31]; /* |x|=|y| return x*0*/
+ }
+ }
+
+ /* determine ix = ilogb(x) */
+ if(hx<0x00100000) { /* subnormal x */
+ if(hx==0) {
+ for (ix = -1043, i=lx; i>0; i<<=1) ix -=1;
+ } else {
+ for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1;
+ }
+ } else ix = (hx>>20)-1023;
+
+ /* determine iy = ilogb(y) */
+ if(hy<0x00100000) { /* subnormal y */
+ if(hy==0) {
+ for (iy = -1043, i=ly; i>0; i<<=1) iy -=1;
+ } else {
+ for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1;
+ }
+ } else iy = (hy>>20)-1023;
+
+ /* set up {hx,lx}, {hy,ly} and align y to x */
+ if(ix >= -1022)
+ hx = 0x00100000|(0x000fffff&hx);
+ else { /* subnormal x, shift x to normal */
+ n = -1022-ix;
+ if(n<=31) {
+ hx = (hx<<n)|(lx>>(32-n));
+ lx <<= n;
+ } else {
+ hx = lx<<(n-32);
+ lx = 0;
+ }
+ }
+ if(iy >= -1022)
+ hy = 0x00100000|(0x000fffff&hy);
+ else { /* subnormal y, shift y to normal */
+ n = -1022-iy;
+ if(n<=31) {
+ hy = (hy<<n)|(ly>>(32-n));
+ ly <<= n;
+ } else {
+ hy = ly<<(n-32);
+ ly = 0;
+ }
+ }
+
+ /* fix point fmod */
+ n = ix - iy;
+ q = 0;
+ while(n--) {
+ hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
+ if(hz<0){hx = hx+hx+(lx>>31); lx = lx+lx;}
+ else {hx = hz+hz+(lz>>31); lx = lz+lz; q++;}
+ q <<= 1;
+ }
+ hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
+ if(hz>=0) {hx=hz;lx=lz;q++;}
+
+ /* convert back to floating value and restore the sign */
+ if((hx|lx)==0) { /* return sign(x)*0 */
+ *quo = (sxy ? -q : q);
+ return Zero[(u_int32_t)sx>>31];
+ }
+ while(hx<0x00100000) { /* normalize x */
+ hx = hx+hx+(lx>>31); lx = lx+lx;
+ iy -= 1;
+ }
+ if(iy>= -1022) { /* normalize output */
+ hx = ((hx-0x00100000)|((iy+1023)<<20));
+ } else { /* subnormal output */
+ n = -1022 - iy;
+ if(n<=20) {
+ lx = (lx>>n)|((u_int32_t)hx<<(32-n));
+ hx >>= n;
+ } else if (n<=31) {
+ lx = (hx<<(32-n))|(lx>>n); hx = sx;
+ } else {
+ lx = hx>>(n-32); hx = sx;
+ }
+ }
+fixup:
+ INSERT_WORDS(x,hx,lx);
+ y = fabs(y);
+ if (y < 0x1p-1021) {
+ if (x+x>y || (x+x==y && (q & 1))) {
+ q++;
+ x-=y;
+ }
+ } else if (x>0.5*y || (x==0.5*y && (q & 1))) {
+ q++;
+ x-=y;
+ }
+ GET_HIGH_WORD(hx,x);
+ SET_HIGH_WORD(x,hx^sx);
+ q &= 0x7fffffff;
+ *quo = (sxy ? -q : q);
+ return x;
+}
--- /dev/null
+/* @(#)e_fmod.c 1.3 95/01/18 */
+/*-
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+
+#include <math.h>
+#include "math_private.h"
+
+static const float Zero[] = {0.0, -0.0,};
+
+/*
+ * Return the IEEE remainder and set *quo to the last n bits of the
+ * quotient, rounded to the nearest integer. We choose n=31 because
+ * we wind up computing all the integer bits of the quotient anyway as
+ * a side-effect of computing the remainder by the shift and subtract
+ * method. In practice, this is far more bits than are needed to use
+ * remquo in reduction algorithms.
+ */
+float
+remquof(float x, float y, int *quo)
+{
+ int32_t n,hx,hy,hz,ix,iy,sx,i;
+ u_int32_t q,sxy;
+
+ GET_FLOAT_WORD(hx,x);
+ GET_FLOAT_WORD(hy,y);
+ sxy = (hx ^ hy) & 0x80000000;
+ sx = hx&0x80000000; /* sign of x */
+ hx ^=sx; /* |x| */
+ hy &= 0x7fffffff; /* |y| */
+
+ /* purge off exception values */
+ if(hy==0||hx>=0x7f800000||hy>0x7f800000) /* y=0,NaN;or x not finite */
+ return (x*y)/(x*y);
+ if(hx<hy) {
+ q = 0;
+ goto fixup; /* |x|<|y| return x or x-y */
+ } else if(hx==hy) {
+ *quo = 1;
+ return Zero[(u_int32_t)sx>>31]; /* |x|=|y| return x*0*/
+ }
+
+ /* determine ix = ilogb(x) */
+ if(hx<0x00800000) { /* subnormal x */
+ for (ix = -126,i=(hx<<8); i>0; i<<=1) ix -=1;
+ } else ix = (hx>>23)-127;
+
+ /* determine iy = ilogb(y) */
+ if(hy<0x00800000) { /* subnormal y */
+ for (iy = -126,i=(hy<<8); i>0; i<<=1) iy -=1;
+ } else iy = (hy>>23)-127;
+
+ /* set up {hx,lx}, {hy,ly} and align y to x */
+ if(ix >= -126)
+ hx = 0x00800000|(0x007fffff&hx);
+ else { /* subnormal x, shift x to normal */
+ n = -126-ix;
+ hx <<= n;
+ }
+ if(iy >= -126)
+ hy = 0x00800000|(0x007fffff&hy);
+ else { /* subnormal y, shift y to normal */
+ n = -126-iy;
+ hy <<= n;
+ }
+
+ /* fix point fmod */
+ n = ix - iy;
+ q = 0;
+ while(n--) {
+ hz=hx-hy;
+ if(hz<0) hx = hx << 1;
+ else {hx = hz << 1; q++;}
+ q <<= 1;
+ }
+ hz=hx-hy;
+ if(hz>=0) {hx=hz;q++;}
+
+ /* convert back to floating value and restore the sign */
+ if(hx==0) { /* return sign(x)*0 */
+ *quo = (sxy ? -q : q);
+ return Zero[(u_int32_t)sx>>31];
+ }
+ while(hx<0x00800000) { /* normalize x */
+ hx <<= 1;
+ iy -= 1;
+ }
+ if(iy>= -126) { /* normalize output */
+ hx = ((hx-0x00800000)|((iy+127)<<23));
+ } else { /* subnormal output */
+ n = -126 - iy;
+ hx >>= n;
+ }
+fixup:
+ SET_FLOAT_WORD(x,hx);
+ y = fabsf(y);
+ if (y < 0x1p-125f) {
+ if (x+x>y || (x+x==y && (q & 1))) {
+ q++;
+ x-=y;
+ }
+ } else if (x>0.5f*y || (x==0.5f*y && (q & 1))) {
+ q++;
+ x-=y;
+ }
+ GET_FLOAT_WORD(hx,x);
+ SET_FLOAT_WORD(x,hx^sx);
+ q &= 0x7fffffff;
+ *quo = (sxy ? -q : q);
+ return x;
+}
* is preserved.
* ====================================================
*
- * $NetBSD: s_rint.c,v 1.11 2002/05/26 22:01:58 wiz Exp $
- * $DragonFly: src/lib/libm/src/s_rint.c,v 1.1 2005/07/26 21:15:20 joerg Exp $
+ * $NetBSD: s_rint.c,v 1.12 2008/04/25 22:21:53 christos Exp $
*/
/*
double
rint(double x)
{
- int32_t i0,j0,sx;
+ int32_t i0,jj0,sx;
u_int32_t i,i1;
double w,t;
EXTRACT_WORDS(i0,i1,x);
sx = (i0>>31)&1;
- j0 = ((i0>>20)&0x7ff)-0x3ff;
- if(j0<20) {
- if(j0<0) {
+ jj0 = ((i0>>20)&0x7ff)-0x3ff;
+ if(jj0<20) {
+ if(jj0<0) {
if(((i0&0x7fffffff)|i1)==0) return x;
i1 |= (i0&0x0fffff);
i0 &= 0xfffe0000;
SET_HIGH_WORD(t,(i0&0x7fffffff)|(sx<<31));
return t;
} else {
- i = (0x000fffff)>>j0;
+ i = (0x000fffff)>>jj0;
if(((i0&i)|i1)==0) return x; /* x is integral */
i>>=1;
if(((i0&i)|i1)!=0) {
- if(j0==19) i1 = 0x40000000; else
- i0 = (i0&(~i))|((0x20000)>>j0);
+ if(jj0==19) i1 = 0x40000000; else
+ i0 = (i0&(~i))|((0x20000)>>jj0);
}
}
- } else if (j0>51) {
- if(j0==0x400) return x+x; /* inf or NaN */
+ } else if (jj0>51) {
+ if(jj0==0x400) return x+x; /* inf or NaN */
else return x; /* x is integral */
} else {
- i = ((u_int32_t)(0xffffffff))>>(j0-20);
+ i = ((u_int32_t)(0xffffffff))>>(jj0-20);
if((i1&i)==0) return x; /* x is integral */
i>>=1;
- if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20));
+ if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(jj0-20));
}
INSERT_WORDS(x,i0,i1);
w = TWO52[sx]+x;
* is preserved.
* ====================================================
*
- * $NetBSD: s_rintf.c,v 1.8 2006/08/01 20:14:35 drochner Exp $
- * $DragonFly: src/lib/libm/src/s_rintf.c,v 1.2 2007/06/17 02:27:53 pavalos Exp $
+ * $NetBSD: s_rintf.c,v 1.9 2008/04/25 22:21:53 christos Exp $
*/
#include <math.h>
float
rintf(float x)
{
- int32_t i0,j0,sx;
+ int32_t i0,jj0,sx;
u_int32_t i,i1;
#ifdef __i386__ /* XXX gcc4 will omit the rounding otherwise */
volatile
float t;
GET_FLOAT_WORD(i0,x);
sx = (i0>>31)&1;
- j0 = ((i0>>23)&0xff)-0x7f;
- if(j0<23) {
- if(j0<0) {
+ jj0 = ((i0>>23)&0xff)-0x7f;
+ if(jj0<23) {
+ if(jj0<0) {
if((i0&0x7fffffff)==0) return x;
i1 = (i0&0x07fffff);
i0 &= 0xfff00000;
SET_FLOAT_WORD(t,(i0&0x7fffffff)|(sx<<31));
return t;
} else {
- i = (0x007fffff)>>j0;
+ i = (0x007fffff)>>jj0;
if((i0&i)==0) return x; /* x is integral */
i>>=1;
- if((i0&i)!=0) i0 = (i0&(~i))|((0x100000)>>j0);
+ if((i0&i)!=0) i0 = (i0&(~i))|((0x100000)>>jj0);
}
} else {
- if(j0==0x80) return x+x; /* inf or NaN */
+ if(jj0==0x80) return x+x; /* inf or NaN */
else return x; /* x is integral */
}
SET_FLOAT_WORD(x,i0);
* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*
- * $FreeBSD: src/lib/msun/src/s_round.c,v 1.1 2004/06/07 08:05:36 das Exp $
- * $NetBSD: s_round.c,v 1.1 2004/07/10 13:49:10 junyoung Exp $
- * $DragonFly: src/lib/libm/src/s_round.c,v 1.1 2005/07/26 21:15:20 joerg Exp $
+ * $NetBSD: s_round.c,v 1.2 2007/08/21 20:10:27 drochner Exp $
*/
#include <math.h>
return (x);
if (x >= 0.0) {
- t = ceil(x);
- if (t - x > 0.5)
- t -= 1.0;
+ t = floor(x);
+ if (x - t >= 0.5)
+ t += 1.0;
return (t);
} else {
- t = ceil(-x);
- if (t + x > 0.5)
- t -= 1.0;
+ x = -x;
+ t = floor(x);
+ if (x - t >= 0.5)
+ t += 1.0;
return (-t);
}
}
* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*
- * $FreeBSD: src/lib/msun/src/s_roundf.c,v 1.1 2004/06/07 08:05:36 das Exp $
- * $NetBSD: s_roundf.c,v 1.1 2004/07/10 13:49:10 junyoung Exp $
- * $DragonFly: src/lib/libm/src/s_roundf.c,v 1.1 2005/07/26 21:15:20 joerg Exp $
+ * $NetBSD: s_roundf.c,v 1.2 2007/08/21 20:10:27 drochner Exp $
*/
#include <math.h>
return (x);
if (x >= 0.0) {
- t = ceilf(x);
- if (t - x > 0.5)
- t -= 1.0;
+ t = floorf(x);
+ if (x - t >= 0.5)
+ t += 1.0;
return (t);
} else {
- t = ceilf(-x);
- if (t + x > 0.5)
- t -= 1.0;
+ x = -x;
+ t = floorf(x);
+ if (x - t >= 0.5)
+ t += 1.0;
return (-t);
}
}
* is preserved.
* ====================================================
*
- * $NetBSD: s_scalbn.c,v 1.12 2002/05/26 22:01:58 wiz Exp $
- * $DragonFly: src/lib/libm/src/s_scalbn.c,v 1.1 2005/07/26 21:15:20 joerg Exp $
+ * $NetBSD: s_scalbn.c,v 1.15 2011/07/26 16:10:16 joerg Exp $
*/
/*
{
int32_t k,hx,lx;
EXTRACT_WORDS(hx,lx,x);
- k = (hx&0x7ff00000)>>20; /* extract exponent */
+ k = ((uint32_t)hx&0x7ff00000)>>20; /* extract exponent */
if (k==0) { /* 0 or subnormal x */
if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */
x *= two54;
GET_HIGH_WORD(hx,x);
- k = ((hx&0x7ff00000)>>20) - 54;
+ k = (((uint32_t)hx&0x7ff00000)>>20) - 54;
if (n< -50000) return tiny*x; /*underflow*/
}
if (k==0x7ff) return x+x; /* NaN or Inf */
* is preserved.
* ====================================================
*
- * $NetBSD: s_scalbnf.c,v 1.8 2002/05/26 22:01:58 wiz Exp $
- * $DragonFly: src/lib/libm/src/s_scalbnf.c,v 1.1 2005/07/26 21:15:20 joerg Exp $
+ * $NetBSD: s_scalbnf.c,v 1.9 2010/04/23 19:17:07 drochner Exp $
*/
#include <math.h>
--- /dev/null
+/* $NetBSD: s_scalbnl.c,v 1.1 2011/07/26 16:10:16 joerg Exp $ */
+
+/*-
+ * Copyright (c) 2011 The NetBSD Foundation, Inc.
+ * All rights reserved.
+ *
+ * This code is derived from software contributed to The NetBSD Foundation
+ * by Joerg Sonnenberger.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions and the following disclaimer in the
+ * documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
+ * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+ * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+ * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR CONTRIBUTORS
+ * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#include <float.h>
+#include <math.h>
+#include <machine/ieee.h>
+
+#ifdef __HAVE_LONG_DOUBLE
+
+#if LDBL_MANT_DIG == 64
+#define FROM_UNDERFLOW 0x1p65L
+#define TO_UNDERFLOW 0x1p-65L
+#elif LDBL_MANT_DIG == 113
+#define FROM_UNDERFLOW 0x1p114L
+#define TO_UNDERFLOW 0x1p-114L
+#else
+#error Unsupported long double format
+#endif
+
+long double
+scalbnl(long double x, int n)
+{
+ union ieee_ext_u u;
+
+ /* Trivial cases first */
+ if (n == 0 || x == 0.0L)
+ return x;
+
+ u.extu_ld = x;
+
+ /* NaN and infinite don't change either, but trigger exception */
+ if (u.extu_ext.ext_exp == EXT_EXP_INFNAN)
+ return x + x;
+
+ /* Protect against integer overflow in calculation of new exponent */
+ if (n > LDBL_MAX_EXP - LDBL_MIN_EXP + LDBL_MANT_DIG)
+ goto overflow;
+ if (n < LDBL_MAX_EXP - LDBL_MIN_EXP + LDBL_MANT_DIG)
+ goto underflow;
+
+ /* Scale denormalized numbers slightly, so that they are normal */
+ if (u.extu_ext.ext_exp == 0) {
+ u.extu_ld *= FROM_UNDERFLOW;
+ n -= LDBL_MANT_DIG + 1;
+ }
+
+ n += u.extu_ext.ext_exp;
+ if (n >= LDBL_MAX_EXP + EXT_EXP_BIAS)
+ goto overflow;
+ /* Positive exponent (incl. bias) means normal result */
+ if (n > 0) {
+ u.extu_ext.ext_exp = n;
+ return u.extu_ld;
+ }
+ /* Shift the exponent and let the multiply below handle subnormal */
+ n += LDBL_MANT_DIG + 1;
+ if (n <= 0)
+ goto underflow;
+ u.extu_ext.ext_exp = n;
+ return u.extu_ld * TO_UNDERFLOW;
+
+underflow:
+ return LDBL_MIN * copysignl(LDBL_MIN, x);
+
+overflow:
+ return LDBL_MAX * copysignl(LDBL_MAX, x);
+}
+
+#endif
* is preserved.
* ====================================================
*
- * $NetBSD: s_sin.c,v 1.10 2002/05/26 22:01:58 wiz Exp $
- * $DragonFly: src/lib/libm/src/s_sin.c,v 1.1 2005/07/26 21:15:20 joerg Exp $
+ * $NetBSD: s_sin.c,v 1.11 2007/08/20 16:01:39 drochner Exp $
*/
/* sin(x)
* is preserved.
* ====================================================
*
- * $NetBSD: s_sinf.c,v 1.7 2002/05/26 22:01:58 wiz Exp $
- * $DragonFly: src/lib/libm/src/s_sinf.c,v 1.1 2005/07/26 21:15:20 joerg Exp $
+ * $NetBSD: s_sinf.c,v 1.8 2007/08/20 16:01:39 drochner Exp $
*/
#include <math.h>
* is preserved.
* ====================================================
*
- * $NetBSD: s_trunc.c,v 1.2 2007/01/17 23:24:22 hubertf Exp $
- * $DragonFly: src/lib/libm/src/s_trunc.c,v 1.1 2007/06/17 06:26:18 pavalos Exp $
+ * $NetBSD: s_trunc.c,v 1.3 2008/04/25 22:21:53 christos Exp $
*/
/*
double
trunc(double x)
{
- int32_t i0,i1,j00;
+ int32_t i0,i1,jj0;
uint32_t i;
EXTRACT_WORDS(i0,i1,x);
- j00 = ((i0>>20)&0x7ff)-0x3ff;
- if(j00<20) {
- if(j00<0) { /* raise inexact if x != 0 */
+ jj0 = ((i0>>20)&0x7ff)-0x3ff;
+ if(jj0<20) {
+ if(jj0<0) { /* raise inexact if x != 0 */
if(huge+x>0.0) {/* |x|<1, so return 0*sign(x) */
i0 &= 0x80000000U;
i1 = 0;
}
} else {
- i = (0x000fffff)>>j00;
+ i = (0x000fffff)>>jj0;
if(((i0&i)|i1)==0) return x; /* x is integral */
if(huge+x>0.0) { /* raise inexact flag */
i0 &= (~i); i1=0;
}
}
- } else if (j00>51) {
- if(j00==0x400) return x+x; /* inf or NaN */
+ } else if (jj0>51) {
+ if(jj0==0x400) return x+x; /* inf or NaN */
else return x; /* x is integral */
} else {
- i = ((u_int32_t)(0xffffffff))>>(j00-20);
+ i = ((u_int32_t)(0xffffffff))>>(jj0-20);
if((i1&i)==0) return x; /* x is integral */
if(huge+x>0.0) /* raise inexact flag */
i1 &= (~i);
* is preserved.
* ====================================================
*
- * $NetBSD: s_truncf.c,v 1.3 2007/01/17 23:24:22 hubertf Exp $
- * $DragonFly: src/lib/libm/src/s_truncf.c,v 1.1 2007/06/17 06:26:18 pavalos Exp $
+ * $NetBSD: s_truncf.c,v 1.4 2008/04/25 22:21:53 christos Exp $
*/
/*
float
truncf(float x)
{
- int32_t i0,j00;
+ int32_t i0,jj0;
uint32_t i;
GET_FLOAT_WORD(i0,x);
- j00 = ((i0>>23)&0xff)-0x7f;
- if(j00<23) {
- if(j00<0) { /* raise inexact if x != 0 */
+ jj0 = ((i0>>23)&0xff)-0x7f;
+ if(jj0<23) {
+ if(jj0<0) { /* raise inexact if x != 0 */
if(huge+x>0.0F) /* |x|<1, so return 0*sign(x) */
i0 &= 0x80000000;
} else {
- i = (0x007fffff)>>j00;
+ i = (0x007fffff)>>jj0;
if((i0&i)==0) return x; /* x is integral */
if(huge+x>0.0F) /* raise inexact flag */
i0 &= (~i);
}
} else {
- if(j00==0x80) return x+x; /* inf or NaN */
+ if(jj0==0x80) return x+x; /* inf or NaN */
else return x; /* x is integral */
}
SET_FLOAT_WORD(x,i0);
* A NaN is a `signalling NaN' if its QUIETNAN bit is clear in its
* high fraction; if the bit is set, it is a `quiet NaN'.
*/
-#define EXT_EXP_INFNAN 32767
+#define EXT_EXP_INFNAN 0x7fff
+#define EXT_EXP_INF 0x7fff
+#define EXT_EXP_NAN 0x7fff
#if 0
#define EXT_QUIETNAN (1 << 30)
struct ieee_ext extu_ext;
};
+#define extu_exp extu_ext.ext_exp
+#define extu_sign extu_ext.ext_sign
+#define extu_fracl extu_ext.ext_fracl
+#define extu_frach extu_ext.ext_frach
+
+#define LDBL_NBIT 0x80000000
+#define mask_nbit_l(u) ((u).extu_frach &= ~LDBL_NBIT)
+
#endif
* A NaN is a `signalling NaN' if its QUIETNAN bit is clear in its
* high fraction; if the bit is set, it is a `quiet NaN'.
*/
-#define EXT_EXP_INFNAN 32767
+#define EXT_EXP_INFNAN 0x7fff
+#define EXT_EXP_INF 0x7fff
+#define EXT_EXP_NAN 0x7fff
#if 0
#define EXT_QUIETNAN (1 << 30)
long double extu_ld;
struct ieee_ext extu_ext;
};
+
+#define extu_exp extu_ext.ext_exp
+#define extu_sign extu_ext.ext_sign
+#define extu_fracl extu_ext.ext_fracl
+#define extu_frach extu_ext.ext_frach
+
+#define LDBL_NBIT 0x80000000
+#define mask_nbit_l(u) ((u).extu_frach &= ~LDBL_NBIT)
projects / dragonfly.git / commitdiff