To improve efficiency, you could consider partitionning your problem:
Decompose your area of interest into several zones, run the process on these zones and re-unite the result.
For your case, you could consider decomposing the input square into 4, for each of them, select the land polygons intesecting them (with a spatial index why not), running the union/...
An option to use difference function with point geometries is using QGIS expressions with Geometry generator or Geometry by expression (see here for details) with this expression:
aggregate('points_2', 'collect', $geometry
Expression used on layer points_1 (10 red dots) with difference from points_2 (6 blue triangles) ...
No real solution has been found yet, but I got a workaround that worked quite nicely.
I got rid of slivers around polygons in QGIS after a Difference tool by first buffering the feature by 5cm. For my case a slight inaccuracy in area was less important than getting rid of the "Difference slivers".
The Difference tool also created some unwanted ...