Update. cvs/libc-ud-970810

1997-08-10 18:48  Ulrich Drepper  <drepper@cygnus.com>

	* Makeconfig: Define build-programs to no if not set and $(config-os)
	is none.
	* Makerules (+depfiles): Don't add $(others) if not build-programs.
	* Rules (others): Don't depend on $(install-bin), $(install-rootbin)
	and $(install-sbin) if not build-programs.
	* Makefile (subdirs): Remove login.
	* sysdeps/mach/hurd/Subdirs: Add login.
	* sysdeps/unix/Subdirs: Add login.

	* sysdeps/generic/sys/mman.h: Test feature macros before defining
	non-POSIX things.  Add MAP_FAILED.
	* sysdeps/unix/bsd/osf/sys/mman.h: Likewise.
	* sysdeps/unix/bsd/sun/sunos4/sys/mman.h: Likewise.
	* sysdeps/unix/bsd/ultrix4/sys/mman.h: Likewise.
	* sysdeps/unix/sysv/irix4/sys/mman.h: Likewise.
	* sysdeps/unix/sysv/linux/sys/mman.h: Rewrite to not use kernel header
	but bits/mman.h.  Test feature macros before definitions non-POSIX
	functions.
	* sysdeps/unix/sysv/linux/alpha/bits/mman.h: New file.
	* sysdeps/unix/sysv/linux/i386/bits/mman.h: New file.
	* sysdeps/unix/sysv/linux/m68k/bits/mman.h: New file.
	* sysdeps/unix/sysv/linux/mips/bits/mman.h: New file.
	* sysdeps/unix/sysv/linux/powerpc/bits/mman.h: New file.
	* sysdeps/unix/sysv/linux/sparc/bits/mman.h: New file.
	* sysdeps/unix/sysv/linux/sparc64/bits/mman.h: New file.
	* sysdeps/unix/sysv/linux/alpha/Dist: Add bits/mman.h.
	* sysdeps/unix/sysv/linux/i386/Dist: Likewise.
	* sysdeps/unix/sysv/linux/m68k/Dist: Likewise.
	* sysdeps/unix/sysv/linux/mips/Dist: Likewise.
	* sysdeps/unix/sysv/linux/powerpc/Dist: Likewise.
	* sysdeps/unix/sysv/linux/sparc/Dist: Likewise.
	* sysdeps/unix/sysv/linux/sparc64/Dist: New file.  Likewise.

	* catgets/open_catalog.c (__open_catalog): Compare result of mmap
	with MAP_FAILED and not -1.
	* elf/dl-load.c (_dl_map_object_from_fd): Likewise.
	* elf/dl-minimal.c (malloc): Likewise.
	* elf/dl-misc.c (_dl_sysdep_read_whole_file): Likewise.
	* locale/loadlocale.c (_nl_load_locale): Likewise.
	* locale/programs/localedef.c: Likewise.
	* malloc/malloc.c (mmap_chunk): Likewise.
	(new_heap): Likewise.
	* string/stratcliff.c: Likewise.
	* string/tester.c: Add tests for stpcpy.

	* gmon/gmon.c (__monstartup): Use calloc instead of malloc+memset.

	* gmon/sys/gmon.h: Add prototype for write_profiling.

	* elf/dl-profile.c: Add decsriptions and implement reading of
	existing data in profiling file.
	* gmon/sys/gmon_out.h (GMON_SHOBJ_VERSION): New macro for shared
	lib profiling data.

	* sysdeps/stub/bits/libc-lock.h (__libc_once_define): Add new first
	parameter for class.

	* inet/inet_ntoa.c: Make once control variable static.
	Correctly handle case where key cannot be retrieved.
	* inet/Makefile (tests): Add tst-ntoa.
	* inet/tst-ntoa.c: New file.

	* manual/time.texi: Document localtime_r, gmtime_r, asctime_r, and
	ctime_r.
	Mention that tzname values must not be changed.

	* manual/users.texi: Document UTMP/WTMP handling functions.
	Mainly written by Mark Kettenis.

	* math/libm-test.c (cbrt_test): Add test for number with aboslute
	value < 1.0.

	* nss/nss_files/files-XXX.c (internal_setent): If opening of file
	failed with EGAIN return NSS_STATUS_TRYAGAIN.
	* nss/nss_files/files-alias.c (internal_setent): Likewise.

	* sysdeps/generic/bits/utmp.h: Pretty print.

	* sysdeps/i386/i586/strcpy.S: New file.  Much faster implementation.
	* sysdeps/i386/i586/stpcpy.S: New file.
	* sysdeps/i386/i686/Implies: Don't fall back on i586 versions since
	these are less efficient than i486/i386 versions.

	* sysdeps/libm-i387/e_rem_pio2.c: Add empty file to prevent the
	version from libm-ieee754 be compiled since the later is not needed.
	* sysdeps/libm-i387/e_rem_pio2f.c: Likewise.
	* sysdeps/libm-i387/e_rem_pio2l.c: Likewise.
	* sysdeps/libm-i387/k_rem_pio2.c: Likewise.
	* sysdeps/libm-i387/k_rem_pio2f.c: Likewise.
	* sysdeps/libm-i387/k_rem_pio2l.c: Likewise.
	* sysdeps/m68k/fpu/e_rem_pio2.c: Likewise.
	* sysdeps/m68k/fpu/e_rem_pio2f.c: Likewise.
	* sysdeps/m68k/fpu/e_rem_pio2l.c: Likewise.
	* sysdeps/m68k/fpu/k_rem_pio2.c: Likewise.
	* sysdeps/m68k/fpu/k_rem_pio2f.c: Likewise.
	* sysdeps/m68k/fpu/k_rem_pio2l.c: Likewise.

	* sysdeps/libm-i387/s_cbrt.S: Fix several bugs.
	* sysdeps/libm-i387/s_cbrtf.S: Likewise.
	* sysdeps/libm-i387/s_cbrtl.S: Likewise.

	* sysdeps/unix/sysv/linux/if_index.c (if_indextoname): Correct
	little bug.

	* sysdeps/unix/sysv/linux/bits/socket.h (struct msghdr): Make field
	msg_iovlen of type int.

	* time/tzfile.c: Change return type of __tzstring to char *.
	* time/tzset.c: Likewise.  Change definition of __tzstring.

	* time/tzset.c: Interpret no DST information in user provided time
	zone specification as it is meant: no DST.
	Patch by Paul Eggert <eggert@twinsun.com>.

1997-07-20 03:01  Geoff Keating  <geoffk@ozemail.com.au>

	* sysdeps/elf/dl-load.c (add_name_to_object): New procedure.
	(_dl_map_object_from_fd): Use add_name_to_object.
	(_dl_map_object): If an object's SONAME is used to resolve a
	dependency, add it to the list of the object's names.

	* inet/getnameinfo.c: Make `domain' non-const.

	* sysdeps/unix/sysv/linux/powerpc/kernel_termios.c: Clean up.

	* math/test-fenv.c (feenv_nomask_test): Don't execute if
	FE_NOMASK_ENV is not defined, or if fesetenv(FE_NOMASK_ENV)
	sets errno to ENOSYS.

	* sysdeps/powerpc/dl-machine.h: Print proper error message on
	unknown reloc type (rather than using assert).

	* sysdeps/unix/sysv/linux/powerpc/profil-counter.h: New file.
	* sysdeps/unix/sysv/linux/powerpc/profil.c: Deleted.

1997-07-16 12:47  Geoff Keating  <geoffk@ozemail.com.au>

	* sysdeps/powerpc/bits/fenv.h (feraiseexcept): New optimising macro.
	(feclearexcept): New optimising macro.
	(FE_NOMASK_ENV): Change to procedure so it can enable exceptions.
	(FE_ENABLED_ENV): Add.
	(FE_NONIEEE_ENV): Add.
	* sysdeps/powerpc/bits/fenv_const.c: Add __fe_enabled_env,
	__fe_nonieee_env; delete __fe_nomask_env.
	* sysdeps/powerpc/bits/fe_nomask.c: New file (stub warning until it
	gets implemented in the kernel).

	* sysdeps/powerpc/fraiseenv.c: Deal with chips that don't have
	FE_INVALID_SOFTWARE implemented. Use macros for bit names for clarity.
	* sysdeps/powerpc/fsetexcptflag.c: Likewise.

	* io/ftw.c: Don't compare holes in structures.
	* sysdeps/unix/sysv/linux/sys/sysmacros.h: Cast the result of the
	macros to `int', because otherwise it might be `long long' which
	the calling code is probably not expecting.

	* sysdeps/libm-ieee754/s_lround.c [NO_LONG_DOUBLE]: Fix a few bugs,
	document the existence of some more.

	* sysdeps/powerpc/s_llrint.c: New file.
	* sysdeps/powerpc/s_lrint.c: New file.
	* sysdeps/powerpc/s_llround.c: New file.
	* sysdeps/powerpc/s_lround.c: New file.

	* sysdeps/powerpc/s_sqrt.c: New file.
	* sysdeps/powerpc/s_sqrtf.c: New file.
	* sysdeps/powerpc/w_sqrt.s: New empty file.
	* sysdeps/powerpc/w_sqrtf.s: New empty file.
	* sysdeps/powerpc/t_sqrt.c: New file.
	* sysdeps/powerpc/test-arithf.c: New file.
	* sysdeps/powerpc/Makefile [subdir=math]: Add t_sqrt to support
	routines.  Add test-arithf to test programs.

	* sysdeps/powerpc/bits/mathdef.h: Add FP_ILOGB0, FP_ILOGBNAN.

	* sysdeps/powerpc/strcmp.s: Simplify drastically. Now much neater,
	and possibly faster (or possibly slower, depending on input).

1997-06-08 22:55  Geoff Keating  <geoffk@ozemail.com.au>

	* sysdeps/powerpc/fenvbits.h: Correct FE_DFL_ENV and FE_NOMASK_ENV
	macros.
	* sysdeps/powerpc/s_rint.c: New file.

1997-05-22 08:47  Geoff Keating  <geoffk@ozemail.com.au>

	* sysdeps/powerpc/Makefile [subdir=math]: Add q_* routines.
	* sysdeps/powerpc/Dist: Add quad_float.h.
	* sysdeps/powerpc/q_dtoq.c: New file.
	* sysdeps/powerpc/q_itoq.c: New file.
	* sysdeps/powerpc/q_lltoq.c: New file.
	* sysdeps/powerpc/q_neg.c: New file.
	* sysdeps/powerpc/q_qtoi.c: New file.
	* sysdeps/powerpc/q_qtoll.c: New file.
	* sysdeps/powerpc/q_qtos.c: New file.
	* sysdeps/powerpc/q_qtou.c: New file.
	* sysdeps/powerpc/q_qtoull.c: New file.
	* sysdeps/powerpc/q_stoq.c: New file.
	* sysdeps/powerpc/q_ulltoq.c: New file.
	* sysdeps/powerpc/q_utoq.c: New file.
	* sysdeps/powerpc/quad_float.h: New file.
	* sysdeps/powerpc/test-arith.c: New file.

	* sysdeps/powerpc/fpu_control.h: Fix _FPU_GETCW.
	* sysdeps/powerpc/fegetround.c: Use mcrfs to be faster and not
	require a stack frame.
	* sysdeps/powerpc/bits/fenv.h: Include inline macro for fegetround.

1997-05-18 05:55  Geoff Keating  <geoffk@ozemail.com.au>

	* sysdeps/powerpc/fenv_libc.h (fegetenv_register,
	fesetenv_register): Add 'volatile'.
	(set_fpscr_bit, reset_fpscr_bit): New macros, FPSCR_*
	constants to use with them.
	* sysdeps/powerpc/s_copysign.S: New file.
	* sysdeps/powerpc/s_copysignf.s: New file.
	* sysdeps/powerpc/s_fabs.S: New file.
	* sysdeps/powerpc/s_fabsf.s: New file.
	* sysdeps/powerpc/s_isnan.c: New file.
	* sysdeps/powerpc/s_isnanf.s: New file.
	* sysdeps/powerpc/s_rintf.c: New file.

	* sysdeps/powerpc/fenvbits.h: Make FE_INVALID the summary
	bit in the FPSCR, not the enable bit.
	* sysdeps/powerpc/fraiseexcpt.c: Consequent change to the above.
	* sysdeps/powerpc/fclrexcpt.c: Correct.
	* sysdeps/powerpc/fsetexcptflag.c: Correct.
	* sysdeps/powerpc/ftestexcpt.c: Is now much simpler.
	* sysdeps/powerpc/fgetexcptflg.c: Simplify.

	* sysdeps/powerpc/strlen.s: Schedule better, save 3 clocks :-).

	* sysdeps/powerpc/dl-machine.h (elf_machine_rela): Add comment
	explaining some reentrancy issues with lazy PLT entries.

1997-08-09 13:04  Mark Kettenis  <kettenis@phys.uva.nl>

	* login/logout.c (logout): utmpname returns -1 on error.

	* login/libutil.map: Remove updwtmp.

	* login/getutline.c: Rename getutline to __getutline and make
	getutline a weak alias.  Make getutxline a weak alias for
	__getutline.
	* login/getutid.c: Rename getutid to __getutid and make getutid a
	weak alias.  Make getutxid a weak alias for __getutid.
	* libc.map: Add getutxid, getutxline.

	* login/utmpname.c (__utmpname): Reset backend right after backend
	endutent call.
	* login/utmp_file.c: Reordered functions.  Remove unecessary
	header files.
	(getutent_r_file, pututline_file): Do not call setutent_file.  At this
	point the file is guaranteed to be open (assert!).
	(getutid_r_file, getutline_r_file): Remove check for validity of file
	descriptor.  At this point the file is guaranteed to be open.
	(getutent_r_file, internal_getut_r, getutline_r_file, pututline_file,
	updwtmp_file): Do not wait when unlocking file.
	* login/utmp_daemon.c: General cleanup and a bit of reordering.
	(getutent_r_daemon, pututline_daemon): Do not call setutent_daemon.
	At this point the socket is guaranteed to be open (assert!).
	(getutid_r_daemon, getutline_r_daemon): Do not check if daemon_sock is
	valid.  At this point the socket is guaranteed to be open (assert!).
	* login/getutline_r.c: Remove unnecessary header files.
	(__getutline_r): Do not call backend setutent.
	* login/getutid_r.c: Remove unnecessary header files.
	(__getutid_r): Do not call backend setutent.
	* login/getutent_r.c: Remove unneccesary header files.
	(__libc_utmp_unknown_functions): Added getutid_r_unknown,
	getutline_r_unknown.
	(setutent_unknown): Only set file backend if setutent for the file
	backend was successful.
	(getutent_r_unknown, pututline_unknown): Call setutent_unknown instead
	of __setutent.  Report failure if no backend was selected.
	(getutid_r_unknown): New function.
	(getutline_r_unknown): New function.
	(__endutent): Reset backend.  This makes sure all backends are checked
	on the next setutent call.

1997-08-08 20:20  Thorsten Kukuk  <kukuk@vt.uni-paderborn.de>

	* nis_cache.c: Replace dummy functions.

	* libc.map: Add xdr_sizeof symbol.

	* sunrpc/Makefile: Add xdr_sizeof to routines.
	* sunrpc/rpc/xdr.h: Add xdr_sizeof prototype.
	* sunrpc/xdr_sizeof.c: New, from tirpc 2.3.

1997-08-08  Thomas Bushnell, n/BSG  <thomas@gnu.ai.mit.edu>

	* sysdeps/mach/bits/libc-lock.h (__libc_once): Define correctly.

1997-08-07  Thomas Bushnell, n/BSG  <thomas@gnu.ai.mit.edu>

	* sysdeps/mach/hurd/profil.c (fetch_samples): Put parens in the
	right place.

1997-08-06  Thomas Bushnell, n/BSG  <thomas@gnu.ai.mit.edu>

	* sysdeps/mach/hurd/profil.c (fetch_samples): Do arithmetic on
	PC's in long long to avoid overflow.

1997-08-07  Thomas Bushnell, n/BSG  <thomas@gnu.ai.mit.edu>

	* sysdeps/mach/bits/libc-lock.h (__libc_once, __libc_once_define):
	New macros.

1997-08-06  Andreas Jaeger  <aj@arthur.rhein-neckar.de>

	* nis/Makefile (headers): Remove bits/nislib.h.

1997-08-06 14:54  Ulrich Drepper  <drepper@cygnus.com>

	* sysdeps/mach/hurd/Subdirs: Add login.

1997-08-06 14:23  Klaus Espenlaub  <kespenla@hydra.informatik.uni-ulm.de>

	* db/hash/hash.c (init_hash): Don't use stat() if it doesn't provide
	the preferred block size.

	* login/programs/database.c (store_state_entry): Don't compile if
	there is no ut_type field.
	(store_state_entry, store_process_entry): Use the ut_tv field
	for timestamps if supported.

	* login/programs/utmpdump.c (print_entry): Always use ut_tv field.

	* login/programs/xtmp.c: Fix numerous xtmp/utmp typos.  Use the ut_tv
	field for timestamps if supported.

	* login/programs/xtmp.h: Fix xtmp/utmp typo.

	* sysdeps/posix/defs.c (stdstream): Change (PTR) to (void *).

	* sysdeps/stub/connect.c (connect): Change to __connect, add alias.
	* sysdeps/stub/send.c (send): Likewise.
	* sysdeps/stub/s_exp2f.c: Emit correct stub_warning().
	* sysdeps/stub/statfs.c: Move stub_warning() to the usual place.
	* sysdeps/stub/init-first.c: Add definition of __libc_pid.

1997-08-05 13:28  Philip Blundell  <pb@spring.nexus.co.uk>

	* sysdeps/standalone/arm/bits/errno.h: Add EISDIR, EOPNOTSUPP;
	tidy up formatting.

	* Makefile (subdirs): Remove `login'.
	* sysdeps/unix/Subdirs: New file; build `login' subdirectory
	for Unix systems.

1997-08-05  Andreas Schwab  <schwab@issan.informatik.uni-dortmund.de>

	* sysdeps/generic/bits/utmpx.h: New file.

	* sysdeps/mach/hurd/Dist: Add some files.
	* sysdeps/mips/Dist: Likewise.
	* sysdeps/mips/mips64/Dist: Likewise.
	* sysdeps/sparc/Dist: Likewise.
	* sysdeps/unix/sysv/linux/mips/Dist: Likewise.
	* sysdeps/unix/sysv/linux/sparc/Dist: Likewise.
	* sysdeps/mips/mipsel/Dist: New file.
	* sysdeps/sparc64/elf/Dist: New file.
	* sysdeps/unix/sysv/linux/sparc64/Dist: New file.

1997-08-05  Andreas Schwab  <schwab@issan.informatik.uni-dortmund.de>

	* libc.map: Add missing symbols.

1997-08-05  Andreas Jaeger  <aj@arthur.rhein-neckar.de>

	* manual/socket.texi: Correct typos.

	* manual/math.texi: Correct typos.

	* manual/time.texi (Formatting Date and Time): Likewise.

1997-08-04 13:06  Thomas Bushnell, n/BSG  <thomas@gnu.ai.mit.edu>

	* gmon/gmon.c (write_gmon): New function; guts from _mcleanup.
	(_mcleanup): Use write_gmon.
	(write_profiling): This function can be safely called at any time
	to write a current histogram without interfering with ongoing
	profiling.

	* sysdeps/mach/hurd/profil.c (fetch_samples): Initialize NSAMPLES.

1997-08-01 17:53  Thomas Bushnell, n/BSG  <thomas@gnu.ai.mit.edu>

	* sysdeps/mach/hurd/profil.c (fetch_samples): Sample buffer need
	not be vm_deallocated; it's a stack buffer.
	(profil_vm_deallocate): Delete prototype.
	(#include <../mach/RPC_vm_deallocate_rpc.c>): Drop this
	inclusion.
	* sysdeps/mach/hurd/Makefile
	($(common-objpfx)hurd/../mach/RPC_vm_deallocate_rpc.c): Delete
	this rule.

	* sysdeps/mach/hurd/profil.c (fetch_samples): New function, guts
	from profile_waiter.
	(profile_waiter): Use fetch_samples.
	(profil): When turning off profiling, fetch the last bunch of
	samples before disabling PC sampling.
	(fetch_samples): Add prototype.

1997-07-30 12:53  Thomas Bushnell, n/BSG  <thomas@gnu.ai.mit.edu>

	* sysdeps/mach/hurd/Makefile: Give normal rules for the mach RPC
	source included by profil.c instead of trying to use
	before-compile grot.

1997-07-23 15:04  Thomas Bushnell, n/BSG  <thomas@gnu.ai.mit.edu>

	* sysdeps/mach/hurd/profil.c (profile_waiter): Do vm_deallocate
	after releasing lock.  Use special RPC stubs and avoid assert,
	thus making this function entirely independent of the threadvar
	mechanism and of cthreads.
	(lock): Convert to be a spin lock instead of a mutex, so that the
	waiter thread need not use cthreads.  Use a fork prepare hook
	instead of _hurd_fork_locks to make sure we are clean before
	forking.
	(fork_profil_child): Renamed from fork_profil.
	(profil_reply_port): New variable.
	(profil_vm_deallocate, profil_task_get_sampled_pcs): Special RPC
	stubs made by special hacks.
	* sysdeps/mach/hurd/Makefile (before-compile): Add the mach RPC
	source files that profil.c needs if we are in the gmon directory.

	* mach/setup-thread.c (__mach_setup_thread): Delete avoidance of a
	cthreads bug that hasn't existed for two years.

	* stdio-common/printf_fp.c (__printf_fp): Correct rounding of numbers
	* locale/programs/locale.c (print_escaped): New function.
	* time/australasia: Likewise.
	(open_database, synchronize_database, initialize_database):
	* nss/nss_files/files-netgrp.c (_nss_netgroup_parseline): Don't
	* Makerules (load-map-file): Currectly handle missing map file.
	cannot rely on M_SQRT2 being defined.
	(log1p): Use __M_SQRT2 not M_SQRT2.

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"