128-bit integer division

General FreeBASIC programming questions.
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: 128-bit integer division

Post by dodicat »

I have my other divider (same principal) doing quotient and remainder.

Code: Select all

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

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

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

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

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

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

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

Function divide(Byval n As String,Byval d As String,places As Long=10,f As String="") As String
    Dim As Integer modstop,runcount
    If f="mod" Then 
        If Len(n)<Len(d) Then Return n
        If Len(n)=Len(d) Then
            If n<d Then Return n
        End If
        modstop=Len(n)
    End If
    Dim As String sign
    If Instr(n,"-") xor Instr(d,"-") Then sign="-"
    n=Ltrim(n,"0"):d=Ltrim(d,"0")
    remove(n,"-"):remove(d,"-")
    If Instr(n,".")=0 Then n+="."
    If Instr(d,".")=0 Then d+="."
    Dim As String s(0 To 10)'to hold digit multiples
    Dim As String ans
    Dim As Long k,flag
    Var dm=decmove(d)
    adjustnumerator(n,dm)
    n+=String(places-decmove(n),"0")
    remove(d,".")
    d=Ltrim(d,"0")
    n=Ltrim(n,"0")
    Dim As String top=Mid(n,1,1)
    remove(top,".")
    For m As Long=1 To Len(n)-1
        If n[m]=Asc(".") Then flag=1:Continue For
        For k =0 To 10 'range up the digits until the product > the numerator part
            If k=10 Then
                s(k)=d+"0"
            Else
                s(k)=smult(d,k)
            End If
            If isbigger(s(k),top) Then ans+=Str(k-1):runcount=runcount+1:Exit For
        Next k
        top=minus(top,s(k-1))
        top+=Chr(n[m])'bring down the next digit
        If flag=1 Then ans+=".":flag=0
        If f="mod" Andalso runcount>=modstop Then 
            If top="" Then Return "0"
            Return Mid(top,1,Len(top)-1)
        End If 
    Next m
    ans=Iif(Instr(ans,"."),Rtrim( Ltrim(ans,"0"),"."),ans)
    If Len(ans)=0 Then ans="0"
    Return sign+ ans
End Function

#define mod_(a,b) divide((a),(b),,"mod")
#define div_(a,b) iif(len((a))<len((b)),"0",divide((a),(b),0))
'===========================================================================

Type answer
    As String quotient,remainder
End Type

Function longdivide(Byval numerator As String,Byval denominator As String) As answer
    Dim As answer a
    a.quotient=div_(numerator,denominator)
    a.remainder=mod_(numerator,denominator)
    Return a
End Function
#define Intrange(f,l) int(Rnd*(((l)+1)-(f))+(f))

Randomize


Print "test a few small numbers"
Dim As answer x
For n As Long=1 To 100
    Var numerator=str(intrange(10,10000))
    Var denominator=str(intrange(10,10000))
    Print numerator;" / ";denominator
    x=longdivide(numerator,denominator)
    Print  x.quotient; "  remainder ";x.remainder
    Print "check ";Val(x.quotient)*Val(denominator)+Val(x.remainder)
    If Str(Val(x.quotient)*Val(denominator)+Val(x.remainder))<> numerator Then Print "ERROR":Sleep
    Print "-----------"
Next
Print "press a key"
Sleep

Var numerator=  "999999999999999999998888888888888888"
Var denominator="18446744073709551616"
Print numerator;" / ";denominator
x=longdivide(numerator,denominator)
Print  x.quotient; "  remainder ";x.remainder
print
numerator=string(2000,"9")
denominator=string(2000,"7")
x=longdivide(numerator,denominator)
Print  x.quotient; "  remainder ";x.remainder
print "Done"
Sleep

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

Re: 128-bit integer division

Post by srvaldez »

Hi dodicat :)
that's very good considering that you are using strings, unfortunately it's a bit slow for use in a library
have you considered converting your code to do binary arithmetic?
that way the numbers would be the same as native integer types but with more precision, say 128-bits
the code that I translated from C does binary arithmetic, but unfortunately in a very inefficient way
for example the division is done one bit at a time -- the way you used to do it with an old 8-bit CPU
dafhi
Posts: 1641
Joined: Jun 04, 2005 9:51

Re: 128-bit integer division

Post by dafhi »

that's cool, dodicat. i've seen vedic methods as well as a YT video about an innuit (Alaskan native) system if i remember correctly

i work at a snail's pace but one day i'd like to work that in code
srvaldez
Posts: 3373
Joined: Sep 25, 2005 21:54

Re: 128-bit integer division

Post by srvaldez »

here's the 64-bit compiled lib of fp256 with example by me
the fp256 Git repo is at https://github.com/piggypiggy/fp256

get the FB translated code of tiny-bignum-c from viewtopic.php?p=296718#p296718
and change const BN_ARRAY_SIZE = 16 / WORD_SIZE
to const BN_ARRAY_SIZE = 32 / WORD_SIZE to make a fair comparison with fp256

tiny-bignum-c example, #cmdline "-arch native -gen gcc -O 2"

Code: Select all

#If 1
    Dim As bn x, y, r, q
	Dim As Long i
	Dim As Double t

	valbn(x, "999999999999999999998888888888888888")
	valbn(y, "18446744073709551616")
	Print "x = 999999999999999999998888888888888888"
	Print "y = 18446744073709551616"
	Print "100000 iterations of bignum_divmod(x, y, q, r)"
	t=Timer
	For i=1 To 100000
		bignum_divmod(x, y, q, r)
	Next
	t=timer-t
	Print "quotient  = ";strbn(q)
	Print "remainder = ";strbn(r)
	Print using "elapsed time ####.#### seconds ";t
	Print "press return to exit ";
	Sleep
#Endif

#If 0
    Dim As bn x, y, r
	Dim As Long i
	Dim As Double t
    
	valbn(x, "54210108624275221")
	valbn(y, "18446744073709551616")
	Print "x = 54210108624275221"
	Print "y = 18446744073709551616"
	Print "100000 iterations of bignum_mul(x, y, r)"
	t=Timer
	For i=1 To 100000
		bignum_mul(x, y, r)
	Next
	t=timer-t

	Print "product = ";strbn(r)
	Print using "elapsed time ####.#### seconds ";t
	Print "press return to exit ";
	Sleep
#Endif
tiny-bignum-c time 0.1831 seconds
fp256 time 0.0095 seconds
fp256 is 19.27 times faster than tiny-bignum-c
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: 128-bit integer division

Post by dodicat »

gmp is quite fast

Code: Select all

#Include "gmp.bi"
#cmdline "-gen gcc -Wc -O2"


Dim Shared As ZString * 1000000 outtext
Type mpf_t As __mpf_struct
Function _mod(n1 As String,n2 As String) As String
    Dim As __mpz_struct answer, mn1, mn2
    __gmpz_init(@answer)
    __gmpz_init_set_str(@mn1, n1, 10)
    __gmpz_init_set_str(@mn2, n2, 10)
    __gmpz_mod(@answer, @mn1, @mn2)
    __gmp_sprintf( @outtext,"%Zi", @answer )
    __gmpz_clear(@answer) : __gmpz_clear(@mn1) : __gmpz_clear(@mn2)
    Return Trim(outtext)
End Function

Function _div(n1 As String,n2 As String) As String
    Dim As __mpz_struct answer, mn1, mn2
    __gmpz_init(@answer)
    __gmpz_init_set_str(@mn1, n1, 10)
    __gmpz_init_set_str(@mn2, n2, 10)
    mpz_div(@answer, @mn1, @mn2)
    __gmp_sprintf( @outtext,"%Zi", @answer )
    __gmpz_clear(@answer) : __gmpz_clear(@mn1) : __gmpz_clear(@mn2)
    Return Trim(outtext)
End Function

Var numerator=  "999999999999999999998888888888888888"
Var denominator="18446744073709551616"

dim as string s1=numerator,s2=denominator
dim as string m,d
dim as double t=timer
for n as long=1 to 10000
 m= _mod(s1,s2)
 d= _div(s1,s2)
next
print timer-t
print d; " remainder ";m
sleep 
result

Code: Select all

 0.05074829998193309
54210108624275221 remainder 12918483735999581752
 
srvaldez
Posts: 3373
Joined: Sep 25, 2005 21:54

Re: 128-bit integer division

Post by srvaldez »

Hi dodicat :)
have considered rendering your arithmetic algorithms into a binary 128-bit integer library?
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: 128-bit integer division

Post by dodicat »

Hi srvaldez.
I am not 100% sure what you mean, as you say, string indexing is a bit slow.
I have found some old C code which produces quotient and remainder division along with the other basic maths.
bignum.h

Code: Select all

/* +++Date last modified: 05-Jul-1997 */

/* BIGNUM.H
**
** Header file with definition of BigNum type for Big Number Arithmetic.
**
** Released to the public domain by Clifford Rhodes on June 15, 1995 with
** no guarantees of any sort as to accuracy or fitness of use.
*/

/* Each int in the Big Number array will hold a digit up to MODULUS - 1.
** Choosing 10000 makes testing easy because each digit contains 4 decimal
** places.
*/

#ifndef BIGNUM__H
#define BIGNUM__H

#include "minmax.h"

#define MODULUS  10000    /* Warning: Must be <= USHRT_MAX! */

typedef unsigned short UShort;
typedef unsigned long  ULong;

/* The Big Number is contained in a structure that has a length, nlen, and
** an array, n[], of unsigned shorts to hold the 'digits'. The most
** significant digit of the big number is at n[0]. The least significant
** digit is at n[nlen - 1];
*/

typedef struct {
      int nlen;      /* Number of int's in n */
      UShort *n;     /* Array of unsigned shorts to hold the number */
} BigNum;

/* Arithmetic function prototypes */
BigNum * BigNumAdd(const BigNum *a, const BigNum *b);
BigNum * BigNumSub(const BigNum *a, const BigNum *b);
BigNum * BigNumMul(const BigNum *a, const BigNum *b);
BigNum * BigNumDiv(const BigNum *a, const BigNum *b, BigNum **c);
BigNum * BigNumAlloc(int nlen);
void     BigNumFree(BigNum *b);

#endif /* BIGNUM__H */ 
and the running functions + an example (test.c say)

Code: Select all

/* +++Date last modified: 05-Jul-1997 */

/* BIGTEST.C
**
** Routines to test Big Number Arithmetic functions found in BIGNUM.C
**
*/

#include <stdlib.h>
#include <stdio.h>
#include <conio.h>
#include <time.h>
#if defined(_MSC_VER)
 #include <memory.h>
#elif defined(__TURBOC__)
 #include <mem.h>
#else
 #include <string.h>
#endif

#include "bignum.h"    /* Typedef for BigNum type */



/* +++Date last modified: 05-Jul-1997 */

/* BIGNUM1.C
**
** Routines to do Big Number Arithmetic. These are multi-precision, unsigned
** natural numbers (0, 1, 2, ...). For the storage method, see the BigNum
** typedef in file BIGNUM.H
**
** Released to the public domain by Clifford Rhodes, June 15, 1995, with
** no guarantees of any sort as to accuracy or fitness of use.
*/

#include <stdlib.h>
#include "bignum.h"    /* Typedef for BigNum type */
#include <string.h>


BigNum * BigNumAlloc(int nlen)
{
      /* Allocates memory for a BigNum object with nlen entries. Returns a
      ** pointer to the memory with the data array initialized to zero if
      ** successful. If not successful, NULL is returned.
      */

      BigNum * big;

      big = (BigNum *) malloc(sizeof(BigNum));
      if (big != NULL)
      {
            big->nlen = nlen;
            if (nlen < 1)
                  big->n = NULL;
            else if ((big->n =
                  (UShort *) calloc(nlen, sizeof(UShort))) == NULL)
            {
                  free(big);
                  return NULL;
            }
      }
      else  return NULL;

      return big;
}

void BigNumFree(BigNum * b)
{
      /* Frees memory used by BigNum object b. */

      if (b != NULL)
      {
            if (b->n != NULL)
                  free(b->n);
            free(b);
      }
}

BigNum * BigNumDiv(const BigNum * a, const BigNum * b, BigNum ** c)
{
      /* Big Number a is divided by Big Number b. If successful a pointer to
      ** the quotient is returned. The user must supply a pointer to a Big
      ** Number pointer, c, to receive the remainder.
      ** If unsuccessful, NULL is returned.
      ** Neither a nor b is changed by the call.
      */

      int      i, j, d, bpos;
      UShort   carry, quo;
      long     m1, m2, m3;
      BigNum * quotient, * atemp, * btemp;

      /* Find non-zero MSB of b, make sure b is not 0 */

      for (bpos = 0; bpos < b->nlen && b->n[bpos] == 0; bpos++)
            ;
      if (bpos == b->nlen)  /* Attempt to divide by zero! */
            return NULL;

      /* Create a copy of b, making sure btemp->n[0] != 0 */

      if ((btemp = BigNumAlloc(b->nlen - bpos)) == NULL)
            return NULL;
      memcpy(btemp->n, b->n + bpos, btemp->nlen * sizeof(UShort));

      /* Create a copy of a, making sure atemp->n[0] == 0 */

      carry = (a->n[0] == 0) ? 0 : 1;
      if ((atemp = BigNumAlloc(a->nlen + carry)) == NULL)
      {
            BigNumFree(btemp);
            return NULL;
      }
      memcpy(atemp->n + carry, a->n, a->nlen * sizeof(UShort));

      /* Allocate memory for quotient and remainder */

      if ((quotient = BigNumAlloc(max(1, atemp->nlen - btemp->nlen))) == NULL)
      {
            BigNumFree(atemp);
            BigNumFree(btemp);
            return NULL;
      }
      if ((*c = BigNumAlloc(btemp->nlen)) == NULL)
      {
            BigNumFree(atemp);
            BigNumFree(btemp);
            BigNumFree(quotient);
            return NULL;
      }
      d = MODULUS / (btemp->n[0] + 1);
      for (carry = 0, j = atemp->nlen - 1; j >= 0; j--)   
      {
            m1 = ((long) d * (long) *(atemp->n + j)) + (long) carry;
            *(atemp->n + j) = (UShort) (m1 % (long) MODULUS);
            carry = (UShort) (m1 / (long) MODULUS);
      }
      for (carry = 0, j = btemp->nlen - 1; j >= 0; j--)
      {
            m1 = ((long) d * (long) *(btemp->n + j)) + (long) carry;
            *(btemp->n + j) = (UShort) (m1 % (long) MODULUS);
            carry = (UShort) (m1 / (long) MODULUS);
      }
      for (j = 0; j < (atemp->nlen - btemp->nlen); j++)   
      {
            if (*btemp->n == *(atemp->n + j))
                  quo = MODULUS - 1;
            else
            {
                  m1 = (long) *(atemp->n + j) * (long) MODULUS;
                  m1 = (m1 + (long) *(atemp->n + j + 1)) / (long) *btemp->n;
                  quo = (UShort) m1;
            }
            m3 = (long) *(atemp->n + j) * (long) MODULUS +
                  (long) *(atemp->n + j + 1);
            do
            {
                  if (btemp->nlen > 1)
                        m1 = (long) quo * (long) *(btemp->n + 1);
                  else  m1 = 0l;
                  m2 = (long) quo * (long) *btemp->n;
                  m2 = (m3 - m2) * (long) MODULUS +
                        (long) *(atemp->n + j + 2);
                  if (m1 > m2)
                        quo--;
            } while (m1 > m2);

            bpos = btemp->nlen - 1;
            i = j + btemp->nlen;
            m2 = 0l;
            while (i >= j)
            {
                  if (bpos >= 0)
                        m1 = (long) quo * (long) *(btemp->n + bpos);
                  else  m1 = 0l;
                  m3 = (long) *(atemp->n + i) - m1 + m2;
                  if (m3 < 0l)
                  {
                        m2 = m3 / (long) MODULUS;
                        m3 %= (long) MODULUS;
                        if (m3 < 0l)
                        {
                              m3 += (long) MODULUS;
                              m2--;
                        }
                  }
                  else  m2 = 0l;
                  *(atemp->n + i) = (UShort) (m3);
                  bpos--;
                  i--;
            }
            if (m2 == 0l)
                  *(quotient->n + j) = quo;
            else
            {
                  *(quotient->n + j) = quo - 1;
                  bpos = btemp->nlen - 1;
                  i = j + btemp->nlen;
                  carry = 0;
                  while (i >= j)
                  {
                        long tmp;

                        if (bpos >= 0)
                              carry += *(btemp->n + bpos);
                        tmp = carry + (long) *(atemp->n + i);
                        if (tmp >= (long) MODULUS)
                        {
                              tmp -= (long) MODULUS;
                              carry = 1;
                        }
                        else  carry = 0;
                        *(atemp->n + i) = (UShort) tmp;
                        bpos--;
                        i--;
                  }
            }
      }
      j = atemp->nlen - btemp->nlen;    
      bpos = 0;
      carry = 0;
      while (j < atemp->nlen)            
      {
            m3 = (long) carry * (long) MODULUS + (long) *(atemp->n + j);
            *((*c)->n + bpos) = (UShort) (m3 / d);
            carry = (UShort) (m3 % d);
            j++;
            bpos++;
      }
      BigNumFree(atemp);   /* Free temporary memory */
      BigNumFree(btemp);

      return quotient;
}


BigNum * BigNumAdd(const BigNum * a, const BigNum * b)
{
      /* Big Numbers a and b are added. If successful a pointer to the sum is
      ** returned. If unsuccessful, a NULL is returned.
      ** Neither a nor b is changed by the call.
      */

      UShort  carry = 0;
      int     size, la, lb;
      long    tsum;
      BigNum  * sum;

      /* Allocate memory for Big Number sum to be returned... */

      size = max(a->nlen, b->nlen) + 1;
      if ((sum = BigNumAlloc(size)) == NULL)
            return NULL;

      /* Get indexes to last digits in each number (Least significant) */

      la = a->nlen - 1;
      lb = b->nlen - 1;
      size--;

      while (la >= 0 && lb >= 0)    /* While both a and b have data */
      {
            tsum = carry + (long) *(a->n + la) + (long) *(b->n + lb);
            *(sum->n + size) = (UShort) (tsum % (long) MODULUS);
            carry = (tsum / (long) MODULUS) ? 1 : 0;
            la--;
            lb--;
            size--;
      }
      if (lb < 0)                   /* Must see if a still has data */
      {
            while (carry && la >= 0)
            {
                  tsum = carry + (long) *(a->n + la);
                  *(sum->n + size) = (UShort) (tsum % (long) MODULUS);
                  carry = (tsum / (long) MODULUS) ? 1 : 0;
                  la--;
                  size--;
            }
            while (la >= 0)
            {
                  *(sum->n + size) = *(a->n + la);
                  la--;
                  size--;
            }
      }
      else                          /* See if b still has data */
      {
            while (carry && lb >= 0)
            {
                  tsum = carry + (long) *(b->n + lb);
                  *(sum->n + size) = (UShort) (tsum % (long) MODULUS);
                  carry = (tsum / (long) MODULUS) ? 1 : 0;
                  lb--;
                  size--;
            }
            while (lb >= 0)
            {
                  *(sum->n + size) = *(b->n + lb);
                  lb--;
                  size--;
            }
      }
      *(sum->n + size) = carry;

      return sum;
}

BigNum * BigNumSub(const BigNum * a, const BigNum * b)
{
      /* Big Numbers a and b are subtracted. a must be >= to b. A pointer to
      ** the difference (a - b) is returned if successful. If unsuccessful,
      ** a NULL is returned.
      ** Neither a nor b is changed by the call.
      */

      int      size, la, lb, borrow = 0;
      long     tdiff;
      BigNum * diff;

      /* Allocate memory for Big Number difference to be returned... */

      if ((diff = BigNumAlloc(a->nlen)) == NULL)
            return NULL;

      la = a->nlen - 1;
      size = la;
      lb = b->nlen - 1;

      while (lb >= 0)
      {
            tdiff = (long) *(a->n + la) - (long) *(b->n + lb) - borrow;
            if (tdiff < 0)
            {
                  tdiff += (long) (MODULUS - 1);
                  borrow = 1;
            }
            else  borrow = 0;
            *(diff->n + size) = (UShort) tdiff + borrow;
            la--;
            lb--;
            size--;
      }
      while (la >= 0)          /* Must get rest of a->n */
      {
            tdiff = (long) *(a->n + la) - borrow;
            if (tdiff < 0)
            {
                  tdiff += (long) (MODULUS - 1);
                  borrow = 1;
            }
            else  borrow = 0;
            *(diff->n + size) = (UShort) tdiff + borrow;
            la--;
            size--;
      }
      if (borrow)   /* We've got an underflow here! */
      {
            BigNumFree(diff);
            return NULL;
      }
      return diff;
}

BigNum * BigNumMul(const BigNum * a, const BigNum * b)
{
      /* Big Numbers a and b are multiplied. A pointer to the product
      ** is returned if successful. If unsuccessful, a NULL is returned.
      ** Neither a nor b is changed by the call.
      */

      int      size, la, lb, apos, bpos, ppos;
      UShort   carry;
      BigNum * product;

      /* Allocate memory for Big Number product to be returned... */

      size = a->nlen + b->nlen;
      if ((product = BigNumAlloc(size)) == NULL)
            return NULL;

      la = a->nlen - 1;
      lb = b->nlen - 1;
      size--;

      bpos = lb;

      while (bpos >= 0)             /* We'll multiply by each digit in b */
      {
            ppos = size + (bpos - lb);    /* Position in product */

            if (*(b->n + bpos) == 0) /* If zero multiplier, skip this pass */
                  ppos = ppos - la - 1;
            else                    /* Multiply a by next b digit */
            {
                  apos = la;
                  carry = 0;
                  while (apos >= 0) /* Go a digit at a time through a */
                  {
                        ULong temp;

                        temp = (ULong) *(a->n + apos) *
                              (ULong) *(b->n + bpos) + carry;

                        /*
                        ** We must add any previous product term in
                        ** this 'column' too
                        */

                        temp += (ULong) *(product->n + ppos);

                        /* Now update product term, remembering any carry */

                        *(product->n + ppos) =
                              (UShort) (temp % (ULong) MODULUS);
                        carry = (UShort) (temp / (ULong) MODULUS);

                        apos--;
                        ppos--;
                  }
                  *(product->n + ppos) = carry;
            }
            bpos--;
      }
      return product;
}


//=======================================================//


ULong Big2Long(BigNum * b)
{
      /* Converts the big number in b to an unsigned long. This is a support
      ** routine for test purposes only. The length of b should be 3 or less
      ** to avoid overflowing the long.
      */

      int i;
      ULong res = 0L;

      for (i = 0; i < b->nlen; i++)
      {
            res *= (ULong) MODULUS;
            res += (ULong) b->n[i];
      }
      return res;
}

void BigNumRand(BigNum **b)
{
      /* This routine fills the big number pointed to by *b to a random long.
      ** It is a support routine for test purposes only.
      */

      int  i = (*b)->nlen - 1;

      memset((*b)->n, 0, (*b)->nlen * sizeof(UShort));

      if (i < 0)
            return;
      else if (i == 0)
            (*b)->n[0] = rand() % MODULUS;
      else if (i == 1)
      {
            (*b)->n[0] = rand() % MODULUS;
            (*b)->n[1] = rand() % MODULUS;
      }
      else
      {
            (*b)->n[i--] = rand() % MODULUS;
            (*b)->n[i--] = rand() % MODULUS;
            (*b)->n[i] = rand() % 21;
      }
}

int main(void)
{
      /* Test the big number arithmetic routines using random longs. Runs
      ** until an error is encountered or a key pressed.
      ** Gives examples of how to call the routines.
      */

      BigNum *a, *b, *c, *d = NULL, *p, *s;
      ULong al, bl, cl, dl, sl, cnt = 0L;

      /* 'Randomly' seed the rand() generator */

      srand((unsigned int) time(NULL));

      /* Create two big numbers */

      a = BigNumAlloc(3);
      b = BigNumAlloc(3);

      do    /* Loop until a key is pressed... */
      {
            if ((++cnt % 1000L) == 0)     /* Show every 1000 iterations */
                  printf("%ld\n", cnt);
            if (kbhit())            /* If a key pressed--we're outta here! */
            {
                  getch();
                  break;
            }

            /* Make the two big numbers random... */

            BigNumRand(&a);
            BigNumRand(&b);

            /* Get their 'long' values... */

            al = Big2Long(a);
            bl = Big2Long(b);

            c = BigNumDiv(a, b, &d);        /* Divide a by b */
            if (c == NULL)
            {
                  printf("Error: Divide did not return quotient\n");
                  break;
            }
            if (d == NULL)
            {
                  printf("Error: Divide did not return remainder\n");
                  break;
            }
            cl = Big2Long(c);        /* Get quotient as a long */
            dl = Big2Long(d);        /* Get remainder as a long */

            if ((al / bl) != cl)
            {
                  printf("Error: Quotient %lu / %lu != %lu\n", al, bl, cl);
                  printf("a = %d%04d%04d   b = %d%04d\n",
                        a->n[0], a->n[1], a->n[2], b->n[0], b->n[1]);
                  break;
            }
            if ((al % bl) != dl)
            {
                  printf("Error: Remainder %lu %c %lu != %lu\n",
                        al, '%', bl, dl);
                  printf("a = %d%04d%04d   b = %d%04d\n",
                        a->n[0], a->n[1], a->n[2], b->n[0], b->n[1]);
                  break;
            }

            /* Reconstruct BigNum a by multiplying the quotient by b and
            ** adding the remainder.
            */

            p = BigNumMul(c, b);
            s = BigNumAdd(p, d);
            sl = Big2Long(s);
            if (sl != al)
            {
                  printf("Error: Couldn't reconstruct original "
                        "divisor: %lu != %lu\n", al, sl);
                  break;
            }
            BigNumFree(s);    /* To use s again, must free it first! */
            s = BigNumSub(a, b);
            if (al >= bl)     /* For BigNumSub(), a must be > than b */
            {
                  sl = Big2Long(s);
                  if (sl != (al - bl))
                  {
                        printf("Error: Subtraction error %lu - %lu != %lu\n",
                              al, bl, sl);
                        break;
                  }
            }
            else if (s != NULL)
            {
                  printf("Invalid subtraction %lu - %lu\n", al, bl);
                  break;
            }

            /* Free all the memory allocated by the arithmetic calls */

            BigNumFree(c);
            BigNumFree(d);
            BigNumFree(p);
            BigNumFree(s);

      } while (1);

      BigNumFree(a);
      BigNumFree(b);

      return 0;
} 
Compiles perfectly with 64 bit gcc, probably 32 bit gcc also.
It is very old code, but you know what they say about old fiddles.
I have made a static fb lib, and am fiddling around to get it to work on fb.
srvaldez
Posts: 3373
Joined: Sep 25, 2005 21:54

Re: 128-bit integer division

Post by srvaldez »

thanks dodicat :)
basically what I mean by binary integer is an integer of base 2^n, preferably n should be 32 so it's limbs will fit in a ulong, tiny-bignum-c fits the bill except that the algorithm's it uses are slow
thank you for the tip about bignum by Clifford Rhodes, Google led me to a site that has his files among others, http://www.shedai.net/c/new/ follow the link under "Assorted math functions"
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: 128-bit integer division

Post by dodicat »

Hi srvaldez.
Did you post an example using Clifford Rhodes's code.
Did you make a lib or translate from c?
I did a few (multiply) tests from a lib and found that many were OK, but some failed by one digit.
Obviously I will have to investigate.
srvaldez
Posts: 3373
Joined: Sep 25, 2005 21:54

Re: 128-bit integer division

Post by srvaldez »

yes I had a post about the FB translation of Clifford Rhodes's code but I found that it was buggy so I removed it, however, I re-tranlated and I think that the translation is ok

Code: Select all

#Pragma Once

#Include Once "crt/long.bi"
#Include Once "crt/stdlib.bi"
#Include Once "crt/string.bi"

Extern "C"

#Define BIGNUM__H
#Define max(X, Y) Iif((X) > (Y), (X), (Y))
Const MODULUS = 10000

Type BigNum
	nlen As Short
	n As Ushort Ptr
End Type

Declare Function BigNumAdd(Byval a As Const BigNum Ptr, Byval b As Const BigNum Ptr) As BigNum Ptr
Declare Function BigNumSub(Byval a As Const BigNum Ptr, Byval b As Const BigNum Ptr) As BigNum Ptr
Declare Function BigNumMul(Byval a As Const BigNum Ptr, Byval b As Const BigNum Ptr) As BigNum Ptr
Declare Function BigNumDiv(Byval a As Const BigNum Ptr, Byval b As Const BigNum Ptr, Byval c As BigNum Ptr Ptr) As BigNum Ptr
Declare Function BigNumAlloc(Byval nlen As Short) As BigNum Ptr
Declare Sub BigNumFree(Byval b As BigNum Ptr)

Sub BigCopy(Byref res As BigNum Ptr, Byval num As Const BigNum Ptr)
	Dim As Long i, p
	p=0
	While num->n[p]=0 Andalso p<num->nlen
		p+=1
	Wend
	BigNumFree(res)
	res = BigNumAlloc(num->nlen-p)
	For i=0 To (num->nlen-1)-p
		res->n[i]=num->n[i+p]
	Next
End Sub

Private Function BigNumAdd(Byval a As Const BigNum Ptr, Byval b As Const BigNum Ptr) As BigNum Ptr
	Dim carry As Ushort = 0
	Dim size As Short
	Dim la As Short
	Dim lb As Short
	Dim tsum As clong
	Dim sum As BigNum Ptr
	size = Iif(a->nlen > b->nlen, a->nlen, b->nlen) + 1
	sum = BigNumAlloc(size)
	If (sum = NULL) Then Return NULL
	la = a->nlen - 1
	lb = b->nlen - 1
	size-=1
	While (la >= 0) Andalso (lb >= 0)
		tsum = (carry + Cast(clong, *(a->n + la))) + Cast(clong, *(b->n + lb))
		(*(sum->n + size)) = Cast(Ushort, tsum Mod Cast(clong, MODULUS))
		carry = Iif(tsum \ Cast(clong, MODULUS), 1, 0)
		la-=1
		lb-=1
		size-=1
	Wend
	If lb < 0 Then
		While carry Andalso (la >= 0)
			tsum = carry + Cast(clong, *(a->n + la))
			(*(sum->n + size)) = Cast(Ushort, tsum Mod Cast(clong, MODULUS))
			carry = Iif(tsum \ Cast(clong, MODULUS), 1, 0)
			la-=1
			size-=1
		Wend
		While la >= 0
			(*(sum->n + size)) = *(a->n + la)
			la-=1
			size-=1
		Wend
	Else
		While carry Andalso (lb >= 0)
			tsum = carry + Cast(clong, *(b->n + lb))
			(*(sum->n + size)) = Cast(Ushort, tsum Mod Cast(clong, MODULUS))
			carry = Iif(tsum \ Cast(clong, MODULUS), 1, 0)
			lb-=1
			size-=1
		Wend
		While lb >= 0
			(*(sum->n + size)) = *(b->n + lb)
			lb-=1
			size-=1
		Wend
	End If
	(*(sum->n + size)) = carry
	Return sum
End Function

Private Function BigNumSub(Byval a As Const BigNum Ptr, Byval b As Const BigNum Ptr) As BigNum Ptr
	Dim size As Short
	Dim la As Short
	Dim lb As Short
	Dim borrow As Short = 0
	Dim tdiff As clong
	Dim diff As BigNum Ptr
	diff = BigNumAlloc(a->nlen)
	If (diff = NULL) Then Return NULL
	la = a->nlen - 1
	size = la
	lb = b->nlen - 1
	While lb >= 0
		tdiff = (Cast(clong, *(a->n + la)) - Cast(clong, *(b->n + lb))) - borrow
		If tdiff < 0 Then
			tdiff += Cast(clong, MODULUS - 1)
			borrow = 1
		Else
			borrow = 0
		End If
		(*(diff->n + size)) = Cast(Ushort, tdiff) + borrow
		la-=1
		lb-=1
		size-=1
	Wend
	While la >= 0
		tdiff = Cast(clong, *(a->n + la)) - borrow
		If tdiff < 0 Then
			tdiff += Cast(clong, MODULUS - 1)
			borrow = 1
		Else
			borrow = 0
		End If
		(*(diff->n + size)) = Cast(Ushort, tdiff) + borrow
		la-=1
		size-=1
	Wend
	If borrow Then
		BigNumFree(diff)
		Return NULL
	End If
	Return diff
End Function

Private Function BigNumMul(Byval a As Const BigNum Ptr, Byval b As Const BigNum Ptr) As BigNum Ptr
	Dim size As Short
	Dim la As Short
	Dim lb As Short
	Dim apos As Short
	Dim bpos As Short
	Dim ppos As Short
	Dim carry As Ushort
	Dim As BigNum Ptr product, an, bn

	BigCopy an, a
	BigCopy bn, b
	size = an->nlen-la + bn->nlen-lb
	product = BigNumAlloc(size)
	If (product = NULL) Then Return NULL
	la = an->nlen - 1
	lb = bn->nlen - 1
	size-=1
	bpos = lb
	While bpos >= 0
		ppos = size + (bpos - lb)
		If (*(bn->n + bpos)) = 0 Then
			ppos = (ppos - la) - 1
		Else
			apos = la
			carry = 0
			While apos >= 0
				Dim temp As Ulong
				temp = (Cast(Ulong, *(an->n + apos)) * Cast(Ulong, *(bn->n + bpos))) + carry
				temp += Cast(Ulong, *(product->n + ppos))
				(*(product->n + ppos)) = Cast(Ushort, temp Mod Cast(Ulong, MODULUS))
				carry = Cast(Ushort, temp \ Cast(Ulong, MODULUS))
				apos-=1
				ppos-=1
			Wend
			(*(product->n + ppos)) = carry
		End If
		bpos-=1
	Wend
	Return product
End Function

Private Function BigNumAlloc(Byval nlen As Short) As BigNum Ptr
	Dim big As BigNum Ptr
	big = Cptr(BigNum Ptr, malloc(Sizeof(BigNum)))
	If big <> NULL Then
		big->nlen = nlen
		If (nlen < 1) Then
			big->n = NULL
		Else
			big->n = Cptr(Ushort Ptr, calloc(nlen, Sizeof(Ushort)))
		If (big->n = NULL) Then
			free(big)
			Return NULL
		End If
	End If
	Else
		Return NULL
	End If
	Return big
End Function

Private Sub BigNumFree(Byval b As BigNum Ptr)
	If b <> NULL Then
		If b->n <> NULL Then
			free(b->n)
		End If
		free(b)
	End If
End Sub

Private Function BigNumDiv(Byval a As Const BigNum Ptr, Byval b As Const BigNum Ptr, Byval c As BigNum Ptr Ptr) As BigNum Ptr
	Dim i As Short
	Dim j As Short
	Dim d As Short
	Dim tmp As Short
	Dim bpos As Short
	Dim carry As Ushort
	Dim quo As Ushort
	Dim m1 As clong
	Dim m2 As clong
	Dim m3 As clong
	Dim quotient As BigNum Ptr
	Dim atemp As BigNum Ptr
	Dim btemp As BigNum Ptr
	bpos = 0
	While bpos < b->nlen Andalso b->n[bpos] = 0
		bpos += 1
	Wend
	If bpos = b->nlen Then
		Return NULL
	End If
	btemp = BigNumAlloc(b->nlen - bpos)
	If (btemp = NULL) Then Return NULL
	memcpy(btemp->n, b->n + bpos, btemp->nlen * Sizeof(Ushort))
	carry = Iif(a->n[0] = 0, 0, 1)
	atemp = BigNumAlloc(a->nlen + carry)
	If (atemp = NULL) Then
		BigNumFree(btemp)
		Return NULL
	End If
	memcpy(atemp->n + carry, a->n, a->nlen * Sizeof(Ushort))
	tmp=Iif((1 > (atemp->nlen - btemp->nlen)), 1, (atemp->nlen - btemp->nlen))
	quotient = BigNumAlloc(tmp)
	If (quotient = NULL) Then
		BigNumFree(atemp)
		BigNumFree(btemp)
		Return NULL
	End If
	*c = BigNumAlloc(btemp->nlen)
	If (*c = NULL) Then
		BigNumFree(atemp)
		BigNumFree(btemp)
		BigNumFree(quotient)
		Return NULL
	End If
	d = MODULUS \ (btemp->n[0] + 1)
	carry = 0
	j = atemp->nlen - 1
	While (j >= 0)
		m1 = (Clng(d) * Clng(*(atemp->n + j))) + Clng(carry)
		*(atemp->n + j) = Cushort(m1 Mod Clng(MODULUS))
		carry = Cushort(m1 \ Clng(MODULUS))
		j-=1
	Wend
	carry = 0
	j = btemp->nlen - 1
	While (j >= 0)
		m1 = (Clng(d) * Clng(*(btemp->n + j))) + Clng(carry)
		*(btemp->n + j) = Cushort(m1 Mod Clng(MODULUS))
		carry = Cushort(m1 \ Clng(MODULUS))
		j-=1
	Wend
	j = 0
	While (j < (atemp->nlen - btemp->nlen))
		If (*btemp->n = *(atemp->n + j)) Then
			quo = MODULUS - 1
		Else
			m1 = Clng(*(atemp->n + j)) * Clng(MODULUS)
			m1 = (m1 + Clng(*(atemp->n + j + 1))) \ Clng(*btemp->n)
			quo = Cushort(m1) 
		End If
		m3 = Clng(*(atemp->n + j)) * Clng(MODULUS) + Clng(*(atemp->n + j + 1))
		Do
			If (btemp->nlen > 1) Then
				m1 = Clng(quo) * Clng(*(btemp->n + 1))
			Else
				m1 = 0l
			End If
			m2 = Clng(quo) * Clng(*btemp->n)
			m2 = (m3 - m2) * Clng(MODULUS) + Clng(*(atemp->n + j + 2))
			If (m1 > m2) Then quo-=1
		Loop While (m1 > m2)
		bpos = btemp->nlen - 1
		i = j + btemp->nlen
		m2 = 0l
		While (i >= j)
			If (bpos >= 0) Then
				m1 = Clng(quo) * Clng(*(btemp->n + bpos))
			Else
				m1 = 0l
			End If
			m3 = Clng(*(atemp->n + i)) - m1 + m2
			If (m3 < 0l) Then
				m2 = m3 \ Clng(MODULUS)
				m3 Mod= Clng(MODULUS)
				If (m3 < 0l) Then
					m3 += Clng(MODULUS)
					m2-=1
				End If
			Else
				m2 = 0l
			End If
			*(atemp->n + i) = Cushort(m3)
			bpos-=1
			i-=1
		Wend
		If (m2 = 0l) Then
			*(quotient->n + j) = quo
		Else
			*(quotient->n + j) = quo - 1
			bpos = btemp->nlen - 1
			i = j + btemp->nlen
			carry = 0
			While (i >= j)
				Dim As Long tmp
				If (bpos >= 0) Then carry += *(btemp->n + bpos)
				tmp = carry + Clng(*(atemp->n + i))
				If (tmp >= Clng(MODULUS)) Then
					tmp -= Clng(MODULUS)
					carry = 1
				Else
					carry = 0
				End If
				*(atemp->n + i) = Cushort(tmp)
				bpos-=1
				i-=1 
			Wend
		End If
		j+=1
	Wend
	j = atemp->nlen - btemp->nlen
	bpos = 0
	carry = 0
	While j < atemp->nlen
		m3 = (Cast(clong, carry) * Cast(clong, MODULUS)) + Cast(clong, *(atemp->n + j))
		(*((*c)->n + bpos)) = Cast(Ushort, m3 \ d)
		carry = Cast(Ushort, m3 Mod d)
		j+=1
		bpos+=1
	Wend
	BigNumFree(atemp)
	BigNumFree(btemp)
	Return quotient
End Function


Private Function Big2Long(Byval b As BigNum Ptr) As Ulong
'       Converts the big number in b to an unsigned long. This is a support
'      ** routine for test purposes only. The length of b should be 3 or less
'      ** to avoid overflowing the long.
'      
	Dim i As Short
	Dim res As Ulong = Cast(clong, 0)
	For i = 0 To b->nlen - 1
		res *= Cast(clong, MODULUS)
		res += Culng(b->n[i])
	Next
	Return res
End Function

Private Sub BigNumRand(Byval b As BigNum Ptr Ptr)
	Dim i As Short = (*b)->nlen - 1
	memset((*b)->n, 0, (*b)->nlen * Sizeof(Ushort))
	If i < 0 Then
		Return
	Elseif i = 0 Then
		(*b)->n[0] = rand() Mod MODULUS
	Elseif i = 1 Then
		(*b)->n[0] = rand() Mod MODULUS
		(*b)->n[1] = rand() Mod MODULUS
	Else
		(*b)->n[i] = Rnd Mod MODULUS : i-=1
		(*b)->n[i] = Rnd Mod MODULUS : i-=1
		(*b)->n[i] = Rnd Mod 21
	End If
End Sub

End Extern

Function BigStr(Byval num As BigNum Ptr) As String
	Dim As Long i, ln
	Dim As String s, d
	If num=0 Then
		Exit Function
	End If	
	s=Str(num->n[0])
	For i=1 To num->nlen-1
		d=Ltrim(Str(num->n[i]))
		ln=Len(d)
		If ln<4 Then
			d=String(4-ln, "0")+d
		End If
		s+=d
	Next
	s=Ltrim(s, "0")
	If s="" Then s="0"
	Return s
End Function

Sub BigVal(Byref num As BigNum Ptr, Byref strn As String)
	Dim As String digits, s=strn
	Dim As Long i, ln, nlen
	s=Trim(s)
	s=Ltrim(s, "0")
	ln=Len(s)
	nlen=ln\4
	If (4*nlen)<ln Then nlen+=1
	If nlen < num->nlen Then nlen = num->nlen
	If ln=0 Then
		If num=0 Then num = BigNumAlloc(nlen)
		Exit Sub
	End If
	If num=0 Then
		num = BigNumAlloc(nlen)
	End If
	If num=0 Then
		Exit Sub
	End If
	If num->nlen<nlen Then
		BigNumFree(num)
		num = BigNumAlloc(nlen)
	End If
	While ln>0
		digits=Right(s, 4)
		s=Left(s, ln-4)
		ln-=4
		nlen-=1
		num->n[nlen]=Cushort(digits)
	Wend
	If nlen>0 Then
		For i=0 To nlen-1
			num->n[i]=0
		Next
	End If
End Sub

Sub Long2BigNum(Byref res As BigNum Ptr, Byval num As Ulong)
	Dim As Long i
	If res->nlen<3 Then
		BigNumFree(res)
		res = BigNumAlloc(3)
	End If
	res->n[res->nlen-1]=num Mod MODULUS
	num\=MODULUS
	res->n[res->nlen-2]=num Mod MODULUS
	num\=MODULUS
	res->n[res->nlen-3]=num
End Sub

'main_test

Dim As BigNum Ptr a, b, c, d
Dim As Long i, l
Dim As Double t
Dim As String s, dg

a = BigNumAlloc(3)
b = BigNumAlloc(3)
c = BigNumAlloc(3)
d = BigNumAlloc(3)

BigVal a, "999999999999999999998888888888888888"
BigVal b, "18446744073709551616"

t=Timer
For i=1 To 100000
	c=BigNumDiv(a, b, @d)
Next i
t=timer-t

Print BigStr(c)
Print "54210108624275221"
Print BigStr(d)
Print "12918483735999581752"
Print Using "elapsed time ####.#### seconds ";t
Print "==================================================================="
BigVal(a, "3141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068")

BigCopy c, a

t=Timer

For i=1 To 8
	c=BigNumMul(c, c)
Next

For i=1 To 255
	c=BigNumDiv(c, a, @d)
Next

t=timer-t

Print BigStr(c)
b=BigNumSub(c, a)
Print Using "elapsed time ####.#### seconds ";t


BigNumFree(d)
BigNumFree(c)
BigNumFree(b)
BigNumFree(a)


a = BigNumAlloc(3)
b = BigNumAlloc(3)
c = BigNumAlloc(3)
d = BigNumAlloc(3)

t=Timer
Long2BigNum c, 1
For i=2 To 1000
	Long2BigNum a, i
	c=BigNumMul(c, a)
Next
t=timer-t
Print BigStr(c)

BigNumFree(d)
BigNumFree(c)
BigNumFree(b)
BigNumFree(a)

Print Using "elapsed time ####.#### seconds ";t
Print "press return to exit ";
Sleep
Last edited by srvaldez on Feb 12, 2023 10:43, edited 1 time in total.
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: 128-bit integer division

Post by dodicat »

Thanks srvaldez, your declared functions work perfectly with a static lib from the C code.
I get about .4 seconds for the timed bit.
Post Reply