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.

/Marcus

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?

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.

>
> 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.

This is true, which is why only the initial lookup is crucial.

Terje

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

On Tuesday, March 22, 2022 at 9:24:58 AM UTC-4, Terje Mathisen wrote:
> 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?
> 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.
> >
> > 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.
> This is true, which is why only the initial lookup is crucial.
>
> Terje
>
> --
> - <Terje.Mathisen at tmsw.no>
> "almost all programming can be viewed as an exercise in caching"

Should be included.

I have sigmoid approximate too.

On Tuesday, March 22, 2022 at 3:25:11 AM UTC-5, 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.
<
My feelings are::
a) screw the approximations
b) because doing them to faithful accuracy is easy
c) the right answer is only 14 cycles
d) if you build the FMAC unit correctly
See: USPTO 10,761,806
>
> 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.
>
> /Marcus

On 2022-03-22, Terje Mathisen wrote:
> 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?
>
> 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.

>
>>
>> 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.
>
> This is true, which is why only the initial lookup is crucial.
>
> Terje
>

On 2022-03-22, MitchAlsup wrote:
> On Tuesday, March 22, 2022 at 3:25:11 AM UTC-5, 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.
> <
> My feelings are::
> a) screw the approximations
> b) because doing them to faithful accuracy is easy
> c) the right answer is only 14 cycles
> d) if you build the FMAC unit correctly
> See: USPTO 10,761,806

Thanks, that makes for a good read!

While I see the benefits of fast argument reduction and a generalized
coefficient + FMA engine for implementing transcendental functions, the
reciprocal and reciprocal sqrt functions are simpler and common enough
to possibly warrant a dedicated solution?

E.g. the CRAY could do a full DP division by means of approximate
reciprocal + NR-step + multiply by nominator in just 3-4 clock cycles
(vectorized), which is pretty impressive even by today's standards.

>>
>> 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.
>>
>> /Marcus

On Tuesday, March 22, 2022 at 10:26:47 AM UTC-5, Marcus wrote:
> On 2022-03-22, Terje Mathisen wrote:
> > 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?
> >
> > 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?
<
0.510 ULP or better Double Precision.
<
> 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 (!).
<
There are a couple hundred cases where you have 26-bits of accuracy
do a Newton-Raphson iteration and you still cannot round correctly,
so you have to do a N-R iteration after you get more than 52-bit of
precision.
>
> 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.
<
I do a 7-th order polynomial
> >
> >>
> >> 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.
> >
> > This is true, which is why only the initial lookup is crucial.
> >
> > Terje
> >

On Tuesday, March 22, 2022 at 10:56:06 AM UTC-5, Marcus wrote:
> On 2022-03-22, MitchAlsup wrote:
> > On Tuesday, March 22, 2022 at 3:25:11 AM UTC-5, 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.
> > <
> > My feelings are::
> > a) screw the approximations
> > b) because doing them to faithful accuracy is easy
> > c) the right answer is only 14 cycles
> > d) if you build the FMAC unit correctly
> > See: USPTO 10,761,806
> Thanks, that makes for a good read!
>
> While I see the benefits of fast argument reduction and a generalized
> coefficient + FMA engine for implementing transcendental functions, the
> reciprocal and reciprocal sqrt functions are simpler and common enough
> to possibly warrant a dedicated solution?
<
You get to make those choices.
>
> E.g. the CRAY could do a full DP division by means of approximate
> reciprocal + NR-step + multiply by nominator in just 3-4 clock cycles
> (vectorized), which is pretty impressive even by today's standards.
<
Its OK if you are content with 46-bits of accuracy.
<
> >>
> >> 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.
> >>
> >> /Marcus

MitchAlsup <MitchAlsup@aol.com> schrieb:

> My feelings are::
> a) screw the approximations
> b) because doing them to faithful accuracy is easy
> c) the right answer is only 14 cycles
> d) if you build the FMAC unit correctly
> See: USPTO 10,761,806

Impressive (as I said before). This being a patent means that
you forbid everybody to do what you claim, except outside the US
(I believe Marcus is based in Sweden) and, of course, research is
OK too. It is unclear if an FPGA-based implementation would fall
under this patent ("processor" does not seem to be defined, and as
I have learnt the hard way, you need to define things in patents).

Regarding approximate reciprocals: I would prefer a CPU which has
different widths of floating point data (16, 32, 64 and possibly 128
bits), if the smaller types are faster and more energy efficient
than the larger types. People who need a low-cost, low-accuracy
sqrt could then use the 16-bit version.

On Tuesday, March 22, 2022 at 12:26:03 PM UTC-5, Thomas Koenig wrote:
> MitchAlsup <Mitch...@aol.com> schrieb:
> > My feelings are::
> > a) screw the approximations
> > b) because doing them to faithful accuracy is easy
> > c) the right answer is only 14 cycles
> > d) if you build the FMAC unit correctly
> > See: USPTO 10,761,806
<
> Impressive (as I said before). This being a patent means that
> you forbid everybody to do what you claim, except outside the US
> (I believe Marcus is based in Sweden) and, of course, research is
> OK too. It is unclear if an FPGA-based implementation would fall
> under this patent ("processor" does not seem to be defined, and as
> I have learnt the hard way, you need to define things in patents).
<
I am willing to sell licenses for less money that it takes you to
assign a design engineer and have him/her read the document
thoroughly.
<
What is claimed is the necessary organization. An FPGA with
the necessary organization would infringe.
>
> Regarding approximate reciprocals: I would prefer a CPU which has
> different widths of floating point data (16, 32, 64 and possibly 128
> bits), if the smaller types are faster and more energy efficient
> than the larger types. People who need a low-cost, low-accuracy
> sqrt could then use the 16-bit version.
<
I am nor arguing against that point. From GPU work, the single precision
transcendentals take 5 cycles for 1 ULP accuracy (not quite faithful),
1 cycle throughput. Due to the organization of SP trans, 16-bit takes
just as long as 32-bit <in GPUs>.
<
In order to get this working at the 128-bit level, you are going to need
FP with 128-bits of fraction; just so you can correctly calculate the
coefficients the polynomial is based upon. My estimate is that 128-bit
transcendentals will cost ~40 cycles.

MitchAlsup <MitchAlsup@aol.com> schrieb:
> On Tuesday, March 22, 2022 at 12:26:03 PM UTC-5, Thomas Koenig wrote:
>> MitchAlsup <Mitch...@aol.com> schrieb:
>> > My feelings are::
>> > a) screw the approximations
>> > b) because doing them to faithful accuracy is easy
>> > c) the right answer is only 14 cycles
>> > d) if you build the FMAC unit correctly
>> > See: USPTO 10,761,806
><
>> Impressive (as I said before). This being a patent means that
>> you forbid everybody to do what you claim, except outside the US
>> (I believe Marcus is based in Sweden) and, of course, research is
>> OK too. It is unclear if an FPGA-based implementation would fall
>> under this patent ("processor" does not seem to be defined, and as
>> I have learnt the hard way, you need to define things in patents).
><
> I am willing to sell licenses for less money that it takes you to
> assign a design engineer and have him/her read the document
> thoroughly.

I'm not going to pay for this, my profession is something
quite different :-)

> What is claimed is the necessary organization. An FPGA with
> the necessary organization would infringe.

Hm, if I had been involved in writing that patent, I would probably
have included a description of what "processor" actually means
in this context (which is why stuff I (co)write tends to be much
longer, possibly unnecessarily so; it is, however, in a field
totally unrelated to computer architecture).

But it is of course your patent and your style of writing :-)

>> Regarding approximate reciprocals: I would prefer a CPU which has
>> different widths of floating point data (16, 32, 64 and possibly 128
>> bits), if the smaller types are faster and more energy efficient
>> than the larger types. People who need a low-cost, low-accuracy
>> sqrt could then use the 16-bit version.
><
> I am nor arguing against that point. From GPU work, the single precision
> transcendentals take 5 cycles for 1 ULP accuracy (not quite faithful),
> 1 cycle throughput. Due to the organization of SP trans, 16-bit takes
> just as long as 32-bit <in GPUs>.
><
> In order to get this working at the 128-bit level, you are going to need
> FP with 128-bits of fraction; just so you can correctly calculate the
> coefficients the polynomial is based upon. My estimate is that 128-bit
> transcendentals will cost ~40 cycles.

Certainly faster than doing the same thing in software.

The little test program

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)
do i=1,m
a = sqrt(a)
read (unit=c,fmt=*) k
write (*,*) a(k)
end do
end program main

calculates 10**7 16-byte square roots in around 2.5 seconds
on a CPU with 2.2 GHz, so around 550 cycles per square root,
on a Zen 1 (I have not bothered to filter out the overhead).

On Tuesday, March 22, 2022 at 10:47:31 PM UTC+2, Thomas Koenig wrote:
> MitchAlsup <Mitch...@aol.com> schrieb:
> > On Tuesday, March 22, 2022 at 12:26:03 PM UTC-5, Thomas Koenig wrote:
> >> MitchAlsup <Mitch...@aol.com> schrieb:
> >> > My feelings are::
> >> > a) screw the approximations
> >> > b) because doing them to faithful accuracy is easy
> >> > c) the right answer is only 14 cycles
> >> > d) if you build the FMAC unit correctly
> >> > See: USPTO 10,761,806
> ><
> >> Impressive (as I said before). This being a patent means that
> >> you forbid everybody to do what you claim, except outside the US
> >> (I believe Marcus is based in Sweden) and, of course, research is
> >> OK too. It is unclear if an FPGA-based implementation would fall
> >> under this patent ("processor" does not seem to be defined, and as
> >> I have learnt the hard way, you need to define things in patents).
> ><
> > I am willing to sell licenses for less money that it takes you to
> > assign a design engineer and have him/her read the document
> > thoroughly.
> I'm not going to pay for this, my profession is something
> quite different :-)
> > What is claimed is the necessary organization. An FPGA with
> > the necessary organization would infringe.
> Hm, if I had been involved in writing that patent, I would probably
> have included a description of what "processor" actually means
> in this context (which is why stuff I (co)write tends to be much
> longer, possibly unnecessarily so; it is, however, in a field
> totally unrelated to computer architecture).
>
> But it is of course your patent and your style of writing :-)
> >> Regarding approximate reciprocals: I would prefer a CPU which has
> >> different widths of floating point data (16, 32, 64 and possibly 128
> >> bits), if the smaller types are faster and more energy efficient
> >> than the larger types. People who need a low-cost, low-accuracy
> >> sqrt could then use the 16-bit version.
> ><
> > I am nor arguing against that point. From GPU work, the single precision
> > transcendentals take 5 cycles for 1 ULP accuracy (not quite faithful),
> > 1 cycle throughput. Due to the organization of SP trans, 16-bit takes
> > just as long as 32-bit <in GPUs>.
> ><
> > In order to get this working at the 128-bit level, you are going to need
> > FP with 128-bits of fraction; just so you can correctly calculate the
> > coefficients the polynomial is based upon. My estimate is that 128-bit
> > transcendentals will cost ~40 cycles.
> Certainly faster than doing the same thing in software.
>
> The little test program
>
> 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)
> do i=1,m
> a = sqrt(a)
> read (unit=c,fmt=*) k
> write (*,*) a(k)
> end do
> end program main
>
> calculates 10**7 16-byte square roots in around 2.5 seconds
> on a CPU with 2.2 GHz, so around 550 cycles per square root,
> on a Zen 1 (I have not bothered to filter out the overhead).

If I understand correctly, real(kind=qp) is equivalent of IEEE binary128.
IIRC, my own quad-precision class calculates square roots 2-3 times faster.
And it has better precision that gnu binary128 - 128bit mantissa instead of 113 bits.
So, it seems, for binary128 I should be able to calculate sqrt a little faster yet.
Of course, not in 40 clock cycles, but may be, in 105-110.
Or may be not. My internal format is easier to parse in sw than binary128
so, possibly more complicated parsing will bring the speed back to the same point.

Besides, does not your test end up calculating square root of exactly 1 for majority of loop iterations?
Since we're talking about sw, there is a danger that this case is not representative of typical timing.

On Tuesday, March 22, 2022 at 7:24:58 AM UTC-6, Terje Mathisen wrote:

> Even when you need exact results, having that reciprocal starting point
> means that you will converge on that value significantly faster.

Of course, but if the only use of an approximate reciprocal is to
begin a software division routine... why not just have a divide
instruction instead? One that is fast and efficient.

John Savard

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?

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

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

Thomas Koenig wrote:
> MitchAlsup <MitchAlsup@aol.com> schrieb:
>> On Tuesday, March 22, 2022 at 12:26:03 PM UTC-5, Thomas Koenig wrote:
>>> MitchAlsup <Mitch...@aol.com> schrieb:
>>>> My feelings are::
>>>> a) screw the approximations
>>>> b) because doing them to faithful accuracy is easy
>>>> c) the right answer is only 14 cycles
>>>> d) if you build the FMAC unit correctly
>>>> See: USPTO 10,761,806
>> <
>>> Impressive (as I said before). This being a patent means that
>>> you forbid everybody to do what you claim, except outside the US
>>> (I believe Marcus is based in Sweden) and, of course, research is
>>> OK too. It is unclear if an FPGA-based implementation would fall
>>> under this patent ("processor" does not seem to be defined, and as
>>> I have learnt the hard way, you need to define things in patents).
>> <
>> I am willing to sell licenses for less money that it takes you to
>> assign a design engineer and have him/her read the document
>> thoroughly.
>
> I'm not going to pay for this, my profession is something
> quite different :-)
>
>
>> What is claimed is the necessary organization. An FPGA with
>> the necessary organization would infringe.
>
> Hm, if I had been involved in writing that patent, I would probably
> have included a description of what "processor" actually means
> in this context (which is why stuff I (co)write tends to be much
> longer, possibly unnecessarily so; it is, however, in a field
> totally unrelated to computer architecture).
>
> But it is of course your patent and your style of writing :-)
>
>>> Regarding approximate reciprocals: I would prefer a CPU which has
>>> different widths of floating point data (16, 32, 64 and possibly 128
>>> bits), if the smaller types are faster and more energy efficient
>>> than the larger types. People who need a low-cost, low-accuracy
>>> sqrt could then use the 16-bit version.
>> <
>> I am nor arguing against that point. From GPU work, the single precision
>> transcendentals take 5 cycles for 1 ULP accuracy (not quite faithful),
>> 1 cycle throughput. Due to the organization of SP trans, 16-bit takes
>> just as long as 32-bit <in GPUs>.
>> <
>> In order to get this working at the 128-bit level, you are going to need
>> FP with 128-bits of fraction; just so you can correctly calculate the
>> coefficients the polynomial is based upon. My estimate is that 128-bit
>> transcendentals will cost ~40 cycles.
>
> Certainly faster than doing the same thing in software.
>
> The little test program
>
> 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)
> do i=1,m
> a = sqrt(a)
> read (unit=c,fmt=*) k
> write (*,*) a(k)
> end do
> end program main
>
> calculates 10**7 16-byte square roots in around 2.5 seconds
> on a CPU with 2.2 GHz, so around 550 cycles per square root,
> on a Zen 1 (I have not bothered to filter out the overhead).
>
If those random inputs cover the full 128-bit range, then I'm quite
impressed, 550 cycles really isn't that bad:

I would probably implement sqrtq() by taking the input, extracting the
top 53 mantissa bits and the bottom exponent bit (plus a 1022 or 1023
bias), then take the sqrt() of that. Two NR iterations with 128-bit and
256-bit precision would get _very_ close to perfect rounding.

OTOH, it might make more sense to start with the reciprocal sqrt lookup
(assume at least 9 bits), then do the much simpler InvSqrt NR with
increasing precision, where only the last two iterations would need
extended precision, and the final one is modified to produce sqrt
instead of InvSqrt. To get perfect rounding for the final 112-bit
mantissa means that we need almost 230 bits for that last stage, making
it very hard to do with just double precision variables (needing 3
double variables), while having a 64x64->128 unsigned int multiplier
makes it all work out nicely.

A 128x128->256 MUL needs 4 64x64->128 operations plus a bunch of
carries, I'm guessing it could be done in less than 20 cycles?

;; Inputs in r10:r9 and r12:r11, result in r16:r15:r14:r13
mov rax,r9
mul r11
mov r13,rax
mov r14,rdx

mov rax,r10
mul r11
xor r15,r15
add r14,rax
adc r15,rdx

mov rax,r9
mul r12
xor r16,r16
add r14,rax
adc r15,rdx
adc r16,0

mov rax,r10
mul r12
add r15,rax
adc r16,rdx

With a wide OoO cpu those MULs are all independent so they can overlap,
only the carries are sequential, so in theory it could run in close to
10 cycles using a 4-cycle MUL as the building block, but 20 seems safe.

Terje

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

Quadibloc wrote:
> On Tuesday, March 22, 2022 at 7:24:58 AM UTC-6, Terje Mathisen wrote:
>
>> Even when you need exact results, having that reciprocal starting point
>> means that you will converge on that value significantly faster.
>
> Of course, but if the only use of an approximate reciprocal is to
> begin a software division routine... why not just have a divide
> instruction instead? One that is fast and efficient.

I was more worried about InvSqrt(), which can of course be used as a
building block for both FSQRT and FDIV as well.

Terje

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

Michael S <already5chosen@yahoo.com> schrieb:

>> The little test program
>>
>> 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)
>> do i=1,m
>> a = sqrt(a)
>> read (unit=c,fmt=*) k
>> write (*,*) a(k)
>> end do
>> end program main
>>
>> calculates 10**7 16-byte square roots in around 2.5 seconds
>> on a CPU with 2.2 GHz, so around 550 cycles per square root,
>> on a Zen 1 (I have not bothered to filter out the overhead).
>
> If I understand correctly, real(kind=qp) is equivalent of IEEE binary128.

This is what GNU Fortran uses on most platforms, using the quadmath
library.

> IIRC, my own quad-precision class calculates square roots 2-3 times faster.

That would be interesting. Could you post a benchmarks?

> And it has better precision that gnu binary128 - 128bit mantissa instead of 113 bits.
> So, it seems, for binary128 I should be able to calculate sqrt a little faster yet.
> Of course, not in 40 clock cycles, but may be, in 105-110.
> Or may be not. My internal format is easier to parse in sw than binary128
> so, possibly more complicated parsing will bring the speed back to the same point.
>
> Besides, does not your test end up calculating square root of
> exactly 1 for majority of loop iterations? Since we're talking
> about sw, there is a danger that this case is not representative
> of typical timing.

Not exactly one, but close to one.

Here's an updated program:

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

whose output is

7671809140563249423302811.97887168825
2769803086965.43433961116406142440907
1664272.53987002812039255042083903983
1290.06687418522150052441522368030520
35.9175009457119872785306931352753634
5.99312113557802094243808280085584780
2.44808519777764698252724820741238982
1.56463580355865785822592147498455777
1.25085402967678760971072508986042575
1.11841585721805089654732018717047560
Used 2.37363 seconds to calculate 1.00000E+07 square roots.

so the number of cycles is similar.

I just ran the same program on a POWER with hardware
support for 128-bit floats (using to-be-released support for
-mabi=ieeelongdouble), and the output was

8625926447798841161143265.26897396334
2936992755830.16056596538977565335987
1713765.66537848475066226136266112131
1309.10872939511206616748020241410596
36.1816076120881072462540037278106935
6.01511492925015583581640205312340543
2.45257312413924894901462105871990474
1.56606932290344316961162320054757086
1.25142691472712187692445496527231813
1.11867194240631684707050674605147131
Used 0.31389 seconds to calculate 1.00000E+07 square roots

This machine also uses a 2.2 GHz clock (or claims to - hard to
be sure with a virtual machine, it is likely to lie about its CPU
frequency), so that would come to around 70 cycles including the
loop and load/store overhead.

On Wednesday, March 23, 2022 at 7:54:38 PM UTC+2, Thomas Koenig wrote:
> Michael S <already...@yahoo.com> schrieb:
> >> The little test program
> >>
> >> 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)
> >> do i=1,m
> >> a = sqrt(a)
> >> read (unit=c,fmt=*) k
> >> write (*,*) a(k)
> >> end do
> >> end program main
> >>
> >> calculates 10**7 16-byte square roots in around 2.5 seconds
> >> on a CPU with 2.2 GHz, so around 550 cycles per square root,
> >> on a Zen 1 (I have not bothered to filter out the overhead).
> >
> > If I understand correctly, real(kind=qp) is equivalent of IEEE binary128.
> This is what GNU Fortran uses on most platforms, using the quadmath
> library.
> > IIRC, my own quad-precision class calculates square roots 2-3 times faster.
> That would be interesting. Could you post a benchmarks?

The benchmark is a part of test suite for my quad-precision class.
It is call, unsurprisingly, tst_sqrt.
https://github.com/already5chosen/extfloat/tree/master/128
The code was developed on Windows and tested both with MSVC and gcc (under msys64).
It was never built or tested on Linux or on platforms other than x64, but I expect that it will take
less than hour to make it running on aarch64-linux/gcc or clang or on little-endian POWER-Linux.
The biggest expected difficulty would be in test infrastructure. My timing report uses __rdtsc() that
will have to be replaced.
Big-endian POWER - I am less sure. It should work, but probably not on the first try.

The speed is reported in nominal __rdtsc() cycles so it underestimates the # of actual CPU cycles.
Less so on desktops, more so on laptops or multicore servers.

And it is faster than what I remembered.
On core i5-3450 (Ivy Bridge, nominal freq=3100) it reports 78 cycles.
On Xeon E3-1271 v3 (Haswell, nominal freq=3600) it reports 72 cycles.
On Xeon E-2176G (Skylake, nominal freq=3700) it reports 58 cycles.

For the last machine I know that the actual frequency=4250, so actual # of cycles = 67.
For the other two I am not sure would guess that i5 is 3400-3500 and Xeon-E3 is at max, turbo=4000.
So, actual # of cycles is, respectively, 86-88 and 80 cycles.
Zen1 is very similar to Haswell, so I expect ~80 cycles, too.

> > And it has better precision that gnu binary128 - 128bit mantissa instead of 113 bits.
> > So, it seems, for binary128 I should be able to calculate sqrt a little faster yet.
> > Of course, not in 40 clock cycles, but may be, in 105-110.
> > Or may be not. My internal format is easier to parse in sw than binary128
> > so, possibly more complicated parsing will bring the speed back to the same point.
> >
> > Besides, does not your test end up calculating square root of
> > exactly 1 for majority of loop iterations? Since we're talking
> > about sw, there is a danger that this case is not representative
> > of typical timing.
> Not exactly one, but close to one.
>
> Here's an updated program:
> 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
>
> whose output is
>
> 7671809140563249423302811.97887168825
> 2769803086965.43433961116406142440907
> 1664272.53987002812039255042083903983
> 1290.06687418522150052441522368030520
> 35.9175009457119872785306931352753634
> 5.99312113557802094243808280085584780
> 2.44808519777764698252724820741238982
> 1.56463580355865785822592147498455777
> 1.25085402967678760971072508986042575
> 1.11841585721805089654732018717047560
> Used 2.37363 seconds to calculate 1.00000E+07 square roots.
>
> so the number of cycles is similar.
>
> I just ran the same program on a POWER with hardware
> support for 128-bit floats (using to-be-released support for
> -mabi=ieeelongdouble), and the output was
>
> 8625926447798841161143265.26897396334
> 2936992755830.16056596538977565335987
> 1713765.66537848475066226136266112131
> 1309.10872939511206616748020241410596
> 36.1816076120881072462540037278106935
> 6.01511492925015583581640205312340543
> 2.45257312413924894901462105871990474
> 1.56606932290344316961162320054757086
> 1.25142691472712187692445496527231813
> 1.11867194240631684707050674605147131
> Used 0.31389 seconds to calculate 1.00000E+07 square roots
>
> This machine also uses a 2.2 GHz clock (or claims to - hard to
> be sure with a virtual machine, it is likely to lie about its CPU
> frequency), so that would come to around 70 cycles including the
> loop and load/store overhead.

POWER with 128-bit FP hardware running at 2.2 GHz?
It seems to me, such HW does not exist.
The last POWER CPU that was able to run that so slowly was POWER5+, but it had no quad-precision HW.
If I am not mistaken, for POWER9 a minimal clock is 2.8 GHz, for POWER6/7/8 it's higher.

Michael S <already5chosen@yahoo.com> schrieb:
> On Wednesday, March 23, 2022 at 7:54:38 PM UTC+2, Thomas Koenig wrote:
>> Michael S <already...@yahoo.com> schrieb:
>> >> The little test program
>> >>
>> >> 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)
>> >> do i=1,m
>> >> a = sqrt(a)
>> >> read (unit=c,fmt=*) k
>> >> write (*,*) a(k)
>> >> end do
>> >> end program main
>> >>
>> >> calculates 10**7 16-byte square roots in around 2.5 seconds
>> >> on a CPU with 2.2 GHz, so around 550 cycles per square root,
>> >> on a Zen 1 (I have not bothered to filter out the overhead).
>> >
>> > If I understand correctly, real(kind=qp) is equivalent of IEEE binary128.
>> This is what GNU Fortran uses on most platforms, using the quadmath
>> library.
>> > IIRC, my own quad-precision class calculates square roots 2-3 times faster.
>> That would be interesting. Could you post a benchmarks?
>
> The benchmark is a part of test suite for my quad-precision class.
> It is call, unsurprisingly, tst_sqrt.
> https://github.com/already5chosen/extfloat/tree/master/128
> The code was developed on Windows and tested both with MSVC and gcc (under msys64).
> It was never built or tested on Linux or on platforms other than x64, but I expect that it will take
> less than hour to make it running on aarch64-linux/gcc or clang or on little-endian POWER-Linux.
> The biggest expected difficulty would be in test infrastructure. My timing report uses __rdtsc() that
> will have to be replaced.

Wall-clock timing is probably better because rdtsc is notoriously
unreliable.

> Big-endian POWER - I am less sure. It should work, but probably not on the first try.
>
> The speed is reported in nominal __rdtsc() cycles so it underestimates the # of actual CPU cycles.
> Less so on desktops, more so on laptops or multicore servers.
>
> And it is faster than what I remembered.
> On core i5-3450 (Ivy Bridge, nominal freq=3100) it reports 78 cycles.
> On Xeon E3-1271 v3 (Haswell, nominal freq=3600) it reports 72 cycles.
> On Xeon E-2176G (Skylake, nominal freq=3700) it reports 58 cycles.
>
> For the last machine I know that the actual frequency=4250, so actual # of cycles = 67.
> For the other two I am not sure would guess that i5 is 3400-3500 and Xeon-E3 is at max, turbo=4000.
> So, actual # of cycles is, respectively, 86-88 and 80 cycles.
> Zen1 is very similar to Haswell, so I expect ~80 cycles, too.
>
>> > And it has better precision that gnu binary128 - 128bit mantissa instead of 113 bits.
>> > So, it seems, for binary128 I should be able to calculate sqrt a little faster yet.
>> > Of course, not in 40 clock cycles, but may be, in 105-110.
>> > Or may be not. My internal format is easier to parse in sw than binary128
>> > so, possibly more complicated parsing will bring the speed back to the same point.
>> >
>> > Besides, does not your test end up calculating square root of
>> > exactly 1 for majority of loop iterations? Since we're talking
>> > about sw, there is a danger that this case is not representative
>> > of typical timing.
>> Not exactly one, but close to one.
>>
>> Here's an updated program:
>> 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
>>
>> whose output is
>>
>> 7671809140563249423302811.97887168825
>> 2769803086965.43433961116406142440907
>> 1664272.53987002812039255042083903983
>> 1290.06687418522150052441522368030520
>> 35.9175009457119872785306931352753634
>> 5.99312113557802094243808280085584780
>> 2.44808519777764698252724820741238982
>> 1.56463580355865785822592147498455777
>> 1.25085402967678760971072508986042575
>> 1.11841585721805089654732018717047560
>> Used 2.37363 seconds to calculate 1.00000E+07 square roots.
>>
>> so the number of cycles is similar.
>>
>> I just ran the same program on a POWER with hardware
>> support for 128-bit floats (using to-be-released support for
>> -mabi=ieeelongdouble), and the output was
>>
>> 8625926447798841161143265.26897396334
>> 2936992755830.16056596538977565335987
>> 1713765.66537848475066226136266112131
>> 1309.10872939511206616748020241410596
>> 36.1816076120881072462540037278106935
>> 6.01511492925015583581640205312340543
>> 2.45257312413924894901462105871990474
>> 1.56606932290344316961162320054757086
>> 1.25142691472712187692445496527231813
>> 1.11867194240631684707050674605147131
>> Used 0.31389 seconds to calculate 1.00000E+07 square roots
>>
>> This machine also uses a 2.2 GHz clock (or claims to - hard to
>> be sure with a virtual machine, it is likely to lie about its CPU
>> frequency), so that would come to around 70 cycles including the
>> loop and load/store overhead.
>
> POWER with 128-bit FP hardware running at 2.2 GHz?
> It seems to me, such HW does not exist.

You can always tune down the frequency...

$ head /proc/cpuinfo
processor : 0
cpu : POWER9 (architected), altivec supported
clock : 2200.000000MHz
revision : 2.2 (pvr 004e 1202)

processor : 1
cpu : POWER9 (architected), altivec supported
clock : 2200.000000MHz
revision : 2.2 (pvr 004e 1202)

> The last POWER CPU that was able to run that so slowly was POWER5+, but it had no quad-precision HW.
> If I am not mistaken, for POWER9 a minimal clock is 2.8 GHz, for POWER6/7/8 it's higher.

or, like I said, it may be lying about its core frequency
because it is running under a hypervisor.

On Thursday, March 24, 2022 at 11:46:23 AM UTC+2, Thomas Koenig wrote:
> Michael S <already...@yahoo.com> schrieb:
> > On Wednesday, March 23, 2022 at 7:54:38 PM UTC+2, Thomas Koenig wrote:
> >> Michael S <already...@yahoo.com> schrieb:
> >> >> The little test program
> >> >>
> >> >> 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)
> >> >> do i=1,m
> >> >> a = sqrt(a)
> >> >> read (unit=c,fmt=*) k
> >> >> write (*,*) a(k)
> >> >> end do
> >> >> end program main
> >> >>
> >> >> calculates 10**7 16-byte square roots in around 2.5 seconds
> >> >> on a CPU with 2.2 GHz, so around 550 cycles per square root,
> >> >> on a Zen 1 (I have not bothered to filter out the overhead).
> >> >
> >> > If I understand correctly, real(kind=qp) is equivalent of IEEE binary128.
> >> This is what GNU Fortran uses on most platforms, using the quadmath
> >> library.
> >> > IIRC, my own quad-precision class calculates square roots 2-3 times faster.
> >> That would be interesting. Could you post a benchmarks?
> >
> > The benchmark is a part of test suite for my quad-precision class.
> > It is call, unsurprisingly, tst_sqrt.
> > https://github.com/already5chosen/extfloat/tree/master/128
> > The code was developed on Windows and tested both with MSVC and gcc (under msys64).
> > It was never built or tested on Linux or on platforms other than x64, but I expect that it will take
> > less than hour to make it running on aarch64-linux/gcc or clang or on little-endian POWER-Linux.
> > The biggest expected difficulty would be in test infrastructure. My timing report uses __rdtsc() that
> > will have to be replaced.
> Wall-clock timing is probably better because rdtsc is notoriously

If all you want to know is wall clock then [on relative modern x86] rdtsc is very reliable.
At least when running on physical machines, as I did in all my measurements.
The problem with it is not reliability, but portability to non-x86 architectures.

> unreliable.
> > Big-endian POWER - I am less sure. It should work, but probably not on the first try.
> >
> > The speed is reported in nominal __rdtsc() cycles so it underestimates the # of actual CPU cycles.
> > Less so on desktops, more so on laptops or multicore servers.
> >
> > And it is faster than what I remembered.
> > On core i5-3450 (Ivy Bridge, nominal freq=3100) it reports 78 cycles.
> > On Xeon E3-1271 v3 (Haswell, nominal freq=3600) it reports 72 cycles.
> > On Xeon E-2176G (Skylake, nominal freq=3700) it reports 58 cycles.
> >
> > For the last machine I know that the actual frequency=4250, so actual # of cycles = 67.
> > For the other two I am not sure would guess that i5 is 3400-3500 and Xeon-E3 is at max, turbo=4000.
> > So, actual # of cycles is, respectively, 86-88 and 80 cycles.
> > Zen1 is very similar to Haswell, so I expect ~80 cycles, too.
> >
> >> > And it has better precision that gnu binary128 - 128bit mantissa instead of 113 bits.
> >> > So, it seems, for binary128 I should be able to calculate sqrt a little faster yet.
> >> > Of course, not in 40 clock cycles, but may be, in 105-110.
> >> > Or may be not. My internal format is easier to parse in sw than binary128
> >> > so, possibly more complicated parsing will bring the speed back to the same point.
> >> >
> >> > Besides, does not your test end up calculating square root of
> >> > exactly 1 for majority of loop iterations? Since we're talking
> >> > about sw, there is a danger that this case is not representative
> >> > of typical timing.
> >> Not exactly one, but close to one.
> >>
> >> Here's an updated program:
> >> 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
> >>
> >> whose output is
> >>
> >> 7671809140563249423302811.97887168825
> >> 2769803086965.43433961116406142440907
> >> 1664272.53987002812039255042083903983
> >> 1290.06687418522150052441522368030520
> >> 35.9175009457119872785306931352753634
> >> 5.99312113557802094243808280085584780
> >> 2.44808519777764698252724820741238982
> >> 1.56463580355865785822592147498455777
> >> 1.25085402967678760971072508986042575
> >> 1.11841585721805089654732018717047560
> >> Used 2.37363 seconds to calculate 1.00000E+07 square roots.
> >>
> >> so the number of cycles is similar.
> >>
> >> I just ran the same program on a POWER with hardware
> >> support for 128-bit floats (using to-be-released support for
> >> -mabi=ieeelongdouble), and the output was
> >>
> >> 8625926447798841161143265.26897396334
> >> 2936992755830.16056596538977565335987
> >> 1713765.66537848475066226136266112131
> >> 1309.10872939511206616748020241410596
> >> 36.1816076120881072462540037278106935
> >> 6.01511492925015583581640205312340543
> >> 2.45257312413924894901462105871990474
> >> 1.56606932290344316961162320054757086
> >> 1.25142691472712187692445496527231813
> >> 1.11867194240631684707050674605147131
> >> Used 0.31389 seconds to calculate 1.00000E+07 square roots
> >>
> >> This machine also uses a 2.2 GHz clock (or claims to - hard to
> >> be sure with a virtual machine, it is likely to lie about its CPU
> >> frequency), so that would come to around 70 cycles including the
> >> loop and load/store overhead.
> >
> > POWER with 128-bit FP hardware running at 2.2 GHz?
> > It seems to me, such HW does not exist.
> You can always tune down the frequency...
>
> $ head /proc/cpuinfo
> processor : 0
> cpu : POWER9 (architected), altivec supported
> clock : 2200.000000MHz
> revision : 2.2 (pvr 004e 1202)
>
> processor : 1
> cpu : POWER9 (architected), altivec supported
> clock : 2200.000000MHz
> revision : 2.2 (pvr 004e 1202)
> > The last POWER CPU that was able to run that so slowly was POWER5+, but it had no quad-precision HW.
> > If I am not mistaken, for POWER9 a minimal clock is 2.8 GHz, for POWER6/7/8 it's higher.
> or, like I said, it may be lying about its core frequency
> because it is running under a hypervisor.

Michael S <already5chosen@yahoo.com> schrieb:
> On Thursday, March 24, 2022 at 11:46:23 AM UTC+2, Thomas Koenig wrote:
>> Michael S <already...@yahoo.com> schrieb:
>> > On Wednesday, March 23, 2022 at 7:54:38 PM UTC+2, Thomas Koenig wrote:
>> >> Michael S <already...@yahoo.com> schrieb:
>> >> >> The little test program
>> >> >>
>> >> >> 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)
>> >> >> do i=1,m
>> >> >> a = sqrt(a)
>> >> >> read (unit=c,fmt=*) k
>> >> >> write (*,*) a(k)
>> >> >> end do
>> >> >> end program main
>> >> >>
>> >> >> calculates 10**7 16-byte square roots in around 2.5 seconds
>> >> >> on a CPU with 2.2 GHz, so around 550 cycles per square root,
>> >> >> on a Zen 1 (I have not bothered to filter out the overhead).
>> >> >
>> >> > If I understand correctly, real(kind=qp) is equivalent of IEEE binary128.
>> >> This is what GNU Fortran uses on most platforms, using the quadmath
>> >> library.
>> >> > IIRC, my own quad-precision class calculates square roots 2-3 times faster.
>> >> That would be interesting. Could you post a benchmarks?
>> >
>> > The benchmark is a part of test suite for my quad-precision class.
>> > It is call, unsurprisingly, tst_sqrt.
>> > https://github.com/already5chosen/extfloat/tree/master/128
>> > The code was developed on Windows and tested both with MSVC and gcc (under msys64).
>> > It was never built or tested on Linux or on platforms other than x64, but I expect that it will take
>> > less than hour to make it running on aarch64-linux/gcc or clang or on little-endian POWER-Linux.
>> > The biggest expected difficulty would be in test infrastructure. My timing report uses __rdtsc() that
>> > will have to be replaced.
>> Wall-clock timing is probably better because rdtsc is notoriously
>
> If all you want to know is wall clock then [on relative modern x86] rdtsc is very reliable.

So, what is the wall-clock timing for calculating 10^7 square roots in
128 bit precision, on what sort of system?

On Thursday, March 24, 2022 at 4:59:22 PM UTC+2, Thomas Koenig wrote:
> Michael S <already...@yahoo.com> schrieb:
> > On Thursday, March 24, 2022 at 11:46:23 AM UTC+2, Thomas Koenig wrote:
> >> Michael S <already...@yahoo.com> schrieb:
> >> > On Wednesday, March 23, 2022 at 7:54:38 PM UTC+2, Thomas Koenig wrote:
> >> >> Michael S <already...@yahoo.com> schrieb:
> >> >> >> The little test program
> >> >> >>
> >> >> >> 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)
> >> >> >> do i=1,m
> >> >> >> a = sqrt(a)
> >> >> >> read (unit=c,fmt=*) k
> >> >> >> write (*,*) a(k)
> >> >> >> end do
> >> >> >> end program main
> >> >> >>
> >> >> >> calculates 10**7 16-byte square roots in around 2.5 seconds
> >> >> >> on a CPU with 2.2 GHz, so around 550 cycles per square root,
> >> >> >> on a Zen 1 (I have not bothered to filter out the overhead).
> >> >> >
> >> >> > If I understand correctly, real(kind=qp) is equivalent of IEEE binary128.
> >> >> This is what GNU Fortran uses on most platforms, using the quadmath
> >> >> library.
> >> >> > IIRC, my own quad-precision class calculates square roots 2-3 times faster.
> >> >> That would be interesting. Could you post a benchmarks?
> >> >
> >> > The benchmark is a part of test suite for my quad-precision class.
> >> > It is call, unsurprisingly, tst_sqrt.
> >> > https://github.com/already5chosen/extfloat/tree/master/128
> >> > The code was developed on Windows and tested both with MSVC and gcc (under msys64).
> >> > It was never built or tested on Linux or on platforms other than x64, but I expect that it will take
> >> > less than hour to make it running on aarch64-linux/gcc or clang or on little-endian POWER-Linux.
> >> > The biggest expected difficulty would be in test infrastructure. My timing report uses __rdtsc() that
> >> > will have to be replaced.
> >> Wall-clock timing is probably better because rdtsc is notoriously
> >
> > If all you want to know is wall clock then [on relative modern x86] rdtsc is very reliable.
> So, what is the wall-clock timing for calculating 10^7 square roots in
> 128 bit precision, on what sort of system?

The needed numbers are in a post above.
Core i5-3450: 78/3100e6 *1e7 = 0.252s
Xeon E3-1271 v3: 72/3600e6*1e7= 0.200s
Xeon E-2176G: 58/3700e6*1e7 = 0.157s

On 3/22/2022 11:44 AM, MitchAlsup wrote:
> On Tuesday, March 22, 2022 at 10:56:06 AM UTC-5, Marcus wrote:
>> On 2022-03-22, MitchAlsup wrote:
>>> On Tuesday, March 22, 2022 at 3:25:11 AM UTC-5, 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.
>>> <
>>> My feelings are::
>>> a) screw the approximations
>>> b) because doing them to faithful accuracy is easy
>>> c) the right answer is only 14 cycles
>>> d) if you build the FMAC unit correctly
>>> See: USPTO 10,761,806
>> Thanks, that makes for a good read!
>>
>> While I see the benefits of fast argument reduction and a generalized
>> coefficient + FMA engine for implementing transcendental functions, the
>> reciprocal and reciprocal sqrt functions are simpler and common enough
>> to possibly warrant a dedicated solution?
> <
> You get to make those choices.
>>
>> E.g. the CRAY could do a full DP division by means of approximate
>> reciprocal + NR-step + multiply by nominator in just 3-4 clock cycles
>> (vectorized), which is pretty impressive even by today's standards.
> <
> Its OK if you are content with 46-bits of accuracy.
> <

A lot likely depends on the approximate reciprocal.

Say, for FRCP:
MagicConstant-FpBits

Accuracy can be improved slightly by using a lookup table for the
high-order bits, which gives a bias to add to the high-order bits of the
result.

One still needs a few N-R stages though, so I ended up skipping these as
they didn't buy much over doing it with integer ops. Hardware could
potentially do special cases though, like turning Zero into NaN or similar.

The N-R can be mapped to FMAC, though this case is helped if one has
either multiple versions (with built in sign negation), or the ability
to negate values cheaply.

Eg:
A*B+C
A*B-C
-A*B+C
-A*B-C

For graphics-processing and similar (with FP16), one can potentially
skip the N-R step.

So, say:
V2=V0/V1
Becomes:
V2=V0*(MagicBias_Div-V1)

Usually works because people are not overly fussy about exact values
when it comes to pixels.

Similar, goes for SQRT:
(Val>>1)+MagicBias_Sqrt
Or, InvSqrt:
MagicBias_InvSqrt-(Val>>1)

....

Likewise, for FP16 one could possibly get along OK using a single-stage
approximation as the final result in many cases.

>>>>
>>>> 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.
>>>>
>>>> /Marcus

The itanium had only an fp reciprocal; no division, that had to be done in
software. See https://www.cl.cam.ac.uk/~jrh13/papers/hol00.pdf

On Tue, 22 Mar 2022, 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.
>
> /Marcus
>

On Thursday, March 24, 2022 at 5:31:28 PM UTC+2, Michael S wrote:
> On Thursday, March 24, 2022 at 4:59:22 PM UTC+2, Thomas Koenig wrote:
> > Michael S <already...@yahoo.com> schrieb:
> > > On Thursday, March 24, 2022 at 11:46:23 AM UTC+2, Thomas Koenig wrote:
> > >> Michael S <already...@yahoo.com> schrieb:
> > >> > On Wednesday, March 23, 2022 at 7:54:38 PM UTC+2, Thomas Koenig wrote:
> > >> >> Michael S <already...@yahoo.com> schrieb:
> > >> >> >> The little test program
> > >> >> >>
> > >> >> >> 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)
> > >> >> >> do i=1,m
> > >> >> >> a = sqrt(a)
> > >> >> >> read (unit=c,fmt=*) k
> > >> >> >> write (*,*) a(k)
> > >> >> >> end do
> > >> >> >> end program main
> > >> >> >>
> > >> >> >> calculates 10**7 16-byte square roots in around 2.5 seconds
> > >> >> >> on a CPU with 2.2 GHz, so around 550 cycles per square root,
> > >> >> >> on a Zen 1 (I have not bothered to filter out the overhead).
> > >> >> >
> > >> >> > If I understand correctly, real(kind=qp) is equivalent of IEEE binary128.
> > >> >> This is what GNU Fortran uses on most platforms, using the quadmath
> > >> >> library.
> > >> >> > IIRC, my own quad-precision class calculates square roots 2-3 times faster.
> > >> >> That would be interesting. Could you post a benchmarks?
> > >> >
> > >> > The benchmark is a part of test suite for my quad-precision class.
> > >> > It is call, unsurprisingly, tst_sqrt.
> > >> > https://github.com/already5chosen/extfloat/tree/master/128
> > >> > The code was developed on Windows and tested both with MSVC and gcc (under msys64).
> > >> > It was never built or tested on Linux or on platforms other than x64, but I expect that it will take
> > >> > less than hour to make it running on aarch64-linux/gcc or clang or on little-endian POWER-Linux.
> > >> > The biggest expected difficulty would be in test infrastructure. My timing report uses __rdtsc() that
> > >> > will have to be replaced.
> > >> Wall-clock timing is probably better because rdtsc is notoriously
> > >
> > > If all you want to know is wall clock then [on relative modern x86] rdtsc is very reliable.
> > So, what is the wall-clock timing for calculating 10^7 square roots in
> > 128 bit precision, on what sort of system?
> The needed numbers are in a post above.
> Core i5-3450: 78/3100e6 *1e7 = 0.252s
> Xeon E3-1271 v3: 72/3600e6*1e7= 0.200s
> Xeon E-2176G: 58/3700e6*1e7 = 0.157s

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.
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.