Implementation of the Simpson's rule and the length of an Ellipse(-arc)

Discussion forum for C++ and script developers who are using the QCAD development platform or who are looking to contribute to QCAD (translations, documentation, etc).

Moderator: andrew

Forum rules

Always indicate your operating system and QCAD version.

Attach drawing files, scripts and screenshots.

Post one question per topic.

Post Reply
CVH
Premier Member
Posts: 4917
Joined: Wed Sep 27, 2017 4:17 pm

Implementation of the Simpson's rule and the length of an Ellipse(-arc)

Post by CVH » Mon Sep 22, 2025 7:34 am

Andrew, all,

REllipse::getSimpsonLength(a, b) in REllipse.cpp (Current) is basically the implementation of Calculation of Ellipse Arc Length found online.

Both methods use 20 intervals.
As per Simpson's rule the factor q is 1.0 for the first and last interval, 4.0 when even, otherwise a factor of 2.0 is used.
Minor difference: In QCAD the interval designation is zero based, in the linked reference it starts with 1, odd/even interval factors are swapped.

An important note for the Java applet is Scope : First Quadrant, 0 x R
Where R is the length of the Semi-Axis lying on the X-axis aka the Major Radius.

Another important fact is that yi is a positive root, the sum of the factored y values is positive thus the df / 3.0 part must be positive.
Meaning that df must be positive and that a2 > a1 so that df = (a2-a1)/20 is larger than zero.

The use of getSimpsonLength(a, b) is limited to REllipse::getLength() for returning an estimate of the length of the ellipse(-arc).
It is thus important that the implementation in getLength() complies with the above.


Studied the implementation in depth knowing that the length of an Ellipse-arc may fail in some situations (FS#2693 >> FS#2565).
:arrow: It doesn't comply with the first quadrant limit.
:arrow: Estimates of a half full ellipse circumference are used, then also over no more than 20 intervals.
:arrow: It fails with for example: An ellipse starting at 45°, passing 0° in CW sense and ends at (-)180°... =2Pi is the issue here.

:!: Not that it matters a lot to know the exact length of an ellipse-arc but when returning NaN drawing tools like D2 may fail :!:

4 years ago a 'faster' algorithm was implemented for the circumference of a full ellipse. ('faster' is debatable, elegant too :wink: )
Not referenced but the Padé 3/3 approximation is said to be very practical and more efficient than the earlier used Ramanujan II approximation.
Padé 3/3 doesn't use a root and the overall accuracy is increased by a 1.5 notch on the log scale. (Review of known formulae)
My 2 main reasons for proposing it.

It is not the limit of what can be achieved.
The Cantrell-Ramanujan approximation (C6) has for example a maximum error of 14.5ppm.
Even better approximations exist at the cost of complexity. (e.g. 152ppb Stan's Library)
(Also see notes in wire Centroids but these calculate a nearly exact full circumference by well suited infinite series expansion.)


Knowing the 'almost' exact circumference of a full ellipse one can easily figure out a half or a quarter ellipse circumference.
Because it is symmetrical in regards with the X and/or Y-axis simply divide for example the Padé 3/3 approximation by 2 or 4.

Then we could chop up an arbitrary ellipse-arc and estimate parts or remainders in ... or mirrored to the first quadrant.
Our best option is then to estimate the more flatter part and subtract that from a quarter ellipse circumference if required.
Already limited to call getSimpsonLength(a, b) at most twice in combination with the 'fast' method for a full ellipse.
It will always be an estimated length without implementing elliptic integrals of the first and second kind.


Next would be to implement all this under QCAD in ECMAScript and report back on the success. :wink:
I can then compare the summed length of some arbitrary ellipse-arcs with the nearly exact full circumference from Centroids.
e.g. 3x 1/3 full ellipse-arcs should match with a full ellipse and so on.

Regards,
CVH

CVH
Premier Member
Posts: 4917
Joined: Wed Sep 27, 2017 4:17 pm

Re: Implementation of the Simpson's rule and the length of an Ellipse(-arc)

Post by CVH » Tue Sep 23, 2025 8:36 am

All,

I can already share my first observation concerning the estimation of half the full ellipse circumference.
In detail the use of getSimpsonLength(0.0, Pi) or getSimpsonLength(Pi, 2Pi) in REllipse::getSimpsonLength(a, b).
(Second main issue)

Rather good for nearly round ellipses (ratio ≈1.0 or Eccentricity near zero).
The more the ellipse is elongated, flatter, the poorer the result.

2 times a quarter ellipse (getSimpsonLength(0.0, Pi/2)) is a far better estimation for a half ellipse.
In the mid range with for example ratio = 0.5 it is over 150000 times more accurate.
But that rapidly declines to 6.5 times more accurate for a 0.05 ratio or when almost flat.

For rather elongated ellipses the reason must be obvious.
We are applying the Simpson's rule on sections with high to very high curvature near the ellipse major points.
For a half ellipse twice near both major points and that with coarser steps than 2 times the quarter ellipse approach.

Then 1/4 or 1/2 of the Padé 3/3 approximation of a full ellipse is always more accurate.


For the record: Although never attempted to estimate for in REllipse::getSimpsonLength(a, b)
An ellipse length from Pi/2 to 3Pi/2 exhibits similar poor results as those found for 0-Pi or Pi-2Pi.
:arrow: What supports the statement to avoid the extremities and rather estimate the flatter part if convenient.

Regards,
CVH

Post Reply

Return to “QCAD Programming, Script Programming and Contributing”