Path: utzoo!utgpu!jarvis.csri.toronto.edu!mailrus!cs.utexas.edu!uunet!willett!ForthNet From: ForthNet@willett.UUCP (ForthNet articles from GEnie) Newsgroups: comp.lang.forth Subject: Mathematical routines Message-ID: <251.UUL1.3#5129@willett.UUCP> Date: 11 Jan 90 23:13:38 GMT Organization: Latest Link in ForthNet Chain (Pittsburgh, PA) Lines: 89 Category 3, Topic 10 Message 30 Wed Jan 10, 1990 F.SERGEANT [Frank] at 21:04 CST Jim Matthews on Usenet writes: >I need to take the square root of a number ranging from 0 - >1,596,801,600. I am using MVP-FORTH (79-standard, no floating pt.) with >a few non-standard words that may apply (D*, D/MOD, M*/). I would >appreciate any input. Please reply by sending mail to >floyd@cscwam.umd.edu. Thanx!!! I don't know if this will get back to Usenet or not, but here are some quick thoughts. ( Newton's method for square root ) ( or at least my understanding of it) : ~SQRT ( ud - ud') ( "approximate square root") 32 FOR D2/ 2DUP OR WHILE NEXT THEN ( shift right while number is non-zero ) 2DROP 32 POP ( ie index) - ( power-of-2) 2/ ( power-of-2/2) 1 0 ( now we'll raise it back up to 1/2 the power of original number) ROT ?DUP IF 1- FOR D2* NEXT THEN ( ud') ; ( above gives us our first guess by halving the power of original number) : 2STEP ( ud-sq ud-last-root-guess - ud-sq ud-new-root-guess) 2OVER 2OVER ( sq guess sq guess) 2DUP D* ( sq guess sq guess^2) D- ( sq guess sq-sq') 2OVER D2* ( sq guess sq-sq' derivative) D/ D+ ( sq new-guess) ; ( given the original number and a guess as to its square root, 2STEP finds the difference between the original number and the square of our guess and divides that difference by the derivative dy/dx where y is our original number and x is the square root of y. For y = x^2 ; dy/dx = 2x. This quotient is added to our previous guess of x to get our new guess of x.) : DSQRT ( ud - ud^1/2) 2DUP D0= IF 2DROP 0 0 ( force result to zero if original number=0) ELSE 2DUP ~2ROOT ( original-number 1st-guess ) BEGIN 2DUP 2PUSH ( or >R >R ) ( save guess on return stack ) 2STEP ( improve guess ) 2DUP 2POP D= ( see if new guess is same as old ) UNTIL ( when guess stops changing, quit ) 2SWAP 2OVER ( guess number guess ) 2DUP D* ( guess number guess^2 ) D< IF 1 0 D- THEN ( guess ) ( Make sure that answer isn't too high by one. ) ( This gets around problem of D/ truncating toward zero. ) ( I don't think this would be needed if D/ is floored.) THEN ; I think 2STEP only needs to be done about 5 times at the most. It might be better to leave out the test for the new guess equaling the previous guess and just do 2STEP 5 or 6 times regardless of whether fewer times would be sufficient, e.g. : DSQRT ( ud - ud^1/2) 2DUP D0= IF 2DROP 0 0 ( force result to zero if original number=0) ELSE 2DUP ~2ROOT ( original-number 1st-guess ) 5 FOR 2STEP NEXT ( improve guess 6 times) 2SWAP 2OVER ( guess number guess ) 2DUP D* ( guess number guess^2 ) D< IF 1 0 D- THEN ( guess ) ( Make sure that answer isn't too high by one. ) ( This gets around problem of D/ truncating toward zero. ) ( I don't think this would be needed if D/ is floored.) THEN ; DSQRT gave the correct value for the few numbers I tested, but I'm not sure whether it goes crazy for certain numbers. D/ needs to be signed. You can replace 2PUSH with >R >R and 2POP with R> R>. You can replace 5 FOR ... NEXT with 6 0 DO ... LOOP. You can get the first guess ( ~SQRT) various ways. Actually, the first guess does not need to be very good, if you are willing to repeat 2STEP more times. It might be worth it to hard-code ~SQRT as 7 or 10000 and repeat 2STEP perhaps 20 times ( eg : ~SQRT ( - ud) 7 0 ; ) Good luck. -- Frank ----- This message came from GEnie via willett through a semi-automated process. Report problems to: 'uunet!willett!dwp' or 'willett!dwp@gateway.sei.cmu.edu'