Script 'mail_helper' called by obssrc Hello community, here is the log from the commit of package openspecfun for openSUSE:Factory checked in at 2024-07-08 19:06:56 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Comparing /work/SRC/openSUSE:Factory/openspecfun (Old) and /work/SRC/openSUSE:Factory/.openspecfun.new.2080 (New) ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Package is "openspecfun" Mon Jul 8 19:06:56 2024 rev:5 rq:1185968 version:0.5.6 Changes: -------- --- /work/SRC/openSUSE:Factory/openspecfun/openspecfun.changes 2022-12-19 14:08:11.838643554 +0100 +++ /work/SRC/openSUSE:Factory/.openspecfun.new.2080/openspecfun.changes 2024-07-08 19:07:09.412318910 +0200 @@ -1,0 +2,8 @@ +Tue Jun 25 03:46:31 UTC 2024 - ming li <[email protected]> + +- update to 0.5.6: + * Fix my_copysign bug + * Remove rem_pio2 + * Apply fixes listed in 'Remark on Algorithm 644' + +------------------------------------------------------------------- Old: ---- openspecfun-0.5.5.tar.gz New: ---- openspecfun-0.5.6.tar.gz ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Other differences: ------------------ ++++++ openspecfun.spec ++++++ --- /var/tmp/diff_new_pack.FW7gdB/_old 2024-07-08 19:07:10.340352850 +0200 +++ /var/tmp/diff_new_pack.FW7gdB/_new 2024-07-08 19:07:10.340352850 +0200 @@ -1,7 +1,7 @@ # # spec file for package openspecfun # -# Copyright (c) 2022 SUSE LLC +# Copyright (c) 2024 SUSE LLC # # All modifications and additions to the file contributed by third parties # remain the property of their copyright owners, unless otherwise agreed @@ -16,10 +16,10 @@ # -%define so_ver 1 +%define so_ver 2 %define libname lib%{name}%{so_ver} Name: openspecfun -Version: 0.5.5 +Version: 0.5.6 Release: 0 Summary: A collection of special mathematical functions License: MIT ++++++ openspecfun-0.5.5.tar.gz -> openspecfun-0.5.6.tar.gz ++++++ diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/openspecfun-0.5.5/Faddeeva/Faddeeva.cc new/openspecfun-0.5.6/Faddeeva/Faddeeva.cc --- old/openspecfun-0.5.5/Faddeeva/Faddeeva.cc 2021-06-02 13:28:45.000000000 +0200 +++ new/openspecfun-0.5.6/Faddeeva/Faddeeva.cc 2023-07-06 00:39:43.000000000 +0200 @@ -140,6 +140,7 @@ functions near the origin. Use gnulib functions if GNULIB_NAMESPACE is defined. 18 December 2012: Slight tweaks (remove recomputation of x*x in Dawson) + 12 May 2015: Bugfix for systems lacking copysign function. */ ///////////////////////////////////////////////////////////////////////// @@ -199,7 +200,7 @@ # elif defined(GNULIB_NAMESPACE) // we are using using gnulib <cmath> # define copysign GNULIB_NAMESPACE::copysign # elif (__cplusplus < 201103L) && !defined(HAVE_COPYSIGN) && !defined(__linux__) && !(defined(__APPLE__) && defined(__MACH__)) && !defined(_AIX) -static inline double my_copysign(double x, double y) { return y<0 ? -x : x; } +static inline double my_copysign(double x, double y) { return x<0 != y<0 ? -x : x; } # define copysign my_copysign # endif diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/openspecfun-0.5.5/Make.inc new/openspecfun-0.5.6/Make.inc --- old/openspecfun-0.5.5/Make.inc 2021-06-02 13:28:45.000000000 +0200 +++ new/openspecfun-0.5.6/Make.inc 2023-07-06 00:39:43.000000000 +0200 @@ -4,9 +4,9 @@ # Do not forget to bump SOMINOR when changing VERSION, # and SOMAJOR when breaking ABI in a backward-incompatible way -VERSION = 0.5.4 -SOMAJOR = 1 -SOMINOR = 4 +VERSION = 0.6.0 +SOMAJOR = 2 +SOMINOR = 1 DESTDIR = prefix = /usr/local @@ -84,7 +84,7 @@ ifeq ($(ARCH),x86_64) override ARCH := amd64 -endif +endif ifneq (,$(findstring MINGW,$(OS))) override OS=WINNT diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/openspecfun-0.5.5/Makefile new/openspecfun-0.5.6/Makefile --- old/openspecfun-0.5.5/Makefile 2021-06-02 13:28:45.000000000 +0200 +++ new/openspecfun-0.5.6/Makefile 2023-07-06 00:39:43.000000000 +0200 @@ -1,7 +1,7 @@ OPENLIBM_HOME=$(abspath .) include ./Make.inc -SUBDIRS = amos Faddeeva rem_pio2 +SUBDIRS = amos Faddeeva define INC_template TEST=test diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/openspecfun-0.5.5/amos/dgamln.f new/openspecfun-0.5.6/amos/dgamln.f --- old/openspecfun-0.5.5/amos/dgamln.f 2021-06-02 13:28:45.000000000 +0200 +++ new/openspecfun-0.5.6/amos/dgamln.f 2023-07-06 00:39:43.000000000 +0200 @@ -132,7 +132,7 @@ IERR=0 IF (Z.LE.0.0D0) GO TO 70 IF (Z.GT.101.0D0) GO TO 10 - NZ = INT(SNGL(Z)) + NZ = INT(Z) FZ = Z - FLOAT(NZ) IF (FZ.GT.0.0D0) GO TO 10 IF (NZ.GT.100) GO TO 10 diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/openspecfun-0.5.5/amos/zbknu.f new/openspecfun-0.5.6/amos/zbknu.f --- old/openspecfun-0.5.5/amos/zbknu.f 2021-06-02 13:28:45.000000000 +0200 +++ new/openspecfun-0.5.6/amos/zbknu.f 2023-07-06 00:39:43.000000000 +0200 @@ -58,7 +58,7 @@ STI = -ZI*RCAZ RZR = (STR+STR)*RCAZ RZI = (STI+STI)*RCAZ - INU = INT(SNGL(FNU+0.5D0)) + INU = INT(FNU+0.5D0) DNU = FNU - DBLE(FLOAT(INU)) IF (DABS(DNU).EQ.0.5D0) GO TO 110 DNU2 = 0.0D0 diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/openspecfun-0.5.5/rem_pio2/Make.files new/openspecfun-0.5.6/rem_pio2/Make.files --- old/openspecfun-0.5.5/rem_pio2/Make.files 2021-06-02 13:28:45.000000000 +0200 +++ new/openspecfun-0.5.6/rem_pio2/Make.files 1970-01-01 01:00:00.000000000 +0100 @@ -1 +0,0 @@ -$(CUR_SRCS) += e_rem_pio2.c e_rem_pio2f.c k_rem_pio2.c diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/openspecfun-0.5.5/rem_pio2/e_rem_pio2.c new/openspecfun-0.5.6/rem_pio2/e_rem_pio2.c --- old/openspecfun-0.5.5/rem_pio2/e_rem_pio2.c 2021-06-02 13:28:45.000000000 +0200 +++ new/openspecfun-0.5.6/rem_pio2/e_rem_pio2.c 1970-01-01 01:00:00.000000000 +0100 @@ -1,193 +0,0 @@ - -/* @(#)e_rem_pio2.c 1.4 95/01/18 */ -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunSoft, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - * - * Optimized by Bruce D. Evans. - */ - -//__FBSDID("$FreeBSD: src/lib/msun/src/e_rem_pio2.c,v 1.22 2011/06/19 17:07:58 kargl Exp $"); - -/* __ieee754_rem_pio2(x,y) - * - * return the remainder of x rem pi/2 in y[0]+y[1] - * use __kernel_rem_pio2() - */ - -#include <float.h> - -#ifdef USE_OPENLIBM -#include "openlibm.h" -#include "openlibm_complex.h" -#else -#include <math.h> -#include <complex.h> -#endif - -#include "math_private.h" - -/* - * invpio2: 53 bits of 2/pi - * pio2_1: first 33 bit of pi/2 - * pio2_1t: pi/2 - pio2_1 - * pio2_2: second 33 bit of pi/2 - * pio2_2t: pi/2 - (pio2_1+pio2_2) - * pio2_3: third 33 bit of pi/2 - * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3) - */ - -static const double -zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */ -two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */ -invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */ -pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */ -pio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */ -pio2_2 = 6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */ -pio2_2t = 2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */ -pio2_3 = 2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */ -pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */ - -#ifndef INLINE_REM_PIO2 -extern -#endif -//__inline int -int -__ieee754_rem_pio2(double x, double *y) -{ - double z,w,t,r,fn; - double tx[3],ty[2]; - int32_t e0,i,j,nx,n,ix,hx; - uint32_t low; - - GET_HIGH_WORD(hx,x); /* high word of x */ - ix = hx&0x7fffffff; -#if 0 /* Must be handled in caller. */ - if(ix<=0x3fe921fb) /* |x| ~<= pi/4 , no need for reduction */ - {y[0] = x; y[1] = 0; return 0;} -#endif - if (ix <= 0x400f6a7a) { /* |x| ~<= 5pi/4 */ - if ((ix & 0xfffff) == 0x921fb) /* |x| ~= pi/2 or 2pi/2 */ - goto medium; /* cancellation -- use medium case */ - if (ix <= 0x4002d97c) { /* |x| ~<= 3pi/4 */ - if (hx > 0) { - z = x - pio2_1; /* one round good to 85 bits */ - y[0] = z - pio2_1t; - y[1] = (z-y[0])-pio2_1t; - return 1; - } else { - z = x + pio2_1; - y[0] = z + pio2_1t; - y[1] = (z-y[0])+pio2_1t; - return -1; - } - } else { - if (hx > 0) { - z = x - 2*pio2_1; - y[0] = z - 2*pio2_1t; - y[1] = (z-y[0])-2*pio2_1t; - return 2; - } else { - z = x + 2*pio2_1; - y[0] = z + 2*pio2_1t; - y[1] = (z-y[0])+2*pio2_1t; - return -2; - } - } - } - if (ix <= 0x401c463b) { /* |x| ~<= 9pi/4 */ - if (ix <= 0x4015fdbc) { /* |x| ~<= 7pi/4 */ - if (ix == 0x4012d97c) /* |x| ~= 3pi/2 */ - goto medium; - if (hx > 0) { - z = x - 3*pio2_1; - y[0] = z - 3*pio2_1t; - y[1] = (z-y[0])-3*pio2_1t; - return 3; - } else { - z = x + 3*pio2_1; - y[0] = z + 3*pio2_1t; - y[1] = (z-y[0])+3*pio2_1t; - return -3; - } - } else { - if (ix == 0x401921fb) /* |x| ~= 4pi/2 */ - goto medium; - if (hx > 0) { - z = x - 4*pio2_1; - y[0] = z - 4*pio2_1t; - y[1] = (z-y[0])-4*pio2_1t; - return 4; - } else { - z = x + 4*pio2_1; - y[0] = z + 4*pio2_1t; - y[1] = (z-y[0])+4*pio2_1t; - return -4; - } - } - } - if(ix<0x413921fb) { /* |x| ~< 2^20*(pi/2), medium size */ -medium: - /* Use a specialized rint() to get fn. Assume round-to-nearest. */ - STRICT_ASSIGN(double,fn,x*invpio2+0x1.8p52); - fn = fn-0x1.8p52; -#ifdef HAVE_EFFICIENT_IRINT - n = irint(fn); -#else - n = (int32_t)fn; -#endif - r = x-fn*pio2_1; - w = fn*pio2_1t; /* 1st round good to 85 bit */ - { - uint32_t high; - j = ix>>20; - y[0] = r-w; - GET_HIGH_WORD(high,y[0]); - i = j-((high>>20)&0x7ff); - if(i>16) { /* 2nd iteration needed, good to 118 */ - t = r; - w = fn*pio2_2; - r = t-w; - w = fn*pio2_2t-((t-r)-w); - y[0] = r-w; - GET_HIGH_WORD(high,y[0]); - i = j-((high>>20)&0x7ff); - if(i>49) { /* 3rd iteration need, 151 bits acc */ - t = r; /* will cover all possible cases */ - w = fn*pio2_3; - r = t-w; - w = fn*pio2_3t-((t-r)-w); - y[0] = r-w; - } - } - } - y[1] = (r-y[0])-w; - return n; - } - /* - * all other (large) arguments - */ - if(ix>=0x7ff00000) { /* x is inf or NaN */ - y[0]=y[1]=x-x; return 0; - } - /* set z = scalbn(|x|,ilogb(x)-23) */ - GET_LOW_WORD(low,x); - e0 = (ix>>20)-1046; /* e0 = ilogb(z)-23; */ - INSERT_WORDS(z, ix - ((int32_t)(e0<<20)), low); - for(i=0;i<2;i++) { - tx[i] = (double)((int32_t)(z)); - z = (z-tx[i])*two24; - } - tx[2] = z; - nx = 3; - while(tx[nx-1]==zero) nx--; /* skip zero term */ - n = __kernel_rem_pio2(tx,ty,e0,nx,1); - if(hx<0) {y[0] = -ty[0]; y[1] = -ty[1]; return -n;} - y[0] = ty[0]; y[1] = ty[1]; return n; -} diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/openspecfun-0.5.5/rem_pio2/e_rem_pio2f.c new/openspecfun-0.5.6/rem_pio2/e_rem_pio2f.c --- old/openspecfun-0.5.5/rem_pio2/e_rem_pio2f.c 2021-06-02 13:28:45.000000000 +0200 +++ new/openspecfun-0.5.6/rem_pio2/e_rem_pio2f.c 1970-01-01 01:00:00.000000000 +0100 @@ -1,88 +0,0 @@ -/* e_rem_pio2f.c -- float version of e_rem_pio2.c - * Conversion to float by Ian Lance Taylor, Cygnus Support, [email protected]. - * Debugged and optimized by Bruce D. Evans. - */ - -/* - * ==================================================== - * 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. - * ==================================================== - */ - -//__FBSDID("$FreeBSD: src/lib/msun/src/e_rem_pio2f.c,v 1.32 2009/06/03 08:16:34 ed Exp $"); - -/* __ieee754_rem_pio2f(x,y) - * - * return the remainder of x rem pi/2 in *y - * use double precision for everything except passing x - * use __kernel_rem_pio2() for large x - */ - -#include <float.h> - -#ifdef USE_OPENLIBM -# include <openlibm.h> -#else -# include <math.h> -#endif -#include "math_private.h" - -/* - * invpio2: 53 bits of 2/pi - * pio2_1: first 33 bit of pi/2 - * pio2_1t: pi/2 - pio2_1 - */ - -static const double -invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */ -pio2_1 = 1.57079631090164184570e+00, /* 0x3FF921FB, 0x50000000 */ -pio2_1t = 1.58932547735281966916e-08; /* 0x3E5110b4, 0x611A6263 */ - -#ifndef INLINE_REM_PIO2F -extern -#endif -//__inline int -int -__ieee754_rem_pio2f(float x, double *y) -{ - double w,r,fn; - double tx[1],ty[1]; - float z; - int32_t e0,n,ix,hx; - - GET_FLOAT_WORD(hx,x); - ix = hx&0x7fffffff; - /* 33+53 bit pi is good enough for medium size */ - if(ix<0x4dc90fdb) { /* |x| ~< 2^28*(pi/2), medium size */ - /* Use a specialized rint() to get fn. Assume round-to-nearest. */ - STRICT_ASSIGN(double,fn,x*invpio2+0x1.8p52); - fn = fn-0x1.8p52; -#ifdef HAVE_EFFICIENT_IRINT - n = irint(fn); -#else - n = (int32_t)fn; -#endif - r = x-fn*pio2_1; - w = fn*pio2_1t; - *y = r-w; - return n; - } - /* - * all other (large) arguments - */ - if(ix>=0x7f800000) { /* x is inf or NaN */ - *y=x-x; return 0; - } - /* set z = scalbn(|x|,ilogb(|x|)-23) */ - e0 = (ix>>23)-150; /* e0 = ilogb(|x|)-23; */ - SET_FLOAT_WORD(z, ix - ((int32_t)(e0<<23))); - tx[0] = z; - n = __kernel_rem_pio2(tx,ty,e0,1,0); - if(hx<0) {*y = -ty[0]; return -n;} - *y = ty[0]; return n; -} diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/openspecfun-0.5.5/rem_pio2/fpmath.h new/openspecfun-0.5.6/rem_pio2/fpmath.h --- old/openspecfun-0.5.5/rem_pio2/fpmath.h 2021-06-02 13:28:45.000000000 +0200 +++ new/openspecfun-0.5.6/rem_pio2/fpmath.h 1970-01-01 01:00:00.000000000 +0100 @@ -1,127 +0,0 @@ -/*- - * Copyright (c) 2003 Mike Barcroft <[email protected]> - * Copyright (c) 2002 David Schultz <[email protected]> - * All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * - * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - * - * $FreeBSD: src/lib/libc/include/fpmath.h,v 1.4 2008/12/23 22:20:59 marcel Exp $ - */ -#ifndef _FPMATH_H_ -#define _FPMATH_H_ - -// We need to find out the byte ordering of the system. -// -// This will be stored in the following macros: -// _FPMATH_LITTLE_ENDIAN: An arbitrary constant -// _FPMATH_BIG_ENDIAN: An arbitrary constant -// _FPMATH_BYTE_ORDER: Either _FPMATH_LITTLE_ENDIAN or _FPMATH_BIG_ENDIAN - -// The GNU preprocessor defines a few things by itself -// See http://gcc.gnu.org/onlinedocs/cpp/Common-Predefined-Macros.html#Common-Predefined-Macros -#if defined(__ORDER_LITTLE_ENDIAN__) && defined(__ORDER_BIG_ENDIAN__) && defined(__BYTE_ORDER__) -# define _FPMATH_BIG_ENDIAN __ORDER_BIG_ENDIAN__ -# define _FPMATH_LITTLE_ENDIAN __ORDER_LITTLE_ENDIAN__ -# define _FPMATH_BYTE_ORDER __BYTE_ORDER__ -// If not defined by the preprocessor, let us try to do it ourselves. -// GNU C Library (GLIBC) - Linux, GNU/kFreeBSD, GNU/Hurd, etc. -#elif defined(__GLIBC__) -# include <features.h> -# include <endian.h> -# define _FPMATH_LITTLE_ENDIAN __LITTLE_ENDIAN -# define _FPMATH_BIG_ENDIAN __BIG_ENDIAN -# define _FPMATH_BYTE_ORDER __BYTE_ORDER -// OSX -#elif defined(__APPLE__) -# include <machine/endian.h> -# define _FPMATH_LITTLE_ENDIAN LITTLE_ENDIAN -# define _FPMATH_BIG_ENDIAN BIG_ENDIAN -# define _FPMATH_BYTE_ORDER BYTE_ORDER -// FreeBSD -#elif defined(__FreeBSD__) -# include <machine/endian.h> -# define _FPMATH_LITTLE_ENDIAN _LITTLE_ENDIAN -# define _FPMATH_BIG_ENDIAN _BIG_ENDIAN -# define _FPMATH_BYTE_ORDER _BYTE_ORDER -// Windows -#elif defined(_WIN32) -# define _FPMATH_LITTLE_ENDIAN 1234 -# define _FPMATH_BIG_ENDIAN 4321 -# define _FPMATH_BYTE_ORDER _FPMATH_LITTLE_ENDIAN -// Solaris -#elif defined(__sun) -# define _FPMATH_LITTLE_ENDIAN 1234 -# define _FPMATH_BIG_ENDIAN 4321 -# include <sys/isa_defs.h> -# ifdef _LITTLE_ENDIAN -# define _FPMATH_BYTE_ORDER _FPMATH_LITTLE_ENDIAN -# endif -# ifdef _BIG_ENDIAN -# define _FPMATH_BYTE_ORDER _FPMATH_BIG_ENDIAN -# endif -#endif - -#ifndef _IEEE_WORD_ORDER -#define _IEEE_WORD_ORDER _FPMATH_BYTE_ORDER -#endif - -union IEEEf2bits { - float f; - struct { -#if _FPMATH_BYTE_ORDER == _FPMATH_LITTLE_ENDIAN - unsigned int man :23; - unsigned int exp :8; - unsigned int sign :1; -#else /* _FPMATH_BIG_ENDIAN */ - unsigned int sign :1; - unsigned int exp :8; - unsigned int man :23; -#endif - } bits; -}; - -#define DBL_MANH_SIZE 20 -#define DBL_MANL_SIZE 32 - -union IEEEd2bits { - double d; - struct { -#if _FPMATH_BYTE_ORDER == _FPMATH_LITTLE_ENDIAN -#if _IEEE_WORD_ORDER == _FPMATH_LITTLE_ENDIAN - unsigned int manl :32; -#endif - unsigned int manh :20; - unsigned int exp :11; - unsigned int sign :1; -#if _IEEE_WORD_ORDER == _FPMATH_BIG_ENDIAN - unsigned int manl :32; -#endif -#else /* _FPMATH_BIG_ENDIAN */ - unsigned int sign :1; - unsigned int exp :11; - unsigned int manh :20; - unsigned int manl :32; -#endif - } bits; -}; - -#endif diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/openspecfun-0.5.5/rem_pio2/k_rem_pio2.c new/openspecfun-0.5.6/rem_pio2/k_rem_pio2.c --- old/openspecfun-0.5.5/rem_pio2/k_rem_pio2.c 2021-06-02 13:28:45.000000000 +0200 +++ new/openspecfun-0.5.6/rem_pio2/k_rem_pio2.c 1970-01-01 01:00:00.000000000 +0100 @@ -1,446 +0,0 @@ - -/* @(#)k_rem_pio2.c 1.3 95/01/18 */ -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunSoft, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -//__FBSDID("$FreeBSD: src/lib/msun/src/k_rem_pio2.c,v 1.11 2008/02/25 11:43:20 bde Exp $"); - -/* - * __kernel_rem_pio2(x,y,e0,nx,prec) - * double x[],y[]; int e0,nx,prec; - * - * __kernel_rem_pio2 return the last three digits of N with - * y = x - N*pi/2 - * so that |y| < pi/2. - * - * The method is to compute the integer (mod 8) and fraction parts of - * (2/pi)*x without doing the full multiplication. In general we - * skip the part of the product that are known to be a huge integer ( - * more accurately, = 0 mod 8 ). Thus the number of operations are - * independent of the exponent of the input. - * - * (2/pi) is represented by an array of 24-bit integers in ipio2[]. - * - * Input parameters: - * x[] The input value (must be positive) is broken into nx - * pieces of 24-bit integers in double precision format. - * x[i] will be the i-th 24 bit of x. The scaled exponent - * of x[0] is given in input parameter e0 (i.e., x[0]*2^e0 - * match x's up to 24 bits. - * - * Example of breaking a double positive z into x[0]+x[1]+x[2]: - * e0 = ilogb(z)-23 - * z = scalbn(z,-e0) - * for i = 0,1,2 - * x[i] = floor(z) - * z = (z-x[i])*2**24 - * - * - * y[] ouput result in an array of double precision numbers. - * The dimension of y[] is: - * 24-bit precision 1 - * 53-bit precision 2 - * 64-bit precision 2 - * 113-bit precision 3 - * The actual value is the sum of them. Thus for 113-bit - * precison, one may have to do something like: - * - * long double t,w,r_head, r_tail; - * t = (long double)y[2] + (long double)y[1]; - * w = (long double)y[0]; - * r_head = t+w; - * r_tail = w - (r_head - t); - * - * e0 The exponent of x[0]. Must be <= 16360 or you need to - * expand the ipio2 table. - * - * nx dimension of x[] - * - * prec an integer indicating the precision: - * 0 24 bits (single) - * 1 53 bits (double) - * 2 64 bits (extended) - * 3 113 bits (quad) - * - * External function: - * double scalbn(), floor(); - * - * - * Here is the description of some local variables: - * - * jk jk+1 is the initial number of terms of ipio2[] needed - * in the computation. The minimum and recommended value - * for jk is 3,4,4,6 for single, double, extended, and quad. - * jk+1 must be 2 larger than you might expect so that our - * recomputation test works. (Up to 24 bits in the integer - * part (the 24 bits of it that we compute) and 23 bits in - * the fraction part may be lost to cancelation before we - * recompute.) - * - * jz local integer variable indicating the number of - * terms of ipio2[] used. - * - * jx nx - 1 - * - * jv index for pointing to the suitable ipio2[] for the - * computation. In general, we want - * ( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8 - * is an integer. Thus - * e0-3-24*jv >= 0 or (e0-3)/24 >= jv - * Hence jv = max(0,(e0-3)/24). - * - * jp jp+1 is the number of terms in PIo2[] needed, jp = jk. - * - * q[] double array with integral value, representing the - * 24-bits chunk of the product of x and 2/pi. - * - * q0 the corresponding exponent of q[0]. Note that the - * exponent for q[i] would be q0-24*i. - * - * PIo2[] double precision array, obtained by cutting pi/2 - * into 24 bits chunks. - * - * f[] ipio2[] in floating point - * - * iq[] integer array by breaking up q[] in 24-bits chunk. - * - * fq[] final product of x*(2/pi) in fq[0],..,fq[jk] - * - * ih integer. If >0 it indicates q[] is >= 0.5, hence - * it also indicates the *sign* of the result. - * - */ - - -/* - * 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 <float.h> - -#ifdef USE_OPENLIBM -# include <openlibm.h> -#else -# include <math.h> -#endif -#include "math_private.h" - -static const int init_jk[] = {3,4,4,6}; /* initial value for jk */ - -/* - * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi - * - * integer array, contains the (24*i)-th to (24*i+23)-th - * bit of 2/pi after binary point. The corresponding - * floating value is - * - * ipio2[i] * 2^(-24(i+1)). - * - * NB: This table must have at least (e0-3)/24 + jk terms. - * For quad precision (e0 <= 16360, jk = 6), this is 686. - */ -static const int32_t ipio2[] = { -0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, -0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, -0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, -0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41, -0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8, -0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF, -0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, -0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, -0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, -0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, -0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B, - -#if LDBL_MAX_EXP > 1024 -#if LDBL_MAX_EXP > 16384 -#error "ipio2 table needs to be expanded" -#endif -0x47C419, 0xC367CD, 0xDCE809, 0x2A8359, 0xC4768B, 0x961CA6, -0xDDAF44, 0xD15719, 0x053EA5, 0xFF0705, 0x3F7E33, 0xE832C2, -0xDE4F98, 0x327DBB, 0xC33D26, 0xEF6B1E, 0x5EF89F, 0x3A1F35, -0xCAF27F, 0x1D87F1, 0x21907C, 0x7C246A, 0xFA6ED5, 0x772D30, -0x433B15, 0xC614B5, 0x9D19C3, 0xC2C4AD, 0x414D2C, 0x5D000C, -0x467D86, 0x2D71E3, 0x9AC69B, 0x006233, 0x7CD2B4, 0x97A7B4, -0xD55537, 0xF63ED7, 0x1810A3, 0xFC764D, 0x2A9D64, 0xABD770, -0xF87C63, 0x57B07A, 0xE71517, 0x5649C0, 0xD9D63B, 0x3884A7, -0xCB2324, 0x778AD6, 0x23545A, 0xB91F00, 0x1B0AF1, 0xDFCE19, -0xFF319F, 0x6A1E66, 0x615799, 0x47FBAC, 0xD87F7E, 0xB76522, -0x89E832, 0x60BFE6, 0xCDC4EF, 0x09366C, 0xD43F5D, 0xD7DE16, -0xDE3B58, 0x929BDE, 0x2822D2, 0xE88628, 0x4D58E2, 0x32CAC6, -0x16E308, 0xCB7DE0, 0x50C017, 0xA71DF3, 0x5BE018, 0x34132E, -0x621283, 0x014883, 0x5B8EF5, 0x7FB0AD, 0xF2E91E, 0x434A48, -0xD36710, 0xD8DDAA, 0x425FAE, 0xCE616A, 0xA4280A, 0xB499D3, -0xF2A606, 0x7F775C, 0x83C2A3, 0x883C61, 0x78738A, 0x5A8CAF, -0xBDD76F, 0x63A62D, 0xCBBFF4, 0xEF818D, 0x67C126, 0x45CA55, -0x36D9CA, 0xD2A828, 0x8D61C2, 0x77C912, 0x142604, 0x9B4612, -0xC459C4, 0x44C5C8, 0x91B24D, 0xF31700, 0xAD43D4, 0xE54929, -0x10D5FD, 0xFCBE00, 0xCC941E, 0xEECE70, 0xF53E13, 0x80F1EC, -0xC3E7B3, 0x28F8C7, 0x940593, 0x3E71C1, 0xB3092E, 0xF3450B, -0x9C1288, 0x7B20AB, 0x9FB52E, 0xC29247, 0x2F327B, 0x6D550C, -0x90A772, 0x1FE76B, 0x96CB31, 0x4A1679, 0xE27941, 0x89DFF4, -0x9794E8, 0x84E6E2, 0x973199, 0x6BED88, 0x365F5F, 0x0EFDBB, -0xB49A48, 0x6CA467, 0x427271, 0x325D8D, 0xB8159F, 0x09E5BC, -0x25318D, 0x3974F7, 0x1C0530, 0x010C0D, 0x68084B, 0x58EE2C, -0x90AA47, 0x02E774, 0x24D6BD, 0xA67DF7, 0x72486E, 0xEF169F, -0xA6948E, 0xF691B4, 0x5153D1, 0xF20ACF, 0x339820, 0x7E4BF5, -0x6863B2, 0x5F3EDD, 0x035D40, 0x7F8985, 0x295255, 0xC06437, -0x10D86D, 0x324832, 0x754C5B, 0xD4714E, 0x6E5445, 0xC1090B, -0x69F52A, 0xD56614, 0x9D0727, 0x50045D, 0xDB3BB4, 0xC576EA, -0x17F987, 0x7D6B49, 0xBA271D, 0x296996, 0xACCCC6, 0x5414AD, -0x6AE290, 0x89D988, 0x50722C, 0xBEA404, 0x940777, 0x7030F3, -0x27FC00, 0xA871EA, 0x49C266, 0x3DE064, 0x83DD97, 0x973FA3, -0xFD9443, 0x8C860D, 0xDE4131, 0x9D3992, 0x8C70DD, 0xE7B717, -0x3BDF08, 0x2B3715, 0xA0805C, 0x93805A, 0x921110, 0xD8E80F, -0xAF806C, 0x4BFFDB, 0x0F9038, 0x761859, 0x15A562, 0xBBCB61, -0xB989C7, 0xBD4010, 0x04F2D2, 0x277549, 0xF6B6EB, 0xBB22DB, -0xAA140A, 0x2F2689, 0x768364, 0x333B09, 0x1A940E, 0xAA3A51, -0xC2A31D, 0xAEEDAF, 0x12265C, 0x4DC26D, 0x9C7A2D, 0x9756C0, -0x833F03, 0xF6F009, 0x8C402B, 0x99316D, 0x07B439, 0x15200C, -0x5BC3D8, 0xC492F5, 0x4BADC6, 0xA5CA4E, 0xCD37A7, 0x36A9E6, -0x9492AB, 0x6842DD, 0xDE6319, 0xEF8C76, 0x528B68, 0x37DBFC, -0xABA1AE, 0x3115DF, 0xA1AE00, 0xDAFB0C, 0x664D64, 0xB705ED, -0x306529, 0xBF5657, 0x3AFF47, 0xB9F96A, 0xF3BE75, 0xDF9328, -0x3080AB, 0xF68C66, 0x15CB04, 0x0622FA, 0x1DE4D9, 0xA4B33D, -0x8F1B57, 0x09CD36, 0xE9424E, 0xA4BE13, 0xB52333, 0x1AAAF0, -0xA8654F, 0xA5C1D2, 0x0F3F0B, 0xCD785B, 0x76F923, 0x048B7B, -0x721789, 0x53A6C6, 0xE26E6F, 0x00EBEF, 0x584A9B, 0xB7DAC4, -0xBA66AA, 0xCFCF76, 0x1D02D1, 0x2DF1B1, 0xC1998C, 0x77ADC3, -0xDA4886, 0xA05DF7, 0xF480C6, 0x2FF0AC, 0x9AECDD, 0xBC5C3F, -0x6DDED0, 0x1FC790, 0xB6DB2A, 0x3A25A3, 0x9AAF00, 0x9353AD, -0x0457B6, 0xB42D29, 0x7E804B, 0xA707DA, 0x0EAA76, 0xA1597B, -0x2A1216, 0x2DB7DC, 0xFDE5FA, 0xFEDB89, 0xFDBE89, 0x6C76E4, -0xFCA906, 0x70803E, 0x156E85, 0xFF87FD, 0x073E28, 0x336761, -0x86182A, 0xEABD4D, 0xAFE7B3, 0x6E6D8F, 0x396795, 0x5BBF31, -0x48D784, 0x16DF30, 0x432DC7, 0x356125, 0xCE70C9, 0xB8CB30, -0xFD6CBF, 0xA200A4, 0xE46C05, 0xA0DD5A, 0x476F21, 0xD21262, -0x845CB9, 0x496170, 0xE0566B, 0x015299, 0x375550, 0xB7D51E, -0xC4F133, 0x5F6E13, 0xE4305D, 0xA92E85, 0xC3B21D, 0x3632A1, -0xA4B708, 0xD4B1EA, 0x21F716, 0xE4698F, 0x77FF27, 0x80030C, -0x2D408D, 0xA0CD4F, 0x99A520, 0xD3A2B3, 0x0A5D2F, 0x42F9B4, -0xCBDA11, 0xD0BE7D, 0xC1DB9B, 0xBD17AB, 0x81A2CA, 0x5C6A08, -0x17552E, 0x550027, 0xF0147F, 0x8607E1, 0x640B14, 0x8D4196, -0xDEBE87, 0x2AFDDA, 0xB6256B, 0x34897B, 0xFEF305, 0x9EBFB9, -0x4F6A68, 0xA82A4A, 0x5AC44F, 0xBCF82D, 0x985AD7, 0x95C7F4, -0x8D4D0D, 0xA63A20, 0x5F57A4, 0xB13F14, 0x953880, 0x0120CC, -0x86DD71, 0xB6DEC9, 0xF560BF, 0x11654D, 0x6B0701, 0xACB08C, -0xD0C0B2, 0x485551, 0x0EFB1E, 0xC37295, 0x3B06A3, 0x3540C0, -0x7BDC06, 0xCC45E0, 0xFA294E, 0xC8CAD6, 0x41F3E8, 0xDE647C, -0xD8649B, 0x31BED9, 0xC397A4, 0xD45877, 0xC5E369, 0x13DAF0, -0x3C3ABA, 0x461846, 0x5F7555, 0xF5BDD2, 0xC6926E, 0x5D2EAC, -0xED440E, 0x423E1C, 0x87C461, 0xE9FD29, 0xF3D6E7, 0xCA7C22, -0x35916F, 0xC5E008, 0x8DD7FF, 0xE26A6E, 0xC6FDB0, 0xC10893, -0x745D7C, 0xB2AD6B, 0x9D6ECD, 0x7B723E, 0x6A11C6, 0xA9CFF7, -0xDF7329, 0xBAC9B5, 0x5100B7, 0x0DB2E2, 0x24BA74, 0x607DE5, -0x8AD874, 0x2C150D, 0x0C1881, 0x94667E, 0x162901, 0x767A9F, -0xBEFDFD, 0xEF4556, 0x367ED9, 0x13D9EC, 0xB9BA8B, 0xFC97C4, -0x27A831, 0xC36EF1, 0x36C594, 0x56A8D8, 0xB5A8B4, 0x0ECCCF, -0x2D8912, 0x34576F, 0x89562C, 0xE3CE99, 0xB920D6, 0xAA5E6B, -0x9C2A3E, 0xCC5F11, 0x4A0BFD, 0xFBF4E1, 0x6D3B8E, 0x2C86E2, -0x84D4E9, 0xA9B4FC, 0xD1EEEF, 0xC9352E, 0x61392F, 0x442138, -0xC8D91B, 0x0AFC81, 0x6A4AFB, 0xD81C2F, 0x84B453, 0x8C994E, -0xCC2254, 0xDC552A, 0xD6C6C0, 0x96190B, 0xB8701A, 0x649569, -0x605A26, 0xEE523F, 0x0F117F, 0x11B5F4, 0xF5CBFC, 0x2DBC34, -0xEEBC34, 0xCC5DE8, 0x605EDD, 0x9B8E67, 0xEF3392, 0xB817C9, -0x9B5861, 0xBC57E1, 0xC68351, 0x103ED8, 0x4871DD, 0xDD1C2D, -0xA118AF, 0x462C21, 0xD7F359, 0x987AD9, 0xC0549E, 0xFA864F, -0xFC0656, 0xAE79E5, 0x362289, 0x22AD38, 0xDC9367, 0xAAE855, -0x382682, 0x9BE7CA, 0xA40D51, 0xB13399, 0x0ED7A9, 0x480569, -0xF0B265, 0xA7887F, 0x974C88, 0x36D1F9, 0xB39221, 0x4A827B, -0x21CF98, 0xDC9F40, 0x5547DC, 0x3A74E1, 0x42EB67, 0xDF9DFE, -0x5FD45E, 0xA4677B, 0x7AACBA, 0xA2F655, 0x23882B, 0x55BA41, -0x086E59, 0x862A21, 0x834739, 0xE6E389, 0xD49EE5, 0x40FB49, -0xE956FF, 0xCA0F1C, 0x8A59C5, 0x2BFA94, 0xC5C1D3, 0xCFC50F, -0xAE5ADB, 0x86C547, 0x624385, 0x3B8621, 0x94792C, 0x876110, -0x7B4C2A, 0x1A2C80, 0x12BF43, 0x902688, 0x893C78, 0xE4C4A8, -0x7BDBE5, 0xC23AC4, 0xEAF426, 0x8A67F7, 0xBF920D, 0x2BA365, -0xB1933D, 0x0B7CBD, 0xDC51A4, 0x63DD27, 0xDDE169, 0x19949A, -0x9529A8, 0x28CE68, 0xB4ED09, 0x209F44, 0xCA984E, 0x638270, -0x237C7E, 0x32B90F, 0x8EF5A7, 0xE75614, 0x08F121, 0x2A9DB5, -0x4D7E6F, 0x5119A5, 0xABF9B5, 0xD6DF82, 0x61DD96, 0x023616, -0x9F3AC4, 0xA1A283, 0x6DED72, 0x7A8D39, 0xA9B882, 0x5C326B, -0x5B2746, 0xED3400, 0x7700D2, 0x55F4FC, 0x4D5901, 0x8071E0, -#endif - -}; - -static const double PIo2[] = { - 1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */ - 7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */ - 5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */ - 3.28200341580791294123e-22, /* 0x3B78CC51, 0x60000000 */ - 1.27065575308067607349e-29, /* 0x39F01B83, 0x80000000 */ - 1.22933308981111328932e-36, /* 0x387A2520, 0x40000000 */ - 2.73370053816464559624e-44, /* 0x36E38222, 0x80000000 */ - 2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */ -}; - -static const double -zero = 0.0, -one = 1.0, -two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */ -twon24 = 5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */ - -int -__kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec) -{ - int32_t jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih; - double z,fw,f[20],fq[20],q[20]; - - /* initialize jk*/ - jk = init_jk[prec]; - jp = jk; - - /* determine jx,jv,q0, note that 3>q0 */ - jx = nx-1; - jv = (e0-3)/24; if(jv<0) jv=0; - q0 = e0-24*(jv+1); - - /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */ - j = jv-jx; m = jx+jk; - for(i=0;i<=m;i++,j++) f[i] = (j<0)? zero : (double) ipio2[j]; - - /* compute q[0],q[1],...q[jk] */ - for (i=0;i<=jk;i++) { - for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j]; q[i] = fw; - } - - jz = jk; -recompute: - /* distill q[] into iq[] reversingly */ - for(i=0,j=jz,z=q[jz];j>0;i++,j--) { - fw = (double)((int32_t)(twon24* z)); - iq[i] = (int32_t)(z-two24*fw); - z = q[j-1]+fw; - } - - /* compute n */ - z = scalbn(z,q0); /* actual value of z */ - z -= 8.0*floor(z*0.125); /* trim off integer >= 8 */ - n = (int32_t) z; - z -= (double)n; - ih = 0; - if(q0>0) { /* need iq[jz-1] to determine n */ - i = (iq[jz-1]>>(24-q0)); n += i; - iq[jz-1] -= i<<(24-q0); - ih = iq[jz-1]>>(23-q0); - } - else if(q0==0) ih = iq[jz-1]>>23; - else if(z>=0.5) ih=2; - - if(ih>0) { /* q > 0.5 */ - n += 1; carry = 0; - for(i=0;i<jz ;i++) { /* compute 1-q */ - j = iq[i]; - if(carry==0) { - if(j!=0) { - carry = 1; iq[i] = 0x1000000- j; - } - } else iq[i] = 0xffffff - j; - } - if(q0>0) { /* rare case: chance is 1 in 12 */ - switch(q0) { - case 1: - iq[jz-1] &= 0x7fffff; break; - case 2: - iq[jz-1] &= 0x3fffff; break; - } - } - if(ih==2) { - z = one - z; - if(carry!=0) z -= scalbn(one,q0); - } - } - - /* check if recomputation is needed */ - if(z==zero) { - j = 0; - for (i=jz-1;i>=jk;i--) j |= iq[i]; - if(j==0) { /* need recomputation */ - for(k=1;iq[jk-k]==0;k++); /* k = no. of terms needed */ - - for(i=jz+1;i<=jz+k;i++) { /* add q[jz+1] to q[jz+k] */ - f[jx+i] = (double) ipio2[jv+i]; - for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j]; - q[i] = fw; - } - jz += k; - goto recompute; - } - } - - /* chop off zero terms */ - if(z==0.0) { - jz -= 1; q0 -= 24; - while(iq[jz]==0) { jz--; q0-=24;} - } else { /* break z into 24-bit if necessary */ - z = scalbn(z,-q0); - if(z>=two24) { - fw = (double)((int32_t)(twon24*z)); - iq[jz] = (int32_t)(z-two24*fw); - jz += 1; q0 += 24; - iq[jz] = (int32_t) fw; - } else iq[jz] = (int32_t) z ; - } - - /* convert integer "bit" chunk to floating-point value */ - fw = scalbn(one,q0); - for(i=jz;i>=0;i--) { - q[i] = fw*(double)iq[i]; fw*=twon24; - } - - /* compute PIo2[0,...,jp]*q[jz,...,0] */ - for(i=jz;i>=0;i--) { - for(fw=0.0,k=0;k<=jp&&k<=jz-i;k++) fw += PIo2[k]*q[i+k]; - fq[jz-i] = fw; - } - - /* compress fq[] into y[] */ - switch(prec) { - case 0: - fw = 0.0; - for (i=jz;i>=0;i--) fw += fq[i]; - y[0] = (ih==0)? fw: -fw; - break; - case 1: - case 2: - fw = 0.0; - for (i=jz;i>=0;i--) fw += fq[i]; - STRICT_ASSIGN(double,fw,fw); - y[0] = (ih==0)? fw: -fw; - fw = fq[0]-fw; - for (i=1;i<=jz;i++) fw += fq[i]; - y[1] = (ih==0)? fw: -fw; - break; - case 3: /* painful */ - for (i=jz;i>0;i--) { - fw = fq[i-1]+fq[i]; - fq[i] += fq[i-1]-fw; - fq[i-1] = fw; - } - for (i=jz;i>1;i--) { - fw = fq[i-1]+fq[i]; - fq[i] += fq[i-1]-fw; - fq[i-1] = fw; - } - for (fw=0.0,i=jz;i>=2;i--) fw += fq[i]; - if(ih==0) { - y[0] = fq[0]; y[1] = fq[1]; y[2] = fw; - } else { - y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw; - } - } - return n&7; -} diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/openspecfun-0.5.5/rem_pio2/math_private.h new/openspecfun-0.5.6/rem_pio2/math_private.h --- old/openspecfun-0.5.5/rem_pio2/math_private.h 2021-06-02 13:28:45.000000000 +0200 +++ new/openspecfun-0.5.6/rem_pio2/math_private.h 1970-01-01 01:00:00.000000000 +0100 @@ -1,375 +0,0 @@ -/* - * ==================================================== - * 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. - * ==================================================== - */ - -/* - * from: @(#)fdlibm.h 5.1 93/09/24 - * $FreeBSD: src/lib/msun/src/math_private.h,v 1.34 2011/10/21 06:27:56 das Exp $ - */ - -#ifndef _MATH_PRIVATE_H_ -#define _MATH_PRIVATE_H_ - -#include "types-compat.h" -#include "fpmath.h" -#include <complex.h> -#include <stdint.h> - -/* - * The original fdlibm 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. - */ - -/* - * A union which permits us to convert between a double and two 32 bit - * ints. - */ - -#if _IEEE_WORD_ORDER == _FPMATH_BIG_ENDIAN - -typedef union -{ - double value; - struct - { - uint32_t msw; - uint32_t lsw; - } parts; - struct - { - uint64_t w; - } xparts; -} ieee_double_shape_type; - -#endif - -#if _IEEE_WORD_ORDER == _FPMATH_LITTLE_ENDIAN - -typedef union -{ - double value; - struct - { - uint32_t lsw; - uint32_t msw; - } parts; - struct - { - uint64_t w; - } xparts; -} 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 a 64-bit int from a double. */ -#define EXTRACT_WORD64(ix,d) \ -do { \ - ieee_double_shape_type ew_u; \ - ew_u.value = (d); \ - (ix) = ew_u.xparts.w; \ -} 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 a double from a 64-bit int. */ -#define INSERT_WORD64(d,ix) \ -do { \ - ieee_double_shape_type iw_u; \ - iw_u.xparts.w = (ix); \ - (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; - /* FIXME: Assumes 32 bit int. */ - unsigned int 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) - -/* Get expsign as a 16 bit int from a long double. */ - -#define GET_LDBL_EXPSIGN(i,d) \ -do { \ - union IEEEl2bits ge_u; \ - ge_u.e = (d); \ - (i) = ge_u.xbits.expsign; \ -} while (0) - -/* Set expsign of a long double from a 16 bit int. */ - -#define SET_LDBL_EXPSIGN(d,v) \ -do { \ - union IEEEl2bits se_u; \ - se_u.e = (d); \ - se_u.xbits.expsign = (v); \ - (d) = se_u.e; \ -} while (0) - - -//VBS -#define STRICT_ASSIGN(type, lval, rval) ((lval) = (rval)) - -/* VBS -#ifdef FLT_EVAL_METHOD -// Attempt to get strict C99 semantics for assignment with non-C99 compilers. -#if FLT_EVAL_METHOD == 0 || __GNUC__ == 0 -#define STRICT_ASSIGN(type, lval, rval) ((lval) = (rval)) -#else -#define STRICT_ASSIGN(type, lval, rval) do { \ - volatile type __lval; \ - \ - if (sizeof(type) >= sizeof(double)) \ - (lval) = (rval); \ - else { \ - __lval = (rval); \ - (lval) = __lval; \ - } \ -} while (0) -#endif -#endif -*/ - -/* - * Common routine to process the arguments to nan(), nanf(), and nanl(). - */ -void _scan_nan(uint32_t *__words, int __num_words, const char *__s); - -#ifdef __GNUCLIKE_ASM - -/* Asm versions of some functions. */ - -#ifdef __amd64__ -static __inline int -irint(double x) -{ - int n; - - __asm__("cvtsd2si %1,%0" : "=r" (n) : "x" (x)); - return (n); -} -#define HAVE_EFFICIENT_IRINT -#endif - -#ifdef __i386__ -static __inline int -irint(double x) -{ - int n; - - __asm__("fistl %0" : "=m" (n) : "t" (x)); - return (n); -} -#define HAVE_EFFICIENT_IRINT -#endif - -#endif /* __GNUCLIKE_ASM */ - -/* - * ieee style elementary functions - * - * We rename functions here to improve other sources' diffability - * against fdlibm. - */ -#define __ieee754_sqrt sqrt -#define __ieee754_acos acos -#define __ieee754_acosh acosh -#define __ieee754_log log -#define __ieee754_log2 log2 -#define __ieee754_atanh atanh -#define __ieee754_asin asin -#define __ieee754_atan2 atan2 -#define __ieee754_exp exp -#define __ieee754_cosh cosh -#define __ieee754_fmod fmod -#define __ieee754_pow pow -#define __ieee754_lgamma lgamma -#define __ieee754_gamma gamma -#define __ieee754_lgamma_r lgamma_r -#define __ieee754_gamma_r gamma_r -#define __ieee754_log10 log10 -#define __ieee754_sinh sinh -#define __ieee754_hypot hypot -#define __ieee754_j0 j0 -#define __ieee754_j1 j1 -#define __ieee754_y0 y0 -#define __ieee754_y1 y1 -#define __ieee754_jn jn -#define __ieee754_yn yn -#define __ieee754_remainder remainder -#define __ieee754_scalb scalb -#define __ieee754_sqrtf sqrtf -#define __ieee754_acosf acosf -#define __ieee754_acoshf acoshf -#define __ieee754_logf logf -#define __ieee754_atanhf atanhf -#define __ieee754_asinf asinf -#define __ieee754_atan2f atan2f -#define __ieee754_expf expf -#define __ieee754_coshf coshf -#define __ieee754_fmodf fmodf -#define __ieee754_powf powf -#define __ieee754_lgammaf lgammaf -#define __ieee754_gammaf gammaf -#define __ieee754_lgammaf_r lgammaf_r -#define __ieee754_gammaf_r gammaf_r -#define __ieee754_log10f log10f -#define __ieee754_log2f log2f -#define __ieee754_sinhf sinhf -#define __ieee754_hypotf hypotf -#define __ieee754_j0f j0f -#define __ieee754_j1f j1f -#define __ieee754_y0f y0f -#define __ieee754_y1f y1f -#define __ieee754_jnf jnf -#define __ieee754_ynf ynf -#define __ieee754_remainderf remainderf -#define __ieee754_scalbf scalbf - -/* fdlibm kernel function */ -int __kernel_rem_pio2(double*,double*,int,int,int); - -/* double precision kernel functions */ -#ifdef INLINE_REM_PIO2 -__inline -#endif -int __ieee754_rem_pio2(double,double*); -double __kernel_sin(double,double,int); -double __kernel_cos(double,double); -double __kernel_tan(double,double,int); -double __ldexp_exp(double,int); -#ifdef _COMPLEX_H -double complex __ldexp_cexp(double complex,int); -#endif - -/* float precision kernel functions */ -#ifdef INLINE_REM_PIO2F -__inline -#endif -int __ieee754_rem_pio2f(float,double*); -#ifdef INLINE_KERNEL_SINDF -__inline -#endif -float __kernel_sindf(double); -#ifdef INLINE_KERNEL_COSDF -__inline -#endif -float __kernel_cosdf(double); -#ifdef INLINE_KERNEL_TANDF -__inline -#endif -float __kernel_tandf(double,int); -float __ldexp_expf(float,int); -#ifdef _COMPLEX_H -float complex __ldexp_cexpf(float complex,int); -#endif - -/* long double precision kernel functions */ -long double __kernel_sinl(long double, long double, int); -long double __kernel_cosl(long double, long double); -long double __kernel_tanl(long double, long double, int); - -#ifdef _WIN32 -# ifdef IMPORT_EXPORTS -# define DLLEXPORT __declspec(dllimport) -# else -# define DLLEXPORT __declspec(dllexport) -# endif -#else -#define DLLEXPORT __attribute__ ((visibility("default"))) -#endif - -#endif /* !_MATH_PRIVATE_H_ */ diff -urN '--exclude=CVS' '--exclude=.cvsignore' '--exclude=.svn' '--exclude=.svnignore' old/openspecfun-0.5.5/rem_pio2/types-compat.h new/openspecfun-0.5.6/rem_pio2/types-compat.h --- old/openspecfun-0.5.5/rem_pio2/types-compat.h 2021-06-02 13:28:45.000000000 +0200 +++ new/openspecfun-0.5.6/rem_pio2/types-compat.h 1970-01-01 01:00:00.000000000 +0100 @@ -1,23 +0,0 @@ -#ifndef _TYPES_COMPAT_H_ -#define _TYPES_COMPAT_H_ - -#include <sys/types.h> -#include <limits.h> -#include <stdint.h> - -#ifdef __linux__ -/* Not sure what to do about __pure2 on linux */ -#define __pure2 -#endif - -#if defined(_WIN32) || defined(__sun) -/* Not sure what to do about __pure2 on windows */ -#define __pure2 -typedef uint8_t u_int8_t; -typedef uint16_t u_int16_t; -typedef uint32_t u_int32_t; -typedef uint64_t u_int64_t; -#endif - - -#endif
