Precision Math Software

Mind you String Math is doing allot of processing of strings, = slow as the dickens!

zzz means sleeping, waiting for you to press a key...

Option _Explicit
_Title "String Math 2021-06-14" ' b+ from SM2 (2021 June) a bunch of experiments to fix and improve speeds.
' June 2021 fix some old String Math procedures, better nInverse with new LT frunction, remove experimental procedures.
' Now with decent sqrRoot it works independent of Mr$() = Math Regulator that handles signs and decimals and calls to
' add$(), subtr$, mult$, divide$ (100 significant digits),  add$(), subtr$, mult$ are exact!
' If you need higher precsion divide, I recommend use nInverse on denominator (integer)
' then add sign and decimal and mult$() that number with numerator to get divsion answer in higher precision than 100.
' (See how Mr$() handles division and just call nInverse$ with what precision you need.)
' The final function showDP$() is for displaying these number to a set amount of Decimal Places.

' The main code is sampler of tests performed with these functions.

Randomize Timer 'now that it's seems to be running silent
Screen _NewImage(1200, 700, 32)
_Delay .25
_ScreenMove _Middle

'test new stuff
'Dim j As Long
'For j = 1 To 40
'    Print nInverse$(String$(j, "9"), 100)
'Next
'End

Dim r$, ruler$
ruler$ = "0123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890" + Chr$(10)
ruler$ = ruler$ + "0         1         2         3         4         5         6         7         8         9        10        11        12"

' jack error reported 2021-06-04 confirmed!   fixed
Print mr$(".00000000000000000000000000000000000000000000200000000000000054307978001764", "-", ".000000000000000000000000000000000000000000002")
Print ".00000000000000000000000000000000000000000000000000000000000054307978001764"
' debug tests
Print mr$("-5", "+", "-2100"), " OK if -2105"
Print mr$("." + String$(20, "0") + "1", "+", "1" + String$(40, "0") + "1") ' first real test of add
Print ruler$
Print Val("." + String$(20, "0") + "1") + Val("1" + String$(40, "0") + "1") ' test print of big and small val
Print mr$("-.00071", "+", ".00036" + String$(35, "0") + "9")
'-.00071
' .00036000000000000000000000000000000000009
Print "-.00034999999999999999999999999999999999991"

'testing a different subtract sub
Print mr$("10", "-", "5"), " 5 OK"
Print mr$("-10", "+", "5"), " -5 OK"
Print mr$("-10", "-", "-5"), " -5 OK"
Print mr$("-10", "-", "5"), " -15  OK  added"
Print mr$("10", "-", "-5"), " 15  OK  added"
Print mr$("-.010", "-", "-5"), "4.99 OK"
Print mr$("-.010", "-", "5"), "-5.01  OK just added"
Print mr$(".010", "-", "5"), " -4.99 OK"
' jack error reported 2021-06-04 confirmed!  variation below
r$ = mr$(".00000000200000000000000054307978001764", "-", ".0000000020000000000000001") '16 0 wrong  8 wrong
Print "   mr$ rtnd:"; r$
Print "    compare:.00000000000000000000000044307978001764" ' 2021-06-05 finally!
'                  .0000000020000000000000001
'                  .00000000200000000000000054307978001764

r$ = mr$(".00000000000000000100000000000000000011", "-", ".000000000000000001000000000000000001") ' bad too
Print "   mr$ rtnd:"; r$
'         ".000000000000000001000000000000000001"
Print "    compare:-.00000000000000000000000000000000000089"
r$ = mr$(".00000000000000000100000000000000000111", "-", ".000000000000000002000000000000000001")
'       ".000000000000000002000000000000000001"
'       ".00000000000000000100000000000000000111"
'       ".00000000000000000099999999999999999989"
Print "   mr$ rtnd:"; r$
Print "    compare:-.00000000000000000099999999999999999989"

r$ = mr$(".00000000000000000000000000999", "-", "1")
'-1.00000000000000000000000000000000000000000
'  .00000000000000000000000000999
' -.99999999999999999999999999001
Print "   mr$ rtnd:"; r$
Print "    compare:-.99999999999999999999999999001"

r$ = mr$("1", "+", "-1000000000000000000000000000000000000000")
Print "   mr$ rtnd:"; r$
'1000000000000000000000000000000000000000
'-999999999999999999999999999999999999999
Print "    compare:-999999999999999999999999999999999999999"

' check jack problems with FB translation  2021-06-03  errors must be in FB trans from QB64
Print Mid$(mr$(" .1", "/", "3"), 1, 100) ' too long?
Print Mid$(mr$("1.1", "/", "9"), 1, 100)
Print Mid$(mr$("1.38", "/", "1.2"), 1, 100)

' another error reported by jack 2021-06-06 fixed (same problem as subtr$)
Print mr$("1.000000000000000000000001000000000000000000000001", "*", ".000000000000000000000000000000000000000000000001")
Print ".000000000000000000000000000000000000000000000001"
Print "                                                1000000000000000000000001000000000000000000000001"
Print ruler$
Print "checking inverse of 50 = .02 "; nInverse$("50", 20) 'OK  .020000... length 20
Print
Print "zzz... see inverse of STx number now takes 2.2 secs with fixed subtr$() sub,"
Print "can find 115 Fibonacci terms in it."
Sleep
Cls
Dim inverseSTx$, start, done
start = Timer(.001)
inverseSTx$ = nInverse$("999999999999999999999998999999999999999999999999", 2785) ' now 2.2 secs with new subtr from 19 secs
'inverseSTx$ = nInverse2$("999999999999999999999998999999999999999999999999", 817) ' 13 sec damn added 7 secs! now 19.xx
' 816 in 4.22 secs only 64 terms   817 sigDigits matching in 12.93 secs gets all 115 Fibonacci Terms
done = Timer(.001) - start
Print Mid$(inverseSTx$, 1, 3000)
Print
Print "Inverse time:"; done; "   zzz... press any to search for Fibonacci Terms"
Sleep
Dim As Long startSearch, termN, find
Dim f1$, f2$, searchFor$
f1$ = "1"
f2$ = "1"
startSearch = 1
termN = 2
Do
   searchFor$ = mr$(f1$, "+", f2$)
   find = InStr(startSearch, inverseSTx$, searchFor$)
   If find Then
       termN = termN + 1
       Print "Term Number"; termN; " = "; searchFor$; " found at"; find
       f1$ = f2$
       f2$ = searchFor$
       startSearch = find + Len(searchFor$)
   Else
       Print searchFor$; " not found."
       Exit Do
   End If
Loop
Print
' test factorial speed
Dim fact$, i As _Unsigned _Integer64, refFact$, cont$
_KeyClear
Input "Press y for yes, let's do 10000 factorial test, takes quite a bit of time (3.25 mins) "; fact$
If fact$ = "y" Then
   start = Timer(.001)
   fact$ = "1"
   For i = 2 To 10000
       fact$ = TrimLead0$(mult$(fact$, _Trim$(Str$(i))))
       If i Mod 100 = 0 Then Print i; "factorial length ="; Len(fact$)
   Next
   done = Timer(.001) - start
   Print i, Len(fact$), done
   ' save it
   Open "calc 10000!.txt" For Output As #1
   Print #1, fact$
   Close #1
   Beep
   Print Len(fact$), done, "   zzz... press any to compare to reference 10000!."
   Sleep
   _KeyClear
   If _FileExists("10000!.txt") Then
       Open "10000!.txt" For Input As #1
       Input #1, refFact$
       Close #1
       Print "Comparing fact$ to reference fact$:"
       For i = 1 To Len(fact$)
           If Mid$(fact$, i, 1) <> Mid$(refFact$, i, 1) Then
               Print i, Mid$(fact$, i, 1), Mid$(refFact$, i, 1)
               Beep
               Input "Mismatch! Continue? y for yes "; cont$
               If cont$ <> "y" Then Exit For
           End If
       Next
       Print "Compare finished, mismatchs already noted if any."
   Else
       Print "Can't find 10000!.txt reference file."
   End If
   _KeyClear
   Print "zzz... press any to start sqr estimating"
   Sleep
   _KeyClear
End If

Dim n$, result$
Do
   'remember everything is strings
   Input "Enter a number to find it's square root, just enter to quit "; n$
   If n$ = "" Then End
   result$ = sqrRoot$(n$)
   Print result$
   Print "Length ="; Len(result$)
   Print
Loop

' ==  String Math 2021-06-14 Procedure start here (aprox 412 LOC for copy/paste into your app) ==

Function sqrRoot$ (nmbr$)
   Dim n$, guess$, lastGuess$, other$, sum$, imaginary$, loopcnt
   If Left$(nmbr$, 1) = "-" Then 'handle neg numbers
       imaginary$ = "*i": n$ = Mid$(nmbr$, 2)
   Else
       imaginary$ = "": n$ = nmbr$
   End If
   guess$ = mr$(n$, "/", "2")
   other$ = n$
   Do
       loopcnt = loopcnt + 1
       If (Mid$(guess$, 1, 105) = Mid$(lastGuess$, 1, 105)) Then
           ' go past 100 matching digits for 100 digit precision
           sqrRoot$ = Mid$(other$, 1, 101) + imaginary$
           ' try other factor for guess$  sometimes it nails answer without all digits
           Exit Function
       Else
           lastGuess$ = guess$
           sum$ = mr$(guess$, "+", other$)
           guess$ = mr$(sum$, "/", "2")
           other$ = mr$(n$, "/", guess$)
       End If
   Loop
End Function

Function add$ (a$, b$) 'add 2 positive integers assume a and b are just numbers no - signs
   'set a and b numbers to same length and multiple of 18 so can take 18 digits at a time
   Dim As Long la, lb, m, g
   Dim sa As _Unsigned _Integer64, sb As _Unsigned _Integer64, co As _Unsigned _Integer64
   Dim fa$, fb$, t$, new$, result$
   la = Len(a$): lb = Len(b$)
   If la > lb Then m = Int(la / 18) + 1 Else m = Int(lb / 18) + 1
   fa$ = Right$(String$(m * 18, "0") + a$, m * 18)
   fb$ = Right$(String$(m * 18, "0") + b$, m * 18)

   'now taking 18 digits at a time Thanks Steve McNeill
   For g = 1 To m
       sa = Val(Mid$(fa$, (m - g) * 18 + 1, 18))
       sb = Val(Mid$(fb$, (m - g) * 18 + 1, 18))
       t$ = Right$(String$(36, "0") + _Trim$(Str$(sa + sb + co)), 36)
       co = Val(Mid$(t$, 1, 18))
       new$ = Mid$(t$, 19)
       result$ = new$ + result$
   Next
   If co Then result$ = Str$(co) + result$
   add$ = result$
End Function

' This is used in nInverse$ not by Mr$ because there it saves time!
Function subtr1$ (a$, b$)
   Dim As Long la, lb, lResult, i, ca, cb, w
   Dim result$, fa$, fb$

   la = Len(a$): lb = Len(b$)
   If la > lb Then lResult = la Else lResult = lb
   result$ = Space$(lResult)
   fa$ = result$: fb$ = result$
   Mid$(fa$, lResult - la + 1) = a$
   Mid$(fb$, lResult - lb + 1) = b$
   For i = lResult To 1 Step -1
       ca = Val(Mid$(fa$, i, 1))
       cb = Val(Mid$(fb$, i, 1))
       If cb > ca Then ' borrow 10
           Mid$(result$, i, 1) = Right$(Str$(10 + ca - cb), 1)
           w = i - 1
           While w > 0 And Mid$(fa$, w, 1) = "0"
               Mid$(fa$, w, 1) = "9"
               w = w - 1
           Wend
           Mid$(fa$, w, 1) = Right$(Str$(Val(Mid$(fa$, w, 1)) - 1), 1)
       Else
           Mid$(result$, i, 1) = Right$(Str$(ca - cb), 1)
       End If
   Next
   subtr1$ = result$
End Function

' 2021-06-08 fix up with new mr call that decides the sign and puts the greater number first
Function subtr$ (sum$, minus$) ' assume both numbers are positive all digits
   Dim As Long m, g, p
   Dim VB As _Unsigned _Integer64, vs As _Unsigned _Integer64, tenE18 As _Unsigned _Integer64
   Dim ts$, tm$, sign$, LG$, sm$, t$, result$

   ts$ = _Trim$(sum$): tm$ = _Trim$(minus$) ' fixed subtr$ 2021-06-05
   If trim0(ts$) = trim0$(tm$) Then subtr$ = "0": Exit Function 'proceed knowing not equal
   tenE18 = 1000000000000000000 'yes!!! no dang E's
   sign$ = ""
   m = Int(Len(ts$) / 18) + 1
   LG$ = Right$(String$(m * 18, "0") + ts$, m * 18)
   sm$ = Right$(String$(m * 18, "0") + tm$, m * 18)

   'now taking 18 digits at a time From Steve I learned we can do more than 1 digit at a time
   For g = 1 To m
       VB = Val(Mid$(LG$, m * 18 - g * 18 + 1, 18))
       vs = Val(Mid$(sm$, m * 18 - g * 18 + 1, 18))
       If vs > VB Then
           t$ = Right$(String$(18, "0") + _Trim$(Str$(tenE18 - vs + VB)), 18)
           p = (m - g) * 18
           While p > 0 And Mid$(LG$, p, 1) = "0"
               Mid$(LG$, p, 1) = "9"
               p = p - 1
           Wend
           If p > 0 Then Mid$(LG$, p, 1) = _Trim$(Str$(Val(Mid$(LG$, p, 1)) - 1))
       Else
           t$ = Right$(String$(18, "0") + _Trim$(Str$(VB - vs)), 18)
       End If
       result$ = t$ + result$
   Next
   subtr$ = result$
End Function

Function TrimLead0$ (s$) 'for treating strings as number (pos integers)
   Dim copys$
   Dim As Long i, find
   copys$ = _Trim$(s$) 'might as well remove spaces too
   i = 1: find = 0
   While i < Len(copys$) And Mid$(copys$, i, 1) = "0"
       i = i + 1: find = 1
   Wend
   If find = 1 Then copys$ = Mid$(copys$, i)
   If copys$ = "" Then TrimLead0$ = "0" Else TrimLead0$ = copys$
End Function

' catchy? mr$ for math regulator  cop$ = " + - * / " 1 of 4 basic arithmetics
' Fixed so that add and subtract have signs calc'd in Mr and correct call to add or subtract made
' with bigger minus smaller in subtr$() call
Function mr$ (a$, cop$, b$)
   Dim op$, ca$, cb$, aSgn$, bSgn$, postOp$, sgn$, rtn$
   Dim As Long adp, bdp, dp, lpop, aLTb

   op$ = _Trim$(cop$) 'save fixing each time
   ca$ = _Trim$(a$): cb$ = _Trim$(b$) 'make copies in case we change
   'strip signs and decimals
   If Left$(ca$, 1) = "-" Then
       aSgn$ = "-": ca$ = Mid$(ca$, 2)
   Else
       aSgn$ = ""
   End If
   dp = InStr(ca$, ".")
   If dp > 0 Then
       adp = Len(ca$) - dp
       ca$ = Mid$(ca$, 1, dp - 1) + Mid$(ca$, dp + 1)
   Else
       adp = 0
   End If
   If Left$(cb$, 1) = "-" Then
       bSgn$ = "-": cb$ = Mid$(cb$, 2)
   Else
       bSgn$ = ""
   End If
   dp = InStr(cb$, ".")
   If dp > 0 Then
       bdp = Len(cb$) - dp
       cb$ = Mid$(cb$, 1, dp - 1) + Mid$(cb$, dp + 1)
   Else
       bdp = 0
   End If
   If op$ = "+" Or op$ = "-" Or op$ = "/" Then 'add or subtr  even up strings on right of decimal
       'even up the right sides of decimals if any
       If adp > bdp Then dp = adp Else dp = bdp
       If adp < dp Then ca$ = ca$ + String$(dp - adp, "0")
       If bdp < dp Then cb$ = cb$ + String$(dp - bdp, "0")
   ElseIf op$ = "*" Then
       dp = adp + bdp
   End If
   If op$ = "*" Or op$ = "/" Then
       If aSgn$ = bSgn$ Then sgn$ = "" Else sgn$ = "-"
   End If

   'now according to signs and op$ call add$ or subtr$
   If op$ = "-" Then ' make it adding according to signs because that is done for + next!
       If bSgn$ = "-" Then bSgn$ = "" Else bSgn$ = "-" ' flip bSgn$ with op$
       op$ = "+" ' turn this over to + op already done! below
   End If
   If op$ = "+" Then
       If aSgn$ = bSgn$ Then 'really add
           postOp$ = add$(ca$, cb$)
           sgn$ = aSgn$
       ElseIf aSgn$ <> bSgn$ Then 'have a case of subtraction
           'but which is first and which is 2nd and should final sign be pos or neg
           If TrimLead0$(ca$) = TrimLead0(cb$) Then 'remove case a = b
               mr$ = "0": Exit Function
           Else
               aLTb = LTE(ca$, cb$)
               If aSgn$ = "-" Then
                   If aLTb Then ' b - a = pos
                       postOp$ = subtr$(cb$, ca$)
                       sgn$ = ""
                   Else '  a > b so a - sgn wins  - (a - b)
                       postOp$ = subtr$(ca$, cb$)
                       sgn$ = "-"
                   End If
               Else ' b has the - sgn
                   If aLTb Then ' result is -
                       postOp$ = subtr$(cb$, ca$)
                       sgn$ = "-"
                   Else ' result is pos
                       postOp$ = subtr$(ca$, cb$)
                       sgn$ = ""
                   End If
               End If
           End If
       End If
   ElseIf op$ = "*" Then
       postOp$ = mult$(ca$, cb$)
   ElseIf op$ = "/" Then
       postOp$ = divide$(ca$, cb$)
   End If ' which op
   If op$ <> "/" Then 'put dp back
       lpop = Len(postOp$) ' put decimal back if there is non zero stuff following it
       If Len(Mid$(postOp$, lpop - dp + 1)) Then ' fix 1 extra dot appearing in 10000! ?!
           If TrimLead0$(Mid$(postOp$, lpop - dp + 1)) <> "0" Then '  .0   or .00 or .000  ??
               postOp$ = Mid$(postOp$, 1, lpop - dp) + "." + Mid$(postOp$, lpop - dp + 1)
           End If
       End If
   End If
   rtn$ = trim0$(postOp$) 'trim lead 0's then tack on sign
   If rtn$ <> "0" Then mr$ = sgn$ + rtn$ Else mr$ = rtn$
End Function

Function divide$ (n$, d$) ' goal here is 100 digits precision not 100 digits past decimal
   Dim di$, ndi$
   Dim As Long nD
   If n$ = "0" Then divide$ = "0": Exit Function
   If d$ = "0" Then divide$ = "div 0": Exit Function
   If d$ = "1" Then divide$ = n$: Exit Function

   ' aha! found a bug when d$ gets really huge 100 is no where near enough!!!!
   ' 2021-06-03 fix by adding 100 to len(d$), plus have to go a little past 100 like 200
   di$ = Mid$(nInverse$(d$, Len(d$) + 200), 2) 'chop off decimal point after
   nD = Len(di$)
   ndi$ = mult$(n$, di$)
   ndi$ = Mid$(ndi$, 1, Len(ndi$) - nD) + "." + Right$(String$(nD, "0") + Right$(ndi$, nD), nD)
   divide$ = ndi$
End Function

' This uses Subtr1$ is Positive Integer only!
' DP = Decimal places = says when to quit if don't find perfect divisor before
Function nInverse$ (n$, DP As Long) 'assume decimal at very start of the string of digits returned
   Dim m$(1 To 9), si$, r$, outstr$, d$
   Dim i As Long
   For i = 1 To 9
       si$ = _Trim$(Str$(i))
       m$(i) = mult$(si$, n$)
   Next
   outstr$ = ""
   If n$ = "0" Then nInverse$ = "Div 0": Exit Function
   If n$ = "1" Then nInverse$ = "1": Exit Function
   outstr$ = "." 'everything else n > 1 is decimal 8/17
   r$ = "10"
   Do
       While LT(r$, n$) ' 2021-06-08 this should be strictly Less Than
           outstr$ = outstr$ + "0" '            add 0 to the  output string
           If Len(outstr$) = DP + 1 Then nInverse$ = outstr$: Exit Function 'DP length?
           r$ = r$ + "0"
       Wend
       For i = 9 To 1 Step -1
           If LTE(m$(i), r$) Then d$ = _Trim$(Str$(i)): Exit For
       Next
       outstr$ = outstr$ + d$
       If Len(outstr$) = DP + 1 Then nInverse$ = outstr$: Exit Function
       r$ = subtr1$(r$, mult$(d$, n$)) 'r = r -d*n    ' 2021-06-08 subtr1 works faster
       If TrimLead0$(r$) = "0" Then nInverse$ = outstr$: Exit Function ' add trimlead0$ 6/08
       r$ = r$ + "0" 'add another place
   Loop
End Function

Function mult$ (a$, b$) 'assume both positive integers prechecked as all digits strings
   Dim As Long la, lb, m, g, dp
   Dim As _Unsigned _Integer64 v18, sd, co
   Dim f18$, f1$, t$, build$, accum$

   If a$ = "0" Then mult$ = "0": Exit Function
   If b$ = "0" Then mult$ = "0": Exit Function
   If a$ = "1" Then mult$ = b$: Exit Function
   If b$ = "1" Then mult$ = a$: Exit Function
   'find the longer number and make it a mult of 18 to take 18 digits at a time from it
   la = Len(a$): lb = Len(b$)
   If la > lb Then
       m = Int(la / 18) + 1
       f18$ = Right$(String$(m * 18, "0") + a$, m * 18)
       f1$ = b$
   Else
       m = Int(lb / 18) + 1
       f18$ = Right$(String$(m * 18, "0") + b$, m * 18)
       f1$ = a$
   End If
   For dp = Len(f1$) To 1 Step -1 'dp = digit position of the f1$
       build$ = "" 'line builder
       co = 0
       'now taking 18 digits at a time Thanks Steve McNeill
       For g = 1 To m
           v18 = Val(Mid$(f18$, (m - g) * 18 + 1, 18))
           sd = Val(Mid$(f1$, dp, 1))
           t$ = Right$(String$(19, "0") + _Trim$(Str$(v18 * sd + co)), 19)
           co = Val(Mid$(t$, 1, 1))
           build$ = Mid$(t$, 2) + build$
       Next g
       If co Then build$ = _Trim$(Str$(co)) + build$
       If dp = Len(f1$) Then
           accum$ = build$
       Else
           accum$ = add$(accum$, build$ + String$(Len(f1$) - dp, "0"))
       End If
   Next dp
   mult$ = accum$
End Function

'this function needs TrimLead0$(s$)
Function LTE (a$, b$) ' a$ Less Than or Equal b$  comparison of 2 strings
   Dim ca$, cb$
   Dim As Long la, lb, i
   ca$ = TrimLead0$(a$): cb$ = TrimLead0(b$)
   la = Len(ca$): lb = Len(cb$)
   If ca$ = cb$ Then
       LTE = -1
   ElseIf la < lb Then ' a is smaller
       LTE = -1
   ElseIf la > lb Then ' a is bigger
       LTE = 0
   ElseIf la = lb Then ' equal lengths
       For i = 1 To Len(ca$)
           If Val(Mid$(ca$, i, 1)) > Val(Mid$(cb$, i, 1)) Then
               LTE = 0: Exit Function
           ElseIf Val(Mid$(ca$, i, 1)) < Val(Mid$(cb$, i, 1)) Then
               LTE = -1: Exit Function
           End If
       Next
   End If
End Function

'need this for ninverse faster than subtr$ for sign
Function LT (a$, b$) ' a$ Less Than or Equal b$  comparison of 2 strings
   Dim ca$, cb$
   Dim As Long la, lb, i
   ca$ = TrimLead0$(a$): cb$ = TrimLead0(b$)
   la = Len(ca$): lb = Len(cb$)
   If la < lb Then ' a is smaller
       LT = -1
   ElseIf la > lb Then ' a is bigger
       LT = 0
   ElseIf la = lb Then ' equal lengths
       For i = 1 To Len(ca$)
           If Val(Mid$(ca$, i, 1)) > Val(Mid$(cb$, i, 1)) Then
               LT = 0: Exit Function
           ElseIf Val(Mid$(ca$, i, 1)) < Val(Mid$(cb$, i, 1)) Then
               LT = -1: Exit Function
           End If
       Next
   End If
End Function

Function TrimTail0$ (s$)
   Dim copys$
   Dim As Long dp, i, find
   copys$ = _Trim$(s$) 'might as well remove spaces too
   TrimTail0$ = copys$
   dp = InStr(copys$, ".")
   If dp > 0 Then
       i = Len(copys$): find = 0
       While i > dp And Mid$(copys$, i, 1) = "0"
           i = i - 1: find = 1
       Wend
       If find = 1 Then
           If i = dp Then
               TrimTail0$ = Mid$(copys$, 1, dp - 1)
           Else
               TrimTail0$ = Mid$(copys$, 1, i)
           End If
       End If
   End If
End Function

Function trim0$ (s$)
   Dim cs$, si$
   cs$ = s$
   si$ = Left$(cs$, 1)
   If si$ = "-" Then cs$ = Mid$(cs$, 2)
   cs$ = TrimLead0$(cs$)
   cs$ = TrimTail0$(cs$)
   If Right$(cs$, 1) = "." Then cs$ = Mid$(cs$, 1, Len(cs$) - 1)
   If si$ = "-" Then trim0$ = si$ + cs$ Else trim0$ = cs$
End Function

' for displaying truncated numbers say to 60 digits
Function showDP$ (num$, nDP As Long)
   Dim cNum$
   Dim As Long dp, d, i
   cNum$ = num$ 'since num$ could get changed
   showDP$ = num$
   dp = InStr(num$, ".")
   If dp > 0 Then
       If Len(Mid$(cNum$, dp + 1)) > nDP Then
           d = Val(Mid$(cNum$, dp + nDP + 1, 1))
           If d > 4 Then
               cNum$ = "0" + cNum$ ' tack on another 0 just in case 9's all the way to left
               dp = dp + 1
               i = dp + nDP
               While Mid$(cNum$, i, 1) = "9" Or Mid$(cNum$, i, 1) = "."
                   If Mid$(cNum$, i, 1) = "9" Then
                       Mid$(cNum$, i, 1) = "0"
                   End If
                   i = i - 1
               Wend
               Mid$(cNum$, i, 1) = _Trim$(Str$(Val(Mid$(cNum$, i, 1)) + 1)) 'last non 9 digit
               cNum$ = Mid$(cNum$, 1, dp + nDP) 'chop it
               showDP$ = trim0$(cNum$)
           Else
               showDP$ = Mid$(cNum$, 1, dp + nDP)
           End If
       End If
   End If
End Function

Later version works with powers and Binary number.

Yes the functions Add, subtract,... operators handle integers but the Mr$ function handles the negative sign and decimal points in the 2 operands and calls the appropriate integer function gets the answer and then adds the sign and decimal point.

eg
answer$ = mr$(".123456789", "+", "-123456789")
Print answer$
_Clipboard$ = answer$ ' -123456788.876543211

So these are done in 2 step process, convert numbers to integer, do the math, convert integer back to proper sign and decimal.

SqrRoot$ uses an approximation method to get to Square Root with Add, Subtract, Mult, Divide through Mr$( ) calls.

Val(ans$) is likely going to undo all the work of String Math with super huge or super small numbers. BTW I am not a fan of Scientific Notation which is one main reason I went through all the trouble of String math, the other is arbitrary precision (for division you have to make arrangements to other than 100 digit with nInverse() Function like I did with Stx special number that gets 115 terms of Fibonacci series in a single decimal expansion, it's in the demo.)

Once you go string with a number you pretty much have to stay string, no more Val(). That is why the functions only take String Type Operands.

Also if
Dim as Long a, b, c doesn't make a, b, c Long then your version of QB64.exe is way out of date. After version 2.0 we have to be careful using Function names as variables in Function definitions that was Major Update to way Functions are defined and will not be compatable with old code same as Dim As (Type) list a, b, c,...

"Once you go string with a number you pretty much have to stay string, no more Val()."


I see your point. So for my fractal programming I can see I'll have to write some routines that let me find 'greater than', 'less than','equal to' etc. It's like another mini-language.

Ahh, thanks again. I hadn't noticed those two very eloquently devised functions...

Regards-

Honestly for drawing fractals default single precision should be fine. Pixels round to nearest integer:

_Title "Sierpinski flies a kite by bplus 2017-10-16"
' after playing with Ashish Kite Fractal

Screen _NewImage(1200, 700, 32)
_ScreenMove 100, 20
While 1
   Cls
   drawKite 600, 540, 200, a
   _Display
   _Limit 20
   a = a + _Pi(2 / 360)
Wend
Sleep
Sub drawKite (xx, yy, s, a)
   x = xx: y = yy
   x2 = x + 3 * s * Cos(_Pi(1 / 2) - a / 2): y2 = y + 3 * s * Sin(_Pi(1 / 2) - a / 2)
   x3 = x + 3 * s * Cos(_Pi(1 / 2) + a / 2): y3 = y + 3 * s * Sin(_Pi(1 / 2) + a / 2)
   SierLineTri x, y, x2, y2, x3, y3, 0
   'LINE (x, y)-(x + s * COS(_PI(2) - a / 2), (y - s) + s * SIN(_PI(2) - a / 2))
   'LINE (x, y)-(x + s * COS(_PI + a / 2), (y - s) + s * SIN(_PI + a / 2))

   If s > 10 Then
       drawKite x + 1 * s * Cos(_Pi(2) - a), (y - s) + 1 * s * Sin(_Pi(2) - a), s / 2, a
       drawKite x + 1 * s * Cos(_Pi + a), (y - s) + 1 * s * Sin(_Pi + a), s / 2, a
   End If
End Sub
Sub SierLineTri (x1, y1, x2, y2, x3, y3, depth)
   If depth = 0 Then 'draw out triangle if level 0
       Line (x1, y1)-(x2, y2)
       Line (x2, y2)-(x3, y3)
       Line (x1, y1)-(x3, y3)
   End If
   'find midpoints
   If x2 < x1 Then mx1 = (x1 - x2) / 2 + x2 Else mx1 = (x2 - x1) / 2 + x1
   If y2 < y1 Then my1 = (y1 - y2) / 2 + y2 Else my1 = (y2 - y1) / 2 + y1
   If x3 < x2 Then mx2 = (x2 - x3) / 2 + x3 Else mx2 = (x3 - x2) / 2 + x2
   If y3 < y2 Then my2 = (y2 - y3) / 2 + y3 Else my2 = (y3 - y2) / 2 + y2
   If x3 < x1 Then mx3 = (x1 - x3) / 2 + x3 Else mx3 = (x3 - x1) / 2 + x1
   If y3 < y1 Then my3 = (y1 - y3) / 2 + y3 Else my3 = (y3 - y1) / 2 + y1

   Line (mx1, my1)-(mx2, my2) '  'draw all inner triangles
   Line (mx2, my2)-(mx3, my3)
   Line (mx1, my1)-(mx3, my3)

   If depth < 4 Then 'not done so call me again
       SierLineTri x1, y1, mx1, my1, mx3, my3, depth + 1
       SierLineTri x2, y2, mx1, my1, mx2, my2, depth + 1
       SierLineTri x3, y3, mx3, my3, mx2, my2, depth + 1
   End If
End Sub


Yes, I was recently checking out some programs in the Samples section from this Forum's Home Page.

I was very impressed with this!
'> Merged with Zom-B's smart $include merger 0.51

DefDbl A-Z

'####################################################################################################################
'# Math Library V1.0 (include)
'# By Zom-B
'####################################################################################################################

Const sqrt2 = 1.41421356237309504880168872420969807856967187537695 ' Knuth01
Const sqrt3 = 1.73205080756887729352744634150587236694280525381038 ' Knuth02
Const sqrt5 = 2.23606797749978969640917366873127623544061835961153 ' Knuth03
Const sqrt10 = 3.16227766016837933199889354443271853371955513932522 ' Knuth04
Const cubert2 = 1.25992104989487316476721060727822835057025146470151 ' Knuth05
Const cubert3 = 1.44224957030740838232163831078010958839186925349935 ' Knuth06
Const q2pow025 = 1.18920711500272106671749997056047591529297209246382 ' Knuth07
Const phi = 1.61803398874989484820458683436563811772030917980576 ' Knuth08
Const log2 = 0.69314718055994530941723212145817656807550013436026 ' Knuth09
Const log3 = 1.09861228866810969139524523692252570464749055782275 ' Knuth10
Const log10 = 2.30258509299404568401799145468436420760110148862877 ' Knuth11
Const logpi = 1.14472988584940017414342735135305871164729481291531 ' Knuth12
Const logphi = 0.48121182505960344749775891342436842313518433438566 ' Knuth13
Const q1log2 = 1.44269504088896340735992468100189213742664595415299 ' Knuth14
Const q1log10 = 0.43429448190325182765112891891660508229439700580367 ' Knuth15
Const q1logphi = 2.07808692123502753760132260611779576774219226778328 ' Knuth16
Const pi = 3.14159265358979323846264338327950288419716939937511 ' Knuth17
Const deg2rad = 0.01745329251994329576923690768488612713442871888542 ' Knuth18
Const q1pi = 0.31830988618379067153776752674502872406891929148091 ' Knuth19
Const pisqr = 9.86960440108935861883449099987615113531369940724079 ' Knuth20
Const gamma05 = 1.7724538509055160272981674833411451827975494561224 '  Knuth21
Const gamma033 = 2.6789385347077476336556929409746776441286893779573 '  Knuth22
Const gamma067 = 1.3541179394264004169452880281545137855193272660568 '  Knuth23
Const e = 2.71828182845904523536028747135266249775724709369996 ' Knuth24
Const q1e = 0.36787944117144232159552377016146086744581113103177 ' Knuth25
Const esqr = 7.38905609893065022723042746057500781318031557055185 ' Knuth26
Const eulergamma = 0.57721566490153286060651209008240243104215933593992 ' Knuth27
Const expeulergamma = 1.7810724179901979852365041031071795491696452143034 '  Knuth28
Const exppi025 = 2.19328005073801545655976965927873822346163764199427 ' Knuth29
Const sin1 = 0.84147098480789650665250232163029899962256306079837 ' Knuth30
Const cos1 = 0.54030230586813971740093660744297660373231042061792 ' Knuth31
Const zeta3 = 1.2020569031595942853997381615114499907649862923405 '  Knuth32
Const nloglog2 = 0.36651292058166432701243915823266946945426344783711 ' Knuth33

Const logr10 = 0.43429448190325182765112891891660508229439700580367
Const logr2 = 1.44269504088896340735992468100189213742664595415299
Const pi05 = 1.57079632679489661923132169163975144209858469968755
Const pi2 = 6.28318530717958647692528676655900576839433879875021
Const q05log10 = 0.21714724095162591382556445945830254114719850290183
Const q05log2 = 0.72134752044448170367996234050094606871332297707649
Const q05pi = 0.15915494309189533576888376337251436203445964574046
Const q13 = 0.33333333333333333333333333333333333333333333333333
Const q16 = 0.16666666666666666666666666666666666666666666666667
Const q2pi = 0.63661977236758134307553505349005744813783858296183
Const q2sqrt5 = 0.89442719099991587856366946749251049417624734384461
Const rad2deg = 57.2957795130823208767981548141051703324054724665643
Const sqrt02 = 0.44721359549995793928183473374625524708812367192231
Const sqrt05 = 0.70710678118654752440084436210484903928483593768847
Const sqrt075 = 0.86602540378443864676372317075293618347140262690519
Const y2q112 = 1.05946309435929526456182529494634170077920431749419 ' Chromatic base

'####################################################################################################################
'# Screen mode selector v1.0 (include)
'# By Zom-B
'####################################################################################################################

videoaspect:
Data "all aspect",15
Data "4:3",11
Data "16:10",10
Data "16:9",14
Data "5:4",13
Data "3:2",12
Data "5:3",9
Data "1:1",7
Data "other",8
Data ,

videomodes:
Data 256,256,7
Data 320,240,1
Data 400,300,1
Data 512,384,1
Data 512,512,7
Data 640,480,1
Data 720,540,1
Data 768,576,1
Data 800,480,2
Data 800,600,1
Data 854,480,3
Data 1024,600,8
Data 1024,640,2
Data 1024,768,1
Data 1024,1024,7
Data 1152,768,5
Data 1152,864,1
Data 1280,720,3
Data 1280,768,6
Data 1280,800,2
Data 1280,854,5
Data 1280,960,1
Data 1280,1024,4
Data 1366,768,3
Data 1400,1050,1
Data 1440,900,2
Data 1440,960,5
Data 1600,900,3
Data 1600,1200,1
Data 1680,1050,2
Data 1920,1080,3
Data 1920,1200,2
Data 2048,1152,3
Data 2048,1536,1
Data 2048,2048,7
Data ,,

'####################################################################################################################
'# Ultra Fractal Gradient library v1.0 (include)
'# By Zom-B
'#
'# Smooth Gradient algorithm from Ultra Fractal (www.ultrafractal.com)
'####################################################################################################################

Type GRADIENTPOINT
   index As Single
   r As Single
   g As Single
   b As Single
   rdr As Single
   rdl As Single
   gdr As Single
   gdl As Single
   bdr As Single
   bdl As Single
End Type

'$dynamic

Dim Shared gradientSmooth(1) As _Byte '_BIT <- bugged
Dim Shared gradientPoints(1) As Integer
Dim Shared gradient(1, 1) As GRADIENTPOINT

'####################################################################################################################
'# Ultra Fractal Gradient library v1.0 (include)
'# By Zom-B
'#
'# Smooth Gradient algorithm from Ultra Fractal (www.ultrafractal.com)
'####################################################################################################################

Type OPACITYPOINT
   index As Single
   o As Single
   odr As Single
   odl As Single
End Type

'$dynamic

Dim Shared opacitySmooth(0 To 0) As _Byte '_BIT <- bugged
Dim Shared opacityPoints(0 To 0) As Integer
Dim Shared opacity(0 To 0, 0 To 0) As OPACITYPOINT


'####################################################################################################################
'# InterestingSpiral2
'# By Zom-B
'#
'# Original art by Mark Hammond (markch1@mindspring.com) (Sep 10, 2002)
'####################################################################################################################

Const Doantialias = -1
Const Usegaussian = 0

'####################################################################################################################

_Title "InterestingSpiral2"
Width 80, 40

Color 11
Print
Print Tab(31); "InterestingSpiral2"
Color 7
Print
Print Tab(6); "Original art by Mark Hammond (markch1@mindspring.com) (Sep 10, 2002)"
Print Tab(19); "Converted to Quick Basic and QB64 by Zom-B"

selectScreenMode 6, 32

'####################################################################################################################

Dim Shared sizeX%, sizeY%
Dim Shared maxX%, maxY%
Dim Shared halfX%, halfY%

sizeX% = _Width
sizeY% = _Height
maxX% = sizeX% - 1
maxY% = sizeY% - 1
halfX% = sizeX% \ 2
halfY% = sizeY% \ 2

Dim Shared zx(250), zy(250)
Dim Shared px, py
Dim Shared magn

magn = 0.375 / halfY%

'####################################################################################################################

setNumGradients 2
setNumOpacities 2

addGradientPoint 0, -0.1150, 0.780, 0.553, 0.420
addGradientPoint 0, 0.0450, 0.169, 0.000, 0.000
addGradientPoint 0, 0.2175, 0.992, 0.878, 0.741
addGradientPoint 0, 0.3850, 0.169, 0.000, 0.000
addGradientPoint 0, 0.5525, 0.725, 0.694, 0.600
addGradientPoint 0, 0.7200, 0.169, 0.000, 0.000
addGradientPoint 0, 0.8850, 0.780, 0.553, 0.420
addGradientPoint 0, 1.0450, 0.169, 0.000, 0.000
setGradientSmooth 0, -1

addGradientPoint 1, -0.140, 0.502, 0.251, 0.063
addGradientPoint 1, 0.110, 0.502, 0.251, 0.063
addGradientPoint 1, 0.235, 0.706, 0.502, 0.251
addGradientPoint 1, 0.360, 0.867, 0.749, 0.561
addGradientPoint 1, 0.485, 1.000, 1.000, 1.000
addGradientPoint 1, 0.610, 0.867, 0.749, 0.561
addGradientPoint 1, 0.735, 0.706, 0.502, 0.251
addGradientPoint 1, 0.860, 0.502, 0.251, 0.063
addGradientPoint 1, 1.110, 0.502, 0.251, 0.063
setGradientSmooth 1, -1

addOpacityPoint 0, -0.210, 0
addOpacityPoint 0, 0.285, 1
addOpacityPoint 0, 0.790, 0
addOpacityPoint 0, 1.285, 1
setOpacitySmooth 0, -1

addOpacityPoint 1, -0.0025, 0
addOpacityPoint 1, 0.4975, 1
addOpacityPoint 1, 0.9975, 0
addOpacityPoint 1, 1.4975, 1
setOpacitySmooth 1, -1

renderProgressive 256, 4

i$ = Input$(1)
End

'####################################################################################################################

Sub renderProgressive (startSize%, endSize%)
   pixStep% = startSize%

   pixWidth% = pixStep% - 1
   For y% = 0 To maxY% Step pixStep%
       For x% = 0 To maxX% Step pixStep%
           calcPoint x%, y%, r%, g%, b%
           Line (x%, y%)-Step(pixWidth%, pixWidth%), _RGB(r%, g%, b%), BF
       Next
       If InKey$ = Chr$(27) Then System
   Next

   Do
       pixSize% = pixStep% \ 2
       pixWidth% = pixSize% - 1
       For y% = 0 To maxY% Step pixStep%
           y1% = y% + pixSize%
           For x% = 0 To maxX% Step pixStep%
               x1% = x% + pixSize%

               If x1% < sizeX% Then
                   calcPoint x1%, y%, r%, g%, b%
                   Line (x1%, y%)-Step(pixWidth%, pixWidth%), _RGB(r%, g%, b%), BF
               End If
               If y1% < sizeY% Then
                   calcPoint x%, y1%, r%, g%, b%
                   Line (x%, y1%)-Step(pixWidth%, pixWidth%), _RGB(r%, g%, b%), BF
                   If x1% < sizeX% Then
                       calcPoint x1%, y1%, r%, g%, b%
                       Line (x1%, y1%)-Step(pixWidth%, pixWidth%), _RGB(r%, g%, b%), BF
                   End If
               End If
           Next
           If InKey$ = Chr$(27) Then System
       Next
       pixStep% = pixStep% \ 2
   Loop While pixStep% > 2

   For y% = 0 To maxY% Step 2
       y1% = y% + 1
       For x% = 0 To maxX% Step 2
           x1% = x% + 1

           If x1% < sizeX% Then
               calcPoint x1%, y%, r%, g%, b%
               PSet (x1%, y%), _RGB(r%, g%, b%)
           End If
           If y1% < sizeY% Then
               calcPoint x%, y1%, r%, g%, b%
               PSet (x%, y1%), _RGB(r%, g%, b%)
               If x1% < sizeX% Then
                   calcPoint x1%, y1%, r%, g%, b%
                   PSet (x1%, y1%), _RGB(r%, g%, b%)
               End If
           End If
       Next
       If InKey$ = Chr$(27) Then System
   Next

   If Not Doantialias Then Exit Sub

   endArea% = endSize% * endSize%

   If Usegaussian Then
       For y% = 0 To maxY%
           For x% = 0 To maxX%
               c& = Point(x%, y%)
               r% = _Red(c&)
               g% = _Green(c&)
               b% = _Blue(c&)
               For i% = 2 To endArea%
                   Do 'Marsaglia polar method for random gaussian
                       u! = Rnd * 2 - 1
                       v! = Rnd * 2 - 1
                       s! = u! * u! + v! * v!
                   Loop While s! >= 1 Or s! = 0
                   s! = Sqr(-2 * Log(s!) / s!) * 0.5
                   u! = u! * s!
                   v! = v! * s!

                   calcPoint x% + u!, y% + v!, r1%, g1%, b1%

                   r% = r% + r1%
                   g% = g% + g1%
                   b% = b% + b1%
               Next

               PSet (x%, y%), _RGB(CInt(r% / endArea%), CInt(g% / endArea%), CInt(b% / endArea%))
               If InKey$ = Chr$(27) Then System
           Next
       Next
   Else
       For y% = 0 To maxY%
           For x% = 0 To maxX%
               r% = 0
               g% = 0
               b% = 0
               For v% = 0 To endSize% - 1
                   y1! = y% + v% / endSize%
                   For u% = 0 To endSize% - 1
                       If u% = 0 And v& = 0 Then
                           c& = Point(x%, y%)
                       Else
                           x1! = x% + u% / endSize%
                           calcPoint x1!, y1!, r1%, g1%, b1%
                       End If
                       r% = r% + r1%
                       g% = g% + g1%
                       b% = b% + b1%
                   Next
               Next
               PSet (x%, y%), _RGB(CInt(r% / endArea%), CInt(g% / endArea%), CInt(b% / endArea%))
               If InKey$ = Chr$(27) Then System
           Next
       Next
   End If
End Sub

'####################################################################################################################

Sub calcPoint (screenX!, screenY!, r%, g%, b%)
   applyLocation screenX!, screenY!

   calcFractal numIter%

   If numIter% < 250 Then
       calcOutside numIter%, index1!, index2!, index3!, index4!

       getGradient 0, index1!, r!, g!, b!
       getGradient 0, index2!, r2!, g2!, b2!
       getGradient 0, index3!, r3!, g3!, b3!
       getOpacity 0, index3!, o3!
       getGradient 1, index4!, r4!, g4!, b4!
       getOpacity 1, index4!, o4!

       r! = r! - r2!: If r! < 0 Then r! = 0
       g! = g! - g2!: If g! < 0 Then g! = 0
       b! = b! - b2!: If b! < 0 Then b! = 0

       r2! = r!: g2! = g!: b2! = b!
       mergeHardLight r2!, g2!, b2!, r3!, g3!, b3!
       r! = r! + (r2! - r!) * o3!
       g! = g! + (g2! - g!) * o3!
       b! = b! + (b2! - b!) * o3!

       r2! = r!: g2! = g!: b2! = b!
       mergeSoftLight r2!, g2!, b2!, r4!, g4!, b4!
       r! = r! + (r2! - r!) * o4!
       g! = g! + (g2! - g!) * o4!
       b! = b! + (b2! - b!) * o4!
   Else
       'r! = 0: g! = 0: b! = 0
   End If

   r% = r! * 255
   g% = g! * 255
   b% = b! * 255
End Sub

'####################################################################################################################

Sub applyLocation (inX!, inY!)
   px = (inX! - halfX%) * magn
   py = (halfY% - inY!) * magn
End Sub

'####################################################################################################################

Sub calcFractal (numIter%)
   zx(0) = px: zy(0) = py
   x = px: y = py
   xx = x * x: yy = y * y

   i% = -1
   For numIter% = 1 To 250
       If i% Then
           y = 2 * x * y + 0.5276360735394
           x = xx - yy + 0.5494321598602
       Else
           t = x - xx + yy
           y = y - 2 * y * x
           x = t * 1.5875 + y * .51875
           y = y * 1.5875 - t * .51875
       End If
       i% = Not i%

       zx(numIter%) = x: zy(numIter%) = y

       xx = x * x: yy = y * y

       If xx + yy >= 128 Then
           Outside% = true
           Exit For
       End If
   Next

   numIter% = numIter% - 1
End Sub

'####################################################################################################################

Sub calcOutside (numIter%, index1!, index2!, index3!, index4!)
   dist1 = 1E38: dist2 = 1E38: dist3 = 0: dist4 = 1E+38

   For a% = 1 To numIter%
       x = zx(a%): y = zy(a%)
       r = Sqr(x * x + y * y)

       If r >= 0.1 And r <= 0.75 Then
           a = 4 * r * r
           x2 = x - Cos(a) * r
           y2 = y - Sin(a) * r
           a = x2 * x2 + y2 * y2
           If dist4 > a Then dist4 = a
       End If

       a = atan2(Abs(y), Abs(x))
       If a < 0 Then a = a + pi2
       r = r + a

       If dist2 > r Then dist2 = r

       x2 = x - r * r + 0.25
       y2 = y - r
       r = x2 * x2 + y2 * y2

       If dist1 > r Then dist1 = r

       x2 = px * px + py * py

       y2 = Int(0.5 - (y * px - x * py) / x2)
       x2 = Int(0.5 - (x * px + y * py) / x2)

       x = px * x2 - py * y2 + x
       y = px * y2 + py * x2 + y

       a = x * x + y * y
       If dist3 < a Then dist3 = a
   Next
   cAsin Sqr(dist1) ^ .1, 0, zr, zi
   index1! = (zr * zr + zi * zi) + 0.865

   cAsin dist2, 0, zr, zi
   index2! = Sqr(zr * zr + zi * zi) / 2 + 0.985

   index3! = Sqr(Sqr(dist3)) / 2 + 0.37

   If dist4 = 1E38 Then index4! = 0.388 Else index4! = Sqr(Sqr(dist4)) / 2
End Sub

'####################################################################################################################

Sub cAsin (re, im, outRe, outIm)
   ' = Inverse Sine = LOG(sqrt(1-z^2) + z*i) * -i
   a = 1 - re * re + im * im
   b = -2 * re * im
   c = Sqr(a * a + b * b)
   If b < 0 Then
       b = re - Sqr((c - a) * 0.5)
   Else
       b = re + Sqr((c - a) * 0.5)
   End If
   a = Sqr((c + a) * 0.5) - im
   outRe = atan2(b, a)
   outIm = Log(a * a + b * b) * -0.5
End Sub

'####################################################################################################################
'# Math Library V0.11 (routines)
'# By Zom-B
'####################################################################################################################

'> merger: Skipping unused FUNCTION remainder% (a%, b%)

'> merger: Skipping unused FUNCTION fRemainder (a, b)

'####################################################################################################################

'> merger: Skipping unused FUNCTION safeLog (x)

'####################################################################################################################

'> merger: Skipping unused FUNCTION asin (y)

'> merger: Skipping unused FUNCTION acos (y)

'> merger: Skipping unused FUNCTION safeAcos (y)

Function atan2 (y, x)
   If x > 0 Then
       atan2 = Atn(y / x)
   ElseIf x < 0 Then
       If y > 0 Then
           atan2 = Atn(y / x) + _Pi
       Else
           atan2 = Atn(y / x) - _Pi
       End If
   ElseIf y > 0 Then
       atan2 = _Pi / 2
   Else
       atan2 = -_Pi / 2
   End If
End Function

'####################################################################################################################
'# Screen mode selector v1.0 (routines)
'# By Zom-B
'####################################################################################################################

Sub selectScreenMode (yOffset%, colors%)
   Dim aspectName$(10), aspectCol%(10)
   Restore videoaspect
   For y% = 0 To 10
       Read aspectName$(y%), aspectCol%(y%)
       If aspectCol%(y%) = 0 Then numAspect% = y% - 1: Exit For
   Next

   Dim vidX%(100), vidY%(100), vidA%(100)
   Restore videomodes
   For y% = 1 To 100
       Read vidX%(y%), vidY%(y%), vidA%(y%)
       If (vidX%(y%) <= 0) Then numModes% = y% - 1: Exit For
   Next

   If numModes% > _Height - yOffset% - 1 Then numModes% = _Height - yOffset% - 1

   Def Seg = &HB800
   Locate yOffset% + 1, 1
   Print "Select video mode:"; Tab(61); "Click "
   Poke yOffset% * 160 + 132, 31

   y% = 0
   lastY% = 0
   selectedAspect% = 0
   reprint% = 1
   lastButton% = 0
   Do
       If InKey$ = Chr$(27) Then Cls: System
       If reprint% Then
           reprint% = 0

           For x% = 1 To numModes%
               Locate yOffset% + x% + 1, 1
               Color 7, 0
               Print Using "##:"; x%;
               If selectedAspect% = 0 Then
                   Color aspectCol%(vidA%(x%))
               ElseIf selectedAspect% = vidA%(x%) Then
                   Color 15
               Else
                   Color 8
               End If
               Print Str$(vidX%(x%)); ","; vidY%(x%);
           Next

           For x% = 0 To numAspect%
               If x% > 0 And selectedAspect% = x% Then
                   Color aspectCol%(x%), 3
               Else
                   Color aspectCol%(x%), 0
               End If
               Locate yOffset% + x% + 2, 64
               Print "<"; aspectName$(x%); ">"
           Next
       End If
       If _MouseInput Then
           If lastY% > 0 Then
               For x% = 0 To 159 Step 2
                   Poke lastY% + x%, Peek(lastY% + x%) And &HEF
               Next
           End If

           x% = _MouseX
           y% = _MouseY - yOffset% - 1

           If x% <= 60 Then
               If y% > 0 And y% <= numModes% Then
                   If _MouseButton(1) = 0 And lastButton% Then Exit Do
                   y% = (yOffset% + y%) * 160 + 1
                   For x% = 0 To 119 Step 2
                       Poke y% + x%, Peek(y% + x%) Or &H10
                   Next
               Else
                   y% = 0
               End If
           Else
               If y% > 0 And y% - 1 <= numAspect% Then
                   If _MouseButton(1) Then
                       selectedAspect% = y% - 1
                       reprint% = 1
                   End If
                   y% = (yOffset% + y%) * 160 + 1
                   For x% = 120 To 159 Step 2
                       Poke y% + x%, Peek(y% + x%) Or &H10
                   Next
               Else
                   y% = 0
               End If
           End If
           lastY% = y%
           lastButton% = _MouseButton(1)
       End If
   Loop

   Cls 'bug evasion for small video modes
   Screen _NewImage(vidX%(y%), vidY%(y%), colors%)
End Sub

'####################################################################################################################
'# Ultra Fractal Gradient library v1.0 (routines)
'# By Zom-B
'#
'# Smooth Gradient algorithm from Ultra Fractal (www.ultrafractal.com)
'####################################################################################################################

'> merger: Skipping unused SUB defaultGradient (gi%)

'> merger: Skipping unused SUB grayscaleGradient (gi%)

'####################################################################################################################

Sub setNumGradients (gi%)
   offset% = LBound(gradientPoints) - 1
   ReDim _Preserve gradientSmooth(gi% + offset%) As _Byte '_BIT <- bugged
   ReDim _Preserve gradientPoints(gi% + offset%) As Integer
   ReDim _Preserve gradient(gi% + offset%, 1) As GRADIENTPOINT
End Sub

Sub addGradientPoint (gi%, index!, r!, g!, b!)
   p% = gradientPoints(gi%)

   If UBound(gradient, 2) < p% Then
       ReDim _Preserve gradient(0 To UBound(gradient, 1), 0 To p%) As GRADIENTPOINT
   End If

   gradient(gi%, p%).index = index!
   gradient(gi%, p%).r = r!
   gradient(gi%, p%).g = g!
   gradient(gi%, p%).b = b!
   gradientPoints(gi%) = p% + 1
End Sub

Sub setGradientSmooth (gi%, s%)
   gradientSmooth(gi%) = s%

   If gradientSmooth(0) = 0 Then Exit Sub

   For i% = 1 To gradientPoints(gi%) - 1
       i1% = i% + 1
       If i1% = gradientPoints(gi%) Then i1% = 2

       dxl! = gradient(gi%, i%).index - gradient(gi%, i% - 1).index
       dxr! = gradient(gi%, i1%).index - gradient(gi%, i%).index
       If dxl! < 0 Then dxl! = dxl! + 1
       If dxr! < 0 Then dxr! = dxr! + 1

       d! = (gradient(gi%, i%).r - gradient(gi%, i% - 1).r) * dxr!
       If d! = 0 Then
           gradient(gi%, i%).rdr = 0
           gradient(gi%, i%).rdl = 0
       Else
           d! = (gradient(gi%, i1%).r - gradient(gi%, i%).r) * dxl! / d!
           If d! <= 0 Then
               gradient(gi%, i%).rdr = 0
               gradient(gi%, i%).rdl = 0
           Else
               gradient(gi%, i%).rdr = 1 / (1 + d!)
               gradient(gi%, i%).rdl = gradient(gi%, i%).rdr - 1
           End If
       End If

       d! = (gradient(gi%, i%).g - gradient(gi%, i% - 1).g) * dxr!
       If d! = 0 Then
           gradient(gi%, i%).gdr = 0
           gradient(gi%, i%).gdl = 0
       Else
           d! = (gradient(gi%, i1%).g - gradient(gi%, i%).g) * dxl! / d!
           If d! <= 0 Then
               gradient(gi%, i%).gdr = 0
               gradient(gi%, i%).gdl = 0
           Else
               gradient(gi%, i%).gdr = 1 / (1 + d!)
               gradient(gi%, i%).gdl = gradient(gi%, i%).gdr - 1
           End If
       End If

       d! = (gradient(gi%, i%).b - gradient(gi%, i% - 1).b) * dxr!
       If d! = 0 Then
           gradient(gi%, i%).bdr = 0
           gradient(gi%, i%).bdl = 0
       Else
           d! = (gradient(gi%, i1%).b - gradient(gi%, i%).b) * dxl! / d!
           If d! <= 0 Then
               gradient(gi%, i%).bdr = 0
               gradient(gi%, i%).bdl = 0
           Else
               gradient(gi%, i%).bdr = 1 / (1 + d!)
               gradient(gi%, i%).bdl = gradient(gi%, i%).bdr - 1
           End If
       End If
   Next
End Sub

'####################################################################################################################

Sub getGradient (gi%, index!, red!, green!, blue!)
   If index! < 0 Then x! = 0 Else x! = index! - Int(index!)

   For l% = gradientPoints(gi%) - 2 To 1 Step -1
       If gradient(gi%, l%).index <= x! Then
           Exit For
       End If
   Next

   r% = l% + 1
   u! = (x! - gradient(gi%, l%).index) / (gradient(gi%, r%).index - gradient(gi%, l%).index)

   If gradientSmooth(gi%) Then
       u2! = u! * u!
       u3! = u2! * u!
       ur! = u3! - (u2! + u2!) + u!
       ul! = u2! - u3!

       red! = gradient(gi%, l%).r + (gradient(gi%, r%).r - gradient(gi%, l%).r) * (u3! + 3 * (gradient(gi%, l%).rdr * ur! + (1 + gradient(gi%, r%).rdl) * ul!))
       green! = gradient(gi%, l%).g + (gradient(gi%, r%).g - gradient(gi%, l%).g) * (u3! + 3 * (gradient(gi%, l%).gdr * ur! + (1 + gradient(gi%, r%).gdl) * ul!))
       blue! = gradient(gi%, l%).b + (gradient(gi%, r%).b - gradient(gi%, l%).b) * (u3! + 3 * (gradient(gi%, l%).bdr * ur! + (1 + gradient(gi%, r%).bdl) * ul!))
   Else
       red! = gradient(gi%, l%).r + (gradient(gi%, r%).r - gradient(gi%, l%).r) * u!
       green! = gradient(gi%, l%).g + (gradient(gi%, r%).g - gradient(gi%, l%).g) * u!
       blue! = gradient(gi%, l%).b + (gradient(gi%, r%).b - gradient(gi%, l%).b) * u!
   End If
End Sub

'> merger: Skipping unused SUB testGradient (gi%)

'####################################################################################################################
'# Ultra Fractal Opacity library v1.0 (routines)
'# By Zom-B
'#
'# Smooth Opacity algorithm from Ultra Fractal (www.ultrafractal.com)
'####################################################################################################################

Sub setNumOpacities (gi%)
   offset% = LBound(opacityPoints) - 1
   ReDim _Preserve opacitySmooth(gi% + offset%) As _Byte '_BIT <- bugged
   ReDim _Preserve opacityPoints(gi% + offset%) As Integer
   ReDim _Preserve opacity(gi% + offset%, 1) As OPACITYPOINT
End Sub

Sub addOpacityPoint (gi%, index!, o!)
   p% = opacityPoints(gi%)

   If UBound(opacity, 2) < p% Then
       ReDim _Preserve opacity(0 To UBound(opacity, 1), 0 To p%) As OPACITYPOINT
   End If

   opacity(gi%, p%).index = index!
   opacity(gi%, p%).o = o!
   opacityPoints(gi%) = p% + 1
End Sub

Sub setOpacitySmooth (gi%, s%)
   opacitySmooth(gi%) = s%

   If opacitySmooth(0) = 0 Then Exit Sub

   For i% = 1 To opacityPoints(gi%) - 1
       i1% = i% + 1
       If i1% = opacityPoints(gi%) Then i1% = 2

       dxl! = opacity(gi%, i%).index - opacity(gi%, i% - 1).index
       dxr! = opacity(gi%, i1%).index - opacity(gi%, i%).index
       If dxl! < 0 Then dxl! = dxl! + 1
       If dxr! < 0 Then dxr! = dxr! + 1

       d! = (opacity(gi%, i%).o - opacity(gi%, i% - 1).o) * dxr!
       If d! = 0 Then
           opacity(gi%, i%).odr = 0
           opacity(gi%, i%).odl = 0
       Else
           d! = (opacity(gi%, i1%).o - opacity(gi%, i%).o) * dxl! / d!
           If d! <= 0 Then
               opacity(gi%, i%).odr = 0
               opacity(gi%, i%).odl = 0
           Else
               opacity(gi%, i%).odr = 1 / (1 + d!)
               opacity(gi%, i%).odl = opacity(gi%, i%).odr - 1
           End If
       End If
   Next
End Sub

'####################################################################################################################

Sub getOpacity (gi%, index!, opacity!)
   If index! < 0 Then x! = 0 Else x! = index! - Int(index!)

   For l% = opacityPoints(gi%) - 2 To 1 Step -1
       If opacity(gi%, l%).index <= x! Then
           Exit For
       End If
   Next

   r% = l% + 1
   u! = (x! - opacity(gi%, l%).index) / (opacity(gi%, r%).index - opacity(gi%, l%).index)

   If opacitySmooth(gi%) Then
       u2! = u! * u!
       u3! = u2! * u!
       ur! = u3! - (u2! + u2!) + u!
       ul! = u2! - u3!

       opacity! = opacity(gi%, l%).o + (opacity(gi%, r%).o - opacity(gi%, l%).o) * (u3! + 3 * (opacity(gi%, l%).odr * ur! + (1 + opacity(gi%, r%).odl) * ul!))
   Else
       opacity! = opacity(gi%, l%).o + (opacity(gi%, r%).o - opacity(gi%, l%).o) * u!
   End If
End Sub

'> merger: Skipping unused SUB testOpacity (gi%)

'####################################################################################################################
'# Merge modes library v0.1 (routines)
'# By Zom-B
'####################################################################################################################

'> merger: Skipping unused SUB testMerge

'####################################################################################################################

'> merger: Skipping unused SUB mergeOverlay (br!, bg!, bb!, tr!, tg!, tb!)

Sub mergeHardLight (br!, bg!, bb!, tr!, tg!, tb!)
   If tr! <= 0.5 Then br! = br! * tr! * 2 Else br! = 1 - (1 - br!) * (1 - tr!) * 2
   If tg! <= 0.5 Then bg! = bg! * tg! * 2 Else bg! = 1 - (1 - bg!) * (1 - tg!) * 2
   If tb! <= 0.5 Then bb! = bb! * tb! * 2 Else bb! = 1 - (1 - bb!) * (1 - tb!) * 2
End Sub

Sub mergeSoftLight (br!, bg!, bb!, tr!, tg!, tb!)
   If tr! <= 0.5 Then br! = br! * (tr! + 0.5) Else br! = 1 - (1 - br!) * (1.5 - tr!)
   If tg! <= 0.5 Then bg! = bg! * (tg! + 0.5) Else bg! = 1 - (1 - bg!) * (1.5 - tg!)
   If tb! <= 0.5 Then bb! = bb! * (tb! + 0.5) Else bb! = 1 - (1 - bb!) * (1.5 - tb!)
End Sub

'> merger: Skipping unused SUB mergeColor (r!, g!, b!, r2!, g2!, b2!)

'> merger: Skipping unused SUB mergeHSLAddition (r!, g!, b!, r2!, g2!, b2!)

'####################################################################################################################

'> merger: Skipping unused SUB mergeHue (r!, g!, b!, r2!, g2!, b2!)

'> merger: Skipping unused SUB rgb2hsl (r!, g!, b!, chr!, smallest!, hue!, sat!, lum!)

'> merger: Skipping unused SUB hsl2rgb (hue!, sat!, lum!, r!, g!, b!)

'> merger: Skipping unused SUB hsl2rgb2 (hue!, chr!, smallest!, r!, g!, b!)


It would take plenty of study but the speed of processing is really nice! maybe some tips could be picked up from this example?

I wonder if all that vodeo mode crap can be replaced by setting up your screen with _DesktopWidth and height and using _NewImage to create the screen.