string math

General FreeBASIC programming questions.
srvaldez
Posts: 3373
Joined: Sep 25, 2005 21:54

string math

Post by srvaldez »

just for fun, to compare performance with my decfloat routines

Code: Select all

'BIGNUM.BAS v0.n
'Sep-Dec 1996 by Marc Kummel aka Treebeard.
'Contact mkummel@rain.org, http://www.rain.org/~mkummel/
'
' ** site no longer available, use the link below
' https://web.archive.org/web/20200220020034/http://www.rain.org/~mkummel/tbvault.html  
/'
  Conditions:
-------------

This program and source code are yours to use and modify as you will, but 
they are offered as freeware with no warranty whatsoever.  Give me credit, 
but do not distribute any changes under my name, or attribute such changes 
to me in any way.  You're on your own! 
'/
' adapted to FreeBasic by srvaldez

Const zero = "0", one = "1"
Const dp = ".", neg = "-", asc0 = 48
Const negative = -1, positive = 1 'returned by bComp()
Const maxlongdig = 8 'max digits in long var&

'useful shared stuff
Dim Shared as long digits
Dim shared as string errors

const null = ""
digits = 2785 'digits for division etc

's = |s|
'
SUB bAbs (byref s as string)
	IF LEFT(s, 1) = neg THEN s = MID(s, 2)
END SUB

'return true if s is negative
'
FUNCTION bIsNeg (byref s as string) as long
	bIsNeg = (LEFT(s, 1) = neg)
END FUNCTION

'return sign of number (-1 or +1)
'
FUNCTION bSign (byref s as string) as long
	IF bIsNeg(s) THEN bSign = negative ELSE bSign = positive
END FUNCTION


'return the largest of two integers
'
FUNCTION bMaxInt (byval n1 as long, byval n2 as long) as long
	IF n1 >= n2 THEN bMaxInt = n1 ELSE bMaxInt = n2
END FUNCTION

'Compare two numbers using fast string compares.  This can screw up since it
'uses string length, eg it reports "8"<"8." so watch out.  The practice in
'these routines is no leading or trailing 0s and no final "."  See bClean().
'
'Return 1 if s1 > s2
'       0 if s1 = s2
'      -1 if s1 < s2
'
FUNCTION bComp (byref s1 as string, byref s2 as string) as long
	'dim as string dp
	dim as long sign1, sign2, dp1, dp2, arg, s1flag, s2flag
	'kludge to fix 0<.1
	IF LEFT(s1, 1) = dp THEN s1 = zero + s1: s1flag = true
	IF LEFT(s2, 1) = dp THEN s2 = zero + s2: s2flag = true

	sign1 = (LEFT(s1, 1) = neg)
	sign2 = (LEFT(s2, 1) = neg)
	dp1 = INSTR(s1, dp): IF dp1 = 0 THEN dp1 = LEN(s1) + 1
	dp2 = INSTR(s2, dp): IF dp2 = 0 THEN dp2 = LEN(s2) + 1

	IF sign1 <> sign2 THEN
		IF sign1 THEN arg = -1 ELSE arg = 1
	ELSEIF s1 = s2 THEN
		arg = 0
	ELSEIF (dp1 < dp2) OR ((dp1 = dp2) AND (s1 < s2)) THEN
		arg = -1
	ELSE
		arg = 1
	END IF

	IF sign1 AND sign2 THEN arg = -arg
	IF s1flag THEN s1 = MID(s1, 2)
	IF s2flag THEN s2 = MID(s2, 2)
	bComp = arg

END FUNCTION

'return true if s1 > s2
'
FUNCTION bIsMore (byref s1 as string, byref s2 as string) as long
	bIsMore = (bComp(s1, s2) = 1)
END FUNCTION
'Strip leading 0s and final "." (but leave something)
'
SUB bStripZero (byref s as string)
	dim as long n
	n = 1
	DO WHILE MID(s, n, 1) = zero
		n = n + 1
	LOOP
	IF n > 1 THEN s = MID(s, n)
	IF RIGHT(s, 1) = dp THEN s = LEFT(s, LEN(s) - 1)
	IF LEN(s) = 0 THEN s = zero
END SUB

'Strip trailing 0s to "." (but leave something)
'
SUB bStripTail (byref s as string)
	dim as long n
	n = LEN(s)
	DO WHILE MID(s, n, 1) = zero
		n = n - 1
		IF n <= 1 THEN EXIT DO
	LOOP
	IF n THEN IF MID(s, n, 1) = dp THEN n = n - 1
	s = LEFT(s, n)
	IF LEN(s) = 0 THEN s = zero
END SUB

'Strip s$ to whole number and base 10 integer logarithm and sign.  Decimal
'point is implied after the first digit, and slog% counts places left or
'right.  bLogPut() reverses the process, and bLogDp() gives info on the
'decimals. Tricky, but it works and simplifies dividing and multipling.
'eg s$ -> s$ , slog%
'  660 -> 66 ,  2     (6.6 * 10^ 2)    (or 660,2 if zeroflag%=false)
'  6.6 -> 66 ,  0     (6.6 * 10^ 0)
' .066 -> 66 , -2     (6.6 * 10^-2)
'bDiv(), bMul(), and bSqr() use this to trim unnecessary zeros and to locate
'decimal point.  These set zeroflag% to trim trailing zeros, but bDivIntMod()
'must set it false in order to figure remainder of division.  A kludge.
'
SUB bLogGet (byref s as string, byref slog as long, byref sign as long, byref zeroflag as long)
	dim as long n, dpt
	IF LEFT(s, 1) = neg THEN s = MID(s, 2): sign = negative ELSE sign = positive
	bStripZero s
	dpt = INSTR(s, dp)
	SELECT CASE dpt
		CASE 0
			slog = LEN(s) - 1
		CASE 1
			n = dpt + 1
			DO WHILE MID(s, n, 1) = zero
				n = n + 1
			LOOP
			s = MID(s, n)
			slog = dpt - n
		CASE ELSE
			s = LEFT(s, dpt - 1) + MID(s, dpt + 1)
			slog = dpt - 2
	END SELECT

	'remove trailing 0's if zeroflag
	IF zeroflag THEN bStripTail s
  
END SUB

'Strip a number to "standard form" with no leading or trailing 0s and no
'final "."  All routines should return all arguments in this form.
'
SUB bClean (byref s as string)
	dim as long sign
	IF LEFT(s, 1) = neg THEN s = MID(s, 2): sign = true
	bStripZero s
	IF INSTR(s, dp) THEN bStripTail s
	IF sign AND s <> zero THEN s = neg + s
END SUB

'Restore a number from the integer and log figured in bLogGet(). s$ is taken
'as a number with the decimal after first digit, and decimal is moved slog%
'places left or right, adding 0s as required. Called by bDiv() and bMul().
'
SUB bLogPut (byref s as string, byref slog as long, byref sign as long)
	dim as long last
	last = LEN(s)
	IF LEN(s) = 0 OR s = zero THEN
		s = zero
	ELSEIF slog < 0 THEN
		s = dp + STRING(-slog - 1, zero) + s
	ELSEIF slog > last - 1 THEN
		s = s + STRING(slog - last + 1, zero) + dp
	ELSE
		s = LEFT(s, slog + 1) + dp + MID(s, slog + 2)
	END IF
	bClean s
	IF sign = negative THEN s = neg + s
END SUB

'shift decimal n% digits (minus=left), i.e multiply/divide by 10.
'
SUB bShift (byref s as string, byval n as long)
	dim as long slog, sign
	bLogGet s, slog, sign, 0 'false
	bLogPut s, slog + n, sign
END SUB

's = -s
'
SUB bNeg (byref s as string)
	IF LEFT(s, 1) = neg THEN s = MID(s, 2) ELSE s = neg + s
END SUB

'Take whole number and log from bLogGet() and return number of decimal
'places in the expanded number; OR take string and number of decimal points
'desired and return the log.  It works both ways.
'
FUNCTION bLogDp (byref s as string, byval logdp as long) as long
	bLogDp = LEN(s) - 1 - logdp
END FUNCTION

'out = s1 / s2 using fast long-integer algorithm. s2$ must be <= 8 digits.
's1$ and s2$ must be stripped first, no decimals.
'
SUB bDivLong (byref s1 as string, byref s2 as string, byref quotient as string, byref remain as string)
	dim as long dividend, remainder, divisor, dig, i
	quotient = null
	remainder = 0
	divisor = VAL(s2)

	FOR i = 1 TO digits
		dividend = remainder * 10 + VAL(MID(s1, i, 1))
		dig = dividend \ divisor
		quotient = quotient + CHR(asc0 + dig)
		remainder = dividend - dig * divisor
	NEXT i

	IF LEN(quotient) = 0 THEN quotient = zero
	remain = LTRIM(STR(remainder))

END SUB

'out = s1 / s2 using character algorithm, digit by digit, slow but honest.
's1$ and s2$ must be stripped first, no decimals.
'
SUB bDivChar (byref s1 as string, byref s2 as string, byref quotient as string, byref remainder as string)
	dim as long last1, last2, ldvd, dig, borrow, i, j, n, lrem
	dim as string dvd
	last1 = LEN(s1)     'length of the dividend
	last2 = LEN(s2)     'length of the divisor
	quotient = null
	remainder = null

	FOR i = 1 TO digits
		'get next digit of dividend or zero$ if past end
		IF i <= last1 THEN
			dvd = remainder + MID(s1, i, 1)
		ELSE
			dvd = remainder + zero
		END IF

		'if dividend < divisor then digit%=0 else have to calculate it.
		'do fast compare using string operations. see bComp%()
		bStripZero dvd
		ldvd = LEN(dvd)
		IF (ldvd < last2) OR ((ldvd = last2) AND (dvd < s2)) THEN
			'divisor is bigger, so digit is 0, easy!
			dig = 0
			remainder = dvd

		ELSE
			'dividend is bigger, but no more than 9 times bigger.
			'subtract divisor until we get remainder less than divisor.
			'time hog, average is 5 tries through j% loop.  There's a better way.
			FOR dig = 1 TO 9
				remainder = null
				borrow = 0
				FOR j = 0 TO ldvd - 1
					n = last2 - j
					IF n < 1 THEN n = 0 ELSE n = VAL(MID(s2, n, 1))
					n = VAL(MID(dvd, ldvd - j, 1)) - n - borrow
					IF n >= 0 THEN borrow = 0 ELSE borrow = 1: n = n + 10
					remainder = CHR(asc0 + n) + remainder
				NEXT j

				'if remainder < divisor then exit
				bStripZero remainder
				lrem = LEN(remainder)
				IF (lrem < last2) OR ((lrem = last2) AND (remainder < s2)) THEN EXIT FOR

				dvd = remainder
				ldvd = LEN(dvd)
			NEXT dig

		END IF
		quotient = quotient + CHR(asc0 + dig)
	NEXT i

END SUB

'out = s1 / s2
'
SUB bDiv (byref s1 as string, byref s2 as string, byref outs as string)
	dim as string t
	dim as long slog1, slog2, sign1, sign2, outlog, outsign, olddigits
	'strip divisor
	t = s2
	bLogGet t, slog2, sign2, -1 'true

	'divide by zero?
	IF t = zero THEN
		outs = errors

		'do powers of 10 with shifts
	ELSEIF t = one THEN
		outs = s1
		sign1 = bSign(outs)
		IF sign1 = negative THEN bAbs outs
		bShift outs, -slog2
		IF sign1 <> sign2 THEN bNeg outs

		'the hard way
	ELSE
		'strip all
		s2 = t: t = null
		bLogGet s1, slog1, sign1, -1 'true

		'figure decimal point and sign of answer
		outlog = slog1 + bLogDp(s2, slog2)
		IF sign1 <> sign2 THEN outsign = negative ELSE outsign = positive

		'bump digits past leading zeros and always show whole quotient
		olddigits = digits
		digits = digits + LEN(s2)
		IF digits < outlog + 1 THEN digits = outlog + 1

		'do it, ignore remainder
		IF LEN(s2) <= maxlongdig THEN bDivLong s1, s2, outs, t ELSE bDivChar s1, s2, outs, t

		'clean up
		bLogPut outs, outlog, outsign
		bLogPut s1, slog1, sign1
		bLogPut s2, slog2, sign2
		digits = olddigits
	END IF

END SUB

'Trim leading spaces, add decimal points, eliminate signs.
'Returns last%=length of string, dpt%=decimal place, sign%=-1 or 1.
'Called only by bAdd() and bSub() which needs a final decimal point.
'
SUB bStripDp (byref s as string, byref last as long, byref dpt as long, byref sign as long)
	IF LEFT(s, 1) = neg THEN s = MID(s, 2): sign = negative ELSE sign = positive
	bStripZero s
	IF INSTR(s, dp) = 0 THEN s = s + dp
	IF s = dp THEN s = "0."
	dpt = INSTR(s, dp)
	last = LEN(s)
END SUB

declare SUB bAdd (byref s1 as string, byref s2 as string, byref outs as string)
'out = s1 - s2
'
SUB bSub (byref s1 as string, byref s2 as string, byref outs as string)
	dim as long last, last1, last2, sign1, sign2, dpt, dp1, dp2
	dim as long d1, d2, borrow, swapflag, i, n
	'strip the numbers
	bStripDp s1, last1, dp1, sign1
	bStripDp s2, last2, dp2, sign2

	'treat different signs as addition
	IF sign1 = negative AND sign2 = positive THEN
		bNeg s1
		bNeg s2
		bAdd s1, s2, outs
		bNeg s2
		EXIT SUB
	ELSEIF sign1 = positive AND sign2 = negative THEN
		bAdd s1, s2, outs
		bNeg s2
		EXIT SUB
	END IF

	'align the decimal points and digit pointers
	last = bMaxInt(last1 - dp1, last2 - dp2)
	d1 = last + dp1
	d2 = last + dp2
	dpt = bMaxInt(dp1, dp2)
	last = dpt + last
	outs = SPACE(last)
	borrow = 0

	'always subtract smaller from bigger to avoid complements
	IF bIsMore(s2, s1) THEN SWAP s2, s1: SWAP d2, d1: swapflag = true

	'do the subtraction right to left
	FOR i = last TO 1 STEP -1
		IF i <> dpt THEN
			IF d1 > 0 THEN n = VAL(MID(s1, d1, 1)) ELSE n = 0
			IF d2 > 0 THEN n = n - VAL(MID(s2, d2, 1))
			n = n - borrow
			IF n >= 0 THEN borrow = 0 ELSE borrow = 1: n = n + 10
			MID(outs, i, 1) = CHR(asc0 + n)
		ELSE
			MID(outs, i, 1) = dp
		END IF
		d1 = d1 - 1
		d2 = d2 - 1
	NEXT i

	'clean up
	IF sign1 = negative THEN s1 = neg + s1: s2 = neg + s2
	IF swapflag THEN SWAP s2, s1: sign1 = -sign1
	IF sign1 = negative THEN outs = neg + outs
	bClean s1
	bClean s2
	bClean outs

END SUB

'out = s1 + s2
'
SUB bAdd (byref s1 as string, byref s2 as string, byref outs as string)
	dim as long last1, last2, sign1, sign2, dp1, dp2
	dim as long last, i, d1, d2, dpt, carry, n
	'strip the numbers
	bStripDp s1, last1, dp1, sign1
	bStripDp s2, last2, dp2, sign2

	'treat different signs as subtraction and exit
	IF sign1 = negative AND sign2 = positive THEN
		bSub s2, s1, outs
		bNeg s1
		EXIT SUB
	ELSEIF sign1 = positive AND sign2 = negative THEN
		bSub s1, s2, outs
		bNeg s2
		EXIT SUB
	END IF

	'align the decimal points and digit pointers
	last = bMaxInt(last1 - dp1, last2 - dp2)
	d1 = last + dp1
	d2 = last + dp2
	dpt = bMaxInt(dp1, dp2)
	last = dpt + last
	outs = SPACE(last)
	carry = 0

	'do the addition right to left
	FOR i = last TO 1 STEP -1
		IF i <> dpt THEN
			n = carry
			IF d1 > 0 THEN n = n + VAL(MID(s1, d1, 1))
			IF d2 > 0 THEN n = n + VAL(MID(s2, d2, 1))
			carry = n \ 10
			MID(outs, i, 1) = CHR(asc0 + (n MOD 10))
		ELSE
			MID(outs, i, 1) = dp
		END IF
		d1 = d1 - 1
		d2 = d2 - 1
	NEXT i
	IF carry THEN outs = one + outs

	'clean up
	IF sign1 = negative THEN s1 = neg + s1: s2 = neg + s2: outs = neg + outs
	bClean s1
	bClean s2
	bClean outs
END SUB

'out = s1 * s2 using character algorithm, slow but honest.  Whole numbers
'only.  Inner loop is optimized and hard to understand, but it works.
'
SUB bMulChar (byref s1 as string, byref s2 as string, byref outs as string)
	dim as long last, last1, last2, i, j, k, sj, ej, product
	last1 = LEN(s1)
	last2 = LEN(s2)
	last = last1 + last2
	outs = SPACE(last)
	product = 0
	FOR i = 0 TO last - 1
		K = last1 - i
		sj = 1 - K: IF sj < 0 THEN sj = 0
		ej = last1 - K: IF ej > last2 - 1 THEN ej = last2 - 1
		FOR j = sj TO ej
			product = product + VAL(MID(s1, K + j, 1)) * VAL(MID(s2, last2 - j, 1))
		NEXT j
		MID(outs, last - i, 1) = CHR(asc0 + CINT(product MOD 10))
		product = product \ 10
	NEXT i
	IF product THEN outs = LTRIM(STR(product)) + outs
END SUB

'out = s1 * s2 using fast long-integer algorithm. s2$ must be <= 8 digits.
's1$ and s2$ must be stripped first, whole numbers only.
'
SUB bMulLong (byref s1 as string, byref s2 as string, byref outs as string)
	dim as long last1, s2L, i, product
	last1 = LEN(s1)
	s2L = VAL(s2)
	outs = SPACE(last1)
	FOR i = last1 TO 1 STEP -1
		product = product + VAL(MID(s1, i, 1)) * s2L
		MID(outs, i, 1) = CHR(asc0 + CINT(product MOD 10))
		product = product \ 10
	NEXT i
	IF product THEN outs = LTRIM(STR(product)) + outs
END SUB

'out = s1 * s2
'
Sub bMul (byref s1 as string, byref s2 as string, byref outs as string)
	dim as long slog1, slog2, sign1, sign2, outdp, outsign
	dim as long outlog, swapflag
	dim as string t
    'strip multiplier
    t = s2
    bLogGet t, slog2, sign2, true

    'times 0
    If t = zero Then
        outs = zero

        'do powers of 10 with shifts
    ElseIf t = one Then
        outs = s1
		sign1 = bSign(outs)
		IF sign1 = negative THEN bAbs outs
		bShift outs, slog2
		IF sign1 <> sign2 THEN bNeg outs

	'the hard way
	ELSE
		'strip all
		s2 = t: t = null
		bLogGet s1, slog1, sign1, true

		'figure decimal point and sign of answer
		outdp = bLogDp(s1, slog1) + bLogDp(s2, slog2)
		IF sign1 <> sign2 THEN outsign = negative ELSE outsign = positive

		'always multiply by the shorter number
		IF LEN(s2) > LEN(s1) THEN SWAP s1, s2: swapflag = true

		'do it
		IF LEN(s2) <= maxlongdig THEN bMulLong s1, s2, outs ELSE bMulChar s1, s2, outs

		'clean up
		outlog = bLogDp(outs, outdp)
		bLogPut outs, outlog, outsign
		IF swapflag THEN SWAP s1, s2
		bLogPut s1, slog1, sign1
		bLogPut s2, slog2, sign2

  END IF

END SUB

'========================================================================
dim as double t
dim as string n, m, r
dim as long i, j
t = Timer
n = "1"
m = "999999999999999999999998999999999999999999999999"

r = ""
bDiv(n, m, r)
j = InStr(r, ".")
r = Mid(r, j + 1)
j = 1
For i = 1 To 58
    Print Mid(r, j, 24), Mid(r, j + 24, 24)
    j = j + 48
Next

t = Timer - t
'Print r
Print t
'========================================================================
time .174 seconds
decfloat time .0194 seconds
D.J.Peters
Posts: 8586
Joined: May 28, 2005 3:28
Contact:

Re: string math

Post by D.J.Peters »

Only a hint, put two decimal 0-9 numbers (4bit) in one byte and use the fast BCD arithmetic CPU instructions !
https://en.wikipedia.org/wiki/Intel_BCD_opcode

Joshy
jj2007
Posts: 2326
Joined: Oct 23, 2016 15:28
Location: Roma, Italia
Contact:

Re: string math

Post by jj2007 »

Joshy, the BCD instructions are notoriously old and slow, but they might still be faster than a high level language, of course ;-)
srvaldez
Posts: 3373
Joined: Sep 25, 2005 21:54

Re: string math

Post by srvaldez »

the decimal-adjust instructions are not available in 64-bit, besides I think I had enough of this math stuff for a while
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: string math

Post by dodicat »

Here is high school long division.

Code: Select all


Function divide(n1 As String,n2 As String,decimal_places As integer,dpflag As String="s") As String
          Dim As String number=n1,divisor=n2
          dpflag=lcase(dpflag)
          'For MOD
          dim as integer modstop
          if dpflag="mod" then 
              if len(n1)<len(n2) then return n1
              if len(n1)=len(n2) then
                  if n1<n2 then return n1
                  end if
              modstop=len(n1)-len(n2)+1
              end if
          if dpflag<>"mod" then
     If dpflag<>"s"  Then dpflag="raw" 
     end if
        Dim runcount As integer
        '_______  LOOK UP TABLES ______________
        Dim Qmod(0 To 19) As Ubyte
        Dim bool(0 To 19) As Ubyte
        For z As Integer=0 To 19
    Qmod(z)=(z Mod 10+48)
    bool(z)=(-(10>z))
Next z
Dim answer As String   'THE ANSWER STRING  

'_______ SET THE DECIMAL WHERE IT SHOULD BE AT _______
Dim As String part1,part2
#macro set(decimal)
#macro insert(s,char,position)
If position > 0 And position <=Len(s) Then
part1=Mid$(s,1,position-1)
part2=Mid$(s,position)
s=part1+char+part2
End if
#endmacro
insert(answer,".",decpos)
  answer=thepoint+zeros+answer
If dpflag="raw" Then
    answer=Mid(answer,1,decimal_places)
    End if
#endmacro
'______________________________________________
'__________ SPLIT A STRING ABOUT A CHARACTRR __________
Dim As String var1,var2
    Dim pst As integer
      #macro split(stri,char,var1,var2)
    pst=Instr(stri,char)
    var1="":var2=""
    If pst<>0 Then
    var1=Rtrim(Mid(stri,1,pst),".")
    var2=Ltrim(Mid(stri,pst),".")
Else
    var1=stri
    End if
    #endmacro
    
       #macro Removepoint(s)
       split(s,".",var1,var2)
#endmacro
'__________ GET THE SIGN AND CLEAR THE -ve __________________
Dim sign As String
          If Left(number,1)="-" Xor Left (divisor,1)="-" Then sign="-"
            If Left(number,1)="-" Then  number=Ltrim(number,"-")
            If Left (divisor,1)="-" Then divisor=Ltrim(divisor,"-")
              
'DETERMINE THE DECIMAL POSITION BEFORE THE DIVISION
Dim As integer lennint,lenddec,lend,lenn,difflen
split(number,".",var1,var2)
lennint=Len(var1)
split(divisor,".",var1,var2)
lenddec=Len(var2)

If Instr(number,".") Then 
    Removepoint(number)
    number=var1+var2
    End if
If Instr(divisor,".") Then 
    Removepoint(divisor)
    divisor=var1+var2
    End if
Dim As integer numzeros
numzeros=Len(number)
number=Ltrim(number,"0"):divisor=Ltrim (divisor,"0")
numzeros=numzeros-Len(number)
lend=Len(divisor):lenn=Len(number)
If lend>lenn Then difflen=lend-lenn
Dim decpos As integer=lenddec+lennint-lend+2-numzeros 'THE POSITION INDICATOR
Dim _sgn As Byte=-Sgn(decpos)
If _sgn=0 Then _sgn=1
Dim As String thepoint=String(_sgn,".") 'DECIMAL AT START (IF)
Dim As String zeros=String(-decpos+1,"0")'ZEROS AT START (IF) e.g. .0009
if dpflag<>"mod" then
If Len(zeros) =0 Then dpflag="s"
end if
Dim As integer runlength
If Len(zeros) Then 
     runlength=decimal_places
     answer=String(Len(zeros)+runlength+10,"0")
    If dpflag="raw" Then 
        runlength=1
        answer=String(Len(zeros)+runlength+10,"0")
        If decimal_places>Len(zeros) Then
            runlength=runlength+(decimal_places-Len(zeros))
            answer=String(Len(zeros)+runlength+10,"0")
            End If
            End If

Else
decimal_places=decimal_places+decpos
runlength=decimal_places
answer=String(Len(zeros)+runlength+10,"0")
End if
'___________DECIMAL POSITION DETERMINED  _____________

'SET UP THE VARIABLES AND START UP CONDITIONS
number=number+String(difflen+decimal_places,"0")
        Dim count As integer
        Dim temp As String
        Dim copytemp As String
        Dim topstring As String
        Dim copytopstring As String
        Dim As integer lenf,lens
        Dim As Ubyte takeaway,subtractcarry
        Dim As integer n3,diff
       If Ltrim(divisor,"0")="" Then Return "Error :division by zero"   
        lens=Len(divisor)
         topstring=Left(number,lend)
         copytopstring=topstring
        Do
            count=0
        Do
            count=count+1
            copytemp=temp
    
            Do
'___________________ QUICK SUBTRACTION loop _________________              
            
lenf=Len(topstring)
If  lens<lenf=0 Then 'not
If Lens>lenf Then
temp= "done"
Exit Do
End if
If divisor>topstring Then 
temp= "done"
Exit Do
End if
End if

  diff=lenf-lens
        temp=topstring
        subtractcarry=0
        
        For n3=lenf-1 To diff Step -1
            takeaway= topstring[n3]-divisor[n3-diff]+10-subtractcarry
            temp[n3]=Qmod(takeaway)
            subtractcarry=bool(takeaway)
        Next n3 
        If subtractcarry=0 Then Exit Do
         If n3=-1 Then Exit Do
        For n3=n3 To 0 Step -1 
            takeaway= topstring[n3]-38-subtractcarry
             temp[n3]=Qmod(takeaway)
            subtractcarry=bool(takeaway)
            if subtractcarry=0 then exit do
            Next n3
        Exit Do
        
        Loop 'single run
        temp=Ltrim(temp,"0")
        If temp="" Then temp= "0"
            topstring=temp
        Loop Until temp="done"
     ' INDIVIDUAL CHARACTERS CARVED OFF ________________       
        runcount=runcount+1
       If count=1 Then
           topstring=copytopstring+Mid(number,lend+runcount,1)
           Else
       topstring=copytemp+Mid(number,lend+runcount,1)
   End If
       copytopstring=topstring
       topstring=Ltrim(topstring,"0")
       if dpflag="mod" then
       if runcount=modstop then 
           if topstring="" then return "0"
           return mid(topstring,1,len(topstring)-1)
           end if
       end if
       answer[runcount-1]=count+47
       If topstring="" And runcount>Len(n1)+1 Then
           Exit Do
           End if
   Loop Until runcount=runlength+1
   
   ' END OF RUN TO REQUIRED DECIMAL PLACES
   set(decimal) 'PUT IN THE DECIMAL POINT
  'THERE IS ALWAYS A DECIMAL POINT SOMEWHERE IN THE ANSWER
  'NOW GET RID OF IT IF IT IS REDUNDANT
       answer=Rtrim(answer,"0")
       answer=Rtrim(answer,".")
       answer=Ltrim(answer,"0")
       If answer="" Then Return "0"
   Return sign+answer
End Function

Function GetFibonacci(n As Long) As String
    Static As Ubyte Qmod(0 To 19)={48,49,50,51,52,53,54,55,56,57,48,49,50,51,52,53,54,55,56,57}
    Static As Ubyte Qbool(0 To 19)={0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1}
    Dim As Ubyte addup,addcarry
    Dim As Long diff,n_,LL
    If n=0 Then Return "0"
    If n=1 Or n=2 Then Return "1"
    If n=3 Then
        Return "2"
    Else
        Dim As String sl="1",l="2",term
        For x As Long= 1 To n-3
            LL=Len(l)
            diff=-(LL<>Len(sl))
            term="0"+l
            addcarry=0
            For n_=LL-1 To diff Step -1
                addup=sl[n_-diff]+l[n_]-96
                term[n_+1]=Qmod(addup+addcarry)
                addcarry=Qbool(addup+addcarry)
            Next n_
            If addcarry=0 Then  term=Ltrim(term,"0"):Goto skip
            If n_=-1      Then term[0]=addcarry+48  :Goto skip
            For n_=n_ To 0 Step -1
                addup=l[n_]-48
                term[n_+1]=Qmod(addup+addcarry)
                addcarry=Qbool(addup+addcarry)
                If addcarry=0 Then Exit For
            Next n_
            term[0]=addcarry+48
            If addcarry=0 Then term=Ltrim(term,"0")
            skip:
            sl=l
            l=term
        Next x
        Function =term
    End If
End Function
'200 185000
'150 100000
'100=40000
'50=10000
'20=2000
'10=400
'5=80
'3=30
dim as double t=timer
dim as integer n=100

dim as string f=string(n-1,"9")+"8"+string(n,"9")
dim as integer places= _ 
 8.039594990255599 _
+ 1.338375048882362*n _
+0.8882107460008227*n^2 _
+0.4071177958035269*n^3 _
 -1.335359078473998e-002 *n^4 _
+1.708530082532032e-004 *n^5 _
 -9.296417579629874e-007 *n^6 _
+ 1.8032897354154e-009 *n^7

f= divide("1",f,places)
dim as integer counter
'get the counter
for i as integer=1 to len(f) step n'91
    counter+=1
next i
width 120,counter+len(f)/95+50
print "START, the Fibonacci numbers are held in the division"
print f
print
print "The extracted Fibonnaci numbers"
counter=0
dim as string s
for i as integer=1 to len(f) step n'91
  if counter<>401 then  s+="FIB "+str(counter)+" = "+ltrim(mid(f,i+1,n),"0") +chr(10)
    counter+=1
next i
print s
    print
    print "Check the last FIB number independently"
    print getfibonacci(400)
    print timer-t
    'print counter
    'print places
sleep

 
srvaldez
Posts: 3373
Joined: Sep 25, 2005 21:54

Re: string math

Post by srvaldez »

hi dodicat
that's very fast!
your code can get up-to fibonacci(480)
SamL
Posts: 58
Joined: Nov 23, 2019 17:30
Location: Minnesota

Re: string math

Post by SamL »

Nice
srvaldez
Posts: 3373
Joined: Sep 25, 2005 21:54

Re: string math

Post by srvaldez »

Hi dodicat
I would like to contract you to implement multiple precision division for my decfloat routines, I know you are infinitely more capable than me, if you are interested then you can send me a pm
at the moment I do division by reciprocal which is kind of slow
just in case you are interested then you may use the updated decfloat routines from viewtopic.php?p=283661#p283661
here's what I have in mind

Code: Select all

Function fpdiv(Byref x As decfloat, Byref y As decfloat) As decfloat
    Dim As decfloat fac1, fac2, one, fac3
    Dim As Integer i, count, er, is_power_of_ten
    Dim As Short sign
    Dim As Longint n, d, q, r
    
    fac1=x
    fac2=y

    If fac2.exponent=0 Then ' if fac2 = 0, return
        ' a divide-by-zero error and
        ' bail out.
        For i=0 To NUM_DWORDS
            fac1.mantissa(i)=999999999
        Next
        fac1.exponent=99999+BIAS+1
        er=DIVZ_ERR
        Return fac1
    Elseif fac1.exponent=0 Then 'fact1=0, just return
        er=0
        Return fac1
    Else
        one.mantissa(0)=100000000 'one = 1.0
        one.exponent=BIAS+1
        one.sign=0
        For i=1 To NUM_DWORDS
            one.mantissa(i)=0
        Next
        'check to see if fac2 is a power of ten
        is_power_of_ten=0
        If fac2.mantissa(0)=100000000 Then
            is_power_of_ten=1
            For i=1 To NUM_DWORDS
                If fac2.mantissa(i)<>0 Then
                    is_power_of_ten=0
                    Exit For
                End If
            Next
        End If
        'if fac2 is a power of ten then all we need to do is to adjust the sign and exponent and we are finished
        If is_power_of_ten=1 Then
            fac1.sign=fac1.sign Xor fac2.sign
            fac1.exponent=fac1.exponent-fac2.exponent+BIAS+1
            Return fac1
        End If

        ' here the actual division would take place
        ' let's assume that the result is in fac3
        
    End If

    NORM_FAC1(fac3) 

    'calculate the exponent of the quotient.
    fac3.exponent=x.exponent-y.exponent+BIAS+1
    'calculate the sign of quotient.
    fac3.sign=x.sign Xor y.sign

    Return fac3
End Function
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: string math

Post by dodicat »

I did the division exactly as long division at school.
All that is needed is add and subtract (strings in my case).
I got a pen and paper and created a long division and followed the steps, moving the decimals along, getting each digit to suit (adding up n times until
n marked the value just short ( n is just one digit 0 to 9), then subtracted then bring the next digit of the numerator down and repeated.
This way you are gathering up all the n digits which form the answer.
I see you have multiply by an integer, this would save the adding, at school you would mentally multiply the divisor by a digit to see which one (1 to 9)
got you nearest so you could subtract for the next run.
Maybe I should practise with doubles, I think I might have to use dec2str with your method and just str with doubles, because the numerator has to be split to bring down a digit into play.
Richard and Albert in their big number methods also used reciprocals.
It is a bit of a pain in the neck simulating school long division, just a pedantic plod through the process.
srvaldez
Posts: 3373
Joined: Sep 25, 2005 21:54

Re: string math

Post by srvaldez »

thank you dodicat, I think I will attack this next week, right now I am busy with other things
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: string math

Post by dodicat »

Hi srvaldez.
I have reworked the long division, and used functions to do the tasks.
The long division mechanics are easy enough, getting a correct decimal place is fiddly.
It is not perfect, there will be anomalies if tested to destruction, but maybe it will help your multiple precision division in some way.

Code: Select all


Sub Remove(Text As String,Char As String)
      Var index = 0,asci=Asc(char)
      For i As Integer = 0 To Len(Text) - 1
            If Text[i] <> ASCi Then Text[index] = Text[i] : index =index+ 1
      Next : Text = Left(Text,index)
End Sub

Sub insertdecimal(s As String,position as long,char as string=".")
   Dim As String part1,part2
      If position > 0 And position <=Len(s) Then
            part1=Mid(s,1,position-1) 
            part2=Mid(s,position)
            s=part1+char+part2
      End If
End Sub

Function smult(n As String,d As long) As String 'multiply n by a single digit
      Dim As String ans=String(Len(n)+1,"*")
      Dim As Ubyte carry,main,temp
      For z As Integer=Len(n)-1 To 0 Step -1
            temp=(d)*(n[z]-48)+carry
            main=temp Mod 10 +48
            carry=temp\10
            ans[z+1]=main
      Next z
      ans[0]=carry+48
      Return Iif(Len(Ltrim(ans,"0")),Ltrim(ans,"0"),"0")
End Function

Function minus(num1 As String,num2 As String) As String 'subtract two strings
      Static As Ubyte sm(0 To 19)={48,49,50,51,52,53,54,55,56,57,48,49,50,51,52,53,54,55,56,57}
      Static As Ubyte sb(0 To 19)={1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0}       
      Dim As Integer bigger,swapflag           
      Dim sign As String * 1
      Var lenf=Len(NUM1)
      Var lens=Len(NUM2)
      #macro finishup()
      answer=Ltrim(answer,"0")
      If answer="" Then Return "0"
      If swapflag=1 Then Swap NUM1,NUM2
      Return sign+answer
      #endmacro
      #macro compare()
      If Lens>lenf Then bigger= -1:Goto fin
      If Lens<lenf Then bigger =0:Goto fin
      If NUM2>NUM1 Then 
            bigger=-1
      Else
            bigger= 0
      End If
      fin:
      #endmacro
      
      compare()
      If bigger Then 
            sign="-"
            Swap NUM2,NUM1
            Swap lens,lenf
            swapflag=1
      End If
      Var diff=lenf-lens
      Dim As String answer=NUM1
      Dim As Integer n
      Dim As Ubyte takeaway,subtractcarry
      subtractcarry=0
      For n=lenf-1 To diff Step -1 
            takeaway= num1[n]-num2[n-diff]+10-subtractcarry
            answer[n]=sm(takeaway)
            subtractcarry=sb(takeaway)
      Next n 
      
      If subtractcarry=0 Then:finishup():End If
      If n=-1 Then:finishup():End If
      For n=n To 0 Step -1 
            takeaway= num1[n]-38-subtractcarry 
            answer[n]=sm(takeaway)
            subtractcarry=sb(takeaway)
            If subtractcarry=0 Then Exit For
      Next n
      finishup()
End Function      

Function isbigger(n1 As String,n2 As String) As Long 'is n1>n2
      If Len(n1)>Len(n2) Then Return -1
      If Len(n1)<Len(n2) Then Return 0
      If n1>n2 Then 
            Return -1
      Else
            Return 0
      End If
End Function

Function decmove(s As String) As Long
      If Instr(s,".")=0 Then Return 0
      Return Len(s)-Instr(s,".")
End Function

Sub adjustnumerator(d As String,dm As Long)
      While decmove(d)<= dm
            d+="0"
      Wend
      Var p=Instr(d,".")+dm
      remove(d,".")
      insertdecimal(d,p)
End Sub

Function divide(n As String,d As String,places As Long) As String
    dim as string sign
      if instr(n,"-") xor instr(d,"-") then sign="-"
      n=ltrim(n,"0"):d=ltrim(d,"0")
      remove(n,"-"):remove(d,"-")
      if instr(n,".")=0 then n+="."
      if instr(d,".")=0 then d+="."
      Static As String s(0 To 10)'to hold digit multiples
      Dim As String ans
      Dim As Long k,flag
      Var dm=decmove(d)
      adjustnumerator(n,dm)
      n+=String(places-decmove(n),"0")
      remove(d,".")
      d=ltrim(d,"0")
      n=ltrim(n,"0")
      Dim As String top=Mid(n,1,1)
      'if instr(top,".") then top=rtrim(top,"0"):remove(top,".")
      remove(top,".")
      For m As Long=1 To Len(n)-1

            If n[m]=Asc(".") Then flag=1:goto nxt
           
            For k =0 To 10 'range up the digits until the product > the numerator part
                  If k=10 Then
                        s(k)=d+"0"
                  Else
                        s(k)=smult(d,k)
                  End If
                  If isbigger(s(k),top) Then ans+=Str(k-1):Exit For
            Next k
            top=minus(top,s(k-1))
            
            top+=Chr(n[m])'bring down the next digit
            if flag=1 then ans+=".":flag=0
            nxt:
      Next m
      
      ans=iif(instr(ans,"."),Ltrim(ans,"0"),ans)
      if len(ans)=0 then ans="0"
      
      Return sign+ ans
End Function


width 150,520

dim as double n,d
For m As Long=1 To 500
    do
      n=rnd-rnd
      d=rnd-rnd
      loop until instr(str(n),"e")=0 and instr(str(d),"e")=0
      Var dd=divide(str(n),str(d),30)
      Print m;" "; str(n);"/";str(d);tab(50);n/d;tab(80);dd;tab(120);cbool(Csng(n/d)=Csng(Val(dd)))
     
Next m
print

var t=timer
print divide("1",str(t),500)
print 1/t,"done"
sleep


 
The 500 tests run faster with -O(n) optimisations of course.
Last edited by dodicat on Aug 25, 2021 22:02, edited 2 times in total.
srvaldez
Posts: 3373
Joined: Sep 25, 2005 21:54

Re: string math

Post by srvaldez »

thank you dodicat
I translated Dr. David M. Smith's Mathematica code to FB http://dmsmith.lmu.build/MComp1996.pdf and it works ok but I think that it could be made faster if the base was changed from 10000 to 100000000
the advantage of using floating point is that the decimal place is easily determined by simply subtracting the exponents
my problem is that the long multiplication is not very fast, even with my tweaks
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: string math

Post by dodicat »

Albert had the fastest multiplier, he used base 1e7.
We tested base 1e4 through base 1e9 but 1e7 was the fastest.
He had it in one intact function.
Pity he wasn't here to show the low down.
It was a grid type method at the core.
Here is my expanded high school divider tested with actual bigger than ulongint numbers.
(Fibonacci by division again)

Code: Select all


Sub Remove(Text As String,Char As String)
      Var index = 0,asci=Asc(char)
      For i As Integer = 0 To Len(Text) - 1
            If Text[i] <> ASCi Then Text[index] = Text[i] : index =index+ 1
      Next : Text = Left(Text,index)
End Sub

Sub insertdecimal(s As String,position as long,char as string=".")
   Dim As String part1,part2
      If position > 0 And position <=Len(s) Then
            part1=Mid(s,1,position-1) 
            part2=Mid(s,position)
            s=part1+char+part2
      End If
End Sub

Function smult(n As String,d As long) As String 'multiply n by a single digit
      Dim As String ans=String(Len(n)+1,"*")
      Dim As Ubyte carry,main,temp
      For z As Integer=Len(n)-1 To 0 Step -1
            temp=(d)*(n[z]-48)+carry
            main=temp Mod 10 +48
            carry=temp\10
            ans[z+1]=main
      Next z
      ans[0]=carry+48
      Return Iif(Len(Ltrim(ans,"0")),Ltrim(ans,"0"),"0")
End Function

Function minus(num1 As String,num2 As String) As String 'subtract two strings
      Static As Ubyte sm(0 To 19)={48,49,50,51,52,53,54,55,56,57,48,49,50,51,52,53,54,55,56,57}
      Static As Ubyte sb(0 To 19)={1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0}       
      Dim As Integer bigger,swapflag           
      Dim sign As String * 1
      Var lenf=Len(NUM1)
      Var lens=Len(NUM2)
      #macro finishup()
      answer=Ltrim(answer,"0")
      If answer="" Then Return "0"
      If swapflag=1 Then Swap NUM1,NUM2
      Return sign+answer
      #endmacro
      #macro compare()
      If Lens>lenf Then bigger= -1:Goto fin
      If Lens<lenf Then bigger =0:Goto fin
      If NUM2>NUM1 Then 
            bigger=-1
      Else
            bigger= 0
      End If
      fin:
      #endmacro
      
      compare()
      If bigger Then 
            sign="-"
            Swap NUM2,NUM1
            Swap lens,lenf
            swapflag=1
      End If
      Var diff=lenf-lens
      Dim As String answer=NUM1
      Dim As Integer n
      Dim As Ubyte takeaway,subtractcarry
      subtractcarry=0
      For n=lenf-1 To diff Step -1 
            takeaway= num1[n]-num2[n-diff]+10-subtractcarry
            answer[n]=sm(takeaway)
            subtractcarry=sb(takeaway)
      Next n 
      
      If subtractcarry=0 Then:finishup():End If
      If n=-1 Then:finishup():End If
      For n=n To 0 Step -1 
            takeaway= num1[n]-38-subtractcarry 
            answer[n]=sm(takeaway)
            subtractcarry=sb(takeaway)
            If subtractcarry=0 Then Exit For
      Next n
      finishup()
End Function      

Function isbigger(n1 As String,n2 As String) As Long 'is n1>n2
      If Len(n1)>Len(n2) Then Return -1
      If Len(n1)<Len(n2) Then Return 0
      If n1>n2 Then 
            Return -1
      Else
            Return 0
      End If
End Function

Function decmove(s As String) As Long
      If Instr(s,".")=0 Then Return 0
      Return Len(s)-Instr(s,".")
End Function

Sub adjustnumerator(d As String,dm As Long)
      While decmove(d)<= dm
            d+="0"
      Wend
      Var p=Instr(d,".")+dm
      remove(d,".")
      insertdecimal(d,p)
End Sub

Function divide(n As String,d As String,places As Long) As String
    dim as string sign
      if instr(n,"-") xor instr(d,"-") then sign="-"
      n=ltrim(n,"0"):d=ltrim(d,"0")
      remove(n,"-"):remove(d,"-")
      if instr(n,".")=0 then n+="."
      if instr(d,".")=0 then d+="."
      Static As String s(0 To 10)'to hold digit multiples
      Dim As String ans
      Dim As Long k,flag
      Var dm=decmove(d)
      adjustnumerator(n,dm)
      n+=String(places-decmove(n),"0")
      remove(d,".")
      d=ltrim(d,"0")
      n=ltrim(n,"0")
      Dim As String top=Mid(n,1,1)
      'if instr(top,".") then top=rtrim(top,"0"):remove(top,".")
      remove(top,".")
      For m As Long=1 To Len(n)-1

            If n[m]=Asc(".") Then flag=1:goto nxt
            
            For k =0 To 10 'range up the digits until the product > the numerator part
                  If k=10 Then
                        s(k)=d+"0"
                  Else
                        s(k)=smult(d,k)
                  End If
                  If isbigger(s(k),top) Then ans+=Str(k-1):Exit For
            Next k
            top=minus(top,s(k-1))
            
            top+=Chr(n[m])'bring down the next digit
            if flag=1 then ans+=".":flag=0
            nxt:
      Next m
      
      ans=iif(instr(ans,"."),Ltrim(ans,"0"),ans)
      if len(ans)=0 then ans="0"
      
      Return sign+ ans
End Function


Function GetFibonacci(n As Long) As String
    Static As Ubyte Qmod(0 To 19)={48,49,50,51,52,53,54,55,56,57,48,49,50,51,52,53,54,55,56,57}
    Static As Ubyte Qbool(0 To 19)={0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1}
    Dim As Ubyte addup,addcarry
    Dim As Long diff,n_,LL
    If n=0 Then Return "0"
    If n=1 Or n=2 Then Return "1"
    If n=3 Then
        Return "2"
    Else
        Dim As String sl="1",l="2",term
        For x As Long= 1 To n-3
            LL=Len(l)
            diff=-(LL<>Len(sl))
            term="0"+l
            addcarry=0
            For n_=LL-1 To diff Step -1
                addup=sl[n_-diff]+l[n_]-96
                term[n_+1]=Qmod(addup+addcarry)
                addcarry=Qbool(addup+addcarry)
            Next n_
            If addcarry=0 Then  term=Ltrim(term,"0"):Goto skip
            If n_=-1      Then term[0]=addcarry+48  :Goto skip
            For n_=n_ To 0 Step -1
                addup=l[n_]-48
                term[n_+1]=Qmod(addup+addcarry)
                addcarry=Qbool(addup+addcarry)
                If addcarry=0 Then Exit For
            Next n_
            term[0]=addcarry+48
            If addcarry=0 Then term=Ltrim(term,"0")
            skip:
            sl=l
            l=term
        Next x
        Function =term
    End If
End Function
'200 185000
'150 100000
'100=40000
'50=10000
'20=2000
'10=400
'5=80
'3=30
dim as double t=timer
dim as integer n=100+10

dim as string f=string(n-1,"9")+"8"+string(n,"9")
dim as integer places= _
 8.039594990255599 _
+ 1.338375048882362*n _
+0.8882107460008227*n^2 _
+0.4071177958035269*n^3 _
 -1.335359078473998e-002 *n^4 _
+1.708530082532032e-004 *n^5 _
 -9.296417579629874e-007 *n^6 _
+ 1.8032897354154e-009 *n^7

f= divide("1",f,places)
dim as integer counter
'get the counter
for i as integer=1 to len(f) step n'91
    counter+=1
next i

width 150,1000
print "START, the Fibonacci numbers are held in the division"
print f
print
print "The extracted Fibonnaci numbers"
counter=0
dim as string s,dd
for i as integer=1 to len(f) step n'91
      if ltrim(mid(f,i+1,n),"0")=getfibonacci(counter )then dd="  OK  " else dd="  ?  "
    print  "FIB ";str(counter);" = ";ltrim(mid(f,i+1,n),"0");"  ";dd',getfibonacci(counter )
    counter+=1
next i
print s
    print
    print "Check the last FIB number independently"
    print getfibonacci(487)
    print "time for everything "; timer-t
sleep

 
  
Last edited by dodicat on Aug 25, 2021 21:55, edited 2 times in total.
srvaldez
Posts: 3373
Joined: Sep 25, 2005 21:54

Re: string math

Post by srvaldez »

nice demo dodicat, here's the Fibonacci demo by division using the a fore mentioned algorithm

Code: Select all

' using the algorithm by Dr. David M. Smith
' "A Multiple-Precision Division Algorithm"
' http://dmsmith.lmu.build/MComp1996.pdf
'
' adapted to FreeBASIC by srvaldez
'
' Fibonacci demo

declare Sub divide (result() As Double, n() As Double, d() As Double)

Const dimension = 12096 ' 48380 digits
Dim As Double n(1 To dimension), d(1 To dimension)
Dim As Double result(1 To dimension), j, t
Dim As String s, digit
dim as integer i

' n(1) holds the exponent of n, n(2) holds the first digit of the numerator
' in sci notation it's 10^1 * .1 or .1e1
n(1) = 1: n(2) = 1
' d(1) holds the exponent of d, d(2) holds the first digit of the denominator
' d(3) an onward hold the rest of the denominator
d(1) = 202: d(2) = 9

For j = 3 To 50 + 3
    d(j) = 9999
Next
d(24 + 3) = 9998
d(50 + 3) = 9000

t = Timer
    divide(result(), n(), d())
t = Timer - t

s = ""
For j = 2 To UBound(n)
    digit = Trim(Str(result(j)))
    While Len(digit) < 4
        digit = "0" + digit
    Wend
    s = s + digit
Next
s = String(-result(1)-1, "0") + s
j = 1

For i = 1 To 481
    Print Mid(s, j, 100)
    j = j + 101
Next

Print "elapsed time for the division is "; t; " seconds"

Sleep

Function min (a As Long, b As Long) as long
    If a < b Then min = a Else min = b
End Function

Function RealW (w() As Double, j As Long) as double
    Dim wx As Double
    wx = ((w(j - 1) * 10000 + w(j)) * 10000 + w(j + 1)) * 10000
    If UBound(w) >= (j + 2) Then wx = wx + w(j + 2)
    RealW = wx
End Function

Sub subtract (w() As Double, q As Long, d() As Double, ka As Long, kb As Long)
    Dim As Long j
    For j = ka To kb
        w(j) = w(j) - q * d(j - ka + 2)
    Next
End Sub

Sub normalize (w() As Double, ka As Long, q As Long)
    w(ka) = w(ka) + w(ka - 1) * 10000
    w(ka - 1) = q
End Sub

Sub finalnorm (w() As Double, kb As Long)
    Dim As Long carry, j
    For j = kb To 3 Step -1
        If w(j) < 0 Then
            carry = ((-w(j) - 1) \ 10000) + 1
        Else
            If w(j) >= 10000 Then
                carry = -(w(j) \ 10000)
            Else
                carry = 0
            End If
        End If
        w(j) = w(j) + carry * 10000
        w(j - 1) = w(j - 1) - carry
    Next
End Sub

Sub divide (result() As Double, n() As Double, d() As Double)
    Dim As Long b, j, last, laststep, q, t
    Dim As Long stp
    Dim As Double xd, xn, round
    Dim As Double w(1 To UBound(n) + 4)
    b = 10000
    For j = UBound(n) To UBound(w)
        w(j) = 0
    Next
    t = UBound(n) - 1
    w(1) = n(1) - d(1) + 1
    w(2) = 0
    For j = 2 To UBound(n)
        w(j + 1) = n(j)
    Next
    xd = (d(2) * b + d(3)) * b + d(4) + d(5) / b
    laststep = t + 2
    For stp = 1 To laststep
        xn = RealW(w(), (stp + 2))
        q = Int(xn / xd)
        last = min(stp + t + 1, UBound(W))
        subtract(w(), q, d(), stp + 2, last)
        normalize(w(), stp + 2, q)
    Next
    finalnorm(w(), laststep + 1)
    If w(2) <> 0 Then laststep = laststep - 1
    round = w(laststep + 1) / b
    If round >= 0.5 Then w(laststep) = w(laststep) + 1
    If w(2) = 0 Then
        For j = 1 To t + 1
            result(j) = w(j + 1)
        Next
    Else
        For j = 1 To t + 1
            result(j) = w(j)
        Next
    End If
    If w(2) = 0 Then result(1) = w(1) - 1 Else result(1) = w(1)
End Sub
srvaldez
Posts: 3373
Joined: Sep 25, 2005 21:54

Re: string math

Post by srvaldez »

I mentioned that the multiplication was not fast, however I added an optional precision parameter so as to allow the Newton Raphson method to vary the precision, so for example the first estimate is done using double precision at approximately 16 digits, for the reciprocal that only uses multiplication and addition I use 40 digits for the first iteration, 72 for the second, 136 for the third and so on
for the Fibonacci example above the reciprocal is about 50% to 60% faster than the division, but the speed advantage is only there when you reach a certain degree of high precision, for lower precision, division is much faster than the reciprocal
Post Reply