Jump to content

Calculations with BIG numbers

Costinbos77

By Costinbos77
in AutoLISP, Visual LISP & DCL

 Share

Recommended Posts

Costinbos77 Rising Star

Costinbos77

Posted
Posted (edited)

Hi ,

if p(i) are as (491677.5042  371092.5820) , how will be possible to do the next calculations : 

(setq selP nil  selP (ssGet '((0 . "Point"))) ) ; Select only 3 points !
(if (and (setq e0 (ssName selP 0)) (setq e1 (ssName selP 1)) (setq e2 (ssName selP 2)) )
 (progn
  (setq p1 (trans (cdr (assoc 10 (entGet e0))) 0 1)  p2 (trans (cdr (assoc 10 (entGet e1))) 0 1)  p3 (trans (cdr (assoc 10 (entGet e2))) 0 1) )
;;;  (setq p1 (cdr (assoc 10 (entGet e0)))  p2 (cdr (assoc 10 (entGet e1)))  p3 (cdr (assoc 10 (entGet e2))) )
  (if (and p1 p2 p3) ; Locale
   (progn
    (setq x1 (car p1)  y1 (cadr p1)   x2 (car p2)  y2 (cadr p2)   x3 (car p3)  y3 (cadr p3)
          a (* (+ (* (- y2 y3) x1) (* (- x3 x2) y1) (* x2 y3) (* x3 y2 -1.) ) 2.)
          x (/ (+ (* (+ (* x1 x1) (* y1 y1)) (- y2 y3))  (* (+ (* x2 x2) (* y2 y2)) (- y3 y1))  (* (+ (* x3 x3) (* y3 y3)) (- y1 y2)) ) a)
          y (/ (+ (* (+ (* x1 x1) (* y1 y1)) (- x3 x2))  (* (+ (* x2 x2) (* y2 y2)) (- x1 x3))  (* (+ (* x3 x3) (* y3 y3)) (- x2 x1)) ) a)
          dX (- x x1)  dY (- y y1)  r (sqrt (+ (* dX dX) (* dY dY)))
    ) ; s
    (entMake (list '(0 . "Circle") (cons 8 "Lay1") '(6 . "ByLayer") '(62 . 256) (cons 10 (trans (list x y 0) 1 0) ) (cons 40 r) )  ) ; end of e
   ) ; prog then
   (prinC "\n   NO  3  VALID  Points  ;")
  ) ; if
)) ; if and

With numbers as (1001.123 5001.123) or (61001.123 51001.123) or (100,005.123 100,006.123) it is working fine .

But , with numbers bigger than 100,000 is giving wrong answers of calculations .

 

Regards ,

Costin

Edited by Costinbos77
marko_ribar Veteran

marko_ribar

Posted
Posted

What is solution - CIRCLE in relation to triangle? Picture would be fine... BTW. CAD can't calculate correct with big numbers - that's known issue... You could though try to truncate numbers by scaling coordinates of points through reference point and then finally rescale CIRCLE as result also through that same reference point...

Costinbos77 Rising Star

Costinbos77

Posted
  • Author
Posted

Hi marko_ribar ,

 

Quote

...

What is solution - CIRCLE in relation to triangle? Picture would be fine... BTW

...

 

It is about 3 points on a circle .

 

The real problem is : (* x1 x1)  , (* y1 y1) ... because are generating big numbers .

 

I tried to divide the number by  1000 , is not working , because exactly the last decimals are making the difference .

Will work something if from a numer as 491677.5042 , I will use only the small part , as 677.5042 and at the end , to add 491000.0 .

 

Regards,

Costin

marko_ribar Veteran

marko_ribar

Posted
Posted (edited)

Just a side thought... If solution is circumcircle of triangle, why do you need to use math formula... I'd prefer and suggest that you construct resulting circle... And also good to mention, you have built in CIRCLE command with 3P (3 point) option, so you can directly draw circle through command usage...

Edited by marko_ribar
Costinbos77 Rising Star

Costinbos77

Posted
  • Author
Posted (edited)

 

Quote

...

you have built in CIRCLE command with 3P (3 point) option, so you can directly draw circle through command usage...

...

Good point marko_ribar ! But I think will be slow for 500 circles or more .

 

This is the result :

1839464647_3.13DevC2.jpg.2f475773c84af19cd0447e54e025c23c.jpg

 

I need the center of the circle defined by 3 points and compared with the center of the designed one.

 

So , I reduce the values if are bigger then 100,000 : 

(setq selP nil  selP (ssGet '((0 . "Point"))) ) ; Select only 3 points !
(if (and (setq e0 (ssName selP 0)) (setq e1 (ssName selP 1)) (setq e2 (ssName selP 2)) )
 (progn
  (setq p1 (trans (cdr (assoc 10 (entGet e0))) 0 1)  p2 (trans (cdr (assoc 10 (entGet e1))) 0 1)  p3 (trans (cdr (assoc 10 (entGet e2))) 0 1) )
;;;  (setq p1 (cdr (assoc 10 (entGet e0)))  p2 (cdr (assoc 10 (entGet e1)))  p3 (cdr (assoc 10 (entGet e2))) )
  (if (and p1 p2 p3) ; Locale
   (progn
;;;    (setq x1 (/ (car p1) 1000.)  y1 (/ (cadr p1) 1000.)   x2 (/ (car p2) 1000.)  y2 (/ (cadr p2) 1000.)   x3 (/ (car p3) 1000.)  y3 (/ (cadr p3) 1000.)
    (setq x1 (car p1)  y1 (cadr p1)   x2 (car p2)  y2 (cadr p2)   x3 (car p3)  y3 (cadr p3) )
    
(if (or (> x1 100000) (> x2 100000) (> x3 100000) )
 (setq x0 (* (fix (/ (min x1 x2 x3) 1000.)) 1000.)  x1 (- x1 x0)  x2 (- x2 x0)  x3 (- x3 x0) )
 (setq x0 0)
) ; if

(if (or (> y1 100000) (> y2 100000) (> y3 100000) )
 (setq y0 (* (fix (/ (min y1 y2 y3) 1000.)) 1000.)  y1 (- y1 y0)  y2 (- y2 y0)  y3 (- y3 y0)  )
 (setq y0 0)
) ; if
    
    (setq a (* (+ (* (- y2 y3) x1) (* (- x3 x2) y1) (* x2 y3) (* x3 y2 -1.) ) 2.)
          x (+ (/ (+ (* (+ (* x1 x1) (* y1 y1)) (- y2 y3))  (* (+ (* x2 x2) (* y2 y2)) (- y3 y1))  (* (+ (* x3 x3) (* y3 y3)) (- y1 y2)) ) a) x0)
          y (+ (/ (+ (* (+ (* x1 x1) (* y1 y1)) (- x3 x2))  (* (+ (* x2 x2) (* y2 y2)) (- x1 x3))  (* (+ (* x3 x3) (* y3 y3)) (- x2 x1)) ) a) y0)
;;;          x (* x 1000.)  y (* y 1000.)
          dX (- (+ x0 x1) x)  dY (- (+ y0 y1) y)  r (sqrt (+ (* dX dX) (* dY dY)))
    ) ; s
    (entMake (list '(0 . "Circle") (cons 8 "Lay1") '(6 . "ByLayer") '(62 . 256) (cons 10 (trans (list x y 0) 1 0) ) (cons 40 r) )  ) ; end of e
   ) ; prog then
   (prinC "\n   NO  3  VALID  Points  ;")
  ) ; if
)) ; if and

At the moment , it is working like that . will see next .

 

 

Regards,

 

Costin

Edited by Costinbos77
dlanorh Mentor

dlanorh

Posted
Posted (edited)

Have you tried doing your calculations as logs?  (* a b) = (exp (+ (log a) (log b))). The natural logs are much smaller numbers i.e. (log 491677.5042) = 13.1056 (to 4 decimal places)

 

to : multiply two numbers add their logs as above

      divide two numbers subtract their logs i.e. (/ 4.0 2.0) = (exp (- (log 4.0) (log 2.0)))

      sqrt = (exp (/ (log num) 2.0))

 

and then use exp function (antilog) to convert back to a number

Edited by dlanorh
clarify
lrm Experienced

lrm

Posted
Posted

dlanorh, I do not think using logs gains you any precision. As you know real numbers in AutoCAD are stored as a mantissa and an exponent. The 56 bit mantissa yields about 15 significant figures of precision. The exponent uses the remaining 8 bits of the 64 bit value.   A log in effect represents the mantissa and exponent as one number. It's not clear if additional bits can be used beyond 64 to store the log. 

It's true that with real numbers if, for example, you  multiply a 10 digit number by another 10 digit number you can get a result with 20 digits of precision.  The result will be rounded (truncated?) to 15 digits within AutoCAD.   With logs you would add the two number resulting in a rounding of only the least significant digit. The question becomes what is the precision of the antilog of the number?

I think the solution is to use a geometric solution that avoids multiplication as Lee Mac's code does by finding the intersection of the two perpendicular bisectors of the 3 points.  This leads to the question of the precision of the inters function and its affect on the result.   I think this could become significant if the 3 points are almost collinear as the perpendicular bisectors will be almost parallel.

Costinbos77 Rising Star

Costinbos77

Posted
  • Author
Posted (edited)

1. Interesting idea dlanorh . Thanks !

 

I was trying like this : 

(setq x (exp (- (log 
 (+ (exp (+ (log (+ (exp (+ (log x1) (log x1))) (exp (+ (log y1) (log y1)))) ) (log (- y2 y3))))
    (exp (+ (log (+ (exp (+ (log x2) (log x2))) (exp (+ (log y2) (log y2)))) ) (log (- y3 y1))))
    (exp (+ (log (+ (exp (+ (log x3) (log x3))) (exp (+ (log y3) (log y3)))) ) (log (- y1 y2))))
  ) ) (log a)) ) )

but, sometimes when : a< 0 , (- y2 y3) < 0 ... , log will not work with them .

 

2. Thanks BIGAL ! I will do some research on it .

 

3 .Thanks lrm , for your detailed answer !

 

Regards ,

Costin

Edited by Costinbos77
Stefan BMR Proficient

Stefan BMR

Posted
Posted

Salut Costin

 

Looking at your formula, something seems off. The "a" variable almost cries for a symmetrical formula... Just saying.

Maybe something like: 

(* (+ (* (- y2 y3) x1) (* (- y3 y1) x2) (* (- y1 y2) x3)) 2.)

Works just fine for points' coordinates around 100,000.00 and a circle about 1.00 radius.

Costinbos77 Rising Star

Costinbos77

Posted
  • Author
Posted (edited)

Hi Stefan BMR ,

 

Formulas , are fine and are working with any coordinates (see the bottom coordinates from video) .

It is not a perfect solution , but , for my needs it is perfect !

It is important to work in the Origin (0 0 0) for UCS and somewhere between 300,000.0 - 600,000.0 ( national coordinates system ).

The radius of circles will not be bigger than 1.5 m .

 

My solution will have a problem when some values of coordinates will be below 100,000.0 and others will be above of 100,000.0 , which is not my case .

This is why I changed the test from (> x1 100,000.0) to  (> x1 100,005.0) . I don't know exactly the limit , but I know the values as 100,006.0 , are working fine .

 

An UCS around the working area , will drop down the values of coordinates and will be no problems with big numbers in formulas .

 

Regards ,

Costin

 

Edited by Costinbos77
Stefan BMR Proficient

Stefan BMR

Posted
Posted

So, it's not a problem with big numbers. Using the correct formula you can get the circle you need.

Look at the sample below. That tiny difference in "a" value is really important for big numbers .

For small numbers, it's not, but I don't get it, if you know the formula is wrong, why not change it?

  

X1   299999.929857717
Y1   599999.8083234681
X2   299999.3306365923
Y2   600000.0352543053
X3   299998.7266854801
Y3   599999.749548545

(setq a (* (+ (* (- y2 y3) x1) (* (- x3 x2) y1) (* x2 y3) (* x3 y2 -1.) ) 2.))
(rtos a 2 16) -> 0.6165771484375
(/ x1 a) -> 486557.0

(setq a (* (+ (* (- y2 y3) x1) (* (- y3 y1) x2) (* (- y1 y2) x3)) 2.))
(rtos a 2 16) -> 0.6165121169469785
(/ x1 a) -> 486608.0

 

Costinbos77 Rising Star

Costinbos77

Posted
  • Author
Posted

Stefan BMR , 

the formulas are right !

 

1945122727_3.13DevC3.thumb.jpg.9d378d879b38232ad066895fce7ae321.jpg

 

Use the above code lines , and you will have the same results .

 

Regards ,

Costin

Stefan BMR Proficient

Stefan BMR

Posted
Posted

Well, you are right, both formulas are the same, the terms are just rearranged!

In this case, I guess the magnitude of the partial terms makes the difference.

In "a" expression, (* (- y2 y3) x1) is in the same range as x1 and (* x2 y3) is 10^6 bigger, which leads to a tiny difference.

With your formula in WCS, the circle is always off. With the other formula, I had no errors in my tests.

 

Costinbos77 Rising Star

Costinbos77

Posted
  • Author
Posted (edited)

Hi Stefan BMR ,

 

45 minutes ago, Stefan BMR said:

...

In "a" expression, (* (- y2 y3) x1) is in the same range as x1 and (* x2 y3) is 10^6 bigger, which leads to a tiny difference.

...

 

 

Good point , thanks !

 

Regards,

Costin

Edited by Costinbos77

Join the conversation

You can post now and register later. If you have an account, sign in now to post with your account.
Note: Your post will require moderator approval before it will be visible.

Guest
Unfortunately, your content contains terms that we do not allow. Please edit your content to remove the highlighted words below.
Reply to this topic...

×   Pasted as rich text.   Restore formatting

  Only 75 emoji are allowed.

×   Your link has been automatically embedded.   Display as a link instead

×   Your previous content has been restored.   Clear editor

×   You cannot paste images directly. Upload or insert images from URL.

×
 Share
Go to topic listing
×
×
  • Create New...