original files of m241
This commit is contained in:
neelimap
@@ -0,0 +1,215 @@
|
||||
# Makefile - makefile for target/src/math
|
||||
#
|
||||
# Copyright (c) 1996-2005, 2007-2010, 2013, 2014, 2018 Wind River Systems, Inc.
|
||||
#
|
||||
# The right to copy, distribute, modify or otherwise make use
|
||||
# of this software may be licensed only pursuant to the terms
|
||||
# of an applicable Wind River license agreement.
|
||||
#
|
||||
# modification history
|
||||
# --------------------
|
||||
# 03f,28nov17,mst Fix for CSO-31359 - PENTIUM4 first decimal digit overflows
|
||||
# to ":" in dtoa() conversion
|
||||
# 03e,07aug14,mcm Fix for VXW6-81137; floating point processing for printf has
|
||||
# inaccuracies; using BSD functionality
|
||||
# 03d,22apr13,yjw ARMARCH7M support
|
||||
# 03c,21sep10,sem Add new IA cpus. Remove pentium2/3
|
||||
# 03b,06sep10,fan added floatFuncBind.o for LP64 mode.
|
||||
# 03a,15jun10,m_h Thumb-2 Support
|
||||
# 02z,23jul09,mze removing PAD64
|
||||
# 02y,19jun09,mze replacing PAD64
|
||||
# 02x,29may09,cww Removed fppLib and fppShow
|
||||
# 02w,31mar09,jb Adding PAD64 support
|
||||
# 02w,03mar09,m_h Add support for ARMARCH7
|
||||
# 02v,30sep08,gls removed build of fioLib.o if _WRS_CONFIG_FORMATTED_OUT is
|
||||
# undefined (Defect #136485)
|
||||
# 02u,19sep08,dtr Tool e500v2 should use WRS_SPE_SUPPORT. WIND00134979
|
||||
# 02t,19aug08,mze TOOL changes
|
||||
# 02s,02apr07,pes Correct make warnings about Extraneous text after 'else'
|
||||
# and 'endif' directives.
|
||||
# 02r,12mar07,pes Update for MIPS library restructure.
|
||||
# 02q,22feb07,pes Add support for new MIPS lib structure.
|
||||
# 02p,27sep05,kk make fioBaseLib.c build the way as as fioLib.c. (SPR# 112910)
|
||||
# 02o,15sep05,mcm Fix for SPR#103295 - adding fioLib.c etc. to DOC_FILES
|
||||
# 02n,02aug05,rhe Break up for scailability
|
||||
# 02m,07apr05,jb Adding ARMARCH6
|
||||
# 02l,28oct04,to removed non-PAL support.
|
||||
# 02k,13apr04,rec PPC970 support
|
||||
# 02j,04aug04,scm coprocessor support modifications for ARM.
|
||||
# 02i,21apr04,h_k SH coprocessor support modifications.
|
||||
# 02h,13apr04,zmm Pentium coprocessor support modifications.
|
||||
# 02g,17mar04,pes Added PAL support for MIPS.
|
||||
# 02f,22mar04,mil Added -D_WRS_SPE_SUPPORT for fioLib.c for PPC32 sf$(TOOL).
|
||||
# 02e,13feb04,pch PPC32 floating point adjustments
|
||||
# 02d,19nov03,c_c Do not rely on VX_CPU_FAMILY setting.
|
||||
# 02c,19aug03,dbt Added PAL support for VxSim.
|
||||
# 02b,12dec02,pch Add PPC32; move FP files from src/os
|
||||
# 02a,22apr03,dbt Added SIMLINUX support.
|
||||
# 01z,18dec01,to add ARMARCH5(_T) support
|
||||
# 01y,29nov01,dee add Coldfire libraries
|
||||
# 01x,23aug01,pch Add PPC440
|
||||
# 01w,29oct01,sn moved contents of target/src/gnu/libm here
|
||||
# 01v,12oct01,tam added repackaging support
|
||||
# 01u,02oct01,sn added cabs/hypot support for arches that lacked it
|
||||
# 01t,23jul01,scm change XScale name to conform to coding standards...
|
||||
# 01s,04may01,scm add STRONGARM support
|
||||
# 01r,26oct00,s_m renamed PPC405 cpu types
|
||||
# 01q,06jul00,alp Added PPC405 Support
|
||||
# 01p,11dec00,scm replace references to ARMSA2 with XScale
|
||||
# 01o,31aug00,scm add sa2 support...
|
||||
# 01n,20jan99,cdp removed support for old ARM libraries.
|
||||
# 01m,13nov98,cdp added generic ARMARCH3/4/4_T support.
|
||||
# 01l,20mar01,zl removed obsolete SH CPU types.
|
||||
# 01k,03mar00,zl merged SH support from T1
|
||||
# 01j,28dec00,zmm add support for MIPS32sfgnu
|
||||
# 01i,27dec00,roz add support for MIPS64gnu
|
||||
# 01h,18aug98,tpr added PowerPC EC 603 support.
|
||||
# 01g,09apr98,hdn added support for Pentium.
|
||||
# 01f,07apr97,tpr removed cabs.o cbrt.o rint.o for the I960.
|
||||
# 01e,26mar97,cdp added ARM7TDMI_T support.
|
||||
# 01d,04mar97,kkk added R4650 support
|
||||
# 01c,03dec96,cdp added OBJS_ARM*.
|
||||
# 01b,26nov96,pr added SIMSPARCSOLARIS support
|
||||
# 01a,18jun96,yp created from 02c of MakeSkel
|
||||
#
|
||||
# DESCRIPTION
|
||||
# This file contains the makefile rules for building the non-ANSI portable
|
||||
# math routines.
|
||||
#
|
||||
# Not all modules are used in all architectures.
|
||||
#
|
||||
#*/
|
||||
|
||||
TGT_DIR=$(WIND_BASE)/target
|
||||
|
||||
LIB_BASE_NAME = fp
|
||||
|
||||
DOC_FILES= fioBaseLib.c fioLib.c
|
||||
|
||||
include $(TGT_DIR)/h/make/defs.library
|
||||
|
||||
ifeq ($(VX_CPU_FAMILY),mips)
|
||||
ifeq ($(BUILD_FP_PREFIX),sf)
|
||||
OBJS_RINT = rint.o
|
||||
else
|
||||
OBJS_RINT =
|
||||
endif
|
||||
endif
|
||||
|
||||
|
||||
OBJS_MIPSI3 = fpType.o cabs.o
|
||||
OBJS_MIPSI2 = fpType.o cabs.o $OBJS_RINT
|
||||
OBJS_MIPSI32 = fpType.o cabs.o $OBJS_RINT
|
||||
OBJS_MIPSI32R2 = fpType.o cabs.o $OBJS_RINT
|
||||
OBJS_MIPSI64 = fpType.o cabs.o $OBJS_RINT
|
||||
OBJS_MIPSI64R2 = fpType.o cabs.o $OBJS_RINT
|
||||
|
||||
OBJS_I960KA = fpType.o cabs.o
|
||||
OBJS_I960KB = fpType.o cabs.o
|
||||
OBJS_I960CA = fpType.o cabs.o
|
||||
OBJS_I960JX = fpType.o cabs.o
|
||||
OBJS_I960HX = fpType.o cabs.o
|
||||
|
||||
OBJS_MC68000 = fpType.o
|
||||
OBJS_MC68010 = fpType.o
|
||||
OBJS_MC68020 = fpType.o
|
||||
OBJS_MC68040 = fpType.o
|
||||
OBJS_MC68LC040 = fpType.o
|
||||
OBJS_MC68060 = fpType.o
|
||||
OBJS_CPU32 = fpType.o
|
||||
|
||||
OBJS_R3000 = fpType.o
|
||||
OBJS_R4000 = fpType.o
|
||||
OBJS_R4650 = fpType.o
|
||||
OBJS_R4000_64 = fpType.o
|
||||
|
||||
OBJS_SPARC = fpType.o cabs.o cbrt.o rint.o
|
||||
OBJS_SPARClite = fpType.o cabs.o cbrt.o rint.o
|
||||
OBJS_SIMLINUX = fpType.o cabs.o cbrt.o rint.o
|
||||
OBJS_SIMNT = fpType.o cabs.o cbrt.o rint.o
|
||||
OBJS_SIMSPARCSOLARIS = fpType.o cabs.o cbrt.o rint.o
|
||||
|
||||
OBJS_G100 = fpType.o
|
||||
OBJS_G200 = fpType.o
|
||||
|
||||
OBJS_I80386 = fpType.o
|
||||
OBJS_I80486 = fpType.o
|
||||
OBJS_PENTIUM = fpType.o
|
||||
OBJS_PENTIUM4 = fpType.o
|
||||
OBJS_CORE2 = fpType.o
|
||||
OBJS_X86_64 = fpType.o
|
||||
OBJS_ATOM = fpType.o
|
||||
OBJS_CORE = fpType.o
|
||||
OBJS_NEHALEM = fpType.o
|
||||
|
||||
OBJS_AM29030 = fpType.o
|
||||
OBJS_AM29200 = fpType.o
|
||||
|
||||
OBJS_PPC32 = fpType.o cabs.o
|
||||
OBJS_PPC603 = $(OBJS_PPC32)
|
||||
OBJS_PPC604 = $(OBJS_PPC32)
|
||||
OBJS_PPC970 = $(OBJS_PPC32)
|
||||
|
||||
OBJS_ARMARCH4 = fpType.o cabs.o
|
||||
OBJS_ARMARCH4_T = fpType.o cabs.o
|
||||
OBJS_ARMARCH5 = fpType.o cabs.o
|
||||
OBJS_ARMARCH5_T = fpType.o cabs.o
|
||||
OBJS_ARMARCH6 = fpType.o cabs.o
|
||||
OBJS_ARMARCH6_T = fpType.o cabs.o
|
||||
OBJS_ARMARCH7 = fpType.o cabs.o
|
||||
OBJS_ARMARCH7_T2 = fpType.o cabs.o
|
||||
OBJS_ARMARCH7M = fpType.o cabs.o
|
||||
OBJS_STRONGARM = fpType.o cabs.o
|
||||
OBJS_XSCALE = fpType.o cabs.o
|
||||
|
||||
OBJS_SH7000 = fpType.o cabs.o
|
||||
OBJS_SH7600 = fpType.o cabs.o
|
||||
OBJS_SH7700 = fpType.o cabs.o
|
||||
OBJS_SH7750 = fpType.o cabs.o
|
||||
|
||||
OBJS_MCF5200 = fpType.o
|
||||
OBJS_MCF5400 = fpType.o
|
||||
|
||||
OBJS = $(OBJS_$(CPU)) e_fmod.o e_pow.o e_remainder.o \
|
||||
e_sqrt.o s_rint.o s_scalbn.o fioBaseLib.o floatLib.o \
|
||||
bsd_cvt.o mprec.o dtoa.o floatFuncBind.o
|
||||
|
||||
ifdef _WRS_CONFIG_FORMATTED_OUT
|
||||
OBJS += fioLib.o
|
||||
endif
|
||||
|
||||
include $(TGT_DIR)/h/make/rules.library
|
||||
|
||||
# These must be after the inclusion of rules.library to
|
||||
# overwrite/add to the CFLAGS.
|
||||
|
||||
#PPC604 has _WRS_ALTIVEC_SUPPORT defined in ppc604.h so not needed
|
||||
#ifeq ($(CPU),PPC604)
|
||||
#CFLAGS_fioLib.o = -D_WRS_ALTIVEC_SUPPORT
|
||||
#CFLAGS_fioBaseLib.o = -D_WRS_ALTIVEC_SUPPORT
|
||||
#endif
|
||||
|
||||
ifeq ($(CPU),PPC85XX)
|
||||
CFLAGS_fioLib.o = -D_WRS_SPE_SUPPORT
|
||||
CFLAGS_fioBaseLib.o = -D_WRS_SPE_SUPPORT
|
||||
endif
|
||||
|
||||
ifeq ($(CPU),PPC32)
|
||||
ifeq ($(_WRS_CONFIG_FP),e500v2)
|
||||
CFLAGS_fioLib.o = -D_WRS_SPE_SUPPORT
|
||||
CFLAGS_fioBaseLib.o = -D_WRS_SPE_SUPPORT
|
||||
endif
|
||||
ifeq ($(_WRS_CONFIG_FP),hard)
|
||||
CFLAGS_fioLib.o = $(OPTION_ALTIVEC_C)
|
||||
CFLAGS_fioBaseLib.o = $(OPTION_ALTIVEC_C)
|
||||
endif
|
||||
ifeq ($(_WRS_CONFIG_FP),soft)
|
||||
CFLAGS_fioLib.o = -D_WRS_SPE_SUPPORT
|
||||
CFLAGS_fioBaseLib.o = -D_WRS_SPE_SUPPORT
|
||||
endif
|
||||
endif
|
||||
|
||||
ifeq ($(CPU),PENTIUM4)
|
||||
CFLAGS_dtoa.o = -DCheck_FLT_ROUNDS -O1
|
||||
endif
|
||||
|
||||
@@ -0,0 +1,108 @@
|
||||
/*-
|
||||
* Copyright (c) 1990 The Regents of the University of California.
|
||||
* All rights reserved.
|
||||
*
|
||||
* This code is derived from software contributed to Berkeley by
|
||||
* Chris Torek.
|
||||
*
|
||||
* 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. All advertising materials mentioning features or use of this software
|
||||
* must display the following acknowledgement:
|
||||
* This product includes software developed by the University of
|
||||
* California, Berkeley and its contributors.
|
||||
* 4. 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.
|
||||
*/
|
||||
|
||||
#include <stdio.h>
|
||||
|
||||
/*
|
||||
* * Flags used during conversion.
|
||||
* */
|
||||
#define ALT 0x001 /* alternate form */
|
||||
#define HEXPREFIX 0x002 /* add 0x or 0X prefix */
|
||||
#define LADJUST 0x004 /* left adjustment */
|
||||
#define LONGDBL 0x008 /* long double; unimplemented */
|
||||
#define LONGINT 0x010 /* long integer */
|
||||
#define QUADINT 0x020 /* quad integer */
|
||||
#define SHORTINT 0x040 /* short integer */
|
||||
#define ZEROPAD 0x080 /* zero (as opposed to blank) pad */
|
||||
#define FPT 0x100 /* Floating point number */
|
||||
|
||||
void
|
||||
setAlt(int flags)
|
||||
{
|
||||
flags |= ALT;
|
||||
}
|
||||
|
||||
void
|
||||
setFpt(int flags)
|
||||
{
|
||||
flags |= FPT;
|
||||
}
|
||||
|
||||
extern void __dtoa(double value, int mode, int ndigits, int *decpt, int *dsgn, char **rve, char **buf, int fl_type);
|
||||
|
||||
void
|
||||
bsd_cvt(value, ndigits, flags, sign, decpt, ch, length, buffer)
|
||||
double value;
|
||||
int ndigits, flags, *decpt, ch, *length;
|
||||
char *sign, **buffer;
|
||||
{
|
||||
int mode, dsgn;
|
||||
char *digits, *bp, *rve;
|
||||
digits = *buffer;
|
||||
|
||||
if (ch == 'f') {
|
||||
mode = 3; /* ndigits after the decimal point */
|
||||
} else {
|
||||
/* To obtain ndigits after the decimal point for the 'e'
|
||||
* and 'E' formats, round to ndigits + 1 significant
|
||||
* figures.
|
||||
*/
|
||||
if (ch == 'e' || ch == 'E') {
|
||||
ndigits++;
|
||||
}
|
||||
mode = 2; /* ndigits significant digits */
|
||||
}
|
||||
|
||||
if (value < 0) {
|
||||
value = -value;
|
||||
*sign = '-';
|
||||
} else
|
||||
*sign = '\000';
|
||||
__dtoa(value, mode, ndigits, decpt, &dsgn, &rve, buffer, (ch == 'f' || ch =='F')? 1: 0);
|
||||
if ((ch != 'g' && ch != 'G') || flags & ALT) { /* Print trailing zeros */
|
||||
bp = digits + ndigits;
|
||||
if (ch == 'f') {
|
||||
if (*digits == '0' && value)
|
||||
*decpt = -ndigits + 1;
|
||||
bp += *decpt;
|
||||
}
|
||||
if (value == 0) /* kludge for __dtoa irregularity */
|
||||
rve = bp;
|
||||
while (rve < bp)
|
||||
*rve++ = '0';
|
||||
}
|
||||
*length = (int)(rve - digits);
|
||||
}
|
||||
@@ -0,0 +1,217 @@
|
||||
/* cabs.c - math routines */
|
||||
|
||||
/* Copyright 1992 Wind River Systems, Inc. */
|
||||
|
||||
/*
|
||||
modification history
|
||||
--------------------
|
||||
01d,15feb93,jdi made hypot() NOMANUAL.
|
||||
01c,14oct92,jdi made cabs() NOMANUAL.
|
||||
01b,13oct92,jdi mangen fixes.
|
||||
01a,30jul92,smb documentation.
|
||||
*/
|
||||
|
||||
/*
|
||||
DESCRIPTION
|
||||
* Copyright (c) 1985 Regents of the University of California.
|
||||
* All rights reserved.
|
||||
*
|
||||
* Redistribution and use in source and binary forms are permitted
|
||||
* provided that the above copyright notice and this paragraph are
|
||||
* duplicated in all such forms and that any documentation,
|
||||
* advertising materials, and other materials related to such
|
||||
* distribution and use acknowledge that the software was developed
|
||||
* by the University of California, Berkeley. The name of the
|
||||
* University may not be used to endorse or promote products derived
|
||||
* from this software without specific prior written permission.
|
||||
* THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
|
||||
* IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
|
||||
* WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
|
||||
*
|
||||
* All recipients should regard themselves as participants in an ongoing
|
||||
* research project and hence should feel obligated to report their
|
||||
* experiences (good or bad) with these elementary function codes, using
|
||||
* the sendbug(8) program, to the authors.
|
||||
*
|
||||
#ifndef lint
|
||||
static char sccsid[] = "@(#)cabs.c 5.3 (Berkeley) 6/30/88";
|
||||
#endif * not lint *
|
||||
|
||||
SEE ALSO: American National Standard X3.159-1989
|
||||
|
||||
NOMANUAL
|
||||
*/
|
||||
|
||||
|
||||
#include "vxWorks.h"
|
||||
#include "math.h"
|
||||
|
||||
#if defined(vax)||defined(tahoe) /* VAX D format */
|
||||
#ifdef vax
|
||||
#define _0x(A,B) 0x/**/A/**/B
|
||||
#else /* vax */
|
||||
#define _0x(A,B) 0x/**/B/**/A
|
||||
#endif /* vax */
|
||||
/* static double */
|
||||
/* r2p1hi = 2.4142135623730950345E0 , Hex 2^ 2 * .9A827999FCEF32 */
|
||||
/* r2p1lo = 1.4349369327986523769E-17 , Hex 2^-55 * .84597D89B3754B */
|
||||
/* sqrt2 = 1.4142135623730950622E0 ; Hex 2^ 1 * .B504F333F9DE65 */
|
||||
static long r2p1hix[] = { _0x(8279,411a), _0x(ef32,99fc)};
|
||||
static long r2p1lox[] = { _0x(597d,2484), _0x(754b,89b3)};
|
||||
static long sqrt2x[] = { _0x(04f3,40b5), _0x(de65,33f9)};
|
||||
#define r2p1hi (*(double*)r2p1hix)
|
||||
#define r2p1lo (*(double*)r2p1lox)
|
||||
#define sqrt2 (*(double*)sqrt2x)
|
||||
#else /* defined(vax)||defined(tahoe) */
|
||||
static double
|
||||
r2p1hi = 2.4142135623730949234E0 , /*Hex 2^1 * 1.3504F333F9DE6 */
|
||||
r2p1lo = 1.2537167179050217666E-16 , /*Hex 2^-53 * 1.21165F626CDD5 */
|
||||
sqrt2 = 1.4142135623730951455E0 ; /*Hex 2^ 0 * 1.6A09E667F3BCD */
|
||||
#endif /* defined(vax)||defined(tahoe) */
|
||||
|
||||
/******************************************************************************
|
||||
*
|
||||
* hypot - return the square root of X^2 + Y^2 where Z = X + iY
|
||||
*
|
||||
* This routine returns the square root of X^2 + Y^2 where Z=X+iY
|
||||
* in double precision (VAX D format 56 bits, IEEE double 53 bits).
|
||||
*
|
||||
* Required system support functions: copysign(), finite(), scalb(), sqrt().
|
||||
*
|
||||
* Method:
|
||||
* 1. replace x by |x| and y by |y|, and swap x and
|
||||
* y if y > x (hence x is never smaller than y).
|
||||
* 2. Hypot(x,y) is computed by:
|
||||
* .CS
|
||||
* Case I, x/y > 2
|
||||
*
|
||||
* y
|
||||
* hypot = x + -----------------------------
|
||||
* 2
|
||||
* sqrt ( 1 + [x/y] ) + x/y
|
||||
*
|
||||
* Case II, x/y <= 2
|
||||
* y
|
||||
* hypot = x + --------------------------------------------------
|
||||
* 2
|
||||
* [x/y] - 2
|
||||
* (sqrt(2)+1) + (x-y)/y + -----------------------------
|
||||
* 2
|
||||
* sqrt ( 1 + [x/y] ) + sqrt(2)
|
||||
* .CE
|
||||
*
|
||||
* Special cases:
|
||||
* .CS
|
||||
* hypot(x,y) is INF if x or y is +INF or -INF; else
|
||||
* hypot(x,y) is NAN if x or y is NAN.
|
||||
* .CE
|
||||
*
|
||||
* ACCURACY:
|
||||
* hypot() returns the sqrt(x^2+y^2) with error less than 1 ulps (units
|
||||
* in the last place). See Kahan's "Interval Arithmetic Options in the
|
||||
* Proposed IEEE Floating Point Arithmetic Standard", Interval Mathematics
|
||||
* 1980, Edited by Karl L.E. Nickel, pp 99-128. (A faster but less accurate
|
||||
* code follows in comments.) In a test run with 500,000 random arguments
|
||||
* on a VAX, the maximum observed error was 0.959 ulps.
|
||||
*
|
||||
* CONSTANTS:
|
||||
* The hexadecimal values are the intended ones for the following constants.
|
||||
* The decimal values may be used, provided that the compiler will convert
|
||||
* from decimal to binary accurately enough to produce the hexadecimal values
|
||||
* shown.
|
||||
*
|
||||
* RETURNS: The square root of X^2 + Y^2 where Z = X + iY.
|
||||
*
|
||||
* INTERNAL:
|
||||
* Coded in C by K.C. Ng, 11/28/84;
|
||||
* Revised by K.C. Ng, 7/12/85.
|
||||
*/
|
||||
|
||||
double hypot(x,y)
|
||||
double x, y;
|
||||
{
|
||||
static double zero=0, one=1,
|
||||
small=1.0E-18; /* fl(1+small)==1 */
|
||||
static ibig=30; /* fl(1+2**(2*ibig))==1 */
|
||||
double copysign(),scalb(),logb(),sqrt(),t,r;
|
||||
int finite(), exp;
|
||||
|
||||
if(finite(x))
|
||||
if(finite(y))
|
||||
{
|
||||
x=copysign(x,one);
|
||||
y=copysign(y,one);
|
||||
if(y > x)
|
||||
{ t=x; x=y; y=t; }
|
||||
if(x == zero) return(zero);
|
||||
if(y == zero) return(x);
|
||||
exp= logb(x);
|
||||
if(exp-(int)logb(y) > ibig )
|
||||
/* raise inexact flag and return |x| */
|
||||
{ one+small; return(x); }
|
||||
|
||||
/* start computing sqrt(x^2 + y^2) */
|
||||
r=x-y;
|
||||
if(r>y) { /* x/y > 2 */
|
||||
r=x/y;
|
||||
r=r+sqrt(one+r*r); }
|
||||
else { /* 1 <= x/y <= 2 */
|
||||
r/=y; t=r*(r+2.0);
|
||||
r+=t/(sqrt2+sqrt(2.0+t));
|
||||
r+=r2p1lo; r+=r2p1hi; }
|
||||
|
||||
r=y/r;
|
||||
return(x+r);
|
||||
|
||||
}
|
||||
|
||||
else if(y==y) /* y is +-INF */
|
||||
return(copysign(y,one));
|
||||
else
|
||||
return(y); /* y is NaN and x is finite */
|
||||
|
||||
else if(x==x) /* x is +-INF */
|
||||
return (copysign(x,one));
|
||||
else if(finite(y))
|
||||
return(x); /* x is NaN, y is finite */
|
||||
#if !defined(vax)&&!defined(tahoe)
|
||||
else if(y!=y) return(y); /* x and y is NaN */
|
||||
#endif /* !defined(vax)&&!defined(tahoe) */
|
||||
else return(copysign(y,one)); /* y is INF */
|
||||
}
|
||||
|
||||
/*****************************************************************************
|
||||
*
|
||||
* cabs - return the absolute value of the complex number Z = X + iY
|
||||
*
|
||||
* This routine returns the absolute value of the complex number Z = X + iY,
|
||||
* in double precision (VAX D format 56 bits, IEEE double 53 bits).
|
||||
*
|
||||
* Required kernel function: hypot()
|
||||
*
|
||||
* Method :
|
||||
* .CS
|
||||
* cabs(z) = hypot(x,y) .
|
||||
* .CE
|
||||
*
|
||||
* RETURNS: The absolute value of Z = X + iY.
|
||||
*
|
||||
* INTERNAL:
|
||||
* Coded in C by K.C. Ng, 11/28/84.
|
||||
* Revised by K.C. Ng, 7/12/85.
|
||||
*
|
||||
* NOMANUAL
|
||||
*/
|
||||
|
||||
double cabs(z)
|
||||
struct { double x, y;} z;
|
||||
{
|
||||
return hypot(z.x,z.y);
|
||||
}
|
||||
|
||||
double
|
||||
z_abs(z)
|
||||
struct { double x,y;} *z;
|
||||
{
|
||||
return hypot(z->x,z->y);
|
||||
}
|
||||
@@ -0,0 +1,124 @@
|
||||
/* cbrt.c - software cube-root function */
|
||||
|
||||
/* Copyright 1992 Wind River Systems, Inc. */
|
||||
|
||||
/*
|
||||
modification history
|
||||
--------------------
|
||||
01a,30jul92,smb documentation.
|
||||
*/
|
||||
|
||||
/*
|
||||
DESCRIPTION
|
||||
* Copyright (c) 1985 Regents of the University of California.
|
||||
* All rights reserved.
|
||||
*
|
||||
* Redistribution and use in source and binary forms are permitted
|
||||
* provided that the above copyright notice and this paragraph are
|
||||
* duplicated in all such forms and that any documentation,
|
||||
* advertising materials, and other materials related to such
|
||||
* distribution and use acknowledge that the software was developed
|
||||
* by the University of California, Berkeley. The name of the
|
||||
* University may not be used to endorse or promote products derived
|
||||
* from this software without specific prior written permission.
|
||||
* THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
|
||||
* IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
|
||||
* WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
|
||||
*
|
||||
* All recipients should regard themselves as participants in an ongoing
|
||||
* research project and hence should feel obligated to report their
|
||||
* experiences (good or bad) with these elementary function codes, using
|
||||
* the sendbug(8) program, to the authors.
|
||||
*
|
||||
#ifndef lint
|
||||
|
||||
SEE ALSO: American National Standard X3.159-1989
|
||||
*/
|
||||
|
||||
#include "vxWorks.h"
|
||||
#include "math.h"
|
||||
|
||||
#if !defined(vax)&&!defined(tahoe)
|
||||
|
||||
static unsigned long B1 = 715094163, /* B1 = (682-0.03306235651)*2**20 */
|
||||
B2 = 696219795; /* B2 = (664-0.03306235651)*2**20 */
|
||||
static double
|
||||
C= 19./35.,
|
||||
D= -864./1225.,
|
||||
E= 99./70.,
|
||||
F= 45./28.,
|
||||
G= 5./14.;
|
||||
|
||||
/*****************************************************************************
|
||||
* cbrt -
|
||||
*
|
||||
* kahan's cube root (53 bits IEEE double precision)
|
||||
* for IEEE machines only
|
||||
* coded in C by K.C. Ng, 4/30/85
|
||||
*
|
||||
* Accuracy:
|
||||
* better than 0.667 ulps according to an error analysis. Maximum
|
||||
* error observed was 0.666 ulps in an 1,000,000 random arguments test.
|
||||
*
|
||||
* Warning: this code is semi machine dependent; the ordering of words in
|
||||
* a floating point number must be known in advance. I assume that the
|
||||
* long interger at the address of a floating point number will be the
|
||||
* leading 32 bits of that floating point number (i.e., sign, exponent,
|
||||
* and the 20 most significant bits).
|
||||
* On a National machine, it has different ordering; therefore, this code
|
||||
* must be compiled with flag -DNATIONAL.
|
||||
*/
|
||||
|
||||
double cbrt(x)
|
||||
double x;
|
||||
{
|
||||
double r,s,t=0.0,w;
|
||||
unsigned long *px = (unsigned long *) &x,
|
||||
*pt = (unsigned long *) &t,
|
||||
mexp,sign;
|
||||
|
||||
#ifdef national /* ordering of words in a floating points number */
|
||||
int n0=1,n1=0;
|
||||
#else /* national */
|
||||
int n0=0,n1=1;
|
||||
#endif /* national */
|
||||
|
||||
mexp=px[n0]&0x7ff00000;
|
||||
if(mexp==0x7ff00000) return(x); /* cbrt(NaN,INF) is itself */
|
||||
if(x==0.0) return(x); /* cbrt(0) is itself */
|
||||
|
||||
sign=px[n0]&0x80000000; /* sign= sign(x) */
|
||||
px[n0] ^= sign; /* x=|x| */
|
||||
|
||||
|
||||
/* rough cbrt to 5 bits */
|
||||
if(mexp==0) /* subnormal number */
|
||||
{pt[n0]=0x43500000; /* set t= 2**54 */
|
||||
t*=x; pt[n0]=pt[n0]/3+B2;
|
||||
}
|
||||
else
|
||||
pt[n0]=px[n0]/3+B1;
|
||||
|
||||
|
||||
/* new cbrt to 23 bits, may be implemented in single precision */
|
||||
r=t*t/x;
|
||||
s=C+r*t;
|
||||
t*=G+F/(s+E+D/s);
|
||||
|
||||
/* chopped to 20 bits and make it larger than cbrt(x) */
|
||||
pt[n1]=0; pt[n0]+=0x00000001;
|
||||
|
||||
|
||||
/* one step newton iteration to 53 bits with error less than 0.667 ulps */
|
||||
s=t*t; /* t*t is exact */
|
||||
r=x/s;
|
||||
w=t+t;
|
||||
r=(r-t)/(w+r); /* r-t is exact */
|
||||
t=t+t*r;
|
||||
|
||||
|
||||
/* retore the sign bit */
|
||||
pt[n0] |= sign;
|
||||
return(t);
|
||||
}
|
||||
#endif
|
||||
@@ -0,0 +1,935 @@
|
||||
/****************************************************************
|
||||
*
|
||||
* The author of this software is David M. Gay.
|
||||
*
|
||||
* Copyright (c) 1991, 2006 by AT&T.
|
||||
*
|
||||
* Permission to use, copy, modify, and distribute this software for any
|
||||
* purpose without fee is hereby granted, provided that this entire notice
|
||||
* is included in all copies of any software which is or includes a copy
|
||||
* or modification of this software and in all copies of the supporting
|
||||
* documentation for such software.
|
||||
*
|
||||
* THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
|
||||
* WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY
|
||||
* REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
|
||||
* OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
|
||||
*
|
||||
***************************************************************/
|
||||
|
||||
/* Please send bug reports to
|
||||
David M. Gay
|
||||
AT&T Bell Laboratories, Room 2C-463
|
||||
600 Mountain Avenue
|
||||
Murray Hill, NJ 07974-2070
|
||||
U.S.A.
|
||||
dmg@research.att.com or research!dmg
|
||||
*/
|
||||
|
||||
#include "mprec.h"
|
||||
#include <string.h>
|
||||
#include <stdlib.h>
|
||||
|
||||
static int
|
||||
_DEFUN (quorem,
|
||||
(b, S),
|
||||
_Jv_Bigint * b _AND _Jv_Bigint * S)
|
||||
{
|
||||
int n;
|
||||
__Long borrow, y;
|
||||
__ULong carry, q, ys;
|
||||
__ULong *bx, *bxe, *sx, *sxe;
|
||||
#ifdef Pack_32
|
||||
__Long z;
|
||||
__ULong si, zs;
|
||||
#endif
|
||||
|
||||
n = S->_wds;
|
||||
#ifdef DEBUG
|
||||
/*debug*/ if (b->_wds > n)
|
||||
/*debug*/ Bug ("oversize b in quorem");
|
||||
#endif
|
||||
if (b->_wds < n)
|
||||
return 0;
|
||||
sx = S->_x;
|
||||
sxe = sx + --n;
|
||||
bx = b->_x;
|
||||
bxe = bx + n;
|
||||
q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
|
||||
#ifdef DEBUG
|
||||
/*debug*/ if (q > 9)
|
||||
/*debug*/ Bug ("oversized quotient in quorem");
|
||||
#endif
|
||||
if (q)
|
||||
{
|
||||
borrow = 0;
|
||||
carry = 0;
|
||||
do
|
||||
{
|
||||
#ifdef Pack_32
|
||||
si = *sx++;
|
||||
ys = (si & 0xffff) * q + carry;
|
||||
zs = (si >> 16) * q + (ys >> 16);
|
||||
carry = zs >> 16;
|
||||
y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
|
||||
borrow = y >> 16;
|
||||
Sign_Extend (borrow, y);
|
||||
z = (*bx >> 16) - (zs & 0xffff) + borrow;
|
||||
borrow = z >> 16;
|
||||
Sign_Extend (borrow, z);
|
||||
Storeinc (bx, z, y);
|
||||
#else
|
||||
ys = *sx++ * q + carry;
|
||||
carry = ys >> 16;
|
||||
y = *bx - (ys & 0xffff) + borrow;
|
||||
borrow = y >> 16;
|
||||
Sign_Extend (borrow, y);
|
||||
*bx++ = y & 0xffff;
|
||||
#endif
|
||||
}
|
||||
while (sx <= sxe);
|
||||
if (!*bxe)
|
||||
{
|
||||
bx = b->_x;
|
||||
while (--bxe > bx && !*bxe)
|
||||
--n;
|
||||
b->_wds = n;
|
||||
}
|
||||
}
|
||||
if (cmp (b, S) >= 0)
|
||||
{
|
||||
q++;
|
||||
borrow = 0;
|
||||
carry = 0;
|
||||
bx = b->_x;
|
||||
sx = S->_x;
|
||||
do
|
||||
{
|
||||
#ifdef Pack_32
|
||||
si = *sx++;
|
||||
ys = (si & 0xffff) + carry;
|
||||
zs = (si >> 16) + (ys >> 16);
|
||||
carry = zs >> 16;
|
||||
y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
|
||||
borrow = y >> 16;
|
||||
Sign_Extend (borrow, y);
|
||||
z = (*bx >> 16) - (zs & 0xffff) + borrow;
|
||||
borrow = z >> 16;
|
||||
Sign_Extend (borrow, z);
|
||||
Storeinc (bx, z, y);
|
||||
#else
|
||||
ys = *sx++ + carry;
|
||||
carry = ys >> 16;
|
||||
y = *bx - (ys & 0xffff) + borrow;
|
||||
borrow = y >> 16;
|
||||
Sign_Extend (borrow, y);
|
||||
*bx++ = y & 0xffff;
|
||||
#endif
|
||||
}
|
||||
while (sx <= sxe);
|
||||
bx = b->_x;
|
||||
bxe = bx + n;
|
||||
if (!*bxe)
|
||||
{
|
||||
while (--bxe > bx && !*bxe)
|
||||
--n;
|
||||
b->_wds = n;
|
||||
}
|
||||
}
|
||||
return q;
|
||||
}
|
||||
|
||||
#ifdef DEBUG
|
||||
#include <stdio.h>
|
||||
|
||||
void
|
||||
print (_Jv_Bigint * b)
|
||||
{
|
||||
int i, wds;
|
||||
__Long *x, y;
|
||||
wds = b->_wds;
|
||||
x = b->_x+wds;
|
||||
i = 0;
|
||||
do
|
||||
{
|
||||
x--;
|
||||
fprintf (stderr, "%08x", *x);
|
||||
}
|
||||
while (++i < wds);
|
||||
fprintf (stderr, "\n");
|
||||
}
|
||||
#endif
|
||||
|
||||
/* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
|
||||
*
|
||||
* Inspired by "How to Print Floating-Point Numbers Accurately" by
|
||||
* Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
|
||||
*
|
||||
* Modifications:
|
||||
* 1. Rather than iterating, we use a simple numeric overestimate
|
||||
* to determine k = floor(log10(d)). We scale relevant
|
||||
* quantities using O(log2(k)) rather than O(k) multiplications.
|
||||
* 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
|
||||
* try to generate digits strictly left to right. Instead, we
|
||||
* compute with fewer bits and propagate the carry if necessary
|
||||
* when rounding the final digit up. This is often faster.
|
||||
* 3. Under the assumption that input will be rounded nearest,
|
||||
* mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
|
||||
* That is, we allow equality in stopping tests when the
|
||||
* round-nearest rule will give the same floating-point value
|
||||
* as would satisfaction of the stopping test with strict
|
||||
* inequality.
|
||||
* 4. We remove common factors of powers of 2 from relevant
|
||||
* quantities.
|
||||
* 5. When converting floating-point integers less than 1e16,
|
||||
* we use floating-point arithmetic rather than resorting
|
||||
* to multiple-precision integers.
|
||||
* 6. When asked to produce fewer than 15 digits, we first try
|
||||
* to get by with floating-point arithmetic; we resort to
|
||||
* multiple-precision integer arithmetic only if we cannot
|
||||
* guarantee that the floating-point calculation has given
|
||||
* the correctly rounded result. For k requested digits and
|
||||
* "uniformly" distributed input, the probability is
|
||||
* something like 10^(k-15) that we must resort to the long
|
||||
* calculation.
|
||||
*/
|
||||
|
||||
|
||||
char *
|
||||
_DEFUN (_dtoa_r,
|
||||
(ptr, _d, mode, ndigits, decpt, sign, rve, float_type),
|
||||
struct _Jv_reent *ptr _AND
|
||||
double _d _AND
|
||||
int mode _AND
|
||||
int ndigits _AND
|
||||
int *decpt _AND
|
||||
int *sign _AND
|
||||
char **rve _AND
|
||||
int float_type)
|
||||
{
|
||||
/*
|
||||
float_type == 0 for double precision, 1 for float.
|
||||
|
||||
Arguments ndigits, decpt, sign are similar to those
|
||||
of ecvt and fcvt; trailing zeros are suppressed from
|
||||
the returned string. If not null, *rve is set to point
|
||||
to the end of the return value. If d is +-Infinity or NaN,
|
||||
then *decpt is set to 9999.
|
||||
|
||||
mode:
|
||||
0 ==> shortest string that yields d when read in
|
||||
and rounded to nearest.
|
||||
1 ==> like 0, but with Steele & White stopping rule;
|
||||
e.g. with IEEE P754 arithmetic , mode 0 gives
|
||||
1e23 whereas mode 1 gives 9.999999999999999e22.
|
||||
2 ==> max(1,ndigits) significant digits. This gives a
|
||||
return value similar to that of ecvt, except
|
||||
that trailing zeros are suppressed.
|
||||
3 ==> through ndigits past the decimal point. This
|
||||
gives a return value similar to that from fcvt,
|
||||
except that trailing zeros are suppressed, and
|
||||
ndigits can be negative.
|
||||
4-9 should give the same return values as 2-3, i.e.,
|
||||
4 <= mode <= 9 ==> same return as mode
|
||||
2 + (mode & 1). These modes are mainly for
|
||||
debugging; often they run slower but sometimes
|
||||
faster than modes 2-3.
|
||||
4,5,8,9 ==> left-to-right digit generation.
|
||||
6-9 ==> don't try fast floating-point estimate
|
||||
(if applicable).
|
||||
|
||||
> 16 ==> Floating-point arg is treated as single precision.
|
||||
|
||||
Values of mode other than 0-9 are treated as mode 0.
|
||||
|
||||
Sufficient space is allocated to the return value
|
||||
to hold the suppressed trailing zeros.
|
||||
*/
|
||||
|
||||
int bbits, b2, b5, be, dig, i, ieps, ilim0, j, j1, k, k0,
|
||||
k_check, leftright, m2, m5, s2, s5, try_quick;
|
||||
int ilim = 0, ilim1 = 0, spec_case = 0;
|
||||
union double_union d, d2, eps;
|
||||
__Long L;
|
||||
#ifndef Sudden_Underflow
|
||||
int denorm;
|
||||
__ULong x;
|
||||
#endif
|
||||
_Jv_Bigint *b, *b1, *delta, *mlo = NULL, *mhi, *S;
|
||||
double ds;
|
||||
char *s, *s0;
|
||||
|
||||
d.d = _d;
|
||||
|
||||
if (ptr->_result)
|
||||
{
|
||||
ptr->_result->_k = ptr->_result_k;
|
||||
ptr->_result->_maxwds = 1 << ptr->_result_k;
|
||||
Bfree (ptr, ptr->_result);
|
||||
ptr->_result = 0;
|
||||
}
|
||||
|
||||
if (word0 (d) & Sign_bit)
|
||||
{
|
||||
/* set sign for everything, including 0's and NaNs */
|
||||
*sign = 1;
|
||||
word0 (d) &= ~Sign_bit; /* clear sign bit */
|
||||
}
|
||||
else
|
||||
*sign = 0;
|
||||
|
||||
#if defined(IEEE_Arith) + defined(VAX)
|
||||
#ifdef IEEE_Arith
|
||||
if ((word0 (d) & Exp_mask) == Exp_mask)
|
||||
#else
|
||||
if (word0 (d) == 0x8000)
|
||||
#endif
|
||||
{
|
||||
/* Infinity or NaN */
|
||||
*decpt = 9999;
|
||||
s =
|
||||
#ifdef IEEE_Arith
|
||||
!word1 (d) && !(word0 (d) & 0xfffff) ? "Infinity" :
|
||||
#endif
|
||||
"NaN";
|
||||
if (rve)
|
||||
*rve =
|
||||
#ifdef IEEE_Arith
|
||||
s[3] ? s + 8 :
|
||||
#endif
|
||||
s + 3;
|
||||
return s;
|
||||
}
|
||||
#endif
|
||||
#ifdef IBM
|
||||
d.d += 0; /* normalize */
|
||||
#endif
|
||||
if (!d.d)
|
||||
{
|
||||
*decpt = 1;
|
||||
s = "0";
|
||||
if (rve)
|
||||
*rve = s + 1;
|
||||
return s;
|
||||
}
|
||||
|
||||
b = d2b (ptr, d.d, &be, &bbits);
|
||||
#ifdef Sudden_Underflow
|
||||
i = (int) (word0 (d) >> Exp_shift1 & (Exp_mask >> Exp_shift1));
|
||||
#else
|
||||
if ((i = (int) (word0 (d) >> Exp_shift1 & (Exp_mask >> Exp_shift1))))
|
||||
{
|
||||
#endif
|
||||
d2.d = d.d;
|
||||
word0 (d2) &= Frac_mask1;
|
||||
word0 (d2) |= Exp_11;
|
||||
#ifdef IBM
|
||||
if (j = 11 - hi0bits (word0 (d2) & Frac_mask))
|
||||
d2.d /= 1 << j;
|
||||
#endif
|
||||
|
||||
/* log(x) ~=~ log(1.5) + (x-1.5)/1.5
|
||||
* log10(x) = log(x) / log(10)
|
||||
* ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
|
||||
* log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
|
||||
*
|
||||
* This suggests computing an approximation k to log10(d) by
|
||||
*
|
||||
* k = (i - Bias)*0.301029995663981
|
||||
* + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
|
||||
*
|
||||
* We want k to be too large rather than too small.
|
||||
* The error in the first-order Taylor series approximation
|
||||
* is in our favor, so we just round up the constant enough
|
||||
* to compensate for any error in the multiplication of
|
||||
* (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
|
||||
* and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
|
||||
* adding 1e-13 to the constant term more than suffices.
|
||||
* Hence we adjust the constant term to 0.1760912590558.
|
||||
* (We could get a more accurate k by invoking log10,
|
||||
* but this is probably not worthwhile.)
|
||||
*/
|
||||
|
||||
i -= Bias;
|
||||
#ifdef IBM
|
||||
i <<= 2;
|
||||
i += j;
|
||||
#endif
|
||||
#ifndef Sudden_Underflow
|
||||
denorm = 0;
|
||||
}
|
||||
else
|
||||
{
|
||||
/* d is denormalized */
|
||||
|
||||
i = bbits + be + (Bias + (P - 1) - 1);
|
||||
x = i > 32 ? word0 (d) << (64 - i) | word1 (d) >> (i - 32)
|
||||
: word1 (d) << (32 - i);
|
||||
d2.d = x;
|
||||
word0 (d2) -= 31 * Exp_msk1; /* adjust exponent */
|
||||
i -= (Bias + (P - 1) - 1) + 1;
|
||||
denorm = 1;
|
||||
}
|
||||
#endif
|
||||
ds = (d2.d - 1.5) * 0.289529654602168 + 0.1760912590558 + i * 0.301029995663981;
|
||||
k = (int) ds;
|
||||
if (ds < 0. && ds != k)
|
||||
k--; /* want k = floor(ds) */
|
||||
k_check = 1;
|
||||
if (k >= 0 && k <= Ten_pmax)
|
||||
{
|
||||
if (d.d < tens[k])
|
||||
k--;
|
||||
k_check = 0;
|
||||
}
|
||||
j = bbits - i - 1;
|
||||
if (j >= 0)
|
||||
{
|
||||
b2 = 0;
|
||||
s2 = j;
|
||||
}
|
||||
else
|
||||
{
|
||||
b2 = -j;
|
||||
s2 = 0;
|
||||
}
|
||||
if (k >= 0)
|
||||
{
|
||||
b5 = 0;
|
||||
s5 = k;
|
||||
s2 += k;
|
||||
}
|
||||
else
|
||||
{
|
||||
b2 -= k;
|
||||
b5 = -k;
|
||||
s5 = 0;
|
||||
}
|
||||
if (mode < 0 || mode > 9)
|
||||
mode = 0;
|
||||
try_quick = 1;
|
||||
if (mode > 5)
|
||||
{
|
||||
mode -= 4;
|
||||
try_quick = 0;
|
||||
}
|
||||
leftright = 1;
|
||||
switch (mode)
|
||||
{
|
||||
case 0:
|
||||
case 1:
|
||||
ilim = ilim1 = -1;
|
||||
i = 18;
|
||||
ndigits = 0;
|
||||
break;
|
||||
case 2:
|
||||
leftright = 0;
|
||||
/* no break */
|
||||
case 4:
|
||||
if (ndigits <= 0)
|
||||
ndigits = 1;
|
||||
ilim = ilim1 = i = ndigits;
|
||||
break;
|
||||
case 3:
|
||||
leftright = 0;
|
||||
/* no break */
|
||||
case 5:
|
||||
i = ndigits + k + 1;
|
||||
ilim = i;
|
||||
ilim1 = i - 1;
|
||||
if (i <= 0)
|
||||
i = 1;
|
||||
}
|
||||
j = sizeof (__Long);
|
||||
for (ptr->_result_k = 0; (int) (sizeof (_Jv_Bigint) - sizeof (__Long)) + j <= i;
|
||||
j <<= 1)
|
||||
ptr->_result_k++;
|
||||
ptr->_result = Balloc (ptr, ptr->_result_k);
|
||||
s = s0 = (char *) ptr->_result;
|
||||
|
||||
if (ilim >= 0 && ilim <= Quick_max && try_quick)
|
||||
{
|
||||
/* Try to get by with floating-point arithmetic. */
|
||||
|
||||
i = 0;
|
||||
d2.d = d.d;
|
||||
k0 = k;
|
||||
ilim0 = ilim;
|
||||
ieps = 2; /* conservative */
|
||||
if (k > 0)
|
||||
{
|
||||
ds = tens[k & 0xf];
|
||||
j = k >> 4;
|
||||
if (j & Bletch)
|
||||
{
|
||||
/* prevent overflows */
|
||||
j &= Bletch - 1;
|
||||
d.d /= bigtens[n_bigtens - 1];
|
||||
ieps++;
|
||||
}
|
||||
for (; j; j >>= 1, i++)
|
||||
if (j & 1)
|
||||
{
|
||||
ieps++;
|
||||
ds *= bigtens[i];
|
||||
}
|
||||
d.d /= ds;
|
||||
}
|
||||
else if ((j1 = -k))
|
||||
{
|
||||
d.d *= tens[j1 & 0xf];
|
||||
for (j = j1 >> 4; j; j >>= 1, i++)
|
||||
if (j & 1)
|
||||
{
|
||||
ieps++;
|
||||
d.d *= bigtens[i];
|
||||
}
|
||||
}
|
||||
if (k_check && d.d < 1. && ilim > 0)
|
||||
{
|
||||
if (ilim1 <= 0)
|
||||
goto fast_failed;
|
||||
ilim = ilim1;
|
||||
k--;
|
||||
d.d *= 10.;
|
||||
ieps++;
|
||||
}
|
||||
eps.d = ieps * d.d + 7.;
|
||||
word0 (eps) -= (P - 1) * Exp_msk1;
|
||||
if (ilim == 0)
|
||||
{
|
||||
S = mhi = 0;
|
||||
d.d -= 5.;
|
||||
if (d.d > eps.d)
|
||||
goto one_digit;
|
||||
if (d.d < -eps.d)
|
||||
goto no_digits;
|
||||
goto fast_failed;
|
||||
}
|
||||
#ifndef No_leftright
|
||||
if (leftright)
|
||||
{
|
||||
/* Use Steele & White method of only
|
||||
* generating digits needed.
|
||||
*/
|
||||
eps.d = 0.5 / tens[ilim - 1] - eps.d;
|
||||
for (i = 0;;)
|
||||
{
|
||||
L = d.d;
|
||||
d.d -= L;
|
||||
*s++ = '0' + (int) L;
|
||||
if (d.d < eps.d)
|
||||
goto ret1;
|
||||
if (1. - d.d < eps.d)
|
||||
goto bump_up;
|
||||
if (++i >= ilim)
|
||||
break;
|
||||
eps.d *= 10.;
|
||||
d.d *= 10.;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
#endif
|
||||
/* Generate ilim digits, then fix them up. */
|
||||
eps.d *= tens[ilim - 1];
|
||||
for (i = 1;; i++, d.d *= 10.)
|
||||
{
|
||||
L = d.d;
|
||||
d.d -= L;
|
||||
*s++ = '0' + (int) L;
|
||||
if (i == ilim)
|
||||
{
|
||||
if (d.d > 0.5 + eps.d)
|
||||
goto bump_up;
|
||||
else if (d.d < 0.5 - eps.d)
|
||||
{
|
||||
while (*--s == '0');
|
||||
s++;
|
||||
goto ret1;
|
||||
}
|
||||
break;
|
||||
}
|
||||
}
|
||||
#ifndef No_leftright
|
||||
}
|
||||
#endif
|
||||
fast_failed:
|
||||
s = s0;
|
||||
d.d = d2.d;
|
||||
k = k0;
|
||||
ilim = ilim0;
|
||||
}
|
||||
|
||||
/* Do we have a "small" integer? */
|
||||
|
||||
if (be >= 0 && k <= Int_max)
|
||||
{
|
||||
/* Yes. */
|
||||
ds = tens[k];
|
||||
if (ndigits < 0 && ilim <= 0)
|
||||
{
|
||||
S = mhi = 0;
|
||||
if (ilim < 0 || d.d <= 5 * ds)
|
||||
goto no_digits;
|
||||
goto one_digit;
|
||||
}
|
||||
for (i = 1;; i++)
|
||||
{
|
||||
L = d.d / ds;
|
||||
d.d -= L * ds;
|
||||
#ifdef Check_FLT_ROUNDS
|
||||
/* If FLT_ROUNDS == 2, L will usually be high by 1 */
|
||||
if (d.d < 0)
|
||||
{
|
||||
L--;
|
||||
d.d += ds;
|
||||
}
|
||||
#endif
|
||||
*s++ = '0' + (int) L;
|
||||
if (i == ilim)
|
||||
{
|
||||
d.d += d.d;
|
||||
if (d.d > ds || (d.d == ds && L & 1))
|
||||
{
|
||||
bump_up:
|
||||
while (*--s == '9')
|
||||
if (s == s0)
|
||||
{
|
||||
k++;
|
||||
*s = '0';
|
||||
break;
|
||||
}
|
||||
++*s++;
|
||||
}
|
||||
break;
|
||||
}
|
||||
if (!(d.d *= 10.))
|
||||
break;
|
||||
}
|
||||
goto ret1;
|
||||
}
|
||||
|
||||
m2 = b2;
|
||||
m5 = b5;
|
||||
mhi = mlo = 0;
|
||||
if (leftright)
|
||||
{
|
||||
if (mode < 2)
|
||||
{
|
||||
i =
|
||||
#ifndef Sudden_Underflow
|
||||
denorm ? be + (Bias + (P - 1) - 1 + 1) :
|
||||
#endif
|
||||
#ifdef IBM
|
||||
1 + 4 * P - 3 - bbits + ((bbits + be - 1) & 3);
|
||||
#else
|
||||
1 + P - bbits;
|
||||
#endif
|
||||
}
|
||||
else
|
||||
{
|
||||
j = ilim - 1;
|
||||
if (m5 >= j)
|
||||
m5 -= j;
|
||||
else
|
||||
{
|
||||
s5 += j -= m5;
|
||||
b5 += j;
|
||||
m5 = 0;
|
||||
}
|
||||
if ((i = ilim) < 0)
|
||||
{
|
||||
m2 -= i;
|
||||
i = 0;
|
||||
}
|
||||
}
|
||||
b2 += i;
|
||||
s2 += i;
|
||||
mhi = i2b (ptr, 1);
|
||||
}
|
||||
if (m2 > 0 && s2 > 0)
|
||||
{
|
||||
i = m2 < s2 ? m2 : s2;
|
||||
b2 -= i;
|
||||
m2 -= i;
|
||||
s2 -= i;
|
||||
}
|
||||
if (b5 > 0)
|
||||
{
|
||||
if (leftright)
|
||||
{
|
||||
if (m5 > 0)
|
||||
{
|
||||
mhi = pow5mult (ptr, mhi, m5);
|
||||
b1 = mult (ptr, mhi, b);
|
||||
Bfree (ptr, b);
|
||||
b = b1;
|
||||
}
|
||||
if ((j = b5 - m5))
|
||||
b = pow5mult (ptr, b, j);
|
||||
}
|
||||
else
|
||||
b = pow5mult (ptr, b, b5);
|
||||
}
|
||||
S = i2b (ptr, 1);
|
||||
if (s5 > 0)
|
||||
S = pow5mult (ptr, S, s5);
|
||||
|
||||
/* Check for special case that d is a normalized power of 2. */
|
||||
|
||||
if (mode < 2)
|
||||
{
|
||||
if (!word1 (d) && !(word0 (d) & Bndry_mask)
|
||||
#ifndef Sudden_Underflow
|
||||
&& word0(d) & Exp_mask
|
||||
#endif
|
||||
)
|
||||
{
|
||||
/* The special case */
|
||||
b2 += Log2P;
|
||||
s2 += Log2P;
|
||||
spec_case = 1;
|
||||
}
|
||||
else
|
||||
spec_case = 0;
|
||||
}
|
||||
|
||||
/* Arrange for convenient computation of quotients:
|
||||
* shift left if necessary so divisor has 4 leading 0 bits.
|
||||
*
|
||||
* Perhaps we should just compute leading 28 bits of S once
|
||||
* and for all and pass them and a shift to quorem, so it
|
||||
* can do shifts and ors to compute the numerator for q.
|
||||
*/
|
||||
|
||||
#ifdef Pack_32
|
||||
if ((i = ((s5 ? 32 - hi0bits (S->_x[S->_wds - 1]) : 1) + s2) & 0x1f))
|
||||
i = 32 - i;
|
||||
#else
|
||||
if ((i = ((s5 ? 32 - hi0bits (S->_x[S->_wds - 1]) : 1) + s2) & 0xf))
|
||||
i = 16 - i;
|
||||
#endif
|
||||
if (i > 4)
|
||||
{
|
||||
i -= 4;
|
||||
b2 += i;
|
||||
m2 += i;
|
||||
s2 += i;
|
||||
}
|
||||
else if (i < 4)
|
||||
{
|
||||
i += 28;
|
||||
b2 += i;
|
||||
m2 += i;
|
||||
s2 += i;
|
||||
}
|
||||
if (b2 > 0)
|
||||
b = lshift (ptr, b, b2);
|
||||
if (s2 > 0)
|
||||
S = lshift (ptr, S, s2);
|
||||
if (k_check)
|
||||
{
|
||||
if (cmp (b, S) < 0)
|
||||
{
|
||||
k--;
|
||||
b = multadd (ptr, b, 10, 0); /* we botched the k estimate */
|
||||
if (leftright)
|
||||
mhi = multadd (ptr, mhi, 10, 0);
|
||||
ilim = ilim1;
|
||||
}
|
||||
}
|
||||
if (ilim <= 0 && mode > 2)
|
||||
{
|
||||
if (ilim < 0 || cmp (b, S = multadd (ptr, S, 5, 0)) <= 0)
|
||||
{
|
||||
/* no digits, fcvt style */
|
||||
no_digits:
|
||||
k = -1 - ndigits;
|
||||
goto ret;
|
||||
}
|
||||
one_digit:
|
||||
*s++ = '1';
|
||||
k++;
|
||||
goto ret;
|
||||
}
|
||||
if (leftright)
|
||||
{
|
||||
if (m2 > 0)
|
||||
mhi = lshift (ptr, mhi, m2);
|
||||
|
||||
/* Single precision case, */
|
||||
if (float_type)
|
||||
mhi = lshift (ptr, mhi, 29);
|
||||
|
||||
/* Compute mlo -- check for special case
|
||||
* that d is a normalized power of 2.
|
||||
*/
|
||||
|
||||
mlo = mhi;
|
||||
if (spec_case)
|
||||
{
|
||||
mhi = Balloc (ptr, mhi->_k);
|
||||
Bcopy (mhi, mlo);
|
||||
mhi = lshift (ptr, mhi, Log2P);
|
||||
}
|
||||
|
||||
for (i = 1;; i++)
|
||||
{
|
||||
dig = quorem (b, S) + '0';
|
||||
/* Do we yet have the shortest decimal string
|
||||
* that will round to d?
|
||||
*/
|
||||
j = cmp (b, mlo);
|
||||
delta = diff (ptr, S, mhi);
|
||||
j1 = delta->_sign ? 1 : cmp (b, delta);
|
||||
Bfree (ptr, delta);
|
||||
#ifndef ROUND_BIASED
|
||||
if (j1 == 0 && !mode && !(word1 (d) & 1))
|
||||
{
|
||||
if (dig == '9')
|
||||
goto round_9_up;
|
||||
if (j > 0)
|
||||
dig++;
|
||||
*s++ = dig;
|
||||
goto ret;
|
||||
}
|
||||
#endif
|
||||
if (j < 0 || (j == 0 && !mode
|
||||
#ifndef ROUND_BIASED
|
||||
&& !(word1 (d) & 1)
|
||||
#endif
|
||||
))
|
||||
{
|
||||
if (j1 > 0)
|
||||
{
|
||||
b = lshift (ptr, b, 1);
|
||||
j1 = cmp (b, S);
|
||||
if ((j1 > 0 || (j1 == 0 && dig & 1))
|
||||
&& dig++ == '9')
|
||||
goto round_9_up;
|
||||
}
|
||||
*s++ = dig;
|
||||
goto ret;
|
||||
}
|
||||
if (j1 > 0)
|
||||
{
|
||||
if (dig == '9')
|
||||
{ /* possible if i == 1 */
|
||||
round_9_up:
|
||||
*s++ = '9';
|
||||
goto roundoff;
|
||||
}
|
||||
*s++ = dig + 1;
|
||||
goto ret;
|
||||
}
|
||||
*s++ = dig;
|
||||
if (i == ilim)
|
||||
break;
|
||||
b = multadd (ptr, b, 10, 0);
|
||||
if (mlo == mhi)
|
||||
mlo = mhi = multadd (ptr, mhi, 10, 0);
|
||||
else
|
||||
{
|
||||
mlo = multadd (ptr, mlo, 10, 0);
|
||||
mhi = multadd (ptr, mhi, 10, 0);
|
||||
}
|
||||
}
|
||||
}
|
||||
else
|
||||
for (i = 1;; i++)
|
||||
{
|
||||
*s++ = dig = quorem (b, S) + '0';
|
||||
if (i >= ilim)
|
||||
break;
|
||||
b = multadd (ptr, b, 10, 0);
|
||||
}
|
||||
|
||||
/* Round off last digit */
|
||||
|
||||
b = lshift (ptr, b, 1);
|
||||
j = cmp (b, S);
|
||||
if (j > 0 || (j == 0))
|
||||
{
|
||||
roundoff:
|
||||
while (*--s == '9')
|
||||
if (s == s0)
|
||||
{
|
||||
k++;
|
||||
*s++ = '1';
|
||||
goto ret;
|
||||
}
|
||||
++*s++;
|
||||
}
|
||||
else
|
||||
{
|
||||
while (*--s == '0');
|
||||
s++;
|
||||
}
|
||||
ret:
|
||||
Bfree (ptr, S);
|
||||
if (mhi)
|
||||
{
|
||||
if (mlo && mlo != mhi)
|
||||
Bfree (ptr, mlo);
|
||||
Bfree (ptr, mhi);
|
||||
}
|
||||
ret1:
|
||||
Bfree (ptr, b);
|
||||
*s = 0;
|
||||
*decpt = k + 1;
|
||||
if (rve)
|
||||
*rve = s;
|
||||
return s0;
|
||||
}
|
||||
void free_Bigints(struct _Jv_Bigint *p)
|
||||
{
|
||||
struct _Jv_Bigint *l = p;
|
||||
while (l)
|
||||
{
|
||||
struct _Jv_Bigint *next = l->_next;
|
||||
free (l);
|
||||
l = next;
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
_VOID
|
||||
_DEFUN (_dtoa,
|
||||
(_d, mode, ndigits, decpt, sign, rve, buf, float_type),
|
||||
double _d _AND
|
||||
int mode _AND
|
||||
int ndigits _AND
|
||||
int *decpt _AND
|
||||
int *sign _AND
|
||||
char **rve _AND
|
||||
char **buf _AND
|
||||
int float_type)
|
||||
{
|
||||
struct _Jv_reent reent;
|
||||
char *p;
|
||||
int i, len;
|
||||
|
||||
memset (&reent, 0, sizeof reent);
|
||||
|
||||
p = _dtoa_r (&reent, _d, mode, ndigits, decpt, sign, rve, float_type);
|
||||
len = (int)(*rve - p);
|
||||
strncpy (*buf, p, len);
|
||||
/* map rve to same char in copied string */
|
||||
*rve = *buf + len;
|
||||
|
||||
for (i = 0; i < reent._max_k; ++i)
|
||||
{
|
||||
free_Bigints(reent._freelist[i]);
|
||||
}
|
||||
if (reent._freelist)
|
||||
free (reent._freelist);
|
||||
if (reent._result)
|
||||
free (reent._result);
|
||||
|
||||
free_Bigints(reent._p5s);
|
||||
|
||||
}
|
||||
@@ -0,0 +1,141 @@
|
||||
|
||||
/* @(#)e_fmod.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.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
/*
|
||||
* __ieee754_fmod(x,y)
|
||||
* Return x mod y in exact arithmetic
|
||||
* Method: shift and subtract
|
||||
*/
|
||||
|
||||
#include "fdlibm.h"
|
||||
|
||||
#ifndef _DOUBLE_IS_32BITS
|
||||
|
||||
#ifdef __STDC__
|
||||
static const double one = 1.0, Zero[] = {0.0, -0.0,};
|
||||
#else
|
||||
static double one = 1.0, Zero[] = {0.0, -0.0,};
|
||||
#endif
|
||||
|
||||
#ifdef __STDC__
|
||||
double __ieee754_fmod(double x, double y)
|
||||
#else
|
||||
double __ieee754_fmod(x,y)
|
||||
double x,y ;
|
||||
#endif
|
||||
{
|
||||
__int32_t n,hz,ix,iy,i;
|
||||
__uint32_t hx,hy,sx;
|
||||
__uint32_t lx,ly,lz;
|
||||
|
||||
EXTRACT_WORDS(hx,lx,x);
|
||||
EXTRACT_WORDS(hy,ly,y);
|
||||
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)) return x; /* |x|<|y| return x */
|
||||
if(lx==ly)
|
||||
return Zero[(__uint32_t)sx>>31]; /* |x|=|y| return x*0*/
|
||||
}
|
||||
|
||||
/* determine ix = ilogb(x) */
|
||||
if(hx<0x00100000) { /* subnormal x */
|
||||
if(hx==0) {
|
||||
for (ix = -1043, i=(__int32_t)lx; i>0; i<<=1) ix -=1;
|
||||
} else {
|
||||
for (ix = -1022,i=(__int32_t)(hx<<11); i>0; i<<=1) ix -=1;
|
||||
}
|
||||
} else ix = (__int32_t)(hx>>20)-1023;
|
||||
|
||||
/* determine iy = ilogb(y) */
|
||||
if(hy<0x00100000) { /* subnormal y */
|
||||
if(hy==0) {
|
||||
for (iy = -1043, i=(__int32_t)ly; i>0; i<<=1) iy -=1;
|
||||
} else {
|
||||
for (iy = -1022,i=(__int32_t)(hy<<11); i>0; i<<=1) iy -=1;
|
||||
}
|
||||
} else iy = (__int32_t)(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;
|
||||
while(n--) {
|
||||
hz=(__int32_t)(hx-hy);lz=lx-ly; if(lx<ly) hz -= 1;
|
||||
if(hz<0){hx = hx+hx+(lx>>31); lx = lx+lx;}
|
||||
else {
|
||||
if((hz|lz)==0) /* return sign(x)*0 */
|
||||
return Zero[(__uint32_t)sx>>31];
|
||||
hx = (__uint32_t)(hz+hz)+(lz>>31); lx = lz+lz;
|
||||
}
|
||||
}
|
||||
hz=(__int32_t)(hx-hy);lz=lx-ly; if(lx<ly) hz -= 1;
|
||||
if(hz>=0) {hx=(__uint32_t)hz;lx=lz;}
|
||||
|
||||
/* convert back to floating value and restore the sign */
|
||||
if((hx|lx)==0) /* return sign(x)*0 */
|
||||
return Zero[(__uint32_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));
|
||||
INSERT_WORDS(x,hx|sx,lx);
|
||||
} else { /* subnormal output */
|
||||
n = -1022 - iy;
|
||||
if(n<=20) {
|
||||
lx = (lx>>n)|((__uint32_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;
|
||||
}
|
||||
INSERT_WORDS(x,hx|sx,lx);
|
||||
x *= one; /* create necessary signal */
|
||||
}
|
||||
return x; /* exact output */
|
||||
}
|
||||
|
||||
#endif /* defined(_DOUBLE_IS_32BITS) */
|
||||
@@ -0,0 +1,312 @@
|
||||
|
||||
/* @(#)e_pow.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.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
/* __ieee754_pow(x,y) return x**y
|
||||
*
|
||||
* n
|
||||
* Method: Let x = 2 * (1+f)
|
||||
* 1. Compute and return log2(x) in two pieces:
|
||||
* log2(x) = w1 + w2,
|
||||
* where w1 has 53-24 = 29 bit trailing zeros.
|
||||
* 2. Perform y*log2(x) = n+y' by simulating muti-precision
|
||||
* arithmetic, where |y'|<=0.5.
|
||||
* 3. Return x**y = 2**n*exp(y'*log2)
|
||||
*
|
||||
* Special cases:
|
||||
* 1. (anything) ** 0 is 1
|
||||
* 2. (anything) ** 1 is itself
|
||||
* 3. (anything) ** NAN is NAN
|
||||
* 4. NAN ** (anything except 0) is NAN
|
||||
* 5. +-(|x| > 1) ** +INF is +INF
|
||||
* 6. +-(|x| > 1) ** -INF is +0
|
||||
* 7. +-(|x| < 1) ** +INF is +0
|
||||
* 8. +-(|x| < 1) ** -INF is +INF
|
||||
* 9. +-1 ** +-INF is NAN
|
||||
* 10. +0 ** (+anything except 0, NAN) is +0
|
||||
* 11. -0 ** (+anything except 0, NAN, odd integer) is +0
|
||||
* 12. +0 ** (-anything except 0, NAN) is +INF
|
||||
* 13. -0 ** (-anything except 0, NAN, odd integer) is +INF
|
||||
* 14. -0 ** (odd integer) = -( +0 ** (odd integer) )
|
||||
* 15. +INF ** (+anything except 0,NAN) is +INF
|
||||
* 16. +INF ** (-anything except 0,NAN) is +0
|
||||
* 17. -INF ** (anything) = -0 ** (-anything)
|
||||
* 18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
|
||||
* 19. (-anything except 0 and inf) ** (non-integer) is NAN
|
||||
*
|
||||
* Accuracy:
|
||||
* pow(x,y) returns x**y nearly rounded. In particular
|
||||
* pow(integer,integer)
|
||||
* always returns the correct integer provided it is
|
||||
* representable.
|
||||
*
|
||||
* Constants :
|
||||
* The hexadecimal values are the intended ones for the following
|
||||
* constants. The decimal values may be used, provided that the
|
||||
* compiler will convert from decimal to binary accurately enough
|
||||
* to produce the hexadecimal values shown.
|
||||
*/
|
||||
|
||||
#include "fdlibm.h"
|
||||
|
||||
#ifndef _DOUBLE_IS_32BITS
|
||||
|
||||
#ifdef __STDC__
|
||||
static const double
|
||||
#else
|
||||
static double
|
||||
#endif
|
||||
bp[] = {1.0, 1.5,},
|
||||
dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */
|
||||
dp_l[] = { 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */
|
||||
zero = 0.0,
|
||||
one = 1.0,
|
||||
two = 2.0,
|
||||
two53 = 9007199254740992.0, /* 0x43400000, 0x00000000 */
|
||||
huge = 1.0e300,
|
||||
tiny = 1.0e-300,
|
||||
/* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
|
||||
L1 = 5.99999999999994648725e-01, /* 0x3FE33333, 0x33333303 */
|
||||
L2 = 4.28571428578550184252e-01, /* 0x3FDB6DB6, 0xDB6FABFF */
|
||||
L3 = 3.33333329818377432918e-01, /* 0x3FD55555, 0x518F264D */
|
||||
L4 = 2.72728123808534006489e-01, /* 0x3FD17460, 0xA91D4101 */
|
||||
L5 = 2.30660745775561754067e-01, /* 0x3FCD864A, 0x93C9DB65 */
|
||||
L6 = 2.06975017800338417784e-01, /* 0x3FCA7E28, 0x4A454EEF */
|
||||
P1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */
|
||||
P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */
|
||||
P3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */
|
||||
P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */
|
||||
P5 = 4.13813679705723846039e-08, /* 0x3E663769, 0x72BEA4D0 */
|
||||
lg2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
|
||||
lg2_h = 6.93147182464599609375e-01, /* 0x3FE62E43, 0x00000000 */
|
||||
lg2_l = -1.90465429995776804525e-09, /* 0xBE205C61, 0x0CA86C39 */
|
||||
ovt = 8.0085662595372944372e-0017, /* -(1024-log2(ovfl+.5ulp)) */
|
||||
cp = 9.61796693925975554329e-01, /* 0x3FEEC709, 0xDC3A03FD =2/(3ln2) */
|
||||
cp_h = 9.61796700954437255859e-01, /* 0x3FEEC709, 0xE0000000 =(float)cp */
|
||||
cp_l = -7.02846165095275826516e-09, /* 0xBE3E2FE0, 0x145B01F5 =tail of cp_h*/
|
||||
ivln2 = 1.44269504088896338700e+00, /* 0x3FF71547, 0x652B82FE =1/ln2 */
|
||||
ivln2_h = 1.44269502162933349609e+00, /* 0x3FF71547, 0x60000000 =24b 1/ln2*/
|
||||
ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
|
||||
|
||||
#ifdef __STDC__
|
||||
double __ieee754_pow(double x, double y)
|
||||
#else
|
||||
double __ieee754_pow(x,y)
|
||||
double x, y;
|
||||
#endif
|
||||
{
|
||||
double z,ax,z_h,z_l,p_h,p_l;
|
||||
double y1,t1,t2,r,s,t,u,v,w;
|
||||
__int32_t n;
|
||||
__uint32_t ix,iy,i,j,k,yisint;
|
||||
__uint32_t lx,ly,hx,hy;
|
||||
|
||||
EXTRACT_WORDS(hx,lx,x);
|
||||
EXTRACT_WORDS(hy,ly,y);
|
||||
ix = hx&0x7fffffff; iy = hy&0x7fffffff;
|
||||
|
||||
/* y==zero: x**0 = 1 */
|
||||
if((iy|ly)==0) return one;
|
||||
|
||||
/* +-NaN return x+y */
|
||||
if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) ||
|
||||
iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0)))
|
||||
return x+y;
|
||||
|
||||
/* determine if y is an odd int when x < 0
|
||||
* yisint = 0 ... y is not an integer
|
||||
* yisint = 1 ... y is an odd int
|
||||
* yisint = 2 ... y is an even int
|
||||
*/
|
||||
yisint = 0;
|
||||
if((__int32_t)hx<0) {
|
||||
if(iy>=0x43400000) yisint = 2; /* even integer y */
|
||||
else if(iy>=0x3ff00000) {
|
||||
k = (iy>>20)-0x3ff; /* exponent */
|
||||
if(k>20) {
|
||||
j = ly>>(52-k);
|
||||
if((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);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* special value of y */
|
||||
if(ly==0) {
|
||||
if (iy==0x7ff00000) { /* y is +-inf */
|
||||
if(((ix-0x3ff00000)|lx)==0)
|
||||
return y - y; /* inf**+-1 is NaN */
|
||||
else if (ix >= 0x3ff00000)/* (|x|>1)**+-inf = inf,0 */
|
||||
return ((__int32_t)hy>=0)? y: zero;
|
||||
else /* (|x|<1)**-,+inf = inf,0 */
|
||||
return ((__int32_t)hy<0)?-y: zero;
|
||||
}
|
||||
if(iy==0x3ff00000) { /* y is +-1 */
|
||||
if((__int32_t)hy<0) return one/x; else return x;
|
||||
}
|
||||
if(hy==0x40000000) return x*x; /* y is 2 */
|
||||
if(hy==0x3fe00000) { /* y is 0.5 */
|
||||
if((__int32_t)hx>=0) /* x >= +0 */
|
||||
return __ieee754_sqrt(x);
|
||||
}
|
||||
}
|
||||
|
||||
ax = fabs(x);
|
||||
/* special value of x */
|
||||
if(lx==0) {
|
||||
if(ix==0x7ff00000||ix==0||ix==0x3ff00000){
|
||||
z = ax; /*x is +-0,+-inf,+-1*/
|
||||
if((__int32_t)hy<0) z = one/z; /* z = (1/|x|) */
|
||||
if((__int32_t)hx<0) {
|
||||
if(((ix-0x3ff00000)|yisint)==0) {
|
||||
z = (z-z)/(z-z); /* (-1)**non-int is NaN */
|
||||
} else if(yisint==1)
|
||||
z = -z; /* (x<0)**odd = -(|x|**odd) */
|
||||
}
|
||||
return z;
|
||||
}
|
||||
}
|
||||
|
||||
/* (x<0)**(non-int) is NaN */
|
||||
/* CYGNUS LOCAL: This used to be
|
||||
if((((hx>>31)+1)|yisint)==0) return (x-x)/(x-x);
|
||||
but ANSI C says a right shift of a signed negative quantity is
|
||||
implementation defined. */
|
||||
if(((((__uint32_t)hx>>31)-1)|yisint)==0) return (x-x)/(x-x);
|
||||
|
||||
/* |y| is huge */
|
||||
if(iy>0x41e00000) { /* if |y| > 2**31 */
|
||||
if(iy>0x43f00000){ /* if |y| > 2**64, must o/uflow */
|
||||
if(ix<=0x3fefffff) return ((__int32_t)hy<0)? huge*huge:tiny*tiny;
|
||||
if(ix>=0x3ff00000) return ((__int32_t)hy>0)? huge*huge:tiny*tiny;
|
||||
}
|
||||
/* over/underflow if x is not close to one */
|
||||
if(ix<0x3fefffff) return ((__int32_t)hy<0)? huge*huge:tiny*tiny;
|
||||
if(ix>0x3ff00000) return ((__int32_t)hy>0)? huge*huge:tiny*tiny;
|
||||
/* now |1-x| is tiny <= 2**-20, suffice to compute
|
||||
log(x) by x-x^2/2+x^3/3-x^4/4 */
|
||||
t = x-1; /* t has 20 trailing zeros */
|
||||
w = (t*t)*(0.5-t*(0.3333333333333333333333-t*0.25));
|
||||
u = ivln2_h*t; /* ivln2_h has 21 sig. bits */
|
||||
v = t*ivln2_l-w*ivln2;
|
||||
t1 = u+v;
|
||||
SET_LOW_WORD(t1,0);
|
||||
t2 = v-(t1-u);
|
||||
} else {
|
||||
double s2,s_h,s_l,t_h,t_l;
|
||||
n = 0;
|
||||
/* take care subnormal number */
|
||||
if(ix<0x00100000)
|
||||
{ax *= two53; n -= 53; GET_HIGH_WORD(ix,ax); }
|
||||
n += (__int32_t)((ix)>>20)-0x3ff;
|
||||
j = ix&0x000fffff;
|
||||
/* determine interval */
|
||||
ix = j|0x3ff00000; /* normalize ix */
|
||||
if(j<=0x3988E) k=0; /* |x|<sqrt(3/2) */
|
||||
else if(j<0xBB67A) k=1; /* |x|<sqrt(3) */
|
||||
else {k=0;n+=1;ix -= 0x00100000;}
|
||||
SET_HIGH_WORD(ax,ix);
|
||||
|
||||
/* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
|
||||
u = ax-bp[k]; /* bp[0]=1.0, bp[1]=1.5 */
|
||||
v = one/(ax+bp[k]);
|
||||
s = u*v;
|
||||
s_h = s;
|
||||
SET_LOW_WORD(s_h,0);
|
||||
/* t_h=ax+bp[k] High */
|
||||
t_h = zero;
|
||||
SET_HIGH_WORD(t_h,((ix>>1)|0x20000000)+0x00080000+(k<<18));
|
||||
t_l = ax - (t_h-bp[k]);
|
||||
s_l = v*((u-s_h*t_h)-s_h*t_l);
|
||||
/* compute log(ax) */
|
||||
s2 = s*s;
|
||||
r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6)))));
|
||||
r += s_l*(s_h+s);
|
||||
s2 = s_h*s_h;
|
||||
t_h = 3.0+s2+r;
|
||||
SET_LOW_WORD(t_h,0);
|
||||
t_l = r-((t_h-3.0)-s2);
|
||||
/* u+v = s*(1+...) */
|
||||
u = s_h*t_h;
|
||||
v = s_l*t_h+t_l*s;
|
||||
/* 2/(3log2)*(s+...) */
|
||||
p_h = u+v;
|
||||
SET_LOW_WORD(p_h,0);
|
||||
p_l = v-(p_h-u);
|
||||
z_h = cp_h*p_h; /* cp_h+cp_l = 2/(3*log2) */
|
||||
z_l = cp_l*p_h+p_l*cp+dp_l[k];
|
||||
/* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
|
||||
t = (double)n;
|
||||
t1 = (((z_h+z_l)+dp_h[k])+t);
|
||||
SET_LOW_WORD(t1,0);
|
||||
t2 = z_l-(((t1-t)-dp_h[k])-z_h);
|
||||
}
|
||||
|
||||
s = one; /* s (sign of result -ve**odd) = -1 else = 1 */
|
||||
if(((((__uint32_t)hx>>31)-1)|(yisint-1))==0)
|
||||
s = -one;/* (-ve)**(odd int) */
|
||||
|
||||
/* 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;
|
||||
z = p_l+p_h;
|
||||
EXTRACT_WORDS(j,i,z);
|
||||
if (j>=0x40900000) { /* z >= 1024 */
|
||||
if(((j-0x40900000)|i)!=0) /* if z > 1024 */
|
||||
return s*huge*huge; /* overflow */
|
||||
else {
|
||||
if(p_l+ovt>z-p_h) return s*huge*huge; /* overflow */
|
||||
}
|
||||
} else if((j&0x7fffffff)>=0x4090cc00 ) { /* z <= -1075 */
|
||||
if(((j-0xc090cc00)|i)!=0) /* z < -1075 */
|
||||
return s*tiny*tiny; /* underflow */
|
||||
else {
|
||||
if(p_l<=z-p_h) return s*tiny*tiny; /* underflow */
|
||||
}
|
||||
}
|
||||
/*
|
||||
* compute 2**(p_h+p_l)
|
||||
*/
|
||||
i = j&0x7fffffff;
|
||||
k = (i>>20)-0x3ff;
|
||||
n = 0;
|
||||
if(i>0x3fe00000) { /* if |z| > 0.5, set n = [z+0.5] */
|
||||
n = (__int32_t)j+(__int32_t)(0x00100000>>(k+1));
|
||||
k = (__uint32_t)((n&0x7fffffff)>>20)-0x3ff; /* new k for n */
|
||||
t = zero;
|
||||
SET_HIGH_WORD(t,(__uint32_t)(n&~(0x000fffff>>k)));
|
||||
n = ((n&0x000fffff)|0x00100000)>>(20-k);
|
||||
if((__int32_t)j<0) n = -n;
|
||||
p_h -= t;
|
||||
}
|
||||
t = p_l+p_h;
|
||||
SET_LOW_WORD(t,0);
|
||||
u = t*lg2_h;
|
||||
v = (p_l-(t-p_h))*lg2+t*lg2_l;
|
||||
z = u+v;
|
||||
w = v-(z-u);
|
||||
t = z*z;
|
||||
t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
|
||||
r = (z*t1)/(t1-two)-(w+z*w);
|
||||
z = one-(r-z);
|
||||
GET_HIGH_WORD(j,z);
|
||||
j += (__uint32_t)(n<<20);
|
||||
if((j>>20)<=0) z = scalbn(z,(int)n); /* subnormal output */
|
||||
else SET_HIGH_WORD(z,j);
|
||||
return s*z;
|
||||
}
|
||||
|
||||
#endif /* defined(_DOUBLE_IS_32BITS) */
|
||||
@@ -0,0 +1,80 @@
|
||||
|
||||
/* @(#)e_remainder.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.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
/* __ieee754_remainder(x,p)
|
||||
* Return :
|
||||
* returns x REM p = x - [x/p]*p as if in infinite
|
||||
* precise arithmetic, where [x/p] is the (infinite bit)
|
||||
* integer nearest x/p (in half way case choose the even one).
|
||||
* Method :
|
||||
* Based on fmod() return x-[x/p]chopped*p exactlp.
|
||||
*/
|
||||
|
||||
#include "fdlibm.h"
|
||||
|
||||
#ifndef _DOUBLE_IS_32BITS
|
||||
|
||||
#ifdef __STDC__
|
||||
static const double zero = 0.0;
|
||||
#else
|
||||
static double zero = 0.0;
|
||||
#endif
|
||||
|
||||
|
||||
#ifdef __STDC__
|
||||
double __ieee754_remainder(double x, double p)
|
||||
#else
|
||||
double __ieee754_remainder(x,p)
|
||||
double x,p;
|
||||
#endif
|
||||
{
|
||||
__uint32_t hx,hp;
|
||||
__uint32_t sx,lx,lp;
|
||||
double p_half;
|
||||
|
||||
EXTRACT_WORDS(hx,lx,x);
|
||||
EXTRACT_WORDS(hp,lp,p);
|
||||
sx = hx&0x80000000;
|
||||
hp &= 0x7fffffff;
|
||||
hx &= 0x7fffffff;
|
||||
|
||||
/* purge off exception values */
|
||||
if((hp|lp)==0) return (x*p)/(x*p); /* p = 0 */
|
||||
if((hx>=0x7ff00000)|| /* x not finite */
|
||||
((hp>=0x7ff00000)&& /* p is NaN */
|
||||
(((hp-0x7ff00000)|lp)!=0)))
|
||||
return (x*p)/(x*p);
|
||||
|
||||
|
||||
if (hp<=0x7fdfffff) x = __ieee754_fmod(x,p+p); /* now x < 2p */
|
||||
if (((hx-hp)|(lx-lp))==0) return zero*x;
|
||||
x = fabs(x);
|
||||
p = fabs(p);
|
||||
if (hp<0x00200000) {
|
||||
if(x+x>p) {
|
||||
x-=p;
|
||||
if(x+x>=p) x -= p;
|
||||
}
|
||||
} else {
|
||||
p_half = 0.5*p;
|
||||
if(x>p_half) {
|
||||
x-=p;
|
||||
if(x>=p_half) x -= p;
|
||||
}
|
||||
}
|
||||
GET_HIGH_WORD(hx,x);
|
||||
SET_HIGH_WORD(x,hx^sx);
|
||||
return x;
|
||||
}
|
||||
|
||||
#endif /* defined(_DOUBLE_IS_32BITS) */
|
||||
@@ -0,0 +1,456 @@
|
||||
|
||||
/* @(#)e_sqrt.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.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
/* __ieee754_sqrt(x)
|
||||
* Return correctly rounded sqrt.
|
||||
* ------------------------------------------
|
||||
* | Use the hardware sqrt if you have one |
|
||||
* ------------------------------------------
|
||||
* Method:
|
||||
* Bit by bit method using integer arithmetic. (Slow, but portable)
|
||||
* 1. Normalization
|
||||
* Scale x to y in [1,4) with even powers of 2:
|
||||
* find an integer k such that 1 <= (y=x*2^(2k)) < 4, then
|
||||
* sqrt(x) = 2^k * sqrt(y)
|
||||
* 2. Bit by bit computation
|
||||
* Let q = sqrt(y) truncated to i bit after binary point (q = 1),
|
||||
* i 0
|
||||
* i+1 2
|
||||
* s = 2*q , and y = 2 * ( y - q ). (1)
|
||||
* i i i i
|
||||
*
|
||||
* To compute q from q , one checks whether
|
||||
* i+1 i
|
||||
*
|
||||
* -(i+1) 2
|
||||
* (q + 2 ) <= y. (2)
|
||||
* i
|
||||
* -(i+1)
|
||||
* If (2) is false, then q = q ; otherwise q = q + 2 .
|
||||
* i+1 i i+1 i
|
||||
*
|
||||
* With some algebric manipulation, it is not difficult to see
|
||||
* that (2) is equivalent to
|
||||
* -(i+1)
|
||||
* s + 2 <= y (3)
|
||||
* i i
|
||||
*
|
||||
* The advantage of (3) is that s and y can be computed by
|
||||
* i i
|
||||
* the following recurrence formula:
|
||||
* if (3) is false
|
||||
*
|
||||
* s = s , y = y ; (4)
|
||||
* i+1 i i+1 i
|
||||
*
|
||||
* otherwise,
|
||||
* -i -(i+1)
|
||||
* s = s + 2 , y = y - s - 2 (5)
|
||||
* i+1 i i+1 i i
|
||||
*
|
||||
* One may easily use induction to prove (4) and (5).
|
||||
* Note. Since the left hand side of (3) contain only i+2 bits,
|
||||
* it does not necessary to do a full (53-bit) comparison
|
||||
* in (3).
|
||||
* 3. Final rounding
|
||||
* After generating the 53 bits result, we compute one more bit.
|
||||
* Together with the remainder, we can decide whether the
|
||||
* result is exact, bigger than 1/2ulp, or less than 1/2ulp
|
||||
* (it will never equal to 1/2ulp).
|
||||
* The rounding mode can be detected by checking whether
|
||||
* huge + tiny is equal to huge, and whether huge - tiny is
|
||||
* equal to huge for some floating point number "huge" and "tiny".
|
||||
*
|
||||
* Special cases:
|
||||
* sqrt(+-0) = +-0 ... exact
|
||||
* sqrt(inf) = inf
|
||||
* sqrt(-ve) = NaN ... with invalid signal
|
||||
* sqrt(NaN) = NaN ... with invalid signal for signaling NaN
|
||||
*
|
||||
* Other methods : see the appended file at the end of the program below.
|
||||
*---------------
|
||||
*/
|
||||
|
||||
#include "fdlibm.h"
|
||||
|
||||
#ifndef _DOUBLE_IS_32BITS
|
||||
|
||||
#ifdef __STDC__
|
||||
static const double one = 1.0, tiny=1.0e-300;
|
||||
#else
|
||||
static double one = 1.0, tiny=1.0e-300;
|
||||
#endif
|
||||
|
||||
#ifdef __STDC__
|
||||
double __ieee754_sqrt(double x)
|
||||
#else
|
||||
double __ieee754_sqrt(x)
|
||||
double x;
|
||||
#endif
|
||||
{
|
||||
double z;
|
||||
__int32_t sign = (int)0x80000000;
|
||||
__uint32_t r,t1,s1,ix1,q1;
|
||||
__uint32_t ix0,s0,q,m,t,i;
|
||||
|
||||
EXTRACT_WORDS(ix0,ix1,x);
|
||||
|
||||
/* take care of Inf and NaN */
|
||||
if((ix0&0x7ff00000)==0x7ff00000) {
|
||||
return x*x+x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf
|
||||
sqrt(-inf)=sNaN */
|
||||
}
|
||||
/* take care of zero */
|
||||
if(ix0<=0) {
|
||||
if(((ix0&(~sign))|ix1)==0) return x;/* sqrt(+-0) = +-0 */
|
||||
else {
|
||||
__int32_t tmp;
|
||||
tmp = (__int32_t)ix0;
|
||||
if(tmp<0)
|
||||
return (x-x)/(x-x); /* sqrt(-ve) = sNaN */
|
||||
}
|
||||
}
|
||||
/* normalize x */
|
||||
m = (ix0>>20);
|
||||
if(m==0) { /* subnormal x */
|
||||
while(ix0==0) {
|
||||
m -= 21;
|
||||
ix0 |= (ix1>>11); ix1 <<= 21;
|
||||
}
|
||||
for(i=0;(ix0&0x00100000)==0;i++) ix0<<=1;
|
||||
m -= i-1;
|
||||
ix0 |= (ix1>>(32-i));
|
||||
ix1 <<= i;
|
||||
}
|
||||
m -= 1023; /* unbias exponent */
|
||||
ix0 = (ix0&0x000fffff)|0x00100000;
|
||||
if(m&1){ /* odd m, double x to make it even */
|
||||
ix0 += ix0 + ((ix1&sign)>>31);
|
||||
ix1 += ix1;
|
||||
}
|
||||
m >>= 1; /* m = [m/2] */
|
||||
|
||||
/* generate sqrt(x) bit by bit */
|
||||
ix0 += ix0 + ((ix1&sign)>>31);
|
||||
ix1 += ix1;
|
||||
q = q1 = s0 = s1 = 0; /* [q,q1] = sqrt(x) */
|
||||
r = 0x00200000; /* r = moving bit from right to left */
|
||||
|
||||
while(r!=0) {
|
||||
t = s0+r;
|
||||
if(t<=ix0) {
|
||||
s0 = t+r;
|
||||
ix0 -= t;
|
||||
q += r;
|
||||
}
|
||||
ix0 += ix0 + ((ix1&sign)>>31);
|
||||
ix1 += ix1;
|
||||
r>>=1;
|
||||
}
|
||||
|
||||
r = (__uint32_t)sign;
|
||||
while(r!=0) {
|
||||
t1 = s1+r;
|
||||
t = s0;
|
||||
if((t<ix0)||((t==ix0)&&(t1<=ix1))) {
|
||||
s1 = t1+r;
|
||||
if(((t1&sign)==(__uint32_t)sign)&&(s1&sign)==0) s0 += 1;
|
||||
ix0 -= t;
|
||||
if (ix1 < t1) ix0 -= 1;
|
||||
ix1 -= t1;
|
||||
q1 += r;
|
||||
}
|
||||
ix0 += ix0 + ((ix1&sign)>>31);
|
||||
ix1 += ix1;
|
||||
r>>=1;
|
||||
}
|
||||
|
||||
/* use floating add to find out rounding direction */
|
||||
if((ix0|ix1)!=0) {
|
||||
z = one-tiny; /* trigger inexact flag */
|
||||
if (z>=one) {
|
||||
z = one+tiny;
|
||||
if (q1==(__uint32_t)0xffffffff) { q1=0; q += 1;}
|
||||
else if (z>one) {
|
||||
if (q1==(__uint32_t)0xfffffffe) q+=1;
|
||||
q1+=2;
|
||||
} else
|
||||
q1 += (q1&1);
|
||||
}
|
||||
}
|
||||
ix0 = (q>>1)+0x3fe00000;
|
||||
ix1 = q1>>1;
|
||||
if ((q&1)==1) ix1 |= sign;
|
||||
ix0 += (m <<20);
|
||||
INSERT_WORDS(z,ix0,ix1);
|
||||
return z;
|
||||
}
|
||||
|
||||
#endif /* defined(_DOUBLE_IS_32BITS) */
|
||||
|
||||
/*
|
||||
Other methods (use floating-point arithmetic)
|
||||
-------------
|
||||
(This is a copy of a drafted paper by Prof W. Kahan
|
||||
and K.C. Ng, written in May, 1986)
|
||||
|
||||
Two algorithms are given here to implement sqrt(x)
|
||||
(IEEE double precision arithmetic) in software.
|
||||
Both supply sqrt(x) correctly rounded. The first algorithm (in
|
||||
Section A) uses newton iterations and involves four divisions.
|
||||
The second one uses reciproot iterations to avoid division, but
|
||||
requires more multiplications. Both algorithms need the ability
|
||||
to chop results of arithmetic operations instead of round them,
|
||||
and the INEXACT flag to indicate when an arithmetic operation
|
||||
is executed exactly with no roundoff error, all part of the
|
||||
standard (IEEE 754-1985). The ability to perform shift, add,
|
||||
subtract and logical AND operations upon 32-bit words is needed
|
||||
too, though not part of the standard.
|
||||
|
||||
A. sqrt(x) by Newton Iteration
|
||||
|
||||
(1) Initial approximation
|
||||
|
||||
Let x0 and x1 be the leading and the trailing 32-bit words of
|
||||
a floating point number x (in IEEE double format) respectively
|
||||
|
||||
1 11 52 ...widths
|
||||
------------------------------------------------------
|
||||
x: |s| e | f |
|
||||
------------------------------------------------------
|
||||
msb lsb msb lsb ...order
|
||||
|
||||
|
||||
------------------------ ------------------------
|
||||
x0: |s| e | f1 | x1: | f2 |
|
||||
------------------------ ------------------------
|
||||
|
||||
By performing shifts and subtracts on x0 and x1 (both regarded
|
||||
as integers), we obtain an 8-bit approximation of sqrt(x) as
|
||||
follows.
|
||||
|
||||
k := (x0>>1) + 0x1ff80000;
|
||||
y0 := k - T1[31&(k>>15)]. ... y ~ sqrt(x) to 8 bits
|
||||
Here k is a 32-bit integer and T1[] is an integer array containing
|
||||
correction terms. Now magically the floating value of y (y's
|
||||
leading 32-bit word is y0, the value of its trailing word is 0)
|
||||
approximates sqrt(x) to almost 8-bit.
|
||||
|
||||
Value of T1:
|
||||
static int T1[32]= {
|
||||
0, 1024, 3062, 5746, 9193, 13348, 18162, 23592,
|
||||
29598, 36145, 43202, 50740, 58733, 67158, 75992, 85215,
|
||||
83599, 71378, 60428, 50647, 41945, 34246, 27478, 21581,
|
||||
16499, 12183, 8588, 5674, 3403, 1742, 661, 130,};
|
||||
|
||||
(2) Iterative refinement
|
||||
|
||||
Apply Heron's rule three times to y, we have y approximates
|
||||
sqrt(x) to within 1 ulp (Unit in the Last Place):
|
||||
|
||||
y := (y+x/y)/2 ... almost 17 sig. bits
|
||||
y := (y+x/y)/2 ... almost 35 sig. bits
|
||||
y := y-(y-x/y)/2 ... within 1 ulp
|
||||
|
||||
|
||||
Remark 1.
|
||||
Another way to improve y to within 1 ulp is:
|
||||
|
||||
y := (y+x/y) ... almost 17 sig. bits to 2*sqrt(x)
|
||||
y := y - 0x00100006 ... almost 18 sig. bits to sqrt(x)
|
||||
|
||||
2
|
||||
(x-y )*y
|
||||
y := y + 2* ---------- ...within 1 ulp
|
||||
2
|
||||
3y + x
|
||||
|
||||
|
||||
This formula has one division fewer than the one above; however,
|
||||
it requires more multiplications and additions. Also x must be
|
||||
scaled in advance to avoid spurious overflow in evaluating the
|
||||
expression 3y*y+x. Hence it is not recommended uless division
|
||||
is slow. If division is very slow, then one should use the
|
||||
reciproot algorithm given in section B.
|
||||
|
||||
(3) Final adjustment
|
||||
|
||||
By twiddling y's last bit it is possible to force y to be
|
||||
correctly rounded according to the prevailing rounding mode
|
||||
as follows. Let r and i be copies of the rounding mode and
|
||||
inexact flag before entering the square root program. Also we
|
||||
use the expression y+-ulp for the next representable floating
|
||||
numbers (up and down) of y. Note that y+-ulp = either fixed
|
||||
point y+-1, or multiply y by nextafter(1,+-inf) in chopped
|
||||
mode.
|
||||
|
||||
I := FALSE; ... reset INEXACT flag I
|
||||
R := RZ; ... set rounding mode to round-toward-zero
|
||||
z := x/y; ... chopped quotient, possibly inexact
|
||||
If(not I) then { ... if the quotient is exact
|
||||
if(z=y) {
|
||||
I := i; ... restore inexact flag
|
||||
R := r; ... restore rounded mode
|
||||
return sqrt(x):=y.
|
||||
} else {
|
||||
z := z - ulp; ... special rounding
|
||||
}
|
||||
}
|
||||
i := TRUE; ... sqrt(x) is inexact
|
||||
If (r=RN) then z=z+ulp ... rounded-to-nearest
|
||||
If (r=RP) then { ... round-toward-+inf
|
||||
y = y+ulp; z=z+ulp;
|
||||
}
|
||||
y := y+z; ... chopped sum
|
||||
y0:=y0-0x00100000; ... y := y/2 is correctly rounded.
|
||||
I := i; ... restore inexact flag
|
||||
R := r; ... restore rounded mode
|
||||
return sqrt(x):=y.
|
||||
|
||||
(4) Special cases
|
||||
|
||||
Square root of +inf, +-0, or NaN is itself;
|
||||
Square root of a negative number is NaN with invalid signal.
|
||||
|
||||
|
||||
B. sqrt(x) by Reciproot Iteration
|
||||
|
||||
(1) Initial approximation
|
||||
|
||||
Let x0 and x1 be the leading and the trailing 32-bit words of
|
||||
a floating point number x (in IEEE double format) respectively
|
||||
(see section A). By performing shifs and subtracts on x0 and y0,
|
||||
we obtain a 7.8-bit approximation of 1/sqrt(x) as follows.
|
||||
|
||||
k := 0x5fe80000 - (x0>>1);
|
||||
y0:= k - T2[63&(k>>14)]. ... y ~ 1/sqrt(x) to 7.8 bits
|
||||
|
||||
Here k is a 32-bit integer and T2[] is an integer array
|
||||
containing correction terms. Now magically the floating
|
||||
value of y (y's leading 32-bit word is y0, the value of
|
||||
its trailing word y1 is set to zero) approximates 1/sqrt(x)
|
||||
to almost 7.8-bit.
|
||||
|
||||
Value of T2:
|
||||
static int T2[64]= {
|
||||
0x1500, 0x2ef8, 0x4d67, 0x6b02, 0x87be, 0xa395, 0xbe7a, 0xd866,
|
||||
0xf14a, 0x1091b,0x11fcd,0x13552,0x14999,0x15c98,0x16e34,0x17e5f,
|
||||
0x18d03,0x19a01,0x1a545,0x1ae8a,0x1b5c4,0x1bb01,0x1bfde,0x1c28d,
|
||||
0x1c2de,0x1c0db,0x1ba73,0x1b11c,0x1a4b5,0x1953d,0x18266,0x16be0,
|
||||
0x1683e,0x179d8,0x18a4d,0x19992,0x1a789,0x1b445,0x1bf61,0x1c989,
|
||||
0x1d16d,0x1d77b,0x1dddf,0x1e2ad,0x1e5bf,0x1e6e8,0x1e654,0x1e3cd,
|
||||
0x1df2a,0x1d635,0x1cb16,0x1be2c,0x1ae4e,0x19bde,0x1868e,0x16e2e,
|
||||
0x1527f,0x1334a,0x11051,0xe951, 0xbe01, 0x8e0d, 0x5924, 0x1edd,};
|
||||
|
||||
(2) Iterative refinement
|
||||
|
||||
Apply Reciproot iteration three times to y and multiply the
|
||||
result by x to get an approximation z that matches sqrt(x)
|
||||
to about 1 ulp. To be exact, we will have
|
||||
-1ulp < sqrt(x)-z<1.0625ulp.
|
||||
|
||||
... set rounding mode to Round-to-nearest
|
||||
y := y*(1.5-0.5*x*y*y) ... almost 15 sig. bits to 1/sqrt(x)
|
||||
y := y*((1.5-2^-30)+0.5*x*y*y)... about 29 sig. bits to 1/sqrt(x)
|
||||
... special arrangement for better accuracy
|
||||
z := x*y ... 29 bits to sqrt(x), with z*y<1
|
||||
z := z + 0.5*z*(1-z*y) ... about 1 ulp to sqrt(x)
|
||||
|
||||
Remark 2. The constant 1.5-2^-30 is chosen to bias the error so that
|
||||
(a) the term z*y in the final iteration is always less than 1;
|
||||
(b) the error in the final result is biased upward so that
|
||||
-1 ulp < sqrt(x) - z < 1.0625 ulp
|
||||
instead of |sqrt(x)-z|<1.03125ulp.
|
||||
|
||||
(3) Final adjustment
|
||||
|
||||
By twiddling y's last bit it is possible to force y to be
|
||||
correctly rounded according to the prevailing rounding mode
|
||||
as follows. Let r and i be copies of the rounding mode and
|
||||
inexact flag before entering the square root program. Also we
|
||||
use the expression y+-ulp for the next representable floating
|
||||
numbers (up and down) of y. Note that y+-ulp = either fixed
|
||||
point y+-1, or multiply y by nextafter(1,+-inf) in chopped
|
||||
mode.
|
||||
|
||||
R := RZ; ... set rounding mode to round-toward-zero
|
||||
switch(r) {
|
||||
case RN: ... round-to-nearest
|
||||
if(x<= z*(z-ulp)...chopped) z = z - ulp; else
|
||||
if(x<= z*(z+ulp)...chopped) z = z; else z = z+ulp;
|
||||
break;
|
||||
case RZ:case RM: ... round-to-zero or round-to--inf
|
||||
R:=RP; ... reset rounding mod to round-to-+inf
|
||||
if(x<z*z ... rounded up) z = z - ulp; else
|
||||
if(x>=(z+ulp)*(z+ulp) ...rounded up) z = z+ulp;
|
||||
break;
|
||||
case RP: ... round-to-+inf
|
||||
if(x>(z+ulp)*(z+ulp)...chopped) z = z+2*ulp; else
|
||||
if(x>z*z ...chopped) z = z+ulp;
|
||||
break;
|
||||
}
|
||||
|
||||
Remark 3. The above comparisons can be done in fixed point. For
|
||||
example, to compare x and w=z*z chopped, it suffices to compare
|
||||
x1 and w1 (the trailing parts of x and w), regarding them as
|
||||
two's complement integers.
|
||||
|
||||
...Is z an exact square root?
|
||||
To determine whether z is an exact square root of x, let z1 be the
|
||||
trailing part of z, and also let x0 and x1 be the leading and
|
||||
trailing parts of x.
|
||||
|
||||
If ((z1&0x03ffffff)!=0) ... not exact if trailing 26 bits of z!=0
|
||||
I := 1; ... Raise Inexact flag: z is not exact
|
||||
else {
|
||||
j := 1 - [(x0>>20)&1] ... j = logb(x) mod 2
|
||||
k := z1 >> 26; ... get z's 25-th and 26-th
|
||||
fraction bits
|
||||
I := i or (k&j) or ((k&(j+j+1))!=(x1&3));
|
||||
}
|
||||
R:= r ... restore rounded mode
|
||||
return sqrt(x):=z.
|
||||
|
||||
If multiplication is cheaper then the foregoing red tape, the
|
||||
Inexact flag can be evaluated by
|
||||
|
||||
I := i;
|
||||
I := (z*z!=x) or I.
|
||||
|
||||
Note that z*z can overwrite I; this value must be sensed if it is
|
||||
True.
|
||||
|
||||
Remark 4. If z*z = x exactly, then bit 25 to bit 0 of z1 must be
|
||||
zero.
|
||||
|
||||
--------------------
|
||||
z1: | f2 |
|
||||
--------------------
|
||||
bit 31 bit 0
|
||||
|
||||
Further more, bit 27 and 26 of z1, bit 0 and 1 of x1, and the odd
|
||||
or even of logb(x) have the following relations:
|
||||
|
||||
-------------------------------------------------
|
||||
bit 27,26 of z1 bit 1,0 of x1 logb(x)
|
||||
-------------------------------------------------
|
||||
00 00 odd and even
|
||||
01 01 even
|
||||
10 10 odd
|
||||
10 00 even
|
||||
11 01 even
|
||||
-------------------------------------------------
|
||||
|
||||
(4) Special cases (see (4) of Section A).
|
||||
|
||||
*/
|
||||
@@ -0,0 +1,270 @@
|
||||
|
||||
/* @(#)fdlibm.h 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.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
/* CYGNUS LOCAL: Include files. */
|
||||
#include <math.h>
|
||||
#include "ieeefp.h"
|
||||
|
||||
/* CYGNUS LOCAL: Default to XOPEN_MODE. */
|
||||
#define _XOPEN_MODE
|
||||
|
||||
#ifndef __P
|
||||
#ifdef __STDC__
|
||||
#define __P(p) p
|
||||
#else
|
||||
#define __P(p) ()
|
||||
#endif
|
||||
#endif
|
||||
|
||||
#define HUGE ((float)3.40282346638528860e+38)
|
||||
|
||||
/*
|
||||
* set X_TLOSS = pi*2**52, which is possibly defined in <values.h>
|
||||
* (one may replace the following line by "#include <values.h>")
|
||||
*/
|
||||
|
||||
#define X_TLOSS 1.41484755040568800000e+16
|
||||
|
||||
/* glue typedefs for VxWorks */
|
||||
|
||||
typedef int __int32_t;
|
||||
typedef unsigned int __uint32_t;
|
||||
|
||||
/* Functions that are not documented, and are not in <math.h>. */
|
||||
|
||||
extern double logb __P((double));
|
||||
#ifdef _SCALB_INT
|
||||
extern double scalb __P((double, int));
|
||||
#else
|
||||
extern double scalb __P((double, double));
|
||||
#endif
|
||||
extern double scalbn __P((double x, int n));
|
||||
extern double significand __P((double));
|
||||
extern double copysign __P((double x, double y));
|
||||
|
||||
/* ieee style elementary functions */
|
||||
extern double __ieee754_sqrt __P((double));
|
||||
extern double __ieee754_acos __P((double));
|
||||
extern double __ieee754_acosh __P((double));
|
||||
extern double __ieee754_log __P((double));
|
||||
extern double __ieee754_atanh __P((double));
|
||||
extern double __ieee754_asin __P((double));
|
||||
extern double __ieee754_atan2 __P((double,double));
|
||||
extern double __ieee754_exp __P((double));
|
||||
extern double __ieee754_cosh __P((double));
|
||||
extern double __ieee754_fmod __P((double,double));
|
||||
extern double __ieee754_pow __P((double,double));
|
||||
extern double __ieee754_lgamma_r __P((double,int *));
|
||||
extern double __ieee754_gamma_r __P((double,int *));
|
||||
extern double __ieee754_log10 __P((double));
|
||||
extern double __ieee754_sinh __P((double));
|
||||
extern double __ieee754_hypot __P((double,double));
|
||||
extern double __ieee754_j0 __P((double));
|
||||
extern double __ieee754_j1 __P((double));
|
||||
extern double __ieee754_y0 __P((double));
|
||||
extern double __ieee754_y1 __P((double));
|
||||
extern double __ieee754_jn __P((int,double));
|
||||
extern double __ieee754_yn __P((int,double));
|
||||
extern double __ieee754_remainder __P((double,double));
|
||||
extern __int32_t __ieee754_rem_pio2 __P((double,double*));
|
||||
#ifdef _SCALB_INT
|
||||
extern double __ieee754_scalb __P((double,int));
|
||||
#else
|
||||
extern double __ieee754_scalb __P((double,double));
|
||||
#endif
|
||||
|
||||
/* fdlibm kernel function */
|
||||
extern double __kernel_standard __P((double,double,int));
|
||||
extern double __kernel_sin __P((double,double,int));
|
||||
extern double __kernel_cos __P((double,double));
|
||||
extern double __kernel_tan __P((double,double,int));
|
||||
extern int __kernel_rem_pio2 __P((double*,double*,int,int,int,const __int32_t*));
|
||||
|
||||
/* Undocumented float functions. */
|
||||
extern float logbf __P((float));
|
||||
#ifdef _SCALB_INT
|
||||
extern float scalbf __P((float, int));
|
||||
#else
|
||||
extern float scalbf __P((float, float));
|
||||
#endif
|
||||
extern float significandf __P((float));
|
||||
|
||||
/* ieee style elementary float functions */
|
||||
extern float __ieee754_sqrtf __P((float));
|
||||
extern float __ieee754_acosf __P((float));
|
||||
extern float __ieee754_acoshf __P((float));
|
||||
extern float __ieee754_logf __P((float));
|
||||
extern float __ieee754_atanhf __P((float));
|
||||
extern float __ieee754_asinf __P((float));
|
||||
extern float __ieee754_atan2f __P((float,float));
|
||||
extern float __ieee754_expf __P((float));
|
||||
extern float __ieee754_coshf __P((float));
|
||||
extern float __ieee754_fmodf __P((float,float));
|
||||
extern float __ieee754_powf __P((float,float));
|
||||
extern float __ieee754_lgammaf_r __P((float,int *));
|
||||
extern float __ieee754_gammaf_r __P((float,int *));
|
||||
extern float __ieee754_log10f __P((float));
|
||||
extern float __ieee754_sinhf __P((float));
|
||||
extern float __ieee754_hypotf __P((float,float));
|
||||
extern float __ieee754_j0f __P((float));
|
||||
extern float __ieee754_j1f __P((float));
|
||||
extern float __ieee754_y0f __P((float));
|
||||
extern float __ieee754_y1f __P((float));
|
||||
extern float __ieee754_jnf __P((int,float));
|
||||
extern float __ieee754_ynf __P((int,float));
|
||||
extern float __ieee754_remainderf __P((float,float));
|
||||
extern __int32_t __ieee754_rem_pio2f __P((float,float*));
|
||||
#ifdef _SCALB_INT
|
||||
extern float __ieee754_scalbf __P((float,int));
|
||||
#else
|
||||
extern float __ieee754_scalbf __P((float,float));
|
||||
#endif
|
||||
|
||||
/* float versions of fdlibm kernel functions */
|
||||
extern float __kernel_sinf __P((float,float,int));
|
||||
extern float __kernel_cosf __P((float,float));
|
||||
extern float __kernel_tanf __P((float,float,int));
|
||||
extern int __kernel_rem_pio2f __P((float*,float*,int,int,int,const __int32_t*));
|
||||
|
||||
/* The original code used statements like
|
||||
n0 = ((*(int*)&one)>>29)^1; * index of high word *
|
||||
ix0 = *(n0+(int*)&x); * high word of x *
|
||||
ix1 = *((1-n0)+(int*)&x); * low word of x *
|
||||
to dig two 32 bit words out of the 64 bit IEEE floating point
|
||||
value. That is non-ANSI, and, moreover, the gcc instruction
|
||||
scheduler gets it wrong. We instead use the following macros.
|
||||
Unlike the original code, we determine the endianness at compile
|
||||
time, not at run time; I don't see much benefit to selecting
|
||||
endianness at run time. */
|
||||
|
||||
#ifndef __IEEE_BIG_ENDIAN
|
||||
#ifndef __IEEE_LITTLE_ENDIAN
|
||||
#error Must define endianness
|
||||
#endif
|
||||
#endif
|
||||
|
||||
/* A union which permits us to convert between a double and two 32 bit
|
||||
ints. */
|
||||
|
||||
#ifdef __IEEE_BIG_ENDIAN
|
||||
|
||||
typedef union
|
||||
{
|
||||
double value;
|
||||
struct
|
||||
{
|
||||
__uint32_t msw;
|
||||
__uint32_t lsw;
|
||||
} parts;
|
||||
} ieee_double_shape_type;
|
||||
|
||||
#endif
|
||||
|
||||
#ifdef __IEEE_LITTLE_ENDIAN
|
||||
|
||||
typedef union
|
||||
{
|
||||
double value;
|
||||
struct
|
||||
{
|
||||
__uint32_t lsw;
|
||||
__uint32_t msw;
|
||||
} parts;
|
||||
} ieee_double_shape_type;
|
||||
|
||||
#endif
|
||||
|
||||
/* Get two 32 bit ints from a double. */
|
||||
|
||||
#define EXTRACT_WORDS(ix0,ix1,d) \
|
||||
do { \
|
||||
ieee_double_shape_type ew_u; \
|
||||
ew_u.value = (d); \
|
||||
(ix0) = ew_u.parts.msw; \
|
||||
(ix1) = ew_u.parts.lsw; \
|
||||
} while (0)
|
||||
|
||||
/* Get the more significant 32 bit int from a double. */
|
||||
|
||||
#define GET_HIGH_WORD(i,d) \
|
||||
do { \
|
||||
ieee_double_shape_type gh_u; \
|
||||
gh_u.value = (d); \
|
||||
(i) = gh_u.parts.msw; \
|
||||
} while (0)
|
||||
|
||||
/* Get the less significant 32 bit int from a double. */
|
||||
|
||||
#define GET_LOW_WORD(i,d) \
|
||||
do { \
|
||||
ieee_double_shape_type gl_u; \
|
||||
gl_u.value = (d); \
|
||||
(i) = gl_u.parts.lsw; \
|
||||
} while (0)
|
||||
|
||||
/* Set a double from two 32 bit ints. */
|
||||
|
||||
#define INSERT_WORDS(d,ix0,ix1) \
|
||||
do { \
|
||||
ieee_double_shape_type iw_u; \
|
||||
iw_u.parts.msw = (ix0); \
|
||||
iw_u.parts.lsw = (ix1); \
|
||||
(d) = iw_u.value; \
|
||||
} while (0)
|
||||
|
||||
/* Set the more significant 32 bits of a double from an int. */
|
||||
|
||||
#define SET_HIGH_WORD(d,v) \
|
||||
do { \
|
||||
ieee_double_shape_type sh_u; \
|
||||
sh_u.value = (d); \
|
||||
sh_u.parts.msw = (v); \
|
||||
(d) = sh_u.value; \
|
||||
} while (0)
|
||||
|
||||
/* Set the less significant 32 bits of a double from an int. */
|
||||
|
||||
#define SET_LOW_WORD(d,v) \
|
||||
do { \
|
||||
ieee_double_shape_type sl_u; \
|
||||
sl_u.value = (d); \
|
||||
sl_u.parts.lsw = (v); \
|
||||
(d) = sl_u.value; \
|
||||
} while (0)
|
||||
|
||||
/* A union which permits us to convert between a float and a 32 bit
|
||||
int. */
|
||||
|
||||
typedef union
|
||||
{
|
||||
float value;
|
||||
__uint32_t word;
|
||||
} ieee_float_shape_type;
|
||||
|
||||
/* Get a 32 bit int from a float. */
|
||||
|
||||
#define GET_FLOAT_WORD(i,d) \
|
||||
do { \
|
||||
ieee_float_shape_type gf_u; \
|
||||
gf_u.value = (d); \
|
||||
(i) = gf_u.word; \
|
||||
} while (0)
|
||||
|
||||
/* Set a float from a 32 bit int. */
|
||||
|
||||
#define SET_FLOAT_WORD(d,i) \
|
||||
do { \
|
||||
ieee_float_shape_type sf_u; \
|
||||
sf_u.word = (i); \
|
||||
(d) = sf_u.value; \
|
||||
} while (0)
|
||||
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
@@ -0,0 +1,44 @@
|
||||
/* floatFuncBind.c - float point io management function binding library */
|
||||
|
||||
/*
|
||||
* Copyright (c) 2010 Wind River Systems, Inc.
|
||||
*
|
||||
* The right to copy, distribute, modify or otherwise make use
|
||||
* of this software may be licensed only pursuant to the terms
|
||||
* of an applicable Wind River license agreement.
|
||||
*/
|
||||
|
||||
/*
|
||||
modification history
|
||||
--------------------
|
||||
|
||||
*/
|
||||
|
||||
/* Includes */
|
||||
|
||||
#include <vxWorks.h>
|
||||
#include <fioLib.h>
|
||||
#include <stdarg.h>
|
||||
|
||||
int (* _func_fioFltFormatRtn)
|
||||
(
|
||||
va_list * pVaList, /* vararg list */
|
||||
int precision, /* precision */
|
||||
BOOL doAlt, /* doAlt boolean */
|
||||
int fmtSym, /* format symbol */
|
||||
BOOL * pDoSign, /* where to fill in doSign result */
|
||||
char * pBufStart, /* buffer start */
|
||||
char * pBufEnd /* buffer end */
|
||||
) = NULL;
|
||||
|
||||
|
||||
BOOL (* _func_fioFltScanRtn)
|
||||
(
|
||||
FAST int * pReturn,
|
||||
int returnSpec, /* 0, 'l', 'L' */
|
||||
int fieldwidth,
|
||||
FAST GET_NEXT_CHAR_FUNCPTR getRtn,
|
||||
FAST void *getArg,
|
||||
int * pCh,
|
||||
int * pNchars
|
||||
) = NULL;
|
||||
@@ -0,0 +1,727 @@
|
||||
/* floatLib.c - floating-point formatting and scanning library */
|
||||
|
||||
/*
|
||||
* Copyright (c) 1991-1994, 1996, 1999, 2001, 2005, 2009-2010, 2014, 2016, 2018, 2021 Wind River Systems, Inc.
|
||||
*
|
||||
* The right to copy, distribute or otherwise make use of this software
|
||||
* may be licensed only pursuant to the terms of an applicable Wind River
|
||||
* license agreement.
|
||||
*/
|
||||
|
||||
/*
|
||||
modification history
|
||||
--------------------
|
||||
02t,14dec21,ume Refactor the cvt for prec (VXW6-87362)
|
||||
02s,15ian18,mst Fix for VXW6-86417
|
||||
02r,25may16,mst Fixes for VXW6-85514/VXW6-85490/VXW6-85496
|
||||
(floatlib/dtoa 64-bit support)
|
||||
02q,03mar16,mst Fixes for VXW6-85306 and VXW6-85240.
|
||||
02p,07sep14,mcm Updated fix for VXW6-81137; fixed PAD loop.
|
||||
02o,07aug14,mcm Fix for VXW6-81137; floating point processing for printf has
|
||||
inaccuracies; using BSD functionality
|
||||
02n,17may14,mcm Fix for VXW6-82645; making sscanf process floats like strtod
|
||||
done.
|
||||
02m,06sep10,fan removed funcBindP.h
|
||||
02l,16apr10,mcm Fixing LP64 floating point scanf and a warning.
|
||||
02k,02mar10,mcm Fix for WIND00188900. Also removed warnings (1. removed
|
||||
ununsed floatFpTypeGet declaration and 2. added parantheses
|
||||
to a logical expression in function cvt).
|
||||
02j,23sep09,fan _WRS_KEEP_USING_FUNCPTR in floatLib.c removed
|
||||
02i,03jul09,fan updated for LP64 data model
|
||||
02h,27sep05,kk moved fioFltInstall from fioLib.c to here
|
||||
02g,13may05,gvm Fix for SPR109279: adding BSD copyright notice in-line
|
||||
02f,19oct01,dcb Fix the routine title line to match the coding conventions.
|
||||
02e,14mar99,jdi doc: removed refs to config.h and/or configAll.h (SPR 25663).
|
||||
02d,30jul96,dbt fixed floatScan to avoid loose of precision (SPR #3829).
|
||||
Updated copyright.
|
||||
02c,01dec94,ism fixed alt 'g' format bug. Added in-line docs.
|
||||
02b,06oct94,ism fixed %.XXg outputting more than XX sig figs (SPR#3695)
|
||||
02a,27jan93,jdi documentation cleanup for 5.1; fixed SPR 1399.
|
||||
01z,30jul92,kdl Restored 01w changes (clobbered by 01x); removed fpTypeGet()
|
||||
and isInf, isNan, & isZero macros.
|
||||
01y,30jul92,rrr Removed decl of DOUBLE (now in mathLibP.h)
|
||||
01x,30jul92,smb (Accidental checkin.)
|
||||
01w,30jul92,kdl Significantly reworked for new ansi fio; moved frexp(),
|
||||
+jcf ldexp(), modf() to src/libc/math; moved DOUBLE union
|
||||
definition to mathP.h; removed special handling for i960.
|
||||
01v,18jul92,smb Changed errno.h to errnoLib.h.
|
||||
01u,13jul92,smb frexp, ldexp, modf now in libc for the SPARC (temporary change)
|
||||
01t,26may92,rrr the tree shuffle
|
||||
01s,04mar92,wmd placed declaration of atof() outside of 960 conditional,
|
||||
included string.h.
|
||||
01r,10dec91,gae added includes for ANSI.
|
||||
01q,28oct91,wmd added in type casting for *pCh in floatScan as per Intel mods.
|
||||
01p,04oct91,rrr passed through the ansification filter
|
||||
-changed functions to ansi style
|
||||
-changed includes to have absolute path from h/
|
||||
-fixed #else and #endif
|
||||
-changed VOID to void
|
||||
-changed copyright notice
|
||||
01o,09jun91,del integration of mods by intel for interfacing to gnu960
|
||||
libraries and fix bugs from 01m.
|
||||
01n,05apr91,jdi documentation -- removed header parens and x-ref numbers;
|
||||
doc review by dnw.
|
||||
01m,25mar91,del interfaces to gnu960 flt. pt. libraries for I960CA
|
||||
01l,07feb91,jaa documentation cleanup.
|
||||
01k,23aug90,elr Documentation
|
||||
01j,01aug90,dnw fixed bug in floatScan scanning <nn>e<nn> format numbers.
|
||||
changed floatFormat to take pVaList instead of vaList.
|
||||
01i,04jun90,dnw changed floatFormat to take a vararg va_list instead of a
|
||||
double arg, so that fioLib can be free of all flt pt refs.
|
||||
replaced floatAtoF with new routine floatScan, written
|
||||
from scratch, with our own algorithm for scanning flt pt
|
||||
numbers, designed to work with scanField() in fioLib.
|
||||
fixed bug incorrectly scanning large numbers that are too big.
|
||||
fixed documentation.
|
||||
make ecvtb, fcvtb, gcvtb be "no manual".
|
||||
01h,07mar90,jdi documentation cleanup.
|
||||
01g,20apr89,dab fixed precision bug in cvtb().
|
||||
documentation touchup in fcvtb() and cvtb().
|
||||
01f,05jun88,dnw changed name from fltLib to floatLib.
|
||||
01e,30may88,dnw changed to v4 names.
|
||||
01d,28may88,dnw changed atof() to fltAtoF().
|
||||
01c,17mar88,gae added 'E' and 'G' format specifiers.
|
||||
01b,05nov87,jlf documentation
|
||||
01a,01aug87,gae written/extracted from fioLib.c.
|
||||
*/
|
||||
|
||||
/*
|
||||
DESCRIPTION
|
||||
This library provides the floating-point I/O formatting and scanning
|
||||
support routines.
|
||||
|
||||
The floating-point formatting and scanning support routines are not
|
||||
directly callable; they are connected to call-outs in the printf()/scanf()
|
||||
family of functions in fioLib. This is done dynamically by the routine
|
||||
floatInit(), which is called by the root task, usrRoot(), in usrConfig.c
|
||||
when the configuration macro INCLUDE_FLOATING_POINT is defined.
|
||||
If this option is omitted (i.e., floatInit() is not called), floating-point
|
||||
format specifications in printf() and sscanf() are not supported.
|
||||
|
||||
INCLUDE FILES: math.h
|
||||
|
||||
SEE ALSO: fioLib
|
||||
*/
|
||||
|
||||
#include <vxWorks.h>
|
||||
#include <vwModNum.h>
|
||||
#include <fioLib.h>
|
||||
#include <string.h>
|
||||
#include <stdlib.h>
|
||||
#include <ctype.h>
|
||||
#include <string.h>
|
||||
#include <errnoLib.h>
|
||||
#include <math.h>
|
||||
#include <float.h>
|
||||
#include <limits.h>
|
||||
#include <stdarg.h>
|
||||
#include <fioLib.h>
|
||||
#include <private/floatioP.h>
|
||||
#include <private/mathP.h>
|
||||
|
||||
/* macros */
|
||||
|
||||
#define to_digit(c) ((c) - '0')
|
||||
#define is_digit(c) ((unsigned)to_digit(c) <= 9)
|
||||
#define to_char(n) ((n) + '0')
|
||||
|
||||
#define LARGE_EXP (INT_MAX/100)
|
||||
|
||||
#define F0 1.25
|
||||
#define F1 F0*F0
|
||||
#define F2 F1*F1
|
||||
#define F3 F2*F2
|
||||
#define F4 F3*F3
|
||||
#define F5 F4*F4
|
||||
#define F6 F5*F5
|
||||
#define F7 F6*F6
|
||||
#define F8 F7*F7
|
||||
#define F9 F8*F8
|
||||
#define F10 F9*F9
|
||||
#define F11 F10*F10
|
||||
#define F12 F11*F11
|
||||
#define F13 F12*F12
|
||||
#define F14 F13*F13
|
||||
|
||||
static const double f125[] = {
|
||||
F0, F1, F2, F3, F4, F5, F6, F7, F8,
|
||||
F9, 1.0
|
||||
};
|
||||
|
||||
#define F125MAX (&f125[sizeof(f125)/sizeof(double) - 1])
|
||||
|
||||
typedef int (*FLOAT_FUNCPTR) (void *);
|
||||
|
||||
/* forward static functions */
|
||||
|
||||
LOCAL int floatFormat (va_list *pVaList, int precision, BOOL doAlt,
|
||||
int fmtSym, BOOL *pDoSign, char *pBufStart,
|
||||
char *pBufEnd);
|
||||
LOCAL BOOL floatScan (int *pReturn, int returnSpec, int fieldwidth,
|
||||
FLOAT_FUNCPTR getRtn, void *getArg, int *pCh, int *pNchars);
|
||||
LOCAL int cvt (double number, int prec, BOOL doAlt, int fmtch,
|
||||
BOOL *pDoSign, char *startp, char *endp);
|
||||
LOCAL char * exponentCvt (char *p, int exp, int fmtch);
|
||||
|
||||
/*******************************************************************************
|
||||
*
|
||||
* fioFltInstall - install routines to print and scan floating-point numbers
|
||||
*
|
||||
* This routine is a hook for floatLib to install its special floating-point
|
||||
* formatting & scanning routines. It should only be called by floatInit().
|
||||
*
|
||||
* NOMANUAL
|
||||
*/
|
||||
|
||||
void fioFltInstall
|
||||
(
|
||||
FUNCPTR formatRtn, /* routine to format floats for output */
|
||||
FUNCPTR scanRtn /* routine to scan input for floats */
|
||||
)
|
||||
{
|
||||
_func_fioFltFormatRtn = formatRtn;
|
||||
_func_fioFltScanRtn = scanRtn;
|
||||
}
|
||||
|
||||
/*******************************************************************************
|
||||
*
|
||||
* floatInit - initialize floating-point I/O support
|
||||
*
|
||||
* This routine must be called if floating-point format specifications are
|
||||
* to be supported by the printf()/scanf() family of routines.
|
||||
* If the configuration macro INCLUDE_FLOATING_POINT is defined, it is called
|
||||
* by the root task, usrRoot(), in usrConfig.c or in prjConfig.c for project
|
||||
* builds.
|
||||
*
|
||||
* RETURNS: N/A
|
||||
*/
|
||||
void floatInit (void)
|
||||
{
|
||||
fioFltInstall ((FLOAT_FUNCPTR)floatFormat, (FLOAT_FUNCPTR)floatScan);
|
||||
}
|
||||
|
||||
/*******************************************************************************
|
||||
*
|
||||
* floatFormat - format arg for output
|
||||
*
|
||||
* This routine converts from a floating-point value to an ASCII representation.
|
||||
* The <fmtSym> parameter indicates the format type desired; the actual work
|
||||
* is done in cvt.
|
||||
*
|
||||
* RETURNS: number of characters placed in buffer
|
||||
*/
|
||||
|
||||
/* extern void setAlt(int); // no effect. Need to set ALT=1 in flags directly. */
|
||||
/* extern void setFpt(int); // no effect and superflous. Would need to set FPT=0x100 directly, but not checked */
|
||||
extern void bsd_cvt(double value, int ndigits, int flags, char *sign, int *decpt, int ch, int *ndig, char **buf);
|
||||
|
||||
LOCAL int floatFormat
|
||||
(
|
||||
va_list * pVaList, /* vararg list */
|
||||
int precision, /* precision */
|
||||
BOOL doAlt, /* doAlt boolean */
|
||||
int fmtSym, /* format symbol */
|
||||
BOOL * pDoSign, /* where to fill in doSign result */
|
||||
char * pBufStart, /* buffer start */
|
||||
char * pBufEnd /* buffer end */
|
||||
)
|
||||
{
|
||||
/* signbit() not available in VxWorks kernel.
|
||||
* gcc 4+ has a polymorphic __builtin_signbit(), but here not performance
|
||||
* critical and all supported CPUs are IEEE 754 compliant,
|
||||
* so accessing the bit directly sufficient and generic enough.
|
||||
*/
|
||||
union s_u { unsigned long long n; double r; } u;
|
||||
u.r = (double) va_arg (*pVaList, double);
|
||||
|
||||
*pDoSign = FALSE; /* assume no sign needed */
|
||||
|
||||
if (isInf(u.r)) /* infinite? */
|
||||
{
|
||||
strcpy (pBufStart, "Inf"); /* fill in the string */
|
||||
if (u.n >> 63) /* less than 0.0, (u.r < 0.0 would also work) */
|
||||
*pDoSign = TRUE; /* we need a sign */
|
||||
return (-3); /* length will be three */
|
||||
}
|
||||
|
||||
if (isNan(u.r)) /* not a number? */
|
||||
{
|
||||
strcpy (pBufStart, "NaN"); /* fill in the string */
|
||||
#undef NEGATIVE_NAN
|
||||
#ifdef NEGATIVE_NAN /* For compatibility to existing test case results maybe don't be too accurate */
|
||||
if (u.n >> 63) /* less than 0.0 (double comparison does not work with NaN */
|
||||
*pDoSign = TRUE; /* we need a sign also for NaN */
|
||||
#endif
|
||||
return (-3); /* length will be three */
|
||||
}
|
||||
|
||||
/* pDoSign must(!) be pre-initialized to FALSE on entering cvt */
|
||||
|
||||
return (cvt (u.r, precision, doAlt, fmtSym, pDoSign, pBufStart, pBufEnd));
|
||||
}
|
||||
|
||||
#define PAD(to, n, p) { \
|
||||
int padN = (n); \
|
||||
while (padN-- > 0) \
|
||||
*to++ = p; \
|
||||
}
|
||||
|
||||
#define PRINT(to, from, num) { \
|
||||
char *from_tmp = from; \
|
||||
int i; \
|
||||
for (i = num; i > 0; i--) \
|
||||
*to++ = *from_tmp++; \
|
||||
}
|
||||
|
||||
#define ALT 0x001 /* alternate form */
|
||||
#define FPT 0x100 /* Floating point number */
|
||||
|
||||
/******************************************************************************
|
||||
*
|
||||
* cvt - helper function for floatFormat
|
||||
*
|
||||
* RETURNS:
|
||||
*/
|
||||
|
||||
LOCAL int cvt
|
||||
(
|
||||
double number,
|
||||
FAST int prec,
|
||||
BOOL doAlt,
|
||||
int fmtch,
|
||||
BOOL * pDoSign,
|
||||
char * startp,
|
||||
char * endp
|
||||
)
|
||||
{
|
||||
char * t;
|
||||
int flags = FPT;
|
||||
/* fioFormatV, the only caller of _func_fioFltFormatRtn, on "%e" or "%g" clips prec to MAXFRACT=39 sufficient for 128bit mantissa,
|
||||
* even on attempting like printf("%.1000g",v). (for comparision: Linux clips prec to 54.)
|
||||
* However e.g. by printf("%300f" a buffer size of more as 300 is needed. startp buffer with BUF=400 is big enough.
|
||||
* Here and inside bsd_cvt/__dtoa there is no protection agains buffer overwrite.
|
||||
* If this would occure, malloc would just replace the stack corruption by a heap+stack corruption.
|
||||
* Normally only little buffer is used. Avoid time overhead and managment of tiny short living malloc.
|
||||
* char *retBuf = malloc((int)(endp - startp)); //by fioFormatV, startp size is BUF=400, but MAXFRACT=39 is enough.
|
||||
* Don't but a second large buffer on stack. char retBufMem[400]; char *retBuf = retBufMem; // Must be a pointer
|
||||
* Share retBuf with startp buffer. Then at front leave space for usual prefix like "0."
|
||||
* If more needed, need to copy retBuf from back to front to final position before filling earlier parts.
|
||||
* See at decpt <= 0) below.
|
||||
*/
|
||||
char *retBuf = startp+3;
|
||||
int ndig = 0;
|
||||
int nexp = 0;
|
||||
int decpt = 0;
|
||||
char expstr[7];
|
||||
char softSign;
|
||||
|
||||
if(doAlt) {
|
||||
flags |= ALT;
|
||||
}
|
||||
|
||||
bsd_cvt(number, prec, flags, &softSign, &decpt, fmtch, &ndig, &retBuf);
|
||||
|
||||
/* The version below was that full of bugs that almost nothing of the import is left over */
|
||||
|
||||
/****************************************************************
|
||||
*
|
||||
* The author of this software is David M. Gay.
|
||||
*
|
||||
* Copyright (c) 1991, 2006 by AT&T.
|
||||
*
|
||||
* Permission to use, copy, modify, and distribute this software for any
|
||||
* purpose without fee is hereby granted, provided that this entire notice
|
||||
* is included in all copies of any software which is or includes a copy
|
||||
* or modification of this software and in all copies of the supporting
|
||||
* documentation for such software.
|
||||
*
|
||||
* THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
|
||||
* WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY
|
||||
* REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
|
||||
* OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
|
||||
*
|
||||
***************************************************************/
|
||||
|
||||
if(softSign) {
|
||||
*pDoSign = TRUE;
|
||||
}
|
||||
|
||||
if(fmtch == 'g' || fmtch == 'G') {
|
||||
if(prec == 0) {
|
||||
prec = 1;
|
||||
}
|
||||
if(decpt <= -4 || decpt > prec) {
|
||||
fmtch = (fmtch == 'g') ? 'e' : 'E';
|
||||
}
|
||||
}
|
||||
|
||||
if(fmtch == 'e' || fmtch == 'E') {
|
||||
char *ret;
|
||||
char *p = expstr;
|
||||
--decpt;
|
||||
ret = exponentCvt(p, decpt, fmtch);
|
||||
nexp = ret - p;
|
||||
}
|
||||
|
||||
fmtch = tolower(fmtch);
|
||||
/* ++ legacy to allow caller to prepending sign, auto skipped by caller if not touched (stays EOS) */
|
||||
t = ++startp; /* get an extra slot for rounding. */
|
||||
if(fmtch >= 'f') { /* 'f' or 'g' */
|
||||
if(fabs(number) == 0.0) { /* On 0.0 1st digit is 0, only one digit and ndig undefined */
|
||||
PRINT(t, "0", 1);
|
||||
/* Special case: on 0.0 even without doAlt fill trailing 0 - if(doAlt) { ... } */
|
||||
if (fmtch == 'f' || doAlt) {
|
||||
if (prec) {
|
||||
PRINT(t, ".", 1);
|
||||
if (fmtch == 'f') {
|
||||
PAD(t, prec, '0');
|
||||
} else {
|
||||
PAD(t, prec - 1, '0');
|
||||
}
|
||||
}
|
||||
}
|
||||
} else
|
||||
if(decpt <= 0) {
|
||||
/* On not using overlapped in/out buffer could replace by
|
||||
* PRINT(t, "0", 1); PRINT(t, ".", 1); PAD(t, -decpt, '0'); PRINT(t, retBuf, ndig);
|
||||
* //if (doAlt) { PAD(t, prec - ndig + decpt, '0'); } // superflous if flags|=ALT properly set on bsd_cvt call
|
||||
*/
|
||||
int i; retBuf+=ndig;
|
||||
|
||||
if (prec)
|
||||
t += 2+ndig-decpt;
|
||||
else
|
||||
t += 1+ndig-decpt;
|
||||
|
||||
for (i = ndig; i > 0; i--) { *--t = *--retBuf; }
|
||||
for (i = -decpt; i>0; i--) { *--t = '0'; }
|
||||
|
||||
if (prec)
|
||||
*--t = '.';
|
||||
|
||||
*--t = '0';
|
||||
if (prec)
|
||||
t += 2+ndig-decpt;
|
||||
else
|
||||
t += 1+ndig-decpt;
|
||||
/*if(prec) {
|
||||
PRINT(t, ".", 1);
|
||||
PAD(t, -decpt, '0');
|
||||
PRINT(t, retBuf, ndig);
|
||||
} */
|
||||
/* if (doAlt) { t += prec-ndig; } // superflous if flags|=ALT properly set on bsd_cvt call */
|
||||
} else
|
||||
if(decpt >= ndig) {
|
||||
PRINT(t, retBuf, ndig);
|
||||
PAD(t, decpt - ndig, '0');
|
||||
if(doAlt) {
|
||||
PRINT(t, ".", 1);
|
||||
/* PAD(t, prec - decpt, '0'); // superflous if flags|=ALT properly set on bsd_cvt call */
|
||||
}
|
||||
}
|
||||
else {
|
||||
char *s = retBuf;
|
||||
PRINT(t, s, decpt);
|
||||
s += decpt;
|
||||
PRINT(t, ".", 1);
|
||||
PRINT(t, s, ndig - decpt);
|
||||
/* if(doAlt) { PAD(t, prec - ndig, '0'); } // superflous if flags|=ALT properly set on bsd_cvt call */
|
||||
}
|
||||
}
|
||||
else {
|
||||
if(fabs(number) == 0.0) { /* handle special independent of if (doAlt) 0.[0..] */
|
||||
PRINT(t, "0", 1);
|
||||
PRINT(t, ".", 1);
|
||||
PAD(t, prec, '0'); /* Always pad, even without doAlt */
|
||||
} else
|
||||
if(ndig > 1 || doAlt) {
|
||||
char *s = retBuf;
|
||||
PRINT(t, s, 1);
|
||||
s++;
|
||||
PRINT(t, ".", 1);
|
||||
PRINT(t, s, ndig - 1);
|
||||
/* ex. 3.10000E+06 : prec=6, ndig=2*/
|
||||
/* if(doAlt) { PAD(t, prec - ndig, '0'); } // superflous if flags|=ALT properly set on bsd_cvt call */
|
||||
/* strip off the trailing visible 0, if not padding to full precision,
|
||||
* (despite possible not 0 digits beyond the visible precision)
|
||||
* if then the last char is a '.', strip it off also.
|
||||
* if not added needlessly no need to strip
|
||||
* if (!doAlt) { while (*(t-1) == '0') { --t; } if (*(t-1) == '.') { --t; } }
|
||||
*/
|
||||
}
|
||||
else {
|
||||
PRINT(t, retBuf, 1);
|
||||
}
|
||||
PRINT(t, expstr, nexp);
|
||||
}
|
||||
|
||||
/***************************************************************************
|
||||
* End of import of David M. Gay.
|
||||
***************************************************************************/
|
||||
|
||||
/*free(retBuf); // no longer allocating, see above */
|
||||
return t - startp;
|
||||
}
|
||||
|
||||
/******************************************************************************
|
||||
*
|
||||
* exponentCvt - helper function for floatFormat
|
||||
*
|
||||
* RETURNS:
|
||||
*/
|
||||
|
||||
LOCAL char *exponentCvt
|
||||
(
|
||||
FAST char * p,
|
||||
FAST int exp,
|
||||
int fmtch
|
||||
)
|
||||
{
|
||||
FAST char *t;
|
||||
char expbuf[MAXEXP];
|
||||
|
||||
*p++ = fmtch;
|
||||
if (exp < 0)
|
||||
{
|
||||
exp = -exp;
|
||||
*p++ = '-';
|
||||
}
|
||||
else
|
||||
*p++ = '+';
|
||||
|
||||
t = expbuf + MAXEXP;
|
||||
|
||||
if (exp > 9)
|
||||
{
|
||||
do
|
||||
{
|
||||
*--t = to_char(exp % 10);
|
||||
} while ((exp /= 10) > 9);
|
||||
|
||||
*--t = to_char(exp);
|
||||
|
||||
for (; t < expbuf + MAXEXP; *p++ = *t++)
|
||||
;
|
||||
}
|
||||
else
|
||||
{
|
||||
*p++ = '0';
|
||||
*p++ = to_char(exp);
|
||||
}
|
||||
|
||||
return (p);
|
||||
}
|
||||
|
||||
/*******************************************************************************
|
||||
*
|
||||
* floatScan - scan ASCII input to a floating-point number
|
||||
*
|
||||
* This routine scans the ASCII input into a floating-point number.
|
||||
*
|
||||
* [+-]<digs>.<digs>e[+-]<digs>
|
||||
*
|
||||
* RETURNS: TRUE if successful, FALSE if unsuccessful
|
||||
*/
|
||||
LOCAL BOOL floatScan
|
||||
(
|
||||
FAST int * pReturn,
|
||||
int returnSpec, /* 0, 'l', 'L' */
|
||||
int fieldwidth,
|
||||
FAST GET_NEXT_CHAR_FUNCPTR getRtn,
|
||||
FAST void *getArg,
|
||||
int * pCh,
|
||||
int * pNchars
|
||||
)
|
||||
{
|
||||
|
||||
#define GET_CHAR(ch, ix) ((ix)++, (ch) = getRtn (getArg))
|
||||
|
||||
static double power10 [] =
|
||||
{1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8};
|
||||
|
||||
FAST int ch = *pCh; /* current char */
|
||||
FAST int ix = 0; /* number of chars consumed */
|
||||
FAST int long1 = 0; /* 1st part of integer representation */
|
||||
FAST int long2 = 0; /* 2nd part of integer representation */
|
||||
int ndigs1 = 0; /* num of digits in long1 */
|
||||
int ndigs2 = 0; /* num of digits in long2 */
|
||||
int ndigsPrec = 0; /* num of digits of precision total */
|
||||
BOOL fracPart = FALSE; /* TRUE = doing digits after '.' */
|
||||
int ndigsFrac = 0; /* number of digits after '.' */
|
||||
int exp = 0; /* exponent */
|
||||
double num; /* double representation of number */
|
||||
BOOL negnum = FALSE; /* number is negative */
|
||||
BOOL negexp = FALSE; /* exponent is negative */
|
||||
BOOL digitsSeen = FALSE; /* set if any digits are encountered */
|
||||
double adj_exp;
|
||||
const double *fp;
|
||||
int exp1;
|
||||
|
||||
/* check for sign */
|
||||
|
||||
if (ix < fieldwidth)
|
||||
{
|
||||
if ((negnum = ((char)ch == '-')) || ((char)ch == '+'))
|
||||
GET_CHAR (ch, ix);
|
||||
}
|
||||
|
||||
/* scan integer and fraction parts */
|
||||
|
||||
for (; (ch != EOF) && (ix < fieldwidth); GET_CHAR (ch, ix))
|
||||
{
|
||||
if (!fracPart && ch == '.')
|
||||
{
|
||||
fracPart = TRUE;
|
||||
continue;
|
||||
}
|
||||
|
||||
if ((ch < '0') || (ch > '9'))
|
||||
{
|
||||
switch (ch)
|
||||
{
|
||||
case 'i':
|
||||
case 'I':
|
||||
GET_CHAR (ch, ix);
|
||||
if (ch != EOF && ix < fieldwidth && (ch == 'n' || ch == 'N'))
|
||||
{
|
||||
GET_CHAR (ch, ix);
|
||||
if (ch != EOF && ix < fieldwidth && (ch == 'f' || ch == 'F'))
|
||||
{
|
||||
/* Infinity */
|
||||
((DOUBLE *) &num)->ldat.l1 = 0x7ff00000;
|
||||
((DOUBLE *) &num)->ldat.l2 = 0;
|
||||
goto returnToCaller;
|
||||
}
|
||||
}
|
||||
break;
|
||||
case 'n':
|
||||
case 'N':
|
||||
GET_CHAR (ch, ix);
|
||||
if (ch != EOF && ix < fieldwidth && (ch == 'a' || ch == 'A'))
|
||||
{
|
||||
GET_CHAR (ch, ix);
|
||||
if (ch != EOF && ix < fieldwidth && (ch == 'n' || ch == 'N'))
|
||||
{
|
||||
/* Not a Number */
|
||||
((DOUBLE *) &num)->ldat.l1 = 0x7ff80000;
|
||||
((DOUBLE *) &num)->ldat.l2 = 1;
|
||||
goto returnToCaller;
|
||||
}
|
||||
}
|
||||
break;
|
||||
default:
|
||||
break;
|
||||
}
|
||||
break;
|
||||
}
|
||||
|
||||
digitsSeen = TRUE;
|
||||
|
||||
ch -= '0'; /* turn ch into digit */
|
||||
|
||||
if (fracPart) /* count digits after '.' */
|
||||
ndigsFrac++;
|
||||
|
||||
if ((ndigsPrec != 0) || (ch != 0)) /* skip leading 0s */
|
||||
{
|
||||
ndigsPrec++;
|
||||
|
||||
/* check if another digit will fit in long1 */
|
||||
|
||||
if (long1 < ((INT_MAX - 9) / 10))
|
||||
{
|
||||
long1 = 10 * long1 + ch;
|
||||
ndigs1++;
|
||||
}
|
||||
else if (ndigsPrec <= DBL_DIG)
|
||||
{
|
||||
long2 = 10 * long2 + ch;
|
||||
ndigs2++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* scan for explicit exponent */
|
||||
|
||||
if ((ix < fieldwidth) && (ch == 'e' || ch == 'E'))
|
||||
{
|
||||
GET_CHAR (ch, ix);
|
||||
|
||||
/* skip over sign */
|
||||
|
||||
if (ix < fieldwidth)
|
||||
{
|
||||
if ((negexp = (ch == '-')) || (ch == '+'))
|
||||
GET_CHAR (ch, ix);
|
||||
}
|
||||
|
||||
for (; (ch != EOF) && (ix < fieldwidth); GET_CHAR (ch, ix))
|
||||
{
|
||||
if ((ch < '0') || (ch > '9'))
|
||||
break;
|
||||
|
||||
exp = 10 * exp + (ch - '0');
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* check that we scanned at least one character */
|
||||
|
||||
if (!digitsSeen)
|
||||
return (FALSE);
|
||||
|
||||
/* put significant digits into double representation */
|
||||
|
||||
num = long1;
|
||||
if (ndigs2 != 0)
|
||||
num = num * power10 [ndigs2] + long2;
|
||||
|
||||
|
||||
/* apply exponent to number; total effective exponent =
|
||||
* + number of digits of precision that were scanned but not
|
||||
* represented in the number (ndigsPrec - ndigs1 - ndigs2)
|
||||
* - number of digits scanned after the '.' (ndigsFrac)
|
||||
* + the scanned explicit exponent
|
||||
*/
|
||||
|
||||
exp = (ndigsPrec - ndigs1 - ndigs2) - ndigsFrac + (negexp ? -exp : exp);
|
||||
|
||||
/* We now have result = num * 10**exp. Change that to
|
||||
num1 * 2**exp1 => num * 1.25**exp * 2**(exp*3) */
|
||||
|
||||
adj_exp = 1.;
|
||||
exp1 = exp;
|
||||
if (exp1 < 0) exp1 = -exp1;
|
||||
fp = f125;
|
||||
while (exp1)
|
||||
{ /* compensate mantissa */
|
||||
if (exp1 & 1) adj_exp *= *fp;
|
||||
exp1 >>= 1;
|
||||
if (fp < F125MAX) fp++;
|
||||
}
|
||||
if (exp >= 0) num *= adj_exp;
|
||||
else num /= adj_exp;
|
||||
exp *= 3;
|
||||
num = ldexp(num, exp);
|
||||
|
||||
/* return value to caller */
|
||||
|
||||
returnToCaller:
|
||||
if (negnum)
|
||||
num = -num;
|
||||
|
||||
if (pReturn != NULL)
|
||||
{
|
||||
switch (returnSpec)
|
||||
{
|
||||
case 'l': * (double *) pReturn = num; break;
|
||||
default: * (float *) pReturn = num; break;
|
||||
}
|
||||
}
|
||||
|
||||
*pCh = ch;
|
||||
*pNchars += ix;
|
||||
|
||||
return (ix != 0);
|
||||
}
|
||||
@@ -0,0 +1,143 @@
|
||||
/* fpType.c - floating point type functions */
|
||||
|
||||
/* Copyright 1992-2001 Wind River Systems, Inc. */
|
||||
|
||||
/*
|
||||
modification history
|
||||
--------------------
|
||||
01i,28jul09,fan remove warning for LP64 data model
|
||||
01h,03oct01,to use IEEE little endian for ARM
|
||||
01g,16may00,pai Implemented fixes for incorrect mask usage in fpTypeGet
|
||||
(SPR #7515), (SPR #8515), (SPR #20508).
|
||||
01f,23jan97,cdp reverse words in double for ARM.
|
||||
01e,19mar95,dvs removed tron reference - no longer supported.
|
||||
01d,01sep93,hdn added byte order condition for I80X86.
|
||||
01c,17aug92,kdl added dummy pointer declaration to work around tron gcc bug.
|
||||
01b,30jul92,kdl removed unused static variables.
|
||||
01a,30jul92,kdl created.
|
||||
*/
|
||||
|
||||
/*
|
||||
DESCRIPTION
|
||||
|
||||
This file contains a routine, fpTypeGet(), which will return the
|
||||
type of a floating point number, where the type is one of NULL,
|
||||
NAN, INTEGER, INF, or REAL.
|
||||
|
||||
NOMANUAL
|
||||
*/
|
||||
|
||||
#include "vxWorks.h"
|
||||
#include "math.h"
|
||||
#include "private/mathP.h"
|
||||
|
||||
#if (_BYTE_ORDER == _BIG_ENDIAN)
|
||||
|
||||
#define LOW 0 /* array indices for mantissa */
|
||||
#define HIGH 1
|
||||
|
||||
#else
|
||||
|
||||
#define LOW 1 /* array indices for mantissa */
|
||||
#define HIGH 0
|
||||
|
||||
#endif /* _BYTE_ORDER == _BIG_ENDIAN */
|
||||
|
||||
|
||||
static unsigned map[32] =
|
||||
{
|
||||
0xffffffff, 0xfffffffe, 0xfffffffc, 0xfffffff8,
|
||||
0xfffffff0, 0xffffffe0, 0xffffffc0, 0xffffff80,
|
||||
0xffffff00, 0xfffffe00, 0xfffffc00, 0xfffff800,
|
||||
0xfffff000, 0xffffe000, 0xffffc000, 0xffff8000,
|
||||
0xffff0000, 0xfffe0000, 0xfffc0000, 0xfff80000,
|
||||
0xfff00000, 0xffe00000, 0xffc00000, 0xff800000,
|
||||
0xff000000, 0xfe000000, 0xfc000000, 0xf8000000,
|
||||
0xf0000000, 0xe0000000, 0xc0000000, 0x80000000,
|
||||
};
|
||||
|
||||
/*******************************************************************************
|
||||
*
|
||||
* fpTypeGet - determine type of fp number
|
||||
*
|
||||
* This routine examines a double floating point value and returns
|
||||
* a code indicating the type (i.e. ZERO, INF, NAN, REAL, INT).
|
||||
*
|
||||
* RETURNS: encoded type (value).
|
||||
*
|
||||
* NOMANUAL
|
||||
*/
|
||||
|
||||
int fpTypeGet
|
||||
(
|
||||
double v, /* value */
|
||||
double * r /* pointer to integral */
|
||||
)
|
||||
{
|
||||
__cv_type cv;
|
||||
int exp;
|
||||
BOOL any;
|
||||
|
||||
|
||||
cv.p_double = v;
|
||||
|
||||
if (((cv.p_mant[LOW] & 0x7fffffff) == 0) && (cv.p_mant[HIGH] == 0))
|
||||
{
|
||||
if (r != NULL)
|
||||
*r = v;
|
||||
return (ZERO);
|
||||
|
||||
}
|
||||
|
||||
exp = (int)((cv.p_mant[LOW] & 0x7ff00000) >> 20);
|
||||
|
||||
if (exp == 2047)
|
||||
{
|
||||
if (r != NULL)
|
||||
*r = v;
|
||||
return (((cv.p_mant[LOW] & 0x000fffff) == 0) &&
|
||||
(cv.p_mant[HIGH] == 0)) ? INF : NAN;
|
||||
}
|
||||
|
||||
if (r == NULL)
|
||||
return (REAL);
|
||||
|
||||
exp -= 1023;
|
||||
any = FALSE;
|
||||
|
||||
if (exp < 0)
|
||||
{
|
||||
any = TRUE;
|
||||
cv.p_mant[LOW] &= 0x80000000;
|
||||
cv.p_mant[HIGH] = 0;
|
||||
}
|
||||
else if (exp < 20)
|
||||
{
|
||||
if ((cv.p_mant[HIGH] != 0) || ((cv.p_mant[LOW] & ~(map[20 - exp])) != 0))
|
||||
{
|
||||
any = TRUE;
|
||||
cv.p_mant[HIGH] = 0;
|
||||
cv.p_mant[LOW] &= map[20 - exp];
|
||||
}
|
||||
}
|
||||
else if (exp == 20)
|
||||
{
|
||||
if ((cv.p_mant[HIGH] & ~(0x00000000)) != 0)
|
||||
{
|
||||
any = TRUE;
|
||||
cv.p_mant[HIGH] &= 0x00000000;
|
||||
}
|
||||
}
|
||||
else if (exp < 52)
|
||||
{
|
||||
if ((cv.p_mant[HIGH] & ~(map [52 - exp])) != 0)
|
||||
{
|
||||
any = TRUE;
|
||||
cv.p_mant[HIGH] &= map [52 - exp];
|
||||
}
|
||||
}
|
||||
|
||||
*r = cv.p_double;
|
||||
|
||||
return (any ? REAL : INTEGER);
|
||||
}
|
||||
@@ -0,0 +1,44 @@
|
||||
/* ieeefp.h */
|
||||
|
||||
/* Copyright 1992-2001 Wind River Systems, Inc. */
|
||||
|
||||
/*
|
||||
modification history
|
||||
--------------------
|
||||
01c,08may03,pes PAL conditional compilation cleanup. Phase 2.
|
||||
01b,14nov01,to use IEEE little endian for ARM
|
||||
01a,31oct01,sn better deduction of endianness
|
||||
added this mod history
|
||||
*/
|
||||
|
||||
#ifndef __INCieeefph
|
||||
#define __INCieeefph
|
||||
|
||||
#include "vxWorks.h"
|
||||
|
||||
#if defined (__H8300__) || defined (__H8300H__)
|
||||
#define __SMALL_BITFIELDS
|
||||
#define _DOUBLE_IS_32BITS
|
||||
#endif
|
||||
|
||||
#ifdef __H8500__
|
||||
#define __SMALL_BITFIELDS
|
||||
#define _DOUBLE_IS_32BITS
|
||||
#endif
|
||||
|
||||
#ifdef __W65__
|
||||
#define __SMALL_BITFIELDS
|
||||
#define _DOUBLE_IS_32BITS
|
||||
#endif
|
||||
|
||||
#ifndef __IEEE_BIG_ENDIAN
|
||||
#ifndef __IEEE_LITTLE_ENDIAN
|
||||
#if _BYTE_ORDER == _LITTLE_ENDIAN
|
||||
#define __IEEE_LITTLE_ENDIAN
|
||||
#else
|
||||
#define __IEEE_BIG_ENDIAN
|
||||
#endif /* _BYTE_ORDER == _LITTLE_ENDIAN */
|
||||
#endif /* not __IEEE_LITTLE_ENDIAN */
|
||||
#endif /* not __IEEE_BIG_ENDIAN */
|
||||
|
||||
#endif /* __INCieeefph */
|
||||
File diff suppressed because it is too large
Load Diff
@@ -0,0 +1,391 @@
|
||||
/****************************************************************
|
||||
*
|
||||
* The author of this software is David M. Gay.
|
||||
*
|
||||
* Copyright (c) 1991, 2000 by AT&T.
|
||||
*
|
||||
* Permission to use, copy, modify, and distribute this software for any
|
||||
* purpose without fee is hereby granted, provided that this entire notice
|
||||
* is included in all copies of any software which is or includes a copy
|
||||
* or modification of this software and in all copies of the supporting
|
||||
* documentation for such software.
|
||||
*
|
||||
* THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
|
||||
* WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY
|
||||
* REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
|
||||
* OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
|
||||
*
|
||||
***************************************************************/
|
||||
|
||||
/* Please send bug reports to
|
||||
David M. Gay
|
||||
AT&T Bell Laboratories, Room 2C-463
|
||||
600 Mountain Avenue
|
||||
Murray Hill, NJ 07974-2070
|
||||
U.S.A.
|
||||
dmg@research.att.com or research!dmg
|
||||
*/
|
||||
|
||||
#ifndef __CLASSPATH_MPREC_H__
|
||||
#define __CLASSPATH_MPREC_H__
|
||||
|
||||
#if 0
|
||||
#include <config.h>
|
||||
#include "config-int.h"
|
||||
#endif
|
||||
#include "ieeefp.h"
|
||||
/* CLASSPATH LOCAL */
|
||||
/* #include "namespace.h" */
|
||||
|
||||
#if defined HAVE_SYS_TYPES_H
|
||||
#include <sys/types.h>
|
||||
#endif
|
||||
|
||||
#if defined HAVE_SYS_CONFIG_H
|
||||
#include <sys/config.h>
|
||||
#endif
|
||||
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
|
||||
/* These typedefs are true for the targets running Java. */
|
||||
|
||||
#ifdef __IEEE_LITTLE_ENDIAN
|
||||
#define IEEE_8087
|
||||
#endif
|
||||
|
||||
#ifdef __IEEE_BIG_ENDIAN
|
||||
#define IEEE_MC68k
|
||||
#endif
|
||||
|
||||
#ifdef __Z8000__
|
||||
#define Just_16
|
||||
#endif
|
||||
|
||||
#ifdef DEBUG
|
||||
#include "stdio.h"
|
||||
#include <stdlib.h>
|
||||
#define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
|
||||
#endif
|
||||
|
||||
|
||||
#ifdef Unsigned_Shifts
|
||||
#define Sign_Extend(a,b) if (b < 0) a |= (uint32_t)0xffff0000;
|
||||
#else
|
||||
#define Sign_Extend(a,b) /*no-op*/
|
||||
#endif
|
||||
|
||||
#if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
|
||||
Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
|
||||
#endif
|
||||
|
||||
#define __Long int32_t
|
||||
#define __ULong uint32_t
|
||||
|
||||
/* If we are going to examine or modify specific bits in a double using
|
||||
the word0 and/or word1 macros, then we must wrap the double inside
|
||||
a union. This is necessary to avoid undefined behavior according to
|
||||
the ANSI C spec. */
|
||||
union double_union
|
||||
{
|
||||
double d;
|
||||
uint32_t i[2];
|
||||
};
|
||||
|
||||
#ifdef IEEE_8087
|
||||
#define word0(x) (x.i[1])
|
||||
#define word1(x) (x.i[0])
|
||||
#else
|
||||
#define word0(x) (x.i[0])
|
||||
#define word1(x) (x.i[1])
|
||||
#endif
|
||||
|
||||
/* The following definition of Storeinc is appropriate for MIPS processors.
|
||||
* An alternative that might be better on some machines is
|
||||
* #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
|
||||
*/
|
||||
#if defined(IEEE_8087) + defined(VAX)
|
||||
#define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
|
||||
((unsigned short *)a)[0] = (unsigned short)c, a++)
|
||||
#else
|
||||
#define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
|
||||
((unsigned short *)a)[1] = (unsigned short)c, a++)
|
||||
#endif
|
||||
|
||||
/* #define P DBL_MANT_DIG */
|
||||
/* Ten_pmax = floor(P*log(2)/log(5)) */
|
||||
/* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
|
||||
/* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
|
||||
/* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
|
||||
|
||||
#if defined(IEEE_8087) + defined(IEEE_MC68k)
|
||||
#if defined (_DOUBLE_IS_32BITS)
|
||||
#define Exp_shift 23
|
||||
#define Exp_shift1 23
|
||||
#define Exp_msk1 ((uint32_t)0x00800000L)
|
||||
#define Exp_msk11 ((uint32_t)0x00800000L)
|
||||
#define Exp_mask ((uint32_t)0x7f800000L)
|
||||
#define P 24
|
||||
#define Bias 127
|
||||
#if 0
|
||||
#define IEEE_Arith /* it is, but the code doesn't handle IEEE singles yet */
|
||||
#endif
|
||||
#define Emin (-126)
|
||||
#define Exp_1 ((uint32_t)0x3f800000L)
|
||||
#define Exp_11 ((uint32_t)0x3f800000L)
|
||||
#define Ebits 8
|
||||
#define Frac_mask ((uint32_t)0x007fffffL)
|
||||
#define Frac_mask1 ((uint32_t)0x007fffffL)
|
||||
#define Ten_pmax 10
|
||||
#define Sign_bit ((uint32_t)0x80000000L)
|
||||
#define Ten_pmax 10
|
||||
#define Bletch 2
|
||||
#define Bndry_mask ((uint32_t)0x007fffffL)
|
||||
#define Bndry_mask1 ((uint32_t)0x007fffffL)
|
||||
#define _LSB 1
|
||||
#define Sign_bit ((uint32_t)0x80000000L)
|
||||
#define Log2P 1
|
||||
#define Tiny0 0
|
||||
#define Tiny1 1
|
||||
#define Quick_max 5
|
||||
#define Int_max 6
|
||||
#define Infinite(x) (word0(x) == ((uint32_t)0x7f800000L))
|
||||
#undef word0
|
||||
#undef word1
|
||||
|
||||
#define word0(x) (x.i[0])
|
||||
#define word1(x) 0
|
||||
#else
|
||||
|
||||
#define Exp_shift 20
|
||||
#define Exp_shift1 20
|
||||
#define Exp_msk1 ((uint32_t)0x100000L)
|
||||
#define Exp_msk11 ((uint32_t)0x100000L)
|
||||
#define Exp_mask ((uint32_t)0x7ff00000L)
|
||||
#define P 53
|
||||
#define Bias 1023
|
||||
#define IEEE_Arith
|
||||
#define Emin (-1022)
|
||||
#define Exp_1 ((uint32_t)0x3ff00000L)
|
||||
#define Exp_11 ((uint32_t)0x3ff00000L)
|
||||
#define Ebits 11
|
||||
#define Frac_mask ((uint32_t)0xfffffL)
|
||||
#define Frac_mask1 ((uint32_t)0xfffffL)
|
||||
#define Ten_pmax 22
|
||||
#define Bletch 0x10
|
||||
#define Bndry_mask ((uint32_t)0xfffffL)
|
||||
#define Bndry_mask1 ((uint32_t)0xfffffL)
|
||||
#define _LSB 1
|
||||
#define Sign_bit ((uint32_t)0x80000000L)
|
||||
#define Log2P 1
|
||||
#define Tiny0 0
|
||||
#define Tiny1 1
|
||||
#define Quick_max 14
|
||||
#define Int_max 14
|
||||
#define Infinite(x) (word0(x) == ((uint32_t)0x7ff00000L)) /* sufficient test for here */
|
||||
#endif
|
||||
|
||||
#else
|
||||
#undef Sudden_Underflow
|
||||
#define Sudden_Underflow
|
||||
#ifdef IBM
|
||||
#define Exp_shift 24
|
||||
#define Exp_shift1 24
|
||||
#define Exp_msk1 ((uint32_t)0x1000000L)
|
||||
#define Exp_msk11 ((uint32_t)0x1000000L)
|
||||
#define Exp_mask ((uint32_t)0x7f000000L)
|
||||
#define P 14
|
||||
#define Bias 65
|
||||
#define Exp_1 ((uint32_t)0x41000000L)
|
||||
#define Exp_11 ((uint32_t)0x41000000L)
|
||||
#define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
|
||||
#define Frac_mask ((uint32_t)0xffffffL)
|
||||
#define Frac_mask1 ((uint32_t)0xffffffL)
|
||||
#define Bletch 4
|
||||
#define Ten_pmax 22
|
||||
#define Bndry_mask ((uint32_t)0xefffffL)
|
||||
#define Bndry_mask1 ((uint32_t)0xffffffL)
|
||||
#define _LSB 1
|
||||
#define Sign_bit ((uint32_t)0x80000000L)
|
||||
#define Log2P 4
|
||||
#define Tiny0 ((uint32_t)0x100000L)
|
||||
#define Tiny1 0
|
||||
#define Quick_max 14
|
||||
#define Int_max 15
|
||||
#else /* VAX */
|
||||
#define Exp_shift 23
|
||||
#define Exp_shift1 7
|
||||
#define Exp_msk1 0x80
|
||||
#define Exp_msk11 ((uint32_t)0x800000L)
|
||||
#define Exp_mask ((uint32_t)0x7f80L)
|
||||
#define P 56
|
||||
#define Bias 129
|
||||
#define Exp_1 ((uint32_t)0x40800000L)
|
||||
#define Exp_11 ((uint32_t)0x4080L)
|
||||
#define Ebits 8
|
||||
#define Frac_mask ((uint32_t)0x7fffffL)
|
||||
#define Frac_mask1 ((uint32_t)0xffff007fL)
|
||||
#define Ten_pmax 24
|
||||
#define Bletch 2
|
||||
#define Bndry_mask ((uint32_t)0xffff007fL)
|
||||
#define Bndry_mask1 ((uint32_t)0xffff007fL)
|
||||
#define _LSB ((uint32_t)0x10000L)
|
||||
#define Sign_bit ((uint32_t)0x8000L)
|
||||
#define Log2P 1
|
||||
#define Tiny0 0x80
|
||||
#define Tiny1 0
|
||||
#define Quick_max 15
|
||||
#define Int_max 15
|
||||
#endif
|
||||
#endif
|
||||
|
||||
#ifndef IEEE_Arith
|
||||
#define ROUND_BIASED
|
||||
#endif
|
||||
|
||||
#ifdef RND_PRODQUOT
|
||||
#define rounded_product(a,b) a = rnd_prod(a, b)
|
||||
#define rounded_quotient(a,b) a = rnd_quot(a, b)
|
||||
#ifdef KR_headers
|
||||
extern double rnd_prod(), rnd_quot();
|
||||
#else
|
||||
extern double rnd_prod(double, double), rnd_quot(double, double);
|
||||
#endif
|
||||
#else
|
||||
#define rounded_product(a,b) a *= b
|
||||
#define rounded_quotient(a,b) a /= b
|
||||
#endif
|
||||
|
||||
#define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
|
||||
#define Big1 ((uint32_t)0xffffffffL)
|
||||
|
||||
#ifndef Just_16
|
||||
/* When Pack_32 is not defined, we store 16 bits per 32-bit long.
|
||||
* This makes some inner loops simpler and sometimes saves work
|
||||
* during multiplications, but it often seems to make things slightly
|
||||
* slower. Hence the default is now to store 32 bits per long.
|
||||
*/
|
||||
|
||||
#ifndef Pack_32
|
||||
#if SIZEOF_VOID_P != 8
|
||||
#define Pack_32
|
||||
#endif
|
||||
#endif
|
||||
#endif
|
||||
|
||||
|
||||
#define MAX_BIGNUMS 16
|
||||
#ifdef Pack_32
|
||||
#define MAX_BIGNUM_WDS 32
|
||||
#else
|
||||
/* Note that this is a workaround for */
|
||||
#define MAX_BIGNUM_WDS 128
|
||||
#endif
|
||||
|
||||
struct _Jv_Bigint
|
||||
{
|
||||
struct _Jv_Bigint *_next;
|
||||
int _k, _maxwds, _sign, _wds;
|
||||
__ULong _x[1];
|
||||
};
|
||||
|
||||
|
||||
#define _PTR void *
|
||||
#define _AND ,
|
||||
#define _NOARGS void
|
||||
#define _CONST const
|
||||
#define _VOLATILE volatile
|
||||
#define _SIGNED signed
|
||||
#define _DOTS , ...
|
||||
#define _VOID void
|
||||
#define _EXFUN(name, proto) name proto
|
||||
#define _DEFUN(name, arglist, args) name(args)
|
||||
#define _DEFUN_VOID(name) name(_NOARGS)
|
||||
#define _CAST_VOID (void)
|
||||
|
||||
|
||||
struct _Jv_reent
|
||||
{
|
||||
/* local copy of errno */
|
||||
int _errno;
|
||||
|
||||
/* used by mprec routines */
|
||||
struct _Jv_Bigint *_result;
|
||||
int _result_k;
|
||||
struct _Jv_Bigint *_p5s;
|
||||
|
||||
struct _Jv_Bigint **_freelist;
|
||||
int _max_k;
|
||||
};
|
||||
|
||||
|
||||
typedef struct _Jv_Bigint _Jv_Bigint;
|
||||
|
||||
#define Balloc _Jv_Balloc
|
||||
#define Bfree _Jv_Bfree
|
||||
#define multadd _Jv_multadd
|
||||
#define s2b _Jv_s2b
|
||||
#define lo0bits _Jv_lo0bits
|
||||
#define hi0bits _Jv_hi0bits
|
||||
#define i2b _Jv_i2b
|
||||
#define mult _Jv_mult
|
||||
#define pow5mult _Jv_pow5mult
|
||||
#define lshift _Jv_lshift
|
||||
#define cmp _Jv__mcmp
|
||||
#define diff _Jv__mdiff
|
||||
#define ulp _Jv_ulp
|
||||
#define b2d _Jv_b2d
|
||||
#define d2b _Jv_d2b
|
||||
#define ratio _Jv_ratio
|
||||
|
||||
#define tens _Jv__mprec_tens
|
||||
#define bigtens _Jv__mprec_bigtens
|
||||
#define tinytens _Jv__mprec_tinytens
|
||||
|
||||
/* #define _dtoa _Jv_dtoa */
|
||||
#define _dtoa __dtoa
|
||||
#define _dtoa_r _Jv_dtoa_r
|
||||
#define _strtod_r _Jv_strtod_r
|
||||
|
||||
extern double _EXFUN(_strtod_r, (struct _Jv_reent *ptr, const char *s00, char **se));
|
||||
extern char* _EXFUN(_dtoa_r, (struct _Jv_reent *ptr, double d,
|
||||
int mode, int ndigits, int *decpt, int *sign,
|
||||
char **rve, int float_type));
|
||||
/* void _EXFUN(_dtoa, (double d, int mode, int ndigits, int *decpt, int *sign,
|
||||
char **rve, char *buf, int float_type)); */
|
||||
|
||||
double _EXFUN(ulp,(double x));
|
||||
double _EXFUN(b2d,(_Jv_Bigint *a , int *e));
|
||||
_Jv_Bigint * _EXFUN(Balloc,(struct _Jv_reent *p, int k));
|
||||
void _EXFUN(Bfree,(struct _Jv_reent *p, _Jv_Bigint *v));
|
||||
_Jv_Bigint * _EXFUN(multadd,(struct _Jv_reent *p, _Jv_Bigint *, int, int));
|
||||
_Jv_Bigint * _EXFUN(s2b,(struct _Jv_reent *, const char*, int, int, __ULong));
|
||||
_Jv_Bigint * _EXFUN(i2b,(struct _Jv_reent *,int));
|
||||
_Jv_Bigint * _EXFUN(mult, (struct _Jv_reent *, _Jv_Bigint *, _Jv_Bigint *));
|
||||
_Jv_Bigint * _EXFUN(pow5mult, (struct _Jv_reent *, _Jv_Bigint *, int k));
|
||||
int _EXFUN(hi0bits,(__ULong));
|
||||
int _EXFUN(lo0bits,(__ULong*));
|
||||
_Jv_Bigint * _EXFUN(d2b,(struct _Jv_reent *p, double d, int *e, int *bits));
|
||||
_Jv_Bigint * _EXFUN(lshift,(struct _Jv_reent *p, _Jv_Bigint *b, int k));
|
||||
_Jv_Bigint * _EXFUN(diff,(struct _Jv_reent *p, _Jv_Bigint *a, _Jv_Bigint *b));
|
||||
int _EXFUN(cmp,(_Jv_Bigint *a, _Jv_Bigint *b));
|
||||
|
||||
double _EXFUN(ratio,(_Jv_Bigint *a, _Jv_Bigint *b));
|
||||
#define Bcopy(x,y) memcpy((char *)&x->_sign, (char *)&y->_sign, y->_wds*sizeof(__Long) + 2*sizeof(int))
|
||||
|
||||
#if defined(_DOUBLE_IS_32BITS) && defined(__v800)
|
||||
#define n_bigtens 2
|
||||
#else
|
||||
#define n_bigtens 5
|
||||
#endif
|
||||
|
||||
extern _CONST double tinytens[];
|
||||
extern _CONST double bigtens[];
|
||||
extern _CONST double tens[];
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
||||
#endif /* __CLASSPATH_MPREC_H__ */
|
||||
@@ -0,0 +1,74 @@
|
||||
/* rint.c - round-to-integer routine */
|
||||
|
||||
/* Copyright 1992 Wind River Systems, Inc. */
|
||||
|
||||
/*
|
||||
modification history
|
||||
--------------------
|
||||
01a,30jul92,smb documentation.
|
||||
*/
|
||||
|
||||
/*
|
||||
DESCRIPTION
|
||||
* algorithm for rint(x) in pseudo-pascal form ...
|
||||
*
|
||||
* real rint(x): real x;
|
||||
* ... delivers integer nearest x in direction of prevailing rounding
|
||||
* ... mode
|
||||
* const L = (last consecutive integer)/2
|
||||
* = 2**55; for VAX D
|
||||
* = 2**52; for IEEE 754 Double
|
||||
* real s,t;
|
||||
* begin
|
||||
* if x != x then return x; ... NaN
|
||||
* if |x| >= L then return x; ... already an integer
|
||||
* s := copysign(L,x);
|
||||
* t := x + s; ... = (x+s) rounded to integer
|
||||
* return t - s
|
||||
* end;
|
||||
*
|
||||
* Note: Inexact will be signaled if x is not an integer, as is
|
||||
* customary for IEEE 754. No other signal can be emitted.
|
||||
|
||||
SEE ALSO: American National Standard X3.159-1989
|
||||
*/
|
||||
|
||||
#include "vxWorks.h"
|
||||
#include "math.h"
|
||||
#include "stddef.h"
|
||||
#include "errno.h"
|
||||
#include "private/mathP.h"
|
||||
|
||||
#if defined(vax)||defined(tahoe)
|
||||
#ifdef vax
|
||||
#define _0x(A,B) 0x/**/A/**/B
|
||||
#else /* vax */
|
||||
#define _0x(A,B) 0x/**/B/**/A
|
||||
#endif /* vax */
|
||||
static long Lx[] = {_0x(0000,5c00),_0x(0000,0000)}; /* 2**55 */
|
||||
#define L *(double *) Lx
|
||||
#else /* defined(vax)||defined(tahoe) */
|
||||
static double L = 4503599627370496.0E0; /* 2**52 */
|
||||
#endif /* defined(vax)||defined(tahoe) */
|
||||
|
||||
/********************************************************************************
|
||||
* rint - needed for sparc math support
|
||||
*
|
||||
* RETURNS:
|
||||
*/
|
||||
|
||||
double
|
||||
rint(x)
|
||||
double x;
|
||||
{
|
||||
double s,t,one = 1.0,copysign();
|
||||
#if !defined(vax)&&!defined(tahoe)
|
||||
if (x != x) /* NaN */
|
||||
return (x);
|
||||
#endif /* !defined(vax)&&!defined(tahoe) */
|
||||
if (copysign(x,one) >= L) /* already an integer */
|
||||
return (x);
|
||||
s = copysign(L,x);
|
||||
t = x + s; /* x+s rounded to integer */
|
||||
return (t - s);
|
||||
}
|
||||
@@ -0,0 +1,86 @@
|
||||
|
||||
/* @(#)s_rint.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.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
/*
|
||||
* rint(x)
|
||||
* Return x rounded to integral value according to the prevailing
|
||||
* rounding mode.
|
||||
* Method:
|
||||
* Using floating addition.
|
||||
* Exception:
|
||||
* Inexact flag raised if x not equal to rint(x).
|
||||
*/
|
||||
|
||||
#include "fdlibm.h"
|
||||
|
||||
#ifndef _DOUBLE_IS_32BITS
|
||||
|
||||
#ifdef __STDC__
|
||||
static const double
|
||||
#else
|
||||
static double
|
||||
#endif
|
||||
TWO52[2]={
|
||||
4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
|
||||
-4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */
|
||||
};
|
||||
|
||||
#ifdef __STDC__
|
||||
double rint(double x)
|
||||
#else
|
||||
double rint(x)
|
||||
double x;
|
||||
#endif
|
||||
{
|
||||
__int32_t j0;
|
||||
__uint32_t i,i0,i1,sx;
|
||||
double w,t;
|
||||
EXTRACT_WORDS(i0,i1,x);
|
||||
sx = (i0>>31)&1;
|
||||
j0 = (__int32_t)((i0>>20)&0x7ff)-0x3ff;
|
||||
if(j0<20) {
|
||||
if(j0<0) {
|
||||
if(((i0&0x7fffffff)|i1)==0) return x;
|
||||
i1 |= (i0&0x0fffff);
|
||||
i0 &= 0xfffe0000;
|
||||
i0 |= ((i1|-i1)>>12)&0x80000;
|
||||
SET_HIGH_WORD(x,i0);
|
||||
w = TWO52[sx]+x;
|
||||
t = w-TWO52[sx];
|
||||
GET_HIGH_WORD(i0,t);
|
||||
SET_HIGH_WORD(t,(i0&0x7fffffff)|(sx<<31));
|
||||
return t;
|
||||
} else {
|
||||
i = (__uint32_t)(0x000fffff)>>j0;
|
||||
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);
|
||||
}
|
||||
}
|
||||
} else if (j0>51) {
|
||||
if(j0==0x400) return x+x; /* inf or NaN */
|
||||
else return x; /* x is integral */
|
||||
} else {
|
||||
i = ((__uint32_t)(0xffffffff))>>(j0-20);
|
||||
if((i1&i)==0) return x; /* x is integral */
|
||||
i>>=1;
|
||||
if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20));
|
||||
}
|
||||
INSERT_WORDS(x,i0,i1);
|
||||
w = TWO52[sx]+x;
|
||||
return w-TWO52[sx];
|
||||
}
|
||||
|
||||
#endif /* _DOUBLE_IS_32BITS */
|
||||
@@ -0,0 +1,105 @@
|
||||
|
||||
/* @(#)s_scalbn.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.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
/*
|
||||
FUNCTION
|
||||
<<scalbn>>, <<scalbnf>>---scale by integer
|
||||
INDEX
|
||||
scalbn
|
||||
INDEX
|
||||
scalbnf
|
||||
|
||||
ANSI_SYNOPSIS
|
||||
#include <math.h>
|
||||
double scalbn(double <[x]>, int <[y]>);
|
||||
float scalbnf(float <[x]>, int <[y]>);
|
||||
|
||||
TRAD_SYNOPSIS
|
||||
#include <math.h>
|
||||
double scalbn(<[x]>,<[y]>)
|
||||
double <[x]>;
|
||||
int <[y]>;
|
||||
float scalbnf(<[x]>,<[y]>)
|
||||
float <[x]>;
|
||||
int <[y]>;
|
||||
|
||||
DESCRIPTION
|
||||
<<scalbn>> and <<scalbnf>> scale <[x]> by <[n]>, returning <[x]> times
|
||||
2 to the power <[n]>. The result is computed by manipulating the
|
||||
exponent, rather than by actually performing an exponentiation or
|
||||
multiplication.
|
||||
|
||||
RETURNS
|
||||
<[x]> times 2 to the power <[n]>.
|
||||
|
||||
PORTABILITY
|
||||
Neither <<scalbn>> nor <<scalbnf>> is required by ANSI C or by the System V
|
||||
Interface Definition (Issue 2).
|
||||
|
||||
*/
|
||||
|
||||
/*
|
||||
* scalbn (double x, int n)
|
||||
* scalbn(x,n) returns x* 2**n computed by exponent
|
||||
* manipulation rather than by actually performing an
|
||||
* exponentiation or a multiplication.
|
||||
*/
|
||||
|
||||
#include "fdlibm.h"
|
||||
|
||||
#ifndef _DOUBLE_IS_32BITS
|
||||
|
||||
#ifdef __STDC__
|
||||
static const double
|
||||
#else
|
||||
static double
|
||||
#endif
|
||||
two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
|
||||
twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
|
||||
huge = 1.0e+300,
|
||||
tiny = 1.0e-300;
|
||||
|
||||
#ifdef __STDC__
|
||||
double scalbn (double x, int n)
|
||||
#else
|
||||
double scalbn (x,n)
|
||||
double x; int n;
|
||||
#endif
|
||||
{
|
||||
__int32_t k;
|
||||
__uint32_t hx,lx;
|
||||
EXTRACT_WORDS(hx,lx,(unsigned int)x);
|
||||
k = (__int32_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 = (__int32_t)((hx&0x7ff00000)>>20) - 54;
|
||||
if (n< -50000) return tiny*x; /*underflow*/
|
||||
}
|
||||
if (k==0x7ff) return x+x; /* NaN or Inf */
|
||||
k = k+n;
|
||||
if (k > 0x7fe) return huge*copysign(huge,x); /* overflow */
|
||||
if (k > 0) /* normal result */
|
||||
{SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x;}
|
||||
if (k <= -54) {
|
||||
if (n > 50000) /* in case integer overflow in n+k */
|
||||
return huge*copysign(huge,x); /*overflow*/
|
||||
else return tiny*copysign(tiny,x); /*underflow*/
|
||||
}
|
||||
k += 54; /* subnormal result */
|
||||
SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20));
|
||||
return x*twom54;
|
||||
}
|
||||
|
||||
#endif /* _DOUBLE_IS_32BITS */
|
||||
Reference in New Issue
Block a user