Introduction to BBC BASIC

[DavidS]

DavidS,

I am going to see about giving it a try in BBC BASIC on RISC OS (hopefully not to much work). Though where is the library with PROCbigmul and PROCbigadd? Or did you not get far enogh for those two to work?

See Richard's post here: https://www.raspberrypi.org/forums/view ... 7#p1500070

Did you get a result from the big fibo BASIC code that linked you to (Also see my previous post)? We know that one works well.

Unfortunately that one is not BBC BASIC text format. The line numbers are not the first thing on the line, and the keywords are lower case, so I have been spending a few minutes here and there translating it to BBC BASIC. Why did you not include the tokenized form, then it would have been correct.

BBC BASIC is case sensitive, as well as type local with everything. thus for example you have the keyword TO and you can have the real variables to, To, and tO, at the same time as you have the integer variables to%, To%, and tO%, and the string variables to$, To$, and tO$, so the listing you provided having the case wrong is a big deal. Taking care of the line number issue was easy.

And that post of Richard only stated that he did not finish his bigint library (so he felt not good enough for entry), it did not say that the two functions used therein do or do not work.

[DavidS]

Though where is the library with PROCbigmul and PROCbigadd? Or did you not get far enogh for those two to work?

As I've said, I developed the library only far enough to allow the Fibo application to run (which it does, so obviously those procedures work!). Specifically there is no support for negetive numbers nor for division, which are pretty serious limitations for a 'bigint' library! I'm happy to publish the library in its current, incomplete, state but I'm not sure that there would be much point; it makes heavy use of 64-bit integers and array pointers, so is not at all portable to ARM BASIC V.

Here's a snippet of the code:

Code: Select all

      REM Unsigned addition
      DEF PROCbiguadd(c%%, a%%, b%%)
      LOCAL A%, F%, r%%, s%%, a%%(), b%%(), c%%() : r%% = 10^18 : s%% = !c%% AND &7FFFFF
      IF a%%!4 < b%%!4 SWAP a%%,b%%
      A% = a%%!4 : a%%!4 = b%%!4
      c%%!4 = a%%!4 : IF c%%!4 > s%% ERROR 100, "PROCbiguadd: BigInt overflow"
      IF INKEY(-256)=&73 IF @platform% AND &40 THEN ]^a%%() = a%%+3 : ]^b%%() = b%%+3 : ]^c%%() = c%%+3 ELSE !^a%%() = a%%+3 : !^b%%() = b%%+3 : !^c%%() = c%%+3
      F% = FNbigadd(c%%(), a%%(), b%%()) <> 0
      a%%!4 = A%
      WHILE F% OR c%%!4 < A%
        IF c%%!4 >= s%% ERROR 100, "PROCbiguadd: BigInt overflow"
        c%%!4 += 1
        IF c%%!4 <= A% c%%(c%%!4 - 1) = a%%(c%%!4 - 1) - F% ELSE c%%(c%%!4 - 1) = - F%
        F% = c%%(c%%!4 - 1) >= r%% : IF F% c%%(c%%!4 - 1) -= r%%
      ENDWHILE
      ENDPROC

      DEF FNbigadd(c%%(), a%%(), b%%())
      LOCAL F%,N%,r%%,d%%() : r%% = 10^18 : N% = DIM(a%%(),1) : DIM d%%(N%)
      c%%() = a%%() + b%%()
      REPEAT
        d%%() = c%%() DIV r%%
        IF SUM(d%%()) = 0 EXIT REPEAT
        c%%() MOD= r%%
        F% OR= d%%(N%)
        IF N% SYS MemMove%%, ^d%%(1), ^d%%(0), N%*8
        d%%(0) = FALSE
        c%%() += d%%()
      UNTIL FALSE
      = F%

I see that the usage of 64-bit integers would be an issue for sure.

though as you know pointers have been around in ARM BASIC for as long as we have had ARM BASIC (though done a bit differently than you show doing it).

I guess I am stuck with using the one that is all wrong (though works correctly) that was provided by request of heater.

[Heater]

DavidS,

Unfortunately that one is not BBC BASIC text format

There is nothing wrong with it. That code runs under RichardRussell's BBC BASIC and a bunch of other BASICs we have tried it with. Your BASIC is, shall we say, dated.

Besides, as far as I can tell BASIC supports celebrate and promote the incompatibilities between BASIC dialects. You should be used to this.

I'm sure it will only take a moment to to massage it into shape for your BASIC. You could use "tr" and/or "sed".

Code: Select all

Why did you not include the tokenized form, then it would have been correct.

It is correct. Tokenized form is not source code. We can only post source code here and I'm not putting build products in a source code repository unless I really have to. That is just wrong.

If you find a confusion in there between "To", "tO", "to" and such like then that is a problem with your BASIC. But I will grant you that I think it is bad programming practice. Do please point them out if you do find any and I'll correct them in the repository.

[ejolson]

Though where is the library with PROCbigmul and PROCbigadd? Or did you not get far enogh for those two to work?

As I've said, I developed the library only far enough to allow the Fibo application to run (which it does, so obviously those procedures work!). Specifically there is no support for negetive numbers nor for division, which are pretty serious limitations for a 'bigint' library! I'm happy to publish the library in its current, incomplete, state but I'm not sure that there would be much point; it makes heavy use of 64-bit integers and array pointers, so is not at all portable to ARM BASIC V.

Here's a snippet of the code:

Code: Select all

      REM Unsigned addition
      DEF PROCbiguadd(c%%, a%%, b%%)
      LOCAL A%, F%, r%%, s%%, a%%(), b%%(), c%%() : r%% = 10^18 : s%% = !c%% AND &7FFFFF
      IF a%%!4 < b%%!4 SWAP a%%,b%%
      A% = a%%!4 : a%%!4 = b%%!4
      c%%!4 = a%%!4 : IF c%%!4 > s%% ERROR 100, "PROCbiguadd: BigInt overflow"
      IF INKEY(-256)=&73 IF @platform% AND &40 THEN ]^a%%() = a%%+3 : ]^b%%() = b%%+3 : ]^c%%() = c%%+3 ELSE !^a%%() = a%%+3 : !^b%%() = b%%+3 : !^c%%() = c%%+3
      F% = FNbigadd(c%%(), a%%(), b%%()) <> 0
      a%%!4 = A%
      WHILE F% OR c%%!4 < A%
        IF c%%!4 >= s%% ERROR 100, "PROCbiguadd: BigInt overflow"
        c%%!4 += 1
        IF c%%!4 <= A% c%%(c%%!4 - 1) = a%%(c%%!4 - 1) - F% ELSE c%%(c%%!4 - 1) = - F%
        F% = c%%(c%%!4 - 1) >= r%% : IF F% c%%(c%%!4 - 1) -= r%%
      ENDWHILE
      ENDPROC

      DEF FNbigadd(c%%(), a%%(), b%%())
      LOCAL F%,N%,r%%,d%%() : r%% = 10^18 : N% = DIM(a%%(),1) : DIM d%%(N%)
      c%%() = a%%() + b%%()
      REPEAT
        d%%() = c%%() DIV r%%
        IF SUM(d%%()) = 0 EXIT REPEAT
        c%%() MOD= r%%
        F% OR= d%%(N%)
        IF N% SYS MemMove%%, ^d%%(1), ^d%%(0), N%*8
        d%%(0) = FALSE
        c%%() += d%%()
      UNTIL FALSE
      = F%

I see that the usage of 64-bit integers would be an issue for sure.

though as you know pointers have been around in ARM BASIC for as long as we have had ARM BASIC (though done a bit differently than you show doing it).

I guess I am stuck with using the one that is all wrong (though works correctly) that was provided by request of heater.

I think there was a version of the same code for 32-bit, because earlier in this thread were posted benchmark scores for both along with the remark that this is one of the few algorithms that ran faster on the C-based 64-bit interpreter compared to the better-optimised 32-bit interpreter.

[DavidS]

DavidS,

Unfortunately that one is not BBC BASIC text format

There is nothing wrong with it. That code runs under RichardRussell's BBC BASIC and a bunch of other BASICs we have tried it with. Your BASIC is, shall we say, dated.

So you are saying that languages that are case sensitive are dated??

So C, C++, Java, JavaScript, etc are all dated by your definition?

It is clearly specified in the manuals of BBC BASIC that it is case sensitive.

Besides, as far as I can tell BASIC supports celebrate and promote the incompatibilities between BASIC dialects. You should be used to this.

I am not as much a fan of BASIC as I am a fan of a particular 3 languages that fall in the catagory of BASIC (BBC BASIC, MS-QuickBASIC [unfortunately died with the Macintosh going 32-bit clean], and MS-AmigaBASIC (very close to MS-QuickBASIC).

I'm sure it will only take a moment to to massage it into shape for your BASIC. You could use "tr" and/or "sed".

Or any number of other tools. Though doing it by hand gives me a chance to fully grock the code, so that I may straiten it out after I give you the run time results of the submitted version.

Code: Select all

Why did you not include the tokenized form, then it would have been correct.

It is correct. Tokenized form is not source code. We can only post source code here and I'm not putting build products in a source code repository unless I really have to. That is just wrong.

That makes no sense. You see even text formats are not standardized. There is a reason that many text editors on many platforms support verious forms of tokenized BASIC (I have actually had a looked at this since we last talked about it). There is a reason it is called Tokenized Source Code.

If you find a confusion in there between "To", "tO", "to" and such like then that is a problem with your BASIC. But I will grant you that I think it is bad programming practice. Do please point them out if you do find any and I'll correct them in the repository.

They are not in your repo. Though they are used in some BBC BASIC programs. And that is the reason that all keywords being upercase is so important.

[RichardRussell]

The line numbers are not the first thing on the line

That has never been a requirement in BBC BASIC. Every implementation of the language I know, including ARM BASIC V, is happy to accept spaces before the line number.

the keywords are lower case,

I agree that Heater should not list a program with lowercase keywords and call it BBC BASIC; at best it's non-standard. But it's easily rectified if you have BBC BASIC for SDL 2.0 installed: select Lowercase Keywords in the Options menu, load the program in, now deselect that option and re-save the program as plain text.

Why did you not include the tokenized form, then it would have been correct.

The tokenised form is not standardised, because it is simply an image of how the program is stored in memory (so for example in Acorn's tokenised format the line numbers are big-endian but in my tokenised format they are little-endian, as they need to be for efficient handling on an x86 CPU). Plain text is the only 'universal' format that everybody can read.

And that post of Richard only stated that he did not finish his bigint library (so he felt not good enough for entry), it did not say that the two functions used therein do or do not work.

It was not necessary to say, since if those functions did not work neither could the Fibo program! I assumed it was evident that the library works sufficiently well to perform the task it was written for.

[Heater]

DavidS,

So you are saying that languages that are case sensitive are dated??

No. I'm saying your BASIC is dated w.r.t the other BASICs. See the Fibo Challenge repository.

Actually I agree with you. I would expect and much prefer languages to be case sensitive. And operating systems. And everything else.

(OMG, I'm agreeing with you and it's on;y Monday!)

Though doing it by hand gives me a chance to fully grock the code, so that I may straiten it out after I give you the run time results of the submitted version.

I like that idea as well.

That makes no sense. You see even text formats are not standardized.

Sure they are. I have never worked in anything that was not ASCII, or one of it's variants for other human languages, or now a days Unicode. Unicode has a few different binary representations but UTF-8 is the thing now a days.

There is a reason that many text editors on many platforms support verious forms of tokenized BASIC

Interesting. I don't recall I have ever used or even seen one. And I have been around. Can you give an example or two?

There is a reason it is called Tokenized Source Code.

Call it what you like. If it's tokens it is not source code.

Except...and I'm not totally sure on this, the Sinclair Spectrum had BASIC keywords on it's keyboard. I always imagined that hitting those keys produced their tokens directly. Never had a Speccy so I don't know for sure. But even then what you see on the screen is source text not tokens.

[Heater]

RichardRussell,

I agree that Heater should not list a program with lowercase keywords and call it BBC BASIC; at best it's non-standard.

Sounds reasonable. I could have a go at fixing it, but to be sure it get's done properly would you mind doing the honors. Pretty please.

[RichardRussell]

I think there was a version of the same code for 32-bit...

Ah, what you're missing is that my 32-bit interpreters still support 64-bit integer variables: you don't need a 64-bit CPU for that! It may seem superficially strange, but it's no different from C which supports the long long data type even when compiled for a 32-bit CPU.

So there are (at least) three varieties of BBC BASIC that you have to account for: ARM BASIC V that DavidS is using (32-bit code, no support for 64-bit integers), BBC BASIC for Windows and the 32-bit editions of BBC BASIC for SDL 2.0 (32-bit code, support for 64-bit integers included) and the 64-bit editions of BBC BASIC for SDL 2.0 (64-bit code, both 32-bit and 64-bit integers).

Sorry if it's confusing, but it's what almost inevitably happens when a language is developed over a period of 38 years!

[DavidS]

RichardRussell,

I agree that Heater should not list a program with lowercase keywords and call it BBC BASIC; at best it's non-standard.

Sounds reasonable. I could have a go at fixing it, but to be sure it get's done properly would you mind doing the honors. Pretty please.

If someone with RichardRussells BBC BASIC would please verify that this version works, then we have it corrected for that version (may need to comment out line 10, I do not know how error handling works in Richard Russells version).

It does not yet work in ARM BASIC (BBC BASIC V on RISC OS), do to getting an integer greater than 32-bits in line 8550. Though I am sure that is not a problem with Richard Russels version.

This is my hand translation of the version that is apparently based on a strange MS-BASIC version:

Code: Select all

       10 : ON ERROR PRINT "LINE : " + STR$(ERL) + " ERROR : " + REPORT$ : END
      11 :     REM fibo_bbc.bas
      20 :     REM This is the classic.bas Fibonacci(4784969) calculator
      30 :    REM with slight modifications to run under Richard Russell's BBCBasic
      40 : 
      50 :  REM Output to file
      60 : *spool fibo/dat
      70 : 
      80 :    REM Initialization required by BBCBasic
      90 :     HIMEM = PAGE + 20000000
      95 :     q1 = 0 : q2 = 0
      98 : 
     100 :  REM classic.bas -- Compute the nth Fibonacci Number
     110 :  REM Written December 25, 2018 by Eric Olson
     120 :  REM
     130 :  REM This program demonstrates the expressiveness of the original
     140 :  REM versions of Microsoft BASIC as measured by explicitly coding
     150 :  REM Karatsuba multiplication for big-number arithmetic and then
     160 :  REM using the doubling formula
     170 :  REM
     180 :  REM     F(2k) = F(k)[2F(k+1)-F(k)]
     190 :  REM   F(2k+1) = F(k+1)[2F(k+1)-F(k)]+(-1)^(k+1)
     200 :  REM   F(2k+1) = F(k)[2F(k)+F(k+1)]+(-1)^(k)
     210 :  REM   F(2k+2) = F(k+1)[2F(k)+F(k+1)]
     220 :  REM
     230 :  REM to compute the nth Fibonacci number.
     231 :  REM
     232 :  REM Version 2:  Minor changes to optimize the O(n^2) multiply and
     233 : : REM prevent overflow when running on MITS BASIC.
     235 : : REM
     236 : : REM To run this program on early versions of Microsoft BASIC please
     237 : : REM remove line 241 to set the default back to single precision
     238 : : REM and change n so the resulting Fibonacci number fits in memory.
     240 :  REM
     241 :  REMdefdbl a-z
     250 :   n=4784969
     251 :  REM n=7499
     252 :  GOSUB 8500
     270 :  d9=INT(n*ln((1+SQR(5))/2)/LN(b8)+7):m9=14
     280 :  DIM m(d9*m9):m0=1
     300 :  a=m0:m0=m0+1+d9:b=m0:m0=m0+1+d9
     310 :  t1=m0:m0=m0+1+d9:t2=m0:m0=m0+1+d9
     400 :  r0=n:GOSUB 1000
     420 :  r7=b:GOSUB 7000
     430 :  STOP
    1000 :  REM Compute nth Fibonacci number
    1010 :  REM    inputs: r0 the value of n
    1020 :  REM   outputs: b the value of F(n)
    1040 :  IF r0<2 THEN m(b)=1:m(b+1)=r0:RETURN
    1060 :  n1=r0:r0=INT((n1-1)/2):GOSUB 1500
    1070 :  p1=n1-4*INT(n1/4)
    1080 :  IF p1=1 OR p1=3 GOTO 1200
    1090 :  r1=t1:r2=a:r3=a:GOSUB 2000
    1110 :  r1=t2:r2=t1:r3=b:GOSUB 2000
    1120 :  r1=t1:r2=b:r3=t2:GOSUB 4000
    1170 :  r1=b:r2=t1:GOSUB 5000:RETURN
    1200 :  r1=t1:r2=b:r3=b:GOSUB 2000
    1210 :  r1=t2:r2=t1:r3=a:GOSUB 3000
    1220 :  r1=t1:r2=b:r3=t2:GOSUB 4000
    1230 :  IF p1=3 THEN 1250
    1240 :  r1=t1:GOSUB 6000:GOTO 1260
    1250 :  r1=t1:GOSUB 6500
    1260 :  r1=b:r2=t1:GOSUB 5000:RETURN
    1500 :  REM Recursive work for nth Fibonacci number
    1510 :  REM    inputs: r0 the value of n
    1520 :  REM   outputs: a the value of F(n)
    1530 :  REM   outputs: b the value of F(n+1)
    1540 :  IF r0=0 THEN m(a)=1:m(a+1)=0:m(b)=1:m(b+1)=1:RETURN
    1600 :  m(m0)=r0:m0=m0+1:r0=INT(r0/2):INT 1500
    1610 :  m0=m0-1:r0=m(m0)
    1620 :  p1=r0-4*INT(r0/4)
    1630 :  IF p1=1 OR p1=3 THEN 1720
    1640 :  r1=t1:r2=b:r3=b:GOSUB 2000
    1650 :  r1=t2:r2=t1:r3=a:GOSUB 3000
    1660 :  r1=t1:r2=a:r3=t2:GOSUB 4000
    1670 :  r1=a:r2=t1:GOSUB 5000
    1680 :  r1=t1:r2=b:r3=t2:GOSUB 4000
    1690 :  IF p1=2 THEN 1710
    1700 :  r1=t1:GOSUB 6000:GOTO 1711
    1710 :  r1=t1:GOSUB 6500
    1711 :  r1=b:r2=t1:GOSUB 5000:RETURN
    1720 :  r1=t1:r2=a:r3=a:GOSUB 2000
    1730 :  r1=t2:r2=t1:r3=b:GOSUB 2000
    1740 :  r1=t1:r2=b:r3=t2:GOSUB 4000
    1750 :  r1=b:r2=t1:GOSUB 5000
    1760 :  r1=t1:r2=a:r3=t2:GOSUB 4000
    1770 :  IF p1=3 THEN 1790
    1780 :  r1=t1:GOSUB 6500:GOTO 1800
    1790 :  r1=t1:GOSUB 6000
    1800 :  r1=a:r2=t1:GOSUB 5000:RETURN
    2000 :  REM Big-number addition
    2010 :  REM  outputs: r1 the value of a+b
    2020 :  REM   inputs: r2 the value of a
    2030 :  REM   inputs: r3 the value of b
    2050 :  IF m(r3)>m(r2) THEN i9=m(r3) ELSE i9=m(r2)
    2060 :  FOR i=1 TO i9+1:m(r1+i)=0:NEXT i
    2070 :  FOR i=1 TO i9:c=0:t=m(r1+i)
    2080 :    IF i<=m(r2) THEN t=t+m(r2+i)
    2090 :    IF i<=m(r3) THEN t=t+m(r3+i)
    2110 :    IF t>=b8 THEN c=1:t=t-b8
    2120 :    m(r1+i)=t:m(r1+i+1)=m(r1+i+1)+c:NEXT i
    2130 :  m(r1)=i9+1
    2140 :  r4=r1:GOSUB 7500
    2150 :  RETURN
    3000 :  REM Big-number subtraction
    3010 :  REM  outputs: r1 the value of a-b
    3020 :  REM   inputs: r2 the value of a
    3030 :  REM   inputs: r3 the value of b
    3050 :  FOR i=1 TO m(r2):m(r1+i)=0:next i
    3060 :  FOR i=1 TO m(r3):t=m(r1+i)+m(r2+i)-m(r3+i)
    3070 :    IF t<0 THEN t=t+b8:m(r1+i+1)=m(r1+i+1)-1
    3080 :    m(r1+i)=t:NEXT i
    3090 :  FOR i=m(r3)+1 TO m(r2):t=m(r1+i)+m(r2+i)
    3100 :    IF t<0 THEN t=t+b8:m(r1+i+1)=m(r1+i+1)-1
    3110 :    m(r1+i)=t:NEXT i
    3120 :  m(r1)=m(r2)
    3130 :  r4=r1:GOSUB 7500
    3150 :  RETURN
    4000 :  REM Big-number multiplication
    4010 :  REM  outputs: r1 the value of a*b
    4020 :  REM   inputs: r2 the value of a
    4030 :  REM   inputs: r3 the value of b
    4040 :  IF m(r2)>80 AND m(r3)>80 THEN 4300
    4050 :  i9=m(r2)+m(r3):FOR i=1 TO i9:m(r1+i)=0:NEXT i
    4070 :  FOR i=1 to m(r2):for j=1 to m(r3)
    4080 :      t=m(r1+i+j-1)+m(r2+i)*m(r3+j)
    4090 :      IF t<b7 THEN 4120
    4100 :      m(r1+i+j-1)=t-b7
    4110 :      m(r1+i+j)=m(r1+i+j)+b6:GOTO 4130
    4120 :      m(r1+i+j-1)=t
    4130 :    NEXT j:NEXT i
    4140 :  c=0:FOR i=1 TO i9:t=m(r1+i)+c
    4150 :    IF t<b8 THEN c=0:GOTO 4170
    4160 :    c=INT(t/b8):t=t-b8*c
    4170 :    m(r1+i)=t:NEXT i
    4180 :  m(r1)=i9
    4190 :  r4=r1:GOSUB 7500
    4230 :  RETURN
    4300 :  REM Big-number Karatsuba algorithm
    4310 :  IF m(r2)<m(r3) THEN i8=m(r3) ELSE i8=m(r2)
    4320 :  i8=int(i8/2)
    4330 :  z0=m0:m0=m0+1+2*i8+1
    4332 :  z2=m0:m0=m0+1+2*i8+3
    4334 :  z1=m0:m0=m0+1+2*i8+5
    4340 :  z3=m0:m0=m0+1+i8+2
    4350 :  z4=m0:m0=m0+1+i8+2
    4360 :  r5=z4:r6=r3:GOSUB 4500
    4370 :  r5=z3:r6=r2:GOSUB 4500
    4380 :  GOSUB 4600:r1=z1:r2=z3:r3=z4:GOSUB 4000:GOSUB 4700
    4400 :  q1=m(r2):IF i8<q1 THEN m(r2)=i8
    4405 :  q2=m(r3):IF i8<q2 THEN m(r3)=i8
    4410 :  GOSUB 4600:r1=z0:GOSUB 4000:GOSUB 4700
    4420 :  m(r2)=q1:m(r3)=q2
    4430 :  q3=q1-i8:q4=q2-i8:IF q3<0 OR q4<0 THEN m(z2)=0:GOTO 8000
    4440 :  q1=m(r2+i8):m(r2+i8)=q3:q2=m(r3+i8):m(r3+i8)=q4
    4450 :  GOSUB 4600:r1=z2:r2=r2+i8:r3=r3+i8:GOSUB 4000:GOSUB 4700
    4460 :  m(r2+i8)=q1:m(r3+i8)=q2
    4470 :  GOTO 8000
    4500 :  REM Add high to low
    4510 :  REM  outputs: r5 the sum of high(a)+low(a)
    4520 :  REM   inputs: r6 the value of a
    4530 :  REM   inputs: i8 the split point
    4540 :  c=0:FOR i=1 TO i8+1:t=c
    4545 :    IF i<=i8 AND i<=m(r6) THEN t=t+m(r6+i)
    4550 :    IF i+i8<=m(r6) THEN t=t+m(r6+i+i8)
    4560 :    IF t>=b8 THEN c=1:t=t-b8 ELSE c=0
    4570 :    m(r5+i)=t:NEXT i:m(r5+i8+2)=c
    4590 :  m(r5)=i8+2:r4=r5:GOSUB 7500
    4595 :  RETURN
    4600 :  REM Save frame
    4610 :  m(m0)=z1:m0=m0+1:m(m0)=z2:m0=m0+1:m(m0)=z0:m0=m0+1
    4620 :  m(m0)=i8:m0=m0+1:m(m0)=q1:m0=m0+1:m(m0)=q2:m0=m0+1
    4630 :  REM Save parameters
    4640 :  m(m0)=r1:m0=m0+1:m(m0)=r2:m0=m0+1:m(m0)=r3:m0=m0+1
    4650 :  RETURN
    4700 :  REM Restore frame
    4710 :  GOSUB 4750
    4720 :  m0=m0-1:q2=m(m0):m0=m0-1:q1=m(m0):m0=m0-1:i8=m(m0)
    4730 :  m0=m0-1:z0=m(m0):m0=m0-1:z2=m(m0):m0=m0-1:z1=m(m0)
    4740 :  RETURN
    4750 :  REM Restore parameters
    4760 :  m0=m0-1:r3=m(m0):m0=m0-1:r2=m(m0):m0=m0-1:r1=m(m0)
    4770 :  RETURN
    5000 :  REM Big-number copy
    5010 :  REM  outputs: r1 the value of a
    5020 :  REM   inputs: r2 the value of a
    5030 :  r4=r2:GOSUB 7500
    5040 :  FOR i=1 TO m(r2):m(r1+i)=m(r2+i):NEXT i
    5050 :  FOR i=m(r2)+1 TO m(r1):m(r1+i)=0:NEXT i
    5060 :  m(r1)=m(r2)
    5070 :  RETURN
    6000 :  REM Big-number decrement
    6010 :  REM   inputs: r1 the value of a
    6020 :  REM  outputs: r1 the value of a-1
    6040 :  i=1:c=1
    6050 :  if c=0 THEN 6080
    6060 :  if m(r1+i)<1 THEN m(r1+i)=b8-1 ELSE m(r1+i)=m(r1+i)-1:c=0
    6070 :  i=i+1:GOTO 6050
    6080 :  r4=r1:GOSUB 7500
    6100 :  RETURN
    6500 :  REM Big-number increment
    6510 :  REM   inputs: r1 the value of a
    6520 :  REM  outputs: r1 the value of a+1
    6540 :  m(r1)=m(r1)+1:m(r1+m(r1))=0:i=1:c=1
    6550 :  if c=0 THEN 6590
    6560 :  t=m(r1+i)+1
    6570 :  if t>=b8 THEN m(r1+i)=t-b8 ELSE m(r1+i)=t:c=0
    6580 :  i=i+1:GOTO 6550
    6590 :  r4=r1:GOSUB 7500
    6600 :  RETURN
    7000 :  REM Big-number print
    7010 :  REM   inputs: r7 the value to print
    7020 :  IF m(r7)=0 THEN PRINT "0":RETURN
    7030 :  FOR i=m(r7) TO 1 STEP -1
    7040 :    s$=STR$(m(r7+i))
    7045 :    IF MID$(s$,1,1)=" " THEN s$=MID$(s$,2,LEN(s$)-1)
    7050 :    IF i=m(r7) OR b9<=LEN(s$) THEN 7070
    7060 :    s$=MID$("0000000000000",1,b9-LEN(s$))+s$
    7070 :    PRINT s$;:NEXT i:PRINT:RETURN
    7500 :  REM Big-number trim
    7510 :  REM    inputs: r4 the value of a
    7520 :  REM   outputs: r4 the trimmed value of a
    7530 :  IF m(r4+m(r4))=0 AND m(r4)>1 THEN m(r4)=m(r4)-1:GOTO 7530
    7540 :  RETURN
    8000 :  REM Tail of Karatsuba
    8003 :  t4=m0-1-2*i8-5
    8005 :  GOSUB 4630:r1=t4:r2=z1:r3=z0:GOSUB 3000
    8010 :  r1=z1:r2=t4:r3=z2:GOSUB 3000:GOSUB 4750
    8020 :  r4=z1:GOSUB 7500
    8030 :  i9=m(r2)+m(r3):i7=2*i8:c=0:FOR i=1 TO i9:t=c
    8040 :    IF i<=m(z0) THEN t=t+m(z0+i)
    8050 :    IF i>i8 AND i-i8<=m(z1) THEN t=t+m(z1+i-i8)
    8060 :    IF i>i7 AND i-i7<=m(z2) THEN t=t+m(z2+i-i7)
    8070 :    IF t<b8 THEN c=0:GOTO 8090
    8080 :    c=INT(t/b8):t=t-b8*c
    8090 :    m(r1+i)=t:NEXT i
    8092 :  m(r1)=i9:r4=r1:GOSUB 7500
    8100 :  m0=z0
    8120 :  RETURN
    8500 :  REM Get dynamic range
    8510 :  REM    outputs: f0 the largest integer
    8511 :  REM    outputs: b7 carry threshold
    8512 :  REM    outputs: b8 radix of limbs
    8513 :  REM    outputs: b9 exponent to b8=10^b9
    8520 :  f9=1
    8530 :  f7=f9+1:IF f7>f9 THEN f9=f9*2:GOTO 8530
    8540 :  f0=f9
    8550 : f7=f0+1:IF f7=f0 THEN f0=INT(f0/2):GOTO 8550
    8560 : f4=0
    8565 : f5=INT(0.5*(f0+f9)+0.5)
    8570 : f7=f5+1:IF f7=f5 THEN f9=f5 ELSE f0=f5
    8580 : IF f5=f4 THEN 8600
    8590 : f4=f5:GOTO 8565
    8600 : f0=f0/2
    8610 : b9=INT(LN(SQR(f0))/LN(10))
    8620 : b8=INT(EXP(LN(10)*b9)+0.5)
    8630 : b7=f0-2*b8^2:IF b7/b8^2<4 THEN b9=b9-1:GOTO 8620
    8631 : b6=int(b7/b8):b7=b8*b6
    8700 : RETURN
    8701 : *spool
    9999 : END

As normal with BBC BASIC simple copy paste (or drag out to text editor) may introduce extra characters from the forum, so take care.

Next up is to attempt to get this version to actually work on BBC BASIC V/ARM BASIC.

The integer size should be a problem for any MS-BASIC crude enough to require such a messed up structure, as it should only allow for 16-bit integers. I think that it was done in a strange way here.

[RichardRussell]

Sounds reasonable. I could have a go at fixing it, but to be sure it get's done properly would you mind doing the honors. Pretty please.

I think this is it:

Code: Select all

      REM fibo_bbc.bas
      REM This is the classic.bas Fibonacci(4784969) calculator
      REM with slight modifications to run under Richard Russell's BBCBasic

      REM Output to file
      *spool fibo.dat

      REM Initialization required by BBCBasic
      HIMEM = PAGE + 20000000
      q1 = 0 : q2 = 0

  100 REM classic.bas -- Compute the nth Fibonacci Number
  110 REM Written December 25, 2018 by Eric Olson
  120 REM
  130 REM This program demonstrates the expressiveness of the original
  140 REM versions of Microsoft BASIC as measured by explicitly coding
  150 REM Karatsuba multiplication for big-number arithmetic and then
  160 REM using the doubling formula
  170 REM
  180 REM     F(2k) = F(k)[2F(k+1)-F(k)]
  190 REM   F(2k+1) = F(k+1)[2F(k+1)-F(k)]+(-1)^(k+1)
  200 REM   F(2k+1) = F(k)[2F(k)+F(k+1)]+(-1)^(k)
  210 REM   F(2k+2) = F(k+1)[2F(k)+F(k+1)]
  220 REM
  230 REM to compute the nth Fibonacci number.
  231 REM
  232 REM Version 2:  Minor changes to optimize the O(n^2) multiply and
  233 REM prevent overflow when running on MITS BASIC.
  235 REM
  236 REM To run this program on early versions of Microsoft BASIC please
  237 REM remove line 241 to set the default back to single precision
  238 REM and change n so the resulting Fibonacci number fits in memory.
  240 REM
  241 DEFdbl a-z
  250 n=4784969
  251 REM n=7499
  252 GOSUB 8500
  270 d9=INT(n*LN((1+SQR(5))/2)/LN(b8)+7):m9=14
  280 DIM m(d9*m9):m0=1
  300 a=m0:m0=m0+1+d9:b=m0:m0=m0+1+d9
  310 t1=m0:m0=m0+1+d9:t2=m0:m0=m0+1+d9
  400 r0=n:GOSUB 1000
  420 r7=b:GOSUB 7000
  430 STOP
 1000 REM Compute nth Fibonacci number
 1010 REM    inputs: r0 the value of n
 1020 REM   outputs: b the value of F(n)
 1040 IF r0<2 THEN m(b)=1:m(b+1)=r0:RETURN
 1060 n1=r0:r0=INT((n1-1)/2):GOSUB 1500
 1070 p1=n1-4*INT(n1/4)
 1080 IF p1=1 OR p1=3 GOTO 1200
 1090 r1=t1:r2=a:r3=a:GOSUB 2000
 1110 r1=t2:r2=t1:r3=b:GOSUB 2000
 1120 r1=t1:r2=b:r3=t2:GOSUB 4000
 1170 r1=b:r2=t1:GOSUB 5000:RETURN
 1200 r1=t1:r2=b:r3=b:GOSUB 2000
 1210 r1=t2:r2=t1:r3=a:GOSUB 3000
 1220 r1=t1:r2=b:r3=t2:GOSUB 4000
 1230 IF p1=3 THEN 1250
 1240 r1=t1:GOSUB 6000:GOTO 1260
 1250 r1=t1:GOSUB 6500
 1260 r1=b:r2=t1:GOSUB 5000:RETURN
 1500 REM Recursive work for nth Fibonacci number
 1510 REM    inputs: r0 the value of n
 1520 REM   outputs: a the value of F(n)
 1530 REM   outputs: b the value of F(n+1)
 1540 IF r0=0 THEN m(a)=1:m(a+1)=0:m(b)=1:m(b+1)=1:RETURN
 1600 m(m0)=r0:m0=m0+1:r0=INT(r0/2):GOSUB 1500
 1610 m0=m0-1:r0=m(m0)
 1620 p1=r0-4*INT(r0/4)
 1630 IF p1=1 OR p1=3 THEN 1720
 1640 r1=t1:r2=b:r3=b:GOSUB 2000
 1650 r1=t2:r2=t1:r3=a:GOSUB 3000
 1660 r1=t1:r2=a:r3=t2:GOSUB 4000
 1670 r1=a:r2=t1:GOSUB 5000
 1680 r1=t1:r2=b:r3=t2:GOSUB 4000
 1690 IF p1=2 THEN 1710
 1700 r1=t1:GOSUB 6000:GOTO 1711
 1710 r1=t1:GOSUB 6500
 1711 r1=b:r2=t1:GOSUB 5000:RETURN
 1720 r1=t1:r2=a:r3=a:GOSUB 2000
 1730 r1=t2:r2=t1:r3=b:GOSUB 2000
 1740 r1=t1:r2=b:r3=t2:GOSUB 4000
 1750 r1=b:r2=t1:GOSUB 5000
 1760 r1=t1:r2=a:r3=t2:GOSUB 4000
 1770 IF p1=3 THEN 1790
 1780 r1=t1:GOSUB 6500:GOTO 1800
 1790 r1=t1:GOSUB 6000
 1800 r1=a:r2=t1:GOSUB 5000:RETURN
 2000 REM Big-number addition
 2010 REM  outputs: r1 the value of a+b
 2020 REM   inputs: r2 the value of a
 2030 REM   inputs: r3 the value of b
 2050 IF m(r3)>m(r2) THEN i9=m(r3) ELSE i9=m(r2)
 2060 FOR i=1 TO i9+1:m(r1+i)=0:NEXT i
 2070 FOR i=1 TO i9:c=0:t=m(r1+i)
 2080   IF i<=m(r2) THEN t=t+m(r2+i)
 2090   IF i<=m(r3) THEN t=t+m(r3+i)
 2110   IF t>=b8 THEN c=1:t=t-b8
 2120   m(r1+i)=t:m(r1+i+1)=m(r1+i+1)+c:NEXT i
 2130 m(r1)=i9+1
 2140 r4=r1:GOSUB 7500
 2150 RETURN
 3000 REM Big-number subtraction
 3010 REM  outputs: r1 the value of a-b
 3020 REM   inputs: r2 the value of a
 3030 REM   inputs: r3 the value of b
 3050 FOR i=1 TO m(r2):m(r1+i)=0:NEXT i
 3060 FOR i=1 TO m(r3):t=m(r1+i)+m(r2+i)-m(r3+i)
 3070   IF t<0 THEN t=t+b8:m(r1+i+1)=m(r1+i+1)-1
 3080   m(r1+i)=t:NEXT i
 3090 FOR i=m(r3)+1 TO m(r2):t=m(r1+i)+m(r2+i)
 3100   IF t<0 THEN t=t+b8:m(r1+i+1)=m(r1+i+1)-1
 3110   m(r1+i)=t:NEXT i
 3120 m(r1)=m(r2)
 3130 r4=r1:GOSUB 7500
 3150 RETURN
 4000 REM Big-number multiplication
 4010 REM  outputs: r1 the value of a*b
 4020 REM   inputs: r2 the value of a
 4030 REM   inputs: r3 the value of b
 4040 IF m(r2)>80 AND m(r3)>80 THEN 4300
 4050 i9=m(r2)+m(r3):FOR i=1 TO i9:m(r1+i)=0:NEXT i
 4070 FOR i=1 TO m(r2):FOR j=1 TO m(r3)
 4080     t=m(r1+i+j-1)+m(r2+i)*m(r3+j)
 4090     IF t<b7 THEN 4120
 4100     m(r1+i+j-1)=t-b7
 4110     m(r1+i+j)=m(r1+i+j)+b6:GOTO 4130
 4120     m(r1+i+j-1)=t
 4130   NEXT j:NEXT i
 4140 c=0:FOR i=1 TO i9:t=m(r1+i)+c
 4150   IF t<b8 THEN c=0:GOTO 4170
 4160   c=INT(t/b8):t=t-b8*c
 4170   m(r1+i)=t:NEXT i
 4180 m(r1)=i9
 4190 r4=r1:GOSUB 7500
 4230 RETURN
 4300 REM Big-number Karatsuba algorithm
 4310 IF m(r2)<m(r3) THEN i8=m(r3) ELSE i8=m(r2)
 4320 i8=INT(i8/2)
 4330 z0=m0:m0=m0+1+2*i8+1
 4332 z2=m0:m0=m0+1+2*i8+3
 4334 z1=m0:m0=m0+1+2*i8+5
 4340 z3=m0:m0=m0+1+i8+2
 4350 z4=m0:m0=m0+1+i8+2
 4360 r5=z4:r6=r3:GOSUB 4500
 4370 r5=z3:r6=r2:GOSUB 4500
 4380 GOSUB 4600:r1=z1:r2=z3:r3=z4:GOSUB 4000:GOSUB 4700
 4400 q1=m(r2):IF i8<q1 THEN m(r2)=i8
 4405 q2=m(r3):IF i8<q2 THEN m(r3)=i8
 4410 GOSUB 4600:r1=z0:GOSUB 4000:GOSUB 4700
 4420 m(r2)=q1:m(r3)=q2
 4430 q3=q1-i8:q4=q2-i8:IF q3<0 OR q4<0 THEN m(z2)=0:GOTO 8000
 4440 q1=m(r2+i8):m(r2+i8)=q3:q2=m(r3+i8):m(r3+i8)=q4
 4450 GOSUB 4600:r1=z2:r2=r2+i8:r3=r3+i8:GOSUB 4000:GOSUB 4700
 4460 m(r2+i8)=q1:m(r3+i8)=q2
 4470 GOTO 8000
 4500 REM Add high to low
 4510 REM  outputs: r5 the sum of high(a)+low(a)
 4520 REM   inputs: r6 the value of a
 4530 REM   inputs: i8 the split point
 4540 c=0:FOR i=1 TO i8+1:t=c
 4545   IF i<=i8 AND i<=m(r6) THEN t=t+m(r6+i)
 4550   IF i+i8<=m(r6) THEN t=t+m(r6+i+i8)
 4560   IF t>=b8 THEN c=1:t=t-b8 ELSE c=0
 4570   m(r5+i)=t:NEXT i:m(r5+i8+2)=c
 4590 m(r5)=i8+2:r4=r5:GOSUB 7500
 4595 RETURN
 4600 REM Save frame
 4610 m(m0)=z1:m0=m0+1:m(m0)=z2:m0=m0+1:m(m0)=z0:m0=m0+1
 4620 m(m0)=i8:m0=m0+1:m(m0)=q1:m0=m0+1:m(m0)=q2:m0=m0+1
 4630 REM Save parameters
 4640 m(m0)=r1:m0=m0+1:m(m0)=r2:m0=m0+1:m(m0)=r3:m0=m0+1
 4650 RETURN
 4700 REM Restore frame
 4710 GOSUB 4750
 4720 m0=m0-1:q2=m(m0):m0=m0-1:q1=m(m0):m0=m0-1:i8=m(m0)
 4730 m0=m0-1:z0=m(m0):m0=m0-1:z2=m(m0):m0=m0-1:z1=m(m0)
 4740 RETURN
 4750 REM Restore parameters
 4760 m0=m0-1:r3=m(m0):m0=m0-1:r2=m(m0):m0=m0-1:r1=m(m0)
 4770 RETURN
 5000 REM Big-number copy
 5010 REM  outputs: r1 the value of a
 5020 REM   inputs: r2 the value of a
 5030 r4=r2:GOSUB 7500
 5040 FOR i=1 TO m(r2):m(r1+i)=m(r2+i):NEXT i
 5050 FOR i=m(r2)+1 TO m(r1):m(r1+i)=0:NEXT i
 5060 m(r1)=m(r2)
 5070 RETURN
 6000 REM Big-number decrement
 6010 REM   inputs: r1 the value of a
 6020 REM  outputs: r1 the value of a-1
 6040 i=1:c=1
 6050 IF c=0 THEN 6080
 6060 IF m(r1+i)<1 THEN m(r1+i)=b8-1 ELSE m(r1+i)=m(r1+i)-1:c=0
 6070 i=i+1:GOTO 6050
 6080 r4=r1:GOSUB 7500
 6100 RETURN
 6500 REM Big-number increment
 6510 REM   inputs: r1 the value of a
 6520 REM  outputs: r1 the value of a+1
 6540 m(r1)=m(r1)+1:m(r1+m(r1))=0:i=1:c=1
 6550 IF c=0 THEN 6590
 6560 t=m(r1+i)+1
 6570 IF t>=b8 THEN m(r1+i)=t-b8 ELSE m(r1+i)=t:c=0
 6580 i=i+1:GOTO 6550
 6590 r4=r1:GOSUB 7500
 6600 RETURN
 7000 REM Big-number print
 7010 REM   inputs: r7 the value to print
 7020 IF m(r7)=0 THEN PRINT "0":RETURN
 7030 FOR i=m(r7) TO 1 STEP -1
 7040   s$=STR$(m(r7+i))
 7045   IF MID$(s$,1,1)=" " THEN s$=MID$(s$,2,LEN(s$)-1)
 7050   IF i=m(r7) OR b9<=LEN(s$) THEN 7070
 7060   s$=MID$("0000000000000",1,b9-LEN(s$))+s$
 7070   PRINT s$;:NEXT i:PRINT:RETURN
 7500 REM Big-number trim
 7510 REM    inputs: r4 the value of a
 7520 REM   outputs: r4 the trimmed value of a
 7530 IF m(r4+m(r4))=0 AND m(r4)>1 THEN m(r4)=m(r4)-1:GOTO 7530
 7540 RETURN
 8000 REM Tail of Karatsuba
 8003 t4=m0-1-2*i8-5
 8005 GOSUB 4630:r1=t4:r2=z1:r3=z0:GOSUB 3000
 8010 r1=z1:r2=t4:r3=z2:GOSUB 3000:GOSUB 4750
 8020 r4=z1:GOSUB 7500
 8030 i9=m(r2)+m(r3):i7=2*i8:c=0:FOR i=1 TO i9:t=c
 8040   IF i<=m(z0) THEN t=t+m(z0+i)
 8050   IF i>i8 AND i-i8<=m(z1) THEN t=t+m(z1+i-i8)
 8060   IF i>i7 AND i-i7<=m(z2) THEN t=t+m(z2+i-i7)
 8070   IF t<b8 THEN c=0:GOTO 8090
 8080   c=INT(t/b8):t=t-b8*c
 8090   m(r1+i)=t:NEXT i
 8092 m(r1)=i9:r4=r1:GOSUB 7500
 8100 m0=z0
 8120 RETURN
 8500 REM Get dynamic range
 8510 REM    outputs: f0 the largest integer
 8511 REM    outputs: b7 carry threshold
 8512 REM    outputs: b8 radix of limbs
 8513 REM    outputs: b9 exponent to b8=10^b9
 8520 f9=1
 8530 f7=f9+1:IF f7>f9 THEN f9=f9*2:GOTO 8530
 8540 f0=f9
 8550 f7=f0+1:IF f7=f0 THEN f0=INT(f0/2):GOTO 8550
 8560 f4=0
 8565 f5=INT(0.5*(f0+f9)+0.5)
 8570 f7=f5+1:IF f7=f5 THEN f9=f5 ELSE f0=f5
 8580 IF f5=f4 THEN 8600
 8590 f4=f5:GOTO 8565
 8600 f0=f0/2
 8610 b9=INT(LN(SQR(f0))/LN(10))
 8620 b8=INT(EXP(LN(10)*b9)+0.5)
 8630 b7=f0-2*b8^2:IF b7/b8^2<4 THEN b9=b9-1:GOTO 8620
 8640 b6=INT(b7/b8):b7=b8*b6
 8700 RETURN
 9999 END
 

[DavidS]

Sounds reasonable. I could have a go at fixing it, but to be sure it get's done properly would you mind doing the honors. Pretty please.

I think this is it:

Code: Select all

      REM fibo_bbc.bas
      REM This is the classic.bas Fibonacci(4784969) calculator
      REM with slight modifications to run under Richard Russell's BBCBasic

      REM Output to file
      *spool fibo.dat

      REM Initialization required by BBCBasic
      HIMEM = PAGE + 20000000
      q1 = 0 : q2 = 0

  100 REM classic.bas -- Compute the nth Fibonacci Number
  110 REM Written December 25, 2018 by Eric Olson
  120 REM
  130 REM This program demonstrates the expressiveness of the original
  140 REM versions of Microsoft BASIC as measured by explicitly coding
  150 REM Karatsuba multiplication for big-number arithmetic and then
  160 REM using the doubling formula
  170 REM
  180 REM     F(2k) = F(k)[2F(k+1)-F(k)]
  190 REM   F(2k+1) = F(k+1)[2F(k+1)-F(k)]+(-1)^(k+1)
  200 REM   F(2k+1) = F(k)[2F(k)+F(k+1)]+(-1)^(k)
  210 REM   F(2k+2) = F(k+1)[2F(k)+F(k+1)]
  220 REM
  230 REM to compute the nth Fibonacci number.
  231 REM
  232 REM Version 2:  Minor changes to optimize the O(n^2) multiply and
  233 REM prevent overflow when running on MITS BASIC.
  235 REM
  236 REM To run this program on early versions of Microsoft BASIC please
  237 REM remove line 241 to set the default back to single precision
  238 REM and change n so the resulting Fibonacci number fits in memory.
  240 REM
  241 DEFdbl a-z
  250 n=4784969
  251 REM n=7499
  252 GOSUB 8500
  270 d9=INT(n*LN((1+SQR(5))/2)/LN(b8)+7):m9=14
  280 DIM m(d9*m9):m0=1
  300 a=m0:m0=m0+1+d9:b=m0:m0=m0+1+d9
  310 t1=m0:m0=m0+1+d9:t2=m0:m0=m0+1+d9
  400 r0=n:GOSUB 1000
  420 r7=b:GOSUB 7000
  430 STOP
 1000 REM Compute nth Fibonacci number
 1010 REM    inputs: r0 the value of n
 1020 REM   outputs: b the value of F(n)
 1040 IF r0<2 THEN m(b)=1:m(b+1)=r0:RETURN
 1060 n1=r0:r0=INT((n1-1)/2):GOSUB 1500
 1070 p1=n1-4*INT(n1/4)
 1080 IF p1=1 OR p1=3 GOTO 1200
 1090 r1=t1:r2=a:r3=a:GOSUB 2000
 1110 r1=t2:r2=t1:r3=b:GOSUB 2000
 1120 r1=t1:r2=b:r3=t2:GOSUB 4000
 1170 r1=b:r2=t1:GOSUB 5000:RETURN
 1200 r1=t1:r2=b:r3=b:GOSUB 2000
 1210 r1=t2:r2=t1:r3=a:GOSUB 3000
 1220 r1=t1:r2=b:r3=t2:GOSUB 4000
 1230 IF p1=3 THEN 1250
 1240 r1=t1:GOSUB 6000:GOTO 1260
 1250 r1=t1:GOSUB 6500
 1260 r1=b:r2=t1:GOSUB 5000:RETURN
 1500 REM Recursive work for nth Fibonacci number
 1510 REM    inputs: r0 the value of n
 1520 REM   outputs: a the value of F(n)
 1530 REM   outputs: b the value of F(n+1)
 1540 IF r0=0 THEN m(a)=1:m(a+1)=0:m(b)=1:m(b+1)=1:RETURN
 1600 m(m0)=r0:m0=m0+1:r0=INT(r0/2):GOSUB 1500
 1610 m0=m0-1:r0=m(m0)
 1620 p1=r0-4*INT(r0/4)
 1630 IF p1=1 OR p1=3 THEN 1720
 1640 r1=t1:r2=b:r3=b:GOSUB 2000
 1650 r1=t2:r2=t1:r3=a:GOSUB 3000
 1660 r1=t1:r2=a:r3=t2:GOSUB 4000
 1670 r1=a:r2=t1:GOSUB 5000
 1680 r1=t1:r2=b:r3=t2:GOSUB 4000
 1690 IF p1=2 THEN 1710
 1700 r1=t1:GOSUB 6000:GOTO 1711
 1710 r1=t1:GOSUB 6500
 1711 r1=b:r2=t1:GOSUB 5000:RETURN
 1720 r1=t1:r2=a:r3=a:GOSUB 2000
 1730 r1=t2:r2=t1:r3=b:GOSUB 2000
 1740 r1=t1:r2=b:r3=t2:GOSUB 4000
 1750 r1=b:r2=t1:GOSUB 5000
 1760 r1=t1:r2=a:r3=t2:GOSUB 4000
 1770 IF p1=3 THEN 1790
 1780 r1=t1:GOSUB 6500:GOTO 1800
 1790 r1=t1:GOSUB 6000
 1800 r1=a:r2=t1:GOSUB 5000:RETURN
 2000 REM Big-number addition
 2010 REM  outputs: r1 the value of a+b
 2020 REM   inputs: r2 the value of a
 2030 REM   inputs: r3 the value of b
 2050 IF m(r3)>m(r2) THEN i9=m(r3) ELSE i9=m(r2)
 2060 FOR i=1 TO i9+1:m(r1+i)=0:NEXT i
 2070 FOR i=1 TO i9:c=0:t=m(r1+i)
 2080   IF i<=m(r2) THEN t=t+m(r2+i)
 2090   IF i<=m(r3) THEN t=t+m(r3+i)
 2110   IF t>=b8 THEN c=1:t=t-b8
 2120   m(r1+i)=t:m(r1+i+1)=m(r1+i+1)+c:NEXT i
 2130 m(r1)=i9+1
 2140 r4=r1:GOSUB 7500
 2150 RETURN
 3000 REM Big-number subtraction
 3010 REM  outputs: r1 the value of a-b
 3020 REM   inputs: r2 the value of a
 3030 REM   inputs: r3 the value of b
 3050 FOR i=1 TO m(r2):m(r1+i)=0:NEXT i
 3060 FOR i=1 TO m(r3):t=m(r1+i)+m(r2+i)-m(r3+i)
 3070   IF t<0 THEN t=t+b8:m(r1+i+1)=m(r1+i+1)-1
 3080   m(r1+i)=t:NEXT i
 3090 FOR i=m(r3)+1 TO m(r2):t=m(r1+i)+m(r2+i)
 3100   IF t<0 THEN t=t+b8:m(r1+i+1)=m(r1+i+1)-1
 3110   m(r1+i)=t:NEXT i
 3120 m(r1)=m(r2)
 3130 r4=r1:GOSUB 7500
 3150 RETURN
 4000 REM Big-number multiplication
 4010 REM  outputs: r1 the value of a*b
 4020 REM   inputs: r2 the value of a
 4030 REM   inputs: r3 the value of b
 4040 IF m(r2)>80 AND m(r3)>80 THEN 4300
 4050 i9=m(r2)+m(r3):FOR i=1 TO i9:m(r1+i)=0:NEXT i
 4070 FOR i=1 TO m(r2):FOR j=1 TO m(r3)
 4080     t=m(r1+i+j-1)+m(r2+i)*m(r3+j)
 4090     IF t<b7 THEN 4120
 4100     m(r1+i+j-1)=t-b7
 4110     m(r1+i+j)=m(r1+i+j)+b6:GOTO 4130
 4120     m(r1+i+j-1)=t
 4130   NEXT j:NEXT i
 4140 c=0:FOR i=1 TO i9:t=m(r1+i)+c
 4150   IF t<b8 THEN c=0:GOTO 4170
 4160   c=INT(t/b8):t=t-b8*c
 4170   m(r1+i)=t:NEXT i
 4180 m(r1)=i9
 4190 r4=r1:GOSUB 7500
 4230 RETURN
 4300 REM Big-number Karatsuba algorithm
 4310 IF m(r2)<m(r3) THEN i8=m(r3) ELSE i8=m(r2)
 4320 i8=INT(i8/2)
 4330 z0=m0:m0=m0+1+2*i8+1
 4332 z2=m0:m0=m0+1+2*i8+3
 4334 z1=m0:m0=m0+1+2*i8+5
 4340 z3=m0:m0=m0+1+i8+2
 4350 z4=m0:m0=m0+1+i8+2
 4360 r5=z4:r6=r3:GOSUB 4500
 4370 r5=z3:r6=r2:GOSUB 4500
 4380 GOSUB 4600:r1=z1:r2=z3:r3=z4:GOSUB 4000:GOSUB 4700
 4400 q1=m(r2):IF i8<q1 THEN m(r2)=i8
 4405 q2=m(r3):IF i8<q2 THEN m(r3)=i8
 4410 GOSUB 4600:r1=z0:GOSUB 4000:GOSUB 4700
 4420 m(r2)=q1:m(r3)=q2
 4430 q3=q1-i8:q4=q2-i8:IF q3<0 OR q4<0 THEN m(z2)=0:GOTO 8000
 4440 q1=m(r2+i8):m(r2+i8)=q3:q2=m(r3+i8):m(r3+i8)=q4
 4450 GOSUB 4600:r1=z2:r2=r2+i8:r3=r3+i8:GOSUB 4000:GOSUB 4700
 4460 m(r2+i8)=q1:m(r3+i8)=q2
 4470 GOTO 8000
 4500 REM Add high to low
 4510 REM  outputs: r5 the sum of high(a)+low(a)
 4520 REM   inputs: r6 the value of a
 4530 REM   inputs: i8 the split point
 4540 c=0:FOR i=1 TO i8+1:t=c
 4545   IF i<=i8 AND i<=m(r6) THEN t=t+m(r6+i)
 4550   IF i+i8<=m(r6) THEN t=t+m(r6+i+i8)
 4560   IF t>=b8 THEN c=1:t=t-b8 ELSE c=0
 4570   m(r5+i)=t:NEXT i:m(r5+i8+2)=c
 4590 m(r5)=i8+2:r4=r5:GOSUB 7500
 4595 RETURN
 4600 REM Save frame
 4610 m(m0)=z1:m0=m0+1:m(m0)=z2:m0=m0+1:m(m0)=z0:m0=m0+1
 4620 m(m0)=i8:m0=m0+1:m(m0)=q1:m0=m0+1:m(m0)=q2:m0=m0+1
 4630 REM Save parameters
 4640 m(m0)=r1:m0=m0+1:m(m0)=r2:m0=m0+1:m(m0)=r3:m0=m0+1
 4650 RETURN
 4700 REM Restore frame
 4710 GOSUB 4750
 4720 m0=m0-1:q2=m(m0):m0=m0-1:q1=m(m0):m0=m0-1:i8=m(m0)
 4730 m0=m0-1:z0=m(m0):m0=m0-1:z2=m(m0):m0=m0-1:z1=m(m0)
 4740 RETURN
 4750 REM Restore parameters
 4760 m0=m0-1:r3=m(m0):m0=m0-1:r2=m(m0):m0=m0-1:r1=m(m0)
 4770 RETURN
 5000 REM Big-number copy
 5010 REM  outputs: r1 the value of a
 5020 REM   inputs: r2 the value of a
 5030 r4=r2:GOSUB 7500
 5040 FOR i=1 TO m(r2):m(r1+i)=m(r2+i):NEXT i
 5050 FOR i=m(r2)+1 TO m(r1):m(r1+i)=0:NEXT i
 5060 m(r1)=m(r2)
 5070 RETURN
 6000 REM Big-number decrement
 6010 REM   inputs: r1 the value of a
 6020 REM  outputs: r1 the value of a-1
 6040 i=1:c=1
 6050 IF c=0 THEN 6080
 6060 IF m(r1+i)<1 THEN m(r1+i)=b8-1 ELSE m(r1+i)=m(r1+i)-1:c=0
 6070 i=i+1:GOTO 6050
 6080 r4=r1:GOSUB 7500
 6100 RETURN
 6500 REM Big-number increment
 6510 REM   inputs: r1 the value of a
 6520 REM  outputs: r1 the value of a+1
 6540 m(r1)=m(r1)+1:m(r1+m(r1))=0:i=1:c=1
 6550 IF c=0 THEN 6590
 6560 t=m(r1+i)+1
 6570 IF t>=b8 THEN m(r1+i)=t-b8 ELSE m(r1+i)=t:c=0
 6580 i=i+1:GOTO 6550
 6590 r4=r1:GOSUB 7500
 6600 RETURN
 7000 REM Big-number print
 7010 REM   inputs: r7 the value to print
 7020 IF m(r7)=0 THEN PRINT "0":RETURN
 7030 FOR i=m(r7) TO 1 STEP -1
 7040   s$=STR$(m(r7+i))
 7045   IF MID$(s$,1,1)=" " THEN s$=MID$(s$,2,LEN(s$)-1)
 7050   IF i=m(r7) OR b9<=LEN(s$) THEN 7070
 7060   s$=MID$("0000000000000",1,b9-LEN(s$))+s$
 7070   PRINT s$;:NEXT i:PRINT:RETURN
 7500 REM Big-number trim
 7510 REM    inputs: r4 the value of a
 7520 REM   outputs: r4 the trimmed value of a
 7530 IF m(r4+m(r4))=0 AND m(r4)>1 THEN m(r4)=m(r4)-1:GOTO 7530
 7540 RETURN
 8000 REM Tail of Karatsuba
 8003 t4=m0-1-2*i8-5
 8005 GOSUB 4630:r1=t4:r2=z1:r3=z0:GOSUB 3000
 8010 r1=z1:r2=t4:r3=z2:GOSUB 3000:GOSUB 4750
 8020 r4=z1:GOSUB 7500
 8030 i9=m(r2)+m(r3):i7=2*i8:c=0:FOR i=1 TO i9:t=c
 8040   IF i<=m(z0) THEN t=t+m(z0+i)
 8050   IF i>i8 AND i-i8<=m(z1) THEN t=t+m(z1+i-i8)
 8060   IF i>i7 AND i-i7<=m(z2) THEN t=t+m(z2+i-i7)
 8070   IF t<b8 THEN c=0:GOTO 8090
 8080   c=INT(t/b8):t=t-b8*c
 8090   m(r1+i)=t:NEXT i
 8092 m(r1)=i9:r4=r1:GOSUB 7500
 8100 m0=z0
 8120 RETURN
 8500 REM Get dynamic range
 8510 REM    outputs: f0 the largest integer
 8511 REM    outputs: b7 carry threshold
 8512 REM    outputs: b8 radix of limbs
 8513 REM    outputs: b9 exponent to b8=10^b9
 8520 f9=1
 8530 f7=f9+1:IF f7>f9 THEN f9=f9*2:GOTO 8530
 8540 f0=f9
 8550 f7=f0+1:IF f7=f0 THEN f0=INT(f0/2):GOTO 8550
 8560 f4=0
 8565 f5=INT(0.5*(f0+f9)+0.5)
 8570 f7=f5+1:IF f7=f5 THEN f9=f5 ELSE f0=f5
 8580 IF f5=f4 THEN 8600
 8590 f4=f5:GOTO 8565
 8600 f0=f0/2
 8610 b9=INT(LN(SQR(f0))/LN(10))
 8620 b8=INT(EXP(LN(10)*b9)+0.5)
 8630 b7=f0-2*b8^2:IF b7/b8^2<4 THEN b9=b9-1:GOTO 8620
 8640 b6=INT(b7/b8):b7=b8*b6
 8700 RETURN
 9999 END
 

And how is my upercased RISC OS formated text version on your BBC BASIC?

[RichardRussell]

It does not yet work in ARM BASIC (BBC BASIC V on RISC OS), do to getting an integer greater than 32-bits in line 8550. Though I am sure that is not a problem with Richard Russels version.

I believe it requires at least double (64-bit) precision floats, so if that's the only issue try running it in ARM BASIC 6 rather than 5. I'm slightly surprised that you're not using BASIC 6 routinely, is that simply because of the inconvenience of soft-loading it in RISC OS?

[DavidS]

It does not yet work in ARM BASIC (BBC BASIC V on RISC OS), do to getting an integer greater than 32-bits in line 8550. Though I am sure that is not a problem with Richard Russels version.

I believe it requires at least double (64-bit) precision floats, so if that's the only issue try running it in ARM BASIC 6 rather than 5. I'm slightly surprised that you're not using BASIC 6 routinely, is that simply because of the inconvenience of soft-loading it in RISC OS?

Actually I am using BASIC 6 (I just think of it as BASIC V, as it is the same langage). By default running a program in the desktop uses BASIC 6, and more specifically BASICVFP in the case of the current versions (uses VFP for HW floating point), as BASIC 6 is just an alias of BASICVFP.

No need to softload, it is in the ROM as a ROMModule, and has been for quite a long time.

Code: Select all

*ROMModules
No. Position    Module Name             Version Status
  1 System ROM  UtilityModule           5.24    Active
  2 System ROM  PCI                     0.15    Active
  3 System ROM  FileSwitch              2.86    Active
  4 System ROM  ResourceFS              0.26    Active
  5 System ROM  TerritoryManager        0.57    Active
  6 System ROM  Messages                1.12    Active
  7 System ROM  MessageTrans            0.49    Active
  8 System ROM  UK                      0.61    Active
  9 System ROM  WindowManager           5.61    Active
 10 System ROM  TaskManager             1.50    Active
 11 System ROM  Desktop                 2.76    Active
 12 System ROM  SharedCLibrary          5.95    Active
 13 System ROM  BASIC                   1.75    Active
 14 System ROM  BASIC64                 1.75    Active
 15 System ROM  BASICVFP                1.75    Active
 16 System ROM  BlendTable              0.02    Active
 17 System ROM  BufferManager           0.39    Active
 18 System ROM  ColourTrans             1.95    Active
 19 System ROM  Debugger                2.02    Active
 20 System ROM  DeviceFS                0.73    Active
 21 System ROM  BCMSupport              0.07    Active
 22 System ROM  Portable                0.81    Active
 23 System ROM  RTSupport               0.17    Active
 24 System ROM  USBDriver               1.28    Active
 25 System ROM  DWCDriver               0.34    Active
 26 System ROM  VCHIQ                   0.10    Active
 27 System ROM  BCMSound                0.04    Active
 28 System ROM  ScreenModes             0.68    Active
 29 System ROM  BCMVideo                0.50    Active
 30 System ROM  DisplayManager          0.44    Active
 31 System ROM  DMAManager              0.23    Active
 32 System ROM  DragASprite             0.21    Active
 33 System ROM  DragAnObject            0.09    Active
 34 System ROM  Draw                    1.22    Active
 35 System ROM  BBCEconet               0.20    Dormant
 36 System ROM  FileCore                3.75    Active
 37 System ROM  RamFS                   2.32    Active
 38 System ROM  Filer                   2.43    Active
 39 System ROM  FilerSWIs               0.05    Active
 40 System ROM  FSLock                  1.24    Active
 41 System ROM  FontManager             3.79    Active
 42 System ROM  FPEmulator              4.36    Active
 43 System ROM  VFPSupport              0.13    Active
 44 System ROM  Free                    0.42    Active
 45 System ROM  Hourglass               2.19    Active
 46 System ROM  IIC                     0.18    Active
 47 System ROM  International           1.69    Active
 48 System ROM  InternationalKeyboard   0.99    Active
 49 System ROM  InverseTable            0.18    Active
 50 System ROM  NetFS                   6.00    Dormant
 51 System ROM  NetFiler                0.86    Active
 52 System ROM  NetPrint                5.61    Dormant
 53 System ROM  NetStatus               2.09    Active
 54 System ROM  Obey                    0.40    Active
 55 System ROM  Pinboard                1.03    Active
 56 System ROM  PipeFS                  0.23    Active
 57 System ROM  RAMFSFiler              0.40    Active
 58 System ROM  ResourceFiler           0.20    Active
 59 System ROM  ROMFonts                0.77    Active
 60 System ROM  RTC                     0.01    Dormant
 61 System ROM  ScreenBlanker           2.34    Active
 62 System ROM  ScrSaver                0.14    Active
 63 System ROM  Serial                  0.25    Active
 64 System ROM  SerialDeviceSupport     0.26    Active
 65 System ROM  ShellCLI                0.37    Active
 66 System ROM  SoundDMA                1.85    Active
 67 System ROM  SoundControl            1.06    Active
 68 System ROM  SoundChannels           1.41    Active
 69 System ROM  SoundScheduler          1.32    Active
 70 System ROM  SpriteExtend            1.82    Active
 71 System ROM  SpriteUtils             1.13    Active
 72 System ROM  Squash                  0.30    Active
 73 System ROM  BootFX                  0.08    Active
 74 System ROM  SuperSample             0.16    Active
 75 System ROM  SystemDevices           1.33    Active
 76 System ROM  TaskWindow              0.80    Running
 77 System ROM  WindowUtils             2.53    Active
 78 System ROM  FilterManager           0.28    Active
 79 System ROM  WaveSynth               1.25    Active
 80 System ROM  StringLib               1.19    Active
 81 System ROM  Percussion              1.18    Active
 82 System ROM  SharedSound             1.20    Active
 83 System ROM  Filer_Action            0.62    Active
 84 System ROM  DOSFS                   1.13    Active
 85 System ROM  SCSIDriver              2.14    Active
 86 System ROM  SCSISoftUSB             0.21    Active
 87 System ROM  SCSIFS                  1.32    Active
 88 System ROM  SCSIFiler               1.30    Active
 89 System ROM  SDIODriver              0.13    Active
 90 System ROM  SDFS                    0.10    Active
 91 System ROM  SDFSFiler               1.12    Active
 92 System ROM  SDCMOS                  0.18    Active
 93 System ROM  ColourPicker            0.56    Active
 94 System ROM  DrawFile                1.59    Active
 95 System ROM  BootCommands            1.49    Active
 96 System ROM  AUNMsgs                 0.30    Active
 97 System ROM  MbufManager             0.30    Active
 98 System ROM  Internet                5.63    Active
 99 System ROM  Resolver                0.72    Active
100 System ROM  Net                     6.26    Dormant
101 System ROM  BootNet                 0.97    Active
102 System ROM  Freeway                 0.41    Unplugged
103 System ROM  ShareFS                 3.59    Unplugged
104 System ROM  MimeMap                 0.19    Active
105 System ROM  LanManFS                2.61    Active
106 System ROM  EtherUSB                0.41    Active
107 System ROM  DHCP                    0.27    Active
108 System ROM  AcornURI                0.26    Active
109 System ROM  !Edit                   1.73    Active
110 System ROM  !Draw                   1.30    Active
111 System ROM  !Paint                  2.23    Active
112 System ROM  !Alarm                  2.93    Active
113 System ROM  !Chars                  2.04    Active
114 System ROM  !Help                   3.25    Active
115 System ROM  Toolbox                 1.58    Active
116 System ROM  Window                  1.79    Active
117 System ROM  ToolAction              0.38    Active
118 System ROM  Menu                    0.40    Active
119 System ROM  Iconbar                 1.23    Active
120 System ROM  ColourDbox              0.22    Active
121 System ROM  ColourMenu              0.22    Active
122 System ROM  DCS                     1.14    Active
123 System ROM  FileInfo                0.20    Active
124 System ROM  FontDbox                0.19    Active
125 System ROM  FontMenu                0.25    Active
126 System ROM  PrintDbox               0.18    Active
127 System ROM  ProgInfo                0.19    Active
128 System ROM  SaveAs                  0.20    Active
129 System ROM  Scale                   0.16    Active
130 System ROM  TextGadgets             0.44    Active
131 System ROM  CDFSDriver              2.41    Active
132 System ROM  CDFSSoftSCSI            0.06    Active
133 System ROM  CDFS                    2.63    Active
134 System ROM  CDFSFiler               2.51    Active
135 System ROM  UnSqueezeAIF            0.10    Active
136 System ROM  GPIO                    1.00    Active
*ShellCLI_TaskQuit

Note ROM Modules number 14 and 15 in that list. BASIC64 (BASIC 6) just loads BASICVFP.

[RichardRussell]

And how is my upercased RISC OS formated text version on your BBC BASIC?

I've not tried it, because it's got colons at the beginning of every line. Where have they come from? Typically they will break BBC BASIC:

Code: Select all

10 : IF FALSE THEN
20 : ENDIF

This reports 'Missing ENDIF in line 10'.

[RichardRussell]

Actually I am using BASIC 6 (I just think of it as BASIC V, as it is the same langage).

Right. In that case I don't think the 'classic BASIC' Fibo program should fail, because 64-bit floats can hold integers up to about 53 bits. If it is failing I don't understand why. What's the error message?

[DavidS]

Actually I am using BASIC 6 (I just think of it as BASIC V, as it is the same langage).

Right. In that case I don't think the 'classic BASIC' Fibo program should fail, because 64-bit floats can hold integers up to about 53 bits. If it is failing I don't understand why. What's the error message?

Code: Select all

*BASICVFP -quit basdFIBO
LINE : 8550 ERROR : Invalid arithmetic operation
*BASIC64 -quit basdFIBO
LINE : 8550 ERROR : Invalid arithmetic operation
*BASIC -quit basdFIBO
LINE : 8550 ERROR : Number too big
*

Also note that the use of the INT function requires that the result be stored as an integer at least for a moment in an ARM register.

[ejolson]

Actually I am using BASIC 6 (I just think of it as BASIC V, as it is the same langage).

Right. In that case I don't think the 'classic BASIC' Fibo program should fail, because 64-bit floats can hold integers up to about 53 bits. If it is failing I don't understand why. What's the error message?

The classic BASIC program should automatically detect the floating-point precision and proceed accordingly. It has been tested on the original MITS/Altair extended Basic written by Microsoft (with a smaller size n and 32-bit floats) as well as FreeBasic, bwBasic, Matrix Brandy Basic, Richard Russell's BBC Basic, Lennart Benschop's gplBasic and ScriptBasic. Minor changes were needed in some cases.

Probably the easiest starting point for RISC OS would be the Matrix Brandy Basic version which I'll include in a forthcoming post.

[RichardRussell]

Also note that the use of the INT function requires that the result be stored as an integer at least for a moment in an ARM register.

If that means you can't use INT() to truncate a value greater than 2^31 to an integer (despite BASIC 6's 64-bit float variables being able to hold much bigger integers) I would call that a bug, and a serious one. In all my BASICs there is no such limitation.

[DavidS]

Also note that the use of the INT function requires that the result be stored as an integer at least for a moment in an ARM register.

If that means you can't use INT() to truncate a value greater than 2^31 to an integer (despite BASIC 6's 64-bit float variables being able to hold much bigger integers) I would call that a bug, and a serious one. In all my BASICs there is no such limitation.

No bug, all documentation of the INT function (parens optional) explicitly states for int real that real must be in the range of possible integers. The same documentation specifies that The range of possible integers ... ... values that may be stored in an integer variable.

So behaviour that allows larger values to be returned from INT would be the bug, as the language definition is fairly clear on that.

[ejolson]

Also note that the use of the INT function requires that the result be stored as an integer at least for a moment in an ARM register.

If that means you can't use INT() to truncate a value greater than 2^31 to an integer (despite BASIC 6's 64-bit float variables being able to hold much bigger integers) I would call that a bug, and a serious one. In all my BASICs there is no such limitation.

No bug, all documentation of the INT function (parens optional) explicitly states for int real that real must be in the range of possible integers. The same documentation specifies that The range of possible integers ... ... values that may be stored in an integer variable.

So behaviour that allows larger values to be returned from INT would be the bug, as the language definition is fairly clear on that.

With Matrix Brandy Basic one gets

Code: Select all

$ sbrandy

Matrix Brandy BASIC V version 1.21.19 (Linux) 01 Mar 2019

Starting with 651516 bytes free

>PRINT INT(1E3)
      1000
>PRINT INT(1E4)
     10000
>PRINT INT(1E5)
    100000
>PRINT INT(1E6)
   1000000
>PRINT INT(1E7)
  10000000
>PRINT INT(1E8)
 100000000
>PRINT INT(1E9)
1000000000
>PRINT INT(1E10)
2147483647

which seems to agree with the convert to a 32-bit integer and back theory. At the same time, Matrix Brandy Basic works fine with the classic line-numbered Fibonacci code. After the necessary uppercase conversion and switching LOG for LN the code looks like

Code: Select all

100 REM CLASSIC.BAS -- COMPUTE THE NTH FIBONACCI NUMBER
110 REM WRITTEN DECEMBER 25, 2018 BY ERIC OLSON
120 REM
130 REM THIS PROGRAM DEMONSTRATES THE EXPRESSIVENESS OF THE ORIGINAL
140 REM VERSIONS OF MICROSOFT BASIC AS MEASURED BY EXPLICITLY CODING
150 REM KARATSUBA MULTIPLICATION FOR BIG-NUMBER ARITHMETIC AND THEN
160 REM USING THE DOUBLING FORMULA
170 REM
180 REM     F(2K) = F(K)[2F(K+1)-F(K)]
190 REM   F(2K+1) = F(K+1)[2F(K+1)-F(K)]+(-1)^(K+1)
200 REM   F(2K+1) = F(K)[2F(K)+F(K+1)]+(-1)^(K)
210 REM   F(2K+2) = F(K+1)[2F(K)+F(K+1)]
220 REM
230 REM TO COMPUTE THE NTH FIBONACCI NUMBER.
231 REM
232 REM VERSION 2:  MINOR CHANGES TO OPTIMIZE THE O(N^2) MULTIPLY AND
233 REM PREVENT OVERFLOW WHEN RUNNING ON MITS BASIC.
235 REM
236 REM TO RUN THIS PROGRAM ON EARLY VERSIONS OF MICROSOFT BASIC PLEASE
237 REM REMOVE LINE 241 TO SET THE DEFAULT BACK TO SINGLE PRECISION
238 REM AND CHANGE N SO THE RESULTING FIBONACCI NUMBER FITS IN MEMORY.
240 REM
241 DEFDBL A-Z
242 Q1=0:Q2=0
250 N=4784969
251 REM N=7499
252 GOSUB 8500
270 D9=INT(N*LN((1+SQR(5))/2)/LN(B8)+7):M9=14
280 DIM M(D9*M9)
290 M0=1:INPUT "? " N:IF N=0 THEN STOP
300 A=M0:M0=M0+1+D9:B=M0:M0=M0+1+D9
310 T1=M0:M0=M0+1+D9:T2=M0:M0=M0+1+D9
400 R0=N:GOSUB 1000
420 R7=B:GOSUB 7000
430 GOTO 290
1000 REM COMPUTE NTH FIBONACCI NUMBER
1010 REM    INPUTS: R0 THE VALUE OF N
1020 REM   OUTPUTS: B THE VALUE OF F(N)
1040 IF R0<2 THEN M(B)=1:M(B+1)=R0:RETURN
1060 N1=R0:R0=INT((N1-1)/2):GOSUB 1500
1070 P1=N1-4*INT(N1/4)
1080 IF P1=1 OR P1=3 THEN 1200
1090 R1=T1:R2=A:R3=A:GOSUB 2000
1110 R1=T2:R2=T1:R3=B:GOSUB 2000
1120 R1=T1:R2=B:R3=T2:GOSUB 4000
1170 R1=B:R2=T1:GOSUB 5000:RETURN
1200 R1=T1:R2=B:R3=B:GOSUB 2000
1210 R1=T2:R2=T1:R3=A:GOSUB 3000
1220 R1=T1:R2=B:R3=T2:GOSUB 4000
1230 IF P1=3 THEN 1250
1240 R1=T1:GOSUB 6000:GOTO 1260
1250 R1=T1:GOSUB 6500
1260 R1=B:R2=T1:GOSUB 5000:RETURN
1500 REM RECURSIVE WORK FOR NTH FIBONACCI NUMBER
1510 REM    INPUTS: R0 THE VALUE OF N
1520 REM   OUTPUTS: A THE VALUE OF F(N)
1530 REM   OUTPUTS: B THE VALUE OF F(N+1)
1540 IF R0=0 THEN M(A)=1:M(A+1)=0:M(B)=1:M(B+1)=1:RETURN
1600 M(M0)=R0:M0=M0+1:R0=INT(R0/2):GOSUB 1500
1610 M0=M0-1:R0=M(M0)
1620 P1=R0-4*INT(R0/4)
1630 IF P1=1 OR P1=3 THEN 1720
1640 R1=T1:R2=B:R3=B:GOSUB 2000
1650 R1=T2:R2=T1:R3=A:GOSUB 3000
1660 R1=T1:R2=A:R3=T2:GOSUB 4000
1670 R1=A:R2=T1:GOSUB 5000
1680 R1=T1:R2=B:R3=T2:GOSUB 4000
1690 IF P1=2 THEN 1710
1700 R1=T1:GOSUB 6000:GOTO 1711
1710 R1=T1:GOSUB 6500
1711 R1=B:R2=T1:GOSUB 5000:RETURN
1720 R1=T1:R2=A:R3=A:GOSUB 2000
1730 R1=T2:R2=T1:R3=B:GOSUB 2000
1740 R1=T1:R2=B:R3=T2:GOSUB 4000
1750 R1=B:R2=T1:GOSUB 5000
1760 R1=T1:R2=A:R3=T2:GOSUB 4000
1770 IF P1=3 THEN 1790
1780 R1=T1:GOSUB 6500:GOTO 1800
1790 R1=T1:GOSUB 6000
1800 R1=A:R2=T1:GOSUB 5000:RETURN
2000 REM BIG-NUMBER ADDITION
2010 REM  OUTPUTS: R1 THE VALUE OF A+B
2020 REM   INPUTS: R2 THE VALUE OF A
2030 REM   INPUTS: R3 THE VALUE OF B
2050 IF M(R3)>M(R2) THEN I9=M(R3) ELSE I9=M(R2)
2060 FOR I=1 TO I9+1:M(R1+I)=0:NEXT I
2070 FOR I=1 TO I9:C=0:T=M(R1+I)
2080 IF I<=M(R2) THEN T=T+M(R2+I)
2090 IF I<=M(R3) THEN T=T+M(R3+I)
2110 IF T>=B8 THEN C=1:T=T-B8
2120 M(R1+I)=T:M(R1+I+1)=M(R1+I+1)+C:NEXT I
2130 M(R1)=I9+1
2140 R4=R1:GOSUB 7500
2150 RETURN
3000 REM BIG-NUMBER SUBTRACTION
3010 REM  OUTPUTS: R1 THE VALUE OF A-B
3020 REM   INPUTS: R2 THE VALUE OF A
3030 REM   INPUTS: R3 THE VALUE OF B
3050 FOR I=1 TO M(R2):M(R1+I)=0:NEXT I
3060 FOR I=1 TO M(R3):T=M(R1+I)+M(R2+I)-M(R3+I)
3070 IF T<0 THEN T=T+B8:M(R1+I+1)=M(R1+I+1)-1
3080 M(R1+I)=T:NEXT I
3090 FOR I=M(R3)+1 TO M(R2):T=M(R1+I)+M(R2+I)
3100 IF T<0 THEN T=T+B8:M(R1+I+1)=M(R1+I+1)-1
3110 M(R1+I)=T:NEXT I
3120 M(R1)=M(R2)
3130 R4=R1:GOSUB 7500
3150 RETURN
4000 REM BIG-NUMBER MULTIPLICATION
4010 REM  OUTPUTS: R1 THE VALUE OF A*B
4020 REM   INPUTS: R2 THE VALUE OF A
4030 REM   INPUTS: R3 THE VALUE OF B
4040 IF M(R2)>80 AND M(R3)>80 THEN 4300
4050 I9=M(R2)+M(R3):FOR I=1 TO I9:M(R1+I)=0:NEXT I
4070 FOR I=1 TO M(R2):FOR J=1 TO M(R3)
4080 T=M(R1+I+J-1)+M(R2+I)*M(R3+J)
4090 IF T<B7 THEN 4120
4100 M(R1+I+J-1)=T-B7
4110 M(R1+I+J)=M(R1+I+J)+B6:GOTO 4130
4120 M(R1+I+J-1)=T
4130 NEXT J:NEXT I
4140 C=0:FOR I=1 TO I9:T=M(R1+I)+C
4150 IF T<B8 THEN C=0:GOTO 4170
4160 C=INT(T/B8):T=T-B8*C
4170 M(R1+I)=T:NEXT I
4180 M(R1)=I9
4190 R4=R1:GOSUB 7500
4230 RETURN
4300 REM BIG-NUMBER KARATSUBA ALGORITHM
4310 IF M(R2)<M(R3) THEN I8=M(R3) ELSE I8=M(R2)
4320 I8=INT(I8/2)
4330 Z0=M0:M0=M0+1+2*I8+1
4332 Z2=M0:M0=M0+1+2*I8+3
4334 Z1=M0:M0=M0+1+2*I8+5
4340 Z3=M0:M0=M0+1+I8+2
4350 Z4=M0:M0=M0+1+I8+2
4360 R5=Z4:R6=R3:GOSUB 4500
4370 R5=Z3:R6=R2:GOSUB 4500
4380 GOSUB 4600:R1=Z1:R2=Z3:R3=Z4:GOSUB 4000:GOSUB 4700
4400 Q1=M(R2):IF I8<Q1 THEN M(R2)=I8
4405 Q2=M(R3):IF I8<Q2 THEN M(R3)=I8
4410 GOSUB 4600:R1=Z0:GOSUB 4000:GOSUB 4700
4420 M(R2)=Q1:M(R3)=Q2
4430 Q3=Q1-I8:Q4=Q2-I8:IF Q3<0 OR Q4<0 THEN M(Z2)=0:GOTO 8000
4440 Q1=M(R2+I8):M(R2+I8)=Q3:Q2=M(R3+I8):M(R3+I8)=Q4
4450 GOSUB 4600:R1=Z2:R2=R2+I8:R3=R3+I8:GOSUB 4000:GOSUB 4700
4460 M(R2+I8)=Q1:M(R3+I8)=Q2
4470 GOTO 8000
4500 REM ADD HIGH TO LOW
4510 REM  OUTPUTS: R5 THE SUM OF HIGH(A)+LOW(A)
4520 REM   INPUTS: R6 THE VALUE OF A
4530 REM   INPUTS: I8 THE SPLIT POINT
4540 C=0:FOR I=1 TO I8+1:T=C
4545 IF I<=I8 AND I<=M(R6) THEN T=T+M(R6+I)
4550 IF I+I8<=M(R6) THEN T=T+M(R6+I+I8)
4560 IF T>=B8 THEN C=1:T=T-B8 ELSE C=0
4570 M(R5+I)=T:NEXT I:M(R5+I8+2)=C
4590 M(R5)=I8+2:R4=R5:GOSUB 7500
4595 RETURN
4600 REM SAVE FRAME
4610 M(M0)=Z1:M0=M0+1:M(M0)=Z2:M0=M0+1:M(M0)=Z0:M0=M0+1
4620 M(M0)=I8:M0=M0+1:M(M0)=Q1:M0=M0+1:M(M0)=Q2:M0=M0+1
4630 REM SAVE PARAMETERS
4640 M(M0)=R1:M0=M0+1:M(M0)=R2:M0=M0+1:M(M0)=R3:M0=M0+1
4650 RETURN
4700 REM RESTORE FRAME
4710 GOSUB 4750
4720 M0=M0-1:Q2=M(M0):M0=M0-1:Q1=M(M0):M0=M0-1:I8=M(M0)
4730 M0=M0-1:Z0=M(M0):M0=M0-1:Z2=M(M0):M0=M0-1:Z1=M(M0)
4740 RETURN
4750 REM RESTORE PARAMETERS
4760 M0=M0-1:R3=M(M0):M0=M0-1:R2=M(M0):M0=M0-1:R1=M(M0)
4770 RETURN
5000 REM BIG-NUMBER COPY
5010 REM  OUTPUTS: R1 THE VALUE OF A
5020 REM   INPUTS: R2 THE VALUE OF A
5030 R4=R2:GOSUB 7500
5040 FOR I=1 TO M(R2):M(R1+I)=M(R2+I):NEXT I
5050 FOR I=M(R2)+1 TO M(R1):M(R1+I)=0:NEXT I
5060 M(R1)=M(R2)
5070 RETURN
6000 REM BIG-NUMBER DECREMENT
6010 REM   INPUTS: R1 THE VALUE OF A
6020 REM  OUTPUTS: R1 THE VALUE OF A-1
6040 I=1:C=1
6050 IF C=0 THEN 6080
6060 IF M(R1+I)<1 THEN M(R1+I)=B8-1 ELSE M(R1+I)=M(R1+I)-1:C=0
6070 I=I+1:GOTO 6050
6080 R4=R1:GOSUB 7500
6100 RETURN
6500 REM BIG-NUMBER INCREMENT
6510 REM   INPUTS: R1 THE VALUE OF A
6520 REM  OUTPUTS: R1 THE VALUE OF A+1
6540 M(R1)=M(R1)+1:M(R1+M(R1))=0:I=1:C=1
6550 IF C=0 THEN 6590
6560 T=M(R1+I)+1
6570 IF T>=B8 THEN M(R1+I)=T-B8 ELSE M(R1+I)=T:C=0
6580 I=I+1:GOTO 6550
6590 R4=R1:GOSUB 7500
6600 RETURN
7000 REM BIG-NUMBER PRINT
7010 REM   INPUTS: R7 THE VALUE TO PRINT
7020 IF M(R7)=0 THEN PRINT "0":RETURN
7030 FOR I=M(R7) TO 1 STEP -1
7040 S$=STR$(M(R7+I))
7045 IF MID$(S$,1,1)=" " THEN S$=MID$(S$,2,LEN(S$)-1)
7050 IF I=M(R7) OR B9<=LEN(S$) THEN 7070
7060 S$=MID$("0000000000000",1,B9-LEN(S$))+S$
7070 PRINT S$;:NEXT I:PRINT:RETURN
7500 REM BIG-NUMBER TRIM
7510 REM    INPUTS: R4 THE VALUE OF A
7520 REM   OUTPUTS: R4 THE TRIMMED VALUE OF A
7530 IF M(R4+M(R4))=0 AND M(R4)>1 THEN M(R4)=M(R4)-1:GOTO 7530
7540 RETURN
8000 REM TAIL OF KARATSUBA
8003 T4=M0-1-2*I8-5
8005 GOSUB 4630:R1=T4:R2=Z1:R3=Z0:GOSUB 3000
8010 R1=Z1:R2=T4:R3=Z2:GOSUB 3000:GOSUB 4750
8020 R4=Z1:GOSUB 7500
8030 I9=M(R2)+M(R3):I7=2*I8:C=0:FOR I=1 TO I9:T=C
8040 IF I<=M(Z0) THEN T=T+M(Z0+I)
8050 IF I>I8 AND I-I8<=M(Z1) THEN T=T+M(Z1+I-I8)
8060 IF I>I7 AND I-I7<=M(Z2) THEN T=T+M(Z2+I-I7)
8070 IF T<B8 THEN C=0:GOTO 8090
8080 C=INT(T/B8):T=T-B8*C
8090 M(R1+I)=T:NEXT I
8092 M(R1)=I9:R4=R1:GOSUB 7500
8100 M0=Z0
8120 RETURN
8500 REM GET DYNAMIC RANGE
8510 REM    OUTPUTS: F0 THE LARGEST INTEGER
8511 REM    OUTPUTS: B7 CARRY THRESHOLD
8512 REM    OUTPUTS: B8 RADIX OF LIMBS
8513 REM    OUTPUTS: B9 EXPONENT TO B8=10^B9
8520 F9=1
8530 F7=F9+1:IF F7>F9 THEN F9=F9*2:GOTO 8530
8540 F0=F9
8550 F7=F0+1:IF F7=F0 THEN F0=INT(F0/2):GOTO 8550
8560 F4=0
8565 F5=INT(0.5*(F0+F9)+0.5)
8570 F7=F5+1:IF F7=F5 THEN F9=F5 ELSE F0=F5
8580 IF F5=F4 THEN 8600
8590 F4=F5:GOTO 8565
8600 F0=F0/2
8610 B9=INT(LN(SQR(F0))/LN(10))
8620 B8=INT(EXP(LN(10)*B9)+0.5)
8630 B7=F0-2*B8^2:IF B7/B8^2<4 THEN B9=B9-1:GOTO 8620
8640 B6=INT(B7/B8):B7=B8*B6
8700 RETURN
9999 END

I suspect use of the wrong base logarithm is what leads to your overflows. A typical run looks like

Code: Select all

$ sbrandy -size 48000000 brandy.bas
? 10
55
? 50
12586269025
? 100
354224848179261915075
? 1000
43466557686937456435688527675040625802564660517371780402481729089536555417949051890403879840079255169295922593080322634775209689623239873322471161642996440906533187938298969649928516003704476137795166849228875
? 0

STOP at line 290
$ 

[ejolson]

Also note that the use of the INT function requires that the result be stored as an integer at least for a moment in an ARM register.

If that means you can't use INT() to truncate a value greater than 2^31 to an integer (despite BASIC 6's 64-bit float variables being able to hold much bigger integers) I would call that a bug, and a serious one. In all my BASICs there is no such limitation.

No bug, all documentation of the INT function (parens optional) explicitly states for int real that real must be in the range of possible integers. The same documentation specifies that The range of possible integers ... ... values that may be stored in an integer variable.

So behaviour that allows larger values to be returned from INT would be the bug, as the language definition is fairly clear on that.

One obtains a different and almost satisfying result with the C program

Code: Select all

#include <stdio.h>
#include <math.h>

int main(){
    for(double x=1e3;x<=1e32;x*=10){
        double y=floor(x);
        printf("floor(%g)=%f %s\n",
            x,y,y==x?"equal":"not equal");
    }
    return 0;
}

which produces the output

Code: Select all

$ ./a.out 
floor(1000)=1000.000000 equal
floor(10000)=10000.000000 equal
floor(100000)=100000.000000 equal
floor(1e+06)=1000000.000000 equal
floor(1e+07)=10000000.000000 equal
floor(1e+08)=100000000.000000 equal
floor(1e+09)=1000000000.000000 equal
floor(1e+10)=10000000000.000000 equal
floor(1e+11)=100000000000.000000 equal
floor(1e+12)=1000000000000.000000 equal
floor(1e+13)=10000000000000.000000 equal
floor(1e+14)=100000000000000.000000 equal
floor(1e+15)=1000000000000000.000000 equal
floor(1e+16)=10000000000000000.000000 equal
floor(1e+17)=100000000000000000.000000 equal
floor(1e+18)=1000000000000000000.000000 equal
floor(1e+19)=10000000000000000000.000000 equal
floor(1e+20)=100000000000000000000.000000 equal
floor(1e+21)=1000000000000000000000.000000 equal
floor(1e+22)=10000000000000000000000.000000 equal
floor(1e+23)=99999999999999991611392.000000 equal
floor(1e+24)=999999999999999983222784.000000 equal
floor(1e+25)=9999999999999998758486016.000000 equal
floor(1e+26)=99999999999999987584860160.000000 equal
floor(1e+27)=999999999999999875848601600.000000 equal
floor(1e+28)=9999999999999999583119736832.000000 equal
floor(1e+29)=99999999999999991433150857216.000000 equal
floor(1e+30)=999999999999999879147136483328.000000 equal
floor(1e+31)=9999999999999998509996388122624.000000 equal
floor(1e+32)=99999999999999987351763694911488.000000 equal

I wonder if Richard's BBC Basic does the same.

[DavidS]

Same results using either LOG or LN.

[ejolson]

Same results using either LOG or LN.

Did you try the Matrix Brandy Basic code just posted above?

[jahboater]

One obtains a different but still not fully satisfying result with the C program

Code: Select all

#include <stdio.h>
#include <math.h>

int main(){
    for(double x=1e3;x<=1e32;x*=10){
        double y=floor(x);
        printf("floor(%g)=%f %s\n",
            x,y,y==x?"equal":"not equal");
    }
    return 0;
}

wonder if Richard's BBC Basic does the same.

The results for this C code are correct and as expected.

floor(x) returns largest integral value not greater than x. (towards negative infinity).

Therefore if x is an integer, floor(x) will return x. And in your example x is always an integer.

For the example I think: long n = lrint(x)

If the rounded value of x cannot be stored in a long, it will report a domain error (and raise the FE_INVALID exception).