Q:

Should be 1.001 us equal to 1001 ns in VHDL?

I'm currently writing a test to check whether the type time is synthesized / simulated properly in various FPGA vendor tools. One corner case is the usage of real literals as the abstract literal for time values, e.g.: 1.001 us.
The IEEE Std. 1076-2008, Section 5.2.4.1 or IEEE Std. 1076-1993, Section 3.1.3, Para 8 states:

There is a position number corresponding to each value of a physical
  type. The position number of the value corresponding to a unit name is
  the number of primary units represented by that unit name. The
  position number of the value corresponding to a physical literal with
  an abstract literal part is the largest integer that is not greater
  than the product of the value of the abstract literal and the position
  number of the accompanying unit name.

As per definition of TIME (see Section 16.2 or 14.2), the primary unit is femto seconds, so that, the position number of 1 us is 1,000,000,000. Thus, the position number of 1.001 us should be 1,001,000,000 which is also the position number of 1001 ns. Thus, both values should be equal.
But, I get different results when I try to synthesize or simulate the following stripped down unit. I have checked, that the minimum time resolution is either 1 fs or 1 ps.
entity physical_test is
    generic (
        C1 : time := 1001 ns;
        C2 : time := 1.001 us;
        C3 : time := TIME'val(integer(real(TIME'pos(1 us)) * 1.001))
    );
    port (
        y : out bit);
end entity physical_test;

architecture rtl of physical_test is
    function f return boolean is
    begin
        report "C1 = " & TIME'image(C1) severity note;
        report "C2 = " & TIME'image(C2) severity note;
        report "C3 = " & TIME'image(C3) severity note;
        return false;
    end f;

    constant C : boolean := f;
begin  -- architecture rtl
    y <= '0';
end architecture rtl;

QuestaSim (ModelSim) was the only tool which reported the expected result:
# ** Note: C1 = 1001000000 fs
# ** Note: C2 = 1001000000 fs
# ** Note: C3 = 1001000000 fs

But, the actual result when synthesizing with Quartus 15.0 or ISE 14.7 is:
Note: "C1 = 1001000000 fs"
Note: "C2 = 1000999999 fs"
Note: "C3 = 1001000000 fs"

Thus, the value of C2 is not as expected. If I write down the cited text as an equation, then I get the expected result in the constant C3.
When I use the integrated simulators of ISE 14.7 or Vivado 2015.4, I get a similar result:
Note: "C1 = 1001000 ps"
Note: "C2 = 1000999 ps"
Note: "C3 = 1001000 ps"

So, should be the behaviour of Quartus / ISE / Vivado considered as a bug? Or is it allowed by the VHDL standard that 1.001 us is not equal to 1001 ns?
EDIT: The bug happens also when I compare 1.001 ps to 1001 fs as well as when I compare 1.001 ns to 1001 ps. As also the manual computation with C3 was correct, it should not be a problem with the accuracy of real.
Just to note, the synthesizer of Vivado reports weird results:
Parameter C1 bound to: 32'b00111011101010100000110001000000    -- 1001000000
Parameter C2 bound to: 32'b10010011011101001011110001101010    -- 2473901162
Parameter C3 bound to: 32'sb00000000000000000000000000000001   -- 1

A:

This is an issue with floating point numbers and has very little to do VHDL. Integer decimal numbers can be converted to binary and back to decimals without loosing information (unless the numbers are too big or too small). This is not the case for decimal fractions. The decimal number 1.001 converted to binary is irrational. When converting back from binary to decimal, there will be rounding errors. 
Quartus and ISE show expected behaviour.
In the Vivado case something happened with the MSB of C2. It appears that some conversion between signed and unsigned integers took place, where it should not have. And C3 was apparently rounded of.
Your examples could be used to support the rule that using the primary unit is the best choice.

