Author Topic: fibonacci(4784969)  (Read 1228 times)

Offline jalih

  • Advocate
  • Posts: 52
fibonacci(4784969)
« on: April 14, 2019, 09:25:39 PM »
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

Offline John

  • Forum Support / SB Dev
  • Posts: 2685
    • ScriptBasic Open Source Project
Re: fibonacci(4784969)
« Reply #1 on: April 15, 2019, 08:38:20 AM »
It seems you have the best time so far.

Thanks for taking on the challenge!

Offline jalih

  • Advocate
  • Posts: 52
Re: fibonacci(4784969)
« Reply #2 on: April 15, 2019, 08:50:36 AM »
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 ;
« Last Edit: May 25, 2019, 03:27:18 AM by jalih »

Offline jalih

  • Advocate
  • Posts: 52
Re: fibonacci(4784969)
« Reply #3 on: April 19, 2019, 09:05:47 AM »
I packaged a binary for RPI. Could someone test it with real Rasberry Pi hardware? I would love to know how fast it is.

Offline John

  • Forum Support / SB Dev
  • Posts: 2685
    • ScriptBasic Open Source Project
Re: fibonacci(4784969)
« Reply #4 on: April 19, 2019, 09:41:04 AM »
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 $


Offline John

  • Forum Support / SB Dev
  • Posts: 2685
    • ScriptBasic Open Source Project
Re: fibonacci(4784969)
« Reply #5 on: April 19, 2019, 02:21:04 PM »
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: Script BASIC
  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?
« Last Edit: April 19, 2019, 02:40:58 PM by John »

Offline John

  • Forum Support / SB Dev
  • Posts: 2685
    • ScriptBasic Open Source Project
Re: fibonacci(4784969)
« Reply #6 on: April 19, 2019, 02:50:13 PM »
Seems I had a typo in my code.

Code: Script BASIC
  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$
« Last Edit: April 19, 2019, 09:11:13 PM by John »

Offline AIR

  • BASIC Developer
  • Posts: 660
Re: fibonacci(4784969)
« Reply #7 on: April 19, 2019, 04:52:01 PM »
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.

Offline John

  • Forum Support / SB Dev
  • Posts: 2685
    • ScriptBasic Open Source Project
Re: fibonacci(4784969)
« Reply #8 on: April 19, 2019, 06:00:43 PM »
The best ScriptBasic can do on 64 bit is Fibonacci(92) without overflowing  MAXINT.

Time:  0.009 seconds.

Offline AIR

  • BASIC Developer
  • Posts: 660
Re: fibonacci(4784969)
« Reply #9 on: April 20, 2019, 10:04:47 AM »
Jali,

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

Offline jalih

  • Advocate
  • Posts: 52
Re: fibonacci(4784969)
« Reply #10 on: April 20, 2019, 10:55:55 AM »
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.

Offline jalih

  • Advocate
  • Posts: 52
Re: fibonacci(4784969)
« Reply #11 on: April 20, 2019, 11:04:43 AM »
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)))

Offline John

  • Forum Support / SB Dev
  • Posts: 2685
    • ScriptBasic Open Source Project
Re: fibonacci(4784969)
« Reply #12 on: April 20, 2019, 12:58:01 PM »
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.

« Last Edit: April 20, 2019, 01:53:18 PM by John »

Offline AIR

  • BASIC Developer
  • Posts: 660
Re: fibonacci(4784969)
« Reply #13 on: April 20, 2019, 02:18:49 PM »
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.

Offline John

  • Forum Support / SB Dev
  • Posts: 2685
    • ScriptBasic Open Source Project
Re: fibonacci(4784969)
« Reply #14 on: April 20, 2019, 02:39:47 PM »
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: Script BASIC
  1. declare sub fibo alias "fibo" lib "gmp"
  2. PRINT fibo(24),"\n"
  3.  

besRETURN_STRING(&res);     Doesn't work either.
« Last Edit: April 20, 2019, 02:44:53 PM by John »