need help translating from FreeBASIC to BCX

[jbk]

it's an interesting exercise to implement your own floating point arithmetic, but sadly to realize that it's too slow except for a calculator application
I wanted to test how difficult/easy it would be to port my decfloat arithmetic  to BCX and also since BCX is so close to C whether it would be faster than FB
that's why I translated the division routine to BCX, it performed about the same with -O2 but about twice as fast with -Ofast, although I prefer -O2 because in my opinion it's more safer/reliable

[MrBcx]

Thanks Robert ...

I sometimes remove chaff while trying to keep the forum coherent and useful.
Multitasking has a way of screwing up my intentions. 

My FP arithmetic code is back now, ready to be ignored for years to come. 

[Robert]

JB,

I thought I uploaded the following code some months ago but I can't find it now.

25 years ago when I was dabbling in MSDOS assembly language, I created
the following library that uses strings to perform long multiplication, addition, and
subtraction.  The code was written for my Basic to Asm Translator (BASM) but it also
works with BCX and any C\C++ compiler.

Hi MrBCX:

The topic seems to have been erased. I bookmarked it. It was at

https://bcxbasiccoders.com/smf/index.php?topic=524.msg2328#msg2328

posted on: June 1 2021 at 10:00:12 AM

[jbk]

thank you

[MrBcx]

JB,

I thought I uploaded the following code some months ago but I can't find it now.

25 years ago when I was dabbling in MSDOS assembly language, I created
the following library that uses strings to perform long multiplication, addition, and
subtraction.  The code was written for my Basic to Asm Translator (BASM) but it also
works with BCX and any C\C++ compiler.

Included in the code below is an example that multiplies

999999999999999999.123456789  *  999999999999999999.123456789

yielding a product of:

999999999999999998246913578000000000.768328000750190521

which has been independently verified using this online calculator:

https://keisan.casio.com/calculator

Maybe you will find it useful.

Code Select Expand


'****************************************************************************
'*  This library was designed specially for the BASM compiler (version 2).  *
'*  Copyrighted 1995 by Kevin T. Diggins.  This code may be distributed     *
'*  freely; if you make improvements to it, send me a copy.  Thanks...      *
'****************************************************************************
'* These FUNCTIONS take 2 strings of ascii digits as their parameters. You  *
'* should add code which limits the lengths of the input.  These FUNCTIONS  *
'* maintain excellent precision for most day to day floating point calc's.  *
'* However, these are not trapped for length overflow.  They work within    *
'* the restrictions of the max string length in your program. What this     *
'* means is if the max string length in your program is 80 chars, you must  *
'* insure that the result from these routines does not exceed that length.  *
'* In fact,  you should probably limit the result to something like:        *
'*                                                                          *
'*           max string length - 10, just to be on the safe side.           *
'*                                                                          *
'* * Remember *  BASM gives you control over the max string length with     *
'*               the $STRING meta-command.  See the help file for more...   *
'****************************************************************************
'* Function fpAdd$ (a$,b$)    addition                                      *
'* Function fpSub$ (a$,b$)    subtraction                                   *
'* Function fpMul$ (a$,b$)    multiplication                                *
'****************************************************************************
'* ** NOTES ** fpadd$ & fpSub$ use each other.  fpMul$ uses fpAdd$ and the  *
'*             functions Lstrip$ and Rstrip$, also in this this file.       *
'****************************************************************************

? fpAdd$ ("999.123456789","100")

? fpSub$ ("500.987654321","100")

? fpMul$ ("999999999999999999.123456789","999999999999999999.123456789")


FUNCTION fpAdd$ (X$,Y$)
  '--------------------------------
  LOCAL Xpoint,Ypoint
  LOCAL SizeA,Carry
  LOCAL Digit1,Digit2,Digit3
  LOCAL A,B,C
  LOCAL Integer1$,Integer2$
  LOCAL Dec1$,Dec2$
  LOCAL Num1$,Num2$
  LOCAL Tmp1$,Tmp2$
  LOCAL Accum$,Result$,Tmp$
  LOCAL Special
  '--------------------------------
  A=INSTR(X$,"-")
  B=INSTR(Y$,"-")

  IF A>0 THEN
    IF B=0 THEN
      Special=1
    END IF
  END IF

  IF A=0 THEN
    IF B>0 THEN
      Special=2
    END IF
  END IF

  IF A>0 THEN
    IF B>0 THEN
      Special = 3
    END IF
  END IF


  IF Special = 1 THEN
    Tmp1$=REMOVE$(X$,"-")
    Tmp2$=REMOVE$(Y$,"-")
    FUNCTION=fpSub$(Tmp2$,Tmp1$)
  END IF


  IF Special = 2 THEN
    Tmp1$=REMOVE$(X$,"-")
    Tmp2$=REMOVE$(Y$,"-")
    FUNCTION=fpSub$(Tmp1$,Tmp2$)
  END IF

  '--------------------------------
  Xpoint = INSTR(X$,".")
  Ypoint = INSTR(Y$,".")

  Integer1$ = EXTRACT$(X$,".")
  Integer2$ = EXTRACT$(Y$,".")

  Integer1$ = REMOVE$(Integer1$,"-")
  Integer2$ = REMOVE$(Integer2$,"-")

  A = LEN(Integer1$)
  B = LEN(Integer2$)

  DO
    C = LEN(Integer1$)
    IF C < B THEN
      Integer1$ = "0"  +  Integer1$
    ELSE
      EXIT LOOP
    END IF
  LOOP

  DO
    C = LEN(Integer2$)
    IF C < A THEN
      Integer2$ = "0"  +  Integer2$
    ELSE
      EXIT LOOP
    END IF
  LOOP

  'the length of both integers are equal at this point

  IF Xpoint = 0 THEN Dec1$ = "."
  IF Ypoint = 0 THEN Dec2$ = "."

  IF Xpoint THEN Dec1$ = Dec1$ + MID$(X$,Xpoint)
  IF Ypoint THEN Dec2$ = Dec2$ + MID$(Y$,Ypoint)

  A = LEN(Dec1$)
  B = LEN(Dec2$)

  DO
    C = LEN(Dec1$)
    IF C < B THEN
      Dec1$ = Dec1$  +  "0"
    ELSE
      EXIT LOOP
    END IF
  LOOP

  DO
    C = LEN(Dec2$)
    IF C < A THEN
      Dec2$ = Dec2$  +  "0"
    ELSE
      EXIT LOOP
    END IF
  LOOP

  'the length of both decimals are equal at this point

  Num1$ = Integer1$ + Dec1$
  Num2$ = Integer2$ + Dec2$

  'the length of both Numbers are equal at this point

  SizeA = LEN(Num1$)

  FOR A = SizeA TO 1 STEP -1
    Tmp1$ = MID$(Num1$,A,1)
    Tmp2$ = MID$(Num2$,A,1)

    IF Tmp1$  < > "." THEN
      Digit1 = VAL(Tmp1$)
      Digit2 = VAL(Tmp2$)
      Digit3 = Digit1  +  Digit2  +  Carry

      IF Digit3 > =  10 THEN
        Digit3  =  Digit3 - 10
        Carry = 1
      ELSE
        Carry  =  0
      END IF
      Accum$  =  Accum$  +  STR$(Digit3)
    ELSE
      Accum$ = Accum$  +  "."
    END IF
  NEXT

  Accum$ = REMOVE$(Accum$," ")

  A = LEN(Accum$)

  FOR B = A TO 1 STEP -1
    Tmp$ = MID$(Accum$,B,1)
    Result$ = Result$ + Tmp$
  NEXT

  IF Carry THEN Result$ ="1" + Result$

  IF Special = 3 THEN Result$ ="-" + Result$

  FUNCTION = Result$
END FUNCTION




FUNCTION fpSub$ (X$,Y$)
  '--------------------------------
  LOCAL Xpoint,Ypoint
  LOCAL SizeA,Borrow
  LOCAL Digit1,Digit2,Digit3
  LOCAL A,B,C
  LOCAL Integer1$,Integer2$
  LOCAL Dec1$,Dec2$
  LOCAL Num1$,Num2$
  LOCAL Tmp1$,Tmp2$
  LOCAL Accum$,Result$,Tmp$
  LOCAL Special,Negated
  '--------------------------------

  A=INSTR(X$,"-")
  B=INSTR(Y$,"-")

  IF A>0 THEN
    IF B=0 THEN
      Special=1
    END IF
  END IF

  IF A=0 THEN
    IF B>0 THEN
      Special=2
    END IF
  END IF

  IF Special = 1 THEN
    Tmp1$=REMOVE$(X$,"-")
    Tmp2$=REMOVE$(Y$,"-")
    Tmp1$=fpAdd$(Tmp1$,Tmp2$)
    Tmp1$="-"+Tmp1$
    REPLACE "--" WITH "-" IN Tmp1$
    FUNCTION=Tmp1$
  END IF

  IF Special = 2 THEN
    Tmp1$=REMOVE$(X$,"-")
    Tmp2$=REMOVE$(Y$,"-")
    FUNCTION=fpAdd$(Tmp1$,Tmp2$)
  END IF

  '--------------------------------

  Xpoint = INSTR(X$,".")
  Ypoint = INSTR(Y$,".")

  Integer1$ = EXTRACT$(X$,".")
  Integer2$ = EXTRACT$(Y$,".")

  A = LEN(Integer1$)
  B = LEN(Integer2$)

  DO
    C = LEN(Integer1$)
    IF C < B THEN
      Integer1$ = "0"  +  Integer1$
    ELSE
      EXIT LOOP
    END IF
  LOOP

  DO
    C = LEN(Integer2$)
    IF C < A THEN
      Integer2$ = "0"  +  Integer2$
    ELSE
      EXIT LOOP
    END IF
  LOOP

  'the length of both integers are equal at this point

  IF Xpoint = 0 THEN
    Dec1$ = "."
  END IF

  IF Ypoint = 0 THEN
    Dec2$ = "."
  END IF

  IF Xpoint THEN
    Dec1$ = Dec1$ + MID$(X$,Xpoint)
  END IF

  IF Ypoint THEN
    Dec2$ = Dec2$ + MID$(Y$,Ypoint)
  END IF

  A = LEN(Dec1$)
  B = LEN(Dec2$)

  DO
    C = LEN(Dec1$)
    IF C < B THEN
      Dec1$ = Dec1$  +  "0"
    ELSE
      EXIT LOOP
    END IF
  LOOP

  DO
    C = LEN(Dec2$)
    IF C < A THEN
      Dec2$ = Dec2$  +  "0"
    ELSE
      EXIT LOOP
    END IF
  LOOP

  'the length of both decimals are equal at this point

  Num1$ = Integer1$ + Dec1$
  Num2$ = Integer2$ + Dec2$

  'the length of both Numbers are equal at this point

  SizeA = LEN(Num1$)

  IF Num1$<Num2$ THEN
    SWAP Num1$,Num2$
    Negated=1
  END IF


  FOR A = SizeA TO 1 STEP -1

    Tmp1$ = MID$(Num1$,A,1)
    Tmp2$ = MID$(Num2$,A,1)

    IF Tmp1$  < > "." THEN
      Digit1 = VAL(Tmp1$)
      Digit2 = VAL(Tmp2$)

      Digit3 = Digit1 -  Digit2 - Borrow

      IF Digit3 < 0 THEN
        Digit3  =  Digit3 + 10
        Borrow = 1
      ELSE
        Borrow  =  0
      END IF
      Accum$  =  Accum$  +  STR$(Digit3)
    ELSE
      Accum$ = Accum$  +  "."
    END IF
  NEXT

  Accum$ = REMOVE$(Accum$," ")

  A = LEN(Accum$)

  FOR B = A TO 1 STEP -1
    Tmp$ = MID$(Accum$,B,1)
    Result$ = Result$ + Tmp$
  NEXT

  A=1

  DO
    Tmp$ = MID$(Result$,A,1)
    IF Tmp$ <> "0" THEN
      Result$ = MID$(Result$,A)
      EXIT LOOP
    END IF
    INCR A
  LOOP
  IF Negated THEN
    Result$ ="-" + Result$
  END IF

  FUNCTION=Result$
END FUNCTION




FUNCTION fpMul$(X$,Y$)
  '---------------------------------------------------------------------------
  ' This function requires the use of Lstrip$, Rstrip$ and FPadd$.  This
  ' function performs minimal checking.  It seems to treat negative and
  ' positive numbers (both input and output) correctly.  For routine calc's,
  ' this should do a good job of maintaining precision.  It can handle much
  ' larger numbers but you will need to use the $STRING statement to increase
  ' the memory needed by the routines to handle it.
  '---------------------------------------------------------------------------
  LOCAL AAA$,BBB$,CCC$,DDD$,Lt$,Rt$,Row$,Temp$
  LOCAL A,B,C
  LOCAL Tmp,Size1,Size2
  LOCAL Dec1,Dec2,Dec3
  LOCAL Top,Bottom,Zeros
  LOCAL Sign

  A=INSTR(X$,"-")
  B=INSTR(Y$,"-")

  IF A THEN Sign = 1
  IF B THEN Sign = Sign +1      '   Sign=1=neg    Sign=2=pos    Sign=0=pos

  AAA$ = REMOVE$(X$,"-")
  BBB$ = REMOVE$(Y$,"-")

  Tmp = INSTR(AAA$,".")
  IF NOT Tmp THEN AAA$ = AAA$ + "."
  Tmp = INSTR(BBB$,".")
  IF NOT Tmp THEN BBB$ = BBB$ + "."
  '-----------------------
  AAA$ = Rstrip$(AAA$,"0")
  BBB$ = Rstrip$(BBB$,"0")
  AAA$ = Lstrip$(AAA$,"0")
  BBB$ = Lstrip$(BBB$,"0")
  '-----------------------
  Size1 = LEN(AAA$)
  Size2 = LEN(BBB$)
  Dec1 = INSTR(AAA$,".")
  Dec1 = Size1 - Dec1
  Dec2 = INSTR(BBB$,".")
  Dec2 = Size2 - Dec2
  Dec3 = Dec1 + Dec2           'final number of Decimal places in the result
  '---------------------------------------------------------------------------
  AAA$ = REMOVE$(AAA$,".")
  BBB$ = REMOVE$(BBB$,".")
  '--------------------------
  IF BBB$ > AAA$ THEN SWAP AAA$,BBB$
  '--------------------------
  Top    = LEN(BBB$)
  Bottom = LEN(AAA$)
  '--------------------------
  FOR A = Bottom TO 1 STEP -1
    Row$ = "0"
    DDD$ = MID$(AAA$,A,1)
    B = VAL(DDD$)
    FOR C = 1 TO B
      Row$ = fpAdd$(Row$,BBB$)
    NEXT
    Row$ = Rstrip$(Row$,".")
    '-----------------------
    Temp$ = REPEAT$(Zeros,"0")
    Row$ = Row$ + Temp$            'this adds the right 0's to each level
    INCR Zeros                    'effectively multiplying by a power of ten
    '-----------------------
    CCC$ = fpAdd$(CCC$,Row$)
  NEXT
  CCC$ = Rstrip$(CCC$,".")
  '------------------------------------
  'now work in the Decimal point -- we subtract Dec3 from the LEN(CCC$)
  '------------------------------------
  Tmp  =  LEN(CCC$)

  WHILE Tmp < Dec3
    CCC$ = "0" + CCC$
    Tmp  =  LEN(CCC$)
  WEND

  Tmp  =  Tmp-Dec3
  Lt$  =  MID$(CCC$,1,Tmp)
  INCR    Tmp
  Rt$  =  MID$(CCC$,Tmp)
  CCC$   =  Lt$ + "." + Rt$
  CCC$   =  Lstrip$(CCC$,"0")
  CCC$   =  Rstrip$(CCC$,".")

  IF Sign = 1 THEN CCC$ = "-" + CCC$     'Sign=1=neg    Sign=2=pos    Sign=0=pos

  FUNCTION  =  CCC$        'and we return it as a string function result
END FUNCTION






FUNCTION Lstrip$(AAA$,B$)
  '----------------------------------------------------------------
  'Returns a copy of AAA$ with leading B$ Chars trimmed from the left
  '----------------------------------------------------------------
  LOCAL A,Tmp1$
  Tmp1$ = B$
  WHILE Tmp1$ = B$
    INCR A
    Tmp1$ = MID$(AAA$,A,1)
  WEND
  FUNCTION=MID$(AAA$,A)
END FUNCTION




FUNCTION Rstrip$(AAA$,B$)
  '------------------------------------------------------------------
  'Returns a copy of AAA$ with trailing B$ Chars trimmed from the right
  '------------------------------------------------------------------
  LOCAL A,Tmp1$
  A = LEN(AAA$)
  INCR A
  Tmp1$ = B$
  WHILE Tmp1$ = B$
    DECR A
    Tmp1$ = MID$(AAA$,A,1)
  WEND
  FUNCTION = MID$(AAA$,1,A)
END FUNCTION



[jbk]

I was not aware of the function round, but I plan on becoming more familiar with BCX
btw the time taken to divide 1 by a 200 digit number to 48380 digits took less than .05 seconds

[MrBcx]

MrBcx,

I managed to make it work, there a number of differences between FB and BCX

Good solution, fewer changes needed.

I still advise against using reserved keywords like "round" as variables.

https://bcxbasiccoders.com/webhelp/html/roundfunction.htm

[jbk]

MrBcx
I managed to make it work, there a number of differences between FB and BCX

Code Select Expand



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

'declare Sub divide (result[] As Double, n[] As Double, d[] As Double)

Const dimension = 12096 ' 48380 digits
Dim n[dimension+1] As Double, d[dimension+1] As Double
Dim result[dimension+1] As Double
dim as single t
Dim s As String, digit As String
dim i as integer, j as integer

' 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[0], &n[0], &d[0])
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 = "0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000" + s

j = 103
Print "0"
For i = 2 To 481
    digit=Mid$(s, j, 100)
    digit=ltrim$(digit, 48)
    Print digit
    j = j + 101
Next

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

End


'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, byval j As Long) as double
    Dim wx As Double
    wx = ((w[j - 1] * 10000 + w[j]) * 10000 + w[j + 1]) * 10000
    If (dimension+4) >= (j + 2) Then wx = wx + w[j + 2]
    function = wx
End Function

Sub subtract (w[] As Double, byval q As Long, d[] As Double, byval ka As Long, byval 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, byval ka As Long, byval q As Long)
    w[ka] = w[ka] + w[ka - 1] * 10000
    w[ka - 1] = q
End Sub

Sub finalnorm (w[] As Double, byval 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 w[dimension + 4+1] As Double

    b = 10000
    For j = 0 To dimension + 4
        w[j] = 0
    Next
    t = dimension - 1
    w[1] = n[1] - d[1] + 1
    w[2] = 0
    For j = 2 To dimension
        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[0], (stp + 2))
        q = Int(xn / xd)
        last = min(stp + t + 1, dimension+4)
        subtract(&w[0], q, &d[0], stp + 2, last)
        normalize(&w[0], stp + 2, q)
    Next
    finalnorm(&w[0], 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


[MrBcx]

MrBcx
unfortunately my code uses Ubound in a number of places which doesn't work on arrays passed to a sub, you simply get 0
I tried passing the array like this n[] but that does not work

You have more work to do ...

For example, this

SUB divide (result[] AS DOUBLE, n[] AS DOUBLE, d[] AS DOUBLE)


will need to be changed to something like this:


SUB divide (result[] AS DOUBLE, n[] AS DOUBLE, d[] AS DOUBLE, UBR%, UBN%, UBD%)



and change stuff like this:

DIM AS DOUBLE w[UBOUND(n) + 4]

to

DIM AS DOUBLE w[UBN% + 4]


and change stuff like this:

  FOR j = UBOUND(n) TO UBOUND(w)
    w[j] = 0
  NEXT

to

FOR j = UBN TO UBOUND(w)
    w[j] = 0
NEXT

[MrBcx]

JB - I'm curious - Do you know if FreeBasic version -worked- correctly?

From the PDF that you linked, the original algorithm was written in FORTRAN,  then
translated to MATHMATICA, then to pascal, then to FreeBasic, then to BCX then to C\C++.

What could possibly go wrong?   

Seriously,  I once wrote a special tool that relied on FORTRAN code that I converted
to C using the AT&T Fortran to C converter.  BCX user Don Baldwin helped me iron out
the final wrinkle to get the entire translation working correctly.  Don programmed in
FORTRAN many years ago, so I was very lucky to have his assistance.

[jbk]

MrBcx
unfortunately my code uses Ubound in a number of places which doesn't work on arrays passed to a sub, you simply get 0
I tried passing the array like this n[] but that does not work

[MrBcx]

I got a successful compile but the result is 0

Some warnings by various compilers

You have global and local variables using the same name.
That is perfectly legal but could cause problems depending
on the circumstance of their use.

declaration of 'i' hides global declaration   
declaration of 'n' hides global declaration
declaration of 'j' hides global declaration
declaration of 'd' hides global declaration
declaration of 't' hides global declaration
declaration of 'result' hides global declaration

CLANG reports 5 instances where you are attempting to take the UBOUND value
of an array passed by reference.  That will return the size of the POINTER, not
the size of the array. 

The solution is to pass the size of the array as another function/sub argument.

[MrBcx]

My changes:

The variable "ROUND" was changed to ROUNDD - ROUND() is a built-in BCX BASIC function
Arrays in FUNCTION Args must always be passed by reference 
Scalars in FUNCTION Args are ALWAYS by value by default but can also be passed by reference
Pay attention to variable names - they are CASE SENSITIVE : w[] and W[] are not the same thing
Pass the address of variables and arrays to FUNCTIONS/SUBS using the & operator and  for  arrays, you must pass

Code Select Expand


[0], not  []


EXAMPLE:  divide(&result[0], &n[0], &d[0])

The following code translates with BCX 777 and compiles with
Lcc-Win32, Pelles C, MSVC, Embarcadero, Mingw, and Clang

Code Select Expand


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

'declare Sub divide (result[] As Double, n[] As Double, d[] As Double)

CONST dimension = 12096 ' 48380 digits
DIM n[dimension] AS DOUBLE, d[dimension] AS DOUBLE
DIM result[dimension] AS DOUBLE, t AS DOUBLE
DIM s AS STRING, digit AS STRING
DIM i AS INTEGER, j AS INTEGER

' 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[0], &n[0], &d[0])
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 = 102
PRINT "0"
FOR i = 2 TO 481
  digit=MID$(s, j, 100)
  digit=LTRIM$(digit, 0)
  PRINT digit
  j = j + 101
NEXT

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

END


'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]
FUNCTION = 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, ROUNDD
  DIM AS DOUBLE w[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[0], (stp + 2))
    q = INT(xn / xd)
    last = MIN(stp + t + 1, UBOUND(w))
    subtract(&w[0], q, &d[0], stp + 2, last)
    normalize(&w[0], stp + 2, q)
  NEXT
  finalnorm(&w[0], laststep + 1)
  IF w[2] <> 0 THEN laststep = laststep - 1
  ROUNDD = w[laststep + 1] / b
  IF ROUNDD >= 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



[jbk]

I got a successful compile but the result is 0

Code Select Expand



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

'declare Sub divide (result[] As Double, n[] As Double, d[] As Double)

Const dimension = 12096 ' 48380 digits
Dim n[dimension] As Double, d[dimension] As Double
Dim result[dimension] As Double
dim as single t
Dim s As String, digit As String
dim i as integer, j as integer

' 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[1], &n[1], &d[1])
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 = 102
Print "0"
For i = 2 To 481
    digit=Mid$(s, j, 100)
    digit=ltrim$(digit, 0)
    Print digit
    j = j + 101
Next

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

End


'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, byval 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]
    function = wx
End Function

Sub subtract (w[] As Double, byval q As Long, d[] As Double, byval ka As Long, byval 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, byval ka As Long, byval q As Long)
    w[ka] = w[ka] + w[ka - 1] * 10000
    w[ka - 1] = q
End Sub

Sub finalnorm (w[] As Double, byval 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 w[UBound(n) + 4] As Double
    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[1], (stp + 2))
        q = Int(xn / xd)
        last = min(stp + t + 1, UBound(w))
        subtract(&w[1], q, &d[1], stp + 2, last)
        normalize(&w[1], stp + 2, q)
    Next
    finalnorm(&w[1], 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


[jbk]

hello
trying to familiarize myself with BCX by translating some FB code to BCX, first the FB code then my attempt at BCX translation

Code Select Expand


' using the algorithm by Dr. David M. Smith
' "A Multiple-Precision Division Algorithm"
' http://dmsmith.lmu.build/MComp1996.pdf
'
' 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 = 102
Print "0"
For i = 2 To 481
    digit=Mid(s, j, 100)
    digit=ltrim(digit, "0")
    Print digit
    j = j + 101
Next

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

End

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

BCX translation

Code Select Expand


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

'declare Sub divide (result[] As Double, n[] As Double, d[] As Double)

Const dimension = 12096 ' 48380 digits
Dim n[dimension] As Double, d[dimension] As Double
Dim result[dimension] As Double, t As Double
Dim s As String, digit As String
dim i as integer, j as integer

' 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 = 102
Print "0"
For i = 2 To 481
    digit=Mid$(s, j, 100)
    digit=ltrim$(digit, 0)
    Print digit
    j = j + 101
Next

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

End


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

Function RealW (byval w[] As Double, byval 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]
    function = wx
End Function

Sub subtract (byval w[] As Double, byval q As Long, byval d[] As Double, byval ka As Long, byval 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 (byval w[] As Double, byval ka As Long, byval q As Long)
    w[ka] = w[ka] + w[ka - 1] * 10000
    w[ka - 1] = q
End Sub

Sub finalnorm (byval w[] As Double, byval 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 (byval result[] As Double, byval n[] As Double, byval 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[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

the errors

Code Select Expand


BCX translated [Divide0.Bas] to [Divide0.Cpp] for a C++ Compiler
Compiling divide0 with MingW64
divide0.cpp: In function 'void divide(double*, double*, double*)':
divide0.cpp:397:18: error: expected primary-expression before ']' token
  397 |       xn=RealW(w[],(stp+2));
      |                  ^
divide0.cpp:400:18: error: expected primary-expression before ']' token
  400 |       subtract(w[],q,d[],stp+2,last);
      |                  ^
divide0.cpp:400:24: error: expected primary-expression before ']' token
  400 |       subtract(w[],q,d[],stp+2,last);
      |                        ^
divide0.cpp:401:19: error: expected primary-expression before ']' token
  401 |       normalize(w[],stp+2,q);
      |                   ^
divide0.cpp:403:15: error: expected primary-expression before ']' token
  403 |   finalnorm(w[],laststep+1);
      |               ^
divide0.cpp: In function 'int main(int, char**)':
divide0.cpp:452:17: error: expected primary-expression before ']' token
  452 |   divide(result[],n[],d[]);
      |                 ^
divide0.cpp:452:21: error: expected primary-expression before ']' token
  452 |   divide(result[],n[],d[]);
      |                     ^
divide0.cpp:452:25: error: expected primary-expression before ']' token
  452 |   divide(result[],n[],d[]);
      |                         ^