Power function

[bplus]

For the math fans, I stumbled upon a way to do power function that might be interesting:

Code: QB64: [Select]

1. _TITLE "Power Function by bplus"
2. 'QB64 X 64 version 1.2 20180228/86  from git b301f92
3.  
4. ' started 2018-07-23  Naalaa has no power function (or operator), so I wrote a power function for it.
5. ' ''Power pack the short version.txt
6. ' ''written for Naalaa 6 by bplus posted 2018-07-23
7. ' ''extracted from large test of fractions, random number functions, string... in Power test.txt
8.  
9. ' OMG the crazy thing worked! It produced decent estimates of roots given the limitations of precision...
10. ' Now I want to see how well it works with far greater precision available. So here we are, looking to see
11. ' how this function compares to the regualar ^ operator in QB64.
12.  
13. 'from Naalaa comments:
14. ' The main purpose of this code: to demo a power function for real numbers,
15.  
16. ' I had an idea for how real numbers to the power of real numbers might be done ie x ^ y = ?
17. ' This means that not only can you take SQR of a number, you can get cube or cube root, quartic, 5th 6th... roots and any multiple
18.  
19. ' It came from this simple idea
20. ' 2 ^ 3.5 = 2 ^ 3 * 2 ^ .5 = 8 * Sqr(2)
21. ' 3 ^ 3.5 = 3 ^ 3 * 3 ^ .5 = 27 * Sqr(3)
22.  
23. ' so 2 ^ 3.25 = 2 ^ 3 * 2 ^ .25
24. ' what is 2 ^ .25 ?  It is sqr(sqr(2)) !
25.  
26. ' likewise 2 ^ 3.125 = 2 ^ 3 * 2 ^ 1/8
27. ' what is 2 ^ 1/8 ? It is sqr(sqr(sqr(2))) !
28.  
29. ' any decimal can be written as a sum of fraction powers of 2 ie 1/2^n, as any integer can be written in powers of 2.
30. ' in binary expansions
31. ' 1/2 = .1       or .5 base 10
32. ' 1/4 = .01      or .25 base 10
33. ' 1/8 = .001     or .125 base 10
34. ' 1/16 = .0001   or .0625 base 10
35.  
36. ' So with binary expansion of decimal, we can SQR and multiply our way to an estimate
37. ' of any real number to the power of another real number using binary expansion of
38. ' the decimal parts as long as we stay in Naalaa integer limits and are mindful of precision.
39.  
40. CONST wW = 800
41. CONST wH = 600
42. SCREEN _NEWIMAGE(wW, wH, 32)
43. _SCREENMOVE 360, 60
44. _DEFINE A-Z AS _FLOAT
45.  
46. DO
47.     PRINT "Testing the power(x, pow) function:"
48.     INPUT "(nothing) quits, Please enter a real number to raise to some power. x = "; x
49.     IF x = 0 THEN EXIT DO
50.     INPUT "(nothing) quits, Please enter a real number for the power. pow = ", pw
51.     IF pw = 0 THEN EXIT DO
52.     result = power(x, pw)
53.     PRINT result; " is what we estimate for"; x; " raised to power of"; pw
54.     PRINT x ^ pw; " is what the ^ operator gives us."
55.     PRINT
56. LOOP
57. PRINT
58. PRINT "This is matching the ^ operator very well! This code is clear proof of concept!"
59. PRINT " OMG, it worked!!!"
60. SLEEP
61.  
62.  
63. ' A power function for real numbers (small ones but still!)
64. ' x to the power of pow
65. FUNCTION power## (x AS _FLOAT, pow AS _FLOAT)
66.     'this sub needs 2 other subs
67.     'bExpand20$
68.     'split though for this split is overkill
69.     DIM build AS _FLOAT
70.     REDIM s(0) AS STRING
71.     r$ = "0" + STR$(pow) 'in case pow starts with decimal
72.     Split s$(), r$, "."
73.     integer$ = s$(0)
74.     build = 1.0
75.     IF integer$ <> "0" THEN
76.         p = VAL(integer$)
77.         FOR i = 1 TO p
78.             build = build * x
79.         NEXT
80.     END IF
81.     IF UBOUND(s$) = 0 THEN power = build: EXIT SUB
82.     'that takes care of integer part,
83.     'now for the fraction part convert decimal to fraction
84.     n$ = s$(1)
85.     ld = LEN(n$)
86.     WHILE RIGHT$(n$, 1) = "0"
87.         n$ = LEFT$(n$, ld - 1)
88.         ld = LEN(n$)
89.     WEND
90.     denom& = 10
91.     FOR i = 2 TO ld
92.         denom& = denom& * 10
93.     NEXT
94.     numer& = VAL(n$)
95.     'OK for bExpand20$ don't have to simplify and that saves us having to extract n and d again from n/d
96.     bs$ = bExpand60$(numer&, denom&)
97.     'at moment we haven't taken any sqr of x
98.     runningXSQR = x
99.     'run through all the 0's and 1's in the bianry expansion of the fraction part of the power float
100.     FOR i = 1 TO LEN(bs$)
101.         'this is the matching sqr of the sqr of the sqr... of x
102.         runningXSQR = SQR(runningXSQR)
103.         'for every 1 in the expansion, multiple our build with the running sqr of ... sqr of x
104.         IF MID$(bs$, i, 1) = "1" THEN build = build * runningXSQR
105.     NEXT
106.     'our build should be a estimate or x to power of pow
107.     power = build
108. END FUNCTION
109.  
110. 'write a series of 1s and 0s that represent the decimal fraction n/d in binary 60 places long
111. FUNCTION bExpand60$ (nOver&, d&)
112.     DIM b AS _FLOAT, r AS _FLOAT
113.     ' b for base
114.     b = 0.5
115.     ' r for remainder
116.     r = nOver& / d&
117.     ' s for string$ 0's and 1's that we will build and return for function value
118.     s$ = ""
119.     ' f for flag to stop
120.     f% = 0
121.     ' c for count to track how far we are, don't want to go past 20
122.     c% = 0
123.     WHILE f% = 0
124.         IF r < b THEN
125.             s$ = s$ + "0"
126.         ELSE
127.             s$ = s$ + "1"
128.             IF r > b THEN
129.                 r = r - b
130.             ELSE
131.                 f% = 1
132.             END IF
133.         END IF
134.         c% = c% + 1
135.         IF c% >= 60 THEN f% = 1
136.         b = b * 0.5
137.     WEND
138.     bExpand60$ = s$
139. END FUNCTION
140.  
141. SUB Split (arr() AS STRING, mystr AS STRING, delim AS STRING)
142.     DIM c AS STRING * 1, curpos AS LONG, arrpos AS LONG, le AS LONG
143.     curpos = 1: arrpos = 0
144.     le = LEN(mystr)
145.     DO UNTIL (curpos > le)
146.         c = MID$(mystr, curpos, 1) ' Get 1 ctr
147.         IF (c <> delim) THEN
148.             arr(arrpos) = arr(arrpos) + c
149.         ELSE
150.             arrpos = arrpos + 1
151.             REDIM _PRESERVE arr(arrpos + 1) AS STRING 'AS STRING required or it assumes the default type which is SINGLE
152.         END IF
153.         curpos = curpos + 1
154.     LOOP
155.     REDIM _PRESERVE arr(arrpos) AS STRING
156. END SUB
157.  

[bplus]

Update Naalaa does have a pow#(x#, y#) function, just not documented (nor keyword highlighted).

[bplus]

I have replaced split with two new (for my toolbox) very handy string functions, leftOf$ and rightOf$.

Like [banned user], constantly working on ever better code.

Code: QB64: [Select]

1. _TITLE "Power Function 2 by bplus"
2. 'QB64 X 64 version 1.2 20180228/86  from git b301f92
3.  
4. ' started 2018-07-23  Naalaa has no power function (or operator), so I wrote a power function for it.
5. ' ''Power pack the short version.txt
6. ' ''written for Naalaa 6 by bplus posted 2018-07-23
7. ' ''extracted from large test of fractions, random number functions, string... in Power test.txt
8. ' 2018-07-24 Power Function 2 split is replaced with two much smaller and very handy string functions.
9.  
10. ' OMG the crazy thing worked! It produced decent estimates of roots given the limitations of precision...
11. ' Now I want to see how well it works with far greater precision available. So here we are, looking to see
12. ' how this function compares to the regualar ^ operator in QB64.
13.  
14. 'from Naalaa comments:
15. ' The main purpose of this code: to demo a power function for real numbers,
16.  
17. ' I had an idea for how real numbers to the power of real numbers might be done ie x ^ y = ?
18. ' This means that not only can you take SQR of a number, you can get cube or cube root, quartic, 5th 6th... roots and any multiple
19.  
20. ' It came from this simple idea
21. ' 2 ^ 3.5 = 2 ^ 3 * 2 ^ .5 = 8 * Sqr(2)
22. ' 3 ^ 3.5 = 3 ^ 3 * 3 ^ .5 = 27 * Sqr(3)
23.  
24. ' so 2 ^ 3.25 = 2 ^ 3 * 2 ^ .25
25. ' what is 2 ^ .25 ?  It is sqr(sqr(2)) !
26.  
27. ' likewise 2 ^ 3.125 = 2 ^ 3 * 2 ^ 1/8
28. ' what is 2 ^ 1/8 ? It is sqr(sqr(sqr(2))) !
29.  
30. ' any decimal can be written as a sum of fraction powers of 2 ie 1/2^n, as any integer can be written in powers of 2.
31. ' in binary expansions
32. ' 1/2 = .1       or .5 base 10
33. ' 1/4 = .01      or .25 base 10
34. ' 1/8 = .001     or .125 base 10
35. ' 1/16 = .0001   or .0625 base 10
36.  
37. ' So with binary expansion of decimal, we can SQR and multiply our way to an estimate
38. ' of any real number to the power of another real number using binary expansion of
39. ' the decimal parts as long as we stay in Naalaa integer limits and are mindful of precision.
40.  
41. CONST wW = 800
42. CONST wH = 600
43. SCREEN _NEWIMAGE(wW, wH, 32)
44. _SCREENMOVE 360, 60
45. _DEFINE A-Z AS _FLOAT
46.  
47. DO
48.     PRINT "Testing the power(x, pow) function:"
49.     INPUT "(nothing) quits, Please enter a real number to raise to some power. x = "; x
50.     IF x = 0 THEN EXIT DO
51.     INPUT "(nothing) quits, Please enter a real number for the power. pow = ", pw
52.     IF pw = 0 THEN EXIT DO
53.     result = power(x, pw)
54.     PRINT result; " is what we estimate for"; x; " raised to power of"; pw
55.     PRINT x ^ pw; " is what the ^ operator gives us."
56.     PRINT
57. LOOP
58. PRINT
59. PRINT "This is matching the ^ operator very well! This code is clear proof of concept!"
60. PRINT " OMG, it worked!!!"
61. SLEEP
62.  
63.  
64. ' A power function for real numbers (small ones but still!)
65. ' x to the power of pow
66. FUNCTION power## (x AS _FLOAT, pow AS _FLOAT)
67.     'this sub needs: bExpand60$, leftOf$, rightOf$
68.  
69.     DIM build AS _FLOAT
70.     r$ = "0" + STR$(pow) 'in case pow starts with decimal
71.     integer$ = leftOf$(r$, ".")
72.     build = 1.0
73.     IF integer$ <> "0" THEN
74.         p = VAL(integer$)
75.         FOR i = 1 TO p
76.             build = build * x
77.         NEXT
78.     END IF
79.     'that takes care of integer part
80.  
81.     n$ = rightOf$(r$, ".")
82.     IF n$ = "" THEN power = build: EXIT SUB
83.  
84.     'remove 0's to right of main digits
85.     ld = LEN(n$)
86.     WHILE RIGHT$(n$, 1) = "0"
87.         n$ = LEFT$(n$, ld - 1)
88.         ld = LEN(n$)
89.     WEND
90.  
91.     'note: we are pretending that the ^ operator is not available, so this is hand made integer power
92.     denom& = 10
93.     FOR i = 2 TO ld
94.         denom& = denom& * 10
95.     NEXT
96.  
97.     'OK for bExpand60$ don't have to simplify fraction and that saves us having to extract n and d again from n/d
98.     bs$ = bExpand60$(VAL(n$), denom&)
99.  
100.     'at moment we haven't taken any sqr of x
101.     runningXSQR = x
102.  
103.     'run through all the 0's and 1's in the bianry expansion, bs$, the fraction part of the power float
104.     FOR i = 1 TO LEN(bs$)
105.         'this is the matching sqr of the sqr of the sqr... of x
106.         runningXSQR = SQR(runningXSQR)
107.         'for every 1 in the expansion, multiple our build with the running sqr of ... sqr of x
108.         IF MID$(bs$, i, 1) = "1" THEN build = build * runningXSQR
109.     NEXT
110.  
111.     'our build should now be an estimate or x to power of pow
112.     power = build
113. END FUNCTION
114.  
115. 'write a series of 1s and 0s that represent the decimal fraction n/d in binary 60 places long
116. FUNCTION bExpand60$ (nOver&, d&)
117.     DIM b AS _FLOAT, r AS _FLOAT
118.     ' b for base
119.     b = 0.5
120.     ' r for remainder
121.     r = nOver& / d&
122.     ' s for string$ 0's and 1's that we will build and return for function value
123.     s$ = ""
124.     ' f for flag to stop
125.     f% = 0
126.     ' c for count to track how far we are, don't want to go past 20
127.     c% = 0
128.     WHILE f% = 0
129.         IF r < b THEN
130.             s$ = s$ + "0"
131.         ELSE
132.             s$ = s$ + "1"
133.             IF r > b THEN
134.                 r = r - b
135.             ELSE
136.                 f% = 1
137.             END IF
138.         END IF
139.         c% = c% + 1
140.         IF c% >= 60 THEN f% = 1
141.         b = b * 0.5
142.     WEND
143.     bExpand60$ = s$
144. END FUNCTION
145.  
146. FUNCTION leftOf$ (source$, of$)
147.     posOf = INSTR(source$, of$)
148.     IF posOf > 0 THEN leftOf$ = MID$(source$, 1, posOf - 1)
149. END FUNCTION
150.  
151. FUNCTION rightOf$ (source$, of$)
152.     posOf = INSTR(source$, of$)
153.     IF posOf > 0 THEN rightOf$ = MID$(source$, posOf + LEN(of$))
154. END FUNCTION
155.  

I have plans for these 2 stringy's, possible for generic properties  keyword = value, or more decimal work or work with fractions with notations _ dividing integer part with fraction part and / dividing numerator and denominator. Let's see what else...?

[Ashish]

Hi bplus! I like it.
However, I can do same thing in just 5 lines of function.

Code: QB64: [Select]

1. FUNCTION power2## (x AS _FLOAT, pow AS _FLOAT)
2.     DIM y AS _FLOAT
3.     y = _CEIL(pow)
4.     IF y = pow THEN power2## = x ^ y ELSE power2## = (x ^ (y - 1)) * (x ^ (pow - (y - 1)))
5. END FUNCTION
6.  

EDIT:
Oh! I noticed, You didn't use even a single '^' in your function. That's the reason you are using string functions.
Amazing!

[bplus]

Yes, Ashish but are you missing the point?

Can you do it without using ^?

That is the challenge!

I mean someone had to figure out how to code the ^ operator.

I just did yesterday and am pretty impressed with myself! ;D

[bplus]

Hi bplus! I like it.
However, I can do same thing in just 5 lines of function.

Code: QB64: [Select]

1. FUNCTION power2## (x AS _FLOAT, pow AS _FLOAT)
2.     DIM y AS _FLOAT
3.     y = _CEIL(pow)
4.     IF y = pow THEN power2## = x ^ y ELSE power2## = (x ^ (y - 1)) * (x ^ (pow - (y - 1)))
5. END FUNCTION
6.  

EDIT:
Oh! I noticed, You didn't use even a single '^' in your function. That's the reason you are using string functions.
Amazing!

OK now you get it! (maybe starting to anyway) ;)

Hey why did you use _CEIL? Wouldn't that turn it to integer and miss out doing roots?

[bplus]

Math lesson:

If x ^ 7 = 1000 what does x = ?

x = 1000 ^ (1/7)

if x ^ 3 = 343 what does x = ?

x = 343 ^ (1/3) the cube root of 343 which is 7 if my multiplying was correct.
x = 343 ^ (.33333...) same thing

Did you know?  SQR(n) = n ^ .5

EDIT: my multiplication was incorrect, sorry 7 ^ 3 <> 483

[bplus]

Code: QB64: [Select]

1. _TITLE "Tower of power and Circle of power by bplus"
2. 'started 2018-07-24 bplus
3. 'using QB64 X 64 version 1.2 20180228/86  from git b301f92
4. CONST wW = 800
5. CONST wH = 600
6. DIM SHARED pi
7. pi = _PI
8. SCREEN _NEWIMAGE(wW, wH, 32)
9. _SCREENMOVE 360, 60
10.  
11. ' if x ^ 7 = wH, what does x = ?
12. x = wH ^ (1 / 7)
13. 'so if I multiply x times itself 7 times I will go the height of this screen
14. cx = wW / 2
15. FOR power = 50 TO 1 STEP -1
16.     CLS
17.     LOCATE 1, 1: PRINT power
18.     x = (wH - 20) ^ (1 / power)
19.     FOR i = 0 TO power - 1
20.         LINE (cx - .5 * x ^ i, 10 + x ^ i)-(cx + .5 * x ^ i, 10 + x ^ (i + 1)), , B
21.     NEXT
22.     _DISPLAY
23.     _LIMIT 10
24. NEXT
25. CLS
26. _DELAY .03
27. cy = wH / 2
28. 'dividing the radius a number of times...
29. FOR power = 50 TO 1 STEP -1
30.     CLS
31.     LOCATE 1, 1: PRINT power
32.     x = (wH - 20) ^ (1 / power)
33.     FOR i = 1 TO power
34.         CIRCLE (cx, cy), .5 * (wH - 20 - x ^ i)
35.     NEXT
36.     _DISPLAY
37.     _LIMIT 10
38. NEXT
39. CLS
40.  
41.  

[johnno56]

Well... You just destroyed my right retina and the left one took a heavy beating... I hope you're happy!! Nah. Kidding. Great job!

J

[bplus]

Ha! burning retinas was not my intension, this time. ;)

I was hoping to illustrate a possible use of power for drawing. Maybe I should keep looking... :)

I did have an idea for planting trees with the rows getting ever closer and the trees ever smaller as they approached the horizon line.

Perhaps an army of space aliens instead of trees would be more dramatic, or rats maybe!

[STxAxTIC]

I like it bplus. Unfortunately I'm not near a QB64-enabled machine for a little while and thus didn't give it an explicit test run - so forgive my ignorance if you've answered this/these:

Can this method extend to very high precision numbers, like 2.7189876543212345678998765432123456789  ^ 3.1498765432345678998765432123456789 ?

And

By a similar-sounding but entirely different token, can it handle very large powers, like 999 ^ 999?

[bplus]

Hi STxAxTIC,

The QB64 code would fail as soon as STR$(the power) starts including letters D or E or whatever for scientific notation.

Code: QB64: [Select]

1. r$ = "0" + STR$(pow) 'in case pow starts with decimal

That would be the first thing (in code sequence) to fail, the next would be _FLOAT's precision.

I am sure the method would be possible for string math with arbitrary long decimal or binary expressions and sqr but without practical need I don't know if want to do all that work.

[STxAxTIC]

Heya bplus,

Thanks for the swift answer. Your post inspired me to waste 20 minutes fiddling with trigonometric expansions of real numbers. It happened to be a total waste of paper and ink, but hardly a waste of time. Say, want any of your polished work(s) to be immortalized in Samples? It wouldn't bother Odin a bit.

[bplus]

Hi STxAxTIC,

Sure Odin may use what please's them.

Can you tell I am particularly pleased with this little gem despite it's limits?

I searched for a Taylor expansion but could find no direct method for powers. I suppose one might be worked out from log and e... but then I found this little clue about 2 ^ 3.5 and got the insight.

I had played around with binary expansion of fractions for years, awed by how patterns of 1's and 0' for rational numbers interweave for instance patterns for n/15 include same patterns for n/3 and n/5 (of course, because they reduce to the same thing but I got to see it played out in graphic diagrams white boxes for 1 and black boxes for 0 a beautiful flow of patterns).

[STxAxTIC]

Hm, for applying Taylor series, first thing that comes to mind is logs and exponentials, exactly like you said. For those who aren't convinced, it goes like this:

i) Suppose we want to compute A = X^Y without explicitly doing the X^Y part.
ii) Take the log of both sides to get: ln A = Y ln X
iii) Exponentiate to land at: A = exp(Y ln X)

... and realize the answer is one Taylor series embedded in another. There are a handful of expansions for ln X in terms of infinite series, but of course these do involve doing integer powers ^2, ^3, ^4, etc. (Easily reduced to trivial multiplication, but you still may accuse me of chicken-and-egging the problem a little bit). Finally, the exp() series is the famous 1 + x + x^2/2 etc.