Big Integer Library

[p_m21987]

Richard Russell has made this post on the raspberry pi forum. I thought it would be appropriate to cross post it here.

https://www.raspberrypi.org/forums/view ... 8#p1466828

Re: Introduction to BBC BASIC
Sun May 12, 2019 11:08 am

I mentioned in the other thread that I was thinking of writing a native code BigInt library for BBC BASIC; it would be an interesting exercise and hopefully provide a route to an efficient solution of the Fibonacci challenge. I've started to look at this, but I'm facing a dilemma: should I use a 10^n radix for the limbs (as the classic BASIC Fibo program does) or a more conventional 2^n radix? The advantage of the former is that the conversion back from the BigInt to decimal becomes trivial (and fast) which is of particular value when outputting a million-digit result! The disadvantage - and it's a big one - is that the BigInt computations themselves become significantly less efficient.

For a general-purpose BigInt library I'm in little doubt that the 2^n limb radix is better. In a typical application the formatting of the output as a decimal number is likely to be a small proportion of the total work (it might not even always be needed, if the output value isn't intended for 'human consumption') and efficiency of the BigInt arithmetic will be the most important factor (GMP uses a 2^n radix incidentally). But I suspect that in the specific case of the million-digit Fibonacci challenge it would give poor results because of the time taken to convert the result to decimal (which I assume is supposed to be part of the challenge).

What to do?! Does anybody have some BASIC code for multiple-precision binary-to-decimal conversion that I could benchmark?

[DDRM]

To follow up on this, Richard has made substantial progress in developing a limited BigInt library, and would value input from others to extend it. I copy, with permission, his comments to me:

Things have moved on since then: I've written a barebones library with sufficient functionality to solve the Fibonacci challenge, meaning only unsigned integers and just add, subtract and multiply operations (increment and decrement added today). I'm not particularly motivated to extend it to become a general purpose bigint library, by the addition of signed numbers and division etc., hence my desire to find out whether somebody else might want to take it on.

Anyone interested? As far as I'm aware the current code isn't available, but if you contact Richard privately he will presumably be happy to share it as a foundation for collaboration - hopefully he will release it anyway - it is clearly already useful in some contexts!

Best wishes,

D

[hellomike]

I sincerely hope that Richard will make it available for BB4W 6.x as well.

[DDRM]

As far as I understand, it is written purely in BBC BASIC, so it should run fine with BB4W. That also makes it more accessible to anyone wanting to help out, since I suspect that there are a limited number of us (not including me!) who might be able to write the routines in assembler...

Best wishes,

D

[Richard Russell]

In principle the Big Integer library (bigint.bbc) could be substantially simplified, and its performance improved, with the availability of array slicing. The existing library has to do nasty things like use SYS "memove" to shift numbers (which are represented internally as arrays), but with array slicing this can be done directly in BASIC - although some care is necessary when the destination and source slices overlap.

For example to shift an array of 64-bit integers left one element one can do:

Code: Select all

      N% = DIM(array%%(),1)
      array%%(0 TO N%-1) = array%%(1 TO N%)

I wonder if anybody would be interested in tackling this. It would probably mean re-writing it from scratch to get the greatest benefit

[Richard Russell]

I wonder if anybody would be interested in tackling this.

Seemingly not, so here's a (hopefully) simpler challenge that somebody might pick up: to write an arbitrary-precision integer division routine to add to the bigint library. The library has integer addition, subtraction and multiplication, but not division.

The most obvious and straightforward approach is to implement the 'grade school' long division routine, formalised as Knuth's Algorithm D. To remind you, the bigint library represents a big integer as an array of 'limbs' each with a radix of 10^18.

It would be really helpful if somebody would offer to attempt this. I did try asking DeepSeek, but I was insufficiently precise with the specification and it used 128-bit temporaries in the solution, which BBC BASIC doesn't have. I should probably have asked it to fix that, but I didn't.

[Richard Russell]

Seemingly not, so here's a (hopefully) simpler challenge that somebody might pick up: to write an arbitrary-precision integer division routine to add to the bigint library. The library has integer addition, subtraction and multiplication, but not division.

No takers even for that simpler task, so I wrote this naïve long-division code (not based on Knuth's Algorithm D) which works but is very slow:

Code: Select all

      REM Unsigned integer division, q%% = a%% DIV b%%:
      DEF PROCbigudiv(q%%,a%%,b%%)
      LOCAL U%,V%,q%%(),a%%(),b%%()
      PTR(q%%()) = q%% : PTR(a%%()) = a%% : PTR(b%%()) = b%%
      REM Remove leading zeroes:
      U% = ABS(a%%(0)) : V% = ABS(b%%(0))
      WHILE a%%(U%) = 0 AND U% > 0 U% -= 1 : ENDWHILE
      WHILE b%%(V%) = 0 AND V% > 0 V% -= 1 : ENDWHILE
      DEF PROC_bigudiv(q%%(),a%%(),b%%(),U%,V%)
      LOCAL J%,d%%,h%%,l%%,t%%(),u%%(),v%%()

      REM Check for division by zero
      IF V% = 0 THEN ERROR 100,"Division by zero"

      REM Handle trivial cases
      IF U% < V% THEN
        q%%(0) = 0 : q%%(1) = 0
        ENDPROC
      ENDIF

      REM Dimension arrays
      DIM t%%(U%+1),u%%(U%+2),v%%(V%)

      REM Normalise
      d%% = 1E18 DIV (b%%(V%)+1) : t%%() = 1,d%%
      PROC_bigumul(u%%(),a%%(),t%%(),U%,1)
      PROC_bigumul(v%%(),b%%(),t%%(),V%,1)
      V% = v%%(0)
      U% = u%%(0) - V%
      IF U% < DIM(q%%(),1) U% += 1
      q%%(0) = U%

      FOR J% = U% TO 1 STEP -1
        q%%(J%) = (u%%(J%+V%) * 1E18+u%%(J%+V%-1)) / v%%(V%)
        l%% = q%%(J%)-256 : h%% = q%%(J%)+256
        REPEAT
          q%%(J%) = (l%%+h%%) DIV 2
          PROC_bigumul(t%%(),q%%(J%-1 TO J%),v%%(),1,V%)
          IF FN_bigucmp(t%%(),u%%(J%-1 TO J%+V%),t%%(0),V%+1) > 0 h%% = q%%(J%) ELSE l%% = q%%(J%)
        UNTIL h%%-l%% <= 1 AND l%% = q%%(J%)
        PROC_bigusub(t%%(),u%%(J%-1 TO J%+V%),t%%(),V%+1,t%%(0))
        u%%(J% TO J%+V%) = t%%(1 TO V%+1)
      NEXT J%

      WHILE q%%(q%%(0)) = 0 AND q%%(0) > 0 q%%(0) -= 1 : ENDWHILE
      ENDPROC

Because it's so slow I asked DeepSeek to review the code, and it was not very complimentary.

But on the other hand none of its suggestions for improvements actually worked, so on that basis I'm not going to take its criticism too seriously. Perhaps a human (with a brain working better than mine is) can spot a way of speeding it up which DeepSeek failed to.