Rounding Errors
In the OpenSCAD manual, rounding is described as
round
The "round" operator returns the greatest or least integer part,
respectively, if the numeric input is positive or negative.
This rounding practice is not in line with IEEE 754
https://en.wikipedia.org/wiki/IEEE_754-1985#Rounding_floating-point_numbers
. There is a phrase overlooked, and thus not implemented, by OpenSCAD: It
should read: "...if the numeric input is positive or negative, so that the
least significant digit is an even number."
The IEEE754 correct rounding result is round(7/2)=4 and round(9/2)=4,
whereas OpenSCAD does round(9/2)=5 !
IEEE 754 choices have been made to result in the least bias when dealing
with aggregates of rounded numbers, a smaller bias than what OpenSCAD
achieves.
Because of rounding errors, the following code produces rather odd shapes
whenever the value of "Resolution" is not even:
Radius=1.5; // radius of Sphere
Resolution=7;
function Sphere(Theta,Phi) = Radius*[cos(Theta)sin(Phi),
sin(Theta)sin(Phi), cos(Phi)]; // a simple sphere
// store all information needed by the computer to create the shape in one
suitably ordered list as a stack of slices!
//AnyShape=(let (I_Step=360/Resolution, K_Step=180/round(Resolution/2)) [
for (k=[0 : 1 : round(Resolution/2)] ) [ for( i=[1 : 1 : Resolution] )
Sphere(iI_Step,kK_Step) ] ] );
AnyShape=(let (Step=360/Resolution) [ for (k=[0 : Step : 180] ) [ for(
i=[Step :Step : 360] ) Sphere(i,k) ] ] );
// Start shape generating list
Points= [ for (a=AnyShape) for (b=a) b ]; // to obtain the
list of points, remove one layer of brackets from "AnyShape"
Top=TopCone();
Side1=( [for( i=[1:1:len(AnyShape)-3] , k=[0:1:len(AnyShape[i])-1 ] ) let
(S=len(AnyShape[i])) [ (i+1)S+k, iS+k , i*S+((k+1) % S) ] ] );
Side2=( [for( i=[1:1:len(AnyShape)-3] , k=[0:1:len(AnyShape[i])-1 ] ) let
(S=len(AnyShape[i])) [ (i+1)S+k, iS+((k+1) % S), (i+1)*S+((k+1) % S) ] ]
);
Bottom=BottomCone();
function TopCone() = (let (S=len(AnyShape[1])) [for( i=[0:1:S-1] ) [ 0,
S+((i+1) % S), S+i ] ] );
function BottomCone() = (let (L=len(AnyShape), S=len(AnyShape[L-1])) [for(
i=[0:1:S-1] ) [ (L-2)*S+i, (L-2)S+((i+1) % S), LS-1 ] ] );
Faces=concat(Top,Side1,Side2,Bottom); // list of all faces needed to
close shape
// End shape generating list
polyhedron(Points,Faces);
The moment the rather naive implementation of the for loop in line 6, which
uses floating point numbers, is replaced by a for loop that enforces the use
of integers instead (line 5) the problem goes away.
Doug.moen has proposed
http://forum.openscad.org/for-statement-doesn-t-do-all-the-steps-sometimes-td19591.html#a19624
another approach. My question is: Does that approach produce a proper sphere
using line 6? and, Is it being implemented? Is it needful to change the
implementation of the for loop if rounding is done as per IEEE 754?
By the way, I am using OpenSCAD 2015.03-1 and Kubuntu.
wolf
--
View this message in context: http://forum.openscad.org/Rounding-Errors-tp21821.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
What you describe is the Anglo-American way of rounding numbers. It makes
sense if you do not count in decimals, but in half-inches, quarter-inches,
etc.
In Europe, where we rather work with factors of 10 instead of factors of 2,
absolute values of .5 and higher are rounded up and the rest down. I think
that computationally, this is a more fair way to do things, as the outcome
follows from the algorithm instead of the other way around.
Some programming languages have an extra parameter to let the programmer
decide what method to use in rounding.
I can imagine that for OpenSCAD, this would be an enhancement, since you
cannot round a variable, look at the outcome, and then change that same
variable.
--
View this message in context: http://forum.openscad.org/Rounding-Errors-tp21821p21824.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
On 2017-07-14 03:47, wolf wrote:
In the OpenSCAD manual, rounding is described as
round
The "round" operator returns the greatest or least integer part,
respectively, if the numeric input is positive or negative.
This rounding practice is not in line with IEEE 754
I believe the OpenSCAD built-in round function exposes the C/C++
function of the same name.
http://www.cplusplus.com/reference/cmath/round/
"Returns the integral value that is nearest to x, with halfway cases
rounded away from zero."
Carsten Arnholm
The C language provides 4 operations for converting an floating point value
to an integer:
floor -- round towards -infinity
ceil -- round towards +infinity
trunc -- round towards zero
round -- round away from zero
OpenSCAD provides the same operations.
None of these is "round towards even". Wolf's post and the Wikipedia
article on "Rounding" give some good reasons why this is a useful operation.
Some forum members raised objections to my proposed range semantics. I
think the underlying problem is that the proposal could change the
behaviour of existing OpenSCAD programs. Anyway, I'm not implementing the
proposal in OpenSCAD, due to these objections (although I did implement it
in my own geometry language).
If you want to play with my range semantics, try defining this function:
function drange(first,step,last) = let (
n = round((last - first) / step)
)
[first: step: first + stepn + step.1];
and then use drange(a,b,c) in place of [a:b:c].
I tried using drange in Wolf's code example, and the results do look better
for Resolution=7.
On 15 July 2017 at 03:43, arnholm@arnholm.org wrote:
On 2017-07-14 03:47, wolf wrote:
In the OpenSCAD manual, rounding is described as
round
The "round" operator returns the greatest or least integer part,
respectively, if the numeric input is positive or negative.
This rounding practice is not in line with IEEE 754
I believe the OpenSCAD built-in round function exposes the C/C++ function
of the same name.
http://www.cplusplus.com/reference/cmath/round/
"Returns the integral value that is nearest to x, with halfway cases
rounded away from zero."
Carsten Arnholm
OpenSCAD mailing list
Discuss@lists.openscad.org
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org
First, I need to correct myself. My original source for the rounding practice
discrepancy was a book," Randal E. Bryant and David R O'Halloran, Computer
Systems, A Programmer's Perspective, Third Edition, Pearson 2017", which
listed four rounding practices. OpenSCAD's round() was missing from them. In
the meantime, I have found another Wikipedia
https://en.wikipedia.org/wiki/IEEE_754#Rounding_rules reference, which
lists five rules, among them OpenSCAD's round().
Of round to even, it says:
"Round to nearest, ties to even – rounds to the nearest value; if the number
falls midway it is rounded to the nearest value with an even (zero) least
significant bit; this is the default for binary floating-point and the
recommended default for decimal."
The discrepancy is thus about default rules. Elsewhere I read that the
majority of hardware does implement round-to-even as default - as it should,
because of marketing pressures.
I have already contacted Randal Bryant for comment on why he omitted one of
the IEEE 754 rounding options, and I will report back his answer once I have
received it. Since Doug.Moen has provided most detail about where OpenSCAD's
rounding rules come from, may I ask him to report back to us
a) which compiler and C version is used to compile OpenSCAD, and
b) if round-to-even is not available from that compiler, when do the
compiler manufacturer's plan to provide it?
Should the latter enquiry be escalated to Marius Kintel as the official
OpenSCAD maintainer?
My interest in rounding rules arises from my research into the origin and
consequences of degenerate triangles. I have in this forum already reported
that rotations by 90 degree are a cause of OpenSCAD failing with a CGAL
assertion error. and there are other causes that are not so easy to pin
down. In executing this rotation, one step is to convert degrees to radians,
i.e. angle=90*PI/180, which is a source of rounding errors. From here, it is
easy to understand that OpenSCAD returns sin(90)=0.999... and not, as it
should, sin(90)=1.
Instead of maintaining tables of special values for e.g. the sine function,
a proper rounding rule may be enough to do the trick - easy to do if the
compiler supports it.
And because C11 - or ISO/IEC 9899:2011, to give it its full name, is three
years younger than IEEE 754-2008, it ought to make available round-to-even.
It only depends on whether compiler writers have seen the urgency to
implement it, and programmers to have learned it exists . . .
wolf
--
View this message in context: http://forum.openscad.org/Rounding-Errors-tp21821p21829.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
The rotate by 90 problem in OpenSCAD is because it uses the C versions of
the trig functions that take radians in the rotate module. It also has trig
functions that work in degrees that it exposes as the built in functions. I
fix the issue by using a user space rotate module.
module rotate(angle) // built-in rotate is inaccurate for 90
degrees, etc
{
a = len(angle) == undef ? [0, 0, angle] : angle;
cx = cos(a[0]);
cy = cos(a[1]);
cz = cos(a[2]);
sx = sin(a[0]);
sy = sin(a[1]);
sz = sin(a[2]);
multmatrix([
[ cy * cz, cz * sx * sy - cx * sz, cx * cz * sy + sx * sz, 0],
[ cy * sz, cx * cz + sx * sy * sz,-cz * sx + cx * sy * sz, 0],
[-sy, cy * sx, cx * cy, 0],
[ 0, 0, 0, 1]
]) children();
}
This works because sin(90) is 1 in OpenSCAD but not in C because it is
impossible to represent 90 degrees in radians with floating point.
OpenSCAD could be fixed to use its degree trig functions everywhere instead
of converting to radians and using the C functions.
On 16 July 2017 at 02:22, wolf wv99999@gmail.com wrote:
First, I need to correct myself. My original source for the rounding
practice
discrepancy was a book," Randal E. Bryant and David R O'Halloran, Computer
Systems, A Programmer's Perspective, Third Edition, Pearson 2017", which
listed four rounding practices. OpenSCAD's round() was missing from them.
In
the meantime, I have found another Wikipedia
https://en.wikipedia.org/wiki/IEEE_754#Rounding_rules reference, which
lists five rules, among them OpenSCAD's round().
Of round to even, it says:
"Round to nearest, ties to even – rounds to the nearest value; if the
number
falls midway it is rounded to the nearest value with an even (zero) least
significant bit; this is the default for binary floating-point and the
recommended default for decimal."
The discrepancy is thus about default rules. Elsewhere I read that the
majority of hardware does implement round-to-even as default - as it
should,
because of marketing pressures.
I have already contacted Randal Bryant for comment on why he omitted one of
the IEEE 754 rounding options, and I will report back his answer once I
have
received it. Since Doug.Moen has provided most detail about where
OpenSCAD's
rounding rules come from, may I ask him to report back to us
a) which compiler and C version is used to compile OpenSCAD, and
b) if round-to-even is not available from that compiler, when do the
compiler manufacturer's plan to provide it?
Should the latter enquiry be escalated to Marius Kintel as the official
OpenSCAD maintainer?
My interest in rounding rules arises from my research into the origin and
consequences of degenerate triangles. I have in this forum already reported
that rotations by 90 degree are a cause of OpenSCAD failing with a CGAL
assertion error. and there are other causes that are not so easy to pin
down. In executing this rotation, one step is to convert degrees to
radians,
i.e. angle=90*PI/180, which is a source of rounding errors. From here, it
is
easy to understand that OpenSCAD returns sin(90)=0.999... and not, as it
should, sin(90)=1.
Instead of maintaining tables of special values for e.g. the sine function,
a proper rounding rule may be enough to do the trick - easy to do if the
compiler supports it.
And because C11 - or ISO/IEC 9899:2011, to give it its full name, is three
years younger than IEEE 754-2008, it ought to make available round-to-even.
It only depends on whether compiler writers have seen the urgency to
implement it, and programmers to have learned it exists . . .
wolf
--
View this message in context: http://forum.openscad.org/
Rounding-Errors-tp21821p21829.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
OpenSCAD mailing list
Discuss@lists.openscad.org
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org
To report back what I learnt:
Software rounding:
The GCC library libc provides rounding modes for all IEEE754-2008 rounding
functions, not just the four doug.moen listed. There is the page
https://www.gnu.org/software/libc/manual/html_node/Rounding.html on
rounding modes, and there is a page
https://www.gnu.org/software/libc/manual/html_node/Rounding-Functions.html#Rounding-Functions
on rounding functions. Round-to-even is available from function rint, via
setting a rounding mode, and also from function roundeven.
Hardware rounding:
Intel hardware
<https://software.intel.com/sites/default/files/managed/39/c5/325462-sdm-vol
-1-2abcd-3abcd.pdf> , page 106, table 4.8 rounds to even as default. It
does not support OpenSCAD-type rounding at all. The reason for that is that
round-to-even average rounding error is only half that of OpenSCAD-type
rounding, and Intel ties to provide accurate results whenever the number
format permits it.
With regard to doug.moen's drange function, here are the results:
http://forum.openscad.org/file/n21934/Sphere-1.png using
AnyShape=(let (Step=360/Resolution) [ for (k=[0 : Step : 180] ) [ for(
i=[Step :Step : 360] ) Sphere(i,k) ] ] );
http://forum.openscad.org/file/n21934/Sphere-2.png using
AnyShape=(let (Step=360/Resolution) [ for (k=drange(0,Step,180)) [ for(
i=drange(Step,Step,360) ) Sphere(i,k) ] ] );
http://forum.openscad.org/file/n21934/Sphere-3.png using
AnyShape=(let (I_Step=360/Resolution, K_Step=180/round(Resolution/2)) [ for
(k=[0 : 1 : round(Resolution/2)] ) [ for( i=[1 : 1 : Resolution] )
Sphere(iI_Step,kK_Step) ] ] );
All three will approximate well sphere at higher resolutions. Which one does
it at the low resolution?
On nophead's code I have one thing to say: it works. Why? All I can say is
that nophead's reasoning is most likely wrong, as all rotations will have to
rely on either fsin, fcos or fsincos assembler functions to access the
hardware, and all three are based upon radians, and not degrees. But since a
programmer or his/her compiler may choose between fsin and fcos on the one
hand and fsincos on the other, without the ability to access the
implementation at assembly language level, it's hard to say what happens.
wolf
--
View this message in context: http://forum.openscad.org/Rounding-Errors-tp21821p21934.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
On 07/25/2017 11:52 PM, wolf wrote:
On nophead's code I have one thing to say: it works. Why? All I can say is
that nophead's reasoning is most likely wrong, as all rotations will have to
rely on either fsin, fcos or fsincos assembler functions to access the
hardware, and all three are based upon radians, and not degrees.>
No, that's the whole point, for some special values there is
no need to call the low level function. Those cases just have
known results.
See sin():
https://github.com/openscad/openscad/blob/master/src/func.cc#L256
and cos():
https://github.com/openscad/openscad/blob/master/src/func.cc#L300
The problem is simply that this is not used in all places, like
nophead said. The code works because it explicitly uses those
functions with the special cases, but the internal rotate() does
not. Fixing that would remove the need for the user space function.
ciao,
Torsten.
tp3 wrote
Fixing that would remove the need for the user space function.
So fix it?
Admin - PM me if you need anything, or if I've done something stupid...
Unless specifically shown otherwise above, my contribution is in the Public Domain; to the extent possible under law, I have waived all copyright and related or neighbouring rights to this work. Obviously inclusion of works of previous authors is not included in the above.
The TPP is no simple “trade agreement.” Fight it! http://www.ourfairdeal.org/ time is running out!
View this message in context: http://forum.openscad.org/Rounding-Errors-tp21821p21936.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
On 07/26/2017 12:36 AM, MichaelAtOz wrote:
tp3 wrote
Fixing that would remove the need for the user space function.
So fix it?
That needs both time and motivation. I have none of the first,
comments like that don't help with the latter...
ciao,
Torsten.