Michael S wrote:
> But one thing is sure - [on Intel] my quad-precision sqrt() is close
> to 15 times faster than the one supplied with GCC Quad-Precision Math
> Library.

That seems almost impossible, how is it even possible to generate
something that slow?

Assuming a naive/classic SQRT NR iteration with a division might explain
it, but it seems like a stretch even so.

If your version is approximate (0.5x ulp but not exactly rounded for all
inputs?), then that might explain a factor of 2, but not 15!

Terje

--
- <Terje.Mathisen at tmsw.no>
"almost all programming can be viewed as an exercise in caching"

On Friday, March 25, 2022 at 12:36:10 PM UTC-5, Terje Mathisen wrote:
> Michael S wrote:
> > But one thing is sure - [on Intel] my quad-precision sqrt() is close
> > to 15 times faster than the one supplied with GCC Quad-Precision Math
> > Library.
> That seems almost impossible, how is it even possible to generate
> something that slow?
>
> Assuming a naive/classic SQRT NR iteration with a division might explain
> it, but it seems like a stretch even so.
<
Convert quad to double:: dozen cycles :: zero cycles in HW
SQRT(double() :: 2 dozen cycles .......... :: 20-24 cycles in HW
Convert double to quad:: ½ dozen cycles :: zero cycles in HW
2 Newton-Raphson iterations (4 quad FMAC delays) :: ½ in HW N-R*
<
(*) in HW one can perform Newton-Raphson iteration without having to
round/normalize except for the final result delivery.
<
>
> If your version is approximate (0.5x ulp but not exactly rounded for all
> inputs?), then that might explain a factor of 2, but not 15!
<
I only see 4×-5× slower.
<
> Terje
>
> --
> - <Terje.Mathisen at tmsw.no>
> "almost all programming can be viewed as an exercise in caching"

MitchAlsup wrote:
> On Friday, March 25, 2022 at 12:36:10 PM UTC-5, Terje Mathisen wrote:
>> Michael S wrote:
>>> But one thing is sure - [on Intel] my quad-precision sqrt() is close
>>> to 15 times faster than the one supplied with GCC Quad-Precision Math
>>> Library.
>> That seems almost impossible, how is it even possible to generate
>> something that slow?
>>
>> Assuming a naive/classic SQRT NR iteration with a division might explain
>> it, but it seems like a stretch even so.
> <
> Convert quad to double:: dozen cycles :: zero cycles in HW
> SQRT(double() :: 2 dozen cycles .......... :: 20-24 cycles in HW
> Convert double to quad:: ½ dozen cycles :: zero cycles in HW
> 2 Newton-Raphson iterations (4 quad FMAC delays) :: ½ in HW N-R*
> <
> (*) in HW one can perform Newton-Raphson iteration without having to
> round/normalize except for the final result delivery.
> <
>>
>> If your version is approximate (0.5x ulp but not exactly rounded for all
>> inputs?), then that might explain a factor of 2, but not 15!
> <
> I only see 4×-5× slower.

I agree.

What about a classic bit/iteration (or CORDIC?), that could eeasily take
700-1000 cycles?

Terje

--
- <Terje.Mathisen at tmsw.no>
"almost all programming can be viewed as an exercise in caching"

On Friday, March 25, 2022 at 4:37:32 PM UTC-5, Terje Mathisen wrote:
> MitchAlsup wrote:
> > On Friday, March 25, 2022 at 12:36:10 PM UTC-5, Terje Mathisen wrote:
> >> Michael S wrote:
> >>> But one thing is sure - [on Intel] my quad-precision sqrt() is close
> >>> to 15 times faster than the one supplied with GCC Quad-Precision Math
> >>> Library.
> >> That seems almost impossible, how is it even possible to generate
> >> something that slow?
> >>
> >> Assuming a naive/classic SQRT NR iteration with a division might explain
> >> it, but it seems like a stretch even so.
> > <
> > Convert quad to double:: dozen cycles :: zero cycles in HW
> > SQRT(double() :: 2 dozen cycles .......... :: 20-24 cycles in HW
> > Convert double to quad:: ½ dozen cycles :: zero cycles in HW
> > 2 Newton-Raphson iterations (4 quad FMAC delays) :: ½ in HW N-R*
> > <
> > (*) in HW one can perform Newton-Raphson iteration without having to
> > round/normalize except for the final result delivery.
> > <
> >>
> >> If your version is approximate (0.5x ulp but not exactly rounded for all
> >> inputs?), then that might explain a factor of 2, but not 15!
> > <
> > I only see 4×-5× slower.
> I agree.
>
> What about a classic bit/iteration (or CORDIC?), that could eeasily take
> 700-1000 cycles?
<
Never more than bit-width(fraction)×2+5 (RCP, RSQRT)
<
The CORDIC you remember took 250-ish cycles to do (effecitvely) Payne-
Hanek argument reduction, followed by 2×fraction CORDIC--giving 300-400
SIN, COS,... Give an argument of 2.7E+709 to COS and you can easily take
700-1400 cycles in argument reduction, followed by a rather brisk 100 cycles
for the actual transcendental calculation. That is basically why we don't do
CORDIC anymore.
<
> Terje
>
>
> --
> - <Terje.Mathisen at tmsw.no>
> "almost all programming can be viewed as an exercise in caching"

Michael S <already5chosen@yahoo.com> schrieb:

> For reference, I compiled your program (MSYS2, gfortran 11.2.0, -O2) and run it on all tree systems.
> Core i5-3450: 3.75962 seconds
> Xeon E3-1271 v3: 2.96402 seconds
> Xeon E-2176G: 2.70312 seconds
>
> So, it seems, either your Zen1 CPU runs at much higher frequency (over 4.5 GHz?) or your version of quadmath library
> is much better than the one, supplied with MSYS2 (mingw-w64-x86_64-gcc-libgfortran 11.2.0-10) or the library
> likes AMD CPUs and hates Intel's.

Or something is going away in the calling sequence of the quadmath
library.

If you have Linux installed, a timing comparison would be of interest.

> Another, not very probable possibility is that your cpu_time() call is lying, but that's easily verifiable with running
> binary under 'time' utility.
> And yet another possibility is that your compiler managed to parallelize a line 'a = sqrt(a)'. Crazy suggestion, I know.
>
> But one thing is sure - [on Intel] my quad-precision sqrt()
> is close to 15 times faster than the one supplied with GCC
> Quad-Precision Math Library.

Have you considered submitting a patch?

MitchAlsup wrote:
> On Friday, March 25, 2022 at 4:37:32 PM UTC-5, Terje Mathisen wrote:
>> MitchAlsup wrote:
>>> On Friday, March 25, 2022 at 12:36:10 PM UTC-5, Terje Mathisen wrote:
>>>> Michael S wrote:
>>>>> But one thing is sure - [on Intel] my quad-precision sqrt() is close
>>>>> to 15 times faster than the one supplied with GCC Quad-Precision Math
>>>>> Library.
>>>> That seems almost impossible, how is it even possible to generate
>>>> something that slow?
>>>>
>>>> Assuming a naive/classic SQRT NR iteration with a division might explain
>>>> it, but it seems like a stretch even so.
>>> <
>>> Convert quad to double:: dozen cycles :: zero cycles in HW
>>> SQRT(double() :: 2 dozen cycles .......... :: 20-24 cycles in HW
>>> Convert double to quad:: ½ dozen cycles :: zero cycles in HW
>>> 2 Newton-Raphson iterations (4 quad FMAC delays) :: ½ in HW N-R*
>>> <
>>> (*) in HW one can perform Newton-Raphson iteration without having to
>>> round/normalize except for the final result delivery.
>>> <
>>>>
>>>> If your version is approximate (0.5x ulp but not exactly rounded for all
>>>> inputs?), then that might explain a factor of 2, but not 15!
>>> <
>>> I only see 4×-5× slower.
>> I agree.
>>
>> What about a classic bit/iteration (or CORDIC?), that could eeasily take
>> 700-1000 cycles?
> <
> Never more than bit-width(fraction)×2+5 (RCP, RSQRT)
> <
> The CORDIC you remember took 250-ish cycles to do (effecitvely) Payne-
> Hanek argument reduction, followed by 2×fraction CORDIC--giving 300-400
> SIN, COS,... Give an argument of 2.7E+709 to COS and you can easily take
> 700-1400 cycles in argument reduction, followed by a rather brisk 100 cycles
> for the actual transcendental calculation. That is basically why we don't do
> CORDIC anymore.

This was purely in a QSQRT() context, i.e. no need for any argument
reduction, just grab the bottom exponent bit and the mantissa, pick a
starting estimate from the top bits and then start churning away,
conditionally adding bits as we go along while maintaining the current
sqrt estimate and its square, and the residue?

Going from x = y*y to x' = (y+1)*(y+1) = x + 2y + 1 which is very easy
in fixed-point, and when we end up we have done one extra (guard) bit,
at which point the residue is the sticky term, making any rounding mode
easy.

I haven't tried to actually write this algorithm, but it seems like it
could be done with LEA (to generate the 2y+1 term), a subtract and a
CMOV to pick either the new smaller residue (adding a '1' bit at the
end) or skipping it. Probably smart to do so branchless since the bits
will be unpredictable in the general case.

If we start with the HW DSQRT, do a little bit of setup work and then
run the algorithm for the next 60 iterations then we need to work with
256-bit variables for the residue and the squared bottom bit, so each
add/sub/cmp will take ~4 cycles, then the CMOVs add 3-4 more, but
definitely doable in 10 cycles/iteration which would land us at slightly
below my 700 cycle estimate.

Terje

Terje

--
- <Terje.Mathisen at tmsw.no>
"almost all programming can be viewed as an exercise in caching"

On Friday, March 25, 2022 at 8:36:10 PM UTC+3, Terje Mathisen wrote:
> Michael S wrote:
> > But one thing is sure - [on Intel] my quad-precision sqrt() is close
> > to 15 times faster than the one supplied with GCC Quad-Precision Math
> > Library.
> That seems almost impossible, how is it even possible to generate
> something that slow?
>
> Assuming a naive/classic SQRT NR iteration with a division might explain
> it, but it seems like a stretch even so.
>
> If your version is approximate (0.5x ulp but not exactly rounded for all
> inputs?), then that might explain a factor of 2, but not 15!

As far as I remember, my version is exact.
Indeed, it was 5 years ago, but looking at the code it seems that I calculate
a low estimate of result with few bits of extra precision then if estimate
is slightly lower than mid-point between two representable numbers then
and only then I calculate a square of the mid-point which tells me where
to round. So, for a common case it costs me only one or two cycles and
for the "bad" case it costs branch mispredict plus, may be 10-15 cycles.

> Terje
>
> --
> - <Terje.Mathisen at tmsw.no>
> "almost all programming can be viewed as an exercise in caching"

On Saturday, March 26, 2022 at 12:45:45 PM UTC+3, Thomas Koenig wrote:
> Michael S <already...@yahoo.com> schrieb:
> > For reference, I compiled your program (MSYS2, gfortran 11.2.0, -O2) and run it on all tree systems.
> > Core i5-3450: 3.75962 seconds
> > Xeon E3-1271 v3: 2.96402 seconds
> > Xeon E-2176G: 2.70312 seconds
> >
> > So, it seems, either your Zen1 CPU runs at much higher frequency (over 4.5 GHz?) or your version of quadmath library
> > is much better than the one, supplied with MSYS2 (mingw-w64-x86_64-gcc-libgfortran 11.2.0-10) or the library
> > likes AMD CPUs and hates Intel's.
> Or something is going away in the calling sequence of the quadmath
> library.
>
> If you have Linux installed, a timing comparison would be of interest.

No, not at home.
Even at work, I am not sure that I have Linux with Fortran compiler installed.

Could you run you program under 'time' ?
At very least it would prove that no autopar involved.

> > Another, not very probable possibility is that your cpu_time() call is lying, but that's easily verifiable with running
> > binary under 'time' utility.
> > And yet another possibility is that your compiler managed to parallelize a line 'a = sqrt(a)'. Crazy suggestion, I know.
> >
> > But one thing is sure - [on Intel] my quad-precision sqrt()
> > is close to 15 times faster than the one supplied with GCC
> > Quad-Precision Math Library.
> Have you considered submitting a patch?

My quad-precision format is different from IEEE binary128==GNU __float128.
So, it would take some work.
And since I am not a great fun of GPL, I'd rather prefer my own work to be given away under
more anarchistic license.

But I don't mind if you or your friends take my source and adapte it to quadmath.
Can put a needed license statement in my public github repository if somebody tells
me a relevant lawyerspeak.

On 2022-03-23, Terje Mathisen wrote:
> Marcus wrote:
>> On 2022-03-22, Terje Mathisen wrote:
>>> Marcus wrote:
>>>> What are your feelings towards including such instructions in an ISA?
>>>
>>> Stupid not to do it?
>>>
>>>>
>>>> My own feelings are mixed.
>>>>
>>>> Pros:
>>>> * Easy to implement in hardware.
>>>> * Can provide a significant speedup for certain workloads, especially
>>>>    if limited accuracy is acceptable.
>>>
>>> Even when you need exact results, having that reciprocal starting
>>> point means that you will converge on that value significantly faster.
>>>
>>> I.e. starting with a ~12-bit approximation gives you float with
>>> fractional ulp accuracy after one NR stage and exact float with one
>>> more iteration, still significantly faster than an SRT divider/sqrt
>>> circuit.
>>
>> Follow-up question 1: What approximation accuracy should one aim for?
>> I've seen 8 bits (yielding SP precision in 2 NR steps), and you mention
>> 12 bits. With 13 bits, would you not get full SP precision in 1 NR step?
>> CRAY gave a full 30 bits in the first approximation (!).
>>
>> Follow-up question 2: Implementation-wise, is it better to do a raw
>> nearest-neighbor ROM lookup, or more sophisticated ROM + linear
>> interpolation? It feels like the latter would give you more bang for
>> the buck, possibly at the cost of added latency.
>
> If you aim for single-cycle SIMD lookup, then I really can't see any
> room for any form of interpolation, while accessing the same ROM table 4
> times is doable, right?

I was initially thinking about having multiple LUTs (one per
FU/"lane"), but my current design is probably simpler since it only
supports the following levels of parallelism:

- 1 x 32 bit FP (scalar), 23 fractional bits
- 2 x 16 bit FP (SIMD), 10 fractional bits
- 4 x 8 bit FP (SIMD), 3 fractional bits

So, I'd probably end up using something like 2 8-bit LUTs and 2 3-bit
LUTs, or something similar. Alternatively, if better accuracy is
required for the initial estimation, 1 12-bit LUT, 1 10-bit LUT and
2 3-bit LUTs. For a wider implementation these figures would have to
increase. But I suspect that it's highly dependent on the implementation
technology. E.g. FPGA vs ASIC would differ I suppose.

The interpolation I was thinking about would be something like a 4-bit
multiplier + an integer adder, just to buy a few extra bits without
having to expand the LUTs too much. The slope factor would be stored as
a few extra bits in the LUT.

>
> A 12 (or 13?) bit lookup giving 12-bit results, along with the required
> exp trickery looks like about 6-7 Kbits of table space.
>
> Terje
>

On Wednesday, March 23, 2022 at 4:26:30 AM UTC-5, Terje Mathisen wrote:
> If you aim for single-cycle SIMD lookup, then I really can't see any
> room for any form of interpolation, while accessing the same ROM table 4
> times is doable, right?
>
> A 12 (or 13?) bit lookup giving 12-bit results, along with the required
> exp trickery looks like about 6-7 Kbits of table space.
<
This table is too big for single cycle lookup except for custom designs.
And besides, 12-bit index to 12-bit data followed by a multiply (of some
sort) only gives you 9.5-bits of accurate result.
<
The other thing to note is to get N-R or Goldschmidt to converge with negative
error under the desired result, wastes 1.5 bits of output from the table. This
makes rounding correctly possible.
<
> Terje
>
> --
> - <Terje.Mathisen at tmsw.no>
> "almost all programming can be viewed as an exercise in caching"

Michael S <already5chosen@yahoo.com> schrieb:
> On Saturday, March 26, 2022 at 12:45:45 PM UTC+3, Thomas Koenig wrote:
>> Michael S <already...@yahoo.com> schrieb:
>> > For reference, I compiled your program (MSYS2, gfortran 11.2.0, -O2) and run it on all tree systems.
>> > Core i5-3450: 3.75962 seconds
>> > Xeon E3-1271 v3: 2.96402 seconds
>> > Xeon E-2176G: 2.70312 seconds
>> >
>> > So, it seems, either your Zen1 CPU runs at much higher frequency (over 4.5 GHz?) or your version of quadmath library
>> > is much better than the one, supplied with MSYS2 (mingw-w64-x86_64-gcc-libgfortran 11.2.0-10) or the library
>> > likes AMD CPUs and hates Intel's.
>> Or something is going away in the calling sequence of the quadmath
>> library.
>>
>> If you have Linux installed, a timing comparison would be of interest.
>
> No, not at home.
> Even at work, I am not sure that I have Linux with Fortran compiler installed.
>
> Could you run you program under 'time' ?
> At very least it would prove that no autopar involved.

As a maintainer of gfortran, I can assure you that there is no
autoparallelization :-).

However, here it is:

$ cat a.f90
program main
implicit none
integer, parameter :: qp = selected_real_kind(30)
integer, parameter :: n = 10**6, m = 10
real :: t1, t2
character(len=20) :: c
integer :: i,k
real(kind=qp), dimension(:), allocatable :: a
allocate (a(n))
c = '10'
call random_number(a)
a = a * 1e50_qp
call cpu_time(t1)
do i=1,m
a = sqrt(a)
read (unit=c,fmt=*) k
write (*,*) a(k)
end do
call cpu_time(t2)
write (*,'(A,F12.5,A,ES15.5,A)') "Used ",t2-t1, " seconds to calculate ", &
real(n*m,qp), " square roots."
end program main
$ gfortran -O3 a.f90
$ time ./a.out
7308311667482999228141838.93337335438
2703388922719.59256619755162850493042
1644198.56547790254000409031824041555
1282.26306406988988805467367099802872
35.8087009547943513227570444215179866
5.98403717859392575839410139231541218
2.44622917540322165097869738482238220
1.56404257467730770560212117881721192
1.25061687765570623839109221766572087
1.11830983079632553071493517444861311
Used 2.35760 seconds to calculate 1.00000E+07 square roots.

real 0m2,475s
user 0m2,445s
sys 0m0,008s
$ head /proc/cpuinfo
processor : 0
vendor_id : AuthenticAMD
cpu family : 23
model : 1
model name : AMD Ryzen 7 1700X Eight-Core Processor
stepping : 1
microcode : 0x8001129
cpu MHz : 2056.406
cache size : 512 KB
physical id : 0

>
>> > Another, not very probable possibility is that your cpu_time() call is lying, but that's easily verifiable with running
>> > binary under 'time' utility.
>> > And yet another possibility is that your compiler managed to parallelize a line 'a = sqrt(a)'. Crazy suggestion, I know.
>> >
>> > But one thing is sure - [on Intel] my quad-precision sqrt()
>> > is close to 15 times faster than the one supplied with GCC
>> > Quad-Precision Math Library.
>> Have you considered submitting a patch?
>
> My quad-precision format is different from IEEE binary128==GNU __float128.
> So, it would take some work.
> And since I am not a great fun of GPL, I'd rather prefer my own work to be given away under
> more anarchistic license.

Unfortunately, this is not what I (usually) do - I am mainly
concerned with the gfortran front end.

> But I don't mind if you or your friends take my source and adapte it to quadmath.
> Can put a needed license statement in my public github repository if somebody tells
> me a relevant lawyerspeak.

I can at least submit a PR. If there is a description of the
algorithm in the source, maybe somebody can take this up and run
with it.

On Saturday, March 26, 2022 at 10:53:09 PM UTC+3, Thomas Koenig wrote:
> Michael S <already...@yahoo.com> schrieb:
> > On Saturday, March 26, 2022 at 12:45:45 PM UTC+3, Thomas Koenig wrote:
> >> Michael S <already...@yahoo.com> schrieb:
> >> > For reference, I compiled your program (MSYS2, gfortran 11.2.0, -O2) and run it on all tree systems.
> >> > Core i5-3450: 3.75962 seconds
> >> > Xeon E3-1271 v3: 2.96402 seconds
> >> > Xeon E-2176G: 2.70312 seconds
> >> >
> >> > So, it seems, either your Zen1 CPU runs at much higher frequency (over 4.5 GHz?) or your version of quadmath library
> >> > is much better than the one, supplied with MSYS2 (mingw-w64-x86_64-gcc-libgfortran 11.2.0-10) or the library
> >> > likes AMD CPUs and hates Intel's.
> >> Or something is going away in the calling sequence of the quadmath
> >> library.
> >>
> >> If you have Linux installed, a timing comparison would be of interest.
> >
> > No, not at home.
> > Even at work, I am not sure that I have Linux with Fortran compiler installed.
> >
> > Could you run you program under 'time' ?
> > At very least it would prove that no autopar involved.
> As a maintainer of gfortran, I can assure you that there is no
> autoparallelization :-).
>
> However, here it is:
>
> $ cat a.f90
> program main
> implicit none
> integer, parameter :: qp = selected_real_kind(30)
> integer, parameter :: n = 10**6, m = 10
> real :: t1, t2
> character(len=20) :: c
> integer :: i,k
> real(kind=qp), dimension(:), allocatable :: a
> allocate (a(n))
> c = '10'
> call random_number(a)
> a = a * 1e50_qp
> call cpu_time(t1)
> do i=1,m
> a = sqrt(a)
> read (unit=c,fmt=*) k
> write (*,*) a(k)
> end do
> call cpu_time(t2)
> write (*,'(A,F12.5,A,ES15.5,A)') "Used ",t2-t1, " seconds to calculate ", &
> real(n*m,qp), " square roots."
> end program main
> $ gfortran -O3 a.f90
> $ time ./a.out
> 7308311667482999228141838.93337335438
> 2703388922719.59256619755162850493042
> 1644198.56547790254000409031824041555
> 1282.26306406988988805467367099802872
> 35.8087009547943513227570444215179866
> 5.98403717859392575839410139231541218
> 2.44622917540322165097869738482238220
> 1.56404257467730770560212117881721192
> 1.25061687765570623839109221766572087
> 1.11830983079632553071493517444861311
> Used 2.35760 seconds to calculate 1.00000E+07 square roots.
>
> real 0m2,475s
> user 0m2,445s
> sys 0m0,008s
> $ head /proc/cpuinfo
> processor : 0
> vendor_id : AuthenticAMD
> cpu family : 23
> model : 1
> model name : AMD Ryzen 7 1700X Eight-Core Processor
> stepping : 1
> microcode : 0x8001129
> cpu MHz : 2056.406
> cache size : 512 KB
> physical id : 0

1700X is listed as base frequency=3400MHz, Turbo=3800MHz and "extra turbo"=3900MHz.
Even at 3900MHz I'd expect that it wouldn't run the variant of sqrtq(), supplied with MSYS2
faster than ~3.0 sec.
So, it has to be better version, somehow.
Unless the library does something that heavily favors AMD. I can't even think what it could be.
May be, when back at work, I'd test on some AMD CPU, but it seems that all AMDs at work are
of Zen3 variety.

Could you test on Intel?

> >
> >> > Another, not very probable possibility is that your cpu_time() call is lying, but that's easily verifiable with running
> >> > binary under 'time' utility.
> >> > And yet another possibility is that your compiler managed to parallelize a line 'a = sqrt(a)'. Crazy suggestion, I know.
> >> >
> >> > But one thing is sure - [on Intel] my quad-precision sqrt()
> >> > is close to 15 times faster than the one supplied with GCC
> >> > Quad-Precision Math Library.
> >> Have you considered submitting a patch?
> >
> > My quad-precision format is different from IEEE binary128==GNU __float128.
> > So, it would take some work.
> > And since I am not a great fun of GPL, I'd rather prefer my own work to be given away under
> > more anarchistic license.
> Unfortunately, this is not what I (usually) do - I am mainly
> concerned with the gfortran front end.
> > But I don't mind if you or your friends take my source and adapte it to quadmath.
> > Can put a needed license statement in my public github repository if somebody tells
> > me a relevant lawyerspeak.
> I can at least submit a PR. If there is a description of the
> algorithm in the source, maybe somebody can take this up and run
> with it.

Description (comments) is minimal, unfortunately.
But code should be easily understood by person, skilled in the field.

On Saturday, March 26, 2022 at 11:30:44 PM UTC+3, Michael S wrote:
> On Saturday, March 26, 2022 at 10:53:09 PM UTC+3, Thomas Koenig wrote:
> > Michael S <already...@yahoo.com> schrieb:
> > > On Saturday, March 26, 2022 at 12:45:45 PM UTC+3, Thomas Koenig wrote:
> > >> Michael S <already...@yahoo.com> schrieb:
> > >> > For reference, I compiled your program (MSYS2, gfortran 11.2.0, -O2) and run it on all tree systems.
> > >> > Core i5-3450: 3.75962 seconds
> > >> > Xeon E3-1271 v3: 2.96402 seconds
> > >> > Xeon E-2176G: 2.70312 seconds
> > >> >
> > >> > So, it seems, either your Zen1 CPU runs at much higher frequency (over 4.5 GHz?) or your version of quadmath library
> > >> > is much better than the one, supplied with MSYS2 (mingw-w64-x86_64-gcc-libgfortran 11.2.0-10) or the library
> > >> > likes AMD CPUs and hates Intel's.
> > >> Or something is going away in the calling sequence of the quadmath
> > >> library.
> > >>
> > >> If you have Linux installed, a timing comparison would be of interest.
> > >
> > > No, not at home.
> > > Even at work, I am not sure that I have Linux with Fortran compiler installed.
> > >
> > > Could you run you program under 'time' ?
> > > At very least it would prove that no autopar involved.
> > As a maintainer of gfortran, I can assure you that there is no
> > autoparallelization :-).
> >
> > However, here it is:
> >
> > $ cat a.f90
> > program main
> > implicit none
> > integer, parameter :: qp = selected_real_kind(30)
> > integer, parameter :: n = 10**6, m = 10
> > real :: t1, t2
> > character(len=20) :: c
> > integer :: i,k
> > real(kind=qp), dimension(:), allocatable :: a
> > allocate (a(n))
> > c = '10'
> > call random_number(a)
> > a = a * 1e50_qp
> > call cpu_time(t1)
> > do i=1,m
> > a = sqrt(a)
> > read (unit=c,fmt=*) k
> > write (*,*) a(k)
> > end do
> > call cpu_time(t2)
> > write (*,'(A,F12.5,A,ES15.5,A)') "Used ",t2-t1, " seconds to calculate ", &
> > real(n*m,qp), " square roots."
> > end program main
> > $ gfortran -O3 a.f90
> > $ time ./a.out
> > 7308311667482999228141838.93337335438
> > 2703388922719.59256619755162850493042
> > 1644198.56547790254000409031824041555
> > 1282.26306406988988805467367099802872
> > 35.8087009547943513227570444215179866
> > 5.98403717859392575839410139231541218
> > 2.44622917540322165097869738482238220
> > 1.56404257467730770560212117881721192
> > 1.25061687765570623839109221766572087
> > 1.11830983079632553071493517444861311
> > Used 2.35760 seconds to calculate 1.00000E+07 square roots.
> >
> > real 0m2,475s
> > user 0m2,445s
> > sys 0m0,008s
> > $ head /proc/cpuinfo
> > processor : 0
> > vendor_id : AuthenticAMD
> > cpu family : 23
> > model : 1
> > model name : AMD Ryzen 7 1700X Eight-Core Processor
> > stepping : 1
> > microcode : 0x8001129
> > cpu MHz : 2056.406
> > cache size : 512 KB
> > physical id : 0
> 1700X is listed as base frequency=3400MHz, Turbo=3800MHz and "extra turbo"=3900MHz.
> Even at 3900MHz I'd expect that it wouldn't run the variant of sqrtq(), supplied with MSYS2
> faster than ~3.0 sec.
> So, it has to be better version, somehow.
> Unless the library does something that heavily favors AMD. I can't even think what it could be.
> May be, when back at work, I'd test on some AMD CPU, but it seems that all AMDs at work are
> of Zen3 variety.
>
> Could you test on Intel?

On my system(s) I see another very strange timing effect:
The speed of sqrtq() strongly depends on the range of the input and does it in unexpected manner:
it is much faster (1.5x to 2x) when input is *outside* the range of double precision!

On Saturday, March 26, 2022 at 5:14:44 PM UTC-5, Michael S wrote:
> On Saturday, March 26, 2022 at 11:30:44 PM UTC+3, Michael S wrote:
> > On Saturday, March 26, 2022 at 10:53:09 PM UTC+3, Thomas Koenig wrote:
> > > Michael S <already...@yahoo.com> schrieb:
> > > > On Saturday, March 26, 2022 at 12:45:45 PM UTC+3, Thomas Koenig wrote:
> > > >> Michael S <already...@yahoo.com> schrieb:
> > > >> > For reference, I compiled your program (MSYS2, gfortran 11.2.0, -O2) and run it on all tree systems.
> > > >> > Core i5-3450: 3.75962 seconds
> > > >> > Xeon E3-1271 v3: 2.96402 seconds
> > > >> > Xeon E-2176G: 2.70312 seconds
> > > >> >
> > > >> > So, it seems, either your Zen1 CPU runs at much higher frequency (over 4.5 GHz?) or your version of quadmath library
> > > >> > is much better than the one, supplied with MSYS2 (mingw-w64-x86_64-gcc-libgfortran 11.2.0-10) or the library
> > > >> > likes AMD CPUs and hates Intel's.
> > > >> Or something is going away in the calling sequence of the quadmath
> > > >> library.
> > > >>
> > > >> If you have Linux installed, a timing comparison would be of interest.
> > > >
> > > > No, not at home.
> > > > Even at work, I am not sure that I have Linux with Fortran compiler installed.
> > > >
> > > > Could you run you program under 'time' ?
> > > > At very least it would prove that no autopar involved.
> > > As a maintainer of gfortran, I can assure you that there is no
> > > autoparallelization :-).
> > >
> > > However, here it is:
> > >
> > > $ cat a.f90
> > > program main
> > > implicit none
> > > integer, parameter :: qp = selected_real_kind(30)
> > > integer, parameter :: n = 10**6, m = 10
> > > real :: t1, t2
> > > character(len=20) :: c
> > > integer :: i,k
> > > real(kind=qp), dimension(:), allocatable :: a
> > > allocate (a(n))
> > > c = '10'
> > > call random_number(a)
> > > a = a * 1e50_qp
> > > call cpu_time(t1)
> > > do i=1,m
> > > a = sqrt(a)
> > > read (unit=c,fmt=*) k
> > > write (*,*) a(k)
> > > end do
> > > call cpu_time(t2)
> > > write (*,'(A,F12.5,A,ES15.5,A)') "Used ",t2-t1, " seconds to calculate ", &
> > > real(n*m,qp), " square roots."
> > > end program main
> > > $ gfortran -O3 a.f90
> > > $ time ./a.out
> > > 7308311667482999228141838.93337335438
> > > 2703388922719.59256619755162850493042
> > > 1644198.56547790254000409031824041555
> > > 1282.26306406988988805467367099802872
> > > 35.8087009547943513227570444215179866
> > > 5.98403717859392575839410139231541218
> > > 2.44622917540322165097869738482238220
> > > 1.56404257467730770560212117881721192
> > > 1.25061687765570623839109221766572087
> > > 1.11830983079632553071493517444861311
> > > Used 2.35760 seconds to calculate 1.00000E+07 square roots.
> > >
> > > real 0m2,475s
> > > user 0m2,445s
> > > sys 0m0,008s
> > > $ head /proc/cpuinfo
> > > processor : 0
> > > vendor_id : AuthenticAMD
> > > cpu family : 23
> > > model : 1
> > > model name : AMD Ryzen 7 1700X Eight-Core Processor
> > > stepping : 1
> > > microcode : 0x8001129
> > > cpu MHz : 2056.406
> > > cache size : 512 KB
> > > physical id : 0
> > 1700X is listed as base frequency=3400MHz, Turbo=3800MHz and "extra turbo"=3900MHz.
> > Even at 3900MHz I'd expect that it wouldn't run the variant of sqrtq(), supplied with MSYS2
> > faster than ~3.0 sec.
> > So, it has to be better version, somehow.
> > Unless the library does something that heavily favors AMD. I can't even think what it could be.
> > May be, when back at work, I'd test on some AMD CPU, but it seems that all AMDs at work are
> > of Zen3 variety.
> >
> > Could you test on Intel?
<
> On my system(s) I see another very strange timing effect:
> The speed of sqrtq() strongly depends on the range of the input and does it in unexpected manner:
> it is much faster (1.5x to 2x) when input is *outside* the range of double precision!
<
Can't rely on :: CVT Q->D when Q would overflow D.
<
But the algorithm only needs the LoB of the exponent and 8-12 MsBs of the fraction to get started.
<
First 3 iterations are SP, next iteration is DP, next 2 are QP.

Michael S <already5chosen@yahoo.com> schrieb:

> On my system(s) I see another very strange timing effect:
> The speed of sqrtq() strongly depends on the range of the input and does it in unexpected manner:
> it is much faster (1.5x to 2x) when input is *outside* the range of double precision!

Source for the sqrt routine in libquadmath is here:

https://gcc.gnu.org/git/?p=gcc.git;a=blob;f=libquadmath/math/sqrtq.c;h=56ea5d3243c06df2276689849b448eebacae6367;hb=refs/heads/master

Hmm.. they first do a test if it is within double precision range,
with two Newton iterations, then a test if it is within long double
range with a single Newton iteration, and if it is outside then
they pick apart the number and run two Newton iterations.

Were the numbers outside the range of double inside the range
of long double? That could explain things (and suggest
an improvement).

On Sunday, March 27, 2022 at 1:25:29 AM UTC+3, Thomas Koenig wrote:
> Michael S <already...@yahoo.com> schrieb:
> > On my system(s) I see another very strange timing effect:
> > The speed of sqrtq() strongly depends on the range of the input and does it in unexpected manner:
> > it is much faster (1.5x to 2x) when input is *outside* the range of double precision!
> Source for the sqrt routine in libquadmath is here:
>
> https://gcc.gnu.org/git/?p=gcc.git;a=blob;f=libquadmath/math/sqrtq.c;h=56ea5d3243c06df2276689849b448eebacae6367;hb=refs/heads/master
>
> Hmm.. they first do a test if it is within double precision range,
> with two Newton iterations, then a test if it is within long double
> range with a single Newton iteration, and if it is outside then
> they pick apart the number and run two Newton iterations.
>
> Were the numbers outside the range of double inside the range
> of long double?

Yes.
I didn't test "outside long double".
If "long double" is 80-bit x87 format then "outside" would be mostly sub-normals, right?

> That could explain things (and suggest
> an improvement).

I think, the whole routine needs rewrite, rather than improvements.
The original author probably was not thinking very hard.

On 2022-03-22, Marcus wrote:
> Hello group!
>
> A class of instructions that is very tempting to include in an ISA is
> approximate floating-point reciprocals (1/x & 1/sqrt(x)), and
> possibly specialized instructions for improving precision
> (a Newton-Raphson step).
>
> Many ISAs have such instructions:
>
> CRAY:
> * 070ijx - Floating reciprocal approximation
> * 067ijk - Reciprocal iteration
>
> ARM:
> * FRECPE - Floating-point reciprocal estimate
> * FRECPS - Floating-point reciprocal step
> * FRSQRTE - Floating-point reciprocal square root estimate
> * FRSQRTS - Floating-point reciprocal square root step
>
> POWER:
> * FRES - Floating reciprocal estimate single
> * FRSQRTE - Floating reciprocal square root estimate
>
> x86:
> * RSQRTSS - Approximate reciprocal square root
>
> TI C67x:
> * RCPSP - Floating-Point reciprocal approximation
> * RSQRSP - Floating-Point reciprocal square root approximation
>
> ...and there are probably others.
>
> What are your feelings towards including such instructions in an ISA?
>
> My own feelings are mixed.
>
> Pros:
> * Easy to implement in hardware.
> * Can provide a significant speedup for certain workloads, especially
>   if limited accuracy is acceptable.
>
> Cons:
> * Hard to specify exact operation.
> * Likely a source of poor portability (borderline undefined behavior).
> * The "step/iteration" instructions can usually be replaced by FMA.

As an aside...

Today I found a pretty nice looking 2nd order polynomial for
approximating 1/x in the interval [1.0, 2.0):

(2*x^2 - 9*x + 13) / 6

Maybe this is known stuff, and I'm not sure that it would be useful for
anything, but I thought I'd share it anyway.

The accuracy that you get is roughly 6-7 bits, so it's not better than a
LUT.

The coefficients are "nice" from an implementation perspective, if you
were to hard-wire the polynomial. It should be possible to implement the
dividend with a single multiplier and a few integer adders. The 3 in the
divisor is muddying the waters, though. Is there a fast way to divide by
three?

/Marcus

Marcus <m.delete@this.bitsnbites.eu> schrieb:

> Today I found a pretty nice looking 2nd order polynomial for
> approximating 1/x in the interval [1.0, 2.0):
>
> (2*x^2 - 9*x + 13) / 6

Interesting. Could be useful for the start of a Newton iteration
for a reciprocal, where the iteration formula would be
x_n+1 = 2*x_n - a * x_n^2 (division-free), or x_n*(2-a*x_n).

>
> Maybe this is known stuff, and I'm not sure that it would be useful for
> anything, but I thought I'd share it anyway.
>
> The accuracy that you get is roughly 6-7 bits, so it's not better than a
> LUT.
>
> The coefficients are "nice" from an implementation perspective, if you
> were to hard-wire the polynomial. It should be possible to implement the
> dividend with a single multiplier and a few integer adders. The 3 in the
> divisor is muddying the waters, though. Is there a fast way to divide by
> three?

The standard trick of dividing by three in fixed point: For
32-bit numbers, multiply the number by the magic number 2863311531
(0xAAAAAAAB) and shift right by 33 bits (in practice, use a
"multiply high" and shift right one bit). For floating point
numbers, this would have to be adjusted somewhat.

See Hacker's Delight, "Unsigned Division by Divisors >= 1".

Or just encode the constants.

On Monday, March 28, 2022 at 12:26:48 AM UTC+3, Thomas Koenig wrote:
> Marcus <m.de...@this.bitsnbites.eu> schrieb:
> > Today I found a pretty nice looking 2nd order polynomial for
> > approximating 1/x in the interval [1.0, 2.0):
> >
> > (2*x^2 - 9*x + 13) / 6
> Interesting. Could be useful for the start of a Newton iteration
> for a reciprocal, where the iteration formula would be
> x_n+1 = 2*x_n - a * x_n^2 (division-free), or x_n*(2-a*x_n).
> >
> > Maybe this is known stuff, and I'm not sure that it would be useful for
> > anything, but I thought I'd share it anyway.
> >
> > The accuracy that you get is roughly 6-7 bits, so it's not better than a
> > LUT.
> >
> > The coefficients are "nice" from an implementation perspective, if you
> > were to hard-wire the polynomial. It should be possible to implement the
> > dividend with a single multiplier and a few integer adders. The 3 in the
> > divisor is muddying the waters, though. Is there a fast way to divide by
> > three?
> The standard trick of dividing by three in fixed point: For
> 32-bit numbers, multiply the number by the magic number 2863311531
> (0xAAAAAAAB) and shift right by 33 bits (in practice, use a
> "multiply high" and shift right one bit). For floating point
> numbers, this would have to be adjusted somewhat.
>

And that causes the whole calculation to be of approximately the same
computational complexity as 2nd-order polynomials with arbitrary 32-bit coefficients,
potentially much better coefficients than those suggested.
The suggestion of Marcus is very naive.
Any suggestion to do sqrt or rsqrt approximation without lookup table is a loose proposition.
Even very small table at the first step, like 32 or 64 entries, helps a lot and saves many steps
down the road.
The only exception to that could be when you do very wide SIMD that has no tools for even small lookups.
Or when you don't care about speed, but that case is sort of off topic.

>
> See Hacker's Delight, "Unsigned Division by Divisors >= 1".
>
> Or just encode the constants.

On Sunday, March 27, 2022 at 2:47:45 PM UTC-5, Marcus wrote:
> On 2022-03-22, Marcus wrote:
> > Hello group!
> >
> > A class of instructions that is very tempting to include in an ISA is
> > approximate floating-point reciprocals (1/x & 1/sqrt(x)), and
> > possibly specialized instructions for improving precision
> > (a Newton-Raphson step).
> >
> > Many ISAs have such instructions:
> >
> > CRAY:
> > * 070ijx - Floating reciprocal approximation
> > * 067ijk - Reciprocal iteration
> >
> > ARM:
> > * FRECPE - Floating-point reciprocal estimate
> > * FRECPS - Floating-point reciprocal step
> > * FRSQRTE - Floating-point reciprocal square root estimate
> > * FRSQRTS - Floating-point reciprocal square root step
> >
> > POWER:
> > * FRES - Floating reciprocal estimate single
> > * FRSQRTE - Floating reciprocal square root estimate
> >
> > x86:
> > * RSQRTSS - Approximate reciprocal square root
> >
> > TI C67x:
> > * RCPSP - Floating-Point reciprocal approximation
> > * RSQRSP - Floating-Point reciprocal square root approximation
> >
> > ...and there are probably others.
> >
> > What are your feelings towards including such instructions in an ISA?
> >
> > My own feelings are mixed.
> >
> > Pros:
> > * Easy to implement in hardware.
> > * Can provide a significant speedup for certain workloads, especially
> > if limited accuracy is acceptable.
> >
> > Cons:
> > * Hard to specify exact operation.
> > * Likely a source of poor portability (borderline undefined behavior).
> > * The "step/iteration" instructions can usually be replaced by FMA.
> As an aside...
>
> Today I found a pretty nice looking 2nd order polynomial for
> approximating 1/x in the interval [1.0, 2.0):
>
> (2*x^2 - 9*x + 13) / 6
<
Equivalent to::
<
.33333×x^2 - 1.5×x + 2.1666
<
The Chebychev 2nd order Coefficients on the interval [1..2) are:
<
0.32323232×x^2 -0.48484848×x + 0.66666667
<
with 6.63 bits of accuracy. BTW the 9th order cheby polynomial only has 21.89
bits of accuracy. Clearly, this is not the right way to start.
<
Splitting the interval into [1..1.5)..[1.5..2] gives 8.9 bits of accuracy
Splitting into 4 intervals gives 11.5 bits of accuracy.
Splitting into 8 intervals gives 14.25 bits of accuracy.
>
> Maybe this is known stuff, and I'm not sure that it would be useful for
> anything, but I thought I'd share it anyway.
>
> The accuracy that you get is roughly 6-7 bits, so it's not better than a
> LUT.
>
> The coefficients are "nice" from an implementation perspective, if you
> were to hard-wire the polynomial. It should be possible to implement the
> dividend with a single multiplier and a few integer adders. The 3 in the
> divisor is muddying the waters, though. Is there a fast way to divide by
> three?
>
> /Marcus

"MitchAlsup" wrote in message
news:600fb4d3-5e5f-490f-b7b6-4301a93f656en@googlegroups.com...

> First 3 iterations are SP, next iteration is DP, next 2 are QP.

Have you investigated the potential of using variable numbers of
terms from the series

x/(1-(1-x*D)) = x*sum([((1-x*D)**n,n=0,∞)])
x/sqrt(1-(1-x**2*D)) =
x*sum([(gamma(n+0.5)/(sqrt(pi)*gamma(n+1.0))*(1-x**2*D)**n,n=0,∞)])

to perhaps remove an iteration somewhere?

Michael S wrote:
> On Monday, March 28, 2022 at 12:26:48 AM UTC+3, Thomas Koenig wrote:
>> Marcus <m.de...@this.bitsnbites.eu> schrieb:
>>> Today I found a pretty nice looking 2nd order polynomial for
>>> approximating 1/x in the interval [1.0, 2.0):
>>>
>>> (2*x^2 - 9*x + 13) / 6
>> Interesting. Could be useful for the start of a Newton iteration
>> for a reciprocal, where the iteration formula would be
>> x_n+1 = 2*x_n - a * x_n^2 (division-free), or x_n*(2-a*x_n).
>>>
>>> Maybe this is known stuff, and I'm not sure that it would be useful for
>>> anything, but I thought I'd share it anyway.
>>>
>>> The accuracy that you get is roughly 6-7 bits, so it's not better than a
>>> LUT.
>>>
>>> The coefficients are "nice" from an implementation perspective, if you
>>> were to hard-wire the polynomial. It should be possible to implement the
>>> dividend with a single multiplier and a few integer adders. The 3 in the
>>> divisor is muddying the waters, though. Is there a fast way to divide by
>>> three?
>> The standard trick of dividing by three in fixed point: For
>> 32-bit numbers, multiply the number by the magic number 2863311531
>> (0xAAAAAAAB) and shift right by 33 bits (in practice, use a
>> "multiply high" and shift right one bit). For floating point
>> numbers, this would have to be adjusted somewhat.
>>
>
> And that causes the whole calculation to be of approximately the same
> computational complexity as 2nd-order polynomials with arbitrary 32-bit coefficients,
> potentially much better coefficients than those suggested.
> The suggestion of Marcus is very naive.
> Any suggestion to do sqrt or rsqrt approximation without lookup table is a loose proposition.
> Even very small table at the first step, like 32 or 64 entries, helps a lot and saves many steps
> down the road.
> The only exception to that could be when you do very wide SIMD that has no tools for even small lookups.
> Or when you don't care about speed, but that case is sort of off topic.

If you are doing SIMD, but without a reciprocal/reciprocal square root
lookup opcode, then I would use the infamous invsqrt() trick. With the
latest Cheby-style tweaks to all the constants, this delivers 10+ bits
after a single iteration.

This is significantly better than an in-register 4-bit/16-term table lookup.

Anyway, I am sure Mitch has much better polynomials to do it all
directly. :-)

Terje

--
- <Terje.Mathisen at tmsw.no>
"almost all programming can be viewed as an exercise in caching"

On Monday, March 28, 2022 at 9:54:31 AM UTC+3, Terje Mathisen wrote:
> Michael S wrote:
> > On Monday, March 28, 2022 at 12:26:48 AM UTC+3, Thomas Koenig wrote:
> >> Marcus <m.de...@this.bitsnbites.eu> schrieb:
> >>> Today I found a pretty nice looking 2nd order polynomial for
> >>> approximating 1/x in the interval [1.0, 2.0):
> >>>
> >>> (2*x^2 - 9*x + 13) / 6
> >> Interesting. Could be useful for the start of a Newton iteration
> >> for a reciprocal, where the iteration formula would be
> >> x_n+1 = 2*x_n - a * x_n^2 (division-free), or x_n*(2-a*x_n).
> >>>
> >>> Maybe this is known stuff, and I'm not sure that it would be useful for
> >>> anything, but I thought I'd share it anyway.
> >>>
> >>> The accuracy that you get is roughly 6-7 bits, so it's not better than a
> >>> LUT.
> >>>
> >>> The coefficients are "nice" from an implementation perspective, if you
> >>> were to hard-wire the polynomial. It should be possible to implement the
> >>> dividend with a single multiplier and a few integer adders. The 3 in the
> >>> divisor is muddying the waters, though. Is there a fast way to divide by
> >>> three?
> >> The standard trick of dividing by three in fixed point: For
> >> 32-bit numbers, multiply the number by the magic number 2863311531
> >> (0xAAAAAAAB) and shift right by 33 bits (in practice, use a
> >> "multiply high" and shift right one bit). For floating point
> >> numbers, this would have to be adjusted somewhat.
> >>
> >
> > And that causes the whole calculation to be of approximately the same
> > computational complexity as 2nd-order polynomials with arbitrary 32-bit coefficients,
> > potentially much better coefficients than those suggested.
> > The suggestion of Marcus is very naive.
> > Any suggestion to do sqrt or rsqrt approximation without lookup table is a loose proposition.
> > Even very small table at the first step, like 32 or 64 entries, helps a lot and saves many steps
> > down the road.
> > The only exception to that could be when you do very wide SIMD that has no tools for even small lookups.
> > Or when you don't care about speed, but that case is sort of off topic.
> If you are doing SIMD, but without a reciprocal/reciprocal square root
> lookup opcode, then I would use the infamous invsqrt() trick.

Magic constant and such?
I don't think it's a good idea if you seek top performance.
One thing to realize is that if you went to trouble of using rsqrt() NR steps
based on custom polynomials, with different coefficients on each step, then it
is better to keep error one-sided on all steps except possibly the last one.
Apart from of being slightly more precise, such strategy has advantage of being
simpler to do with unsigned integer math.

> With the
> latest Cheby-style tweaks to all the constants, this delivers 10+ bits
> after a single iteration.
>
> This is significantly better than an in-register 4-bit/16-term table lookup.
>
> Anyway, I am sure Mitch has much better polynomials to do it all
> directly. :-)

Mitch appears to like Chebychev polynomials.
Likely, because they are easy to find out :-)

They are very close to optimal (in sense of minimizing two-sided error, but
probably can be easily offseted to minimize one-sided error too) when the
interval of input error is already small. But when the interval is large, as it
is a case on the first step without lookup table and, may be, on the second
step, when the 1st step was of low order, then Chebychev polynomials are
measurably suboptimal vs true equiripple polys. The later can be found by Remez
exchange algorithm.
The reason for this is that Chebychev poly of order N is defined as the best
approximation for poly of order N+1. So, when term N+1 in infinite Tailor
series is much bigger than sum of terms N+2 to Infinity then Cheby is good, but
when term N+1 is not dominant it is less good.

> Terje
>
> --
> - <Terje.Mathisen at tmsw.no>
> "almost all programming can be viewed as an exercise in caching"

Michael S wrote:
> On Monday, March 28, 2022 at 9:54:31 AM UTC+3, Terje Mathisen wrote:
>> Michael S wrote:
>>> On Monday, March 28, 2022 at 12:26:48 AM UTC+3, Thomas Koenig wrote:
>>>> Marcus <m.de...@this.bitsnbites.eu> schrieb:
>>>>> Today I found a pretty nice looking 2nd order polynomial for
>>>>> approximating 1/x in the interval [1.0, 2.0):
>>>>>
>>>>> (2*x^2 - 9*x + 13) / 6
>>>> Interesting. Could be useful for the start of a Newton iteration
>>>> for a reciprocal, where the iteration formula would be
>>>> x_n+1 = 2*x_n - a * x_n^2 (division-free), or x_n*(2-a*x_n).
>>>>>
>>>>> Maybe this is known stuff, and I'm not sure that it would be useful for
>>>>> anything, but I thought I'd share it anyway.
>>>>>
>>>>> The accuracy that you get is roughly 6-7 bits, so it's not better than a
>>>>> LUT.
>>>>>
>>>>> The coefficients are "nice" from an implementation perspective, if you
>>>>> were to hard-wire the polynomial. It should be possible to implement the
>>>>> dividend with a single multiplier and a few integer adders. The 3 in the
>>>>> divisor is muddying the waters, though. Is there a fast way to divide by
>>>>> three?
>>>> The standard trick of dividing by three in fixed point: For
>>>> 32-bit numbers, multiply the number by the magic number 2863311531
>>>> (0xAAAAAAAB) and shift right by 33 bits (in practice, use a
>>>> "multiply high" and shift right one bit). For floating point
>>>> numbers, this would have to be adjusted somewhat.
>>>>
>>>
>>> And that causes the whole calculation to be of approximately the same
>>> computational complexity as 2nd-order polynomials with arbitrary 32-bit coefficients,
>>> potentially much better coefficients than those suggested.
>>> The suggestion of Marcus is very naive.
>>> Any suggestion to do sqrt or rsqrt approximation without lookup table is a loose proposition.
>>> Even very small table at the first step, like 32 or 64 entries, helps a lot and saves many steps
>>> down the road.
>>> The only exception to that could be when you do very wide SIMD that has no tools for even small lookups.
>>> Or when you don't care about speed, but that case is sort of off topic.
>> If you are doing SIMD, but without a reciprocal/reciprocal square root
>> lookup opcode, then I would use the infamous invsqrt() trick.
>
> Magic constant and such?
> I don't think it's a good idea if you seek top performance.
> One thing to realize is that if you went to trouble of using rsqrt() NR steps
> based on custom polynomials, with different coefficients on each step, then it
> is better to keep error one-sided on all steps except possibly the last one.
> Apart from of being slightly more precise, such strategy has advantage of being
> simpler to do with unsigned integer math.
>
>> With the
>> latest Cheby-style tweaks to all the constants, this delivers 10+ bits
>> after a single iteration.
>>
>> This is significantly better than an in-register 4-bit/16-term table lookup.
>>
>> Anyway, I am sure Mitch has much better polynomials to do it all
>> directly. :-)
>
> Mitch appears to like Chebychev polynomials.
> Likely, because they are easy to find out :-)
>
> They are very close to optimal (in sense of minimizing two-sided error, but
> probably can be easily offseted to minimize one-sided error too) when the
> interval of input error is already small. But when the interval is large, as it
> is a case on the first step without lookup table and, may be, on the second
> step, when the 1st step was of low order, then Chebychev polynomials are
> measurably suboptimal vs true equiripple polys. The later can be found by Remez
> exchange algorithm.

OK, I might have been misusing the "Cheby" term! I've always assumed
that an optimal N-term poly will have N+1 ripples, all equally far
(typically using a relative measure) away from the true solution at
those points. Sometimes I'll additional restrictions, like wanting the
error at both ends to be identical in sign and magnitude which is a good
idea if I use range reduction and individual coefficients for each
sub-range.

> The reason for this is that Chebychev poly of order N is defined as the best
> approximation for poly of order N+1. So, when term N+1 in infinite Tailor
> series is much bigger than sum of terms N+2 to Infinity then Cheby is good, but
> when term N+1 is not dominant it is less good.

OK,
thanks!

Terje

--
- <Terje.Mathisen at tmsw.no>
"almost all programming can be viewed as an exercise in caching"

On Monday, March 28, 2022 at 3:19:46 AM UTC-5, Michael S wrote:
> On Monday, March 28, 2022 at 9:54:31 AM UTC+3, Terje Mathisen wrote:
> > Michael S wrote:
> > > On Monday, March 28, 2022 at 12:26:48 AM UTC+3, Thomas Koenig wrote:
> > >> Marcus <m.de...@this.bitsnbites.eu> schrieb:
> > >>> Today I found a pretty nice looking 2nd order polynomial for
> > >>> approximating 1/x in the interval [1.0, 2.0):
> > >>>
> > >>> (2*x^2 - 9*x + 13) / 6
> > >> Interesting. Could be useful for the start of a Newton iteration
> > >> for a reciprocal, where the iteration formula would be
> > >> x_n+1 = 2*x_n - a * x_n^2 (division-free), or x_n*(2-a*x_n).
> > >>>
> > >>> Maybe this is known stuff, and I'm not sure that it would be useful for
> > >>> anything, but I thought I'd share it anyway.
> > >>>
> > >>> The accuracy that you get is roughly 6-7 bits, so it's not better than a
> > >>> LUT.
> > >>>
> > >>> The coefficients are "nice" from an implementation perspective, if you
> > >>> were to hard-wire the polynomial. It should be possible to implement the
> > >>> dividend with a single multiplier and a few integer adders. The 3 in the
> > >>> divisor is muddying the waters, though. Is there a fast way to divide by
> > >>> three?
> > >> The standard trick of dividing by three in fixed point: For
> > >> 32-bit numbers, multiply the number by the magic number 2863311531
> > >> (0xAAAAAAAB) and shift right by 33 bits (in practice, use a
> > >> "multiply high" and shift right one bit). For floating point
> > >> numbers, this would have to be adjusted somewhat.
> > >>
> > >
> > > And that causes the whole calculation to be of approximately the same
> > > computational complexity as 2nd-order polynomials with arbitrary 32-bit coefficients,
> > > potentially much better coefficients than those suggested.
> > > The suggestion of Marcus is very naive.
> > > Any suggestion to do sqrt or rsqrt approximation without lookup table is a loose proposition.
> > > Even very small table at the first step, like 32 or 64 entries, helps a lot and saves many steps
> > > down the road.
> > > The only exception to that could be when you do very wide SIMD that has no tools for even small lookups.
> > > Or when you don't care about speed, but that case is sort of off topic.
> > If you are doing SIMD, but without a reciprocal/reciprocal square root
> > lookup opcode, then I would use the infamous invsqrt() trick.
<
> Magic constant and such?
> I don't think it's a good idea if you seek top performance.
> One thing to realize is that if you went to trouble of using rsqrt() NR steps
> based on custom polynomials, with different coefficients on each step, then it
> is better to keep error one-sided on all steps except possibly the last one.
<
Notice that the error term for N-R iterations flip the sign every iteration.
<
> Apart from of being slightly more precise, such strategy has advantage of being
> simpler to do with unsigned integer math.
<
You need to choose the approximat such that the error is uniformly positive after
the last iteration. This allows you to choose from {r+0, r+1, and r+2} for the properly
rounded result; which is easier than choosing {r-1, r+0, r+1} for the properly rounded
result.
<
> > With the
> > latest Cheby-style tweaks to all the constants, this delivers 10+ bits
> > after a single iteration.
> >
> > This is significantly better than an in-register 4-bit/16-term table lookup.
> >
> > Anyway, I am sure Mitch has much better polynomials to do it all
> > directly. :-)
> Mitch appears to like Chebychev polynomials.
> Likely, because they are easy to find out :-)
<
Once you figure out how to calculate them, they are, indeed, easy to determine.
>
> They are very close to optimal (in sense of minimizing two-sided error, but
> probably can be easily offseted to minimize one-sided error too) when the
> interval of input error is already small. But when the interval is large, as it
> is a case on the first step without lookup table and, may be, on the second
> step, when the 1st step was of low order, then Chebychev polynomials are
> measurably suboptimal vs true equiripple polys. The later can be found by Remez
> exchange algorithm.
<
Cheby polynomials are within 1 ULP of optimal, where Remez is better, it is
a lot harder to calculate, and if you are working to the last bit of precision,
even calculating Remez requires more precision. Cheby does not have this
property.
<
> The reason for this is that Chebychev poly of order N is defined as the best
> approximation for poly of order N+1. So, when term N+1 in infinite Tailor
> series is much bigger than sum of terms N+2 to Infinity then Cheby is good, but
> when term N+1 is not dominant it is less good.
<
I calculated my SIN(), COS() Cheby polynomials by starting with a Taylor series
3 terms longer and then removed N+3, N+2, and N+1 to give me the N-th order
cheby polynomial in x^2
<
> > Terje
> >
> > --
> > - <Terje.Mathisen at tmsw.no>
> > "almost all programming can be viewed as an exercise in caching"