fibonacci(4784969)

[jalih]

After following John on Rasberry Pi Forums and seen some people obsessed with fibonacci(4784969), I decided to try calculating that million digit number using 8th. Actually it was really easy and result is really fast. I used locals, so getting rid of those and utilizing stack should give even better performance.

Code: [Select]


\
\ Fibonacci with million digits
\

: fibo-loop
  "a" w:@ "b" w:@ 2 n:* "a" w:@ n:- n:* "d" w:!
  "a" w:@ dup n:* "b" w:@ dup n:* n:+ "e" w:!
  "d" w:@ "a" w:!
  "e" w:@ "b" w:!
  r@ swap n:shr 1 n:band if
    "a" w:@ "b" w:@ n:+ "c" w:!
    "b" w:@ "a" w:!
    "c" w:@ "b" w:!
  then ;


locals:
: fibo
  >r
  0 "a" w:!
  1 "b" w:!
  ' fibo-loop 0 31 loop-
  rdrop
  "a" w:@ ;


: app:main
  d:msec >r
  4784969 fibo
  d:msec r> n:- >r
  "%.f" strfmt \ Convert result to just an 'int' string.
  s:len . " digits:" . cr
  dup 40 s:lsub . " upper 40 digits" . cr
  40 s:rsub . " lower 40 digits" . cr
  r> . " msec" . cr
  bye ;

Result is:

Code: [Select]

C:\temp>8th fibo.8th
1000000 digits:
1072739564180047722936481359622500432190 upper 40 digits
3407167474856539211500699706378405156269 lower 40 digits
279 msec


[John]

It seems you have the best time so far.

Thanks for taking on the challenge!

[jalih]

It seems you have the best time so far.

Here is a little bit faster (a couple of msec) version. I eliminated all local variables and used just stack.

Code: [Select]

\
\ Fibonacci with million digits
\

: fibo-loop
  >r                             \ Put loop counter on r-stack
  \ a b
  2dup 2 n:*                     \ a b a 2b
  over n:-                       \ a b a (2b-a)
  n:* -rot                       \ d=(2b-a)*a a b
  dup n:* swap dup n:* n:+       \ d=(2b-a)*a e=b*b+a*a
  r> r@ swap n:shr 1 n:band if
    dup  \ a b b
    rot  \ b b a
    n:+  \ b c=b+a
  then ;


: fibo  \  n -- fibo(n)
  >r    \ Put n on r-stack
  F0 F1   \ a b
  ' fibo-loop 0 31 loop-
  rdrop
  drop ;  \ a


: app:main
  d:msec >r
  4784969 fibo
  d:msec r> n:- >r
  n:bfloat "%.f" strfmt \ Convert result to just an 'int' string.
  s:len . " digits:" . cr
  dup 40 s:lsub . " upper 40 digits" . cr
  40 s:rsub . " lower 40 digits" . cr
  r> . " msec" . cr
  bye ;


[jalih]

I packaged a binary for RPI. Could someone test it with real Rasberry Pi hardware? I would love to know how fast it is.

[John]

Works on my Raspberry Pi 3 B.


pi@RPi3B:~/rpi8th/rpi32 $ time ./fibo
1000000 digits:
1072739564180047722936481359622500432190 upper 40 digits
3407167474856539211500699706378405156269 lower 40 digits
1946 msec

real   0m2.090s
user   0m2.044s
sys   0m0.040s
pi@RPi3B:~/rpi8th/rpi32 $

[John]

I tied to convert a Python version of the Fibonacci code challenge to ScriptBasic but I keep getting a seg fault and not sure why.

Code: Python

1. fiboA = 0
2. fiboB = 0
3. def Fibonacci(n):
4.   global fiboA
5.   global fiboB
6.   if n == 0:
7.     fiboA = 0
8.     fiboB = 1
9.     return 0
10.   Fibonacci(n // 2)
11.   if (n % 2) == 0 :
12.     t = fiboB + fiboB - fiboA
13.     fiboA = fiboA * t
14.     fiboB = fiboB * t
15.     if ( n % 4 ) == 0 : fiboB = fiboB - 1
16.     else              : fiboB = fiboB + 1
17.   else:
18.     t = fiboA + fiboA + fiboB
19.     fiboA = fiboA * t
20.     fiboB = fiboB * t
21.     if (n % 4) == 1 : fiboA = fiboA + 1
22.     else            : fiboA = fiboA - 1
23.   return fiboA
24.  
25. print Fibonacci(24) # => 46368
26.  

jrs@jrs-laptop:~/sb/examples/test$ time python hippy.py
46368

real   0m0.025s
user   0m0.020s
sys   0m0.005s
jrs@jrs-laptop:~/sb/examples/test$

print(Fast_Fibonacci(4784969))

real   0m35.497s
user   0m35.186s
sys   0m0.025s
jrs@jrs-laptop:~/sb/examples/test$

Code: ScriptBasic

1. fiboA = 0
2. fiboB = 0
3. function Fast_Fibonacci(n)
4.   if n = 0 then
5.     fiboA = 0
6.     fiboB = 1
7.     Fast_Fibonacci = 0
8.   end if
9.   Fast_Fibonacci(n\2)
10.   if (n % 2) = 0 then
11.     t = fiboB + fiboB - fiboA
12.     fiboA *= t
13.     fiboB *= t
14.     if (n % 4) = 0 then
15.       fiboB -= 1
16.     else
17.       fiboB += 1
18.     end if
19.   else
20.     t = fiboA + fiboA + fiboB
21.     fiboA *= t
22.     fiboB *= t
23.     if (n % 4) = 1 then
24.       fiboA += 1
25.     else
26.       fiboA -= 1
27.     end if
28.   end if
29.   Fast_Fibonacci = fiboA
30. end function
31.  
32. ' print(Fast_Fibonacci(4784969))
33.  
34. print Fast_Fibonacci(24),"\n"
35.  

Any ideas?

[John]

Seems I had a typo in my code.

Code: ScriptBasic

1. fiboA = 0
2. fiboB = 0
3. function Fibonacci(n)
4.   if n then
5.     Fibonacci(n \ 2)
6.     if n and 1 then
7.       t = fiboA + fiboA + fiboB
8.       fiboA *= t
9.       fiboB *= t
10.       if (n and 3) = 1 then
11.         fiboA += 1
12.       else
13.         fiboA -= 1
14.       end if
15.     else
16.       t = fiboB + fiboB - fiboA
17.       fiboA *= t
18.       fiboB *= t
19.       if n and 3 then
20.         fiboB += 1
21.       else
22.         fiboB -= 1
23.       end if
24.     end if
25.   else
26.     fiboA = 0
27.     fiboB = 1
28.   end if
29.   Fibonacci = fiboA
30. end function
31.  
32. print Fibonacci(24),"\n"
33.  

jrs@jrs-laptop:~/sb/examples/test$ time scriba hippy.sb
46368

real   0m0.011s
user   0m0.007s
sys   0m0.004s
jrs@jrs-laptop:~/sb/examples/test$ time scriba hippy.sb    <-- Fibonacci(4784969)
-nan

real   0m0.010s
user   0m0.004s
sys   0m0.006s
jrs@jrs-laptop:~/sb/examples/test$

[AIR]

Another way of doing it in Python.  I didn't write this, and don't fully understand all the calculations (I think it uses a MATRIX), but it's the fastest implementation I found.

Code: Python

1. #!/usr/bin/env python
2.  
3. def fib(n):
4.     v1, v2, v3 = 1, 1, 0  
5.     for rec in bin(n)[3:]:
6.         calc = v2*v2
7.         v1, v2, v3 = v1*v1+calc, (v1+v3)*v2, calc+v3*v3
8.         if rec=='1':    v1, v2, v3 = v1+v2, v1, v2
9.     return v2
10.  
11. res = fib(4784969)
12.  

On my Macbook Pro, 2.5Ghz 4 Core i7:
$ time ./fibo.py

real   0m1.100s
user   0m1.080s
sys   0m0.016s

On my RasPi:
$ time ./fibo.py

real   0m25.947s
user   0m25.894s
sys   0m0.050s

AIR.

[John]

The best ScriptBasic can do on 64 bit is Fibonacci(92) without overflowing  MAXINT.

Time:  0.009 seconds.

[AIR]

Jali,

Can you code this so that it outputs the actual Fibonacci number for 1000000?

[jalih]

Can you code this so that it outputs the actual Fibonacci number for 1000000?

Sure, here is a version that takes numbers as command-line parameter to calculate fibonacci number. Included is source + binaries for most systems.

[jalih]

Here is a REXX version that might be easier to read than my 8th version:

Code: [Select]

numeric digits 1000000
say fibo(4784969)
exit


fibo:
  parse arg n

  a = 0
  b = 1
  do i = 1 to length(strip(d2b(n), 'L','0'))
    d = a * (b * 2 - a)
    e = a * a + b * b
    a = d
    b = e
    if substr(strip(d2b(n), 'L','0'), i, 1) = '1' then
      do
        c = a + b
        a = b
        b = c
      end
  end
  return a


d2b: return x2b(d2x(arg(1)))


[John]

AIR,

Attached is the output to a file of the fibonacci(4784969) challenge. (from the following C code - gmp assisted)

I tried to add this as a function to the tools extension module. I have gmp-dev installed. The desired result of the function would be a SB string.

Code: C

1. #include <stdio.h>
2. #include <gmp.h>
3. #include <stdlib.h>
4.  
5. void fibo(int n) {
6.     mpz_t res;
7.     mpz_init(res);
8.  
9.     mpz_fib_ui(res, n);
10.  
11.     gmp_printf( "%Zd\n", res );
12. }
13.  
14. int main(int argc, char* argv[] )
15. {
16.     int n = 4784969;   // The first Fibonacci number with a million digits
17.  
18.     if (argc >= 2) {
19.         n = atol(argv[1]);
20.     }
21.  
22.     fibo(n);
23.     return (0);
24. }
25.  

Can you give me a hand with this? Maybe creating a new gmp extention module is a better way to go.

[AIR]

The challenge is getting the output into a string that SB can understand.

This should give you an idea:

Code: C

1. void fibo(int n) {
2.     char buf[n+1];
3.     memset(buf,0,1);
4.     mpz_t res;
5.     mpz_init(res);
6.  
7.     mpz_fib_ui(res, n);
8.  
9.     gmp_snprintf( buf,sizeof(buf),"%Zd\n", res );
10.     printf("%s\n",buf);
11. }

So now the output is in the 'buf' char array (string).  Work on returning that and you should be okay.

AIR.

[John]

This is what I have so far but doesn't return anything.

Code: C

1. /* GMP Extension Module
2. UXLIBS: -lgmp
3. */
4.  
5. #include <stdio.h>
6. #include <stdlib.h>
7. #include <string.h>
8. #include <unistd.h>
9. #include <gmp.h>
10. #include "../../basext.h"
11.  
12. /**************************
13.  Extension Module Functions
14. **************************/
15.  
16. typedef struct _ModuleObject {
17.   void *HandleArray;
18. }ModuleObject,*pModuleObject;
19.  
20.  
21. besVERSION_NEGOTIATE
22.   return (int)INTERFACE_VERSION;
23. besEND
24.  
25.  
26. besSUB_START
27.   pModuleObject p;
28.  
29.   besMODULEPOINTER = besALLOC(sizeof(ModuleObject));
30.   if( besMODULEPOINTER == NULL )return 0;
31.  
32.   p = (pModuleObject)besMODULEPOINTER;
33.   return 0;
34. besEND
35.  
36.  
37. besSUB_FINISH
38.   pModuleObject p;
39.  
40.   p = (pModuleObject)besMODULEPOINTER;
41.   if( p == NULL )return 0;
42.   return 0;
43. besEND
44.  
45.  
46. /*************
47.  GMP Functions
48. *************/
49.  
50.  
51. besFUNCTION(fibo)
52.   mpz_t res;
53.   mpz_init(res);
54.   int fval;
55.  
56.   besARGUMENTS("i")
57.     &fval
58.   besARGEND
59.  
60.   mpz_fib_ui(res, fval);
61.  
62.   besRETURN_STRING(res);
63. besEND
64.  

Code: ScriptBasic

1. declare sub fibo alias "fibo" lib "gmp"
2. PRINT fibo(24),"\n"
3.  

besRETURN_STRING(&res);     Doesn't work either.