RFC Power PC G3 optimized sqrtf function.
Conn Clark
clark@esteem.com
Thu Dec 14 21:19:00 GMT 2006
Hi everybody,
This is my 1st post and attempt at contributing to glibc
I have written a sqrtf function that is much faster on a PowerPC G3 than
the original one used. It uses the frsqrte instruction and Newton
Goldschmidt iterations to get a result. I have had testing done on G3 and G4
processors and it has results that conform to ieee . Testing on G5's has
shown it gets within at least 2 bits of the correct answer. This shouldn't
matter because the G5 has a hardware sqrtf instruction. It may work on 603,
603e, 604, and 604e processors as well but I have not tested on them. It
will not work on a 601 processor.
The limiting factor on ieee conformance is the frsqrte instruction must
produce a result that is within 1/59th of the correct value. A timing test
on all valid values using the current glibc function takes about 26 minutes
on a iMac g3 400MHz machine. With my implementation it takes about 21
minutes.
Please read the header for more details and give me some feedback.
P.S. Do I need to file copyright assignment papers for this?
Thanks,
Conn
---------------------------------------
Conn Clark
Electronic Systems Technology
415 N. Quay Street Building B1 (509)-735-9092 ext 117
Kennewick, WA. 99336
Gentoo Linux RU13$!!!
-------------- next part --------------
/* Single-precision floating point square root Optimized for G3 Powerpc32
processors written By Conn Clark.
Copyright (C) 2004, 2006 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 Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 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
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, write to the Free
Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston MA
02110-1301 USA.
Notes:
This is sqrtf function utilizes the frsqrte instruction
which is present on the 603, 603e, 604, 604e, G3, G4 and G5
processors. It generates ieee compliant results as long as
the frsqrte returns an estimate that is within 1/59th the
correct number. As far as I know the G5 is the only
processor that has less accuracy than this. The 603, 603e,
604, and 604e may have less accuracy but I don't have a
system to test on. This function won't work on a 601 or
other processor that doesn't support the frsqrte
instruction. If the result of the frsqrte instruction is
less than 1/59th the correct value but with in 1/32 this
function will return a value that is at most 1/4.22e+06 off.
Algorithm Notes:
This takes the reciprocal square root estimate and feeds
it into three Newton Goldschmidt iterations. The result is
a reciprocal square root. This value is then multiplied by
the input value to get the square root.
Care has been taken to avoid loading constants from the
data segment of Glibc. This was done to avoid a data cache
load.
Optimization emphasis has been on dealing with valid
inputs ( 0 <= x < +infinity). Dealing with invalid inputs
could be done faster or better.
*/
.section ".text"
.align 2
.globl __slow_ieee754_sqrtf
.type __slow_ieee754_sqrtf, @function
__slow_ieee754_sqrtf:
b __ieee754_sqrtf
.size __slow_ieee754_sqrtf, .-__slow_ieee754_sqrtf
.section ".text"
.align 2
.globl __ieee754_sqrtf
.type __ieee754_sqrtf, @function
__ieee754_sqrtf:
/* start loading some constants for integer comparison */
lis 3,0x3f00 /* 0.5F equiv as an integer */
lis 4,0x3FC0 /* 1.5F equiv as an integer */
stfsu 1,-12(1) /* store original and get 12 bytes of stack space. */
stw 3,4(1) /* store 0.5F on the stack to be loaded by fpu */
mffs 6 /* store fpu configuration */
stw 4,8(1) /* store 1.5F on the stack to be loaded by fpu */
lis 7, 0x7F80 /* load NAN for testing */
lfs 2,4(1) /* load 0.5F into fpu reg 2 */
/*mtfsfi 7,0 */ /* relax fpu (not needed) */
lwz 5,0(1) /* load original value as an int for testing */
fmr 9,1 /* copy original value into fpu reg 9 */
cmpi 0,5,0x0000 /* test for positive zero */
rlwinm 12,5,0,0,1 /* mask off sign bit and store in reg 12*/
lfs 7,8(1) /* load 1.5F into fpu reg 7 */
cmpl 1,12,7 /* test for NAN results in cr1 */
ble neg_number_or_zero /* branch if less than or equal to zero */
frsqrte 1,1 /* get recip sqrt estimate (does no harm if input is NaN */
beq cr1, not_a_number /* branch if original value was not a number */
fmul 2,2,9 /* begin Goldschmidt */
fmuls 3,1,1 /* saves one clock without effecting accuracy */
fmul 4,2,3
fnmsubs 3,2,3,7 /* saves one clock without effecting accuracy */
fmul 5,3,3
fmul 1,3,1
fnmsubs 3,4,5,7 /* saves one clock without effecting accuracy */
fmul 4,4,5
fmul 1,3,1
fmul 5,3,3
fnmsub 3,4,5,7
fmul 1,3,1
fmul 1,9,1 /* get sqrt from recip sqrt */
addi 1,1,12 /* clean up stack */
mtfsf 0xff,6 /* restore fpu state */
frsp 1,1 /* round result to a float */
blr /* return */
neg_number_or_zero:
lis 9,0x3F80 /* load up equiv of 1.0F */
stw 9,8(1) /* store 1.0F to where fpu can load it */
lis 6, 0x8000 /* negative */
lfs 2,8(1) /* load 1.0F into fpu reg 2 */
beq its_zero /* branch if zero */
cmpl 0,5,6 /* test for negative zero */
be its_zero /* branch if zero */
stfs 6,0(1) /* store fpu status */
lis 3,0x7FC0 /* load aNaN */
lis 4,0x2000 /* load FE_INVALID flag */
lwz 5,0(1) /* load fpu status in gp register */
ori 4,4,0x0200 /* load INV_SQRT flag */
or 5,5,4 /* or FE_INVALID and INV_SQRT flags with fpu status */
stw 3,4(1) /* store aNaN on stack */
stw 5,0(1) /* store fpu status */
lfs 6,0(1) /* load new fpu status */
lfs 1,4(1) /* load aNaN to be returned */
mtfsf 0xff,6 /* update fpu status to new value */
its_zero:
fmuls 1,1,2 /* multiply by 1.0 to set appropriate status bits */
addi 1,1,12 /* clean up stack */
blr /* return */
not_a_number:
lis 6,0x3F80 /* load up equiv of 1.0F */
stw 6,8(1) /* store 1.0F to where fpu can load it */
lfs 2,8(1) /* load 1.0F into fpu reg 2 */
mtfsf 0xff,6 /* restore fpu status */
fmuls 1,1,2 /* multiply by 1.0 to set appropriate status bits */
addi 1,1,12 /* clean up stack */
blr /* return */
.size __ieee754_sqrtf, .-__ieee754_sqrtf
More information about the Libc-alpha
mailing list