decimal floating point in PB

[jack]

here's my implementation of a decimal floating point lib
decfloat.pbi

Code: Select all

; functions that take number of digits instead of number of dwords as optional parameter
; fp2str
; fp2str_exp
; fp2str_fix
; fpeps
; pi_brent_salamin
; the rest take number of dwords as optional parameter

EnableExplicit

CompilerIf Not Defined(number_of_digits,#PB_Constant)
   #number_of_digits = 80
CompilerEndIf
#num_digits = #number_of_digits
CompilerIf ((#num_digits / 8)*8)<#num_digits
   #num_dwords = #num_digits / 8 +1
CompilerElse
   #num_dwords = #num_digits / 8
CompilerEndIf
#bias = 1073741824 ;2 ^ 30

Structure decfloat
   sign.l
   exponent.l
   mantissa.l[#num_dwords+1]
EndStructure

; error definitions

#divz_err = 1 ;divide by zero
#expo_err = 2 ;exponent overflow error

Procedure.s String(n.l, s.s)
   Define.s sn
   sn=Space(n)
   ReplaceString(sn," ",s,#PB_String_InPlace)
   ProcedureReturn sn
EndProcedure

CompilerIf Not Defined(maximum_Pi_Dec_precision,#PB_Constant)
   #maximum_Pi_Dec_precision = 1024
CompilerEndIf
Declare pi_brent_salamin (*pi_bs.decfloat, digits_in.l= #number_of_digits)
Global.decfloat Pi_dec
pi_brent_salamin (Pi_dec, #maximum_Pi_Dec_precision)

Procedure copydec(*n.decfloat, *result.decfloat)
   Define.l i
   *result\sign=*n\sign
   *result\exponent=*n\exponent
   For i=0 To #num_dwords
      *result\mantissa[i]=*n\mantissa[i]
   Next
EndProcedure

Procedure.s fp2str_exp (*n.decfloat, places.l=#number_of_digits-2)
   Define.l i, ex
   Define.s v, f, ts
   If *n\exponent > 0
      ex = (*n\exponent & $7fffffff) - #bias - 1
   Else
      ex = 0
   EndIf
   If *n\sign
      v = "-"
   Else
      v = " "
   EndIf
   ts = Str(*n\mantissa[0])
   ts = Trim(ts)
   If Len(ts) < 8
      ts = ts + String(8 - Len(ts), "0")
   EndIf
   v = v + Left(ts, 1) + "." + Mid(ts, 2)
   For i = 1 To #num_dwords - 1
      ts = Str(*n\mantissa[i])
      ts = Trim(ts)
      If Len(ts) < 8
         ts = String(8 - Len(ts), "0") + ts
      EndIf
      v = v + ts
   Next
   v = Left(v, places + 3)
   f = Trim(Str(Abs(ex)))
   
   ;f = String(9 - Len(f), "0") + f
   If ex < 0
      v = v + "e-"
   Else
      v = v + "e+"
   EndIf
   v = v + f
   ProcedureReturn v
EndProcedure

;long part of num
Procedure fptrunc (*result.decfloat, *num.decfloat)
   Define.decfloat ip
   Define.l ex, ex2, j, k
   
   ex = (*num\exponent & $7fffffff) - #bias
   If ex < 1
      copydec(ip, *result) ;result = ip
      ProcedureReturn
   EndIf
   If ex >= (#num_digits)
      copydec(*num, *result) ;result = num
      ProcedureReturn
   EndIf
   ex2 = ex / 8
   k = ex2
   j = ex % 8
   While ex2 > 0
      ex2 = ex2 - 1
      ip\mantissa[ex2]=*num\mantissa[ex2]
   Wend
   If j = 1
      ip\mantissa[k]= 10000000 * (*num\mantissa[k] / 10000000)
   ElseIf j = 2
      ip\mantissa[k]= 1000000 * (*num\mantissa[k] / 1000000)
   ElseIf j = 3
      ip\mantissa[k]= 100000 * (*num\mantissa[k] / 100000)
   ElseIf j = 4
      ip\mantissa[k]= 10000 * (*num\mantissa[k] / 10000)
   ElseIf j = 5
      ip\mantissa[k]= 1000 * (*num\mantissa[k] / 1000)
   ElseIf j = 6
      ip\mantissa[k]= 100 * (*num\mantissa[k] / 100)
   ElseIf j = 7
      ip\mantissa[k]= 10 * (*num\mantissa[k] / 10)
   ElseIf j = 8
      ip\mantissa[k]= *num\mantissa[k]
   EndIf
   ip\exponent = ex + #bias
   ip\sign = *num\sign
   copydec(ip, *result) ;result = ip
EndProcedure

Procedure.l fptrunc_is_odd (*num.decfloat)
   Define.l ex, j, k
   
   ex = (*num\exponent & $7fffffff) - #bias
   If ex < 1
      ProcedureReturn 0
   EndIf
   If ex >= (#num_digits)
      PrintN( "error in function fptrunc_is_odd")
      ProcedureReturn 99999999
   EndIf
   k = ex / 8
   j = ex % 8
   
   If j = 1
      ProcedureReturn (*num\mantissa[k] / 10000000) & 1
   ElseIf j = 2
      ProcedureReturn (*num\mantissa[k] / 1000000) & 1
   ElseIf j = 3
      ProcedureReturn (*num\mantissa[k] / 100000) & 1
   ElseIf j = 4
      ProcedureReturn (*num\mantissa[k] / 10000) & 1
   ElseIf j = 5
      ProcedureReturn (*num\mantissa[k] / 1000) & 1
   ElseIf j = 6
      ProcedureReturn (*num\mantissa[k] / 100) & 1
   ElseIf j = 7
      ProcedureReturn (*num\mantissa[k] / 10) & 1
   ElseIf j = 8
      ProcedureReturn *num\mantissa[k] & 1
   EndIf
   ProcedureReturn 0
EndProcedure

Procedure.d fp2dbl (*n.decfloat)
   Define.l ex
   Define.s v, f, ts
   If *n\exponent > 0
      ex = (*n\exponent & $7fffffff) - #bias - 1
   Else
      ex = 0
   EndIf
   If *n\sign
      v = "-"
   Else
      v = " "
   EndIf
   ts = Trim(Str(*n\mantissa[0]))
   If Len(ts) < 8
      ts = ts + String(8 - Len(ts), "0")
   EndIf
   v = v + Left(ts, 1) + "." + Mid(ts, 2)
   ts = Trim(Str(*n\mantissa[1]))
   If Len(ts) < 8
      ts = String(8 - Len(ts), "0") + ts
   EndIf
   v = v + ts
   
   f = Str(Abs(ex))
   f = String(5 - Len(f), "0") + f
   If ex < 0
      v = v + "e-"
   Else
      v = v + "e+"
   EndIf
   v = v + f
   ProcedureReturn ValD(v)
EndProcedure

Procedure.s fp2str_fix (*n.decfloat, places.l=#number_of_digits-2)
   Define.l i, ex
   Define.s v, ts, s
   If *n\exponent > 0
      ex = (*n\exponent & $7fffffff) - #bias - 1
   Else
      ex = 0
   EndIf
   If *n\sign
      s = "-"
   Else
      s = " "
   EndIf
   ts = Trim(Str(*n\mantissa[0]))
   If Len(ts) < 8
      ts = ts + String(8 - Len(ts), "0")
   EndIf
   v = ts
   For i = 1 To #num_dwords - 1
      ts = Trim(Str(*n\mantissa[i]))
      If Len(ts) < 8
         ts = String(8 - Len(ts), "0") + ts
      EndIf
      v = v + ts
   Next
   If places < #num_digits
      v = Left(v, places)
   EndIf
   If ex = 0
      v = Left(v, 1) + "." + Mid(v, 2)
   ElseIf ex < 0
      v = "0." + String(Abs(ex) - 1, "0") + v
   ElseIf ex > 0
      v = Left(v, ex + 1) + "." + Mid(v, ex + 2)
   EndIf
   ProcedureReturn s + v
EndProcedure

Declare fpadd (*result.decfloat, *x.decfloat, *y.decfloat, dwords_in.l=#num_dwords)
Declare fpsub (*result.decfloat, *x.decfloat, *y.decfloat, dwords_in.l=#num_dwords)
Declare si2fp (*result.decfloat, m.q, dwords_in.l=#num_dwords)

Procedure.s fp2str (*n.decfloat, places_in.l=#number_of_digits)
   Define.l ex, e, sgn, places=places_in
   Define.decfloat z
   If places>#number_of_digits-1
      places=#number_of_digits-2
   EndIf
   If *n\exponent <> 0
      ex = (*n\exponent & $7fffffff) - #bias - 1
   Else
      ex = 0
   EndIf
   If *n\exponent=0
      ProcedureReturn " 0"
   EndIf
   e=ex
   If e<0
      e=-e
   EndIf
   si2fp (z, 5)
   sgn=*n\sign
   If (ex > (-5)) And (e < places)
      z\exponent = (-(places+e) + #bias + 1)
      If sgn=0
         fpadd(z, z, *n)
      Else
         fpsub(z, *n, z)
      EndIf
      ProcedureReturn fp2str_fix(z, places)
   Else
      z\exponent = (-(places+e+1) + #bias + 1)
      If sgn=0
         fpadd(z, z, *n)
      Else
         fpsub(z, *n, z)
      EndIf
      ProcedureReturn fp2str_exp(z, places)
   EndIf
EndProcedure

; in case the above doesn't work right use this instead
; Procedure.s fp2str (*n.decfloat, places.l=#number_of_digits-2)
;    Define.l ex
;    If *n\exponent <> 0
;       ex = (*n\exponent & $7fffffff) - #bias - 1
;    Else
;       ex = 0
;    EndIf
;    If (Abs(ex) < places) And (places<(#num_digits-1))
;       ProcedureReturn fp2str_fix(*n, places)
;    Else
;       ProcedureReturn fp2str_exp(*n, places)
;    EndIf
; EndProcedure

Procedure.s fp2strT (*n.decfloat, places_in.l=#number_of_digits)
   Define.l ex, e, sgn, places=places_in, k
   Define.decfloat z
   Define.s s, s2
   If places>#number_of_digits-1
      places=#number_of_digits-2
   EndIf
   If *n\exponent <> 0
      ex = (*n\exponent & $7fffffff) - #bias - 1
   Else
      ex = 0
   EndIf
   If *n\exponent=0
      ProcedureReturn " 0"
   EndIf
   e=ex
   If e<0
      e=-e
   EndIf
   si2fp (z, 5)
   sgn=*n\sign
   If (ex > (-5)) And (e < places)
      z\exponent = (-(places+e) + #bias + 1)
      If sgn=0
         fpadd(z, z, *n)
      Else
         fpsub(z, *n, z)
      EndIf
      s=fp2str_fix(z, places)
      s=RTrim(s, "0")
      If Right(s,1)="."
         s=Left(s,Len(s)-1)
      EndIf
      ProcedureReturn s
   Else
      z\exponent = (-(places+e+1) + #bias + 1)
      If sgn=0
         fpadd(z, z, *n)
      Else
         fpsub(z, *n, z)
      EndIf
      s=fp2str_exp(z, places)
      k=FindString(s, "e")
      s2=Left(s, k-1)
      s2=RTrim(s2, "0")
      If Right(s2,1)="."
         s2=Left(s2,Len(s2)-1)
      EndIf
      s=s2+Right(s,Len(s)-k+1)
      ProcedureReturn s
   EndIf
EndProcedure

Procedure str2fp (*result.decfloat, value_in.s)
   Define.l j, s, d, e, ep, ulng
   Define.l ex, es, i, f, fp, fln
   Define.s c, f1, f2, f3, ts, value
   Define.decfloat n
   
   value=value_in
   j = 1
   s = 1
   d = 0
   e = 0
   ep = 0
   ex = 0
   es = 1
   i = 0
   f = 0
   fp = 0
   f1 = ""
   f2 = ""
   f3 = ""
   value = UCase(value)
   fln = Len(value)
   
   While j <= fln
      c = Mid(value, j, 1)
      If ep = 1
         If c = " "
            j = j + 1
            Goto skip_while
         EndIf
         If c = "-"
            es = -es
            c = ""
         EndIf
         If c = "+"
            j = j + 1
            Goto skip_while
         EndIf
         If (c = "0") And (f3 = "")
            j = j + 1
            Goto skip_while
         EndIf
         If (c > "/") And (c < ":") ;c is digit between 0 and 9
            f3 = f3 + c
            ex = 10 * ex + (Asc(c) - 48)
            j = j + 1
            Goto skip_while
         EndIf
      EndIf
      
      If c = " "
         j = j + 1
         Goto skip_while
      EndIf
      If c = "-"
         s = -s
         j = j + 1
         Goto skip_while
      EndIf
      If c = "+"
         j = j + 1
         Goto skip_while
      EndIf
      If c = "."
         If d = 1
            j = j + 1
            Goto skip_while
         EndIf
         d = 1
      EndIf
      If (c > "/") And (c < ":") ;c is digit between 0 and 9
         If ((c = "0") And (i = 0))
            If d = 0
               j = j + 1
               Goto skip_while
            EndIf
            If (d = 1) And (f = 0)
               e = e - 1
               j = j + 1
               Goto skip_while
            EndIf
         EndIf
         If d = 0
            f1 = f1 + c
            i = i + 1
         Else
            If (c > "0")
               fp = 1
            EndIf
            f2 = f2 + c
            f = f + 1
         EndIf
      EndIf
      If c = "E" Or c = "D"
         ep = 1
      EndIf
      j = j + 1
      skip_while:
   Wend
   If fp = 0
      f = 0
      f2 = ""
   EndIf
   
   If s = -1
      s = $8000
   Else
      s = 0
   EndIf
   n\sign = s
   ex = es * ex - 1 + i + e
   f1 = f1 + f2
   f1 = Mid(f1, 1, 1) + Right(f1, Len(f1) - 1)
   fln = Len(f1)
   If Len(f1) > (#num_digits + 1 + 8)
      f1 = Mid(f1, 1, (#num_digits + 1 + 8))
   EndIf
   While Len(f1) < (#num_digits + 1 + 8)
      f1 = f1 + "0"
   Wend
   j = 1
   For i = 0 To #num_dwords
      ts = Mid(f1, j, 8)
      ulng = Val(ts)
      n\mantissa[i]= ulng
      If ulng <> 0
         fp = 1
      EndIf
      j = j + 8
   Next
   If fp
      n\exponent = (ex + #bias + 1)
   Else
      n\exponent = 0
   EndIf
   copydec(n, *result) ;result = n
EndProcedure

Procedure si2fp (*result.decfloat, m.q, dwords_in.l=#num_dwords)
   Define.l dwords, i, ndw
   Define.q n
   Define.decfloat fac1
   dwords = dwords_in
   If dwords > #num_dwords
      dwords = #num_dwords
   EndIf
   
   n = Abs(m)
   ndw = n
   If n > 9999999999999999
      str2fp (fac1, Str(m))
      copydec(fac1, *result) ;result = fac1
      ProcedureReturn
   EndIf
   
   For i = 1 To dwords
      fac1\mantissa[i] = 0
   Next
   
   If m = 0
      fac1\exponent = 0
      fac1\sign = 0
      copydec(fac1, *result) ;result = fac1
      ProcedureReturn
   EndIf
   
   fac1\exponent = #bias
   If n < 100000000
      If n < 10
         fac1\mantissa[0] = n * 10000000
         fac1\exponent = fac1\exponent + 1
      ElseIf n < 100
         fac1\mantissa[0] = n * 1000000
         fac1\exponent = fac1\exponent + 2
      ElseIf n < 1000
         fac1\mantissa[0] = n * 100000
         fac1\exponent = fac1\exponent + 3
      ElseIf n < 10000
         fac1\mantissa[0] = n * 10000
         fac1\exponent = fac1\exponent + 4
      ElseIf n < 100000
         fac1\mantissa[0] = n * 1000
         fac1\exponent = fac1\exponent + 5
      ElseIf n < 1000000
         fac1\mantissa[0] = n * 100
         fac1\exponent = fac1\exponent + 6
      ElseIf n < 10000000
         fac1\mantissa[0] = n * 10
         fac1\exponent = fac1\exponent + 7
      ElseIf n < 100000000
         fac1\mantissa[0] = ndw
         fac1\exponent = fac1\exponent + 8
      EndIf
   EndIf
   If n > 99999999
      fac1\exponent = fac1\exponent + 8
      If n < 1000000000
         fac1\mantissa[0] = n / 10
         fac1\mantissa[1] = (n % 10) * 10000000
         fac1\exponent = fac1\exponent + 1
      ElseIf n < 100000000000
         fac1\mantissa[0] = n / 100
         fac1\mantissa[1] = (n % 100) * 1000000
         fac1\exponent = fac1\exponent + 2
      ElseIf n < 1000000000000
         fac1\mantissa[0] = n / 1000
         fac1\mantissa[1] = (n % 1000) * 100000
         fac1\exponent = fac1\exponent + 3
      ElseIf n < 10000000000000
         fac1\mantissa[0] = n / 10000
         fac1\mantissa[1] = (n % 10000) * 10000
         fac1\exponent = fac1\exponent + 4
      ElseIf n < 100000000000000
         fac1\mantissa[0] = n / 100000
         fac1\mantissa[1] = (n % 100000) * 1000
         fac1\exponent = fac1\exponent + 5
      ElseIf n < 1000000000000000
         fac1\mantissa[0] = n / 1000000
         fac1\mantissa[1] = (n % 1000000) * 100
         fac1\exponent = fac1\exponent + 6
      ElseIf n < 10000000000000000
         fac1\mantissa[0] = n / 10000000
         fac1\mantissa[1] = (n % 10000000) * 10
         fac1\exponent = fac1\exponent + 7
      ElseIf n < 100000000000000000
         fac1\mantissa[0] = n / 100000000
         fac1\mantissa[1] = n % 100000000
         fac1\exponent = fac1\exponent + 8
      EndIf
   EndIf
   If m < 0
      fac1\sign = $8000
   Else
      fac1\sign = 0
   EndIf
   copydec(fac1, *result) ;result = fac1
EndProcedure

Procedure rshift_1 (*mantissa.decfloat, dwords_in.l=#num_dwords)
   Define.l dwords, v1, v2, i
   dwords = dwords_in
   If dwords > #num_dwords
      dwords = #num_dwords
   EndIf
   For i = dwords To 1 Step -1
      v1 = *mantissa\mantissa[i] / 10
      v2 = *mantissa\mantissa[i - 1] % 10
      v2 = v2 * 10000000 + v1
      *mantissa\mantissa[i] = v2
   Next
   *mantissa\mantissa[0] = *mantissa\mantissa[0] / 10
EndProcedure

Procedure lshift_1 (*mantissa.decfloat, dwords_in.l=#num_dwords)
   Define.l dwords, v1, v2, i
   dwords = dwords_in
   If dwords > #num_dwords
      dwords = #num_dwords
   EndIf
   For i = 0 To dwords - 1
      v1 = *mantissa\mantissa[i] % 10000000
      v2 = *mantissa\mantissa[i + 1] / 10000000
      *mantissa\mantissa[i] = v1 * 10 + v2
      *mantissa\mantissa[i + 1] = *mantissa\mantissa[i + 1] % 10000000
   Next
   *mantissa\mantissa[dwords] = 10 * (*mantissa\mantissa[dwords] % 10000000)
EndProcedure

Procedure rshift_2 (*mantissa.decfloat, dwords_in.l=#num_dwords)
   Define.l dwords, v1, v2, i
   dwords = dwords_in
   If dwords > #num_dwords
      dwords = #num_dwords
   EndIf
   For i = dwords To 1 Step -1
      v1 = *mantissa\mantissa[i] / 100
      v2 = *mantissa\mantissa[i - 1] % 100
      v2 = v2 * 1000000 + v1
      *mantissa\mantissa[i] = v2
   Next
   *mantissa\mantissa[0] = *mantissa\mantissa[0] / 100
EndProcedure

Procedure lshift_2 (*mantissa.decfloat, dwords_in.l=#num_dwords)
   Define.l dwords, v1, v2, i
   dwords = dwords_in
   If dwords > #num_dwords
      dwords = #num_dwords
   EndIf
   For i = 0 To dwords - 1
      v1 = *mantissa\mantissa[i] % 1000000
      v2 = *mantissa\mantissa[i + 1] / 1000000
      *mantissa\mantissa[i] = v1 * 100 + v2
      *mantissa\mantissa[i + 1] = *mantissa\mantissa[i + 1] % 1000000
   Next
   *mantissa\mantissa[dwords] = 100 * (*mantissa\mantissa[dwords] % 1000000)
EndProcedure

Procedure rshift_3 (*mantissa.decfloat, dwords_in.l=#num_dwords)
   Define.l dwords, v1, v2, i
   dwords = dwords_in
   If dwords > #num_dwords
      dwords = #num_dwords
   EndIf
   For i = dwords To 1 Step -1
      v1 = *mantissa\mantissa[i] / 1000
      v2 = *mantissa\mantissa[i - 1] % 1000
      v2 = v2 * 100000 + v1
      *mantissa\mantissa[i] = v2
   Next
   *mantissa\mantissa[0] = *mantissa\mantissa[0] / 1000
EndProcedure

Procedure lshift_3 (*mantissa.decfloat, dwords_in.l=#num_dwords)
   Define.l dwords, v1, v2, i
   dwords = dwords_in
   If dwords > #num_dwords
      dwords = #num_dwords
   EndIf
   For i = 0 To dwords - 1
      v1 = *mantissa\mantissa[i] % 100000
      v2 = *mantissa\mantissa[i + 1] / 100000
      *mantissa\mantissa[i] = v1 * 1000 + v2
      *mantissa\mantissa[i + 1] = *mantissa\mantissa[i + 1] % 100000
   Next
   *mantissa\mantissa[dwords] = 1000 * (*mantissa\mantissa[dwords] % 100000)
EndProcedure

Procedure rshift_4 (*mantissa.decfloat, dwords_in.l=#num_dwords)
   Define.l dwords, v1, v2, i
   dwords = dwords_in
   If dwords > #num_dwords
      dwords = #num_dwords
   EndIf
   For i = dwords To 1 Step -1
      v1 = *mantissa\mantissa[i] / 10000
      v2 = *mantissa\mantissa[i - 1] % 10000
      v2 = v2 * 10000 + v1
      *mantissa\mantissa[i] = v2
   Next
   *mantissa\mantissa[0] = *mantissa\mantissa[0] / 10000
EndProcedure

Procedure lshift_4 (*mantissa.decfloat, dwords_in.l=#num_dwords)
   Define.l dwords, v1, v2, i
   dwords = dwords_in
   If dwords > #num_dwords
      dwords = #num_dwords
   EndIf
   For i = 0 To dwords - 1
      v1 = *mantissa\mantissa[i] % 10000
      v2 = *mantissa\mantissa[i + 1] / 10000
      *mantissa\mantissa[i] = v1 * 10000 + v2
      *mantissa\mantissa[i + 1] = *mantissa\mantissa[i + 1] % 10000
   Next
   *mantissa\mantissa[dwords] = 10000 * (*mantissa\mantissa[dwords] % 10000)
EndProcedure

Procedure rshift_5 (*mantissa.decfloat, dwords_in.l=#num_dwords)
   Define.l dwords, v1, v2, i
   dwords = dwords_in
   If dwords > #num_dwords
      dwords = #num_dwords
   EndIf
   For i = dwords To 1 Step -1
      v1 = *mantissa\mantissa[i] / 100000
      v2 = *mantissa\mantissa[i - 1] % 100000
      v2 = v2 * 1000 + v1
      *mantissa\mantissa[i] = v2
   Next
   *mantissa\mantissa[0] = *mantissa\mantissa[0] / 100000
EndProcedure

Procedure lshift_5 (*mantissa.decfloat, dwords_in.l=#num_dwords)
   Define.l dwords, v1, v2, i
   dwords = dwords_in
   If dwords > #num_dwords
      dwords = #num_dwords
   EndIf
   For i = 0 To dwords - 1
      v1 = *mantissa\mantissa[i] % 1000
      v2 = *mantissa\mantissa[i + 1] / 1000
      *mantissa\mantissa[i] = v1 * 100000 + v2
      *mantissa\mantissa[i + 1] = *mantissa\mantissa[i + 1] % 1000
   Next
   *mantissa\mantissa[dwords] = 100000 * (*mantissa\mantissa[dwords] % 1000)
EndProcedure

Procedure rshift_6 (*mantissa.decfloat, dwords_in.l=#num_dwords)
   Define.l dwords, v1, v2, i
   dwords = dwords_in
   If dwords > #num_dwords
      dwords = #num_dwords
   EndIf
   For i = dwords To 1 Step -1
      v1 = *mantissa\mantissa[i] / 1000000
      v2 = *mantissa\mantissa[i - 1] % 1000000
      v2 = v2 * 100 + v1
      *mantissa\mantissa[i] = v2
   Next
   *mantissa\mantissa[0] = *mantissa\mantissa[0] / 1000000
EndProcedure

Procedure lshift_6 (*mantissa.decfloat, dwords_in.l=#num_dwords)
   Define.l dwords, v1, v2, i
   dwords = dwords_in
   If dwords > #num_dwords
      dwords = #num_dwords
   EndIf
   For i = 0 To dwords - 1
      v1 = *mantissa\mantissa[i] % 100
      v2 = *mantissa\mantissa[i + 1] / 100
      *mantissa\mantissa[i] = v1 * 1000000 + v2
      *mantissa\mantissa[i + 1] = *mantissa\mantissa[i + 1] % 100
   Next
   *mantissa\mantissa[dwords] = 1000000 * (*mantissa\mantissa[dwords] % 100)
EndProcedure

Procedure rshift_7 (*mantissa.decfloat, dwords_in.l=#num_dwords)
   Define.l dwords, v1, v2, i
   dwords = dwords_in
   If dwords > #num_dwords
      dwords = #num_dwords
   EndIf
   For i = dwords To 1 Step -1
      v1 = *mantissa\mantissa[i] / 10000000
      v2 = *mantissa\mantissa[i - 1] % 10000000
      v2 = v2 * 10 + v1
      *mantissa\mantissa[i] = v2
   Next
   *mantissa\mantissa[0] = *mantissa\mantissa[0] / 10000000
EndProcedure

Procedure lshift_7 (*mantissa.decfloat, dwords_in.l=#num_dwords)
   Define.l dwords, v1, v2, i
   dwords = dwords_in
   If dwords > #num_dwords
      dwords = #num_dwords
   EndIf
   For i = 0 To dwords - 1
      v1 = *mantissa\mantissa[i] % 10
      v2 = *mantissa\mantissa[i + 1] / 10
      *mantissa\mantissa[i] = v1 * 10000000 + v2
      *mantissa\mantissa[i + 1] = *mantissa\mantissa[i + 1] % 10
   Next
   *mantissa\mantissa[dwords] = 10000000 * (*mantissa\mantissa[dwords] % 10)
EndProcedure

Procedure rshift_8 (*mantissa.decfloat, dwords_in.l=#num_dwords)
   Define.l dwords, i
   dwords = dwords_in
   If dwords > #num_dwords
      dwords = #num_dwords
   EndIf
   For i = dwords To 1 Step -1
      *mantissa\mantissa[i] = *mantissa\mantissa[i - 1]
   Next
   *mantissa\mantissa[0] = 0
EndProcedure

Procedure lshift_8 (*mantissa.decfloat, dwords_in.l=#num_dwords)
   Define.l dwords, i
   dwords.l = dwords_in
   If dwords > #num_dwords
      dwords = #num_dwords
   EndIf
   For i = 0 To dwords - 1
      *mantissa\mantissa[i] = *mantissa\mantissa[i + 1]
   Next
   *mantissa\mantissa[dwords] = 0
EndProcedure

Procedure.l fpcmp (*x.decfloat, *y.decfloat, dwords_in.l=#num_dwords)
   Define.l dwords
   dwords = dwords_in
   
   If dwords > #num_dwords
      dwords = #num_dwords
   EndIf
   Define.l c, i
   If *x\sign < *y\sign
      ProcedureReturn -1
   EndIf
   If *x\sign > *y\sign
      ProcedureReturn 1
   EndIf
   If *x\exponent<*y\exponent
      If *x\sign=0
         ProcedureReturn -1
      Else
         ProcedureReturn 1
      EndIf
   EndIf
   If *x\exponent>*y\exponent
      If *x\sign=0
         ProcedureReturn 1
      Else
         ProcedureReturn -1
      EndIf
   EndIf
   
   For i=0 To dwords
      c=*x\mantissa[i]-*y\mantissa[i] 
      If c<>0
         Break
      EndIf
   Next
   If c=0
      ProcedureReturn 0
   EndIf
   If c<0
      If *x\sign=0
         ProcedureReturn -1
      Else
         ProcedureReturn 1
      EndIf
   EndIf
   If c>0
      If *x\sign=0
         ProcedureReturn 1
      Else
         ProcedureReturn -1
      EndIf
   EndIf
EndProcedure

Procedure norm_fac1 (*fac1.decfloat, dwords_in.l=#num_dwords)
   Define.l dwords, i, er, f
   dwords = dwords_in
   If dwords > #num_dwords
      dwords = #num_dwords
   EndIf
   ; normalize the number in fac1
   ; all routines exit through this one.
   
   ;see if the mantissa is all zeros.
   ;if so, sit the exponent and sign equal to 0.
   er = 0: f = 0
   For i = 0 To dwords
      If *fac1\mantissa[i] > 0
         f = 1
      EndIf
   Next
   If f = 0
      *fac1\exponent = 0
      *fac1\sign = 0
      ProcedureReturn
      ;if the highmost digit in fac1_man is nonzero,
      ;shift the mantissa right 1 digit and
      ;increment the exponent
   ElseIf *fac1\mantissa[0] > 99999999
      rshift_1 (*fac1, dwords)
      *fac1\exponent = *fac1\exponent + 1
   Else
      ;now shift fac1_man 1 to the left until a
      ;nonzero digit appears in the next-to-highest
      ;digit of fac1_man.  decrement exponent for
      ;each shift.
      While *fac1\mantissa[0] = 0
         lshift_8 (*fac1, dwords)
         *fac1\exponent = *fac1\exponent - 8
         If *fac1\exponent = 0
            PrintN( "norm_fac1=expu_err")
            ProcedureReturn
         EndIf
      Wend
      If *fac1\mantissa[0] < 10
         lshift_7 (*fac1, dwords)
         *fac1\exponent = *fac1\exponent - 7
      ElseIf *fac1\mantissa[0] < 100
         lshift_6 (*fac1, dwords)
         *fac1\exponent = *fac1\exponent - 6
      ElseIf *fac1\mantissa[0] < 1000
         lshift_5 (*fac1, dwords)
         *fac1\exponent = *fac1\exponent - 5
      ElseIf *fac1\mantissa[0] < 10000
         lshift_4 (*fac1, dwords)
         *fac1\exponent = *fac1\exponent - 4
      ElseIf *fac1\mantissa[0] < 100000
         lshift_3 (*fac1, dwords)
         *fac1\exponent = *fac1\exponent - 3
      ElseIf *fac1\mantissa[0] < 1000000
         lshift_2 (*fac1, dwords)
         *fac1\exponent = *fac1\exponent - 2
      ElseIf *fac1\mantissa[0] < 10000000
         lshift_1 (*fac1, dwords)
         *fac1\exponent = *fac1\exponent - 1
      EndIf
   EndIf
   ;check for overflow/underflow
   If *fac1\exponent < 0
      PrintN( "norm_fac1=#expo_err")
   EndIf
EndProcedure

Procedure fpadd_aux (*fac1.decfloat, *fac2.decfloat, dwords_in.l=#num_dwords)
   Define.l dwords, v, c, i
   dwords = dwords_in
   If dwords > #num_dwords
      dwords = #num_dwords
   EndIf
   c = 0
   For i = dwords To 1 Step -1
      v = *fac2\mantissa[i] + *fac1\mantissa[i] + c
      If v > 99999999
         v = v - 100000000
         c = 1
      Else
         c = 0
      EndIf
      *fac1\mantissa[i] = v
   Next
   v = *fac1\mantissa[0] + *fac2\mantissa[0] + c
   *fac1\mantissa[0] = v
   
   norm_fac1 (*fac1, dwords)
EndProcedure

Procedure fpsub_aux (*fac1.decfloat, *fac2.decfloat, dwords_in.l=#num_dwords)
   Define.l dwords, v, c, i
   dwords = dwords_in
   If dwords > #num_dwords
      dwords = #num_dwords
   EndIf
   c = 0
   For i = dwords To 1 Step -1
      v = *fac1\mantissa[i] - *fac2\mantissa[i] - c
      If v < 0
         v = v + 100000000
         c = 1
      Else
         c = 0
      EndIf
      *fac1\mantissa[i] = v
   Next
   v = *fac1\mantissa[0] - *fac2\mantissa[0] - c
   *fac1\mantissa[0] = v
   
   norm_fac1 (*fac1, dwords)
EndProcedure

Procedure fpadd (*result.decfloat, *x.decfloat, *y.decfloat, dwords_in.l=#num_dwords)
   Define.l dwords, i, t, c, xsign, ysign
   Define.decfloat fac1, fac2
   dwords = dwords_in
   If dwords > #num_dwords
      dwords = #num_dwords
   EndIf
   
   xsign = *x\sign: *x\sign = 0
   ysign = *y\sign: *y\sign = 0
   c = fpcmp(*x, *y, dwords_in)
   
   *x\sign = xsign
   *y\sign = ysign
   If c < 0
      copydec(*y, fac1) ;fac1 = y
      copydec(*x, fac2) ;fac2 = x
   Else
      copydec(*x, fac1) ;fac1 = x
      copydec(*y, fac2) ;fac2 = y
   EndIf
   t = fac1\exponent - fac2\exponent
   t = ((fac1\exponent & $7fffffff) - #bias - 1) - ((fac2\exponent & $7fffffff) - #bias - 1)
   
   If t < (dwords*8 + 8)
      ;the difference between the two
      ;exponents indicate how many times
      ;we have to multiply the mantissa
      ;of fac2 by 10 (i.e., shift it right 1 place).
      ;if we have to shift more times than
      ;we have digits, the result is already in fac1.
      t = fac1\exponent - fac2\exponent
      If t > 0 And t < (dwords*8 + 8) ;shift
         
         i = t / 8
         While i > 0
            rshift_8 (fac2, dwords)
            t = t - 8
            i = i - 1
         Wend
         If t = 7
            rshift_7 (fac2, dwords)
         ElseIf t = 6
            rshift_6 (fac2, dwords)
         ElseIf t = 5
            rshift_5 (fac2, dwords)
         ElseIf t = 4
            rshift_4 (fac2, dwords)
         ElseIf t = 3
            rshift_3 (fac2, dwords)
         ElseIf t = 2
            rshift_2 (fac2, dwords)
         ElseIf t = 1
            rshift_1 (fac2, dwords)
         EndIf
      EndIf
      ;see if the signs of the two numbers
      ;are the same.  if so, add; if not, subtract.
      If fac1\sign = fac2\sign ;add
         fpadd_aux (fac1, fac2, dwords)
      Else
         fpsub_aux (fac1, fac2, dwords)
      EndIf
   EndIf
   copydec(fac1, *result) ;result = fac1
EndProcedure

Procedure fpsub (*result.decfloat, *x.decfloat, *y.decfloat, dwords_in.l=#num_dwords)
   Define.l dwords
   Define.decfloat fac1, fac2
   dwords = dwords_in
   If dwords > #num_dwords
      dwords = #num_dwords
   EndIf
   copydec(*x, fac1) ;fac1 = x
   copydec(*y, fac2) ;fac2 = y
   fac2\sign = fac2\sign ! $8000
   fpadd (*result, fac1, fac2, dwords)
EndProcedure

Procedure fpmul (*result.decfloat, *x.decfloat, *y.decfloat, dwords_in.l=#num_dwords)
   Define.l dwords, i, j, ex, er, ndw, den, num
   Define.q digit, carry
   Define.decfloat fac1, fac2
   dwords = dwords_in
   If dwords > #num_dwords
      dwords = #num_dwords
   EndIf
   
   Dim fac3.q(dwords)
   copydec(*x, fac1) ;fac1 = x
   copydec(*y, fac2) ;fac2 = y
                     ;check exponents.  if either is zero,
                     ;the result is zero
   If fac1\exponent = 0 Or fac2\exponent = 0 ;result is zero...clear fac1.
      fac1\sign = 0
      fac1\exponent = 0
      For i = 0 To dwords
         fac1\mantissa[i] = 0
      Next
      copydec(fac1, *result) ;result = fac1
      ProcedureReturn
   Else
      
      If ex < 0
         er = #expo_err
         copydec(fac1, *result) ;result = fac1
         ProcedureReturn
      EndIf
      
      ;clear fac3 mantissa
      For i = 0 To dwords
         fac3(i) = 0
      Next
      
      den = dwords
      While fac2\mantissa[den] = 0
         den = den - 1
      Wend
      num = dwords
      While fac1\mantissa[num] = 0
         num = num - 1
      Wend
      
      If num < den
         ;Swap fac1, fac2
         copydec(*y, fac1) ;;fac1=y
         copydec(*x, fac2) ;;fac2=x
         Swap den, num
      EndIf
      
      For j = den To 0 Step -1
         carry = 0
         digit = fac2\mantissa[j]
         For i = num To 0 Step -1
            fac3(i) = fac3(i) + digit * fac1\mantissa[i]
         Next
         For i = num To 0 Step -1
            fac3(i) = fac3(i) + carry
            carry = fac3(i) / 100000000
            fac3(i) = (fac3(i) % 100000000)
         Next
         
         For i = dwords To 1 Step -1
            fac3(i) = fac3(i - 1)
         Next
         fac3(0) = carry
      Next
      
      For i = 0 To dwords
         ndw = fac3(i)
         fac1\mantissa[i] = ndw
      Next
   EndIf
   ;now determine exponent of result.
   ;as you do...watch for overflow.
   ex = fac2\exponent - #bias + fac1\exponent
   fac1\exponent = ex
   ;determine the sign of the product
   fac1\sign = fac1\sign ! fac2\sign
   norm_fac1 (fac1, dwords)
   copydec(fac1, *result) ;result = fac1
EndProcedure

Procedure fpmul_si (*result.decfloat, *x.decfloat, y.q, dwords_in.l=#num_dwords)
   Define.l dwords, ex, er, i, ndw, count
   Define.decfloat fac1, fac2
   Define.q carry, digit, prod, value
   dwords = dwords_in
   If dwords > #num_dwords
      dwords = #num_dwords
   EndIf
   
   copydec(*x, fac1) ;fac1 = x
   digit = Abs(y)
   If digit > 99999999
      si2fp (fac2, y, dwords)
      fpmul (fac1, fac1, fac2, dwords)
      copydec(fac1, *result) ;result = fac1
      ProcedureReturn
   EndIf
   ;check exponents.  if either is zero,
   ;the result is zero
   If fac1\exponent = 0 Or y = 0 ;result is zero...clear fac1.
      fac1\sign = 0
      fac1\exponent = 0
      For count = 0 To dwords
         fac1\mantissa[count] = 0
      Next
      norm_fac1 (fac1, dwords)
      copydec(fac1, *result) ;result = fac1
      ProcedureReturn
   Else
      If digit = 1
         If y < 0
            fac1\sign = fac1\sign ! $8000
         EndIf
         copydec(fac1, *result) ;result = fac1
         ProcedureReturn
      EndIf
      ;now determine exponent of result.
      ;as you do...watch for overflow.
      
      If ex < 0
         er = #expo_err
         copydec(fac1, *result) ;result = fac1
         ProcedureReturn
      EndIf
      
      carry = 0
      
      For i = dwords To 0 Step -1
         prod = digit * fac1\mantissa[i] + carry
         value = (prod % 100000000)
         ndw = value
         fac1\mantissa[i] = ndw
         carry = prod / 100000000
      Next
      
      If carry < 10
         rshift_1 (fac1, dwords)
         fac1\exponent = fac1\exponent + 1
         fac1\mantissa[0] = fac1\mantissa[0] + carry * 10000000
      ElseIf carry < 100
         rshift_2 (fac1, dwords)
         fac1\exponent = fac1\exponent + 2
         fac1\mantissa[0] = fac1\mantissa[0] + carry * 1000000
      ElseIf carry < 1000
         rshift_3 (fac1, dwords)
         fac1\exponent = fac1\exponent + 3
         fac1\mantissa[0] = fac1\mantissa[0] + carry * 100000
      ElseIf carry < 10000
         rshift_4 (fac1, dwords)
         fac1\exponent = fac1\exponent + 4
         fac1\mantissa[0] = fac1\mantissa[0] + carry * 10000
      ElseIf carry < 100000
         rshift_5 (fac1, dwords)
         fac1\exponent = fac1\exponent + 5
         fac1\mantissa[0] = fac1\mantissa[0] + carry * 1000
      ElseIf carry < 1000000
         rshift_6 (fac1, dwords)
         fac1\exponent = fac1\exponent + 6
         fac1\mantissa[0] = fac1\mantissa[0] + carry * 100
      ElseIf carry < 10000000
         rshift_7 (fac1, dwords)
         fac1\exponent = fac1\exponent + 7
         fac1\mantissa[0] = fac1\mantissa[0] + carry * 10
      ElseIf carry < 100000000
         rshift_8 (fac1, dwords)
         fac1\exponent = fac1\exponent + 8
         fac1\mantissa[0] = fac1\mantissa[0] + carry
      EndIf
      
   EndIf
   
   norm_fac1 (fac1, dwords)
   
   If y < 0
      fac1\sign = fac1\sign ! $8000
   EndIf
   copydec(fac1, *result) ;result = fac1
EndProcedure

Procedure fpmul_ll (*result.decfloat, *x.decfloat, y.q, dwords_in.l=#num_dwords)
   Define.l dwords, count, ex, er, i, ndw
   Define.decfloat fac1, fac2
   Define.q n0, n1, carry, digit, prod, value
   dwords = dwords_in
   If dwords > #num_dwords
      dwords = #num_dwords
   EndIf
   
   copydec(*x, fac1) ;fac1 = x
   digit = Abs(y)
   If digit > 99999999 And digit < 10000000000000000
      n0 = digit / 100000000
      n1 = digit % 100000000
      fpmul_si (fac2, fac1, n1, dwords)
      fac1\exponent = fac1\exponent + 8
      fpmul_si (fac1, fac1, n0, dwords)
      fpadd (fac1, fac1, fac2, dwords)
      If y < 0
         fac1\sign = fac1\sign ! $8000
      EndIf
      copydec(fac1, *result) ;result = fac1
      ProcedureReturn
   EndIf
   If digit > 9999999999999999
      str2fp (fac2, Str(y))
      fpmul (fac1, fac1, fac2, dwords)
      copydec(fac1, *result) ;result = fac1
      ProcedureReturn
   EndIf
   si2fp (fac2, y, dwords)
   
   ;check exponents.  if either is zero,
   ;the result is zero
   If fac1\exponent = 0 Or y = 0 ;result is zero...clear fac1.
      fac1\sign = 0
      fac1\exponent = 0
      For count = 0 To dwords
         fac1\mantissa[count] = 0
      Next
      norm_fac1 (fac1, dwords)
      copydec(fac1, *result) ;result = fac1
      ProcedureReturn
   Else
      If digit = 1
         If y < 0
            fac1\sign = fac1\sign ! $8000
         EndIf
         copydec(fac1, *result) ;result = fac1
         ProcedureReturn
      EndIf
      ;now determine exponent of result.
      ;as you do...watch for overflow.
      
      If ex < 0
         er = #expo_err
         copydec(fac1, *result) ;result = fac1
         ProcedureReturn
      EndIf
      
      carry = 0
      
      For i = dwords To 0 Step -1
         prod = digit * fac1\mantissa[i] + carry
         value = (prod % 100000000)
         ndw = value
         fac1\mantissa[i] = ndw
         carry = prod / 100000000
      Next
      
      If carry < 10
         rshift_1 (fac1, dwords)
         fac1\exponent = fac1\exponent + 1
         fac1\mantissa[0] = fac1\mantissa[0] + carry * 10000000
      ElseIf carry < 100
         rshift_2 (fac1, dwords)
         fac1\exponent = fac1\exponent + 2
         fac1\mantissa[0] = fac1\mantissa[0] + carry * 1000000
      ElseIf carry < 1000
         rshift_3 (fac1, dwords)
         fac1\exponent = fac1\exponent + 3
         fac1\mantissa[0] = fac1\mantissa[0] + carry * 100000
      ElseIf carry < 10000
         rshift_4 (fac1, dwords)
         fac1\exponent = fac1\exponent + 4
         fac1\mantissa[0] = fac1\mantissa[0] + carry * 10000
      ElseIf carry < 100000
         rshift_5 (fac1, dwords)
         fac1\exponent = fac1\exponent + 5
         fac1\mantissa[0] = fac1\mantissa[0] + carry * 1000
      ElseIf carry < 1000000
         rshift_6 (fac1, dwords)
         fac1\exponent = fac1\exponent + 6
         fac1\mantissa[0] = fac1\mantissa[0] + carry * 100
      ElseIf carry < 10000000
         rshift_7 (fac1, dwords)
         fac1\exponent = fac1\exponent + 7
         fac1\mantissa[0] = fac1\mantissa[0] + carry * 10
      ElseIf carry < 100000000
         rshift_8 (fac1, dwords)
         fac1\exponent = fac1\exponent + 8
         fac1\mantissa[0] = fac1\mantissa[0] + carry
      EndIf
      
   EndIf
   
   norm_fac1 (fac1, dwords)
   
   If y < 0
      fac1\sign = fac1\sign ! $8000
   EndIf
   copydec(fac1, *result) ;result = fac1
EndProcedure

Procedure fpreciprocal (*result.decfloat, *m.decfloat, dwords_in.l=#num_dwords)
   Define.l dwords
   Define.d x
   dwords=dwords_in
   If dwords>#NUM_DWORDS
      dwords=#NUM_DWORDS
   EndIf
   Define.l k, l, i, ex, s, prec
   Define.decfloat r, r2, one, n
   copydec(*m, n) ;n=m
   l=Log((dwords*8+8)*0.0625)*1.5
   Define.s v,f, ts
   If n\exponent>0
      ex=(n\exponent & $7FFFFFFF)-#BIAS-1
   Else
      ex=0
   EndIf
   
   If n\sign
      v="-"
   Else
      v=" "
   EndIf
   ts=Str(n\mantissa[0])
   If Len(ts)<8
      ts=ts+string(8-Len(ts),"0")
   EndIf
   v=v+Left(ts,1)+"."+Mid(ts,2)
   ts=Str(n\mantissa[1])
   If Len(ts)<8
      ts=string(8-Len(ts),"0")+ts
   EndIf		
   v=v+ts
   x=ValD(v)
   If x = 0
      PrintN( "Div 0")
      copydec(r, *result)
      ProcedureReturn
   EndIf
   If x = 1 And ex=0
      str2fp(r,"1")
      copydec(r, *result)
      ProcedureReturn
   ElseIf x=1
      x=10
      ex=ex-1
   EndIf
   ex=(-1)-ex
   x=1/x
   str2fp(r,StrD(x,16)) 
   r\exponent=ex+#BIAS+1
   one\mantissa[0]=10000000
   one\exponent=#BIAS+1
   copydec(r, r2) ;r2=r
   prec=3
   For k=1 To l
      prec=2*prec-1
      fpmul(r2,n,r, prec)
      fpsub(r2,one, r2, prec)
      fpmul(r2,r, r2, prec)
      fpadd(r,r,r2, prec)
   Next
   copydec(r, *result)
EndProcedure

Procedure.l min (a.l, b.l)
   If a < b
      ProcedureReturn a
   Else
      ProcedureReturn b
   EndIf
EndProcedure

Procedure.d realw (Array w.d(1), j.l)
   Define.d wx
   wx = ((w(j - 1) * 10000 + w(j)) * 10000 + w(j + 1)) * 10000
   If ArraySize(w()) >= (j + 2)
      wx = wx + w(j + 2)
   EndIf
   ProcedureReturn wx
EndProcedure

Procedure subtract (Array w.d(1), q.l, Array d.d(1), ka.l, kb.l)
   Define.l j
   For j = ka To kb
      w(j) = w(j) - q * d(j - ka + 2)
   Next
EndProcedure

Procedure normalise (Array w.d(1), ka.l, q.l)
   w(ka) = w(ka) + w(ka - 1) * 10000
   w(ka - 1) = q
EndProcedure

Procedure finalnorm (Array w.d(1), kb.l)
   Define.l carry, j
   For j = kb To 3 Step -1
      If w(j) < 0
         carry = ((-Int(w(j)) - 1) / 10000) + 1
      Else
         If w(j) >= 10000
            carry = -(Int(w(j)) / 10000)
         Else
            carry = 0
         EndIf
      EndIf
      w(j) = w(j) + carry * 10000
      w(j - 1) = w(j - 1) - carry
   Next
EndProcedure

Procedure fpdiv (*result.decfloat, *x.decfloat, *y.decfloat, dwords_in.l=#num_dwords)
   Define.l dwords, i, er, is_power_of_ten, tmp, b, j, stp, last, laststep, q, t
   Define.d xd, xn, rund
   Define.decfloat fac1, fac2
   dwords = dwords_in
   If dwords > #num_dwords
      dwords = #num_dwords
   EndIf
   
   copydec(*x, fac1) ;fac1 = x
   copydec(*y, fac2) ;fac2 = y
   
   If fac2\exponent = 0  ; if fac2 = 0, return
                         ; a divide-by-zero error and
                         ; bail out.
      For i = 0 To dwords
         fac1\mantissa[i] = 99999999
      Next
      fac1\exponent = 99999 + #bias + 1
      er = #divz_err
      copydec(fac1, *result) ;result = fac1
      ProcedureReturn
   ElseIf fac1\exponent = 0 ;fact1=0, just return
      er = 0
      copydec(fac1, *result) ;result = fac1
      ProcedureReturn
   Else
      ;check to see if fac2 is a power of ten
      is_power_of_ten = 0
      If fac2\mantissa[0] = 10000000
         is_power_of_ten = 1
         For i = 1 To dwords
            If fac2\mantissa[i] <> 0
               is_power_of_ten = 0
               Break
            EndIf
         Next
      EndIf
      ;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
         fac1\sign = fac1\sign ! fac2\sign
         fac1\exponent = fac1\exponent - fac2\exponent + #bias + 1
         copydec(fac1, *result) ;result = fac1
         ProcedureReturn
      EndIf
      
      Dim result.d(2 * dwords + 3)
      Dim n.d(2 * dwords + 3)
      Dim d.d(2 * dwords + 3)
      b= 10000
      Dim w.d(ArraySize(n()) + 4)
      
      For j = 0 To dwords
         tmp=fac1\mantissa[j] / 10000
         n(2 * j + 2) = tmp
         tmp=fac1\mantissa[j] % 10000
         n(2 * j + 3) = tmp
         tmp=fac2\mantissa[j] / 10000
         d(2 * j + 2) = tmp
         tmp=fac2\mantissa[j] % 10000
         d(2 * j + 3) = tmp
      Next
      tmp=(fac1\exponent & $7fffffff) - #bias - 1
      n(1) = tmp
      tmp=(fac2\exponent & $7fffffff) - #bias - 1
      d(1) = tmp
      For j = ArraySize(n()) To ArraySize(w())
         w(j) = 0
      Next
      t = ArraySize(n()) - 1
      w(1) = n(1) - d(1) + 1
      w(2) = 0
      For j = 2 To ArraySize(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, ArraySize(w()))
         subtract (w(), q, d(), (stp + 2), last)
         normalise (w(), (stp + 2), q)
      Next
      finalnorm (w(), (laststep + 1))
      If w(2) <> 0
         laststep = laststep - 1
      EndIf
      rund = w(laststep + 1) / b
      If rund >= 0.5
         w(laststep) = w(laststep) + 1
      EndIf
      If w(2) = 0
         For j = 1 To t + 1
            result(j) = w(j + 1)
         Next
      Else
         For j = 1 To t + 1
            result(j) = w(j)
         Next
      EndIf
      If w(2) = 0
         result(1) = w(1) - 1
      Else
         result(1) = w(1)
      EndIf
      For j = 0 To dwords
         fac1\mantissa[j] = result(2 * j + 2) * 10000 + result(2 * j + 3)
      Next
      norm_fac1 (fac1, dwords)
      fac1\exponent = (result(1) + #bias)
   EndIf
   fac1\sign = fac1\sign ! fac2\sign
   copydec(fac1, *result) ;result = fac1
EndProcedure

; sqrt(num)
; sqrt(num)
Procedure fpsqr (*result.decfloat, *num.decfloat, dwords_in.l=#num_dwords)
   Define.l dwords, ex, aex, k, l, prec
   Define.decfloat r, r2, tmp, n, half
   Define.s ts, v
   Define.d x
   dwords = dwords_in
   If dwords > #num_dwords
      dwords = #num_dwords
   EndIf
   l = Log((dwords*8 + 8) * 0.0625) * 1.5
   ;l=estimated number of iterations needed
   ;first estimate is accurate to about 16 digits
   ;l is approximatly = to log2((#num_digits+9)/16)
   ;#num_digits+9 because decfloat has an extra 9 guard digits
   copydec(*num, n) ;n = num
   si2fp (tmp, 0, dwords)
   If fpcmp(n, tmp, dwords) = 0
      si2fp (r, 0, dwords)
      copydec(r, *result) ;result = r
      ProcedureReturn
   EndIf
   si2fp (tmp, 1, dwords)
   If fpcmp(n, tmp, dwords) = 0
      si2fp (r, 1, dwords)
      copydec(r, *result) ;result = r
      ProcedureReturn
   EndIf
   si2fp (tmp, 0, dwords)
   If fpcmp(n, tmp, dwords) < 0
      si2fp (r, 0, dwords)
      copydec(r, *result) ;result = r
      ProcedureReturn
   EndIf
   ;=====================================================================
   ;hack to bypass the limitation of double exponent range
   ;in case the number is larger than what a double can handle
   ;for example, if the number is 2e500
   ;we separate the exponent and mantissa in this case 2
   ;if the exponent is odd then multiply the mantissa by 10
   ;take the square root and assign it to decfloat
   ;divide the exponent in half for square root
   ;in this case 1.414213562373095e250
   If n\exponent > 0
      ex = (n\exponent & $7fffffff) - #bias - 1
   Else
      ex = 0
   EndIf
   ts = Trim(Str(n\mantissa[0]))
   If Len(ts) < 8
      ts = ts + String(8 - Len(ts), "0")
   EndIf
   v = v + Left(ts, 1) + "." + Mid(ts, 2)
   ts = Trim(Str(n\mantissa[1]))
   If Len(ts) < 8
      ts = String(8 - Len(ts), "0") + ts
   EndIf
   v = v + ts
   x = ValD(v)
   
   If x = 0
      PrintN( "div 0")
      copydec(r, *result) ;result = r
      ProcedureReturn
   EndIf
   If x = 1 And ex = 0
      si2fp (r, 1, dwords)
      copydec(r, *result) ;result = r
      ProcedureReturn
   EndIf
   aex=ex
   If aex<0
      aex=-aex
   EndIf
   If aex & 1
      x = x * 10
      ex = ex - 1
   EndIf
   x = Sqr(x) ;approximation
   v = Trim(StrD(x,16))
   k = FindString(v, ".")
   str2fp (r, v)
   r\exponent = ex / 2 + #bias + 1
   If Len(v) > 1 And k = 0
      r\exponent = r\exponent + 1
   EndIf
   For k = 1 To dwords
      half\mantissa[k] = 0
   Next
   half\mantissa[0] = 50000000
   half\exponent = #bias
   half\sign = 0
   ;=====================================================================
   ;newton-raphson method
   prec = 3
   
   For k = 1 To l + 1
      prec = 2 * prec - 1
      fpdiv (tmp, n, r, prec)
      fpadd (r2, r, tmp, prec)
      fpmul (r, r2, half, prec)
   Next
   copydec(r, *result) ;result = r
EndProcedure

Procedure fpdiv_si (*result.decfloat, *num.decfloat, den.l, dwords_in.l=#num_dwords)
   Define.l dwords, ndw, i
   Define.q carry, remder, divisor, quotient
   Define.decfloat fac1
   dwords = dwords_in
   If dwords > #num_dwords
      dwords = #num_dwords
   EndIf
   remder = 0
   divisor = Abs(den)
   copydec(*num, fac1) ;fac1 = num
   If divisor = 0
      PrintN( "error: divisor = 0")
      copydec(fac1, *result) ;result = fac1
      ProcedureReturn
   EndIf
   If divisor > 99999999
      PrintN( "error: divisor too large")
      copydec(fac1, *result) ;result = fac1
      ProcedureReturn
   EndIf
   
   For i = 0 To dwords
      ndw = i
      quotient = fac1\mantissa[ndw] + remder * 100000000
      remder = quotient % divisor
      fac1\mantissa[ndw] = quotient / divisor
   Next
   quotient = remder * 100000000
   quotient = quotient / divisor
   carry = fac1\mantissa[0]
   
   If carry = 0
      lshift_8 (fac1, dwords)
      fac1\exponent = fac1\exponent - 8
      fac1\mantissa[dwords] = fac1\mantissa[dwords] + quotient
   ElseIf carry < 10
      lshift_7 (fac1, dwords)
      fac1\exponent = fac1\exponent - 7
      fac1\mantissa[dwords] = fac1\mantissa[dwords] + quotient / 10
   ElseIf carry < 100
      lshift_6 (fac1, dwords)
      fac1\exponent = fac1\exponent - 6
      fac1\mantissa[dwords] = fac1\mantissa[dwords] + quotient / 100
   ElseIf carry < 1000
      lshift_5 (fac1, dwords)
      fac1\exponent = fac1\exponent - 5
      fac1\mantissa[dwords] = fac1\mantissa[dwords] + quotient / 1000
   ElseIf carry < 10000
      lshift_4 (fac1, dwords)
      fac1\exponent = fac1\exponent - 4
      fac1\mantissa[dwords] = fac1\mantissa[dwords] + quotient / 10000
   ElseIf carry < 100000
      lshift_3 (fac1, dwords)
      fac1\exponent = fac1\exponent - 3
      fac1\mantissa[dwords] = fac1\mantissa[dwords] + quotient / 100000
   ElseIf carry < 1000000
      lshift_2 (fac1, dwords)
      fac1\exponent = fac1\exponent - 2
      fac1\mantissa[dwords] = fac1\mantissa[dwords] + quotient / 1000000
   ElseIf carry < 10000000
      lshift_1 (fac1, dwords)
      fac1\exponent = fac1\exponent - 1
      fac1\mantissa[dwords] = fac1\mantissa[dwords] + quotient / 10000000
   EndIf
   
   ;norm_fac1(fac1)
   If den < 0
      fac1\sign = fac1\sign ! $8000
   EndIf
   copydec(fac1, *result) ;result = fac1
EndProcedure

;fractional part of num
Procedure fpfrac (*result.decfloat, *num.decfloat, dwords_in.l=#num_dwords)
   Define.l dwords
   Define.decfloat n
   dwords = dwords_in
   If dwords > #num_dwords
      dwords = #num_dwords
   EndIf
   fptrunc (n, *num)
   fpsub (n, *num, n, dwords)
   copydec(n, *result) ;result = n
EndProcedure

;returns the positive of n
Procedure fpabs (*result.decfloat, *n.decfloat)
   Define.decfloat x
   copydec(*n, x) ;x = n
   x\sign = 0
   copydec(x, *result) ;result = x
EndProcedure

;changes the sign of n, if n is positive then n will be negative & vice versa
Procedure fpneg (*result.decfloat, *n.decfloat)
   Define.decfloat x
   copydec(*n, x) ;x = n
   x\sign = x\sign ! $8000
   copydec(x, *result) ;result = x
EndProcedure

Procedure fpfmod (*quotient.decfloat, *f_mod.decfloat, *num.decfloat, *denom.decfloat, dwords_in.l=#num_dwords)
   Define.l dwords
   Define.decfloat q, fm
   If dwords > #num_dwords
      dwords = #num_dwords
   EndIf
   fpdiv (fm, *num, *denom, dwords)
   fptrunc (q, fm)
   fpmul(q, q, *denom, dwords)
   fpsub (fm, *num, q, dwords)
   copydec(q, *quotient) ;quotient = q
   copydec(fm, *f_mod)   ;f_mod = fm
EndProcedure

Procedure fpeps (*result.decfloat, digits_in.l=#number_of_digits)
   Define.l digits, dwords
   Define.decfloat ep
   digits=digits_in
   If digits > #number_of_digits
      digits = #number_of_digits
   EndIf
   dwords=digits/8
   si2fp (ep, 5, dwords)
   ep\exponent = (-(digits) + #bias + 1)
   copydec(ep, *result) ;result = ep
EndProcedure

Procedure fpipow (*result.decfloat, *x.decfloat, e.q, dwords_in.l=#num_dwords)
   Define.l dwords, i
   Define.decfloat y, z
   Define.q n, c
   dwords = dwords_in
   If dwords > #num_dwords
      dwords = #num_dwords
   EndIf
   ;take x to an long power
   
   c = 0
   copydec(*x, y) ;y = x
   n = Abs(e)
   z\sign = 0
   z\exponent = (#bias + 1)
   z\mantissa[0] = 10000000
   For i = 1 To #num_dwords
      z\mantissa[i] = 0
   Next
   While n > 0
      While (n & 1) = 0
         n = n / 2
         fpmul (y, y, y, dwords)
         c = c + 1
      Wend
      n = n - 1
      fpmul (z, y, z, dwords)
      c = c + 1
   Wend
   If e < 0
      si2fp (y, 1, dwords)
      fpdiv (z, y, z, dwords)
   EndIf
   copydec(z, *result) ;result = z
EndProcedure

Procedure fpfactorial (*result.decfloat, n.l, dwords_in.l=#num_dwords)
   Define.l dwords
   Define.q i
   Define.decfloat f
   dwords = dwords_in
   If dwords > #num_dwords
      dwords = #num_dwords
   EndIf
   If n < 0
      PrintN( "inf")
      copydec(f, *result) ;result = f
      ProcedureReturn
   EndIf
   If n = 0 Or n = 1
      si2fp (f, 1, dwords)
      copydec(f, *result) ;result = f
      ProcedureReturn
   EndIf
   si2fp (f, 2, dwords_in)
   If n = 2
      copydec(f, *result) ;result = f
      ProcedureReturn
   EndIf
   For i = 3 To n
      fpmul_si (f, f, i, dwords)
   Next
   copydec(f, *result) ;result = f
EndProcedure

Procedure fpnroot (*result.decfloat, *x.decfloat, p.l, dwords_in.l=#num_dwords)
   Define.l dwords, i, ex, l, prec
   Define.d t, t2
   Define.decfloat ry, tmp, tmp2
   dwords = dwords_in
   If dwords > #num_dwords
      dwords = #num_dwords
   EndIf
   
   l = Log((dwords*8 + 8) * 0.0625) * 1.5
   ex = (*x\exponent & $7fffffff) - #bias - 1
   t = *x\mantissa[0] + *x\mantissa[1] / 100000000
   t = t / 10000000
   t2 = Log(t) / p + Log(10) * ex / p
   t2 = Exp(t2)
   str2fp (ry, StrD(t2,16))
   prec = 3
   
   fpipow (tmp, ry, p - 1, prec)
   fpdiv (tmp, *x, tmp, prec)
   fpmul_si (tmp2, ry, p - 1, prec)
   fpadd (tmp2, tmp2, tmp, prec)
   fpdiv_si (ry, tmp2, p, prec)
   For i = 1 To l
      prec = 2 * prec - 1
      fpipow (tmp, ry, p - 1, prec)
      fpdiv (tmp, *x, tmp, prec)
      fpmul_si (tmp2, ry, p - 1, prec)
      fpadd (tmp2, tmp2, tmp, prec)
      fpdiv_si (ry, tmp2, p, prec)
   Next
   copydec(ry, *result) ;result = ry
EndProcedure

Procedure pi_brent_salamin (*pi_bs.decfloat, digits_in.l= #number_of_digits)
   Define.l dwords, digits, limit2
   Define.decfloat c0, c1, c2, c05, a, b, sum, ak, bk, ck, ab, asq
   Define.decfloat pow2, tmp
   digits=digits_in
   If digits > #number_of_digits
      digits = #number_of_digits
   EndIf
   dwords = digits/8
   If (dwords*8)<digits
      dwords+1
   EndIf
   
   limit2 = -(digits) + #bias + 1
   
   si2fp (c0, 0, dwords)
   copydec(c0, ak) ;ak = c0
   copydec(c0, bk) ;bk = c0
   copydec(c0, ab) ;ab = c0
   copydec(c0, asq);asq = c0
   si2fp (c1, 1, dwords)
   copydec(c1, a) ;a = c1
   copydec(c1, ck);ck = c1
   copydec(c1, pow2) ;pow2 = c1
   si2fp (c2, 2, dwords)
   copydec(c2, b) ;b = c2
   str2fp (c05, ".5")
   copydec(c05, sum) ;sum = c05
   si2fp (*pi_bs, 3, dwords)
   
   fpsqr (b, b, dwords)
   fpdiv (b, c1, b, dwords)
   While fpcmp(ck, c0, dwords) <> 0 And ck\exponent > limit2
      fpadd (ak, a, b, dwords)
      fpmul (ak, c05, ak, dwords)
      fpmul (ab, a, b, dwords)
      fpsqr (bk, ab, dwords)
      fpmul (asq, ak, ak, dwords)
      fpsub (ck, asq, ab, dwords)
      fpmul (pow2, pow2, c2, dwords)
      fpmul (tmp, pow2, ck, dwords)
      fpsub (sum, sum, tmp, dwords)
      copydec(ak, a) ;a = ak
      copydec(bk, b) ;b = bk
   Wend
   fpdiv (tmp, asq, sum, dwords)
   fpmul (*pi_bs, c2, tmp, dwords)
EndProcedure

decfloat-test.pb

Code: Select all

#number_of_digits = 100 ; need to set the constant before including decfloat.pbi

XIncludeFile "decfloat.pbi"

OpenConsole()
Define.decfloat m, x, y, z, v
Define.s s
Define.d t

str2fp (x, "3."+String(#number_of_digits-1, "3"))
str2fp (y, "4."+String(#number_of_digits-1, "4"))

t=ElapsedMilliseconds()

PrintN("x        = "+fp2str(x))
PrintN("y        = "+fp2str(y))
fpadd(z, x, y)
PrintN("x + y    = "+fp2str(z))
;    c = fpcmp(x, y)
;    PrintN("c="+Str(c))
fpsub(z, x, y)
PrintN("x - y    = "+fp2str(z))

fpmul(z, x, y)
PrintN("x * y    = "+fp2str(z))

fpdiv(z, x, y)
PrintN("x / y    = "+fp2str(z))

si2fp (x, 2)
PrintN("x        = "+fp2str(x))
fpsqr (z, x)
PrintN("sqrt(x)  = "+fp2str(z))

fpipow (z, x, -3)
PrintN("x^(-3)   = "+fp2str(z))

fpnroot (z, x, 3)
PrintN("x^(1/3)  = "+fp2str(z))

fpfactorial (z, 10000)
PrintN("10000!   = "+fp2str(z))

pi_brent_salamin (z)
PrintN("Pi       = "+fp2str(z))

fpfrac(v, z)
PrintN("frac(Pi) = "+fp2str(v))

fptrunc(z, z)
PrintN("trunc(Pi)= "+fp2str(z))

si2fp (x, 1)
;PrintN("x        = "+fp2str(x))
str2fp (y, ".18")
PrintN("y        = "+fp2str(y))
;fpdiv (z, x, y)
fpreciprocal(z, y)
PrintN("reciprocal(x)="+fp2str(z))

t=ElapsedMilliseconds()-t
PrintN( "elapsed time "+StrD(t/1000)+" seconds")
Input()
CloseConsole()

log, exp and trig functions are posted below https://www.purebasic.fr/english/viewto ... 77#p585577
output

x = 3.3333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333
y = 4.4444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444
x + y = 7.7777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777
x - y = -1.1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111
x * y = 14.814814814814814814814814814814814814814814814814814814814814814814814814814814814814814814814814
x / y = 0.75000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
x = 2.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
sqrt(x) = 1.4142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415
x^(-3) = 0.12500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
x^(1/3) = 1.2599210498948731647672106072782283505702514647015079800819751121552996765139594837293965624362550
10000! = 2.84625968091705451890641321211986889014805140170279923079417999427441134000376444377299078675778477e+35659
Pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170
frac(Pi) = 0.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706
trunc(Pi)= 3.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
y = 0.18000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
reciprocal(y)= 5.5555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555
elapsed time 0.004 seconds

corrected a bug with str2fp
corrected fpmod as sugested by STARGÅTE, thank you STARGÅTE
exponent increased to 999999999
added the Procedure.s fp2strT (*n.decfloat, places) for trimmed output - no trailing 0's

[netmaestro]

In procedure fpdiv_si I get the following error: For i = 0 To dwords 'For' doesn't support quad variables.

PB 6 beta 10 x86 and x64

[jack]

thanks for reporting netmaestro, I will fix it right away

[davido]

@jack,

Looks great.
Thank you for sharing.

[jack]

thanks for the feedback
I edited the first post to fix the problem with q loop
there's a limitation when using arrays, I suspect it may be due to arrays using the stack?
anyway it works reasonably well for precision of 1000 digits or so

[STARGÅTE]

I played around with it a bit and here are some issues:
I receive errors when the exponent is lager than 99999 because of your code line:

Code: Select all

f = String(5 - Len(f), "0") + f

If I use Mod(1.0, 0.3) the result is 0.0999999900000000000 instead of 0.1 ...

Code: Select all

OpenConsole()
str2fp (x, "1.0")
str2fp (y, "0.3")
fpfmod(v, z, x, y)
PrintN(fp2str(z))
Input()

The problem is, you calculate: mod = (a/b - trunc(a/b))*b , which result in lost digits.
It is better to use: mod = a - trunc(a/b)*b.

Nevertheless, nice implementation. However, please keep in mind, the performance is not in focus.
Due to the high-level and the 10^8 limb base (instead of full 32/64 bit) implementation, you lose speed.

[jack]

thank you STARGÅTE
fixed as per your suggestion

[jack]

here are fplog an fpexp and trig functions

Code: Select all

Procedure fplogTaylor(*result, *x.decfloat, dwords_in.l=#num_dwords)
   Define.l dwords
   dwords=dwords_in
   If dwords>#num_dwords
      dwords=#num_dwords
   EndIf
   ; taylor series

   Define.decfloat g, zero
   Define.l i
   fpabs(g, *x)
   zero\sign=0
   zero\exponent=0
   For i=0 To dwords
      zero\mantissa[i]=0
   Next
   If fpcmp(g, *x, dwords)<>0
      copydec(zero, *result)
      ProcedureReturn ;Exit Procedure
   EndIf
   If fpcmp(*x, zero, dwords)=0
      copydec(zero, *result)
      ProcedureReturn ;Exit Procedure
   EndIf

   Define.l invflag
   Define.decfloat Inv, XX, Term, Accum, strC, _x, tmp, tmp2, tmp3
   Define.decfloat T, B, one, Q, two
   
   one\sign=0
   one\exponent=(#bias+1)
   one\mantissa[0]=10000000
   two\sign=0
   two\exponent=(#bias+1)
   two\mantissa[0]=20000000
   For i=1 To dwords
      one\mantissa[i]=0
      two\mantissa[i]=0
   Next
   copydec(*x, _x)
   If fpcmp(*x, one, dwords)<0
      invflag=1
      fpdiv(_x, one, _x, dwords)
   EndIf
   fpsub(T, _x, one, dwords)
   fpadd(B, _x, one, dwords)
   fpdiv(accum, T, B, dwords)
   fpdiv(Q, T, B, dwords)
   copydec(Q, tmp)
   fpmul(XX, Q, Q, dwords)
   Define.l c=1
   Repeat
      c=c+2
      copydec(tmp, tmp2)
      fpmul(Q, Q, XX, dwords)
      fpdiv_si(term, Q,c, dwords)
      fpadd(accum, tmp, Term, dwords)
      ; Swap tmp,Accum
      copydec(tmp, tmp3)
      copydec(Accum, tmp)
      copydec(tmp3, Accum)
   Until fpcmp(tmp, tmp2, dwords) = 0
   fpmul_si(accum, accum,2, dwords)
   If invflag
      fpneg(*result, accum)
      ProcedureReturn 
   EndIf
   copydec(Accum, *result)
EndProcedure

Procedure fplog (*result.decfloat, *x.decfloat, dwords_in.l=#num_dwords)
   Define.l dwords, i, factor
   Define.decfloat g, one, zero
   dwords=dwords_in
   If dwords>#num_dwords
      dwords=#num_dwords
   EndIf

   zero\sign=0
   zero\exponent=0
   For i=0 To dwords
      zero\mantissa[i]=0
   Next
   si2fp(one, 1, dwords)
   fpabs(g, *x)
   If fpcmp(g, *x, dwords)<>0
      copydec(zero, *result)
      ProcedureReturn ;Exit Procedure
   EndIf
   If fpcmp(*x, zero, dwords)=0
      copydec(zero, *result)
      ProcedureReturn ;Exit Procedure
   EndIf
   If fpcmp(*x, one, dwords)=0
      copydec(zero, *result)
      ProcedureReturn ;Exit Procedure
   EndIf

   Define.decfloat approx, ans, logx
   copydec(*x, approx)
   factor=4096
   fpsqr(approx, approx, dwords)
   fpsqr(approx, approx, dwords)
   fpsqr(approx, approx, dwords)
   fpsqr(approx, approx, dwords)
   fpsqr(approx, approx, dwords)
   fpsqr(approx, approx, dwords)
   fpsqr(approx, approx, dwords)
   fpsqr(approx, approx, dwords)
   fpsqr(approx, approx, dwords)
   fpsqr(approx, approx, dwords)
   fpsqr(approx, approx, dwords)
   fpsqr(approx, approx, dwords)
   fplogTaylor(logx, approx, dwords)
   fpmul_si(*result, logx, factor, dwords)
EndProcedure

Procedure fpexp (*result, *x.decfloat, dwords_in.l=#num_dwords)
   Define.l dwords
   dwords=dwords_in
   If dwords>#num_dwords
      dwords=#num_dwords
   EndIf
   ;taylor series
   Define.decfloat fac, _x, temp, p, term
   Define.l i, c
   si2fp(temp, 0)
   
   fac\sign=0
   fac\exponent=(#bias+1)
   fac\mantissa[0]=10000000
   For i=1 To dwords
      fac\mantissa[i]=0
   Next
   If fpcmp(*x, temp, dwords)=0
      copydec(fac, *result)
      ProcedureReturn
   EndIf
   c=1
   fpdiv_si(_x, *x, 8192, dwords) ;fpdiv_si(x, 67108864)
   copydec(_x, p)
   fpadd(*result, fac, _x, dwords) ;1 + x
   
   Repeat
      c+1
      copydec(*result, temp)
      fpdiv_si(fac, fac, c, dwords)
      fpmul(p, p, _x, dwords)
      fpmul(term, p,fac, dwords)
      fpadd(*result, temp,Term, dwords)
   Until fpcmp(*result,temp, dwords) = 0
   For i=1 To 13
      fpmul(*result, *result, *result, dwords)
   Next
EndProcedure

Procedure fpsin (*result, *x.decfloat, dwords_in.l=#num_dwords)
   Define.l dwords
   dwords=dwords_in
   If dwords>#num_dwords
      dwords=#num_dwords
   EndIf
   
   Define.decfloat XX, Term, Accum, p, temp2, fac, x_2
   Define.decfloat pi2, circ, Ab
   
   copydec(*x, x_2)
   fpmul_si(pi2, pi_dec, 2)
   
   copydec(pi2, circ)
   fpabs(Ab, *x)
   
   If fpcmp(Ab, circ, dwords) > 0
      ;======== CENTRALIZE ==============
      ;floor/ceil to centralize
      Define.decfloat tmp, tmp2
      fpdiv(tmp, x_2, pi2, dwords)
      fpfrac(tmp, tmp) ;frac part
      fpmul(x_2, tmp, pi2, dwords)
   EndIf
   
   Define.l lm, limit, i
   Define.decfloat factor
   lm = dwords*8
   limit = 1+Int(-0.4534499388609258 + 0.0223330028523981 * lm + 5.046181440833308E-7 * lm * lm - 4.233845303980424E-11 * lm * lm * lm)
   If limit<0
      limit=0
   EndIf
   si2fp(factor, 5, dwords)
   fpipow(factor, factor, limit, dwords)
   fpdiv(x_2, x_2, factor, dwords) ;x_=x_/5^limit

   Dim Sgn.l(3): Sgn(3) = 1
   Define.l c = 1
   
   Accum = x_2
   si2fp(fac, 1, dwords)
   copydec(x_2, p)
   fpmul(XX, x_2, x_2, dwords)
   
   Repeat
      c = c + 2
      copydec(Accum, temp2)
      fpmul_si(fac, fac, c * (c - 1), dwords)
      fpmul(p, p, XX, dwords)
      fpdiv(term, p, fac, dwords)
      If Sgn(c & 3)<>0
         fpsub(Accum, temp2, Term, dwords)
      Else
         fpadd(Accum, temp2, Term, dwords)
      EndIf
   Until fpcmp(Accum, temp2, dwords) = 0
   ;multiply the result by 5^limit
   
   For i = 1 To limit
      fpmul(p, Accum, Accum, dwords)
      fpmul(temp2, Accum, p, dwords)
      ;*** sin(5*x) = 5 * sin(x) - 20 * sin(x)^3 + 16 * sin(x)^5
      fpmul_si(Accum, Accum, 5, dwords)
      fpmul_si(Term, temp2, 20, dwords)
      fpmul_si(XX, temp2, 16, dwords)
      fpmul(XX, XX, p, dwords)
      fpsub(Accum, Accum, Term, dwords)
      fpadd(Accum, Accum, XX, dwords)
   Next i
   copydec(Accum, *result)
EndProcedure

Procedure fpcos (*result, *z.decfloat, dwords_in.l=#num_dwords)
   Define.l dwords
   dwords=dwords_in
   If dwords>#num_dwords
      dwords=#num_dwords
   EndIf
   Static.decfloat pi_dec_half
   copydec(Pi_dec, pi_dec_half)
   fpdiv_si(pi_dec_half, pi_dec_half, 2)
   Define.decfloat x_2
   
   fpsub(x_2, pi_dec_half, *z, dwords)
   fpsin(*result, x_2, dwords)
EndProcedure

Procedure fptan (*result, *z.decfloat, dwords_in.l=#num_dwords)
   Define.l dwords
   dwords=dwords_in
   If dwords>#num_dwords
      dwords=#num_dwords
   EndIf
   Define.decfloat x_2, s, c
   copydec(*z, x_2)
   fpsin(s, x_2, dwords)
   fpcos(c, x_2, dwords)
   fpdiv(*result, s, c, dwords)
EndProcedure

Procedure fpatn (*result.decfloat, *x.decfloat, dwords_in.l=#num_dwords)
   Define.l dwords
   dwords=dwords_in
   If dwords>#num_dwords
      dwords=#num_dwords
   EndIf
   Static.decfloat pi_4
   copydec(Pi_dec, pi_4)
   fpdiv_si(pi_4, pi_4, 4)
   Dim Sgn.l(3): Sgn(3) = 1
   Define.l z, c = 1
   Define.decfloat XX, Term, Accum, strC, x_2, mt, mt2, p
   Define.decfloat decnum, one, decnum2, factor
   copydec(*x, decnum2)
   decnum2\sign = 0
   si2fp(one, 1, dwords)
   If fpcmp(decnum2, one, dwords) = 0
      copydec(pi_4, *result)
      *result\sign = *x\sign
      ProcedureReturn
   EndIf
   decnum2\sign = *x\sign
   Define.l limit = 16
   si2fp(factor, 2 << (limit - 1), dwords)
   
   For z = 1 To limit
      fpmul(decnum , decnum2, decnum2, dwords)
      fpadd(decnum , decnum, one, dwords)
      fpsqr(decnum , decnum, dwords)
      fpadd(decnum , decnum, one, dwords)
      fpdiv(decnum , decnum2, decnum, dwords)
      copydec(decnum, decnum2)
   Next z
   
   copydec(decnum, mt)
   copydec(decnum, x_2)
   copydec(decnum, p)
   fpmul(XX, x_2, x_2, dwords)
   Repeat
      c = c + 2
      copydec(mt, mt2)
      si2fp(strC, c, dwords)
      fpmul(p, p, XX, dwords)
      fpdiv(Term, p, strC, dwords)
      If Sgn(c & 3)
         fpsub(Accum, mt, Term, dwords)
      Else
         fpadd(Accum, mt, Term, dwords)
      EndIf
      copydec(mt, decnum2)
      copydec(Accum, mt)
      copydec(decnum2, Accum)
   Until fpcmp(mt, mt2, dwords) = 0
   fpmul(*result, factor, mt, dwords)
EndProcedure

Procedure fpasin (*result.decfloat, *x.decfloat, dwords_in.l=#num_dwords)
   Define.l dwords
   dwords=dwords_in
   If dwords>#num_dwords
      dwords=#num_dwords
   EndIf
   Define.d num
   Define.decfloat one, T, B, term1, minusone
   ; ARCSIN = ATN(x / SQR(-x * x + 1))

   num = fp2dbl(*x)
   If (num > 1) Or (num < -1)
      copydec(one, *result)
      ProcedureReturn
   EndIf
   
   si2fp(one, 1, dwords)
   si2fp(minusone, -1, dwords)
   copydec(*x, T)
   fpmul(B, *x, *x, dwords) ;x*x
                            ;for 1 and -1
   If fpcmp(B, one, dwords) = 0
      Define.decfloat two, atn1
      si2fp(two, 2, dwords)
      fpatn(atn1, one, dwords)
      If fpcmp(*x, minusone, dwords) = 0
         fpmul(two, two, atn1, dwords)
         fpmul(*result, two, minusone, dwords)
         ProcedureReturn
      Else
         fpmul(*result, two, atn1, dwords)
         ProcedureReturn
      EndIf
   EndIf
   fpsub(B, one, B, dwords) ;1-x*x
   fpsqr(B, B, dwords)      ;sqr(1-x*x)
   fpdiv(term1, T, B, dwords)
   fpatn(*result, term1, dwords)
EndProcedure

Procedure fpacos (*result.decfloat, *x.decfloat, dwords_in.l=#num_dwords)
   Define.l dwords
   dwords=dwords_in
   If dwords>#num_dwords
      dwords=#num_dwords
   EndIf
   Define.decfloat one, minusone, two, atn1, tail, T, B, term1, atnterm1 ;,_x,temp
   Define.d num
   ;ARCCOS = ATN(-x / SQR(-x * x + 1)) + 2 * ATN(1)

   num = fp2dbl(*x)
   If (num > 1) Or (num < -1)
      copydec(one, *result)
      ProcedureReturn
   EndIf

   si2fp(one, 1, dwords)
   si2fp(minusone, -1, dwords)
   si2fp(two, 2, dwords)
   fpatn(atn1, one, dwords)
   fpmul(tail, two, atn1, dwords) ;2*atn(1)
   fpmul(T, minusone, *x, dwords) ;-x
   fpmul(B, *x, *x, dwords)       ;x*x
   If fpcmp(B, one, dwords) = 0
      ;for 1 and -1
      If fpcmp(*x, minusone, dwords) = 0
         Define.decfloat four
         si2fp(four, 4, dwords)
         fpmul(*result, four, atn1, dwords)
         ProcedureReturn
      Else
         si2fp(*result, 0, dwords)
         ProcedureReturn
      EndIf
   EndIf
   fpsub(B, one, B, dwords) ;1-x*x
   fpsqr(B, B, dwords)      ;sqr(1-x*x)
   fpdiv(term1, T, B, dwords)
   fpatn(atnterm1, term1, dwords)
   fpadd(*result, atnterm1, tail, dwords)
EndProcedure

[jack]

fixed a couple of bugs in the code in first post and posted logs and trigs in post above

[STARGÅTE]

Why do you define Pi as a String-Constant, when you already have a procedure pi_brent_salamin() to calculate Pi until the requiert precision.
You can also buffer the first evaluation of Pi in a cache.

BTW: Interesting approches for calculating fpsin and fpatn. Do you have a reference?
I use Taylor around 0 and clip the argument between 0 and Pi/2 in case of sin() and minimize the argument in case of arctan().

[jack]

hello STARGÅTE
I appreciate your feedback and suggestions, in answer to your question no I don't have a reference, the routines were written by a fellow at the FreeBasic forum with the alias dodicat, I tweaked them a bit for my purpose
you are right about calculating Pi instead of the string constant the only problem is if someone sets the precision much above 1000 then it may take a long time to compute, I think that I will make Pi a global but with a limit to the precision, the user can always change it if they want

[jack]

changed the first post to calculate Pi in global variable Pi_Dec, also fp2str now should roundup the value
simple test

Code: Select all

#number_of_digits = 10000
#maximum_Pi_Dec_precision = #number_of_digits ; Pi_Dec = Pi, default precision is 1024

XIncludeFile "decfloat.pbi"

OpenConsole()
Define.decfloat x, z
Define.s s
Define.d t

s="1.23456789"
str2fp(x, s)
t=ElapsedMilliseconds()
fpsin(z, x)
PrintN(fp2str(z))
fpasin(z, z)
PrintN(fp2str(z))
fpcos(z, x)
PrintN(fp2str(z))
fpacos(z, z)
PrintN(fp2str(z))
fptan(z, x)
PrintN(fp2str(z))
fpatn(z, z)
PrintN(fp2str(z))
t=ElapsedMilliseconds()-t
PrintN( "elapsed time "+StrD(t/1000)+" seconds")
Input()
CloseConsole()

compiled with the C backend it took 90 seconds

[STARGÅTE]

I have run you code also on my system, and it took 127 seconds.
I was curious to run the same code in my language Lizard. Interestingly, it took also quite long: 15 seconds.
But I observed, that only my ArcSin() need more than 6 seconds. So i used your definition with
; ARCSIN = ATN(x / SQR(-x * x + 1))
and the time was less than 0.2 seconds. And now, the whole code need just 1.7 seconds in Lizard!
So I am very grateful to you, that to showed me potential of optimization.

[jack]

hello STARGÅTE
I am happy that you found some of my code useful
I realize that this decimal floating point implementation is useful for relatively low precision due to it using arrays which are stored on the stack, you can quickly run out of stack memory and the program will simply crash
@STARGÅTE, another function of lizard that could be sped up is the Sqrt function, Calculate(Sqrt(2), 10000) took .47 seconds whereas the decfloat fpsqr took .15 seconds

[STARGÅTE]

@Jack: Thanks for the comment.
In fact, Sqrt(x) is converted to x^(1/2). The problem is now, Calculate() also "calculates" also the 1/2 to a decimal number and then Lizard do not realizes, that it was a simple square root and Lizard have to use Exp(0.5000... * Log(x)), which is of cause slower.
When you run: "two = Calculate(2, 10000); Sqrt(two)" it is quite fast, because then Lizard also used the Newton–Raphson method.
Again, thanks for the hint, I have to take a closer look.