View Issue Details

IDProjectCategoryView StatusLast Update
0036415FPCCompilerpublic2019-12-23 13:15
ReporterTomas Ledvinka Assigned ToFlorian  
PrioritynormalSeverityminorReproducibilityalways
Status resolvedResolutionwon't fix 
Product Version3.0.4 
Summary0036415: certain divisions in single precision only
DescriptionDear fpc developers,

Let me submit a divison bug:

When a floating point constant which can be represented exactly in single precision is divided by integer variable, the result is assumed to be a single-precision number.

program div5;
var five:integer=5;
begin
    writeln(1.0/five);
    writeln(1.0/five-1/five);
end.

--- OUTPUT ---
 2.000000030E-01
 2.9802322276673010E-009
--------------------

versions where problem spotted:
3.0.4 [2019/10/27] for x86_64 on windows
3.0.4+dfsg-18ubuntu2 [2018/08/29] for x86_64

older version with no problem:
fpc 2.6.4 i386-win32 is OK

Thanks
Steps To Reproducesee the program above
Additional Informationdisassembly shows use of single precision arithmetics

movswl TC_$P$HELLOWORLD_$$_FIVE,%eax ; five is 16-bit integer type
cvtsi2ss %eax,%xmm0 ;
movss _$HELLOWORLD$_Ld1,%xmm1 ; OK, constant 1.0 fits into single precision exactly
divss %xmm0,%xmm1 ; NOT OK (only single precision arithmetics used)
movss %xmm1,-20(%rbp) ; NOT OK
TagsNo tags attached.
Fixed in Revision
FPCOldBugId
FPCTarget-
Attached Files

Activities

Florian

2019-12-09 21:39

administrator   ~0119717

See mailing list discussion why it is handled as it is. There is no reason why <single> op <integer> should result in another type than single as the range of single includes integer.

Tomas Ledvinka

2019-12-10 11:11

reporter   ~0119725

Last edited: 2019-12-10 11:13

View 2 revisions

Please, please, check the following program and its output and, please, do not dismiss this issue immediately:

program fourdiv;
var five : integer = 5;
     one : real = 1.0;
begin
    writeln(1.0 / five);
    writeln(one / five);
    writeln( 1 / five);
    writeln( 1 / 5);
end.

It should output four mosly identical lines, yet the first one is now very different:

 2.000000030E-01 <- too much inexact (see digit three inside)
 2.0000000000000001E-001
 2.0000000000000001E-001
 2.0000000000000000E-001

Due to the nature of division operation already integer/integer yields double in Pascal (while +-* behave differently) and the behavior on the first Writeln is too irregular.

Thank you

  T.

Florian

2019-12-10 18:55

administrator   ~0119732

Yes, this is excatly what has been described on the mailing list and perfectly fine.

Marco van de Voort

2019-12-10 20:02

manager   ~0119733

You are doing a floating point divide (/). Floating point always rounds somewhere, depending on storage format and alu hardware.

Tomas Ledvinka

2019-12-11 12:57

reporter   ~0119748

Last edited: 2019-12-11 13:03

View 3 revisions

(Please note, that I did not succeed to register in mailing list -- no reply, no error -- and thus I use the web interface).

Of course that float point has to round somewhere but there is no obvious reason why 1.0/five should be inferior to 1/five. (It was OK in FPC 2.6.4 and update to 3.0.4 broke my code.)

I do use Pascal for teaching/computing with sometimes long formulas and It is very usual in numeric code to write 1.0 even if integer 1 is OK (see code on page 2 of 1234 of well-known textbook Press et al: Numerical Recipes).

Compiler user would not suspect that using 1.0 instead of 1 will ruin his code:

type tMatrix = array of array of real;
...
function HilbertMatrix(n:integer): tMatrix;
var r,s:integer;
    M:tMatrix;
begin
  SetLength(M,n,n);
  for r:=0 to n-1 do
       for s:=0 to n-1 do
           M[r][s] := 1.0/(1+r+s);
  
  HilbertMatrix := M;
end;

begin
  Writeln( 1 / Determinant( HilbertMatrix(4) ) : -12 : 6 );
end.

and would spend long time debugging the Determinant function to find out why his code outputs

6047865.901211

instead of

6048000.000000

Thank you for considering my suggestions.

T.

P.S. And it is even more usual to write 0.5 instead of 1/2. I think that any use of single precision by compiler which impacts precision of the result should be "on demand" only.

Florian

2019-12-11 18:36

administrator   ~0119756

> I think that any use of single precision by compiler which impacts precision of the result should be "on demand" only.

No. People expect that the compiler does not waste computational resources, in particular as double divisions are more expensive than single divisions. 1.0 fits perfectly into a single. So single is sufficient. There is no reason why it should double, extended, float128 or even a bignum. In particular because there are architectures which support only single. If you do not care about computional resources, just use -CF64 (or -CF80 on i386) in your configuration file.

The behaviour of older versions of FPC was governed by the x87 behaviour which does operations with extended precision if not configured otherwise, but the x87 design is unique and not replicated by other FPUs.

Cyrax

2019-12-11 19:42

reporter   ~0119760

Last edited: 2019-12-11 19:57

View 2 revisions

By using test program at note 0119725 as an base (see attached test project), compiling it with FPC trunk r43670 i386-linux gives this output:
 2.000000000E-01
 2.0000000000000000E-001
 2.00000000000000000003E-0001
 2.00000000000000000003E-0001


EDIT: it gives same output when running at i386-win32 target.

project1.zip (1,171 bytes)

Cyrax

2019-12-11 19:44

reporter   ~0119761

@Florian, FPC compiler complains if I use -CF80 option (Error: (11006) Illegal parameter: -CF80). Is this really supported on i386-linux target?

Martok

2019-12-11 19:50

reporter   ~0119762

I guess it comes down to the "principle of least surprise" thing again: any other pascal compiler and all of the C world defaults to double (or even extended) unless specifically reduced, so it's a bit of a stretch to say people would reasonably expect FPC to do the opposite... the only exception I can think of is Fortan, which I'd take as a sign that something is not a good idea.


> If you do not care about computional resources, just use -CF64 (or -CF80 on i386) in your configuration file.

"Error: Illegal parameter: -CF80"
So there is no way to force extended precision, at the very least.
A shorthand for "highest available on target, error if host can't do that" would be incredibly useful, since that is basically what Delphi does.
Speaking of which: could these be automatically added when compiling in $mode delphi, where no user will expect FPC trying to be clever?

Cyrax

2019-12-11 19:59

reporter   ~0119763

x86_64-linux gives this output:
 2.0000000000000001E-001
 2.0000000000000001E-001
 2.0000000000000001E-001
 2.0000000000000001E-001

Cyrax

2019-12-11 20:02

reporter   ~0119764

x64_64-win64 gives this:
2.000000030E-01
 2.0000000000000001E-001
 2.0000000000000001E-001
 2.0000000000000001E-001

Florian

2019-12-11 20:24

administrator   ~0119766

> I guess it comes down to the "principle of least surprise" thing again: any other pascal compiler and all of the C world defaults to double

I would find it very surprising if 1.0 is parsed as a double while it fits into a single. Next bug report will be then about 255 being handled as byte but people expect it to be an UInt128 for whatever reason.

Martok

2019-12-11 21:48

reporter   ~0119772

That's what a standard is for, in the absence of that, documentation, in the absence of that, convention.

Three observations, in reverse order.
a) Taking the literal "255" as type Byte instead of Integer is already Borland-incompatible: """Other decimal numbers [than reals] denote integer-type constants; they must be within the range -2,147,483,648 to 2,147,483,647.""" But at least it's somewhat consistent within FPC (but undocumented AFAICS), and not relevant to this issue anyway.
b) We're not really talking about parsing here, are we? The issue here is a low-precision intermediate giving a low-precision result... sometimes, and not before two years ago.
c) The typeconvs here are a bit mad anyway.
    writeln(1.0 / five);     // Single / Integer  = Single(2.000000000E-01 on win32, 2.000000030E-01 on win64)
    writeln(one / five);     // Double / Integer  = Double
    writeln( 1 / five);      // Integer / Integer = Extended
    writeln( 1 / 5);         // Integer / Integer (compiletime eval) = Extended   
    writeln( 1.0 / 5.0);     // Single / Single (compiletime eval) = Single(2.000000030E-01)

I would at the very least not expect all of them to be true at once.

Florian

2019-12-11 22:15

administrator   ~0119773

I will not discuss further here because: 1) i386, x86_64-win64 and x86_64-linux behave different due to different FPU modells (x87, SSE, x87+SSE) and this is not taken into account in the discussion; 2) I see no need to change anything, to me it makes perfectly sense. And yes, the decision on the type is made in the parser.

Sven Barth

2019-12-12 12:35

manager   ~0119784

There is nothing mad about the conversions:
- a Single or Double divided by an Integer stays a Single or Double
- a Single divided by a Single stays a Single
- a floating point division of an Integer results in a floating point value; this uses the highest available precision.

Tomas Ledvinka

2019-12-12 18:02

reporter   ~0119800

Let me thank everyone for their arguments.

> i386, x86_64-win64 and x86_64-linux behave different due to different FPU modells (x87, SSE, x87+SSE) and this is not taken into account in the discussion;
Let me suggest an FPU-independent rule:  Single/Integer uses the same arithmetic precision as Integer/Integer (because there are many exact integers in single type).

In more detail: Assume we have a float point literal X which can be exactly represented as rational number P/Q (Q is integer power of two and P is integer). Then the expression with integer denominator K in the form X/K should be evaluated in the same precision as expression P/Q/K. (Yes it would be double in my case and I would be happy.)

>I would find it very surprising if 1.0 is parsed as a double while it fits into a single.

Of course, this is OK, since 1.0 fits exactly into single.

>Next bug report will be then about 255 being handled as byte

Again, no problem to parse 255 as byte literal. If addition instead of division is used, the analogy shows the problem is in the bit-width of the arithmetic operation:
Imagine, the expression (255+one) would evaluate to zero because a byte-wide arithmetic would be used for addition. But your Pascal is kind, so it evaluates to 256 whatever integer type of variable "one" is used.

> a Single or Double divided by an Integer stays a Single or Double
Let's assume byte instead of single and plus instead of divide and Imagine how uncomfortable it would be to have rule, which for "255+one" would provide 0, while "256+one" would be 257.

In other words the division operation has a special feature (precision loss) and it should be allowed to return single only if "single / single" division is performed. (A feature of addition is overflow and it is considered by FPC not insisting on the byte type of the result in the example above.)

To demonstrate that the problem is not in the parser which understands 1.0 as single, see the following code, where the second line gives the right result:

    writeln(1.0/five);
    writeln(1.0*(1/five));

Here the integer division "1/five" produces double and multiplication by 1.0 works perfectly even when 1.0 is parsed as single. (Remember: 1.0 is stored in single because of efficiency not because it is inexact.)

Summarized, the division is a special arithmetic operation and thus in the expression "1.0/five" should not be treated worse that 1/five. This can be achieved converting first, second or both operands to double, since FPC now seems to treat "singleA op doubleB" as "double(singleA) op doubleB".

nanobit

2019-12-13 11:07

reporter   ~0119811

On fpc3.2, {$EXCESSPRECISION ON} can be used to get more precision (extended in win32).
But my test shows, some cases still remain single ((1.0 / 5) is single, (1.0 / 5.0) is extended).

Florian

2019-12-22 22:30

administrator   ~0120031

Sven pointed out the rules. They are very simple, waste no computational resources, are cross platform, easy to understand and if higher precision is needed, they can be worked around. Any other handling adds more exceptions and makes things harder to understand.

Issue History

Date Modified Username Field Change
2019-12-09 17:08 Tomas Ledvinka New Issue
2019-12-09 21:39 Florian Assigned To => Florian
2019-12-09 21:39 Florian Status new => resolved
2019-12-09 21:39 Florian Resolution open => no change required
2019-12-09 21:39 Florian FPCTarget => -
2019-12-09 21:39 Florian Note Added: 0119717
2019-12-10 11:11 Tomas Ledvinka Status resolved => feedback
2019-12-10 11:11 Tomas Ledvinka Resolution no change required => reopened
2019-12-10 11:11 Tomas Ledvinka Note Added: 0119725
2019-12-10 11:13 Tomas Ledvinka Note Edited: 0119725 View Revisions
2019-12-10 18:55 Florian Note Added: 0119732
2019-12-10 18:55 Florian Assigned To Florian =>
2019-12-10 20:02 Marco van de Voort Note Added: 0119733
2019-12-11 12:57 Tomas Ledvinka Note Added: 0119748
2019-12-11 12:57 Tomas Ledvinka Status feedback => new
2019-12-11 12:59 Tomas Ledvinka Note Edited: 0119748 View Revisions
2019-12-11 13:03 Tomas Ledvinka Note Edited: 0119748 View Revisions
2019-12-11 18:36 Florian Note Added: 0119756
2019-12-11 19:42 Cyrax File Added: project1.zip
2019-12-11 19:42 Cyrax Note Added: 0119760
2019-12-11 19:44 Cyrax Note Added: 0119761
2019-12-11 19:50 Martok Note Added: 0119762
2019-12-11 19:57 Cyrax Note Edited: 0119760 View Revisions
2019-12-11 19:59 Cyrax Note Added: 0119763
2019-12-11 20:02 Cyrax Note Added: 0119764
2019-12-11 20:24 Florian Note Added: 0119766
2019-12-11 21:48 Martok Note Added: 0119772
2019-12-11 22:15 Florian Note Added: 0119773
2019-12-12 12:35 Sven Barth Note Added: 0119784
2019-12-12 18:02 Tomas Ledvinka Note Added: 0119800
2019-12-13 11:07 nanobit Note Added: 0119811
2019-12-22 22:30 Florian Assigned To => Florian
2019-12-22 22:30 Florian Status new => resolved
2019-12-22 22:30 Florian Resolution reopened => won't fix
2019-12-22 22:30 Florian Note Added: 0120031