Topic: Error working with Pi

[JasonB]

Hi all,

I've been working on some angle functions and have found what appears to be an error when working with pi. It came up when checking whether an angle is a multiple of 2 pi

For example

Code: [Select]

(/ (* 22.0 pi) (* 2.0 pi))returns 11.0, which is what I expected, however

Code: [Select]

(fix (/ (* 22.0 pi) (* 2.0 pi)))returns 10.

What's going on? I've tried the same in other applications without issue. It is not an isolated incident. Test function below throws up other ones.

Code: [Select]

(defun TEST360 (num / count quotient fixpart remainder)
 (setq count 1)
 (princ "\nCount, Remainder, Quotient, Fixed")
 (repeat num
(princ "\n")
(setq quotient (/ (* (float count) pi) (* 2.0 pi)))
(setq fixpart (fix quotient))
(setq remainder (- quotient fixpart))
(princ (itoa count))
(princ ", ")
(princ (rtos remainder 2 8))
(princ ", ")
(princ (rtos quotient 2 8))
(princ ", ")
(princ (itoa fixpart))
(setq count (1+ count))
 ); end repeat
 (prin1)
); end DEFUN TEST360
I also get the same type issues if I use rem instead of fix. Any ideas on why this is occurring, and how to avoid?

[irneb]

Short answer: It's the fix function which shows up floating point errors.

See this thread: http://www.cadtutor.net/forum/showthread.php?60050-Help-What-wrong-in-Fix-defun

[irneb]

What's going on? I've tried the same in other applications without issue. It is not an isolated incident. Test function below throws up other ones.

Most calculating programs try to account for the error by never converting floats to integers by simply cutting off the fraction part - usually some rounding is used. Or in much more specialist programs they don't even use floating points due to this error, e.g. most accounting proggies use what's known as BCD (Binary Coded Decimal).

Here's a good explanation of the problem (from MS): http://support.microsoft.com/kb/42980

[Lee Mac]

Not the best solution, but I would usually add some tolerance before using fix so as to account for rounding of doubles, e.g.:

Code - Auto/Visual Lisp: [Select]

1. (defun _rem2pi ( x )
2.     (equal 0.0 (rem (+ 1e-8 x) (+ pi pi)) 1e-8)
3. )

Code - Auto/Visual Lisp: [Select]

1. _$ (_rem2pi (* 22.0 pi))
2. T
3. _$ (_rem2pi (* 21.0 pi))
4. nil

[JasonB]

Thanks, forced me to go back and check again. I had been using rem, and new that it was a rounding issue. I had a check for around 0 degs, but hadn't considered a similar check for 180 degs. So I was literally looking in the wrong direction! Had been going round and round on this  .

Code I was working on below for interest. The arccos functions are based on a earlier post http://www.theswamp.org/index.php?topic=575.msg7526#msg7526

Code: [Select]

; ANGLE180
; Adjusts the given angle so that it returns a value between 0 - Pi (0 - 180 degs).
;(angle180 angle); angle has to be in radians
; e.g. (angle180 (getorient '(0.0 0.0) "\nSelect angle: "))
(defun ANGLE180 (angle / checkangle fuzz)
 (setq checkangle (rem angle (+ pi pi))) ; Angle Check to determine whether angle is pi or 2pi
 (setq fuzz 1e-10) ; set accuracy limit when making comparisons
(cond
((equal checkangle 0.0 fuzz) ; if the angle is a multiple of 2Pi (360 deg)
(setq angle 0.0) ; set the angle = 0.0
)
((equal checkangle pi fuzz) ; if the angle is a multiple of Pi (180 deg)
(setq angle pi) ; set the angle = Pi (180 degs)
)
((minusp (sin angle)) ; if the sin of the angle is -ve
(setq angle (+ angle pi)) ; then add Pi (180 degs) to rotate our angle
; Use Cos angle, followed by Arccos to ensure returned angle is between 0 - Pi (0 - 180)
(setq angle (acos (cos angle)))
)
(T  ; Otherwise
; Use Cos angle, followed by Arccos to ensure returned angle is between 0 - Pi (0 - 180)
(setq angle (acos (cos angle)))
)
); end cond
angle ; return angle
); end DEFUN ANGLE180


;INVERSE CIRCULAR FUNCTIONS

;ARCSINE
;Source e 88, pgE7 Engineering Formulaes 7th edtion
;Arcsin(x) = Arctan(x / (1 - x²)^0.5) where x is within -1 & +1
;Note on atan in LISP
;(atan num1 [num2])
;If you supply both num1 and num2 arguments, atan returns the arctangent of num1/num2, in radians.
;You can make use of this feature in the formulae above, in doing so you acheive a more consistent result.
(defun ASIN (x)
  (if (<= -1.0 x 1.0) ; if x is within expected range
         (atan x (sqrt (- 1.0 (* x x)))) ; THEN calculate Arcsine using Arctan
(*error* "\nERROR ACOS(x) not <= -1 x 1.") ; ELSE call error handler
  ); end if
); end DEFUN ASIN

;ARCCOS
;Source e 88, pgE7 Engineering Formulaes 7th edtion
;Arccos(x) = Arctan((1 - x²)^0.5 / x) where x is within -1 & +1
;Note on atan in LISP
;(atan num1 [num2])
;If you supply both num1 and num2 arguments, atan returns the arctangent of num1/num2, in radians.
;You can make use of this feature in the formulae above, in doing so you acheive a more consistent result.
(defun ACOS (x)
  (if (<= -1.0 x 1.0) ; if x is within expected range
         (atan (sqrt (- 1.0 (* x x))) x) ; THEN calculate Arccos using Arctan
(*error* "\nERROR ACOS(x) not <= -1 x 1.") ; ELSE call error handler
  ); end if
); end DEFUN ACOS
used this for testing

Code: [Select]

; RTD
; Convert from Radians to Degrees
(defun RTD (radians)
;180 deg = pi
(* radians (/ 180.0 pi))
)

(defun TEST ( num / angle)
(setq angle 0.0)
(repeat num
(princ "\n")
(princ (rtd angle))
(princ " , ")
(princ (rtd (angle180 angle)))
(princ " , ")
(princ (rem angle (+ pi pi)))
(setq angle (+ angle (/ pi 4.0)))
)
(prin1)
)

[Lee Mac]

How about:

Code - Auto/Visual Lisp: [Select]

1. ;; Ang0-pi  -  Lee Mac
2. ;; Returns the supplied angle in the range 0-pi rads
3.  
4. (defun _ang0-pi ( a )
5.     (if (minusp (setq a (rem a pi)))
6.         (_ang0-pi (+ a pi pi))
7.         (if (or (equal a pi 1e-10) (equal a 0.0 1e-10))
8.             0.0
9.             a
10.         )
11.     )
12. )

[JasonB]

Nice one, wish I had thought of something like that earlier.... Much earlier . Hadn't thought of using recursion. I had gone down the path of cos / arccos, as it will return +ve 0 - Pi no matter that the form the value is supplied in. In a similar way you can use arcsine to return +ve 0 - Pi/2.

Code: [Select]

(defun ANGLE90 (angle / x)
 (asin(abs(sin angle)))
); end DEFUN ANGLE90
getting a function to call itself can be a bit of a mind bender, but man it can be powerful.

Code: [Select]

; LIST2STRING
; converts a list to a string using the provided separator.
; will convert the following data types:
; REAL number -> RTOS using dwg defaults
; INTeger number -> ITOA
; LIST -> calls List2string recursively until all list items are processed.
; For all other data types return the symbol name

; Example of usage
; (list2string '(12 10 ("A" 45 (12.2 54.4)) 13 15 19 "Test" (139 . SYMBOL)) ";")
; returns the following string
; "12;10;A;45;12.2;54.4;13;15;19;Test;139;SYMBOL"

; if it is passed a string instead of a list, it will just pass the string back
; Example
; (list2string "A Text String" ";")
; Returns
; "A Text String"

(defun List2String (listobj sep / listitem itemtype liststr)
(vl-load-com); Make sure the visual lisp extensions are available
(while listobj ;
(if (= 'LIST (type listobj)) ; IF the object provided is a list
(setq
listitem (car listobj) ; THEN retrieve the first item in the list
listobj (cdr listobj) ; and Iterate the list to the remaining items
)
(setq
listitem listobj ; ELSE retrieve the whole object
listobj nil ; set our list to nil
)
) ; end if
(setq itemtype (type listitem)) ; Find the data type for the item retrieved
(cond ; depending on the items data type
(( = 'REAL itemtype)(setq listitem (rtos listitem))) ; Convert the Real number to a string
(( = 'INT itemtype)(setq listitem (itoa listitem))) ; Convert the Integer number to a string
(( = 'LIST itemtype)(setq listitem (list2string listitem sep))) ; Recursively retrieve list items by re-calling list2string
(( = 'STR itemtype)(setq listitem listitem)) ; leave strings as is
(T (setq listitem (vl-symbol-name listitem))) ; for all other data types return the symbol name
) ; end cond

(if liststr ; if our string list already exists
(setq liststr (strcat liststr sep listitem)) ; THEN add our next item onto it, placing the seperator between
(setq liststr listitem); ELSE initialise our string list with the first item
) ; end if
) ; end while
; return processed list
liststr
) ; end list2string

[Lee Mac]

The recursive call was purely to ensure the returned angle was positive since the AutoLISP implementation of the rem function permits negative values to be returned.

Why the 'LIST2STRING' function?

[irneb]

Why the 'LIST2STRING' function?

I think it was as a sample for another use of recursion - note the recursive call for nested lists. A bit off topic but anyhow