Why Avoid BASIC on RPi?

[Heater]

DavidS,

Operator overloading is doable in any language that supports user defined data types, including C

Not in C.

An example of how it is done in some C compilers that bend the rules a little (by allowing operator overloading) is below

Well, they are not C compilers then are they.

Out of interest what compilers do you mean and where can I get them to play with?

Remember that we now have inline functions in C,

Inlining is neither here nor there in this discussion. That is an optimization a compiler may or may not perform.

I have used a few C compilers that allow this,... you can do the exact same thing in C++, and it works (leave out the OO stuff completely).

Again, I say, they were not C compilers. And I'd be interested to know exactly what they were and how to get one to play with.

Now, looking at your example, it seems to me that much the same could be done in C++. But it does not save you any effort.

Instead of "typedef struct {...} PIX" we have "class PIX {...}"

Instead of free standing overload functions with two pointers like "PIX operator+(PIX *src, PIX *trans){...}" we put them into the class and only have one parameter, the other being the "this" pointer of the instance itself, like so: "PIX operator+(PIX *trans){...}"

One could write much the same as you present in C++ but I see no point to it. Either way the actual work you have to write to implement your data type still has to be done.

[ejolson]

With another 300 lines of code behind the operator overloading of "+", "-", "*", "=" etc...

I've always found programming C++ objects and operator overloading quite like building with TinkerToys.



It's great fun and everything seems to be going splendidly until you try to use the thing you've built.

[DavidS]

@Heater:
Actually the example I gave would compile in C++, given that it had some more substance it could be usefull.

It is true that it is non-standard in C, as I already said that is the one thing that C++ has that I would like to see in C as standard. Also remember that there was descusion of including this in C11 though for an unknown reason it was discarded, usually when such things are discused it is because it is already used by enough C implementations.

And to answer the question of what C compilers I know of that do this:

• AnderC
• OCC (Orange C Compiler)
• ExCC

I am sure there are more. Though you can also do it simply by compiling with a C++ compiler. Also the syntax I showed is much more usable than embeding it as a method, and more readable. In C++ it is the way to do it in my opinion.

Of course if you really want to talk about OO, take a look at BOOPSI which is a means of implementing OO purely in C in a way that is much more readable than C++.

I would not hold the use of C++ against you if you did not use classes, class methods, and/or new/delete. If you were just using it to get operator overloading it would be reasonable and an exception I would be willing to make.

Though remember that the challenge was for you to do this in C11, without any extensions or non-standard libs. Then I would convert it to the BASIC of my choice.

[Heater]

@ejolson

I've always found programming C++ objects and operator overloading ... great fun and everything seems to be going splendidly until you try to use the thing you've built.

Ha! That is so true.

[jahboater]

@ejolson

I've always found programming C++ objects and operator overloading ... great fun and everything seems to be going splendidly until you try to use the thing you've built.

Ha! That is so true.

Or until some one else tries to maintain your project.

[Heater]

DavidS,

Certainly one can do operator overloading in C++ without using classes and methods. Just use free standing "operator" functions. I don't know why you think that is better. You still need a struct to hold your data but now you need two parameters to your overload functions instead of just one. How is that preferable?

BOOPSI is interesting but if I really want to talk OPP I don't want to be using C and a bunch of macros.

I would not hold the use of C++ against you if you did not use classes, class methods, and/or new/delete. If you were just using it to get operator overloading it would be reasonable and an exception I would be willing to make.

I'm sicking to my class and methods solution. I see no benefit to splitting those methods out into free functions and having to add the second parameter to all of them.

My intention here is to keep it simple such that stripping it down to regular C would not be too onerous if one wanted to. As for using new/delete, it makes no odds, it's a trivial matter to replace:

Code: Select all

    value = new int64_t[width];

with:

Code: Select all

    value = (int64_t*)malloc(sizeof value[0] * width);

and replace "delete[]" with "free". I do not see any benefit in doing so though.

Also, what is actually required is to not actually do all those new/deletes or malloc/frees. That is horribly inefficient. What I want to do there is replace that with my own memory allocation. Something along the lines of getting a big memory allocation at start up then returning pointers into that memory space in the constructors as the algorithm proceeds. Nothing would actually get freed until the end when it all gets freed in one go.

Though remember that the challenge was for you to do this in C11, without any extensions or non-standard libs. Then I would convert it to the BASIC of my choice.

Yes, sorry about that. Given that ejolson has already presented a fine solution in C I'm probably not going to be stripping my version down to C. There seems little point.

I have not looked at his code yet, I wanted to get my effort working first without peeking, but I'm sure it is very tight and fast code, far more efficient than mine. That should already be a good working model for you to develop a BASIC version from.

I think I mentioned earlier, the promise of C++ is "zero cost abstractions". So my problem now is to keep the abstractions of operating overloading, being able to write "r = x + y" instead of "sum(x, y, r)" or "r = sum(x,y)", but get rid of all the overheads involved.

[jahboater]

Also, what is actually required is to not actually do all those new/deletes or malloc/frees. That is horribly inefficient. What I want to do there is replace that with my own memory allocation. Something along the lines of getting a big memory allocation at start up then returning pointers into that memory space in the constructors as the algorithm proceeds. Nothing would actually get freed until the end when it all gets freed in one go.

Ejolson's C solution does just that.

Note that above 128KB, memory allocation is different, it does a mmap() of /dev/zero and the large block ends up in its own segment (possibly not subject to ulimit).

[ejolson]

I have not looked at his code yet, I wanted to get my effort working first without peeking, but I'm sure it is very tight and fast code, far more efficient than mine. That should already be a good working model for you to develop a BASIC version from.

Thanks for your vote of confidence. The efficiency of my code seems to suffer rather badly from having too many 64-bit divisions, especially on the ARMv6 architecture. To avoid this I could have worked in a power-of-2 base or used 32-bit arithmetic throughout. Another place for improvement might be tuning the switch-over point which determines whether the Karatsuba algorithm or O(n^2) multiplication is used. The value currently used in the code was reasonable for Intel PIII class machines but may not be for Raspberry Pi. Anyway, these are just a few pitfalls you may want to think about when writing a new code.

[Heater]

ejolson,

I'm sure it could be done with a power-of-2 base or the full width of 32/64 bit integers. Problem then is getting a decimal result out. That is something I don't want to fight with and the solutions we have in C, Python and JS all take 10 or more times longer to do the hex to decimal conversion than get the result in the first place.

I leaning towards not bothering with swicthing to "naive" multiply, kind of like the elegance of it as it stands.

Just now got my bint multiply working! Here is the top level mul() method:

Code: Select all

bint bint::mul (bint& a)
{
    // Ensure operands are same width.
    while (this->width > a.width)
    {
        a.grow();
    }
    while (a.width > this->width)
    {
        this->grow();
    }

    // Result will be twice as wide as the operands.
    bint result(a.width * 2);

    // The base case, only one element in value, just do the multiply
    if (width == 1)
    {
        int64_t product = value[0] * a.value[0];
        if (product < BASE) {
            result.value[0] = product;
            result.value[1] = 0;
        }
        else
        {
            result.value[0] = (product % BASE);
            result.value[1] = (product / BASE);
        }
        return result; 
    }

    // Calculates the size of the numbers
    int m = (this->width);
    int m2 = m / 2;

    // Split the numbers in the middle
    bint high1 = this->high();
    bint low1 = this->low();
    bint high2 = a.high();
    bint low2 = a.low();

    // Do da karatsaba shuffle, yabba dabba do. 
    bint z0 = low1 * low2;
    bint z2 = high1 * high2;
    bint s1 = low1 + high1;
    bint s2 = low2 + high2;
    bint z1 = s1 * s2;
    bint t1 = z1 - z2 - z0;
    bint t1Shifted = t1.shift1(m2);
    bint z2Shifted = z2.shift2(m);
    result = z2Shifted + t1Shifted + z0;
    result.shrink();
    return result; 
}

[ejolson]

I'm sure it could be done with a power-of-2 base or the full width of 32/64 bit integers. Problem then is getting a decimal result out. That is something I don't want to fight with and the solutions we have in C, Python and JS all take 10 or more times longer to do the hex to decimal conversion than get the result in the first place.

I leaning towards not bothering with swicthing to "naive" multiply, kind of like the elegance of it as it stands.

The Karatsuba algorithm looks nice after the operator overloading.

The profiling done in this post indicates from 22 to 66 percent of the run time is spent in the library function __udivmoddi4 for my code depending on the exact ARM architecture and compiler used. I wonder if the profiling will look similar for your code. I notice you included a nice optimization that avoids the subsequent divisions when doing a multiply by zero.

[DavidS]

The Karatsuba algorithm looks nice after the operator overloading.

The profiling done in thispost indicates from 22 to 66 percent of the run time is spent in the library function __udivmoddi4 for my code depending on the exact ARM architecture and compiler used. I wonder if the profiling will look similar for your code.

Interesting, makes me wonder if there would be a significant improvment in performance from converting to using 32-bit operations (which has been my primarry stumbling block so far).

I must admit that this is the most difficulty I have ever had converting C code into ARM BASIC.

Could be worth translating back to C once I get the BASIC version done. I am doing it in BBC BASIC V (ARM BASIC), as it is the BASIC of RISC OS, and an extended subset of it is the most widely available dialect of BASIC (Brandy BASIC, and others, for DOS, Linux, Windows, Macintosh System Software, Amiga OS, Android, etc). I figure so long as I stick with the subset of the langage that is implemented on all platforms that have versions based on BBC BASIC V, it should be able to be ran on any system. Though that only leaves out SWI calls and *commands (shell commands in RISC OS) being used, the rest of the language is implemented on other platforms.

To bad most systems do not have a compiler for there implementations of BBC BASIC V work alikes. I only know of the compiler for RISC OS.

Unfortunately for the counter challenge that I am working on, that is a comon case that I feel BASIC is more expresive, not using SWI's means no WIMP front end and no drag and drop for file open/save operations (so can not do that the RISC OS way, have to use the command line way).

[jahboater]

The profiling done in this post indicates from 22 to 66 percent of the run time is spent in the library function __udivmoddi4 for my code depending on the exact ARM architecture and compiler used. I wonder if the profiling will look similar for your code. I notice you included a nice optimization that avoids the subsequent divisions when doing a multiply by zero.

FYI, I think __udivmoddi4 is 128 bit division - returning two 64-bit results.

[ejolson]

FYI, I think __udivmoddi4 is 128 bit division - returning two 64-bit results.

You might be right about that. I thought it was 64-bit division.

According to the Wookieepedia

Han Solo: "Threepio was translating Artoo's beeps into a language we could understand—good old Basic."

but not all enoyed Galactic Basic as shown by

Corran Horn: "Do any of you speak Basic?"
Kotaa Zun-qin: "I speak your infidel tongue. It tastes like the waste excretions of an ill vhlor on my tongue, but I can speak it."

Anyway, given the above Karatsuba C++ code that is missing everything else, here is a Fibonacci program written in FreeBASIC that is only missing the Karatsuba algorithm.

Code: Select all

rem fibo.bas -- Compute the nth Fibonacci Number
rem Written December 6, 2018 by Eric Olson
rem
rem This program demonstrates the expressiveness of FreeBASIC as 
rem measured by explicitly coding big-number arithmetic and then
rem using the doubling formulas
rem
rem     F(2k) = F(k)[2F(k+1)-F(k)]
rem   F(2k+1) = F(k+1)^2+F(k)^2
rem   F(2k+2) = F(k+1)[2F(k)+F(k+1)]
rem
rem to compute the nth Fibonacci number.  For simplicity this code
rem does not use Karatsuba multiplication but instead relies on
rem the O(n^2) algorithm in all cases.  When computing F(4784969)
rem this results in significantly more computational work.

const as integer nmax = 250100
const as integer radix = 10000
const as integer rexp = 4

type bignum
    dim as integer n
    dim as integer d(1 to nmax)
end type

sub bigprint(x as bignum)
    dim as integer i
    if x.n=0 then
        print 0
    else
        for i=x.n to 1 step -1
            dim as string s=str(x.d(i))
            if i<x.n then
                s=mid("0000000000",1,rexp-len(s))+s
            end if
            print s;
        next i
        print
    endif
end sub

sub bigmul(z as bignum,x as bignum,y as bignum)
    dim as integer i,j
    for i=1 to x.n+y.n
        z.d(i)=0
    next i
    for i=1 to x.n
        for j=1 to y.n
            dim as integer t,c
            c=0
            t=z.d(i+j-1)+x.d(i)*y.d(j)
            if t>=radix then
                c=t\radix
                t=t mod radix
            end if
            z.d(i+j-1)=t
            z.d(i+j)=z.d(i+j)+c
        next j
    next i
    z.n=x.n+y.n
    while z.d(z.n)=0 and z.n>1
        z.n=z.n-1
    wend
    return
end sub

sub bigadd(z as bignum,x as bignum,y as bignum)
    dim as integer i,max
    if y.n>x.n then
        max=y.n
    else
        max=x.n
    end if
    for i=1 to max+1
        z.d(i)=0
    next i
    for i=1 to max
        dim as integer t,c
        c=0
        if i<=x.n then
            t=z.d(i)+x.d(i)
            if i<=y.n then
                t=t+y.d(i)
            end if
        else
            t=z.d(i)+y.d(i)
        end if
        if t>=radix then
            c=1
            t=t-radix
        end if
        z.d(i)=t
        z.d(i+1)=z.d(i+1)+c
    next i
    z.n=max+1
    while z.d(z.n)=0 and z.n>1
        z.n=z.n-1
    wend
    return
end sub

sub bigsub(z as bignum,x as bignum,y as bignum)
    dim as integer i
    for i=1 to x.n
        z.d(i)=0
    next i
    for i=1 to y.n
        dim as integer t,c
        c=0
        t=z.d(i)+x.d(i)-y.d(i)
        if t<0 then
            t=t+radix
            c=-1
        end if
        z.d(i)=t
        z.d(i+1)=z.d(i+1)+c
    next i
    for i=y.n+1 to x.n
        dim as integer t,c
        c=0
        t=z.d(i)+x.d(i)
        if t<0 then
            t=t+radix
            c=-1
        end if
        z.d(i)=t
        z.d(i+1)=z.d(i+1)+c
    next i
    z.n=x.n
    while z.d(z.n)=0 and z.n>1
        z.n=z.n-1
    wend
    return
end sub

sub fibo(n as integer,a as bignum,b as bignum)
    if n=0 then
        a.n=0:b.n=1:b.d(1)=1
        return
    end if
    fibo(n\2,a,b)
    static as bignum t1,t2,t3
    bigmul(t1,a,a):bigmul(t2,b,b):bigadd(t3,t1,t2)
    if n mod 2=0 then
        rem [a,b]=[a*(b*2-a),a*a+b*b]
        bigadd(t1,b,b):bigsub(t2,t1,a):bigmul(t1,a,t2)
        a=t1:b=t3
    else
        rem [a,b]=[a*a+b*b,b*(2*a+b)]
        bigadd(t1,a,a):bigadd(t2,t1,b):bigmul(t1,b,t2)
        a=t3:b=t1
    end if
    return
end sub

rem main routine
dim as integer n=4784969
if n<2 then
    print str(n)
else
    static as bignum a,b
    fibo(n-1,a,b)
    bigprint(b)
end if

As this program is using the O(n^2) multiplication throughout it runs noticeably slower than codes using asymptotically superior algorithms. Run times on the original Pi 2B are as follows:

Code: Select all

$ time ./fibo | head -c32; time ./fibo | tail -c32
10727395641800477229364813596225
real    36m41.437s
user    36m40.874s
sys 0m0.149s
4856539211500699706378405156269

real    36m14.462s
user    36m11.622s
sys 0m3.141s

In my opinion, FreeBASIC is better than using Cfront to code in C++ and just as much fun as using Bourne's macros to turn C into Algol.

To get FreeBASIC running on the Raspberry Pi, I downloaded the most recent linux-armv7a-hf-debian binary from here, made the directory /usr/local/fbc-299 for safe keeping and unzipped the archive. After installing libncurses5-dev I could compile with the command

Code: Select all

$ /usr/local/fbc-299/bin/fbc -O 3 fibo.bas

[DavidS]

Anyway, given the above Karatsuba C++ code that is missing everything else, here is a Fibonacci program written in FreeBASIC that is only missing the Karatsuba algorithm.

Code: Select all

rem fibo.bas -- Compute the nth Fibonacci Number
rem Written December 6, 2018 by Eric Olson
rem
rem This program demonstrates the expressiveness of FreeBASIC as 
rem measured by explicitly coding big-number arithmetic and then
rem using the doubling formulas
rem
rem     F(2k) = F(k)[2F(k+1)-F(k)]
rem   F(2k+1) = F(k+1)^2+F(k)^2
rem   F(2k+2) = F(k+1)[2F(k)+F(k+1)]
rem
rem to compute the nth Fibonacci number.  For simplicity this code
rem does not use Karatsuba multiplication but instead relies on
rem the O(n^2) algorithm in all cases.  When computing F(4784969)
rem this resutls in significantly more computation work.

const as integer nmax = 250100
const as integer radix = 10000
const as integer rexp = 4

type bignum
    dim as integer n
    dim as integer d(1 to nmax)
end type

sub bigprint(x as bignum)
    dim as integer i
    if x.n=0 then
        print 0
    else
        for i=x.n to 1 step -1
            dim as string s=str(x.d(i))
            if i<x.n then
                s=mid("0000000000",1,rexp-len(s))+s
            end if
            print s;
        next i
        print
    endif
end sub

sub bigmul(z as bignum,x as bignum,y as bignum)
    dim as integer i,j
    for i=1 to x.n+y.n
        z.d(i)=0
    next i
    for i=1 to x.n
        for j=1 to y.n
            dim as integer t,c
            c=0
            t=z.d(i+j-1)+x.d(i)*y.d(j)
            if t>=radix then
                c=t\radix
                t=t mod radix
            end if
            z.d(i+j-1)=t
            z.d(i+j)=z.d(i+j)+c
        next j
    next i
    z.n=x.n+y.n
    while z.d(z.n)=0 and z.n>1
        z.n=z.n-1
    wend
    return
end sub

sub bigadd(z as bignum,x as bignum,y as bignum)
    dim as integer i,max
    if y.n>x.n then
        max=y.n
    else
        max=x.n
    end if
    for i=1 to max+1
        z.d(i)=0
    next i
    for i=1 to max
        dim as integer t,c
        c=0
        if i<=x.n then
            t=z.d(i)+x.d(i)
            if i<=y.n then
                t=t+y.d(i)
            end if
        else
            t=z.d(i)+y.d(i)
        end if
        if t>=radix then
            c=1
            t=t-radix
        end if
        z.d(i)=t
        z.d(i+1)=z.d(i+1)+c
    next i
    z.n=max+1
    while z.d(z.n)=0 and z.n>1
        z.n=z.n-1
    wend
    return
end sub

sub bigsub(z as bignum,x as bignum,y as bignum)
    dim as integer i
    for i=1 to x.n
        z.d(i)=0
    next i
    for i=1 to y.n
        dim as integer t,c
        c=0
        t=z.d(i)+x.d(i)-y.d(i)
        if t<0 then
            t=t+radix
            c=-1
        end if
        z.d(i)=t
        z.d(i+1)=z.d(i+1)+c
    next i
    for i=y.n+1 to x.n
        dim as integer t,c
        c=0
        t=z.d(i)+x.d(i)
        if t<0 then
            t=t+radix
            c=-1
        end if
        z.d(i)=t
        z.d(i+1)=z.d(i+1)+c
    next i
    z.n=x.n
    while z.d(z.n)=0 and z.n>1
        z.n=z.n-1
    wend
    return
end sub

sub fibo(n as integer,a as bignum,b as bignum)
    if n=0 then
        a.n=0:b.n=1:b.d(1)=1
        return
    end if
    fibo(n\2,a,b)
    static as bignum t1,t2,t3
    bigmul(t1,a,a):bigmul(t2,b,b):bigadd(t3,t1,t2)
    if n mod 2=0 then
        rem [a,b]=[a*(b*2-a),a*a+b*b]
        bigadd(t1,b,b):bigsub(t2,t1,a):bigmul(t1,a,t2)
        a=t1:b=t3
    else
        rem [a,b]=[a*a+b*b,b*(2*a+b)]
        bigadd(t1,a,a):bigadd(t2,t1,b):bigmul(t1,b,t2)
        a=t3:b=t1
    end if
    return
end sub

rem main routine
dim as integer n=4784969
if n<2 then
    print str(n)
else
    static as bignum a,b
    fibo(n-1,a,b)
    bigprint(b)
end if

As this program is using the O(n^2) multiplication throughout it runs noticeably slower than codes using asymptotically superior algorithms. Run times on the original Pi 2B are as follows:

Code: Select all

$ time ./fibo | head -c32; time ./fibo | tail -c32
10727395641800477229364813596225
real    36m41.437s
user    36m40.874s
sys 0m0.149s
4856539211500699706378405156269

real    36m14.462s
user    36m11.622s
sys 0m3.141s

In my opinion, FreeBASIC is better than using Cfront to code in C++ and just as much fun as using Bourne's macros to turn C into Algol.

Well if you are going for FreeBASIC (probably better anyway) then you can use 64-bit integers, the datatype is LONGINT or ULONGINT (signed and unsigned respectively).

see:
https://www.freebasic.net/wiki/wikka.ph ... dDataTypes

By the way you just got me to boot into Linux on my new RPi 3B+, as I already have FreeBASIC there, now to re-translate the original using 64-bit integers in FreeBASIC. Thank you for talking me into doing this in Linux.

FreeBASIC is a bit closer to ANSI Extended BASIC than BBC BASIC is (using the SUB and FUNCTION keywords for example).

Thank you for not using classes in your FreeBASIC Code.

[Paeryn]

FYI, I think __udivmoddi4 is 128 bit division - returning two 64-bit results.

No, it is a combined Unsigned DIV and MOD on Double-Integers (long long, typically 64-bit). It takes two 64-bit integers and a pointer to a 64-bit integer. It returns the quotient, the remainder is stored at the address pointed to by the third parameter (if not NULL). Not totally sure what the 4 stands for without looking, I thought it might refer to the bytes in an integer but I know there are similar functions ending in 3. Well, after looking at the source I'm no wiser to what the 4 refers to. It's an internal function anyway for gcc's use.

[ejolson]

FYI, I think __udivmoddi4 is 128 bit division - returning two 64-bit results.

No, it is a combined Unsigned DIV and MOD on Double-Integers (long long, 64-bit). It takes two 64-bit integers and returns both the division and the remainder. Not totally sure what the 4 stands for without looking, I would have said 4 refers to the bytes in an integer but I know there are similar functions ending in 3.

Maybe the 4 indicates that it makes the Fibonacci program run four times slower than expected. Then a base-2 solution that avoids division except when printing should run four times faster.

[DavidS]

Well I have an initial conversion to ARM BASIS, not working yet and my brain is not catching what is likely a very obvious reason. I was just running it to find out if I missed a variable name or some such the first error that comes up is a "number too big" error in the procedure PROCbigmul.

I am posting the conversion thus far not working as I am having trouble concentrating at the moment.

Code: Select all

 REM >ejfiboslow
ON ERROR PRINT "LINE : " + STR$(ERL) + " ERROR : " + REPORT$ : E

REM fibo.bas -- Compute the nth Fibonacci Number
REM Written December 6, 2018 by Eric Olson
REM Translated to ARM BASIC by David Cagle, December 7th 2018.
REM
REM This program demonstrates the expressiveness of FreeBASIC as
REM measured by explicitly coding big-number arithmetic and then
REM using the doubling formulas
REM
REM     F(2k) = F(k)[2F(k+1)-F(k)]
REM   F(2k+1) = F(k+1)^2+F(k)^2
REM   F(2k+2) = F(k+1)[2F(k)+F(k+1)]
REM
REM to compute the nth Fibonacci number.  For simplicity this co
REM does not use Karatsuba multiplication but instead relies on
REM the O(n^2) algorithm in all cases.  When computing F(4784969
REM this resutls in significantly more computation work.

nmax% = 250100
radix% = 10000
rexp% = 4

REM type bignum
REM    dim as integer n
REM    dim as integer d(1 to nmax)
REM end type

REM main routine
n%=4784969
IF n%<2 THEN
    PRINT n%
ELSE
    REM static as bignum a,b
    REM First word (4 bytes) is as n% member of type.
    REM Member d% of FreeBASIC type is starts after the first
    REM word.
    DIM an% nmax%*4+4 :ad%=an%+4: a%=an%
    DIM bn% nmax%*4+4 :bd%=bn%+4: b%=bn%

    REM static as bignum t1,t2,t3
    DIM t1n% nmax%*4+4 :t1d%=t1n%+4
    DIM t2n% nmax%*4+4 :t2d%=t2n%+4
    DIM t3n% nmax%*4+4 :t3d%=t3n%+4


    PROCfibo(n%-1)
    PROCbigprint(b%)
ENDIF
END


DEFPROCbigprint(xn%)
    LOCAL xd%, i%
    xd%=x%+4

    IF xn%=0 THEN
        PRINT 0
    ELSE
        FOR i%=!xn%<<2 TO 0 STEP -4
            s$=STR(xd%!i%)
            IF i%<!xn% THEN
                s$=MID$("0000000000",1,rexp%-LEN(s$))+s$
            ENDIF
            PRINT s$;
        NEXT
        PRINT
    ENDIF
ENDPROC


DEFPROCbigmul(zn%,xn%,yn%)
    LOCAL zd%, xd%, yd%, i%, j%, t%, c%
    zd%=zn%+4: xd%=xn%+4: yd%=yn%+4

    FOR i%=0 TO (!xn%+!yn%)<<2 STEP 4
        zd%!i% = 0
    NEXT
    FOR i%=0 TO !xn%<<2 STEP 4
        FOR j%=1 TO !yn%<<2 STEP 4
            c% = 0
            t% = zd%!(i%+j%-4) + xd%!i% * yd%!j%
            IF t% >= radix% THEN
                c% = t% DIV radix%
                t% = t% MOD radix%
            ENDIF
            zd%!(i%+j%-4) = t%
            zd%!(i%+j%)=zd%!(i%+j%)+c%
        NEXT
    NEXT
    !zn% = !xn% + !yn%
    WHILE zd%!(!zn%<<2) = 0 AND !zn%
        !zn% -= 1
    WEND
ENDPROC


DEFPROCbigadd(zn%,xn%,yn%)
    LOCAL zd%, xd%, yd%, i%, max%, t%, c%
    zd%=zn%+4: xd%=xn%+4: yd%=yn%+4

    IF !yn%>!xn% THEN
        max% = !yn% << 2
    ELSE
        max% = !xn% << 2
    ENDIF
    FOR i%=1 TO nmax%+4 STEP 4
        zd%!i = 0
    NEXT
    FOR i%=1 TO nmax% STEP 4
        c% = 0
        IF i% <= !xn% << 2 THEN
            t% = zd%!i% + xd%!i%
            IF i% <= !yn% << 2 THEN
                t% += yd%!i%
            ENDIF
        ELSE
            t% = zd%!i% + yd%!i%
        ENDIF
        IF t% >= radix% THEN
            c% = 1
            t% -= radix%
        ENDIF
        zd%!i = t%
        zd%!(i+1) += c%
    NEXT
    !zn% = (nmax% + 4) >> 2
    WHILE zd%!(!zn%<<2) = 0 AND !zn%
        !zn% -= 1
    ENDWILE
ENDPROC


DEFPROCbigsub(zn%,xn%,yn%)
    LOCAL zd%, xd%, yd%, i%, max%, t%, c%
    zd%=zn%+4: xd%=xn%+4: yd%=yn%+4

    FOR i%=0 TO !xn%<<2
        zd%!i% = 0
    NEXT
    FOR i%=0 TO !yn%<<2 STEP 4
        c% = 0
        t% = zd%!i% + xd%!i% - yd%!i%
        IF t%<0 THEN
            t% += radix%
            c% = -1
        ENDIF
        zd%!i% = t%
        zd%!(i%+4) = zd%!(i%+4) + c
    NEXT
    FOR i%=(!yn%+1)<<2 TO !xn%<<2 STEP 4
        c% = 0
        t% = zd%!i% + xd%!i%
        IF t% < 0 THEN
            t% += radix%
            c% = -1
        ENDIf
        zd%!i% = t%
        zd%!(i%+4) = zd%!(i%+4) + c%
    NEXT
    !zn% = !xn%
    WHILE zd%!(!zn% << 2) = 0 AND !zn%
        !zn% -= 1
    ENDWHILE
ENDPROC


DEFPROCfibo(nv%)
  IF nv% THEN
    PROCfibo(nv% DIV 2)
    REM static as bignum t1,t2,t3
    PROCbigmul(t1n%,an%,an%):PROCbigmul(t2,bn%,bn%)
                  PROCbigadd(t3n%,t1n%,t2n%)
    IF n% MOD 2 = 0 THEN
        REM [a,b]=[a*(b*2-a),a*a+b*b]
        PROCbigadd(t1n%,bn%,bn%):PROCbigsub(t2n%,t1n%,an%)
                  PROCbigmul(t1n%,an%,t2m%)
        an%=t1n% : ad%=t1d% : bn%=t3n% : bd%=t3d%
    ELSE
        REM [a,b]=[a*a+b*b,b*(2*a+b)]
        PROCbigadd(t1n%,an%,an%):PROCbigadd(t2n%,t1n%,bn%)
                  PROCbigmul(t1n%,bn%,t2n%)
        an%=t3n% : ad%=t3d% : bn%=t1n% : bd%=t3d%
    ENDIF
  ENDIF
    !an% = 0 : !bn% = 1 : !bd% = 1
ENDPROC

Note that i am using indirected memory blocks instead of arrays, hence the scaling of acceses to by 4 times (integers are 4 bytes each, indirection in ARM basic is always to the byte address).

FYI the line where the error is reported is in bigmul and is th line:
t% = zd%!(i%+j%-4) + xd%!i% * yd%!j%

[DavidS]

And the error is definitely on the first call to PROCbigmul, as I just found an error later in that proc that would prevent it from continuing to the end of the PROC.

[ejolson]

And the error is definitely on the first call to PROCbigmul, as I just found an error later in that proc that would prevent it from continuing to the end of the PROC.

What happens if you set nmax to 250 and try to compute F(11) instead?

[jahboater]

No, it is a combined Unsigned DIV and MOD on Double-Integers (long long, typically 64-bit). It takes two 64-bit integers and a pointer to a 64-bit integer. It returns the quotient, the remainder is stored at the address pointed to by the third parameter (if not NULL).

What are the two 64-bit integers it takes as input? I presumed they were the 128-bit dividend.
I guess that is 64-bit division, in the same way that a 64-bit multiply returns a 128-bit result.

[jahboater]

Off topic, but I see the forthcoming release of GCC will have D included ...
https://dlang.org/overview.html

[Heater]

I am floundering with my DIY karatsuba as well....I begin to suspect there is something about this algorithm I don't understand.

I'm trying to write my karatsuba from the description and pseudo code on wikipedia and the top level of that is shown above.

All works well much of the time. It can multiply big numbers, for example it can calculate 2 to the power 1000. With the following result:

Code: Select all

10715086071862673209484250490600018105614048117055336074437503883703510511249361224931983788156958581275946729175531468251871452856923140435984577574698574803934567774824230985421074605062371141877954182153046474983581941267398767559165543946077062914571196477686542167660429831652624386837205668069376
[/quote]
Which Python tells me is correct.

Problem comes, for example, when I try to square that result. Like so:
[code]
    bint two = "2";
    bint res1 = "1";
    for (int i = 0; i < 1000; i++)
    {
        res1 = res1.mul(two);
        res1.shrink(32);
    }
    std::cout << "res1: " << res1 << std::endl;

    bint res2 = res1 * res1;
    std::cout << "res2: " << res2 << std::endl;

Then I get an assertion failure as the calculation of (z1 - z2 - z0) goes negative. I have no support for negative numbers and had assumed that would not happen.

In general terms without getting into the details of my big integer add, sub, etc, a couple of questions:

1) Is that subtraction expected to go negative sometimes, and I should handle it, or does this indicate I have something else wrong somewhere?

2) Unlike the wiki pseudo code I decided to make all my numbers a power of two in length rather than split them into potentially odd length parts. My reasoning being that this would mean the base case was always just a simple 1 element by 1 element multiply. And well, because, powers of two, I though this might make splitting easier and memory management later.

I now suspect this was a bad idea. It likely means using more memory and doing more work than keeping the numbers whatever width they should be. Perhaps it also leads to my problem with negative intermediate results.

Anyone have any thoughts on this?

I could post the whole code here but it I'm guessing it would be two long and tedious for anyone to be bothered to read.

My numbers are working in base 1000000000000000000 by the way. 18 decimal digits per element.

[Heater]

jahboater,

Interesting. From your link: "What is D?"

"...It doesn't come with a VM, a religion, or an overriding philosophy....

What? A programming language without a religion? Cannot be. It's doomed to failure

Perhaps it will get a religious following later.

[jahboater]

What? A programming language without a religion? Cannot be. It's doomed to failure

It does seem appealing. It appears to offer both safety and C like speed and systems programming ability.
Perhaps a better evolutionary step after C than C++ is.

The section in the overview "Who is D for" starts with:-

Programmers who routinely use lint or similar code analysis tools to eliminate bugs before the code is even compiled.
People who compile with maximum warning levels turned on and who instruct the compiler to treat warnings as errors.

[Heater]

jahboater,

Programmers who routinely use lint or similar code analysis tools...People who compile with maximum warning levels turned on...

Yes, sounds like just the thing for me. Can't live without asan and friends!

Of course the next language after "C" should rightfully be "D". Clearly C++ is is only a prototype, an experiment, a detour on the road between between C and D. (Up a blind alley as it happens).