diff options
Diffstat (limited to 'sysdeps/powerpc')
46 files changed, 3119 insertions, 63 deletions
diff --git a/sysdeps/powerpc/Dist b/sysdeps/powerpc/Dist index ae16c0f8fb..282cf1394e 100644 --- a/sysdeps/powerpc/Dist +++ b/sysdeps/powerpc/Dist @@ -1,2 +1,3 @@ fenv_const.c fenv_libc.h +quad_float.h diff --git a/sysdeps/powerpc/Makefile b/sysdeps/powerpc/Makefile index 4100901d62..0a50956640 100644 --- a/sysdeps/powerpc/Makefile +++ b/sysdeps/powerpc/Makefile @@ -1,3 +1,13 @@ ifeq ($(subdir),math) -libm-support += fenv_const +libm-support += fenv_const fe_nomask t_sqrt + +# These routines have not been tested, so it's probable they don't work; +# and they don't provide the complete list of FP routines. So there's +# no point in compiling them. +#sysdep_routines += q_feq q_fne q_utoq q_dtoq q_itoq q_stoq q_neg q_ulltoq \ +# q_lltoq q_qtou q_qtoi q_qtoull q_qtoll q_qtos + +tests += test-arith test-arithf +LDLIBS-test-arith = libm +LDLIBS-test-arithf = libm endif diff --git a/sysdeps/powerpc/bits/fenv.h b/sysdeps/powerpc/bits/fenv.h index 08d998e118..8c26a2548f 100644 --- a/sysdeps/powerpc/bits/fenv.h +++ b/sysdeps/powerpc/bits/fenv.h @@ -37,9 +37,10 @@ enum /* ... except for FE_INVALID, for which we use bit 31. FE_INVALID actually corresponds to bits 7 through 12 and 21 through 23 in the FPSCR, but we can't use that because the current draft - says that it must be a power of 2. Instead we use bit 24 which - is the enable bit for all the FE_INVALID exceptions. */ - FE_INVALID = 1 << 31-24, + says that it must be a power of 2. Instead we use bit 2 which + is the summary bit for all the FE_INVALID exceptions, which + kind of makes sense. */ + FE_INVALID = 1 << 31-2, #define FE_INVALID FE_INVALID #ifdef __USE_GNU @@ -69,19 +70,22 @@ enum FE_INVALID_IMZ = 1 << 31-11, #define FE_INVALID_IMZ FE_INVALID_IMZ - /* Comparison with NaN or SNaN. */ + /* Comparison with NaN or SNaN. */ FE_INVALID_COMPARE = 1 << 31-12, #define FE_INVALID_COMPARE FE_INVALID_COMPARE - /* Invalid operation flag for software (not set by hardware). */ + /* Invalid operation flag for software (not set by hardware). */ + /* Note that some chips don't have this implemented, presumably + because no-one expected anyone to write software for them %-). */ FE_INVALID_SOFTWARE = 1 << 31-21, #define FE_INVALID_SOFTWARE FE_INVALID_SOFTWARE - /* Square root of negative number (including -Inf). */ + /* Square root of negative number (including -Inf). */ + /* Note that some chips don't have this implemented. */ FE_INVALID_SQRT = 1 << 31-22, #define FE_INVALID_SQRT FE_INVALID_SQRT - /* Conversion-to-integer of a NaN or a number too large or too small. */ + /* Conversion-to-integer of a NaN or a number too large or too small. */ FE_INVALID_INTEGER_CONVERSION = 1 << 31-23, #define FE_INVALID_INTEGER_CONVERSION FE_INVALID_INTEGER_CONVERSION @@ -122,7 +126,53 @@ extern const fenv_t __fe_dfl_env; #define FE_DFL_ENV (&__fe_dfl_env) #ifdef __USE_GNU -/* Floating-point environment where none of the exceptions are masked. */ -extern const fenv_t __fe_nomask_env; -# define FE_NOMASK_ENV (&__fe_nomask_env) +/* Floating-point environment where all exceptions are enabled. Note that + this is not sufficient to give you SIGFPE. */ +extern const fenv_t __fe_enabled_env; +# define FE_ENABLED_ENV (&__fe_enabled_env) + +/* Floating-point environment with (processor-dependent) non-IEEE floating + point. */ +extern const fenv_t __fe_nonieee_env; +# define FE_NONIEEE_ENV (&__fe_nonieee_env) + +/* Floating-point environment with all exceptions enabled. Note that + just evaluating this value will set the processor into 'FPU + exceptions imprecise recoverable' mode, which may cause a significant + performance penalty (but have no other visible effect). */ +extern const fenv_t *__fe_nomask_env __P ((void)); +# define FE_NOMASK_ENV (__fe_nomask_env ()) #endif + +#ifdef __OPTIMIZE__ +/* Inline definition for fegetround. */ +# define fegetround() \ + (__extension__ ({ int __fegetround_result; \ + __asm__ ("mcrfs 7,7 ; mfcr %0" \ + : "=r"(__fegetround_result) : : "cr7"); \ + __fegetround_result & 3; })) + +/* Inline definition for feraiseexcept. */ +# define feraiseexcept(__excepts) \ + (__extension__ ({ if (__builtin_constant_p (__excepts) \ + && ((__excepts) & -(__excepts)) == 0 \ + && (__excepts) != FE_INVALID) { \ + if ((__excepts) != 0) \ + __asm__ __volatile__ \ + ("mtfsb1 %0" \ + : : "i"(32 - __builtin_ffs (__excepts))); \ + } else \ + (feraiseexcept) (__excepts); })) + +/* Inline definition for feclearexcept. */ +# define feclearexcept(__excepts) \ + (__extension__ ({ if (__builtin_constant_p (__excepts) \ + && ((__excepts) & -(__excepts)) == 0 \ + && (__excepts) != FE_INVALID) { \ + if ((__excepts) != 0) \ + __asm__ __volatile__ \ + ("mtfsb0 %0" \ + : : "i"(32 - __builtin_ffs (__excepts))); \ + } else \ + (feclearexcept) (__excepts); })) +#endif /* __OPTIMIZE__ */ diff --git a/sysdeps/powerpc/bits/mathdef.h b/sysdeps/powerpc/bits/mathdef.h index c0e6caa843..9f91863333 100644 --- a/sysdeps/powerpc/bits/mathdef.h +++ b/sysdeps/powerpc/bits/mathdef.h @@ -30,7 +30,7 @@ #ifdef __GNUC__ #if __STDC__ == 1 -/* In GNU or ANSI mode, gcc leaves `float' expressions as-is, I think. */ +/* In GNU or ANSI mode, gcc leaves `float' expressions as-is. */ typedef float float_t; /* `float' expressions are evaluated as `float'. */ typedef double double_t; /* `double' expressions are evaluated as @@ -70,3 +70,7 @@ typedef double double_t; #define INFINITY HUGE_VAL #endif + +/* The values returned by `ilogb' for 0 and NaN respectively. */ +#define FP_ILOGB0 0x80000001 +#define FP_ILOGBNAN 0x7fffffff diff --git a/sysdeps/powerpc/dl-machine.h b/sysdeps/powerpc/dl-machine.h index c9d1d7332a..74436d15c1 100644 --- a/sysdeps/powerpc/dl-machine.h +++ b/sysdeps/powerpc/dl-machine.h @@ -57,10 +57,11 @@ #define OPCODE_LI(rd,simm) OPCODE_ADDI(rd,0,simm) #define OPCODE_SLWI(ra,rs,sh) OPCODE_RLWINM(ra,rs,sh,0,31-sh) -#define PPC_DCBST(where) asm __volatile__ ("dcbst 0,%0" : : "r"(where)) -#define PPC_SYNC asm __volatile__ ("sync") -#define PPC_ISYNC asm __volatile__ ("sync; isync") -#define PPC_ICBI(where) asm __volatile__ ("icbi 0,%0" : : "r"(where)) +#define PPC_DCBST(where) asm volatile ("dcbst 0,%0" : : "r"(where)) +#define PPC_SYNC asm volatile ("sync") +#define PPC_ISYNC asm volatile ("sync; isync") +#define PPC_ICBI(where) asm volatile ("icbi 0,%0" : : "r"(where)) +#define PPC_DIE asm volatile ("tweq 0,0") /* Use this when you've modified some code, but it won't be in the instruction fetch queue (or when it doesn't matter if it is). */ @@ -147,9 +148,9 @@ elf_machine_load_address (void) /* The PLT uses Elf32_Rela relocs. */ #define elf_machine_relplt elf_machine_rela - /* This code is used in dl-runtime.c to call the `fixup' function - and then redirect to the address it returns. It is called - from code built in the PLT by elf_machine_runtime_setup. */ +/* This code is used in dl-runtime.c to call the `fixup' function + and then redirect to the address it returns. It is called + from code built in the PLT by elf_machine_runtime_setup. */ #define ELF_MACHINE_RUNTIME_TRAMPOLINE asm ("\ .section \".text\" .align 2 @@ -511,11 +512,14 @@ elf_machine_lazy_rel (struct link_map *map, const Elf32_Rela *reloc) static inline void elf_machine_rela (struct link_map *map, const Elf32_Rela *reloc, const Elf32_Sym *sym, const struct r_found_version *version, - Elf32_addr *const reloc_addr) + Elf32_Addr *const reloc_addr) { +#ifndef RTLD_BOOTSTRAP const Elf32_Sym *const refsym = sym; +#endif Elf32_Word loadbase, finaladdr; const int rinfo = ELF32_R_TYPE (reloc->r_info); + extern char **_dl_argv; if (rinfo == R_PPC_NONE) return; @@ -598,7 +602,6 @@ elf_machine_rela (struct link_map *map, const Elf32_Rela *reloc, if (sym->st_size > refsym->st_size || (_dl_verbose && sym->st_size < refsym->st_size)) { - extern char **_dl_argv; const char *strtab; strtab = ((void *) map->l_addr @@ -631,11 +634,31 @@ elf_machine_rela (struct link_map *map, const Elf32_Rela *reloc, plt = (Elf32_Word *)((char *)map->l_addr + map->l_info[DT_PLTGOT]->d_un.d_val); index = (reloc_addr - plt - PLT_INITIAL_ENTRY_WORDS)/2; - if (index >= PLT_DOUBLE_SIZE) { /* Slots greater than or equal to 2^13 have 4 words available instead of two. */ + /* FIXME: There are some possible race conditions in this code, + when called from 'fixup'. + + 1) Suppose that a lazy PLT entry is executing, a + context switch between threads (or a signal) occurs, + and the new thread or signal handler calls the same + lazy PLT entry. Then the PLT entry would be changed + while it's being run, which will cause a segfault + (almost always). + + 2) Suppose the reverse: that a lazy PLT entry is + being updated, a context switch occurs, and the new + code calls the lazy PLT entry that is being updated. + Then the half-fixed PLT entry will be executed, which + will also almost always cause a segfault. + + These problems don't happen with the 2-word entries, because + only one of the two instructions are changed when a lazy + entry is retargeted at the actual PLT entry; the li + instruction stays the same (we have to update it anyway, + because we might not be updating a lazy PLT entry). */ reloc_addr[0] = OPCODE_LI (11, finaladdr); reloc_addr[1] = OPCODE_ADDIS (11, 11, finaladdr + 0x8000 >> 16); reloc_addr[2] = OPCODE_MTCTR (11); @@ -648,19 +671,27 @@ elf_machine_rela (struct link_map *map, const Elf32_Rela *reloc, num_plt_entries = (map->l_info[DT_PLTRELSZ]->d_un.d_val / sizeof(Elf32_Rela)); + plt[index+PLT_DATA_START_WORDS (num_plt_entries)] = finaladdr; reloc_addr[0] = OPCODE_LI (11, index*4); reloc_addr[1] = OPCODE_B (-(4*(index*2 + 1 - PLT_LONGBRANCH_ENTRY_WORDS + PLT_INITIAL_ENTRY_WORDS))); - plt[index+PLT_DATA_START_WORDS (num_plt_entries)] = finaladdr; } } MODIFIED_CODE (reloc_addr); } else - assert (! "unexpected dynamic reloc type"); + { +#ifdef RTLD_BOOTSTRAP + PPC_DIE; /* There is no point calling _dl_sysdep_error, it + almost certainly hasn't been relocated properly. */ +#else + _dl_sysdep_error (_dl_argv[0] ?: "<program name unknown>", + ": Unknown relocation type\n", NULL); +#endif + } if (rinfo == R_PPC_ADDR16_LO || rinfo == R_PPC_ADDR16_HI || diff --git a/sysdeps/powerpc/e_sqrt.c b/sysdeps/powerpc/e_sqrt.c new file mode 100644 index 0000000000..df80973f58 --- /dev/null +++ b/sysdeps/powerpc/e_sqrt.c @@ -0,0 +1,141 @@ +/* Single-precision floating point square root. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <math.h> +#include <math_private.h> +#include <fenv_libc.h> +#include <inttypes.h> + +static const double almost_half = 0.5000000000000001; /* 0.5 + 2^-53 */ +static const uint32_t a_nan = 0x7fc00000; +static const uint32_t a_inf = 0x7f800000; +static const float two108 = 3.245185536584267269e+32; +static const float twom54 = 5.551115123125782702e-17; +extern const float __t_sqrt[1024]; + +/* The method is based on a description in + Computation of elementary functions on the IBM RISC System/6000 processor, + P. W. Markstein, IBM J. Res. Develop, 34(1) 1990. + Basically, it consists of two interleaved Newton-Rhapson approximations, + one to find the actual square root, and one to find its reciprocal + without the expense of a division operation. The tricky bit here + is the use of the POWER/PowerPC multiply-add operation to get the + required accuracy with high speed. + + The argument reduction works by a combination of table lookup to + obtain the initial guesses, and some careful modification of the + generated guesses (which mostly runs on the integer unit, while the + Newton-Rhapson is running on the FPU). */ +double +__sqrt(double x) +{ + const float inf = *(const float *)&a_inf; + /* x = f_wash(x); *//* This ensures only one exception for SNaN. */ + if (x > 0) + { + if (x != inf) + { + /* Variables named starting with 's' exist in the + argument-reduced space, so that 2 > sx >= 0.5, + 1.41... > sg >= 0.70.., 0.70.. >= sy > 0.35... . + Variables named ending with 'i' are integer versions of + floating-point values. */ + double sx; /* The value of which we're trying to find the + square root. */ + double sg,g; /* Guess of the square root of x. */ + double sd,d; /* Difference between the square of the guess and x. */ + double sy; /* Estimate of 1/2g (overestimated by 1ulp). */ + double sy2; /* 2*sy */ + double e; /* Difference between y*g and 1/2 (se = e * fsy). */ + double shx; /* == sx * fsg */ + double fsg; /* sg*fsg == g. */ + fenv_t fe; /* Saved floating-point environment (stores rounding + mode and whether the inexact exception is + enabled). */ + uint32_t xi0, xi1, sxi, fsgi; + const float *t_sqrt; + + fe = fegetenv_register(); + EXTRACT_WORDS (xi0,xi1,x); + relax_fenv_state(); + sxi = xi0 & 0x3fffffff | 0x3fe00000; + INSERT_WORDS (sx, sxi, xi1); + t_sqrt = __t_sqrt + (xi0 >> 52-32-8-1 & 0x3fe); + sg = t_sqrt[0]; + sy = t_sqrt[1]; + + /* Here we have three Newton-Rhapson iterations each of a + division and a square root and the remainder of the + argument reduction, all interleaved. */ + sd = -(sg*sg - sx); + fsgi = xi0 + 0x40000000 >> 1 & 0x7ff00000; + sy2 = sy + sy; + sg = sy*sd + sg; /* 16-bit approximation to sqrt(sx). */ + INSERT_WORDS (fsg, fsgi, 0); + e = -(sy*sg - almost_half); + sd = -(sg*sg - sx); + if ((xi0 & 0x7ff00000) == 0) + goto denorm; + sy = sy + e*sy2; + sg = sg + sy*sd; /* 32-bit approximation to sqrt(sx). */ + sy2 = sy + sy; + e = -(sy*sg - almost_half); + sd = -(sg*sg - sx); + sy = sy + e*sy2; + shx = sx * fsg; + sg = sg + sy*sd; /* 64-bit approximation to sqrt(sx), + but perhaps rounded incorrectly. */ + sy2 = sy + sy; + g = sg * fsg; + e = -(sy*sg - almost_half); + d = -(g*sg - shx); + sy = sy + e*sy2; + fesetenv_register (fe); + return g + sy*d; + denorm: + /* For denormalised numbers, we normalise, calculate the + square root, and return an adjusted result. */ + fesetenv_register (fe); + return __sqrt(x * two108) * twom54; + } + } + else if (x < 0) + { +#ifdef FE_INVALID_SQRT + feraiseexcept (FE_INVALID_SQRT); + /* For some reason, some PowerPC processors don't implement + FE_INVALID_SQRT. I guess no-one ever thought they'd be + used for square roots... :-) */ + if (!fetestexcept (FE_INVALID)) +#endif + feraiseexcept (FE_INVALID); +#ifndef _IEEE_LIBM + if (_LIB_VERSION != _IEEE_) + x = __kernel_standard(x,x,26); + else +#endif + x = *(const float*)&a_nan; + } + return f_wash(x); +} + +weak_alias (__sqrt, sqrt) +/* Strictly, this is wrong, but the only places where _ieee754_sqrt is + used will not pass in a negative result. */ +strong_alias(__sqrt,__ieee754_sqrt) diff --git a/sysdeps/powerpc/e_sqrtf.c b/sysdeps/powerpc/e_sqrtf.c new file mode 100644 index 0000000000..804dff3c44 --- /dev/null +++ b/sysdeps/powerpc/e_sqrtf.c @@ -0,0 +1,136 @@ +/* Single-precision floating point square root. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <math.h> +#include <math_private.h> +#include <fenv_libc.h> +#include <inttypes.h> + +static const float almost_half = 0.50000006; /* 0.5 + 2^-24 */ +static const uint32_t a_nan = 0x7fc00000; +static const uint32_t a_inf = 0x7f800000; +static const float two48 = 281474976710656.0; +static const float twom24 = 5.9604644775390625e-8; +extern const float __t_sqrt[1024]; + +/* The method is based on a description in + Computation of elementary functions on the IBM RISC System/6000 processor, + P. W. Markstein, IBM J. Res. Develop, 34(1) 1990. + Basically, it consists of two interleaved Newton-Rhapson approximations, + one to find the actual square root, and one to find its reciprocal + without the expense of a division operation. The tricky bit here + is the use of the POWER/PowerPC multiply-add operation to get the + required accuracy with high speed. + + The argument reduction works by a combination of table lookup to + obtain the initial guesses, and some careful modification of the + generated guesses (which mostly runs on the integer unit, while the + Newton-Rhapson is running on the FPU). */ +float +__sqrtf(float x) +{ + const float inf = *(const float *)&a_inf; + /* x = f_washf(x); *//* This ensures only one exception for SNaN. */ + if (x > 0) + { + if (x != inf) + { + /* Variables named starting with 's' exist in the + argument-reduced space, so that 2 > sx >= 0.5, + 1.41... > sg >= 0.70.., 0.70.. >= sy > 0.35... . + Variables named ending with 'i' are integer versions of + floating-point values. */ + float sx; /* The value of which we're trying to find the square + root. */ + float sg,g; /* Guess of the square root of x. */ + float sd,d; /* Difference between the square of the guess and x. */ + float sy; /* Estimate of 1/2g (overestimated by 1ulp). */ + float sy2; /* 2*sy */ + float e; /* Difference between y*g and 1/2 (note that e==se). */ + float shx; /* == sx * fsg */ + float fsg; /* sg*fsg == g. */ + fenv_t fe; /* Saved floating-point environment (stores rounding + mode and whether the inexact exception is + enabled). */ + uint32_t xi, sxi, fsgi; + const float *t_sqrt; + + GET_FLOAT_WORD (xi, x); + fe = fegetenv_register (); + relax_fenv_state (); + sxi = xi & 0x3fffffff | 0x3f000000; + SET_FLOAT_WORD (sx, sxi); + t_sqrt = __t_sqrt + (xi >> 23-8-1 & 0x3fe); + sg = t_sqrt[0]; + sy = t_sqrt[1]; + + /* Here we have three Newton-Rhapson iterations each of a + division and a square root and the remainder of the + argument reduction, all interleaved. */ + sd = -(sg*sg - sx); + fsgi = xi + 0x40000000 >> 1 & 0x7f800000; + sy2 = sy + sy; + sg = sy*sd + sg; /* 16-bit approximation to sqrt(sx). */ + e = -(sy*sg - almost_half); + SET_FLOAT_WORD (fsg, fsgi); + sd = -(sg*sg - sx); + sy = sy + e*sy2; + if ((xi & 0x7f800000) == 0) + goto denorm; + shx = sx * fsg; + sg = sg + sy*sd; /* 32-bit approximation to sqrt(sx), + but perhaps rounded incorrectly. */ + sy2 = sy + sy; + g = sg * fsg; + e = -(sy*sg - almost_half); + d = -(g*sg - shx); + sy = sy + e*sy2; + fesetenv_register (fe); + return g + sy*d; + denorm: + /* For denormalised numbers, we normalise, calculate the + square root, and return an adjusted result. */ + fesetenv_register (fe); + return __sqrtf(x * two48) * twom24; + } + } + else if (x < 0) + { +#ifdef FE_INVALID_SQRT + feraiseexcept (FE_INVALID_SQRT); + /* For some reason, some PowerPC processors don't implement + FE_INVALID_SQRT. I guess no-one ever thought they'd be + used for square roots... :-) */ + if (!fetestexcept (FE_INVALID)) +#endif + feraiseexcept (FE_INVALID); +#ifndef _IEEE_LIBM + if (_LIB_VERSION != _IEEE_) + x = __kernel_standard(x,x,126); + else +#endif + x = *(const float*)&a_nan; + } + return f_washf(x); +} + +weak_alias (__sqrtf, sqrtf) +/* Strictly, this is wrong, but the only places where _ieee754_sqrt is + used will not pass in a negative result. */ +strong_alias(__sqrtf,__ieee754_sqrtf) diff --git a/sysdeps/powerpc/fclrexcpt.c b/sysdeps/powerpc/fclrexcpt.c index 1e66140c2e..cfcf175637 100644 --- a/sysdeps/powerpc/fclrexcpt.c +++ b/sysdeps/powerpc/fclrexcpt.c @@ -19,6 +19,7 @@ #include <fenv_libc.h> +#undef feclearexcept void feclearexcept (int excepts) { @@ -28,7 +29,8 @@ feclearexcept (int excepts) u.fenv = fegetenv_register (); /* Clear the relevant bits. */ - u.l[1] = u.l[1] & ~FE_to_sticky (excepts); + u.l[1] = u.l[1] & ~(-((excepts) >> 31-FPSCR_VX & 1) & FE_ALL_INVALID + | (excepts) & FPSCR_STICKY_BITS); /* Put the new state in effect. */ fesetenv_register (u.fenv); diff --git a/sysdeps/powerpc/fe_nomask.c b/sysdeps/powerpc/fe_nomask.c new file mode 100644 index 0000000000..bca49055d1 --- /dev/null +++ b/sysdeps/powerpc/fe_nomask.c @@ -0,0 +1,32 @@ +/* Procedure definition for FE_NOMASK_ENV. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <fenv.h> +#include <errno.h> + +/* This is presently a stub, until it's decided how the kernels should + support this. */ + +const fenv_t * +__fe_nomask_env(void) +{ + __set_errno (ENOSYS); + return FE_ENABLED_ENV; +} +stub_warning (__fe_nomask_env) diff --git a/sysdeps/powerpc/fegetround.c b/sysdeps/powerpc/fegetround.c index 05395f0797..3bb9fb49e0 100644 --- a/sysdeps/powerpc/fegetround.c +++ b/sysdeps/powerpc/fegetround.c @@ -19,13 +19,11 @@ #include <fenv_libc.h> +#undef fegetround int fegetround (void) { - fenv_union_t u; - - u.fenv = fegetenv_register (); - - /* The rounding mode is bits 30 and 31 of the FPSCR. */ - return u.l[1] & 3; + int result; + asm ("mcrfs 7,7 ; mfcr %0" : "=r"(result) : : "cr7"); \ + return result & 3; } diff --git a/sysdeps/powerpc/fenv_const.c b/sysdeps/powerpc/fenv_const.c index fa35fbc0cf..506bd43df0 100644 --- a/sysdeps/powerpc/fenv_const.c +++ b/sysdeps/powerpc/fenv_const.c @@ -25,5 +25,9 @@ const unsigned long long __fe_dfl_env __attribute__ ((aligned (8))) = 0xfff8000000000000ULL; /* Floating-point environment where none of the exceptions are masked. */ -const unsigned long long __fe_nomask_env __attribute__ ((aligned (8))) = +const unsigned long long __fe_enabled_env __attribute__ ((aligned (8))) = 0xfff80000000000f8ULL; + +/* Floating-point environment with the NI bit set. */ +const unsigned long long __fe_nonieee_env __attribute__ ((aligned (8))) = +0xfff8000000000004ULL; diff --git a/sysdeps/powerpc/fenv_libc.h b/sysdeps/powerpc/fenv_libc.h index 45d61e1565..343be16fcb 100644 --- a/sysdeps/powerpc/fenv_libc.h +++ b/sysdeps/powerpc/fenv_libc.h @@ -22,23 +22,17 @@ #include <fenv.h> -/* Transform a logical or of the FE_* bits to a bit pattern for the - appropriate sticky bits in the FPSCR. */ -#define FE_to_sticky(excepts) \ - (-(excepts & FE_INVALID) & FE_ALL_INVALID \ - | (excepts) & (FE_ALL_EXCEPT & ~FE_INVALID | FE_ALL_INVALID)) - /* The sticky bits in the FPSCR indicating exceptions have occurred. */ #define FPSCR_STICKY_BITS ((FE_ALL_EXCEPT | FE_ALL_INVALID) & ~FE_INVALID) /* Equivalent to fegetenv, but returns a fenv_t instead of taking a pointer. */ #define fegetenv_register() \ - ({ fenv_t env; asm ("mffs %0" : "=f" (env)); env; }) + ({ fenv_t env; asm volatile ("mffs %0" : "=f" (env)); env; }) /* Equivalent to fesetenv, but takes a fenv_t instead of a pointer. */ #define fesetenv_register(env) \ - ({ double d = (env); asm ("mtfsf 0xff,%0" : : "f" (d)); }) + ({ double d = (env); asm volatile ("mtfsf 0xff,%0" : : "f" (d)); }) /* This very handy macro: - Sets the rounding mode to 'round to nearest'; @@ -48,10 +42,65 @@ functions. */ #define relax_fenv_state() asm ("mtfsfi 7,0") +/* Set/clear a particular FPSCR bit (for instance, + reset_fpscr_bit(FPSCR_VE); + prevents INVALID exceptions from being raised). */ +#define set_fpscr_bit(x) asm volatile ("mtfsb1 %0" : : "i"(x)) +#define reset_fpscr_bit(x) asm volatile ("mtfsb0 %0" : : "i"(x)) + typedef union { fenv_t fenv; unsigned int l[2]; } fenv_union_t; +/* Definitions of all the FPSCR bit numbers */ +enum { + FPSCR_FX = 0, /* exception summary */ + FPSCR_FEX, /* enabled exception summary */ + FPSCR_VX, /* invalid operation summary */ + FPSCR_OX, /* overflow */ + FPSCR_UX, /* underflow */ + FPSCR_ZX, /* zero divide */ + FPSCR_XX, /* inexact */ + FPSCR_VXSNAN, /* invalid operation for SNaN */ + FPSCR_VXISI, /* invalid operation for Inf-Inf */ + FPSCR_VXIDI, /* invalid operation for Inf/Inf */ + FPSCR_VXZDZ, /* invalid operation for 0/0 */ + FPSCR_VXIMZ, /* invalid operation for Inf*0 */ + FPSCR_VXVC, /* invalid operation for invalid compare */ + FPSCR_FR, /* fraction rounded [fraction was incremented by round] */ + FPSCR_FI, /* fraction inexact */ + FPSCR_FPRF_C, /* result class descriptor */ + FPSCR_FPRF_FL, /* result less than (usually, less than 0) */ + FPSCR_FPRF_FG, /* result greater than */ + FPSCR_FPRF_FE, /* result equal to */ + FPSCR_FPRF_FU, /* result unordered */ + FPSCR_20, /* reserved */ + FPSCR_VXSOFT, /* invalid operation set by software */ + FPSCR_VXSQRT, /* invalid operation for square root */ + FPSCR_VXCVI, /* invalid operation for invalid integer convert */ + FPSCR_VE, /* invalid operation exception enable */ + FPSCR_OE, /* overflow exception enable */ + FPSCR_UE, /* underflow exception enable */ + FPSCR_ZE, /* zero divide exception enable */ + FPSCR_XE, /* inexact exception enable */ + FPSCR_NI /* non-IEEE mode (typically, no denormalised numbers) */ + /* the remaining two least-significant bits keep the rounding mode */ +}; + +/* This operation (i) sets the appropriate FPSCR bits for its + parameter, (ii) converts SNaN to the corresponding NaN, and (iii) + otherwise passes its parameter through unchanged (in particular, -0 + and +0 stay as they were). The `obvious' way to do this is optimised + out by gcc. */ +#define f_wash(x) \ + ({ double d; asm volatile ("fmul %0,%1,%2" \ + : "=f"(d) \ + : "f" (x), "f"((float)1.0)); d; }) +#define f_washf(x) \ + ({ float f; asm volatile ("fmuls %0,%1,%2" \ + : "=f"(f) \ + : "f" (x), "f"((float)1.0)); f; }) + #endif /* fenv_libc.h */ diff --git a/sysdeps/powerpc/fgetexcptflg.c b/sysdeps/powerpc/fgetexcptflg.c index d6bd58dc69..adbdfe22b9 100644 --- a/sysdeps/powerpc/fgetexcptflg.c +++ b/sysdeps/powerpc/fgetexcptflg.c @@ -23,11 +23,10 @@ void fegetexceptflag (fexcept_t *flagp, int excepts) { fenv_union_t u; - unsigned int flag; /* Get the current state. */ u.fenv = fegetenv_register (); - /* Return that portion that corresponds to the requested exceptions. */ - *flagp = flag = u.l[1] & FPSCR_STICKY_BITS & FE_to_sticky (excepts); + /* Return (all of) it. */ + *flagp = u.l[1]; } diff --git a/sysdeps/powerpc/fpu_control.h b/sysdeps/powerpc/fpu_control.h index d21987ffd8..c4c478672b 100644 --- a/sysdeps/powerpc/fpu_control.h +++ b/sysdeps/powerpc/fpu_control.h @@ -52,6 +52,7 @@ typedef unsigned int fpu_control_t __attribute__ ((__mode__ (__SI__))); #define _FPU_GETCW(cw) ( { \ fpu_control_t tmp[2] __attribute__ ((__aligned__(8))); \ __asm__ ("mffs 0; stfd 0,%0" : "=m" (*tmp) : : "fr0"); \ + (cw)=tmp[1]; \ tmp[1]; } ) #define _FPU_SETCW(cw) { \ fpu_control_t tmp[2] __attribute__ ((__aligned__(8))); \ diff --git a/sysdeps/powerpc/fraiseexcpt.c b/sysdeps/powerpc/fraiseexcpt.c index 4305c3dcd2..d0c7971836 100644 --- a/sysdeps/powerpc/fraiseexcpt.c +++ b/sysdeps/powerpc/fraiseexcpt.c @@ -19,6 +19,7 @@ #include <fenv_libc.h> +#undef feraiseexcept void feraiseexcept (int excepts) { @@ -36,9 +37,17 @@ feraiseexcept (int excepts) u.l[1] = (u.l[1] | excepts & FPSCR_STICKY_BITS /* Turn FE_INVALID into FE_INVALID_SOFTWARE. */ - | excepts << (31 - 21) - (31 - 24) & FE_INVALID_SOFTWARE); + | (excepts >> (31 - FPSCR_VX) - (31 - FPSCR_VXSOFT) + & FE_INVALID_SOFTWARE)); /* Store the new status word (along with the rest of the environment), triggering any appropriate exceptions. */ fesetenv_register (u.fenv); + + if ((excepts & FE_INVALID) + /* For some reason, some PowerPC chips (the 601, in particular) + don't have FE_INVALID_SOFTWARE implemented. Detect this + case and raise FE_INVALID_SNAN instead. */ + && !fetestexcept (FE_INVALID)) + set_fpscr_bit (FPSCR_VXSNAN); } diff --git a/sysdeps/powerpc/fsetexcptflg.c b/sysdeps/powerpc/fsetexcptflg.c index 4279b74849..b762552dbf 100644 --- a/sysdeps/powerpc/fsetexcptflg.c +++ b/sysdeps/powerpc/fsetexcptflg.c @@ -23,15 +23,26 @@ void fesetexceptflag (const fexcept_t *flagp, int excepts) { fenv_union_t u; + fexcept_t flag; /* Get the current state. */ u.fenv = fegetenv_register (); + /* Ignore exceptions not listed in 'excepts'. */ + flag = *flagp & excepts; + /* Replace the exception status */ - u.l[1] = u.l[1] & FPSCR_STICKY_BITS | *flagp & FE_to_sticky (excepts); + u.l[1] = (u.l[1] & ~(FPSCR_STICKY_BITS & excepts) + | flag & FPSCR_STICKY_BITS + | (flag >> (31 - FPSCR_VX) - (31 - FPSCR_VXSOFT) + & FE_INVALID_SOFTWARE)); /* Store the new status word (along with the rest of the environment). This may cause floating-point exceptions if the restored state requests it. */ fesetenv_register (u.fenv); + + /* Deal with FE_INVALID_SOFTWARE not being implemented on some chips. */ + if (flag & FE_INVALID) + feraiseexcept(FE_INVALID); } diff --git a/sysdeps/powerpc/ftestexcept.c b/sysdeps/powerpc/ftestexcept.c index 52733f7ae7..6cb8dd5e28 100644 --- a/sysdeps/powerpc/ftestexcept.c +++ b/sysdeps/powerpc/ftestexcept.c @@ -23,16 +23,11 @@ int fetestexcept (int excepts) { fenv_union_t u; - int flags; /* Get the current state. */ u.fenv = fegetenv_register (); - /* Find the bits that indicate exceptions have occurred. */ - flags = u.l[1] & FPSCR_STICKY_BITS; - - /* Set the FE_INVALID bit if any of the FE_INVALID_* bits are set. */ - flags |= ((u.l[1] & FE_ALL_INVALID) != 0) << 31-24; - - return flags & excepts; + /* The FE_INVALID bit is dealt with correctly by the hardware, so we can + just: */ + return u.l[1] & excepts; } diff --git a/sysdeps/powerpc/q_addsub.c b/sysdeps/powerpc/q_addsub.c new file mode 100644 index 0000000000..e4ef6d8165 --- /dev/null +++ b/sysdeps/powerpc/q_addsub.c @@ -0,0 +1,549 @@ +/* Add or subtract two 128-bit floating point values. C prototype. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <quad_float.h> + +/* Add 'a' to 'b' and put the result in 'result', but treat a[0]=axx, + b[0]=bxx. bxx differs from b[0] only in the high bit, similarly axx. */ +/* Exceptions to raise: + - Invalid (SNaN) + - Invalid (Inf-Inf) + - Overflow + - Underflow + - Inexact + */ + +/* Handle cases where exponent of a or b is maximum. */ +static void +handle_max_exponent(unsigned result[4], + const unsigned a[4], const unsigned b[4], + const unsigned axx, /* Treat as a[0]. */ + const unsigned bxx, /* Treat as b[0]. */ + const unsigned ax, /* axx >> 16 & 0x7fff. */ + const unsigned bx) /* bxx >> 16 & 0x7fff. */ +{ + int ax_ismax, bx_ismax; + unsigned a1,a2,a3, b1,b2,b3; + int a_zeromant, b_zeromant; + + ax_ismax = ax == 0x7fff; + bx_ismax = bx == 0x7fff; + + assert(ax_ismax || bx_ismax); + + a1 = a[1]; a2 = a[2]; a3 = a[3]; + b1 = b[1]; b2 = b[2]; b3 = b[3]; + + a_zeromant = (axx & 0xffff | a1 | a2 | a3) == 0; + b_zeromant = (bxx & 0xffff | b1 | b2 | b3) == 0; + + /* Deal with SNaNs. */ + if ( ax_ismax && !a_zeromant && (axx & 0x8000) == 0 + || bx_ismax && !b_zeromant && (bxx & 0x8000) == 0) + { + set_fpscr_bit(FPSCR_VXSNAN); + axx |= 0x8000; /* Demote the SNaN to a QNaN (whichever of */ + bxx |= 0x8000; /* a or b it was). */ + } + /* Deal with Inf-Inf. */ + else if (a_zeromant && b_zeromant && (axx ^ bxx) == 0x80000000) + { + set_fpscr_bit(FPSCR_VXISI); + bxx |= 0x8000; /* Return an appropriate QNaN. */ + } + + /* Return the lexicographically larger of a or b, ignoring the sign + bits. */ + if ((axx & 0x7fffffff) > (bxx & 0x7fffffff)) goto return_a; + else if ((axx & 0x7fffffff) < (bxx & 0x7fffffff)) goto return_b; + else if (a1 > b1) goto return_a; + else if (a1 < b1) goto return_b; + else if (a2 > b2) goto return_a; + else if (a2 < b2) goto return_b; + else if (a3 > b3) goto return_a; /* I've clearly been writing too */ + else if (a3 < b3) goto return_b; /* much Fortran... */ + + /* If they are equal except for the sign bits, return 'b'. */ + +return_b: + result[0] = bxx; result[1] = b1; result[2] = b2; result[3] = b3; + return; + +return_a: + result[0] = axx; result[1] = a1; result[2] = a2; result[3] = a3; + return; +} + +/* Renormalise and output a FP number. */ +static void +renormalise_value(unsigned result[4], + const unsigned axx, + unsigned ax, + unsigned r0, + unsigned r1, + unsigned r2, + unsigned r3) +{ + int rshift; + if (r0 != 0 || cntlzw(a1) < 16 || 32 > ax-1) + { + rshift = cntlzw(r0)-15 + (-(cntlzw(r0) >> 5) & cntlzw(a1)); + assert(rshift < 32); + if (rshift > ax-1) + { + ax--; + rshift = ax; + } + + result[0] = (axx & 0x80000000 + | ax-rshift << 16 + | r0 << rshift & 0xffff + | a1 >> 32-rshift & 0xffff); + result[1] = a1 << rshift | a2 >> 32-rshift; + result[2] = a2 << rshift | a3 >> 32-rshift; + result[3] = a3 << rshift; + return; + } + result[3] = 0; + /* Special case for zero. */ + if (a1 == 0 && a2 == 0 && a3 == 0) + { + result[0] = axx & 0x80000000; + result[1] = result[2] = 0; + return; + } + while (a1 != 0 && cntlzw(a2) >= 16 && 64 <= ax-1) + { + ax -= 32; + a1 = a2; a2 = a3; a3 = 0; + } + rshift = cntlzw(a1)-15 + (-(cntlzw(a1) >> 5) & cntlzw(a2)); + assert(rshift < 32); + if (rshift > ax-1-32) + { + ax--; + rshift = ax-32; + } + + result[0] = (axx & 0x80000000 + | ax-rshift-32 << 16 + | a1 << rshift & 0xffff + | a2 >> 32-rshift & 0xffff); + result[1] = a2 << rshift | a3 >> 32-rshift; + result[2] = a3 << rshift; + return; +} + +/* Handle the case where one or both numbers are denormalised or zero. + This case almost never happens, so we don't slow the main code + with it. */ +static void +handle_min_exponent(unsigned result[4], + const unsigned a[4], const unsigned b[4], + const unsigned axx, /* Treat as a[0]. */ + const unsigned bxx, /* Treat as b[0]. */ + const unsigned ax, /* axx >> 16 & 0x7fff. */ + const unsigned bx) /* bxx >> 16 & 0x7fff. */ +{ + int ax_denorm, bx_denorm; + unsigned a1,a2,a3, b1,b2,b3; + int a_zeromant, b_zeromant; + + ax_denorm = ax == 0; + bx_denorm = bx == 0; + + assert(ax_denorm || bx_denorm); + + a1 = a[1]; a2 = a[2]; a3 = a[3]; + b1 = b[1]; b2 = b[2]; b3 = b[3]; + + +} + +/* Add a+b+cin modulo 2^32, put result in 'r' and carry in 'cout'. */ +#define addc(r,cout,a,b,cin) \ + do { \ + unsigned long long addc_tmp = (a)+(b)+(cin); + (cout) = addc_tmp >> 32; + (r) = addc_tmp; + } + +/* Calculate a+~b+cin modulo 2^32, put result in 'r' and carry in 'cout'. */ +#define subc(r,cout,a,b,cin) \ + do { \ + unsigned long long addc_tmp = (a)-(b)+(cin)-1; + (cout) = addc_tmp >> 63; + (r) = addc_tmp; + } + +/* Handle the case where both exponents are the same. This requires quite + a different algorithm than the general case. */ +static void +handle_equal_exponents(unsigned result[4], + const unsigned a[4], const unsigned b[4], + const unsigned axx, /* Treat as a[0]. */ + const unsigned bxx, /* Treat as b[0]. */ + unsigned ax) /* [ab]xx >> 16 & 0x7fff. */ +{ + unsigned a1,a2,a3, b1,b2,b3; + int roundmode; + unsigned carry, r0; + + a1 = a[1]; a2 = a[2]; a3 = a[3]; + b1 = b[1]; b2 = b[2]; b3 = b[3]; + + if ((int)(axx ^ bxx) >= 0) + { + int roundmode; + + /* Adding. */ + roundmode = fegetround(); + + /* What about overflow? */ + if (ax == 0x7ffe) + { + /* Oh no! Too big! */ + /* Result: + rounding result + -------- ------ + nearest return Inf with sign of a,b + zero return nearest possible non-Inf value with + sign of a,b + +Inf return +Inf if a,b>0, otherwise return + value just before -Inf. + -Inf return +Inf if a,b>0, otherwise return + value just before -Inf. + */ + set_fpscr_bit(FPSCR_OX); + /* Overflow always produces inexact result. */ + set_fpscr_bit(FPSCR_XX); + + if ( roundmode == FE_TONEAREST + || roundmode == FE_UPWARD && (int)axx >= 0 + || roundmode == FE_DOWNWARD && (int)axx < 0) + { + result[3] = result[2] = result[1] = 0; + result[0] = axx & 0xffff0000 | 0x7fff0000; + } + else + { + result[3] = result[2] = result[1] = 0xffffffff; + result[0] = axx & 0xfffe0000 | 0x7ffeffff; + } + return; + } + + /* We need to worry about rounding/inexact here. Do it like this: */ + if (a3 + b3 & 1) + { + /* Need to round. Upwards? */ + set_fpscr_bit(FPSCR_XX); + carry = ( roundmode == FE_NEAREST && (a3 + b3 & 2) != 0 + || roundmode == FE_UPWARD && (int)axx >= 0 + || roundmode == FE_DOWNWARD && (int)axx < 0); + } + else + carry = 0; /* Result will be exact. */ + + /* Perform the addition. */ + addc(a3,carry,a3,b3,carry); + addc(a2,carry,a2,b2,carry); + addc(a1,carry,a1,b1,carry); + r0 = (axx & 0xffff) + (bxx & 0xffff) + carry; + + /* Shift right by 1. */ + result[3] = a3 >> 1 | a2 << 31; + result[2] = a2 >> 1 | a1 << 31; + result[1] = a1 >> 1 | r0 << 31; + /* Exponent of result is exponent of inputs plus 1. + Sign of result is common sign of inputs. */ + result[0] = r0 >> 1 & 0xffff | axx + 0x10000 & 0xffff0000; + } + else + { + /* Subtracting. */ + + /* Perform the subtraction, a-b. */ + subc(a3,carry,a3,b3,0); + subc(a2,carry,a2,b2,carry); + subc(a1,carry,a1,b1,carry); + subc(r0,carry,a0&0xffff,b0&0xffff,carry); + + /* Maybe we should have calculated b-a... */ + if (carry) + { + subc(a3,carry,0,a3,0); + subc(a2,carry,0,a2,carry); + subc(a1,carry,0,a1,carry); + subc(r0,carry,0,r0,carry); + axx ^= 0x80000000; + } + + renormalise_value(result, axx, ax, r0, a1, a2, a3); + } +} + + +static void +add(unsigned result[4], const unsigned a[4], const unsigned b[4], + unsigned axx, unsigned bxx) +{ + int ax, bx, diff, carry; + unsigned a0,a1,a2,a3, b0,b1,b2,b3,b4, sdiff; + + ax = axx >> 16 & 0x7fff; + bx = bxx >> 16 & 0x7fff; + + /* Deal with NaNs and Inf. */ + if (ax == 0x7fff || bx == 0x7fff) + { + handle_max_exponent(result, a, b, axx, bxx, ax, bx); + return; + } + /* Deal with denorms and zero. */ + if (ax == 0 || bx == 0) + { + handle_min_exponent(result, a, b, axx, bxx, ax, bx); + return; + } + /* Finally, one special case, when both exponents are equal. */ + if (ax == bx) + { + handle_equal_exponents(result, a, b, axx, bxx, ax); + return; + } + + sdiff = axx ^ bxx; + /* Swap a and b if b has a larger magnitude than a, so that a will have + the larger magnitude. */ + if (ax < bx) + { + const unsigned *t; + t = b; b = a; a = t; + diff = bx - ax; + ax = bx; + axx = bxx; + } + else + diff = ax - bx; + + a0 = a[0] & 0xffff | 0x10000; a1 = a[1]; a2 = a[2]; a3 = a[3]; + b0 = b[0] & 0xffff | 0x10000; b1 = b[1]; b2 = b[2]; b3 = b[3]; + if (diff < 32) + { + b4 = b3 << 32-diff; + b3 = b3 >> diff | b2 << 32-biff; + b2 = b2 >> diff | b1 << 32-diff; + b1 = b1 >> diff | b0 << 32-diff; + b0 = b0 >> diff; + } + else if (diff < 64) + { + diff -= 32; + b4 = b3 & 1 | b3 >> (diff == 32) | b2 << 32-biff; + b3 = b2 >> diff | b1 << 32-diff; + b2 = b1 >> diff | b0 << 32-diff; + b1 = b0 >> diff; + b0 = 0; + } + else if (diff < 96) + { + b4 = b2 | b3 | b1 << 32-diff; + b3 = b1 >> diff | b0 << 32-diff; + b2 = b0 >> diff; + b1 = b0 = 0; + } + else if (diff < 128) + { + b4 = b1 | b2 | b3 | b0 << 32-diff; + b3 = b0 >> diff; + b2 = b1 = b0 = 0; + } + else + { + b4 = b0|b1|b2|b3; + b3 = b2 = b1 = b0 = 0; + } + + /* Now, two cases: one for addition, one for subtraction. */ + if ((int)sdiff >= 0) + { + /* Addition. */ + + /* + + /* Perform the addition. */ + addc(a3,carry,a3,b3,0); + addc(a2,carry,a2,b2,carry); + addc(a1,carry,a1,b1,carry); + addc(a0,carry,a0,b0,carry); + + + + if (a0 & 0x20000) + { + /* Need to renormalise by shifting right. */ + /* Shift right by 1. */ + b4 = b4 | a3 << 31; + a3 = a3 >> 1 | a2 << 31; + a2 = a2 >> 1 | a1 << 31; + result[1] = a1 >> 1 | r0 << 31; + /* Exponent of result is exponent of inputs plus 1. + Sign of result is common sign of inputs. */ + result[0] = r0 >> 1 & 0xffff | axx + 0x10000 & 0xffff0000; + } + + + } + else + { + /* Subtraction. */ + + } +} + +/* Add the absolute values of two 128-bit floating point values, + give the result the sign of one of them. The only exception this + can raise is for SNaN. */ +static void +aadd(unsigned result[4], const unsigned a[4], const unsigned b[4]) +{ + unsigned ax, bx, xd; + const unsigned *sml; + unsigned t0,t1,t2,t3,tx, s0,s1,s2,s3,s4, carry; + int rmode, xdelta, shift; + + ax = a[0] >> 16 & 0x7fff; + bx = b[0] >> 16 & 0x7fff; + + /* Deal with . */ + if (ax == 0x7fff) + { + t0 = a[0]; t1 = a[1]; t2 = a[2]; t3 = a[3]; + /* Check for SNaN. */ + if ((t0 & 0x8000) == 0 + && (t0 & 0x7fff | t1 | t2 | t3) != 0) + set_fpscr_bit(FPSCR_VXSNAN); + /* Return b. */ + result[0] = t0; result[1] = t1; result[2] = t2; result[3] = t3; + return; + } + /* Deal with b==Inf or b==NaN. */ + if (bx == 0x7fff) + { + t0 = b[0]; t1 = b[1]; t2 = b[2]; t3 = b[3]; + /* Check for SNaN. */ + if ((t0 & 0x8000) == 0 + && (t0 & 0x7fff | t1 | t2 | t3) != 0) + set_fpscr_bit(FPSCR_VXSNAN); + /* Return b. */ + result[0] = t0; result[1] = t1; result[2] = t2; result[3] = t3; + return; + } + + /* Choose the larger of the two to be 't', and the smaller to be 's'. */ + if (ax > bx) + { + t0 = a[0] & 0xffff | (ax != 0) << 16; + t1 = a[1]; t2 = a[2]; t3 = a[3]; tx = ax; + s0 = b[0] & 0xffff | (bx != 0) << 16; + s1 = b[1]; s2 = b[2]; s3 = b[3]; + xd = ax-bx; + } + else + { + t0 = b[0] & 0xffff | (bx != 0) << 16; + t1 = b[1]; t2 = b[2]; t3 = b[3]; tx = bx; + s0 = a[0] & 0xffff | (ax != 0) << 16; + s1 = a[1]; s2 = a[2]; s3 = a[3]; + sml = a; + xd = bx-ax; + } + + /* Shift 's2' right by 'xd' bits. */ + switch (xd >> 5) + { + case 0: + s4 = 0; + break; + case 1: + s4 = s3; s3 = s2; s2 = s1; s1 = s0; s0 = 0; + break; + case 2: + s4 = s2 | s3 != 0; + s3 = s1; s2 = s0; s1 = s0 = 0; + break; + case 3: + s4 = s1 | (s3|s2) != 0; + s3 = s0; s2 = s1 = s0 = 0; + break; + default: + s4 = s0 | (s3|s2|s1) != 0; + s3 = s2 = s1 = s0 = 0; + } + xd = xd & 0x1f; + if (xd != 0) + { + s4 = s4 >> xd | (s4 << 32-xd) != 0 | s3 << 32-xd; + s3 = s3 >> xd | s2 << 32-xd; + s2 = s2 >> xd | s1 << 32-xd; + s1 = s1 >> xd | s0 << 32-xd; + s0 = s0 >> xd; + } + + /* Do the addition. */ +#define addc(r,cout,a,b,cin) \ + do { \ + unsigned long long addc_tmp = (a)+(b)+(cin); + (cout) = addc_tmp >> 32; + (r) = addc_tmp; + } + addc(t3,carry,t3,s3,0); + addc(t2,carry,t2,s2,carry); + addc(t1,carry,t1,s1,carry); + t0 = t0 + s0 + carry; + + /* Renormalise. */ + xdelta = 15-cntlzw(t0); + if (tx + xdelta <= 0x7fff) + shift = xdelta; + else + { + } +} + +/* Add two 128-bit floating point values. */ +void +__q_add(unsigned result[4], const unsigned a[4], const unsigned b[4]) +{ + if ((a[0] ^ b[0]) >= 0) + aadd(result, a, b); + else + asubtract(result, a, b); +} + +/* Subtract two 128-bit floating point values. */ +void +__q_sub(unsigned result[4], const unsigned a[4], const unsigned b[4]) +{ + if ((a[0] ^ b[0]) < 0) + aadd(result, a, b); + else + asubtract(result, a, b); +} diff --git a/sysdeps/powerpc/q_dtoq.c b/sysdeps/powerpc/q_dtoq.c new file mode 100644 index 0000000000..8cab9e4689 --- /dev/null +++ b/sysdeps/powerpc/q_dtoq.c @@ -0,0 +1,66 @@ +/* 64-bit floating point to 128-bit floating point. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <quad_float.h> + +/* long double _q_dtoq(double a); + Convert 'a' to long double. Don't raise exceptions. + */ + +void +__q_dtoq(unsigned long long result[2], double a) +{ + unsigned ux, rx; + union { + double d; + unsigned long long u; + } u; + + u.d = a; + result[1] = u.u << 64-4; + result[0] = u.u >> 4; + /* correct exponent bias */ + rx = ((long long)u.u >> 52 & 0x87ff) + 16383-1023; + + ux = u.u >> 52 & 0x7ff; + if (ux == 0) + { + if ((u.u & 0x7fffffffffffffffULL) == 0) + { + /* +0.0 or -0.0. */ + rx &= 0x8000; + } + else + { + /* Denormalised number. Renormalise. */ + unsigned long long um = u.u & 0x000fffffffffffffULL; + int cs = cntlzd(um) - 12 + 1; + rx -= cs; + um <<= cs; + result[0] = um >> 4; + result[1] = um << 64-4; + } + } + else if (ux == 0x7ff) + { + /* Inf or NaN. */ + rx |= 0x7fff; + } + *(unsigned short *)result = rx; +} diff --git a/sysdeps/powerpc/q_feq.c b/sysdeps/powerpc/q_feq.c new file mode 100644 index 0000000000..3f5512a4a0 --- /dev/null +++ b/sysdeps/powerpc/q_feq.c @@ -0,0 +1,65 @@ +/* 128-bit floating point unordered equality. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <fenv_libc.h> + +/* int _q_feq(const long double *a, const long double *b); + Returns nonzero if a==b, 0 otherwise. + + Special cases: + NaNs are never equal to anything; + -0 == +0; + If either argument is a SNaN, we set FPSCR_VXSNAN to cause an + 'invalid operation' exception. + */ + +int +__q_feq(const unsigned a[4], const unsigned b[4]) +{ + unsigned a0, b0; + + a0 = a[0]; b0 = b[0]; + + if ((a0 >> 16 & 0x7fff) != 0x7fff && + (b0 >> 16 & 0x7fff) != 0x7fff) + { + unsigned a3,b3,a2,b2,a1,b1; + unsigned z, result; + + a3 = a[3]; b3 = b[3]; a2 = a[2]; + b2 = b[2]; a1 = a[1]; b1 = b[1]; + result = ~(a3 ^ b3); + result &= ~(a2 ^ b2); + result &= ~(a1 ^ b1); + z = a3 | a2; + z |= a1; + z = (z >> 1 | z | a0) & 0x7fffffff; + /* 'z' is 0 iff a==0.0, otherwise between 1 and 0x7fffffff */ + return (result & ~(a0^b0 & ~(-z & 0x80000000)))+1; + } + else + { + /* Deal with SNaNs */ + if (((a0|b0) & 0x8000) == 0) + { + set_fpscr_bit(FPSCR_VXSNAN); + } + return 0; + } +} diff --git a/sysdeps/powerpc/q_fne.c b/sysdeps/powerpc/q_fne.c new file mode 100644 index 0000000000..8f6429cd55 --- /dev/null +++ b/sysdeps/powerpc/q_fne.c @@ -0,0 +1,70 @@ +/* 128-bit floating point unordered not-equality. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <fenv_libc.h> + +/* int _q_fne(const long double *a, const long double *b); + Returns 0 if a==b + + Special cases: + NaNs return 1 always; + -0 == +0; + If either argument is a SNaN, we set FPSCR_VXSNAN to cause an + 'invalid operation' exception. + */ + +int +__q_fne(const unsigned a[4], const unsigned b[4]) +{ + unsigned a0, b0; + + a0 = a[0]; b0 = b[0]; + + if ((a0 >> 16 & 0x7fff) != 0x7fff && + (b0 >> 16 & 0x7fff) != 0x7fff) + { + unsigned a3,b3,a2,b2,a1,b1; + unsigned z, result; + + a3 = a[3]; b3 = b[3]; a2 = a[2]; + result = a3 ^ b3; + b2 = b[2]; + if (result != 0) + return result; + a1 = a[1]; b1 = b[1]; + result = a2 ^ b2; + result |= a1 ^ b1; + z = a3 | a2; + if (result != 0) + return result; + z |= a1; + z = (z >> 1 | z | a0) & 0x7fffffff; + /* 'z' is 0 iff a==0.0, otherwise between 1 and 0x7fffffff */ + return (a0^b0) & ~(-z & 0x80000000); + } + else + { + /* Deal with SNaNs */ + if (((a0|b0) & 0x8000) == 0) + { + set_fpscr_bit(FPSCR_VXSNAN); + } + return 1; + } +} diff --git a/sysdeps/powerpc/q_itoq.c b/sysdeps/powerpc/q_itoq.c new file mode 100644 index 0000000000..aa737391bb --- /dev/null +++ b/sysdeps/powerpc/q_itoq.c @@ -0,0 +1,50 @@ +/* 32-bit fixed point signed integer to 128-bit floating point. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <quad_float.h> + +/* long double _q_itoq(int a); + Convert 'a' to long double. + */ + +void +__q_itoq(unsigned int result[4], int a) +{ + unsigned rx; + int cs, css; + + /* Note that if a==-2^31, then ua will be 2^31. */ + unsigned int ua = abs(a); + + /* Normalize. */ + cs = cntlzw(ua); + ua <<= cs+1; + + /* Calculate the exponent, in 4 easy instructions. */ + css = 31-cs; + rx = 16383+css << 16 & ~css; + + /* Copy the sign bit from 'a'. */ + asm ("rlwimi %0,%1,0,0,0": "=r"(rx) : "r"(a), "0"(rx)); + + /* Put it all together. */ + result[2] = result[3] = 0; + asm ("rlwimi %0,%1,16,16,31": "=r"(result[0]) : "r"(ua), "0"(rx)); + result[1] = ua << 16; +} diff --git a/sysdeps/powerpc/q_lltoq.c b/sysdeps/powerpc/q_lltoq.c new file mode 100644 index 0000000000..6222eff9fc --- /dev/null +++ b/sysdeps/powerpc/q_lltoq.c @@ -0,0 +1,52 @@ +/* 64-bit fixed point signed integer to 128-bit floating point. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <quad_float.h> + +/* long double _q_ulltoq(long long a); + Convert 'a' to long double. + */ + +void +__q_lltoq(unsigned int result[4], long long a) +{ + unsigned rx; + int cs, css; + + /* Calculate absolute value of 'a'. */ + unsigned long long ua = (a ^ a >> 63) - (a >> 63); + + /* Normalize. */ + cs = cntlzd(a); + ua <<= cs+1; + + /* Calculate the exponent, in 4 easy instructions. */ + css = 63-cs; + rx = 16383+css << 16 & ~css; + + /* Copy the sign bit from 'a'. */ + asm ("rlwimi %0,%1,0,0,0": "=r"(rx) : "r"(a >> 32), "0"(rx)); + + /* Put it all together. */ + result[3] = 0; + asm ("rlwimi %0,%1,16,16,31": "=r"(result[0]) : "r"(ua >> 32), "0"(rx)); + ua <<= 16; + result[2] = ua; + result[1] = ua >> 32; +} diff --git a/sysdeps/powerpc/q_neg.c b/sysdeps/powerpc/q_neg.c new file mode 100644 index 0000000000..8e11afdbc6 --- /dev/null +++ b/sysdeps/powerpc/q_neg.c @@ -0,0 +1,38 @@ +/* Negate a 128-bit floating point value. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <quad_float.h> + +/* long double _q_neg(const long double *a); + Negate 'a'. Don't raise exceptions. + */ + +void +__q_neg(unsigned int result[4], unsigned int a[4]) +{ + unsigned int t1,t2; + t1 = a[0]; + t2 = a[1]; + result[0] = t1 ^ 0x80000000; + t1 = a[2]; + result[1] = t2; + t2 = a[3]; + result[2] = t1; + result[3] = t2; +} diff --git a/sysdeps/powerpc/q_qtoi.c b/sysdeps/powerpc/q_qtoi.c new file mode 100644 index 0000000000..a2b1b3fe56 --- /dev/null +++ b/sysdeps/powerpc/q_qtoi.c @@ -0,0 +1,65 @@ +/* 128-bit floating point to 32-bit fixed point signed integer. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <quad_float.h> + +/* int _q_qtoi(const long double *a); + Convert 'a' to signed int by truncation. + + Special cases: + NaN: raise VXCVI, return 0x80000000 + SNaN: raise VXSNAN, VXCVI, and return 0x80000000 + >= 2^32: raise VXCVI, return 0x7fffffff + <= -2^32: raise VXCVI, return 0x80000000 + */ + +int +__q_qtoi(const unsigned long long a[2]) +{ + int ax, sgn; + unsigned long long a0; + + a0 = a[0]; + /* 'sgn' is -1 if 'a' is negative, 0 if positive. */ + sgn = (long long)a0 >> 63; + ax = (a0 >> 48 & 0x7fff) - 16383; + /* Deal with non-special cases. */ + if (ax <= 30) + /* If we had a '>>' operator that behaved properly with respect to + large shifts, we wouldn't need the '& -(ax >> 31)' term, which + clears 'result' if ax is negative. The '^ sgn) - sgn' part is + equivalent in this case to multiplication by 'sgn', but usually + faster. */ + return (((unsigned)(a0 >> 17 | 0x80000000) >> 31-ax & -(ax >> 31)) + ^ sgn) - sgn; + + /* Check for SNaN, if so raise VXSNAN. */ + if ((a0 >> 32 & 0x7fff8000) == 0x7fff0000 + && (a[1] | a0 & 0x00007fffffffffffULL) != 0) + set_fpscr_bit(FPSCR_VXSNAN); + + /* Treat all NaNs as large negative numbers. */ + sgn |= -cntlzw(~(a0 >> 32) & 0x7fff0000) >> 5; + + /* All the specials raise VXCVI. */ + set_fpscr_bit(FPSCR_VXCVI); + + /* Return 0x7fffffff (if a < 0) or 0x80000000 (otherwise). */ + return sgn ^ 0x7fffffff; +} diff --git a/sysdeps/powerpc/q_qtoll.c b/sysdeps/powerpc/q_qtoll.c new file mode 100644 index 0000000000..77a54f8bf9 --- /dev/null +++ b/sysdeps/powerpc/q_qtoll.c @@ -0,0 +1,64 @@ +/* 128-bit floating point to 64-bit fixed point signed integer. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <quad_float.h> + +/* long long _q_qtoll(const long double *a); + Convert 'a' to signed long long by truncation. + + Special cases: + NaN: raise VXCVI, return 0x8000000000000000 + SNaN: raise VXSNAN, VXCVI, and return 0x8000000000000000 + >= 2^32: raise VXCVI, return 0x7fffffffffffffff + <= -2^32: raise VXCVI, return 0x8000000000000000 + */ + +long long +__q_qtoll(const unsigned long long a[2]) +{ + int ax, sgn; + unsigned long long a0, a1; + + a0 = a[0]; a1 = a[1]; + /* 'sgn' is -1 if 'a' is negative, 0 if positive. */ + sgn = (long long)a0 >> 63; + ax = (a0 >> 48 & 0x7fff) - 16383; + /* Deal with non-special cases. */ + if (ax <= 62) + /* If we had a '>>' operator that behaved properly with respect to + large shifts, we wouldn't need the '& -(ax >> 31)' term, which + clears 'result' if ax is negative. The '^ sgn) - sgn' part is + equivalent in this case to multiplication by 'sgn', but faster. */ + return (((a0 << 15 | a1 >> 64-15 | 0x8000000000000000ULL) >> 63-ax + & -(unsigned long long)(ax >> 31)) + ^ (long long)sgn) - sgn; + + /* Check for SNaN, if so raise VXSNAN. */ + if ((a0 >> 32 & 0x7fff8000) == 0x7fff0000 + && (a1 | a0 & 0x00007fffffffffffULL) != 0) + set_fpscr_bit(FPSCR_VXSNAN); + + /* Treat all NaNs as large negative numbers. */ + sgn |= -cntlzw(~(a0 >> 32) & 0x7fff0000) >> 5; + + /* All the specials raise VXCVI. */ + set_fpscr_bit(FPSCR_VXCVI); + + return (long long)sgn ^ 0x7fffffffffffffffLL; +} diff --git a/sysdeps/powerpc/q_qtos.c b/sysdeps/powerpc/q_qtos.c new file mode 100644 index 0000000000..cd2715f8e9 --- /dev/null +++ b/sysdeps/powerpc/q_qtos.c @@ -0,0 +1,84 @@ +/* 128-bit floating point to 32-bit floating point. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <quad_float.h> + +/* float _q_qtos(const long double *a); + Convert 'a' to float. Round as per current rounding flags. + + Input Rounding Output + +/-0 * +/-0 + +/-Inf * +/-Inf + +/-NaN * +/-NaN (with mantissa truncated) + +/-SNaN * +/-NaN (with mantissa truncated, MSB of mantissa <- 1) + && raise VXSNAN + [Note: just truncating the mantissa may not give you + a SNaN!] + |a|>=2^128 Nearest +/-Inf && raise overflow && raise inexact + |a|>=2^128 Truncate +/-(2^128-2^104) && raise overflow && raise inexact + a>=2^128 +Inf +Inf && raise overflow && raise inexact + a<=-2^128 +Inf -(2^128-2^104) && raise overflow && raise inexact + a>=2^128 -Inf +(2^128-2^104) && raise overflow && raise inexact + a<=-2^128 -Inf -Inf && raise overflow && raise inexact + + We also need to raise 'inexact' if the result will be inexact, which + depends on the current rounding mode. + + To avoid having to deal with all that, we convert to a 'double' + that will round correctly (but is not itself rounded correctly), + and convert that to a float. This makes this procedure much + simpler and much faster. */ + +float +__q_qtos(const unsigned long long a[2]) +{ + unsigned long long a0,d; + union { + double d; + unsigned long long ull; + } u; + + a0 = a[0]; + + /* Truncate the mantissa to 48 bits. */ + d = a0 << 4; + /* Set the low bit in the mantissa if any of the bits we are dropping + were 1. This ensures correct rounding, and also distinguishes + 0 and Inf from denormalised numbers and SNaN (respectively). */ + d |= a[1] != 0; + /* Copy the sign bit. */ + d = d & 0x7fffffffffffffffULL | a0 & 0x8000000000000000ULL; + + /* Now, we need to fix the exponent. If the exponent of a was in + the range +127 to -152, or was +16384 or -16383, it is already + correct in 'd'. Otherwise, we need to ensure that the new + exponent is in the range +1023 to +128, or -153 to -1022, with + the same sign as the exponent of 'a'. We can do this by setting + bits 1-3 (the second through fourth-most significant bit) of 'd' + to 101 if bit 1 of 'a' is 1, or 010 if bit 1 of 'a' is 0. */ + if ((a0 >> 56 & 0x7f) - 0x3f > 1) + { + unsigned t = (a0 >> 32+2 & 2 << 31-1-2)*3 + (2 << 31-2); + d = (d & 0x8fffffffffffffffULL + | (unsigned long long)t<<32 & 0x7000000000000000ULL); + } + + u.ull = d; + return (float)u.d; +} diff --git a/sysdeps/powerpc/q_qtou.c b/sysdeps/powerpc/q_qtou.c new file mode 100644 index 0000000000..a5583cc816 --- /dev/null +++ b/sysdeps/powerpc/q_qtou.c @@ -0,0 +1,64 @@ +/* 128-bit floating point to 32-bit fixed point unsigned integer. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <quad_float.h> + +/* unsigned int _q_qtou(const long double *a); + Convert 'a' to unsigned int by truncation. + + Special cases: + -0: return 0 + NaN: raise VXCVI, return 0 + SNaN: raise VXSNAN, VXCVI, and return 0 + Negative numbers (other than -0, but including -Inf): raise VXCVI, return 0 + >= 2^32: raise VXCVI, return 0xffffffff + */ + +unsigned int +__q_qtou(const unsigned long long a[2]) +{ + int ax; + unsigned result; + unsigned long long a0; + + /* Deal with non-special cases. */ + a0 = a[0]; + ax = (a0 >> 48) - 16383; + if (ax <= 31) + /* If we had a '>>' operator that behaved properly with respect to + large shifts, we wouldn't need the '& -(ax >> 31)' term, which + clears 'result' if ax is negative. */ + return (unsigned)(a0 >> 17 | 0x80000000) >> 31-ax & -(ax >> 31); + + /* 'result' is 1 if a is negative, 0 otherwise. */ + result = a0 >> 63; + /* Check for -0, otherwise raise VXCVI. */ + if (a0 != 0x8000000000000000ULL || a[1] != 0) + { + /* Check for SNaN, if so raise VXSNAN. */ + if ((a0 >> 32 & 0x7fff8000) == 0x7fff0000 + && (a[1] | a0 & 0x00007fffffffffffULL) != 0) + set_fpscr_bit(FPSCR_VXSNAN); + /* Treat all NaNs as large negative numbers. */ + result |= cntlzw(~(a0 >> 32) & 0x7fff0000) >> 5; + + set_fpscr_bit(FPSCR_VXCVI); + } + return result-1; +} diff --git a/sysdeps/powerpc/q_qtoull.c b/sysdeps/powerpc/q_qtoull.c new file mode 100644 index 0000000000..e0b84d3e84 --- /dev/null +++ b/sysdeps/powerpc/q_qtoull.c @@ -0,0 +1,65 @@ +/* 128-bit floating point to 64-bit fixed point unsigned integer. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <quad_float.h> + +/* unsigned long long _q_qtoull(const long double *a); + Convert 'a' to unsigned long long by truncation. + + Special cases: + -0: return 0 + NaN: raise VXCVI, return 0 + SNaN: raise VXSNAN, VXCVI, and return 0 + Negative numbers (other than -0, but including -Inf): raise VXCVI, return 0 + >= 2^32: raise VXCVI, return 0xffffffffffffffff + */ + +unsigned long long +__q_qtoull(const unsigned long long a[2]) +{ + int ax; + unsigned result; + unsigned long long a0, a1; + + /* Deal with non-special cases. */ + a0 = a[0]; a1 = a[1]; + ax = (a0 >> 48) - 16383; + if (ax <= 63) + /* If we had a '>>' operator that behaved properly with respect to + large shifts, we wouldn't need the '& -(ax >> 31)' term, which + clears 'result' if ax is negative. */ + return ((a0 << 15 | a1 >> 64-15 | 0x8000000000000000ULL) >> 63-ax + & -(unsigned long long)(ax >> 31)); + + /* 'result' is 1 if a is negative, 0 otherwise. */ + result = a0 >> 63; + /* Check for -0, otherwise raise VXCVI. */ + if (a0 != 0x8000000000000000ULL || a[1] != 0) + { + /* Check for SNaN, if so raise VXSNAN. */ + if ((a0 >> 32 & 0x7fff8000) == 0x7fff0000 + && (a1 | a0 & 0x00007fffffffffffULL) != 0) + set_fpscr_bit(FPSCR_VXSNAN); + /* Treat all NaNs as large negative numbers. */ + result |= cntlzw(~(a0 >> 32) & 0x7fff0000) >> 5; + + set_fpscr_bit(FPSCR_VXCVI); + } + return result-1LL; +} diff --git a/sysdeps/powerpc/q_stoq.c b/sysdeps/powerpc/q_stoq.c new file mode 100644 index 0000000000..766997b97b --- /dev/null +++ b/sysdeps/powerpc/q_stoq.c @@ -0,0 +1,65 @@ +/* 32-bit floating point to 128-bit floating point. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <quad_float.h> + +/* long double _q_stoq(float a); + Convert 'a' to long double. Don't raise exceptions. + */ + +void +__q_stoq(unsigned long long result[2], float a) +{ + unsigned ux, rx; + union { + float d; + unsigned u; + } u; + + u.d = a; + result[1] = 0; + result[0] = u.u << 32-7; + /* correct exponent bias */ + rx = ((int)u.u >> 23 & 0x80ff) + 16383-127; + + ux = u.u >> 23 & 0xff; + if (ux == 0) + { + if ((u.u & 0x7fffffff) == 0) + { + /* +0.0 or -0.0. */ + rx &= 0x8000; + } + else + { + /* Denormalised number. Renormalise. */ + unsigned um = u.u & 0x007fffff; + int cs = cntlzw(um) - 9 + 1; + rx -= cs; + um <<= cs; + result[0] = um << 32-7; + } + } + else if (ux == 0xff) + { + /* Inf or NaN. */ + rx |= 0x7fff; + } + *(unsigned short *)result = rx; +} diff --git a/sysdeps/powerpc/q_ulltoq.c b/sysdeps/powerpc/q_ulltoq.c new file mode 100644 index 0000000000..caae124781 --- /dev/null +++ b/sysdeps/powerpc/q_ulltoq.c @@ -0,0 +1,46 @@ +/* 64-bit fixed point unsigned integer to 128-bit floating point. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <quad_float.h> + +/* long double _q_ulltoq(unsigned long long a); + Convert 'a' to long double. + */ + +void +__q_ulltoq(unsigned int result[4], unsigned long long a) +{ + unsigned rx; + int cs, css; + + /* Normalize. */ + cs = cntlzd(a); + a <<= cs+1; + + /* Calculate the exponent, in 4 easy instructions. */ + css = 63-cs; + rx = 16383+css << 16 & ~css; + + /* Put it all together. */ + result[3] = 0; + result[0] = rx & 0xffff0000 | a >> 48 & 0x0000ffff; + a <<= 16; + result[2] = a; + result[1] = a >> 32; +} diff --git a/sysdeps/powerpc/q_utoq.c b/sysdeps/powerpc/q_utoq.c new file mode 100644 index 0000000000..84ce215a00 --- /dev/null +++ b/sysdeps/powerpc/q_utoq.c @@ -0,0 +1,44 @@ +/* 32-bit fixed point unsigned integer to 128-bit floating point. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <quad_float.h> + +/* long double _q_utoq(unsigned int a); + Convert 'a' to long double. + */ + +void +__q_utoq(unsigned int result[4], unsigned int a) +{ + unsigned rx; + int cs, css; + + /* Normalize. */ + cs = cntlzw(a); + a <<= cs+1; + + /* Calculate the exponent, in 4 easy instructions. */ + css = 31-cs; + rx = 16383+css << 16 & ~css; + + /* Put it all together. */ + result[2] = result[3] = 0; + asm ("rlwimi %0,%1,16,16,31": "=r"(result[0]) : "r"(a), "0"(rx)); + result[1] = a << 16; +} diff --git a/sysdeps/powerpc/quad_float.h b/sysdeps/powerpc/quad_float.h new file mode 100644 index 0000000000..787e8d9ab9 --- /dev/null +++ b/sysdeps/powerpc/quad_float.h @@ -0,0 +1,44 @@ +/* Internal libc stuff for 128-bit IEEE FP emulation routines. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#ifndef _QUAD_FLOAT_H +#define _QUAD_FLOAT_H 1 + +#include <fenv_libc.h> + +/* Returns the number of leading zero bits in 'x' (between 0 and 32 + inclusive). 'x' is treated as being 32 bits long. */ +#define cntlzw(x) \ + ({ unsigned p = (x); \ + unsigned r; \ + asm ("cntlzw %0,%1" : "=r"(r) : "r"(p)); \ + r; }) + +/* Returns the number of leading zero bits in 'x' (between 0 and 64 + inclusive). 'x' is treated as being 64 bits long. */ +#define cntlzd(x) \ + ({ unsigned long long q = (x); \ + unsigned int c1, c2; \ + c1 = cntlzw(q >> 32); \ + c2 = cntlzw(q); \ + c1 + (-(c1 >> 5) & c2); }) + +#define shift_and_or + +#endif /* quad_float.h */ diff --git a/sysdeps/powerpc/s_copysign.S b/sysdeps/powerpc/s_copysign.S new file mode 100644 index 0000000000..adc7df226a --- /dev/null +++ b/sysdeps/powerpc/s_copysign.S @@ -0,0 +1,60 @@ +/* Copy a sign bit between floating-point values. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +/* This has been coded in assembler because GCC makes such a mess of it + when it's coded in C. */ + + .section ".text" + .align 2 + .globl __copysign + .type __copysign,@function +__copysign: +/* double [f1] copysign (double [f1] x, double [f2] y); + copysign(x,y) returns a value with the magnitude of x and + with the sign bit of y. */ + + stwu %r1,-16(%r1) + stfd %f2,8(%r1) + lwz %r3,8(%r1) + cmpwi %r3,0 + addi %r1,%r1,16 + blt 0f + fabs %f1,%f1 + blr +0: fnabs %f1,%f1 + blr +0: + .size __copysign,0b-__copysign + + .globl copysign + .globl copysignf + .globl __copysignf + .weak copysign + .weak copysignf + .set copysign,__copysign +/* It turns out that it's safe to use this code even for single-precision. */ + .set __copysignf,__copysign + .set copysignf,__copysign +#ifdef NO_LONG_DOUBLE + .globl copysignl + .globl __copysignl + .weak copysignl + .set __copysignl,__copysign + .set copysignl,__copysign +#endif diff --git a/sysdeps/powerpc/s_fabs.S b/sysdeps/powerpc/s_fabs.S new file mode 100644 index 0000000000..a52733568d --- /dev/null +++ b/sysdeps/powerpc/s_fabs.S @@ -0,0 +1,42 @@ +/* Floating-point absolute value. PowerPC version. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + + .section ".text" + .align 2 + .globl __fabs + .type __fabs,@function +__fabs: +/* double [f1] fabs (double [f1] x); */ + fabs %f1,%f1 + blr +0: + .size __fabs,0b-__fabs + + .globl fabs,fabsf,__fabsf + .weak fabs,fabsf + .set fabs,__fabs +/* It turns out that it's safe to use this code even for single-precision. */ + .set __fabsf,__fabs + .set fabsf,__fabs +#ifdef NO_LONG_DOUBLE + .globl fabsl,__fabsl + .weak fabsl + .set __fabsl,__fabs + .set fabsl,__fabs +#endif diff --git a/sysdeps/powerpc/s_isnan.c b/sysdeps/powerpc/s_isnan.c new file mode 100644 index 0000000000..34019fdbdb --- /dev/null +++ b/sysdeps/powerpc/s_isnan.c @@ -0,0 +1,47 @@ +/* Return 1 if argument is a NaN, else 0. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include "math.h" +#include <fenv_libc.h> + +int __isnan(double x) +{ + fenv_t savedstate; + int result; + savedstate = fegetenv_register(); + reset_fpscr_bit(FPSCR_VE); + result = !(x == x); + fesetenv_register(savedstate); + return result; +} +weak_alias (__isnan, isnan) +/* It turns out that the 'double' version will also always work for + single-precision. Use explicit assembler to stop gcc complaining + that 'isnanf' takes a float parameter, not double. */ +asm ("\ + .globl __isnanf + .globl isnanf + .weak isnanf + .set __isnanf,__isnan + .set isnanf,__isnan +"); +#ifdef NO_LONG_DOUBLE +strong_alias (__isnan, __isnanl) +weak_alias (__isnan, isnanl) +#endif diff --git a/sysdeps/powerpc/s_llrint.c b/sysdeps/powerpc/s_llrint.c new file mode 100644 index 0000000000..7ca48c029a --- /dev/null +++ b/sysdeps/powerpc/s_llrint.c @@ -0,0 +1,28 @@ +/* Round a long long floating point value to an integer in the current + rounding mode. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include "math.h" + +long long int +__llrint (long double x) +{ + return (long long int) __rintl (x); +} +weak_alias (__llrint, llrint) diff --git a/sysdeps/powerpc/s_llround.c b/sysdeps/powerpc/s_llround.c new file mode 100644 index 0000000000..fbe3a3217c --- /dev/null +++ b/sysdeps/powerpc/s_llround.c @@ -0,0 +1,46 @@ +/* Round long double value to long int. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <math.h> + +/* I think that what this routine is supposed to do is round a value + to the nearest integer, with values exactly on the boundary rounded + away from zero. */ +/* This routine relies on (long long)x, when x is out of range of a long long, + clipping to MAX_LLONG or MIN_LLONG. */ + +long long int +__llround (long double x) +{ + long double xrf; + long long int xr; + xr = (long long int) x; + xrf = (long double) xr; + if (x >= 0.0) + if (x - xrf >= 0.5 && x - xrf < 1.0 && x+1 > 0) + return x+1; + else + return x; + else + if (xrf - x >= 0.5 && xrf - x < 1.0 && x-1 < 0) + return x-1; + else + return x; +} +weak_alias (__llround, llround) diff --git a/sysdeps/powerpc/s_lrint.c b/sysdeps/powerpc/s_lrint.c new file mode 100644 index 0000000000..647cf30c9c --- /dev/null +++ b/sysdeps/powerpc/s_lrint.c @@ -0,0 +1,44 @@ +/* Round floating-point to integer. PowerPC version. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include "math.h" + +#ifdef NO_LONG_DOUBLE + +long int +__lrint (long double x) +{ + return (long int) __rintl(x); +} + +#else + +long int +__lrint (long double x) +{ + union { + double d; + long int ll[2]; + } u; + asm ("fctiw %0,%1" : "=f"(u.d) : "f"(x)); + return d.ll[1]; +} + +#endif +weak_alias (__lrint, lrint) diff --git a/sysdeps/powerpc/s_lround.c b/sysdeps/powerpc/s_lround.c new file mode 100644 index 0000000000..a6cb6c96d2 --- /dev/null +++ b/sysdeps/powerpc/s_lround.c @@ -0,0 +1,46 @@ +/* Round long double value to long int. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <math.h> + +/* I think that what this routine is supposed to do is round a value + to the nearest integer, with values exactly on the boundary rounded + away from zero. */ +/* This routine relies on (long int)x, when x is out of range of a long int, + clipping to MAX_LONG or MIN_LONG. */ + +long int +__lround (long double x) +{ + long double xrf; + long int xr; + xr = (long int) x; + xrf = (long double) xr; + if (x >= 0.0) + if (x - xrf >= 0.5 && x - xrf < 1.0 && x+1 > 0) + return x+1; + else + return x; + else + if (xrf - x >= 0.5 && xrf - x < 1.0 && x-1 < 0) + return x-1; + else + return x; +} +weak_alias (__lround, lround) diff --git a/sysdeps/powerpc/s_rint.c b/sysdeps/powerpc/s_rint.c new file mode 100644 index 0000000000..d6d1dd2f76 --- /dev/null +++ b/sysdeps/powerpc/s_rint.c @@ -0,0 +1,45 @@ +/* Round a 64-bit floating point value to the nearest integer. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include "math.h" + +double +__rint (double x) +{ + static const float TWO52 = 4503599627370496.0; + + if (fabs (x) < TWO52) + if (x > 0.0) + { + x += TWO52; + x -= TWO52; + } + else if (x < 0.0) + { + x -= TWO52; + x += TWO52; + } + + return x; +} +weak_alias (__rint, rint) +#ifdef NO_LONG_DOUBLE +strong_alias (__rint, __rintl) +weak_alias (__rint, rintl) +#endif diff --git a/sysdeps/powerpc/s_rintf.c b/sysdeps/powerpc/s_rintf.c new file mode 100644 index 0000000000..b95c0540bf --- /dev/null +++ b/sysdeps/powerpc/s_rintf.c @@ -0,0 +1,41 @@ +/* Round a 32-bit floating point value to the nearest integer. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include "math.h" + +float +__rintf (float x) +{ + static const float TWO23 = 8388608.0; + + if (fabsf (x) < TWO23) + if (x > 0.0) + { + x += TWO23; + x -= TWO23; + } + else if (x < 0.0) + { + x -= TWO23; + x += TWO23; + } + + return x; +} +weak_alias (__rintf, rintf) diff --git a/sysdeps/powerpc/strlen.s b/sysdeps/powerpc/strlen.s index 9d671cabc3..ea809772eb 100644 --- a/sysdeps/powerpc/strlen.s +++ b/sysdeps/powerpc/strlen.s @@ -61,8 +61,8 @@ # Some notes on register usage: Under the SVR4 ABI, we can use registers # 0 and 3 through 12 (so long as we don't call any procedures) without # saving them. We can also use registers 14 through 31 if we save them. - # We can't use r1 (it's the stack pointer), nor r2 or r13 because the user - # program may expect them to be hold their usual value if we get sent + # We can't use r1 (it's the stack pointer), r2 nor r13 because the user + # program may expect them to hold their usual value if we get sent # a signal. Integer parameters are passed in r3 through r10. # We can use condition registers cr0, cr1, cr5, cr6, and cr7 without saving # them, the others we must save. @@ -80,14 +80,11 @@ strlen: # r9-12 are temporaries. r0 is used as a temporary and for discarded # results. clrrwi %r4,%r3,2 - lis %r6,0xfeff lis %r7,0x7f7f - rlwinm %r10,%r3,0,29,29 - lwz %r8,0(%r4) - addi %r7,%r7,0x7f7f rlwinm %r5,%r3,3,27,28 - cmpwi %cr1,%r10,0 + lwz %r8,0(%r4) li %r9,-1 + addi %r7,%r7,0x7f7f # That's the setup done, now do the first pair of words. # We make an exception and use method (2) on the first two words, to reduce # overhead. @@ -97,12 +94,14 @@ strlen: add %r0,%r0,%r7 nor %r0,%r10,%r0 and. %r8,%r0,%r9 + mtcrf 0x01,%r3 bne done0 - # Handle second word of pair. Put addi between branches to avoid hurting - # branch prediction. + lis %r6,0xfeff addi %r6,%r6,-0x101 + # Are we now aligned to a doubleword boundary? + bt 29,loop - bne %cr1,loop + # Handle second word of pair. lwzu %r8,4(%r4) and %r0,%r7,%r8 or %r10,%r7,%r8 diff --git a/sysdeps/powerpc/t_sqrt.c b/sysdeps/powerpc/t_sqrt.c new file mode 100644 index 0000000000..c49380c0fd --- /dev/null +++ b/sysdeps/powerpc/t_sqrt.c @@ -0,0 +1,144 @@ +const float __t_sqrt[1024] = { +0.7078,0.7064, 0.7092,0.7050, 0.7106,0.7037, 0.7119,0.7023, 0.7133,0.7010, +0.7147,0.6996, 0.7160,0.6983, 0.7174,0.6970, 0.7187,0.6957, 0.7201,0.6943, +0.7215,0.6930, 0.7228,0.6917, 0.7242,0.6905, 0.7255,0.6892, 0.7269,0.6879, +0.7282,0.6866, 0.7295,0.6854, 0.7309,0.6841, 0.7322,0.6829, 0.7335,0.6816, +0.7349,0.6804, 0.7362,0.6792, 0.7375,0.6779, 0.7388,0.6767, 0.7402,0.6755, +0.7415,0.6743, 0.7428,0.6731, 0.7441,0.6719, 0.7454,0.6708, 0.7467,0.6696, +0.7480,0.6684, 0.7493,0.6672, 0.7507,0.6661, 0.7520,0.6649, 0.7532,0.6638, +0.7545,0.6627, 0.7558,0.6615, 0.7571,0.6604, 0.7584,0.6593, 0.7597,0.6582, +0.7610,0.6570, 0.7623,0.6559, 0.7635,0.6548, 0.7648,0.6537, 0.7661,0.6527, +0.7674,0.6516, 0.7686,0.6505, 0.7699,0.6494, 0.7712,0.6484, 0.7725,0.6473, +0.7737,0.6462, 0.7750,0.6452, 0.7762,0.6441, 0.7775,0.6431, 0.7787,0.6421, +0.7800,0.6410, 0.7812,0.6400, 0.7825,0.6390, 0.7837,0.6380, 0.7850,0.6370, +0.7862,0.6359, 0.7875,0.6349, 0.7887,0.6339, 0.7900,0.6330, 0.7912,0.6320, +0.7924,0.6310, 0.7937,0.6300, 0.7949,0.6290, 0.7961,0.6281, 0.7973,0.6271, +0.7986,0.6261, 0.7998,0.6252, 0.8010,0.6242, 0.8022,0.6233, 0.8034,0.6223, +0.8046,0.6214, 0.8059,0.6205, 0.8071,0.6195, 0.8083,0.6186, 0.8095,0.6177, +0.8107,0.6168, 0.8119,0.6158, 0.8131,0.6149, 0.8143,0.6140, 0.8155,0.6131, +0.8167,0.6122, 0.8179,0.6113, 0.8191,0.6104, 0.8203,0.6096, 0.8215,0.6087, +0.8227,0.6078, 0.8238,0.6069, 0.8250,0.6060, 0.8262,0.6052, 0.8274,0.6043, +0.8286,0.6035, 0.8297,0.6026, 0.8309,0.6017, 0.8321,0.6009, 0.8333,0.6000, +0.8344,0.5992, 0.8356,0.5984, 0.8368,0.5975, 0.8379,0.5967, 0.8391,0.5959, +0.8403,0.5950, 0.8414,0.5942, 0.8426,0.5934, 0.8437,0.5926, 0.8449,0.5918, +0.8461,0.5910, 0.8472,0.5902, 0.8484,0.5894, 0.8495,0.5886, 0.8507,0.5878, +0.8518,0.5870, 0.8530,0.5862, 0.8541,0.5854, 0.8552,0.5846, 0.8564,0.5838, +0.8575,0.5831, 0.8587,0.5823, 0.8598,0.5815, 0.8609,0.5808, 0.8621,0.5800, +0.8632,0.5792, 0.8643,0.5785, 0.8655,0.5777, 0.8666,0.5770, 0.8677,0.5762, +0.8688,0.5755, 0.8700,0.5747, 0.8711,0.5740, 0.8722,0.5733, 0.8733,0.5725, +0.8744,0.5718, 0.8756,0.5711, 0.8767,0.5703, 0.8778,0.5696, 0.8789,0.5689, +0.8800,0.5682, 0.8811,0.5675, 0.8822,0.5667, 0.8833,0.5660, 0.8844,0.5653, +0.8855,0.5646, 0.8866,0.5639, 0.8877,0.5632, 0.8888,0.5625, 0.8899,0.5618, +0.8910,0.5611, 0.8921,0.5605, 0.8932,0.5598, 0.8943,0.5591, 0.8954,0.5584, +0.8965,0.5577, 0.8976,0.5570, 0.8987,0.5564, 0.8998,0.5557, 0.9008,0.5550, +0.9019,0.5544, 0.9030,0.5537, 0.9041,0.5530, 0.9052,0.5524, 0.9062,0.5517, +0.9073,0.5511, 0.9084,0.5504, 0.9095,0.5498, 0.9105,0.5491, 0.9116,0.5485, +0.9127,0.5478, 0.9138,0.5472, 0.9148,0.5465, 0.9159,0.5459, 0.9170,0.5453, +0.9180,0.5446, 0.9191,0.5440, 0.9202,0.5434, 0.9212,0.5428, 0.9223,0.5421, +0.9233,0.5415, 0.9244,0.5409, 0.9254,0.5403, 0.9265,0.5397, 0.9276,0.5391, +0.9286,0.5384, 0.9297,0.5378, 0.9307,0.5372, 0.9318,0.5366, 0.9328,0.5360, +0.9338,0.5354, 0.9349,0.5348, 0.9359,0.5342, 0.9370,0.5336, 0.9380,0.5330, +0.9391,0.5324, 0.9401,0.5319, 0.9411,0.5313, 0.9422,0.5307, 0.9432,0.5301, +0.9442,0.5295, 0.9453,0.5289, 0.9463,0.5284, 0.9473,0.5278, 0.9484,0.5272, +0.9494,0.5266, 0.9504,0.5261, 0.9515,0.5255, 0.9525,0.5249, 0.9535,0.5244, +0.9545,0.5238, 0.9556,0.5233, 0.9566,0.5227, 0.9576,0.5221, 0.9586,0.5216, +0.9596,0.5210, 0.9607,0.5205, 0.9617,0.5199, 0.9627,0.5194, 0.9637,0.5188, +0.9647,0.5183, 0.9657,0.5177, 0.9667,0.5172, 0.9677,0.5167, 0.9687,0.5161, +0.9698,0.5156, 0.9708,0.5151, 0.9718,0.5145, 0.9728,0.5140, 0.9738,0.5135, +0.9748,0.5129, 0.9758,0.5124, 0.9768,0.5119, 0.9778,0.5114, 0.9788,0.5108, +0.9798,0.5103, 0.9808,0.5098, 0.9818,0.5093, 0.9828,0.5088, 0.9838,0.5083, +0.9847,0.5077, 0.9857,0.5072, 0.9867,0.5067, 0.9877,0.5062, 0.9887,0.5057, +0.9897,0.5052, 0.9907,0.5047, 0.9917,0.5042, 0.9926,0.5037, 0.9936,0.5032, +0.9946,0.5027, 0.9956,0.5022, 0.9966,0.5017, 0.9976,0.5012, 0.9985,0.5007, +0.9995,0.5002, +1.0010,0.4995, 1.0029,0.4985, 1.0049,0.4976, 1.0068,0.4966, 1.0088,0.4957, +1.0107,0.4947, 1.0126,0.4938, 1.0145,0.4928, 1.0165,0.4919, 1.0184,0.4910, +1.0203,0.4901, 1.0222,0.4891, 1.0241,0.4882, 1.0260,0.4873, 1.0279,0.4864, +1.0298,0.4855, 1.0317,0.4846, 1.0336,0.4837, 1.0355,0.4829, 1.0374,0.4820, +1.0393,0.4811, 1.0411,0.4802, 1.0430,0.4794, 1.0449,0.4785, 1.0468,0.4777, +1.0486,0.4768, 1.0505,0.4760, 1.0523,0.4751, 1.0542,0.4743, 1.0560,0.4735, +1.0579,0.4726, 1.0597,0.4718, 1.0616,0.4710, 1.0634,0.4702, 1.0653,0.4694, +1.0671,0.4686, 1.0689,0.4678, 1.0707,0.4670, 1.0726,0.4662, 1.0744,0.4654, +1.0762,0.4646, 1.0780,0.4638, 1.0798,0.4630, 1.0816,0.4623, 1.0834,0.4615, +1.0852,0.4607, 1.0870,0.4600, 1.0888,0.4592, 1.0906,0.4585, 1.0924,0.4577, +1.0942,0.4570, 1.0960,0.4562, 1.0978,0.4555, 1.0995,0.4547, 1.1013,0.4540, +1.1031,0.4533, 1.1049,0.4525, 1.1066,0.4518, 1.1084,0.4511, 1.1101,0.4504, +1.1119,0.4497, 1.1137,0.4490, 1.1154,0.4483, 1.1172,0.4476, 1.1189,0.4469, +1.1207,0.4462, 1.1224,0.4455, 1.1241,0.4448, 1.1259,0.4441, 1.1276,0.4434, +1.1293,0.4427, 1.1311,0.4421, 1.1328,0.4414, 1.1345,0.4407, 1.1362,0.4401, +1.1379,0.4394, 1.1397,0.4387, 1.1414,0.4381, 1.1431,0.4374, 1.1448,0.4368, +1.1465,0.4361, 1.1482,0.4355, 1.1499,0.4348, 1.1516,0.4342, 1.1533,0.4335, +1.1550,0.4329, 1.1567,0.4323, 1.1584,0.4316, 1.1600,0.4310, 1.1617,0.4304, +1.1634,0.4298, 1.1651,0.4292, 1.1668,0.4285, 1.1684,0.4279, 1.1701,0.4273, +1.1718,0.4267, 1.1734,0.4261, 1.1751,0.4255, 1.1768,0.4249, 1.1784,0.4243, +1.1801,0.4237, 1.1817,0.4231, 1.1834,0.4225, 1.1850,0.4219, 1.1867,0.4213, +1.1883,0.4208, 1.1900,0.4202, 1.1916,0.4196, 1.1932,0.4190, 1.1949,0.4185, +1.1965,0.4179, 1.1981,0.4173, 1.1998,0.4167, 1.2014,0.4162, 1.2030,0.4156, +1.2046,0.4151, 1.2063,0.4145, 1.2079,0.4139, 1.2095,0.4134, 1.2111,0.4128, +1.2127,0.4123, 1.2143,0.4117, 1.2159,0.4112, 1.2175,0.4107, 1.2192,0.4101, +1.2208,0.4096, 1.2224,0.4090, 1.2239,0.4085, 1.2255,0.4080, 1.2271,0.4075, +1.2287,0.4069, 1.2303,0.4064, 1.2319,0.4059, 1.2335,0.4054, 1.2351,0.4048, +1.2366,0.4043, 1.2382,0.4038, 1.2398,0.4033, 1.2414,0.4028, 1.2429,0.4023, +1.2445,0.4018, 1.2461,0.4013, 1.2477,0.4008, 1.2492,0.4003, 1.2508,0.3998, +1.2523,0.3993, 1.2539,0.3988, 1.2555,0.3983, 1.2570,0.3978, 1.2586,0.3973, +1.2601,0.3968, 1.2617,0.3963, 1.2632,0.3958, 1.2648,0.3953, 1.2663,0.3949, +1.2678,0.3944, 1.2694,0.3939, 1.2709,0.3934, 1.2725,0.3929, 1.2740,0.3925, +1.2755,0.3920, 1.2771,0.3915, 1.2786,0.3911, 1.2801,0.3906, 1.2816,0.3901, +1.2832,0.3897, 1.2847,0.3892, 1.2862,0.3887, 1.2877,0.3883, 1.2892,0.3878, +1.2907,0.3874, 1.2923,0.3869, 1.2938,0.3865, 1.2953,0.3860, 1.2968,0.3856, +1.2983,0.3851, 1.2998,0.3847, 1.3013,0.3842, 1.3028,0.3838, 1.3043,0.3834, +1.3058,0.3829, 1.3073,0.3825, 1.3088,0.3820, 1.3103,0.3816, 1.3118,0.3812, +1.3132,0.3807, 1.3147,0.3803, 1.3162,0.3799, 1.3177,0.3794, 1.3192,0.3790, +1.3207,0.3786, 1.3221,0.3782, 1.3236,0.3778, 1.3251,0.3773, 1.3266,0.3769, +1.3280,0.3765, 1.3295,0.3761, 1.3310,0.3757, 1.3324,0.3753, 1.3339,0.3748, +1.3354,0.3744, 1.3368,0.3740, 1.3383,0.3736, 1.3397,0.3732, 1.3412,0.3728, +1.3427,0.3724, 1.3441,0.3720, 1.3456,0.3716, 1.3470,0.3712, 1.3485,0.3708, +1.3499,0.3704, 1.3514,0.3700, 1.3528,0.3696, 1.3542,0.3692, 1.3557,0.3688, +1.3571,0.3684, 1.3586,0.3680, 1.3600,0.3676, 1.3614,0.3673, 1.3629,0.3669, +1.3643,0.3665, 1.3657,0.3661, 1.3672,0.3657, 1.3686,0.3653, 1.3700,0.3650, +1.3714,0.3646, 1.3729,0.3642, 1.3743,0.3638, 1.3757,0.3634, 1.3771,0.3631, +1.3785,0.3627, 1.3800,0.3623, 1.3814,0.3620, 1.3828,0.3616, 1.3842,0.3612, +1.3856,0.3609, 1.3870,0.3605, 1.3884,0.3601, 1.3898,0.3598, 1.3912,0.3594, +1.3926,0.3590, 1.3940,0.3587, 1.3954,0.3583, 1.3968,0.3580, 1.3982,0.3576, +1.3996,0.3572, 1.4010,0.3569, 1.4024,0.3565, 1.4038,0.3562, 1.4052,0.3558, +1.4066,0.3555, 1.4080,0.3551, 1.4094,0.3548, 1.4108,0.3544, 1.4121,0.3541, +1.4135,0.3537 +}; + + +/* Generated by: */ +#if 0 +#include <math.h> +#include <stdio.h> +#include <assert.h> + +int +main(int argc, char **argv) +{ + int i, j; + + printf ("const float __t_sqrt[1024] = {"); + for (i = 0; i < 2; i++) + { + putchar('\n'); + for (j = 0; j < 256; j++) + { + double mval = j/512.0 + 0.5; + double eval = i==0 ? 1.0 : 2.0; + double ls = sqrt(mval*eval); + double hs = sqrt((mval+1/512.0)*eval); + double as = (ls+hs)*0.5; + double lx = 1/(2.0*ls); + double hx = 1/(2.0*hs); + double ax = (lx+hx)*0.5; + + printf("%.4f,%.4f%s",as,ax, + i*j==255 ? "\n" : j % 5 == 4 ? ",\n" : ", "); + assert((hs-ls)/as < 1/256.0); + assert((hx-lx)/ax < 1/256.0); + } + } + printf ("};\n"); + return 0; +} +#endif /* 0 */ diff --git a/sysdeps/powerpc/test-arith.c b/sysdeps/powerpc/test-arith.c new file mode 100644 index 0000000000..c846b0d0ec --- /dev/null +++ b/sysdeps/powerpc/test-arith.c @@ -0,0 +1,604 @@ +/* Test floating-point arithmetic operations. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ +#ifndef _GNU_SOURCE +#define _GNU_SOURCE +#endif +#include <math.h> +#include <stdio.h> +#include <stdlib.h> +#include <string.h> +#include <fenv.h> +#include <assert.h> + +#ifndef ESIZE +typedef double tocheck_t; +#define ESIZE 11 +#define MSIZE 52 +#define FUNC(x) x +#endif + +#define R_NEAREST 1 +#define R_ZERO 2 +#define R_UP 4 +#define R_DOWN 8 +#define R_ALL (R_NEAREST|R_ZERO|R_UP|R_DOWN) +static fenv_t rmodes[4]; +static const char * const rmnames[4] = +{ "nearest","zero","+Inf","-Inf" }; + +typedef union { + tocheck_t tc; + unsigned char c[sizeof(tocheck_t)]; +} union_t; + +/* Don't try reading these in a font that doesn't distinguish + O and zero. */ +typedef enum { + P_Z = 0x0, /* 00000...0 */ + P_000O = 0x1, /* 00011...1 */ + P_001Z = 0x2, /* 00100...0 */ + P_00O = 0x3, /* 00111...1 */ + P_01Z = 0x4, /* 01000...0 */ + P_010O = 0x5, /* 01011...1 */ + P_011Z = 0x6, /* 01100...0 */ + P_0O = 0x7, /* 01111...1 */ + P_1Z = 0x8, /* 10000...0 */ + P_100O = 0x9, /* 10011...1 */ + P_101Z = 0xa, /* 10100...0 */ + P_10O = 0xb, /* 10111...1 */ + P_11Z = 0xc, /* 11000...0 */ + P_110O = 0xd, /* 11011...1 */ + P_111Z = 0xe, /* 11100...0 */ + P_O = 0xf, /* 11111...1 */ + P_Z1 = 0x11, /* 000...001 */ + P_Z10 = 0x12, /* 000...010 */ + P_Z11 = 0x13, /* 000...011 */ + P_0O00 = 0x14, /* 011...100 */ + P_0O01 = 0x15, /* 011...101 */ + P_0O0 = 0x16, /* 011...110 */ + P_1Z1 = 0x19, /* 100...001 */ + P_1Z10 = 0x1a, /* 100...010 */ + P_1Z11 = 0x1b, /* 100...011 */ + P_O00 = 0x1c, /* 111...100 */ + P_O01 = 0x1d, /* 111...101 */ + P_O0 = 0x1e, /* 111...110 */ + P_R = 0x20, /* rrr...rrr */ /* ('r' means random. ) */ + P_Ro = 0x21, /* rrr...rrr, with odd parity. */ + P_0R = 0x22, /* 0rr...rrr */ + P_1R = 0x23, /* 1rr...rrr */ + P_Rno = 0x24, /* rrr...rrr, but not all ones. */ +} pattern_t; + +static void +pattern_fill(pattern_t ptn, unsigned char *start, int bitoffset, int count) +{ +#define bitset(count, value) \ + start[(count)/8] = (start[(count)/8] & ~(1 << 7-(count)%8) \ + | (value) << 7-(count)%8) + int i; + + if (ptn >= 0 && ptn <= 0xf) + { + /* Patterns between 0 and 0xF have the following format: + The LSBit is used to fill the last n-3 bits of the pattern; + The next 3 bits are the first 3 bits of the pattern. */ + for (i = 0; i < count; i++) + if (i < 3) + bitset((bitoffset+i), ptn >> (3-i) & 1); + else + bitset((bitoffset+i), ptn >> 0 & 1); + } + else if (ptn <= 0x1f) + { + /* Patterns between 0x10 and 0x1F have the following format: + The two LSBits are the last two bits of the pattern; + The 0x8 bit is the first bit of the pattern; + The 0x4 bit is used to fill the remainder. */ + for (i = 0; i < count; i++) + if (i == 0) + bitset((bitoffset+i), ptn >> 3 & 1); + else if (i >= count-2) + bitset((bitoffset+i), ptn >> (count-1-i) & 1); + else + bitset((bitoffset+i), ptn >> 2 & 1); + } + else switch (ptn) + { + case P_0R: case P_1R: + assert(count > 0); + bitset(bitoffset, ptn & 1); + count--; + bitoffset++; + case P_R: + for (; count > 0; count--, bitoffset++) + bitset(bitoffset, rand() & 1); + break; + case P_Ro: + { + int op = 1; + assert(count > 0); + for (; count > 1; count--, bitoffset++) + bitset(bitoffset, op ^= (rand() & 1)); + bitset(bitoffset, op); + break; + } + case P_Rno: + { + int op = 1; + assert(count > 0); + for (; count > 1; count--, bitoffset++) + { + int r = rand() & 1; + op &= r; + bitset(bitoffset, r); + } + bitset(bitoffset, rand() & (op ^ 1)); + break; + } + + default: + assert(0); + } +#undef bitset +} + +static tocheck_t +pattern(int negative, pattern_t exp, pattern_t mant) +{ + union_t result; +#if 0 + int i; +#endif + + pattern_fill(negative ? P_O : P_Z, result.c, 0, 1); + pattern_fill(exp, result.c, 1, ESIZE); + pattern_fill(mant, result.c, ESIZE+1, MSIZE); +#if 0 + printf("neg=%d exp=%02x mant=%02x: ", negative, exp, mant); + for (i = 0; i < sizeof(tocheck_t); i++) + printf("%02x", result.c[i]); + printf("\n"); +#endif + return result.tc; +} + +/* Return the closest different tocheck_t to 'x' in the direction of + 'direction', or 'x' if there is no such value. Assumes 'x' is not + a NaN. */ +static tocheck_t +delta(tocheck_t x, int direction) +{ + union_t xx; + int i; + + xx.tc = x; + if (xx.c[0] & 0x80) + direction = -direction; + if (direction == +1) + { + union_t tx; + tx.tc = pattern(xx.c[0] >> 7, P_O, P_Z); + if (memcmp(tx.c, xx.c, sizeof(tocheck_t)) == 0) + return x; + } + for (i = sizeof(tocheck_t)-1; i > 0; i--) + { + xx.c[i] += direction; + if (xx.c[i] != (direction > 0 ? 0 : 0xff)) + return xx.tc; + } + if (direction < 0 && (xx.c[0] & 0x7f) == 0) + return pattern(~(xx.c[0] >> 7) & 1, P_Z, P_Z1); + else + { + xx.c[0] += direction; + return xx.tc; + } +} + +static int nerrors = 0; + +#ifdef FE_ALL_INVALID +static const int all_exceptions = FE_ALL_INVALID | FE_ALL_EXCEPT; +#else +static const int all_exceptions = FE_ALL_EXCEPT; +#endif + +static void +check_result(int line, const char *rm, tocheck_t expected, tocheck_t actual) +{ + if (memcmp(&expected, &actual, sizeof(tocheck_t)) != 0) + { + unsigned char *ex, *ac; + int i; + + printf("%s:%d:round %s:result failed\n" + " expected result 0x", __FILE__, line, rm); + ex = (unsigned char *)&expected; + ac = (unsigned char *)&actual; + for (i = 0; i < sizeof(tocheck_t); i++) + printf("%02x", ex[i]); + printf(" got 0x"); + for (i = 0; i < sizeof(tocheck_t); i++) + printf("%02x", ac[i]); + printf("\n"); + nerrors++; + } +} + +static const struct { + int except; + const char *name; +} excepts[] = { +#define except_entry(ex) { ex, #ex } , +#ifdef FE_INEXACT + except_entry(FE_INEXACT) +#else +# define FE_INEXACT 0 +#endif +#ifdef FE_DIVBYZERO + except_entry(FE_DIVBYZERO) +#else +# define FE_DIVBYZERO 0 +#endif +#ifdef FE_UNDERFLOW + except_entry(FE_UNDERFLOW) +#else +# define FE_UNDERFLOW 0 +#endif +#ifdef FE_OVERFLOW + except_entry(FE_OVERFLOW) +#else +# define FE_OVERFLOW 0 +#endif +#ifdef FE_INVALID + except_entry(FE_INVALID) +#else +# define FE_INVALID 0 +#endif +#ifdef FE_INVALID_SNAN + except_entry(FE_INVALID_SNAN) +#else +# define FE_INVALID_SNAN FE_INVALID +#endif +#ifdef FE_INVALID_ISI + except_entry(FE_INVALID_ISI) +#else +# define FE_INVALID_ISI FE_INVALID +#endif +#ifdef FE_INVALID_IDI + except_entry(FE_INVALID_IDI) +#else +# define FE_INVALID_IDI FE_INVALID +#endif +#ifdef FE_INVALID_ZDZ + except_entry(FE_INVALID_ZDZ) +#else +# define FE_INVALID_ZDZ FE_INVALID +#endif +#ifdef FE_INVALID_COMPARE + except_entry(FE_INVALID_COMPARE) +#else +# define FE_INVALID_COMPARE FE_INVALID +#endif +#ifdef FE_INVALID_SOFTWARE + except_entry(FE_INVALID_SOFTWARE) +#else +# define FE_INVALID_SOFTWARE FE_INVALID +#endif +#ifdef FE_INVALID_SQRT + except_entry(FE_INVALID_SQRT) +#else +# define FE_INVALID_SQRT FE_INVALID +#endif +#ifdef FE_INVALID_INTEGER_CONVERSION + except_entry(FE_INVALID_INTEGER_CONVERSION) +#else +# define FE_INVALID_INTEGER_CONVERSION FE_INVALID +#endif +}; + +static int excepts_missing = 0; + +static void +check_excepts(int line, const char *rm, int expected, int actual) +{ + if (expected & excepts_missing) + expected = expected & ~excepts_missing | FE_INVALID_SNAN; + if ((expected & all_exceptions) != actual) + { + int i; + printf("%s:%d:round %s:exceptions failed\n" + " expected exceptions ", __FILE__, line,rm); + for (i = 0; i < sizeof(excepts)/sizeof(excepts[0]); i++) + if (expected & excepts[i].except) + printf("%s ",excepts[i].name); + if ((expected & all_exceptions) == 0) + printf("- "); + printf("got"); + for (i = 0; i < sizeof(excepts)/sizeof(excepts[0]); i++) + if (actual & excepts[i].except) + printf(" %s",excepts[i].name); + if ((actual & all_exceptions) == 0) + printf("- "); + printf(".\n"); + nerrors++; + } +} + +typedef enum { + B_ADD, B_SUB, B_MUL, B_DIV, B_NEG, B_ABS, B_SQRT +} op_t; +typedef struct { + int line; + op_t op; + int a_sgn; + pattern_t a_exp, a_mant; + int b_sgn; + pattern_t b_exp, b_mant; + int rmode; + int excepts; + int x_sgn; + pattern_t x_exp, x_mant; +} optest_t; +static const optest_t optests[] = { + /* Additions of zero. */ + {__LINE__,B_ADD, 0,P_Z,P_Z, 0,P_Z,P_Z, R_ALL,0, 0,P_Z,P_Z }, + {__LINE__,B_ADD, 1,P_Z,P_Z, 0,P_Z,P_Z, R_ALL & ~R_DOWN,0, 0,P_Z,P_Z }, + {__LINE__,B_ADD, 1,P_Z,P_Z, 0,P_Z,P_Z, R_DOWN,0, 1,P_Z,P_Z }, + {__LINE__,B_ADD, 1,P_Z,P_Z, 1,P_Z,P_Z, R_ALL,0, 1,P_Z,P_Z }, + + /* Additions with NaN. */ + {__LINE__,B_ADD, 0,P_O,P_101Z, 0,P_Z,P_Z, R_ALL,0, 0,P_O,P_101Z }, + {__LINE__,B_ADD, 0,P_O,P_01Z, 0,P_Z,P_Z, R_ALL, + FE_INVALID | FE_INVALID_SNAN, 0,P_O,P_11Z }, + {__LINE__,B_ADD, 0,P_O,P_Z, 0,P_O,P_0O, R_ALL, + FE_INVALID | FE_INVALID_SNAN, 0,P_O,P_O }, + {__LINE__,B_ADD, 0,P_Z,P_Z, 0,P_O,P_11Z, R_ALL,0, 0,P_O,P_11Z }, + {__LINE__,B_ADD, 0,P_O,P_001Z, 0,P_O,P_001Z, R_ALL, + FE_INVALID | FE_INVALID_SNAN, 0,P_O,P_101Z }, + {__LINE__,B_ADD, 0,P_O,P_1Z, 0,P_Z,P_Z, R_ALL,0, 0,P_O,P_1Z }, + {__LINE__,B_ADD, 0,P_0O,P_Z, 0,P_O,P_10O, R_ALL,0, 0,P_O,P_10O }, + + /* Additions with infinity. */ + {__LINE__,B_ADD, 0,P_O,P_Z, 0,P_Z,P_Z, R_ALL,0, 0,P_O,P_Z }, + {__LINE__,B_ADD, 0,P_O,P_Z, 1,P_Z,P_Z, R_ALL,0, 0,P_O,P_Z }, + {__LINE__,B_ADD, 1,P_O,P_Z, 0,P_Z,P_Z, R_ALL,0, 1,P_O,P_Z }, + {__LINE__,B_ADD, 1,P_O,P_Z, 1,P_Z,P_Z, R_ALL,0, 1,P_O,P_Z }, + {__LINE__,B_ADD, 0,P_O,P_Z, 0,P_O,P_Z, R_ALL,0, 0,P_O,P_Z }, + {__LINE__,B_ADD, 1,P_O,P_Z, 1,P_O,P_Z, R_ALL,0, 1,P_O,P_Z }, + {__LINE__,B_ADD, 0,P_O,P_Z, 1,P_O,P_Z, R_ALL, + FE_INVALID | FE_INVALID_ISI, 0,P_O,P_1Z }, + {__LINE__,B_ADD, 1,P_O,P_Z, 0,P_O,P_Z, R_ALL, + FE_INVALID | FE_INVALID_ISI, 0,P_O,P_1Z }, + {__LINE__,B_ADD, 0,P_O,P_Z, 0,P_0O,P_Z, R_ALL,0, 0,P_O,P_Z }, + {__LINE__,B_ADD, 1,P_O,P_Z, 0,P_0O,P_Z, R_ALL,0, 1,P_O,P_Z }, + {__LINE__,B_ADD, 0,P_O,P_Z, 1,P_0O,P_Z, R_ALL,0, 0,P_O,P_Z }, + {__LINE__,B_ADD, 1,P_O,P_Z, 1,P_0O,P_Z, R_ALL,0, 1,P_O,P_Z }, + + /* Overflow (and zero). */ + {__LINE__,B_ADD, 0,P_O0,P_Z, 0,P_O0,P_Z, R_NEAREST | R_UP, + FE_INEXACT | FE_OVERFLOW, 0,P_O,P_Z }, + {__LINE__,B_ADD, 0,P_O0,P_Z, 0,P_O0,P_Z, R_ZERO | R_DOWN, + FE_INEXACT | FE_OVERFLOW, 0,P_O0,P_O }, + {__LINE__,B_ADD, 1,P_O0,P_Z, 1,P_O0,P_Z, R_NEAREST | R_DOWN, + FE_INEXACT | FE_OVERFLOW, 1,P_O,P_Z }, + {__LINE__,B_ADD, 1,P_O0,P_Z, 1,P_O0,P_Z, R_ZERO | R_UP, + FE_INEXACT | FE_OVERFLOW, 1,P_O0,P_O }, + {__LINE__,B_ADD, 0,P_O0,P_Z, 1,P_O0,P_Z, R_ALL & ~R_DOWN, + 0, 0,P_Z,P_Z }, + {__LINE__,B_ADD, 0,P_O0,P_Z, 1,P_O0,P_Z, R_DOWN, + 0, 1,P_Z,P_Z }, + + /* Negation. */ + {__LINE__,B_NEG, 0,P_Z,P_Z, 0,0,0, R_ALL, 0, 1,P_Z,P_Z }, + {__LINE__,B_NEG, 1,P_Z,P_Z, 0,0,0, R_ALL, 0, 0,P_Z,P_Z }, + {__LINE__,B_NEG, 0,P_O,P_Z, 0,0,0, R_ALL, 0, 1,P_O,P_Z }, + {__LINE__,B_NEG, 1,P_O,P_Z, 0,0,0, R_ALL, 0, 0,P_O,P_Z }, + {__LINE__,B_NEG, 0,P_O,P_1Z, 0,0,0, R_ALL, 0, 1,P_O,P_1Z }, + {__LINE__,B_NEG, 1,P_O,P_1Z, 0,0,0, R_ALL, 0, 0,P_O,P_1Z }, + {__LINE__,B_NEG, 0,P_O,P_01Z, 0,0,0, R_ALL, 0, 1,P_O,P_01Z }, + {__LINE__,B_NEG, 1,P_O,P_01Z, 0,0,0, R_ALL, 0, 0,P_O,P_01Z }, + {__LINE__,B_NEG, 0,P_1Z,P_1Z1, 0,0,0, R_ALL, 0, 1,P_1Z,P_1Z1 }, + {__LINE__,B_NEG, 1,P_1Z,P_1Z1, 0,0,0, R_ALL, 0, 0,P_1Z,P_1Z1 }, + {__LINE__,B_NEG, 0,P_Z,P_Z1, 0,0,0, R_ALL, 0, 1,P_Z,P_Z1 }, + {__LINE__,B_NEG, 1,P_Z,P_Z1, 0,0,0, R_ALL, 0, 0,P_Z,P_Z1 }, + + /* Absolute value. */ + {__LINE__,B_ABS, 0,P_Z,P_Z, 0,0,0, R_ALL, 0, 0,P_Z,P_Z }, + {__LINE__,B_ABS, 1,P_Z,P_Z, 0,0,0, R_ALL, 0, 0,P_Z,P_Z }, + {__LINE__,B_ABS, 0,P_O,P_Z, 0,0,0, R_ALL, 0, 0,P_O,P_Z }, + {__LINE__,B_ABS, 1,P_O,P_Z, 0,0,0, R_ALL, 0, 0,P_O,P_Z }, + {__LINE__,B_ABS, 0,P_O,P_1Z, 0,0,0, R_ALL, 0, 0,P_O,P_1Z }, + {__LINE__,B_ABS, 1,P_O,P_1Z, 0,0,0, R_ALL, 0, 0,P_O,P_1Z }, + {__LINE__,B_ABS, 0,P_O,P_01Z, 0,0,0, R_ALL, 0, 0,P_O,P_01Z }, + {__LINE__,B_ABS, 1,P_O,P_01Z, 0,0,0, R_ALL, 0, 0,P_O,P_01Z }, + {__LINE__,B_ABS, 0,P_1Z,P_1Z1, 0,0,0, R_ALL, 0, 0,P_1Z,P_1Z1 }, + {__LINE__,B_ABS, 1,P_1Z,P_1Z1, 0,0,0, R_ALL, 0, 0,P_1Z,P_1Z1 }, + {__LINE__,B_ABS, 0,P_Z,P_Z1, 0,0,0, R_ALL, 0, 0,P_Z,P_Z1 }, + {__LINE__,B_ABS, 1,P_Z,P_Z1, 0,0,0, R_ALL, 0, 0,P_Z,P_Z1 }, + + /* Square root. */ + {__LINE__,B_SQRT, 0,P_Z,P_Z, 0,0,0, R_ALL, 0, 0,P_Z,P_Z }, + {__LINE__,B_SQRT, 1,P_Z,P_Z, 0,0,0, R_ALL, 0, 1,P_Z,P_Z }, + {__LINE__,B_SQRT, 0,P_O,P_1Z, 0,0,0, R_ALL, 0, 0,P_O,P_1Z }, + {__LINE__,B_SQRT, 1,P_O,P_1Z, 0,0,0, R_ALL, 0, 1,P_O,P_1Z }, + {__LINE__,B_SQRT, 0,P_O,P_01Z, 0,0,0, R_ALL, + FE_INVALID | FE_INVALID_SNAN, 0,P_O,P_11Z }, + {__LINE__,B_SQRT, 1,P_O,P_01Z, 0,0,0, R_ALL, + FE_INVALID | FE_INVALID_SNAN, 1,P_O,P_11Z }, + + {__LINE__,B_SQRT, 0,P_O,P_Z, 0,0,0, R_ALL, 0, 0,P_O,P_Z }, + {__LINE__,B_SQRT, 0,P_0O,P_Z, 0,0,0, R_ALL, 0, 0,P_0O,P_Z }, + + {__LINE__,B_SQRT, 1,P_O,P_Z, 0,0,0, R_ALL, + FE_INVALID | FE_INVALID_SQRT, 0,P_O,P_1Z }, + {__LINE__,B_SQRT, 1,P_1Z,P_1Z1, 0,0,0, R_ALL, + FE_INVALID | FE_INVALID_SQRT, 0,P_O,P_1Z }, + {__LINE__,B_SQRT, 1,P_Z,P_Z1, 0,0,0, R_ALL, + FE_INVALID | FE_INVALID_SQRT, 0,P_O,P_1Z }, + +}; + +static void +check_op(void) +{ + int i, j; + tocheck_t r, a, b, x; + int raised; + + for (i = 0; i < sizeof(optests)/sizeof(optests[0]); i++) + { + a = pattern(optests[i].a_sgn, optests[i].a_exp, + optests[i].a_mant); + b = pattern(optests[i].b_sgn, optests[i].b_exp, + optests[i].b_mant); + x = pattern(optests[i].x_sgn, optests[i].x_exp, + optests[i].x_mant); + for (j = 0; j < 4; j++) + if (optests[i].rmode & 1<<j) + { + fesetenv(rmodes+j); + switch (optests[i].op) + { + case B_ADD: r = a + b; break; + case B_SUB: r = a - b; break; + case B_MUL: r = a * b; break; + case B_DIV: r = a / b; break; + case B_NEG: r = -a; break; + case B_ABS: r = FUNC(fabs)(a); break; + case B_SQRT: r = FUNC(sqrt)(a); break; + } + raised = fetestexcept(all_exceptions); + check_result(optests[i].line,rmnames[j],x,r); + check_excepts(optests[i].line,rmnames[j], + optests[i].excepts,raised); + } + } +} + +static void +fail_xr(int line, const char *rm, tocheck_t x, tocheck_t r, tocheck_t xx, + int xflag) +{ + int i; + unsigned char *cx, *cr, *cxx; + + printf("%s:%d:round %s:fail\n with x=0x", __FILE__, line,rm); + cx = (unsigned char *)&x; + cr = (unsigned char *)&r; + cxx = (unsigned char *)&xx; + for (i = 0; i < sizeof(tocheck_t); i++) + printf("%02x", cx[i]); + printf(" r=0x"); + for (i = 0; i < sizeof(tocheck_t); i++) + printf("%02x", cr[i]); + printf(" xx=0x"); + for (i = 0; i < sizeof(tocheck_t); i++) + printf("%02x", cxx[i]); + printf(" inexact=%d\n", xflag != 0); + nerrors++; +} + +static void +check_sqrt(tocheck_t a) +{ + int j; + tocheck_t r0, r1, r2, x0, x1, x2; + int raised = 0; + int ok; + + for (j = 0; j < 4; j++) + { + int excepts; + + fesetenv(rmodes+j); + r1 = FUNC(sqrt)(a); + excepts = fetestexcept(all_exceptions); + fesetenv(FE_DFL_ENV); + raised |= excepts & ~FE_INEXACT; + x1 = r1 * r1 - a; + if (excepts & FE_INEXACT) + { + r0 = delta(r1,-1); r2 = delta(r1,1); + switch (1 << j) + { + case R_NEAREST: + x0 = r0 * r0 - a; x2 = r2 * r2 - a; + ok = fabs(x0) >= fabs(x1) && fabs(x1) <= fabs(x2); + break; + case R_ZERO: case R_DOWN: + x2 = r2 * r2 - a; + ok = x1 <= 0 && x2 >= 0; + break; + case R_UP: + x0 = r0 * r0 - a; + ok = x1 >= 0 && x0 <= 0; + break; + default: + assert(0); + } + } + else + ok = x1 == 0; + if (!ok) + fail_xr(__LINE__,rmnames[j],a,r1,x1,excepts&FE_INEXACT); + } + check_excepts(__LINE__,"all",0,raised); +} + +int main(int argc, char **argv) +{ + int i; + + _LIB_VERSION = _IEEE_; + + /* Set up environments for rounding modes. */ + fesetenv(FE_DFL_ENV); + fesetround(FE_TONEAREST); + fegetenv(rmodes+0); + fesetround(FE_TOWARDSZERO); + fegetenv(rmodes+1); + fesetround(FE_UPWARD); + fegetenv(rmodes+2); + fesetround(FE_DOWNWARD); + fegetenv(rmodes+3); + +#if defined(FE_INVALID_SOFTWARE) || defined(FE_INVALID_SQRT) + /* There's this really stupid feature of the 601... */ + fesetenv(FE_DFL_ENV); + feraiseexcept(FE_INVALID_SOFTWARE); + if (!fetestexcept(FE_INVALID_SOFTWARE)) + excepts_missing |= FE_INVALID_SOFTWARE; + fesetenv(FE_DFL_ENV); + feraiseexcept(FE_INVALID_SQRT); + if (!fetestexcept(FE_INVALID_SQRT)) + excepts_missing |= FE_INVALID_SQRT; +#endif + + check_op(); + for (i = 0; i < 100000; i++) + check_sqrt(pattern(0, P_Rno, P_R)); + for (i = 0; i < 100; i++) + check_sqrt(pattern(0, P_Z, P_R)); + check_sqrt(pattern(0,P_Z,P_Z1)); + + printf("%d errors.\n", nerrors); + return nerrors == 0 ? 0 : 1; +} diff --git a/sysdeps/powerpc/test-arithf.c b/sysdeps/powerpc/test-arithf.c new file mode 100644 index 0000000000..d78ec49009 --- /dev/null +++ b/sysdeps/powerpc/test-arithf.c @@ -0,0 +1,6 @@ +typedef float tocheck_t; +#define ESIZE 8 +#define MSIZE 23 +#define FUNC(x) x##f + +#include "test-arith.c" |