Rounding Errors

classic Classic list List threaded Threaded
12 messages Options
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Rounding Errors

wolf
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. 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(i*I_Step,k*K_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, i*S+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, i*S+((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), L*S-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 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

Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Rounding Errors

MathLover
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.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Rounding Errors

cacb
In reply to this post by wolf
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
[hidden email]
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Rounding Errors

doug.moen
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 + step*n + 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, <[hidden email]> 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
[hidden email]
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org


_______________________________________________
OpenSCAD mailing list
[hidden email]
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Rounding Errors

wolf
In reply to this post by wolf
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 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
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Rounding Errors

nophead
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 <[hidden email]> 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
[hidden email]
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org


_______________________________________________
OpenSCAD mailing list
[hidden email]
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Rounding Errors

wolf
In reply to this post by doug.moen
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 on rounding modes, and there is a page on rounding functions. Round-to-even is available from function rint, via setting a rounding mode, and also from function roundeven.
Hardware rounding:
Intel <a href="https://software.intel.com/sites/default/files/managed/39/c5/325462-sdm-vol -1-2abcd-3abcd.pdf">hardware , 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:
 using
AnyShape=(let (Step=360/Resolution) [ for (k=[0 : Step : 180] ) [ for( i=[Step :Step : 360] ) Sphere(i,k) ] ] );  
 using
AnyShape=(let (Step=360/Resolution) [ for (k=drange(0,Step,180)) [ for( i=drange(Step,Step,360) ) Sphere(i,k) ] ] );  
 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(i*I_Step,k*K_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

tp3
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Rounding Errors

tp3
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.

_______________________________________________
OpenSCAD mailing list
[hidden email]
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org
-- Torsten
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Rounding Errors

MichaelAtOz
Administrator
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!
tp3
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Rounding Errors

tp3
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.

_______________________________________________
OpenSCAD mailing list
[hidden email]
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org
-- Torsten
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Rounding Errors

wolf
In reply to this post by tp3
I have always assumed, based on what has been written on this forum in the past, that OpenSCAD functions like rotate(), translate(), scale(), etc are just front ends for multmatrix(). Am I mistaken there and the actual implementation of rotate() is independent of multmatrix()? It is always possible to force special values like sin(90)->sin(pi/2)==1, but is it actually done? From the code references Torsten has given, I at least am not able to ascertain that, and thus I remain sceptical. Kudos to nophead for finding a workaround, but I am looking for more.
I am looking for information that is trustworthy enough that I can build an exception handler on it, an exception handler that lays to rest CGAL assertion violations, once and for all. I have learned enough by now that I can do just that manually - at least if the .stl files are small enough. But I also want to do that under program control, and with large files, so that others can benefit as well. So I need more accuracy.
X86 hardware provides accuracy these days for floating point operations better than 1E-19 decimal or not worse than 2^-65 binary. OpenSCAD, on the other hand, for the script known as issue1258, suffers from accuracy not better than 1E-5 decimal, or 2^-18 binary, which is a 47 binary digits difference. It is this big loss in precision that is at the heart of degenerate triangles, and all that flows from it.

wolf
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Rounding Errors

kintel
Administrator
Some short answers to requested info:

* OpenSCAD’s round: https://github.com/openscad/openscad/blob/master/src/func.cc#L387
* OpenSCAD’s linear transformations: https://github.com/openscad/openscad/blob/master/src/transform.cc#L86

OpenSCAD is source code that can be compiled on different C++ compilers. We use the default compiler for Linux distros and usually the most recent platform compiler for OS X. For Windows, we use whatever gcc version http://mxe.cc provides.

The source code is always going to be the authoritative information about how anything is implemented. If you cannot find it, ask here and someone can usually point you to the correct location.

 -Marius


_______________________________________________
OpenSCAD mailing list
[hidden email]
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org
Loading...