QGIS area calculation differs when on the fly CRS transformation enabled

I can confirm that it seems to be a bug.

Create a csv file with the following content:

E N
600000 200000
700000 200000
700000 300000
600000 300000

Import it as delimited text with EPSG:21781, enable snapping and draw a polygon shapefile on the four points.

Without OTF, the result of $area/1000000.0 is 10000 m² (which is obviously correct).

Turning OTF on, and selecting the same EPSG:21781, you get 9988.2338 m².

Choosing a different CRS, like EPSG:4326, delivers 9990.5339 m², because the calculation is done on a different ellipsoid (WGS84 instead of bessel).

Vector --> Geometry Tools --> Export/Add Geometry Columns seems to deliver correct values.

The bug already has some tickets: https://issues.qgis.org/issues/10966 and https://issues.qgis.org/issues/12473

EDIT - Disclaimer: I would like to refer the readers to the discussion with ChrisW below. It might be that getting an area based upon an OTF CRS is not a bug after all; that is, at least, in arcgis it also being used to allow geoprocessing two layers from different CRS.

To elaborate on the issue above. As AndreJ as suggested and show - this is probably a bug in qgis's current version. Yet it should be noted that the problem isn't the wrong area, but that on-the-fly transformation affects in anyway on area calculations.

The purpose of on-the-fly transformation/projection is to align data from different sources and with different CRS. That is mainly for display purposes. E.G. arcmap automatically perform on-the-fly projection in any case a layer CRS does not match the data frame CRS.

Arcmap also provides a possibility to edit data while projected on-the-fly, but also notes that: (source)

However, it is important to note that certain editing operations may produce unexpected alignment or accuracy problems, depending on the coordinate systems being used.

Specific editing operations that may cause issues include changing the shapes of features, snapping to the edge or boundary of features, or extending and trimming features. These problems are more likely to occur when the features you are editing are close to the edge or beyond the area of use of the coordinate system

That is to say: on the fly transformation is less acurate than just project the data to a different CRS (which also introduces its own problems).

Having said that it is not surprising that based on a on-the-fly transformation a wrong area is being calculated, yet it is surprising that the fact that on-the-fly was enabled affects in any way the calculation of geometry, which should be based on the data. Thus, it doesn't matter if on-the-fly transformation is based upon the same or a different CRS, area calculation should be identical each time.

To be more practical, if your aim is to compute the area do not use on-the-fly. If you have the wrong CRS, project your data.