Sign in Join

Live chat for Free Pascal, Lazarus, Delphi and every Pascal programmer.

Mailing lists / fpc-devel / 2018-04

Broken frac function in FPC3.1.1 / Windows x86_64

37 messages · 11 people · –

Reply by email Subscribe to this list Your subscription

Read-only archive. Use j and k to move between messages.

Unknown author

If you compile and run this 64-bit program on Win 64 you get a crash

program tfrac;
var
   x,y: double;
begin
   x := 1e20;
   writeln('x=',x);
   y := frac(x);
   writeln('frac(x)=',y);
end.

D:\Work\DAMath>D:\FPC311\bin\i386-win32\ppcrossx64.exe tfrac.pas Free Pascal Compiler version 3.1.1-r20:38794 [2018/04/22] for x86_64 Copyright (c) 1993-2018 by Florian Klaempfl and others Target OS: Win64 for x64
Compiling tfrac.pas
Linking tfrac.exe
9 lines compiled, 0.1 sec, 37472 bytes code, 2036 bytes data

D:\Work\DAMath>tfrac
x= 1.0000000000000000E+020
Runtime error 207 at $0000000100001BD0
$0000000100001BD0
$0000000100001614
$0000000100001696
$00000001000095D0
$0000000100001580
$00000000773D59CD
$000000007763383D

I recently posted an issue to the bugtracker

https://bugs.freepascal.org/view.php?id=33635

But the only responses are by a reporter named Thaddy de Koning, which are for me totally useless and IMO he did miss the point.

Are there changes to the frac function (or its documentation) in the newer 64-bit windows cross-compiler? Especially for the allowed ranges?

This would not be intuitive because the corresponding int function works as
expected, and meanwhile I use the work-around  frac(x) := int(x) -x

Wolfgang Ehrhardt

Reply by email

Thorsten Engler
to Unknown author
But the only responses are by a reporter named Thaddy de Koning, which are for me totally useless and IMO he did miss the point.

It's not the first time the person has been less than helpful...

After looking through other issues he has commented on, I don't think he has contributed anything beyond, usually wrong, often aggressive and/or insulting, notes on the tracker.

Reply by email

Mark Morgan Lloyd
to Thorsten Engler: It's not the first time the person has been less than helpful... After looking through…
Thorsten Engler wrote
> But the only responses are by a reporter named Thaddy de Koning,> which are for me totally useless and IMO he did miss the point.It's not the first time the person has been less than helpful... After looking through other issues he has commented on, I don't think he has contributed anything beyond, usually wrong, often aggressive and/or insulting, notes on the tracker.

Note that "reporter" means that somebody's not a manager of the bug tracker. The designation covers a lot of people who are regular contributors to- and testers of- the core FPC code, and it's probably unfair for somebody to accuse one of them of trolling without considering that.

Having said that, I think that the OP would possibly have made a bit more headway if he'd reported it as a simple regression against the 2016 version of the compiler. I have to say that I don't think Thaddy has exactly helped things by implying that a 64-bit double is converted into an 80-bit extended on an architecture which doesn't support that type.

Mark Morgan Lloyd markMLl .AT. telemetry.co .DOT. uk [Opinions above are the author's, not those of his employers or colleagues]

Reply by email

Thorsten Engler
Show quoted text 5 lines
Note that "reporter" means that somebody's not a manager of the bug tracker. The designation covers a lot of people who are regular contributors to- and testers of- the core FPC code, and it's probably unfair for somebody to accuse one of them of trolling without considering that.

"reporter" is simply default group that anyone who just created a new account on the bug tracker ends up with. It says absolutely nothing about who that person is or what he has done.

But if you check out the link Ondrej Pokorny's posted in the "Dangerous optimization in CASE..OF" thread about a week ago ( https://bugs.freepascal.org/view.php?id=33603 ) you will see that this isn't the first time Thaddy has been mentioned on this mailing list.

If you then search on the tracker for other issues that Thaddy has commented on, you will find that the majority of his notes don't seem particular helpful, but are often enough using quite aggressive and/or insulting language. (And he has been called out for it in comments by developers, e.g. here: https://bugs.freepascal.org/view.php?id=32911 ).

So, you know, whom the shoe fits...

Reply by email

Bart
to Unknown author: If you compile and run this 64-bit program on Win 64 you get a crash program tfrac; var…
If you compile and run this 64-bit program on Win 64 you get a crash

And AFAICS your analysis of the cause (see bugtracker) is correct as well.

function fpc_frac_real(d: ValReal) : ValReal;compilerproc; assembler; nostackframe;
asm
cvttsd2si %xmm0,%rax

    { Windows defines %xmm4 and %xmm5 as first non-parameter volatile registers;
      on SYSV systems all are considered as such, so use %xmm4 }
    cvtsi2sd %rax,%xmm4
    subsd %xmm4,%xmm0
  end;

CVTTSD2SI — Convert with Truncation Scalar Double-Precision Floating-Point Value to Signed Integer
This should not be used to get a ValReal result.

Bart

Reply by email

Sven Barth
to Bart

Bart schrieb am Fr., 27. Apr. 2018, 13:42:

Show quoted text 21 lines
On Wed, Apr 25, 2018 at 11:04 AM, wrote: > If you compile and run this 64-bit program on Win 64 you get a crash And AFAICS your analysis of the cause (see bugtracker) is correct as well. function fpc_frac_real(d: ValReal) : ValReal;compilerproc; assembler; nostackframe; asm cvttsd2si %xmm0,%rax { Windows defines %xmm4 and %xmm5 as first non-parameter volatile registers; on SYSV systems all are considered as such, so use %xmm4 } cvtsi2sd %rax,%xmm4 subsd %xmm4,%xmm0 end; CVTTSD2SI — Convert with Truncation Scalar Double-Precision Floating-Point Value to Signed Integer This should not be used to get a ValReal result.

The code essentially does the following (instruction by instruction):

=== code begin ===

tmpi := int64(d - trunc(d));
tmpd := double(tmpi);
Result := d - tmpd;

=== code end ===

Though why it fails with the given value is a different topic...

Regards,
Sven

Reply by email

Thorsten Engler
to Sven Barth

Highest integer that fits in a Int64:

9223372036854775808

1e20:

100000000000000000000

Your Int is overflowing.

You can’t implement Frac by going through an Integer, that will never work. Except if you have an integer that can hold 1.8E308 (which would be a 1024 bit integer, or thereabouts).

Earlier message 59 lines
From: fpc-devel On Behalf Of Sven Barth via fpc-devel Sent: Friday, 27 April 2018 23:47 To: FPC developers' list Cc: Sven Barth Subject: *** GMX Spamverdacht *** Re: [fpc-devel] Broken frac function in FPC3.1.1 / Windows x86_64 Bart < > schrieb am Fr., 27. Apr. 2018, 13:42: On Wed, Apr 25, 2018 at 11:04 AM, < > wrote: > If you compile and run this 64-bit program on Win 64 you get a crash And AFAICS your analysis of the cause (see bugtracker) is correct as well. function fpc_frac_real(d: ValReal) : ValReal;compilerproc; assembler; nostackframe; asm cvttsd2si %xmm0,%rax { Windows defines %xmm4 and %xmm5 as first non-parameter volatile registers; on SYSV systems all are considered as such, so use %xmm4 } cvtsi2sd %rax,%xmm4 subsd %xmm4,%xmm0 end; CVTTSD2SI — Convert with Truncation Scalar Double-Precision Floating-Point Value to Signed Integer This should not be used to get a ValReal result. The code essentially does the following (instruction by instruction): === code begin === tmpi := int64(d - trunc(d)); tmpd := double(tmpi); Result := d - tmpd; === code end === Though why it fails with the given value is a different topic... Regards, Sven

Reply by email

Sven Barth
to Thorsten Engler

Thorsten Engler schrieb am Fr., 27. Apr. 2018,
16:09:

Show quoted text 18 lines
Highest integer that fits in a Int64: 9223372036854775808 1e20: 100000000000000000000 Your Int is overflowing. You can’t implement Frac by going through an Integer, that will never work. Except if you have an integer that can hold 1.8E308 (which would be a 1024 bit integer, or thereabouts)

Yes, I saw that now as well, though it's even worse as the cvttsd2si instruction in fact only works with 32-bit integers. That additionally means that Trunc() and Round() are broken for such values as well as they rely on the same kind of functions. Int() works because it doesn't have an SSE version.

Regards,
Sven

Reply by email

Sergei Gorelkin
to Sven Barth
Sven Barth wrote 26 lines
Thorsten Engler < > schrieb am Fr., 27. Apr. 2018, 16:09: Highest integer that fits in a Int64:____ 9223372036854775808____ 1e20:____ 100000000000000000000____ __ __ Your Int is overflowing.____ __ __ You can’t implement Frac by going through an Integer, that will never work. Except if you have an integer that can hold 1.8E308 (which would be a 1024 bit integer, or thereabouts) Yes, I saw that now as well, though it's even worse as the cvttsd2si instruction in fact only works with 32-bit integers. That additionally means that Trunc() and Round() are broken for such values as well as they rely on the same kind of functions. Int() works because it doesn't have an SSE version.

That's true for i386. But on x86_64 cvt(t)sd2si instuctions can deal with int64 range, if destination register is a 64-bit one.

Regards,
Sergei

Reply by email

Sven Barth
to Thorsten Engler

Thorsten Engler schrieb am Fr., 27. Apr. 2018,
17:47:

Show quoted text 5 lines
> That's true for i386. But on x86_64 cvt(t)sd2si instuctions can > deal with int64 range, if destination register is a 64-bit one. You are still going to be at least 960-bit short...

I've disabled the SSE variant for now again till we've decided how we want to proceed with this.

Regards,
Sven

Reply by email

Thorsten Engler

For what it’s worth, Delphi simply decided to give up on doing it correctly and silently fail if the double is too large to fit in an Int64.

WriteLn(Frac(1e15+0.5));
WriteLn(Frac(1e16+0.5));

When executed in 32bit code, returns:

5.00000000000000E-0001

5.00000000000000E-0001

And when executed in 64bit code, returns:

5.00000000000000E-0001

0.00000000000000E+0000

Whoever had the *brilliant* idea to deprecate the FPU in 64bit mode should be shoot and quartered, not necessarily in that order. 64bit code is basically useless if you want to do anything meaningful with floats that are larger than fits into an int64.

Earlier message 24 lines
From: fpc-devel On Behalf Of Sven Barth via fpc-devel Sent: Saturday, 28 April 2018 01:51 To: FPC developers' list Cc: Sven Barth Subject: *** GMX Spamverdacht *** Re: [fpc-devel] *** GMX Spamverdacht *** Re: Broken frac function in FPC3.1.1 / Windows x86_64 Thorsten Engler < > schrieb am Fr., 27. Apr. 2018, 17:47: > That's true for i386. But on x86_64 cvt(t)sd2si instuctions can > deal with int64 range, if destination register is a 64-bit one. You are still going to be at least 960-bit short... I've disabled the SSE variant for now again till we've decided how we want to proceed with this. Regards, Sven

Reply by email

Sven Barth
to Thorsten Engler

Thorsten Engler schrieb am Fr., 27. Apr. 2018,
18:48:

Show quoted text 36 lines
For what it’s worth, Delphi simply decided to give up on doing it correctly and silently fail if the double is too large to fit in an Int64. WriteLn(Frac(1e15+0.5)); WriteLn(Frac(1e16+0.5)); When executed in 32bit code, returns: 5.00000000000000E-0001 5.00000000000000E-0001 And when executed in 64bit code, returns: 5.00000000000000E-0001 0.00000000000000E+0000 Whoever had the **brilliant** idea to deprecate the FPU in 64bit mode should be shoot and quartered, not necessarily in that order. 64bit code is basically useless if you want to do anything meaningful with floats that are larger than fits into an int64.

Ehm, what does the depreciation of the FPU have to do with that the floats are too large for an Int64? Extended would only allow for even larger values that would not fit. And all other functionality "simply" needs to be implemented correctly using the "new" SSE functionality, which however is not without its pitfalls as we have seen.

Regards,
Sven

Reply by email

Mattias Gaertner
to Thorsten Engler: For what it’s worth, Delphi simply decided to give up on doing it correctly and…
Show quoted text 7 lines
For what it’s worth, Delphi simply decided to give up on doing it correctly and silently fail if the double is too large to fit in an Int64. [...] WriteLn(Frac(1e15+0.5)); WriteLn(Frac(1e16+0.5)); When executed in 32bit code, returns:

That means it has the Extended 80bit type, so it can handle the 1e16+0.5.
http://docwiki.embarcadero.com/Libraries/Tokyo/en/System.Extended

Show quoted text 7 lines
5.00000000000000E-0001 5.00000000000000E-0001 And when executed in 64bit code, returns:

That means it has only the 64bit double type, so the +0.5 is lost.

5.00000000000000E-0001 0.00000000000000E+0000

Mattias

Reply by email

Thorsten Engler
to Sven Barth: Thorsten Engler schrieb am Fr., 27. Apr. 2018, 18:48: Ehm, what does the depreciation of…

Actually, that test is wrong, because the 32bit code will use Extended float literals instead of Double. If you force it to use double in 32bit you get the same result.

If you want to use SSE instructions, you will need to add additional checks to see if the value will fit into an integer before you use the SSE instructions.

Earlier message 51 lines
From: fpc-devel On Behalf Of Sven Barth via fpc-devel Sent: Saturday, 28 April 2018 03:33 To: FPC developers' list Cc: Sven Barth Subject: *** GMX Spamverdacht *** Re: [fpc-devel] Broken frac function in FPC3.1.1 / Windows x86_64 Thorsten Engler < > schrieb am Fr., 27. Apr. 2018, 18:48: For what it’s worth, Delphi simply decided to give up on doing it correctly and silently fail if the double is too large to fit in an Int64. WriteLn(Frac(1e15+0.5)); WriteLn(Frac(1e16+0.5)); When executed in 32bit code, returns: 5.00000000000000E-0001 5.00000000000000E-0001 And when executed in 64bit code, returns: 5.00000000000000E-0001 0.00000000000000E+0000 Whoever had the *brilliant* idea to deprecate the FPU in 64bit mode should be shoot and quartered, not necessarily in that order. 64bit code is basically useless if you want to do anything meaningful with floats that are larger than fits into an int64. Ehm, what does the depreciation of the FPU have to do with that the floats are too large for an Int64? Extended would only allow for even larger values that would not fit. And all other functionality "simply" needs to be implemented correctly using the "new" SSE functionality, which however is not without its pitfalls as we have seen. Regards, Sven

Reply by email

Thorsten Engler
to Mattias Gaertner: On Sat, 28 Apr 2018 02:48:14 +1000 That means it has the Extended 80bit type, so it can…

Yeah, I realized that myself and made a post about that a few minutes ago which seems to have crossed ways with yours...

I was just projecting my annoyance about the lack of precision when being forced to do math in double instead of extended...

Show quoted text 45 lines
-----Original Message----- From: fpc-devel On Behalf Of Mattias Gaertner Sent: Saturday, 28 April 2018 03:39 To: Subject: *** GMX Spamverdacht *** Re: [fpc-devel] Broken frac function in FPC3.1.1 / Windows x86_64 On Sat, 28 Apr 2018 02:48:14 +1000 "Thorsten Engler" wrote: > For what it’s worth, Delphi simply decided to give up on doing it correctly and silently fail if the double is too large to fit in an Int64. >[...] > WriteLn(Frac(1e15+0.5)); > > WriteLn(Frac(1e16+0.5)); > > When executed in 32bit code, returns: That means it has the Extended 80bit type, so it can handle the 1e16+0.5. http://docwiki.embarcadero.com/Libraries/Tokyo/en/System.Extended > 5.00000000000000E-0001 > > 5.00000000000000E-0001 > > > > And when executed in 64bit code, returns: That means it has only the 64bit double type, so the +0.5 is lost. > 5.00000000000000E-0001 > > 0.00000000000000E+0000 Mattias _______________________________________________ fpc-devel maillist - http://lists.freepascal.org/cgi-bin/mailman/listinfo/fpc-devel

Reply by email

Thorsten Engler
to Sven Barth: Bart schrieb am Fr., 27. Apr. 2018, 13:42: The code essentially does the following…

I’ve only tested it in Delphi, so you’ll have to convert it to the right syntax for fpc, but either of these should do:

function Frac1(const X: Double): Double;

asm

.align 16

.noframe

movq rdx, xmm0

shl rdx, 32

and edx, $7FF00000

cmp edx, $43300000

jge @@zero

cmp edx, $3FE00000

jbe @@skip

cvttsd2si rax, xmm0

cvtsi2sd xmm4, rax

subsd xmm0, xmm4

jmp @@skip

@@zero:

xorpd xmm0, xmm0

@@skip:

end;
function Frac2(const X: Double): Double;

asm

.align 16

.noframe

movq rdx, xmm0

shl rdx, 48

and dx, $7FF0

cmp dx, $4330

jge @@zero

cmp dx, $3FE0

jbe @@skip

cvttsd2si rax, xmm0

cvtsi2sd xmm4, rax

subsd xmm0, xmm4

jmp @@skip

@@zero:

xorpd xmm0, xmm0

@@skip:

end;
Earlier message 59 lines
From: fpc-devel On Behalf Of Sven Barth via fpc-devel Sent: Friday, 27 April 2018 23:47 To: FPC developers' list Cc: Sven Barth Subject: *** GMX Spamverdacht *** Re: [fpc-devel] Broken frac function in FPC3.1.1 / Windows x86_64 Bart < > schrieb am Fr., 27. Apr. 2018, 13:42: On Wed, Apr 25, 2018 at 11:04 AM, < > wrote: > If you compile and run this 64-bit program on Win 64 you get a crash And AFAICS your analysis of the cause (see bugtracker) is correct as well. function fpc_frac_real(d: ValReal) : ValReal;compilerproc; assembler; nostackframe; asm cvttsd2si %xmm0,%rax { Windows defines %xmm4 and %xmm5 as first non-parameter volatile registers; on SYSV systems all are considered as such, so use %xmm4 } cvtsi2sd %rax,%xmm4 subsd %xmm4,%xmm0 end; CVTTSD2SI — Convert with Truncation Scalar Double-Precision Floating-Point Value to Signed Integer This should not be used to get a ValReal result. The code essentially does the following (instruction by instruction): === code begin === tmpi := int64(d - trunc(d)); tmpd := double(tmpi); Result := d - tmpd; === code end === Though why it fails with the given value is a different topic... Regards, Sven

Reply by email

Thorsten Engler

Oops, small mistake caused by last minute change (I replaced rol with shl): it needs to be shr (or ror or rol, they all perform about the same on my cpu).

And in case anyone wonders, the first cmp and branch returns 0 for numbers that would cause an integer overflow, and the 2nd cmp and branch skips the whole thing if the input is between -1 and 1 (so it already is just a fraction).

Also, at least on my CPU (AMD Phenom II X6, so not exactly the newest) the effect of code alignment on performance is huge. (It affects the branch predictor I think.) I’m not sure what the best alignment is for all CPUs.

.align 16 forces alignment of the function entry point to a multiple of 16. If I add between 0 and 3 nops at the start of the function, the timings for calling it 10 million times are:

In range (1e15+0.5):

0 nop 1266149

1 nop 4260343

2 nop 1369745

3 nop 4469482

Out of range (1e16+0.5):

0 nop 881536

1 nop 896240

2 nop 890805

3 nop 871582

Only fraction (0.5):

0 nop 894850

1 nop 2219469

2 nop 955618

3 nop 2303233

Leaving out the check if it’s already a fraction decreases the time for in-range numbers and increases it for once that are already a faction:

In range (1e15+0.5):

do skip 1306063

no skip 1121395

Out of range (1e16+0.5):

do skip 887081

no skip 888925

Only fraction (0.5):

do skip 903330

no skip 1124026

function FracDoSkip(const X: Double): Double;

asm

.align 16

.noframe

movq rdx, xmm0

rol rdx, 32

and edx, $7FF00000

cmp edx, $43300000

jge @@zero

cmp edx, $3FE00000

jbe @@skip

cvttsd2si rax, xmm0

cvtsi2sd xmm4, rax

subsd xmm0, xmm4

jmp @@skip

@@zero:

xorpd xmm0, xmm0

@@skip:

end;
function FracNoSkip(const X: Double): Double;

asm

.align 16

.noframe

movq rdx, xmm0

rol rdx, 32

and edx, $7FF00000

cmp edx, $43300000

jge @@zero

//  cmp       edx, $3FE00000
//  jbe      @@skip
  cvttsd2si rax, xmm0
  cvtsi2sd  xmm4, rax
  subsd     xmm0, xmm4
  jmp      @@skip

@@zero:

xorpd xmm0, xmm0

@@skip:

end;

Cheers,

Thorsten

Earlier message 152 lines
From: fpc-devel On Behalf Of Thorsten Engler Sent: Saturday, 28 April 2018 15:37 To: 'FPC developers' list' Subject: Re: [fpc-devel] *** GMX Spamverdacht *** Re: Broken frac function in FPC3.1.1 / Windows x86_64 I’ve only tested it in Delphi, so you’ll have to convert it to the right syntax for fpc, but either of these should do: function Frac1(const X: Double): Double; asm .align 16 .noframe movq rdx, xmm0 shl rdx, 32 and edx, $7FF00000 cmp edx, $43300000 jge @@zero cmp edx, $3FE00000 jbe @@skip cvttsd2si rax, xmm0 cvtsi2sd xmm4, rax subsd xmm0, xmm4 jmp @@skip @@zero: xorpd xmm0, xmm0 @@skip: end; function Frac2(const X: Double): Double; asm .align 16 .noframe movq rdx, xmm0 shl rdx, 48 and dx, $7FF0 cmp dx, $4330 jge @@zero cmp dx, $3FE0 jbe @@skip cvttsd2si rax, xmm0 cvtsi2sd xmm4, rax subsd xmm0, xmm4 jmp @@skip @@zero: xorpd xmm0, xmm0 @@skip: end; From: fpc-devel < > On Behalf Of Sven Barth via fpc-devel Sent: Friday, 27 April 2018 23:47 To: FPC developers' list < > Cc: Sven Barth < > Subject: *** GMX Spamverdacht *** Re: [fpc-devel] Broken frac function in FPC3.1.1 / Windows x86_64 Bart < > schrieb am Fr., 27. Apr. 2018, 13:42: On Wed, Apr 25, 2018 at 11:04 AM, < > wrote: > If you compile and run this 64-bit program on Win 64 you get a crash And AFAICS your analysis of the cause (see bugtracker) is correct as well. function fpc_frac_real(d: ValReal) : ValReal;compilerproc; assembler; nostackframe; asm cvttsd2si %xmm0,%rax { Windows defines %xmm4 and %xmm5 as first non-parameter volatile registers; on SYSV systems all are considered as such, so use %xmm4 } cvtsi2sd %rax,%xmm4 subsd %xmm4,%xmm0 end; CVTTSD2SI — Convert with Truncation Scalar Double-Precision Floating-Point Value to Signed Integer This should not be used to get a ValReal result. The code essentially does the following (instruction by instruction): === code begin === tmpi := int64(d - trunc(d)); tmpd := double(tmpi); Result := d - tmpd; === code end === Though why it fails with the given value is a different topic... Regards, Sven

Reply by email

Sven Barth
to Thorsten Engler
Thorsten Engler wrote 19 lines
Oops, small mistake caused by last minute change (I replaced rol with shl): it needs to be shr (or ror or rol, they all perform about the same on my cpu). And in case anyone wonders, the first cmp and branch returns 0 for numbers that would cause an integer overflow, and the 2^nd cmp and branch skips  the whole thing if the input is between -1 and 1 (so it already is just a fraction). Also, at least on my CPU (AMD Phenom II X6, so not exactly the newest) the effect of code alignment on performance is huge. (It affects the branch predictor I think.) I’m not sure what the best alignment is for all CPUs. .align 16 forces alignment of the function entry point to a multiple of 16. If I add between 0 and 3 nops at the start of the function, the timings for calling it 10 million times are:

The question is whether we should really return 0 for numbers that would cause an integer overflow as from the user's perspective of this function integers aren't involved at all (after all one passes in a floating point value and receives one).
Maybe it would be better if we'd provide an optimized Int() function and Frac() then simply stays "d - Int(d)".

Regards,
Sven

Reply by email

Mattias Gaertner
to Sven Barth
Show quoted text 7 lines
[...] The question is whether we should really return 0 for numbers that would cause an integer overflow as from the user's perspective of this function integers aren't involved at all (after all one passes in a floating point value and receives one). Maybe it would be better if we'd provide an optimized Int() function and Frac() then simply stays "d - Int(d)".

Since a float only stores the highest bits of a number, the fraction of a high number is not stored aka treated as 0.

Mattias

Reply by email

Thorsten Engler
to Mattias Gaertner

Basically that.

For doubles that don't fit into an Int64, any fraction is beyond the significant number of digits that the double can store anyway, so 0 is the valid and correct result for Frac in that case.

Show quoted text 31 lines
-----Original Message----- From: fpc-devel On Behalf Of Mattias Gaertner Sent: Saturday, 28 April 2018 19:18 To: Subject: Re: [fpc-devel] Broken frac function in FPC3.1.1 / Windows x86_64 On Sat, 28 Apr 2018 10:55:10 +0200 Sven Barth via fpc-devel wrote: >[...] > The question is whether we should really return 0 for numbers that >would cause an integer overflow as from the user's perspective of this >function integers aren't involved at all (after all one passes in a >floating point value and receives one). > Maybe it would be better if we'd provide an optimized Int() function >and > Frac() then simply stays "d - Int(d)". Since a float only stores the highest bits of a number, the fraction of a high number is not stored aka treated as 0. Mattias _______________________________________________ fpc-devel maillist - http://lists.freepascal.org/cgi-bin/mailman/listinfo/fpc-devel

Reply by email

Sven Barth
to Thorsten Engler
Thorsten Engler wrote 7 lines
Basically that. For doubles that don't fit into an Int64, any fraction is beyond the significant number of digits that the double can store anyway, so 0 is the valid and correct result for Frac in that case. > > Since a float only stores the highest bits of a number, the > fraction of a high number is not stored aka treated as 0.

Right. It's clearly too early in the morning to think about such implications :P
In that case: Thorsten, do you have tests that could go with the new Frac() implementation? Something simply that checks the results and does a Halt() if something is wrong would be sufficent. This way I can be sure that my conversion to FPC's syntax will be correct.

Oh, and before I forget it: would you *please* stop top posting? (Except you have a valid reason for it like being blind and using a screen reader)

Regards,
Sven

Reply by email

Thorsten Engler
to Sven Barth
Show quoted text 7 lines
-----Original Message----- From: fpc-devel On Behalf Of Sven Barth via fpc-devel Thorsten, do you have tests that could go with the new Frac() implementation? Something simply that checks the results and does a Halt() if something is wrong would be sufficent.

I was using the following to verify my implementation while writing it (and if I had remembered to actually run it again instead of making the rol to shl change only when posting the source here, I would have noticed my little oops sooner :/ ):

    Assert(Frac1(1) = 0);
    Assert(Frac1(1e15+0.5) = 0.5);
    Assert(Frac1(1e16+0.5) = 0);
    Assert(Frac1(0.5) = 0.5);
    Assert(Frac1(-1) = 0);
    Assert(Frac1(-1e15-0.5) = -0.5);
    Assert(Frac1(-1e16-0.5) = 0);
    Assert(Frac1(-0.5) = -0.5);

Not sure what exact format you want, but I'm sure you can adapt that as needed.

This way I can be sure that my conversion to FPC's syntax will be correct.

I also used this to verify that my implementation produces the same output as Delphi:

    Assert(Frac1(1) = Frac(1));
    Assert(Frac1(1e15+0.5) = Frac(1e15+0.5) );
    Assert(Frac1(1e16+0.5) = Frac(1e16+0.5) );
    Assert(Frac1(0.5) = Frac(0.5));
    Assert(Frac1(-1) = Frac(-1));
    Assert(Frac1(-1e15-0.5) = Frac(-1e15-0.5));
    Assert(Frac1(-1e16-0.5) = Frac(-1e16-0.5));
    Assert(Frac1(-0.5) = Frac(-0.5));
Oh, and before I forget it: would you *please* stop top posting?

Oh, sure. No worries. Standards about that vary from community to community...

The effects of code alignment beyond a granularity of 16 on such short code is interesting:

Code address:
Frac1: 0000000000536430 (48)
Frac2: 0000000000536480 (0)
Frac3: 00000000005364D0 (80)
Frac4: 0000000000536520 (32)
Frac5: 0000000000536570 (112)
Frac6: 00000000005365C0 (64)
Frac7: 0000000000536610 (16)
Frac8: 0000000000536660 (96)

1st run:
In range (1e15+0.5):
Frac 1872424
Frac1 1292064
Frac2 1302669
Frac3 1315819
Frac4 1301370
Frac5 1305682
Frac6 1288991
Frac7 1312695
Frac8 1294104

Out of range (1e16+0.5):
Frac 1605912
Frac1 1081226
Frac2 892385
Frac3 1078618
Frac4 888650
Frac5 1077251
Frac6 1096165
Frac7 1072769
Frac8 882501

Only fraction (0.5):
Frac 1638980
Frac1 990573
Frac2 893984
Frac3 989789
Frac4 884970
Frac5 975254
Frac6 986453
Frac7 983450
Frac8 877745

The number in () after the address is address mod 128.

It's interesting to see that for this particular code (and CPU) the "good" ones are 0, 32, and 96, but NOT 64. I'm sure the results are highly dependent on the CPU.

Also, you can see that no matter what the alignment, it beats Delphi's 64bit Frac implementation by a good margin.

Reply by email

Unknown author
to Thorsten Engler
Thorsten Engler wrote 12 lines
The effects of code alignment beyond a granularity of 16 on such short code is interesting: Code address: Frac1: 0000000000536430 (48) Frac2: 0000000000536480 (0) Frac3: 00000000005364D0 (80) Frac4: 0000000000536520 (32) Frac5: 0000000000536570 (112) Frac6: 00000000005365C0 (64) Frac7: 0000000000536610 (16) Frac8: 0000000000536660 (96)

[...]

The number in () after the address is address mod 128. It's interesting to see that for this particular code (and CPU) the "good" ones are 0, 32, and 96, but NOT 64. I'm sure the results are highly dependent on the CPU.

why not 64? the pattern looks to be bad,good,bad,good,bad,good,bad,good but i'm very likely missing something...

also, not only highly dependent on the CPU but also on what other processes may be running and consuming some CPU time... i'm not even sure that booting linux to "single" mode would get you a system completely dedicated to one task like in the old DOS world...

NOTE: No off-list assistance is given without prior approval. *Please keep mailing list traffic on the list unless* *a signed and pre-paid contract is in effect with us.*

Reply by email

Thorsten Engler
to Unknown author
Show quoted text 15 lines
-----Original Message----- From: fpc-devel On Behalf Of > Code address: > Frac1: 0000000000536430 (48) > Frac2: 0000000000536480 (0) > Frac3: 00000000005364D0 (80) > Frac4: 0000000000536520 (32) > Frac5: 0000000000536570 (112) > Frac6: 00000000005365C0 (64) > Frac7: 0000000000536610 (16) > Frac8: 0000000000536660 (96) why not 64? the pattern looks to be bad,good,bad,good,bad,good,bad,good but i'm very likely missing something...

Take a look at the pattern again, I've marked it with + and - this time:

Out of range (1e16+0.5):
-Frac1 1081226
+Frac2 892385
-Frac3 1078618
+Frac4 888650
-Frac5 1077251
-Frac6 1096165
-Frac7 1072769
+Frac8 882501

Only fraction (0.5):
-Frac1 990573
+Frac2 893984
-Frac3 989789
+Frac4 884970
-Frac5 975254
-Frac6 986453
-Frac7 983450
+Frac8 877745

Why not 64? No idea...

Some more funny points, if I set CODEALIGN to 1, so the 8 functions get aligned all over the board, I get the following result:

Code address:
Frac1: 0000000000536422 (34)
Frac2: 0000000000536458 (88)
Frac3: 000000000053648E (14)
Frac4: 00000000005364C4 (68)
Frac5: 00000000005364FA (122)
Frac6: 0000000000536530 (48)
Frac7: 0000000000536566 (102)
Frac8: 000000000053659C (28)

1st run:
In range (1e15+0.5):
Frac1 1322656
Frac2 1326933
Frac3 1321368
Frac4 1311924
Frac5 1302156
Frac6 1311717
Frac7 1314821
Frac8 1304493

Out of range (1e16+0.5):
Frac1 796050
Frac2 914107
Frac3 1086412
Frac4 792253
Frac5 999116
Frac6 1091263
Frac7 792960
Frac8 996593

Only fraction (0.5):
Frac1 889606
Frac2 887059
Frac3 978125
Frac4 888643
Frac5 989546
Frac6 990639
Frac7 891451
Frac8 993132

The in range numbers are pretty much all the same, the only fraction numbers are in the same ballpark for the "good" and "bad" cases, but for the out of range cases, for case 1, 4, and 7 (offsets of 34, 68, 102) the result is actually about 8-10% better than for offsets 0, 32, and 96.

The real shocker comes if I run it again after adding a single additional nop in a spacer function that I have above Frac1:

procedure XXX1;
asm
.noframe
nop
nop // added this
end;

function Frac1(const X: Double): Double;
...

Code address:
Frac1: 0000000000536423 (35)
Frac2: 0000000000536459 (89)
Frac3: 000000000053648F (15)
Frac4: 00000000005364C5 (69)
Frac5: 00000000005364FB (123)
Frac6: 0000000000536531 (49)
Frac7: 0000000000536567 (103)
Frac8: 000000000053659D (29)

1st run:
In range (1e15+0.5):
Frac1 4541094
Frac2 4451186
Frac3 4446460
Frac4 4450625
Frac5 4323903
Frac6 4389160
Frac7 4422070
Frac8 1291440

Out of range (1e16+0.5):
Frac1 892037
Frac2 882626
Frac3 1082966
Frac4 785305
Frac5 887439
Frac6 991194
Frac7 789911
Frac8 987714

Only fraction (0.5):
Frac1 2245965
Frac2 2172882
Frac3 2260256
Frac4 2154653
Frac5 2158976
Frac6 2257709
Frac7 2147712
Frac8 979118

Total blowout. Time for In Range have ballooned by factor 3-4, for only fraction by over 2.

Clearly, even if you don't care about the differences between aligning to offset 0 vs 16, these show that proper alignment to at least a multiple of 2 is absolutely essential.

I find it interesting that in this case offset 29 (Frac8) performs pretty well for the in range and only fraction cases, but badly for out of range. And offset 103 (Frac7) performs very well for the out of range case, but just as badly as most of the others for in range or fraction only.

I've attached the source (I'm using Delphi 10.2.3, 64bit to compile it) in case anyone wants to try it out on different cpus and with different alignments (change the {$CODEALIGN 1} and add nops to the XXX1 .. XXX8 procedures to finetune alignment).

Show quoted text 7 lines
also, not only highly dependent on the CPU but also on what other processes may be running and consuming some CPU time... i'm not even sure that booting linux to "single" mode would get you a system completely dedicated to one task like in the old DOS world...

There are two approaches you can use under windows to get pretty reliable timings:

a)

  SetThreadPriority(GetCurrentThread, THREAD_PRIORITY_TIME_CRITICAL);
...
  SetThreadPriority(GetCurrentThread, THREAD_PRIORITY_NORMAL);

As long as the total number of active (not in a wait state) realtime threads is smaller than the number of cores, this does a very good job of keeping your process running continuously. You could use SetProcessAffinityMask on top of that to make sure that it keeps running on one specific core. This very greatly reduces the impact of any other background activity on the timing.

This is what I used in this case because it was very easy to do.

The timing results are very consistent between multiple runs (within 1-3% I would say, definitely much below the difference between the different alignments).

b)

instead of the normal high precision timer, use:

function GetThreadTimes(hThread: THandle;
  var lpCreationTime, lpExitTime, lpKernelTime, lpUserTime: TFileTime): BOOL; stdcall;

lpKernelTime and lpUserTime only count the time while the thread is executing, not any time that it spends in a wait state or waiting to be scheduled to run.

This can give you pretty exact timing values for how long your code actually executed, independent of how often the thread was suspended to allow some other thread to run.

You do still have some variation as every thread context switch ends up resetting the TLBs which can significantly impact performance of memory access. So a combination of a) and b) gives the best possible results while running on a pre-emptive multi-tasking OS.

Reply by email

Florian Klämpfl
to Thorsten Engler
Thorsten Engler wrote 6 lines
procedure XXX1; asm .noframe nop nop // added this end;

I did not look at the code in detail but I suspect this is caused by the two branches:

cmp edx, $43300000
jge @@zero
cmp edx, $3FE00000
jbe @@skip

If two branches are in the same 16 byte block, they share a branch prediction entry on almost all x86 processors and branch prediction suffers very much in this case.

So something like

cmp edx, $43300000
jge @@zero
cmp edx, $3FE00000
.align 16
jbe @@skip

might be much better.

Reply by email

Unknown author
to Thorsten Engler: Take a look at the pattern again, I've marked it with + and - this time: Out of range…
Thorsten Engler wrote
I've attached the source (I'm using Delphi 10.2.3, 64bit to compile it) in case anyone wants to try it out on different cpus and with different alignments (change the {$CODEALIGN 1} and add nops to the XXX1 .. XXX8 procedures to finetune alignment).

FWIW: i tried, with my old lazarus 1.6.1 and fpc 3.0.0, to convert the project to a lazarus project but it aborted with no real indication as to the problem... i didn't try with the (also old) trunk install i have... both are pretty old and i've forgotten how i set it all up with fpcup...

so i tried just compiling the lpr that did result from the conversion attempt... it looks to be valid source but the compile also failed... mainly for not being able to find System.Diagnostics...

maybe i'll muddle about with it later... i was just curious to see what this AMD FX Black FX8350 4Ghz 8-core CPU running linux would do...

NOTE: No off-list assistance is given without prior approval. *Please keep mailing list traffic on the list unless* *a signed and pre-paid contract is in effect with us.*

Reply by email

Thorsten Engler
to Florian Klämpfl: I did not look at the code in detail but I suspect this is caused by the two branches:…
Show quoted text 12 lines
-----Original Message----- From: fpc-devel On Behalf Of Florian Klämpfl So something like cmp edx, $43300000 jge @@zero cmp edx, $3FE00000 .align 16 jbe @@skip might be much better.

That ended up making things worse in some cases.

Here is a branchless version:

function Frac1(const X: Double): Double;
asm
.noframe
movq rdx, xmm0
mov rax, rdx
xor rcx, rcx
shr rdx, 32
and edx, $7FF00000
cmp edx, $43300000
cmovge rax, rcx
movq xmm0, rax
cvttsd2si rax, xmm0
cvtsi2sd xmm4, rax
subsd xmm0, xmm4
end;

It performs slightly slower in the "in range" case, noticeable worse in the other 2 cases (as it's exactly the same for all 3).

I would guess that the "in range" case is the most common (you aren't going to call Frac if you know ahead of time that it's always 0 as the number is too big, or if you know that it already is a value between -1 and 1), so the higher cost for the out of range and only fraction cases is probably less important than it might look.

It IS largely independent of code alignment or predictable patterns in the incoming value:

Code address:
Frac1: 0000000000536430 (48)
Frac2: 0000000000536480 (0)
Frac3: 00000000005364D0 (80)
Frac4: 0000000000536520 (32)
Frac5: 0000000000536570 (112)
Frac6: 00000000005365C0 (64)
Frac7: 0000000000536610 (16)
Frac8: 0000000000536660 (96)

1st run:
In range (1e15+0.5):
Frac1 1431794
Frac2 1429232
Frac3 1463357
Frac4 1475042
Frac5 1446016
Frac6 1472979
Frac7 1443244
Frac8 1467528

Out of range (1e16+0.5):
Frac1 1476556
Frac2 1458534
Frac3 1444431
Frac4 1427287
Frac5 1427326
Frac6 1427472
Frac7 1428914
Frac8 1419654

Only fraction (0.5):
Frac1 1470644
Frac2 1475227
Frac3 1447379
Frac4 1529162
Frac5 1509275
Frac6 1485185
Frac7 1500826
Frac8 1524294

Code address:
Frac1: 0000000000536423 (35)
Frac2: 0000000000536458 (88)
Frac3: 000000000053648D (13)
Frac4: 00000000005364C2 (66)
Frac5: 00000000005364F7 (119)
Frac6: 000000000053652C (44)
Frac7: 0000000000536561 (97)
Frac8: 0000000000536596 (22)

1st run:
In range (1e15+0.5):
Frac1 1349334
Frac2 1429198
Frac3 1447011
Frac4 1436476
Frac5 1477058
Frac6 1496887
Frac7 1431293
Frac8 1435460

Out of range (1e16+0.5):
Frac1 1349939
Frac2 1412543
Frac3 1462295
Frac4 1442081
Frac5 1512579
Frac6 1453593
Frac7 1457510
Frac8 1436533

Only fraction (0.5):
Frac1 1371353
Frac2 1443000
Frac3 1437583
Frac4 1415591
Frac5 1474870
Frac6 1437224
Frac7 1452196
Frac8 1453833

Also, it still outperforms Delphi's Frac in all cases.

Reply by email

J. Gareth Moreton

If I may add a possible solution, would code akin to the following fix the problems? (Pray that the mail list formatting doesn't mess it all up!)

----

 {$WARN 7121 off : Check size of memory operand "$1: memory-operand-size is
$2 bits, but expected [$3 bits]"}
 {$ASMMODE Intel}
 const
   MAX_DOUBLE_WITH_FRAC: Double = 4503599627370496;
   SIGN_MASK: QWord  = $7FFFFFFFFFFFFFFF;

function SafeFrac(X: ValReal): ValReal; assembler; nostackframe; asm
  MOVSD     XMM5, [RIP+SIGN_MASK]
  MOVSD     XMM4, XMM0
  ANDPD     XMM4, XMM5 { Remove sign bit }
  COMISD    XMM4, [RIP+MAX_DOUBLE_WITH_FRAC] { If less than, then ZF = 0, PF = 0 and CF = 1}
  JNC       @IsZero    { If carry flag is clear, the absolute value is too
                         large to contain a fractional component (and
                         potentially too large to fit into RAX) }
  CVTTSD2SI RAX,  XMM0
  CVTSI2SD  XMM4, RAX
  SUBPD     XMM0, XMM4
  RET
@IsZero:
  XORPD     XMM0, XMM0
end;

----

The first MOVSD could be removed and simply be merged with ANDPD to save on a register, but I wrote it this way to take advantage of instruction prefetching and offer a slight speed boost.

There is one minor functional change though... when testing this function, I discovered that if you try to pass infinity or minus infinity, the result is zero, whereas currently with Frac, it triggers SIGFPE.

Gareth aka. Kit

Reply by email

J. Gareth Moreton

I've done some speed and accuracy comparisons between our respective Frac functions.  Initially, my "SafeFrac" was marginally faster than "FracDoSkip", but I managed to optimise Thorsten's routine a little bit into the following:

function FracSkip2(const X: ValReal): ValReal; assembler; nostackframe; asm
  align 16
  movq      rax,  xmm0
  shr       rax,  48
  and       ax,   $7FF0
  cmp       ax,   $4330
  jge       @@zero
  cmp       ax,   $3FE0
  jbe       @@skip
  cvttsd2si rax,  xmm0
  cvtsi2sd  xmm4, rax
  subsd     xmm0, xmm4
  ret
@@zero:
  xorpd     xmm0, xmm0
@@skip:
end;

My test compared Frac, FracDoSkip, SafeFrac and what I call FracSkip2 above, which reworks the comparisons to use only 16 bits, and replaces "jmp @@skip" with "ret". The results are as follows (as you can see... all of them are a great improvement over Frac).  Frac raises SIGFPE if plus or minus infinity is passed in, but our functions return zero.  This may or may not be a desirable change.

Code sizes (alignment will round it up to the nearest 16 bytes): Frac = 49 bytes, FracDoSkip = 52 bytes, SafeFrac = 46 bytes, FracSkip2 = 45 bytes.

Long story short, with a few tweaks, Thorsten's routine is the fastest and also the smallest.

****

My test set was:

DataSet: array[0..14] of Double = (1.5, 0, 2251799813685248, 4503599627370496, 1E300, 0.125, 3.6415926535897932384626433832795, -1.5, -2251799813685248, -4503599627370496, -1E300, -0.125, -3.6415926535897932384626433832795, Infinity, NegInfinity);

For each value, it is tested as is, then DataSet[X] + 0.5, then DataSet[X]
- 0.5 (best way to determine how it handles precision without it being optimised out by the compiler).

****

Frac( 1.5000000000000000E+000) =  5.0000000000000000E-001 - Pass - Time = 124.483 ns
FracDoSkip( 1.5000000000000000E+000) =  5.0000000000000000E-001 - Pass - Time = 47.525 ns
SafeFrac( 1.5000000000000000E+000) =  5.0000000000000000E-001 - Pass - Time = 32.707 ns
FracSkip2( 1.5000000000000000E+000) =  5.0000000000000000E-001 - Pass - Time = 34.904 ns
Frac( 2.0000000000000000E+000) =  0.0000000000000000E+000 - Pass - Time = 126.170 ns
FracDoSkip( 2.0000000000000000E+000) =  0.0000000000000000E+000 - Pass - Time = 51.210 ns
SafeFrac( 2.0000000000000000E+000) =  0.0000000000000000E+000 - Pass - Time = 35.351 ns
FracSkip2( 2.0000000000000000E+000) =  0.0000000000000000E+000 - Pass - Time = 33.911 ns
Frac( 1.0000000000000000E+000) =  0.0000000000000000E+000 - Pass - Time = 125.927 ns
FracDoSkip( 1.0000000000000000E+000) =  0.0000000000000000E+000 - Pass - Time = 49.127 ns
SafeFrac( 1.0000000000000000E+000) =  0.0000000000000000E+000 - Pass - Time = 34.695 ns
FracSkip2( 1.0000000000000000E+000) =  0.0000000000000000E+000 - Pass - Time = 36.139 ns
Frac( 0.0000000000000000E+000) =  0.0000000000000000E+000 - Pass - Time = 119.800 ns
FracDoSkip( 0.0000000000000000E+000) =  0.0000000000000000E+000 - Pass - Time = 40.316 ns
SafeFrac( 0.0000000000000000E+000) =  0.0000000000000000E+000 - Pass - Time = 35.875 ns
FracSkip2( 0.0000000000000000E+000) =  0.0000000000000000E+000 - Pass - Time = 34.046 ns
Frac( 5.0000000000000000E-001) =  5.0000000000000000E-001 - Pass - Time = 118.913 ns
FracDoSkip( 5.0000000000000000E-001) =  5.0000000000000000E-001 - Pass - Time = 40.183 ns
SafeFrac( 5.0000000000000000E-001) =  5.0000000000000000E-001 - Pass - Time = 36.783 ns
FracSkip2( 5.0000000000000000E-001) =  5.0000000000000000E-001 - Pass - Time = 34.976 ns
Frac(-5.0000000000000000E-001) = -5.0000000000000000E-001 - Pass - Time = 127.560 ns
FracDoSkip(-5.0000000000000000E-001) = -5.0000000000000000E-001 - Pass - Time = 41.676 ns
SafeFrac(-5.0000000000000000E-001) = -5.0000000000000000E-001 - Pass - Time = 36.577 ns
FracSkip2(-5.0000000000000000E-001) = -5.0000000000000000E-001 - Pass - Time = 34.714 ns
Frac( 2.2517998136852480E+015) =  0.0000000000000000E+000 - Pass - Time = 126.323 ns
FracDoSkip( 2.2517998136852480E+015) =  0.0000000000000000E+000 - Pass - Time = 49.108 ns
SafeFrac( 2.2517998136852480E+015) =  0.0000000000000000E+000 - Pass - Time = 35.376 ns
FracSkip2( 2.2517998136852480E+015) =  0.0000000000000000E+000 - Pass - Time = 36.373 ns
Frac( 2.2517998136852485E+015) =  5.0000000000000000E-001 - Pass - Time = 131.001 ns
FracDoSkip( 2.2517998136852485E+015) =  5.0000000000000000E-001 - Pass - Time = 54.474 ns
SafeFrac( 2.2517998136852485E+015) =  5.0000000000000000E-001 - Pass - Time = 38.834 ns
FracSkip2( 2.2517998136852485E+015) =  5.0000000000000000E-001 - Pass - Time = 37.139 ns
Frac( 2.2517998136852475E+015) =  5.0000000000000000E-001 - Pass - Time = 131.932 ns
FracDoSkip( 2.2517998136852475E+015) =  5.0000000000000000E-001 - Pass - Time = 52.214 ns
SafeFrac( 2.2517998136852475E+015) =  5.0000000000000000E-001 - Pass - Time = 37.093 ns
FracSkip2( 2.2517998136852475E+015) =  5.0000000000000000E-001 - Pass - Time = 35.674 ns
Frac( 4.5035996273704960E+015) =  0.0000000000000000E+000 - Pass - Time = 82.749 ns
FracDoSkip( 4.5035996273704960E+015) =  0.0000000000000000E+000 - Pass - Time = 38.613 ns
SafeFrac( 4.5035996273704960E+015) =  0.0000000000000000E+000 - Pass - Time = 38.575 ns
FracSkip2( 4.5035996273704960E+015) =  0.0000000000000000E+000 - Pass - Time = 33.970 ns
Frac( 4.5035996273704960E+015) =  0.0000000000000000E+000 - Pass - Time = 86.126 ns
FracDoSkip( 4.5035996273704960E+015) =  0.0000000000000000E+000 - Pass - Time = 38.434 ns
SafeFrac( 4.5035996273704960E+015) =  0.0000000000000000E+000 - Pass - Time = 38.636 ns
FracSkip2( 4.5035996273704960E+015) =  0.0000000000000000E+000 - Pass - Time = 33.747 ns
Frac( 4.5035996273704955E+015) =  5.0000000000000000E-001 - Pass - Time = 131.589 ns
FracDoSkip( 4.5035996273704955E+015) =  5.0000000000000000E-001 - Pass - Time = 53.594 ns
SafeFrac( 4.5035996273704955E+015) =  5.0000000000000000E-001 - Pass - Time = 36.617 ns
FracSkip2( 4.5035996273704955E+015) =  5.0000000000000000E-001 - Pass - Time = 36.509 ns
Frac( 1.0000000000000001E+300) =  0.0000000000000000E+000 - Pass - Time = 82.875 ns
FracDoSkip( 1.0000000000000001E+300) =  0.0000000000000000E+000 - Pass - Time = 39.008 ns
SafeFrac( 1.0000000000000001E+300) =  0.0000000000000000E+000 - Pass - Time = 39.112 ns
FracSkip2( 1.0000000000000001E+300) =  0.0000000000000000E+000 - Pass - Time = 34.195 ns
Frac( 1.0000000000000001E+300) =  0.0000000000000000E+000 - Pass - Time = 85.401 ns
FracDoSkip( 1.0000000000000001E+300) =  0.0000000000000000E+000 - Pass - Time = 38.653 ns
SafeFrac( 1.0000000000000001E+300) =  0.0000000000000000E+000 - Pass - Time = 38.655 ns
FracSkip2( 1.0000000000000001E+300) =  0.0000000000000000E+000 - Pass - Time = 34.408 ns
Frac( 1.0000000000000001E+300) =  0.0000000000000000E+000 - Pass - Time = 84.719 ns
FracDoSkip( 1.0000000000000001E+300) =  0.0000000000000000E+000 - Pass - Time = 39.174 ns
SafeFrac( 1.0000000000000001E+300) =  0.0000000000000000E+000 - Pass - Time = 38.876 ns
FracSkip2( 1.0000000000000001E+300) =  0.0000000000000000E+000 - Pass - Time = 33.570 ns
Frac( 1.2500000000000000E-001) =  1.2500000000000000E-001 - Pass - Time = 123.770 ns
FracDoSkip( 1.2500000000000000E-001) =  1.2500000000000000E-001 - Pass - Time = 41.642 ns
SafeFrac( 1.2500000000000000E-001) =  1.2500000000000000E-001 - Pass - Time = 38.704 ns
FracSkip2( 1.2500000000000000E-001) =  1.2500000000000000E-001 - Pass - Time = 35.399 ns
Frac( 6.2500000000000000E-001) =  6.2500000000000000E-001 - Pass - Time = 128.967 ns
FracDoSkip( 6.2500000000000000E-001) =  6.2500000000000000E-001 - Pass - Time = 42.082 ns
SafeFrac( 6.2500000000000000E-001) =  6.2500000000000000E-001 - Pass - Time = 38.199 ns
FracSkip2( 6.2500000000000000E-001) =  6.2500000000000000E-001 - Pass - Time = 36.072 ns
Frac(-3.7500000000000000E-001) = -3.7500000000000000E-001 - Pass - Time = 128.962 ns
FracDoSkip(-3.7500000000000000E-001) = -3.7500000000000000E-001 - Pass - Time = 40.375 ns
SafeFrac(-3.7500000000000000E-001) = -3.7500000000000000E-001 - Pass - Time = 37.153 ns
FracSkip2(-3.7500000000000000E-001) = -3.7500000000000000E-001 - Pass - Time = 34.515 ns
Frac( 3.6415926535897931E+000) =  6.4159265358979312E-001 - Pass - Time = 129.245 ns
FracDoSkip( 3.6415926535897931E+000) =  6.4159265358979312E-001 - Pass - Time = 53.440 ns
SafeFrac( 3.6415926535897931E+000) =  6.4159265358979312E-001 - Pass - Time = 38.390 ns
FracSkip2( 3.6415926535897931E+000) =  6.4159265358979312E-001 - Pass - Time = 36.963 ns
Frac( 4.1415926535897931E+000) =  1.4159265358979312E-001 - Pass - Time = 132.623 ns
FracDoSkip( 4.1415926535897931E+000) =  1.4159265358979312E-001 - Pass - Time = 52.325 ns
SafeFrac( 4.1415926535897931E+000) =  1.4159265358979312E-001 - Pass - Time = 39.016 ns
FracSkip2( 4.1415926535897931E+000) =  1.4159265358979312E-001 - Pass - Time = 36.818 ns
Frac( 3.1415926535897931E+000) =  1.4159265358979312E-001 - Pass - Time = 128.032 ns
FracDoSkip( 3.1415926535897931E+000) =  1.4159265358979312E-001 - Pass - Time = 49.834 ns
SafeFrac( 3.1415926535897931E+000) =  1.4159265358979312E-001 - Pass - Time = 37.077 ns
FracSkip2( 3.1415926535897931E+000) =  1.4159265358979312E-001 - Pass - Time = 37.099 ns
Frac(-1.5000000000000000E+000) = -5.0000000000000000E-001 - Pass - Time = 132.057 ns
FracDoSkip(-1.5000000000000000E+000) = -5.0000000000000000E-001 - Pass - Time = 53.112 ns
SafeFrac(-1.5000000000000000E+000) = -5.0000000000000000E-001 - Pass - Time = 38.287 ns
FracSkip2(-1.5000000000000000E+000) = -5.0000000000000000E-001 - Pass - Time = 36.849 ns
Frac(-1.0000000000000000E+000) =  0.0000000000000000E+000 - Pass - Time = 130.452 ns
FracDoSkip(-1.0000000000000000E+000) =  0.0000000000000000E+000 - Pass - Time = 51.451 ns
SafeFrac(-1.0000000000000000E+000) =  0.0000000000000000E+000 - Pass - Time = 36.993 ns
FracSkip2(-1.0000000000000000E+000) =  0.0000000000000000E+000 - Pass - Time = 36.110 ns
Frac(-2.0000000000000000E+000) =  0.0000000000000000E+000 - Pass - Time = 131.912 ns
FracDoSkip(-2.0000000000000000E+000) =  0.0000000000000000E+000 - Pass - Time = 52.946 ns
SafeFrac(-2.0000000000000000E+000) =  0.0000000000000000E+000 - Pass - Time = 38.330 ns
FracSkip2(-2.0000000000000000E+000) =  0.0000000000000000E+000 - Pass - Time = 37.156 ns
Frac(-2.2517998136852480E+015) =  0.0000000000000000E+000 - Pass - Time = 131.354 ns
FracDoSkip(-2.2517998136852480E+015) =  0.0000000000000000E+000 - Pass - Time = 53.712 ns
SafeFrac(-2.2517998136852480E+015) =  0.0000000000000000E+000 - Pass - Time = 36.978 ns
FracSkip2(-2.2517998136852480E+015) =  0.0000000000000000E+000 - Pass - Time = 36.262 ns
Frac(-2.2517998136852475E+015) = -5.0000000000000000E-001 - Pass - Time = 127.641 ns
FracDoSkip(-2.2517998136852475E+015) = -5.0000000000000000E-001 - Pass - Time = 52.853 ns
SafeFrac(-2.2517998136852475E+015) = -5.0000000000000000E-001 - Pass - Time = 38.318 ns
FracSkip2(-2.2517998136852475E+015) = -5.0000000000000000E-001 - Pass - Time = 37.286 ns
Frac(-2.2517998136852485E+015) = -5.0000000000000000E-001 - Pass - Time = 130.918 ns
FracDoSkip(-2.2517998136852485E+015) = -5.0000000000000000E-001 - Pass - Time = 52.916 ns
SafeFrac(-2.2517998136852485E+015) = -5.0000000000000000E-001 - Pass - Time = 37.928 ns
FracSkip2(-2.2517998136852485E+015) = -5.0000000000000000E-001 - Pass - Time = 36.701 ns
Frac(-4.5035996273704960E+015) =  0.0000000000000000E+000 - Pass - Time = 82.714 ns
FracDoSkip(-4.5035996273704960E+015) =  0.0000000000000000E+000 - Pass - Time = 37.410 ns
SafeFrac(-4.5035996273704960E+015) =  0.0000000000000000E+000 - Pass - Time = 37.091 ns
FracSkip2(-4.5035996273704960E+015) =  0.0000000000000000E+000 - Pass - Time = 33.130 ns
Frac(-4.5035996273704955E+015) = -5.0000000000000000E-001 - Pass - Time = 131.699 ns
FracDoSkip(-4.5035996273704955E+015) = -5.0000000000000000E-001 - Pass - Time = 52.932 ns
SafeFrac(-4.5035996273704955E+015) = -5.0000000000000000E-001 - Pass - Time = 38.499 ns
FracSkip2(-4.5035996273704955E+015) = -5.0000000000000000E-001 - Pass - Time = 37.341 ns
Frac(-4.5035996273704960E+015) =  0.0000000000000000E+000 - Pass - Time = 85.069 ns
FracDoSkip(-4.5035996273704960E+015) =  0.0000000000000000E+000 - Pass - Time = 38.384 ns
SafeFrac(-4.5035996273704960E+015) =  0.0000000000000000E+000 - Pass - Time = 39.041 ns
FracSkip2(-4.5035996273704960E+015) =  0.0000000000000000E+000 - Pass - Time = 34.266 ns
Frac(-1.0000000000000001E+300) =  0.0000000000000000E+000 - Pass - Time = 81.913 ns
FracDoSkip(-1.0000000000000001E+300) =  0.0000000000000000E+000 - Pass - Time = 37.216 ns
SafeFrac(-1.0000000000000001E+300) =  0.0000000000000000E+000 - Pass - Time = 37.385 ns
FracSkip2(-1.0000000000000001E+300) =  0.0000000000000000E+000 - Pass - Time = 34.328 ns
Frac(-1.0000000000000001E+300) =  0.0000000000000000E+000 - Pass - Time = 85.317 ns
FracDoSkip(-1.0000000000000001E+300) =  0.0000000000000000E+000 - Pass - Time = 38.639 ns
SafeFrac(-1.0000000000000001E+300) =  0.0000000000000000E+000 - Pass - Time = 38.644 ns
FracSkip2(-1.0000000000000001E+300) =  0.0000000000000000E+000 - Pass - Time = 34.293 ns
Frac(-1.0000000000000001E+300) =  0.0000000000000000E+000 - Pass - Time = 85.878 ns
FracDoSkip(-1.0000000000000001E+300) =  0.0000000000000000E+000 - Pass - Time = 38.932 ns
SafeFrac(-1.0000000000000001E+300) =  0.0000000000000000E+000 - Pass - Time = 38.651 ns
FracSkip2(-1.0000000000000001E+300) =  0.0000000000000000E+000 - Pass - Time = 34.316 ns
Frac(-1.2500000000000000E-001) = -1.2500000000000000E-001 - Pass - Time = 128.603 ns
FracDoSkip(-1.2500000000000000E-001) = -1.2500000000000000E-001 - Pass - Time = 41.592 ns
SafeFrac(-1.2500000000000000E-001) = -1.2500000000000000E-001 - Pass - Time = 37.280 ns
FracSkip2(-1.2500000000000000E-001) = -1.2500000000000000E-001 - Pass - Time = 34.995 ns
Frac( 3.7500000000000000E-001) =  3.7500000000000000E-001 - Pass - Time = 124.473 ns
FracDoSkip( 3.7500000000000000E-001) =  3.7500000000000000E-001 - Pass - Time = 42.099 ns
SafeFrac( 3.7500000000000000E-001) =  3.7500000000000000E-001 - Pass - Time = 38.716 ns
FracSkip2( 3.7500000000000000E-001) =  3.7500000000000000E-001 - Pass - Time = 36.194 ns
Frac(-6.2500000000000000E-001) = -6.2500000000000000E-001 - Pass - Time = 129.138 ns
FracDoSkip(-6.2500000000000000E-001) = -6.2500000000000000E-001 - Pass - Time = 42.219 ns
SafeFrac(-6.2500000000000000E-001) = -6.2500000000000000E-001 - Pass - Time = 39.724 ns
FracSkip2(-6.2500000000000000E-001) = -6.2500000000000000E-001 - Pass - Time = 34.307 ns
Frac(-3.6415926535897931E+000) = -6.4159265358979312E-001 - Pass - Time = 129.833 ns
FracDoSkip(-3.6415926535897931E+000) = -6.4159265358979312E-001 - Pass - Time = 51.274 ns
SafeFrac(-3.6415926535897931E+000) = -6.4159265358979312E-001 - Pass - Time = 38.494 ns
FracSkip2(-3.6415926535897931E+000) = -6.4159265358979312E-001 - Pass - Time = 37.459 ns
Frac(-3.1415926535897931E+000) = -1.4159265358979312E-001 - Pass - Time = 132.230 ns
FracDoSkip(-3.1415926535897931E+000) = -1.4159265358979312E-001 - Pass - Time = 53.066 ns
SafeFrac(-3.1415926535897931E+000) = -1.4159265358979312E-001 - Pass - Time = 38.658 ns
FracSkip2(-3.1415926535897931E+000) = -1.4159265358979312E-001 - Pass - Time = 36.351 ns
Frac(-4.1415926535897931E+000) = -1.4159265358979312E-001 - Pass - Time = 126.783 ns
FracDoSkip(-4.1415926535897931E+000) = -1.4159265358979312E-001 - Pass - Time = 51.889 ns
SafeFrac(-4.1415926535897931E+000) = -1.4159265358979312E-001 - Pass - Time = 38.785 ns
FracSkip2(-4.1415926535897931E+000) = -1.4159265358979312E-001 - Pass - Time = 36.711 ns
Frac(+Inf) = EXCEPTION - EInvalidOp raised with message "Invalid floating point operation"
FracDoSkip(+Inf) =  0.0000000000000000E+000 - Pass - Time = 39.849 ns SafeFrac(+Inf) =  0.0000000000000000E+000 - Pass - Time = 38.889 ns FracSkip2(+Inf) =  0.0000000000000000E+000 - Pass - Time = 34.289 ns Frac(+Inf) = EXCEPTION - EInvalidOp raised with message "Invalid floating point operation"
FracDoSkip(+Inf) =  0.0000000000000000E+000 - Pass - Time = 40.781 ns SafeFrac(+Inf) =  0.0000000000000000E+000 - Pass - Time = 37.504 ns FracSkip2(+Inf) =  0.0000000000000000E+000 - Pass - Time = 33.043 ns Frac(+Inf) = EXCEPTION - EInvalidOp raised with message "Invalid floating point operation"
FracDoSkip(+Inf) =  0.0000000000000000E+000 - Pass - Time = 40.993 ns SafeFrac(+Inf) =  0.0000000000000000E+000 - Pass - Time = 39.575 ns FracSkip2(+Inf) =  0.0000000000000000E+000 - Pass - Time = 33.041 ns Frac(-Inf) = EXCEPTION - EInvalidOp raised with message "Invalid floating point operation"
FracDoSkip(-Inf) =  0.0000000000000000E+000 - Pass - Time = 40.414 ns SafeFrac(-Inf) =  0.0000000000000000E+000 - Pass - Time = 37.835 ns FracSkip2(-Inf) =  0.0000000000000000E+000 - Pass - Time = 33.294 ns Frac(-Inf) = EXCEPTION - EInvalidOp raised with message "Invalid floating point operation"
FracDoSkip(-Inf) =  0.0000000000000000E+000 - Pass - Time = 39.871 ns SafeFrac(-Inf) =  0.0000000000000000E+000 - Pass - Time = 37.885 ns FracSkip2(-Inf) =  0.0000000000000000E+000 - Pass - Time = 34.041 ns Frac(-Inf) = EXCEPTION - EInvalidOp raised with message "Invalid floating point operation"
FracDoSkip(-Inf) =  0.0000000000000000E+000 - Pass - Time = 40.437 ns SafeFrac(-Inf) =  0.0000000000000000E+000 - Pass - Time = 38.868 ns FracSkip2(-Inf) =  0.0000000000000000E+000 - Pass - Time = 33.819 ns

****

Gareth aka. Kit

Reply by email

J. Gareth Moreton

As an extra point, removing the 'skip' check (i.e. cmp ax, $3FE0, jbe @@skip) removes 6 bytes from the code size and shaves about 2 to 3 nanoseconds off the execution time in most cases, and it could be argued that it's worth going for the 'no skip' version because using Frac on a value of x where |x| < 1 is rather uncommon compared to when |x| >= 1.

However, when running my timing tests, one thing that's confused me is that when using very large inputs like 10^300, the function is at least 5 nanoseconds slower than FracSkip2, even though the code is less complex. This happens even if I put 'align 16' before the @@zero label.

I did wonder if it being a debug build caused some issues, but when I compiled it with full optimisation, both versions of the functions ran slower for numbers of that size (and the original FracDoSkip took about just as long), and SafeFrac beat them by around 5 nanoseconds.

Nevertheless, I conclude that for most situations, using the improved FracNoSkip gives the best performance and size for typical inputs, but this may depend on an individual machine's architecture.

****

function FracNoSkp2(const X: ValReal): ValReal; assembler; nostackframe; asm
  movq      rax,  xmm0
  shr       rax,  48
  and       ax,   $7FF0
  cmp       ax,   $4330
  jge       @@zero
  cvttsd2si rax,  xmm0
  cvtsi2sd  xmm4, rax
  subsd     xmm0, xmm4
  ret
@@zero:
  xorpd     xmm0, xmm0
end;

****

Note: 'align 16' at the start of a procedure is usually unnecessary, as FPC aligns procedures to 16-byte boundaries automatically.  FracNoSkp2 has a code size of 39 bytes, so will fill 48 bytes (3 blocks), which is a block smaller than the original FracNoSkip and the current Frac function.

I've attached my test project to this e-mail if you wish to look at the figures yourselves (I hope attachments work) and make a more informed decision.  This will currenly only run on Windows due to the use of QueryPerformanceCounter for timing checks. These calls will need to be removed to run this on Linux.

Gareth aka. Kit

Reply by email

Thorsten Engler
to J. Gareth Moreton
Show quoted text 8 lines
From: fpc-devel On Behalf Of J. Gareth Moreton Sent: Sunday, 29 April 2018 12:36 As an extra point, removing the 'skip' check (i.e. cmp ax, $3FE0, jbe @@skip) removes 6 bytes from the code size and shaves about 2 to 3 nanoseconds off the execution time in most cases, and it could be argued that it's worth going for the 'no skip' version because using Frac on a value of x where |x| < 1 is rather uncommon compared to when |x| >= 1.

I agree that calling Frac on values that are already just a fraction is probably not going to happen too often.

However, when running my timing tests, one thing that's confused me is that when using very large inputs like 10^300, the function is at least 5 nanoseconds slower than FracSkip2, even though the code is less complex. This happens even if I put 'align 16' before the @@zero label.

I do not see any noticeable difference between 1e16 and 1e300 as inputs:

Code address:
Frac1: 0000000000536440 (64)
Frac2: 0000000000536490 (16)
Frac3: 00000000005364E0 (96)
Frac4: 0000000000536530 (48)
Frac5: 0000000000536580 (0)
Frac6: 00000000005365D0 (80)
Frac7: 0000000000536620 (32)
Frac8: 0000000000536670 (112)

1st run:
In range (1e15+0.5):
Frac1 923470
Frac2 964422
Frac3 967501
Frac4 1027080
Frac5 1005352
Frac6 1052105
Frac7 1011983
Frac8 1048743

Out of range (1e16+0.5):
Frac1 893526
Frac2 998532
Frac3 894644
Frac4 993987
Frac5 895353
Frac6 994606
Frac7 900848
Frac8 992751

Out of range (1e300):
Frac1 897274
Frac2 986679
Frac3 899123
Frac4 999495
Frac5 899438
Frac6 989588
Frac7 885060
Frac8 985288

Only fraction (0.5):
Frac1 954220
Frac2 1046781
Frac3 993959
Frac4 1015032
Frac5 1013128
Frac6 1043157
Frac7 928712
Frac8 988220

Also, it seems to be relatively resilient against changes in code alignment even if it's not a multiple of 16:

Code address:
Frac1: 0000000000536433 (51)
Frac2: 000000000053645D (93)
Frac3: 0000000000536487 (7)
Frac4: 00000000005364B1 (49)
Frac5: 00000000005364DB (91)
Frac6: 0000000000536505 (5)
Frac7: 000000000053652F (47)
Frac8: 0000000000536559 (89)

1st run:
In range (1e15+0.5):
Frac1 946247
Frac2 904187
Frac3 902870
Frac4 1025163
Frac5 931021
Frac6 895990
Frac7 1050683
Frac8 952305

Out of range (1e16+0.5):
Frac1 883588
Frac2 877412
Frac3 809785
Frac4 831095
Frac5 976555
Frac6 711201
Frac7 791657
Frac8 897085

Out of range (1e300):
Frac1 902103
Frac2 901861
Frac3 802404
Frac4 808002
Frac5 972999
Frac6 710888
Frac7 804050
Frac8 875901

Only fraction (0.5):
Frac1 945212
Frac2 904468
Frac3 915325
Frac4 997584
Frac5 945569
Frac6 898036
Frac7 1071561
Frac8 906152

Nevertheless, I conclude that for most situations, using the improved FracNoSkip gives the best performance and size for typical inputs, but this may depend on an individual machine's architecture.

Seems we got a winner.

I was considering the ret like that, but didn't do it as I was worried because SEH under windows expects function prologues and epilogues that exactly match a specific pattern. But in hindsight, this is a no stack frame leaf function anyway, so I don't think that matters.

Cheers,
Thorsten

Reply by email

Florian Klaempfl
to Thorsten Engler: That ended up making things worse in some cases. Here is a branchless version: function…
Thorsten Engler wrote 14 lines
> -----Original Message----- > From: fpc-devel On Behalf > Of Florian Klämpfl > So something like > > cmp edx, $43300000 > jge @@zero > cmp edx, $3FE00000 > .align 16 > jbe @@skip > > might be much better. That ended up making things worse in some cases.

Can you take a look at the generated machine code if delphi uses proper multi byte nops. If not, the align might make things indeed worse.

Reply by email

Thorsten Engler
to Florian Klaempfl
Show quoted text 10 lines
-----Original Message----- From: fpc-devel On Behalf Of Florian Klaempfl Sent: Monday, 30 April 2018 04:28 > That ended up making things worse in some cases. Can you take a look at the generated machine code if delphi uses proper multi byte nops. If not, the align might make things indeed worse.

It does.

The problem was not the time required by the nops, but that for certain entry point alignments (among them the 16 byte alignments) the presence of this .align triggered the 3-4 times increase in processing time. I didn't look any closer into it as the version that J. Gareth worked out is faster and isn't alignment sensitive.

Reply by email