QB64.org Forum

Active Forums => Programs => Topic started by: jack on June 06, 2021, 09:11:41 am

Title: reciprocal without division
Post by: jack on June 06, 2021, 09:11:41 am

just for fun, a reciprocal function without division using the fast_pow trick http://www.qb64.org/wiki/C_Libraries

Code: QB64: [Select]

1. Declare CustomType Library
2.     Sub memcpy (ByVal dest As _Offset, Byval source As _Offset, Byval bytes As Long)
3. End Declare
4.  
5.  
6. Print recip(3)
7.  
8. Function recip# (x#)
9.     Dim As Double a
10.     Dim As Long tmp
11.     a = x#
12.     memcpy _Offset(tmp), _Offset(a) + 4, 4
13.     tmp = tmp - 1072632447
14.     tmp = 1072632447 - tmp
15.     memcpy _Offset(a) + 4, _Offset(tmp), 4
16.     ' refine the initial approximation using the Newton-Raphson method
17.     a = a * (2 - x# * a)
18.     a = a * (2 - x# * a)
19.     a = a * (2 - x# * a)
20.     a = a * (2 - x# * a)
21.     recip = a
22. End Function
23.  

Quote

.3333333333333333

Title: Re: reciprocal without division
Post by: david_uwi on June 06, 2021, 10:22:34 am

I'm an advocate of keeping it simple. These fancy optimization routines are all good using paper and pencil, but a computer is fast and even faster if the calculations required on each iteration are simple. So I present bisection for reciprocals...

Code: QB64: [Select]

1. PRINT "input a number"
2. INPUT x
3. y1 = x
4. y2 = 0
5. sq0 = 9999
6. DO
7.     sq1 = (y1 + y2) / 2
8.     q1 = sq1 * x
9.     IF q1 > 1! THEN y1 = sq1
10.     IF q1 < 1! THEN y2 = sq1
11.     IF ABS(sq1 - sq0) < .00000001 THEN EXIT DO
12.     sq0 = sq1
13. LOOP
14. PRINT sq1
15.  

Title: Re: reciprocal without division
Post by: jack on June 06, 2021, 11:04:21 am

@david_uwi it doesn't count if you use division, see the thread's title
you could change / 2 by * 0.5, but I am willing to bet that my method is both much faster an more accurate
try x = 300000000000000#

Text Only | Text with Attachments