fibonacci(4784969)

[AIR]

riveraa@dpi:~/Projects/Python/gmp $ time ./fibogmp.py
1000000


real    0m0.362s
user    0m0.362s
sys     0m0.000s

Is this on PC or RPI 3 B+ ?

RPI 3B+


riveraa@dpi: ~/Projects/Python/gmp $ uname -a
Linux dpi 4.14.98-v7+ #1200 SMP Tue Feb 12 20:27:48 GMT 2019 armv7l GNU/Linux


AIR.

[jalih]

RPI 3B+

Very good time! How long does it take to output the whole result?

[AIR]

Printing is always expensive....



real    0m1.895s
user    0m1.817s
sys     0m0.050s

[John]

GMP is the only way to fly. (with unreal numbers)

[jalih]

Printing is always expensive....

real    0m1.895s
user    0m1.817s
sys     0m0.050s

I made a small modification that gave about 26 percent more speed to fibonacci calculation. My current 8th version on my PC calculates and prints 1000000 digit number, if output is redirected into a file in 243 msec. Not sure, how RPI 3 B+ would perform.

[AIR]

I made a small modification that gave about 26 percent more speed to fibonacci calculation. My current 8th version on my PC calculates and prints 1000000 digit number, if output is redirected into a file in 243 msec. Not sure, how RPI 3 B+ would perform.

Are you saving as text or data?

Because if I convert the number to a byte array in GO, and save that, these are the timings I get on my MBP:

real    0m0.240s
user    0m0.232s
sys    0m0.009s

[jalih]

Are you saving as text or data?

I format the result into string as it's a big-float and I want to display just the integer part. After that I print how many digits, followed by a carriage return and the string containing actual million digit number. Output is directed into file from command-line.

[AIR]

Jali, you say you "format" the number into a string rather than say "convert".  Interesting.

[John]

It's getting to the point you can't even say fibo before the results are returned.

[AIR]

Final GO version (I'm fibonacci'd out already!):

Code: Go

1. package main
2.  
3. import (
4.     "fmt"
5.     "io/ioutil"
6.  
7.     big "github.com/ncw/gmp"
8. )
9.  
10. func Mul(x, y *big.Int) *big.Int {
11.     return big.NewInt(0).Mul(x, y)
12. }
13.  
14. func Add(x, y *big.Int) *big.Int {
15.     return big.NewInt(0).Add(x, y)
16. }
17.  
18. func fib(n int) *big.Int {
19.     v1, v2, v3 := big.NewInt(1), big.NewInt(1), big.NewInt(0)
20.     s := fmt.Sprintf("%b", n)
21.     for i := 1; i < len(s); i++ {
22.         calc := Mul(v2, v2)
23.         v1, v2, v3 = Add(Mul(v1, v1), calc), Mul(Add(v1, v3), v2), Add(Mul(v3, v3), calc)
24.         if s[i] == 49 {
25.             v1, v2, v3 = Add(v1, v2), v1, v2
26.         }
27.     }
28.     return v2
29. }
30.  
31. func main() {
32.  
33.     result := fib(4784969).String()
34.     fmt.Println(len(result))
35.  
36.     ioutil.WriteFile("output.txt", []byte(result), 0644)
37. }

[riveraa@Kiara ~/Projects/go/Fib] $ time ./Fib
1000000

real    0m0.130s
user    0m0.122s
sys    0m0.007s


    Hardware Overview:

      Model Name: MacBook Pro
      Model Identifier: MacBookPro14,1
      Processor Name: Intel Core i5
      Processor Speed: 2.3 GHz
      Number of Processors: 1
      Total Number of Cores: 2
      L2 Cache (per Core): 256 KB
      L3 Cache: 4 MB
      Hyper-Threading Technology: Enabled
      Memory: 16 GB


On my RasPi, the same program:

real    0m2.246s
user    0m2.218s
sys    0m0.040s


AIR.

[John]

Nice finish!

[AIR]

Admittedly, that laptop is pretty fast considering it's only dual core, but it's only about a year old.

Here's a more realistic type of run on my Linux server, which is running on an old Optiplex 7010 (7yr old HW, 16GB ram, 8TB storage):

riveraa@nas:~/Projects/go/Fib$ time ./Fib
1000000

real    0m0.235s
user    0m0.216s
sys    0m0.016s

Hardware:

Architecture:          x86_64
CPU op-mode(s):        32-bit, 64-bit
Byte Order:            Little Endian
CPU(s):                4
On-line CPU(s) list:   0-3
Thread(s) per core:    1
Core(s) per socket:    4
Socket(s):             1
NUMA node(s):          1
Vendor ID:             GenuineIntel
CPU family:            6
Model:                 58
Model name:            Intel(R) Core(TM) i5-3470 CPU @ 3.20GHz

[petelomax]

Code: Text

1. include mpfr.e
2.  
3. mpz res = NULL, prev, next
4. integer lastn
5. atom t0 = time()
6.  
7. function fibonampz(integer n) -- resumable, works for -ve numbers, yields mpz
8. integer absn = abs(n)
9.     if res=NULL or absn!=abs(lastn)+1 then
10.         if res=NULL then
11.             prev = mpz_init(0)
12.             res = mpz_init(1)
13.             next = mpz_init()
14.         else
15.             if n==lastn then return res end if
16.         end if
17.         mpz_fib2_ui(res,prev,absn)
18.     else
19.         if lastn<0 and remainder(lastn,2)=0 then
20.             mpz_mul_si(res,res,-1)
21.         end if
22.         mpz_add(next,res,prev)
23.         {prev,res,next} = {res,next,prev}
24.     end if
25.     if n<0 and remainder(n,2)=0 then
26.         mpz_mul_si(res,res,-1)
27.     end if
28.     lastn = n
29.     return res
30. end function
31.  
32. string s = mpz_get_str(fibonampz(4784969))
33. integer l = length(s)
34. s[40..-40] = "..."
35. ?{l,s}
36. ?elapsed(time()-t0)
37.  

[John]

Looks good Pete!

Is this written in your BASIC?

Is the BIGINT library like GMP?

How long does it take for the 1 mil digit Fibo?

[John]

AIR,

I'm trying to get the GMP extension module that I added Heater's routines to allow me to do the 1 mil fibo using the BIGINT math. I'm getting a seg fault with the return of the string. Your FIBO function code works great. Can you give me a hand adapting that same method to the other functions I added?

What has me confused is in your fibo() function you create a buffer the size of n but the resulting length is 1 million bytes. Should  I use the _ui versions of the function calls? I'm still a bit confused about initial buffer size.

Heater's C version

Code: C

1. //
2. // An experiment in doing integer arithmetic using GMP with all numbers represented by strings.
3. //
4. // By heater.
5. // Modified June 11, 2019 to use base 32 strings for intermediate results.
6. //
7. #include <stdio.h>
8. #include <gmp.h>
9. #include <stdlib.h>
10. #include <memory.h>
11. #include <string.h>
12.  
13. // Number base used for internal calculations by GMP.
14. int BASE = 32;
15.  
16. mpz_t op1;
17. mpz_t op2;
18. mpz_t res;
19.  
20. // Functions letis, addis, subis and mulis do large integer arithmetic on integers represented by strings.
21.  
22. void writeis(const char *s) {
23.     mpz_set_str (op1, s, BASE);
24.     char *buf=mpz_get_str (0, 10, op1);
25.     puts(buf);
26.     free(buf);
27. }
28.  
29. char* letis(const char* s) {
30.     return strdup(s);
31. }
32.  
33. char* addis(const char* s1, const char* s2) {
34.     mpz_set_str (op1, s1, BASE);
35.     mpz_set_str (op2, s2, BASE);
36.     mpz_add (res, op1, op2);  // result = x * y
37.     char* res_string = mpz_get_str (0, BASE, res);
38.     return res_string;
39. }
40.  
41. char* subis(const char* s1, const char* s2) {
42.     mpz_set_str (op1, s1, BASE);
43.     mpz_set_str (op2, s2, BASE);
44.     mpz_sub (res, op1, op2);  // result = x * y
45.     char* res_string = mpz_get_str (0, BASE, res);
46.     return res_string;
47. }
48.  
49. char* mulis(const char* s1, const char* s2) {
50.     mpz_set_str (op1, s1, BASE);
51.     mpz_set_str (op2, s2, BASE);
52.     mpz_mul (res, op1, op2);  // result = x * y
53.     char* res_string = mpz_get_str (0, BASE, res);
54.     return res_string;
55. }
56.  
57. char* memo[3];
58.  
59. void init_memo() {
60.     memo[0] = letis("0");
61.     memo[1] = letis("1");
62.     memo[2] = letis("1");
63. }
64.  
65. void clean_memo() {
66.     free(memo[0]);
67.     free(memo[1]);
68.     free(memo[2]);
69. }
70.  
71.  
72. // Return the n'th Fibonacci number as a decimal string for integer n
73. char* fibois (int n) {
74.     char* res;
75.     if (n <= 2) {
76.         return letis(memo[n]);
77.     }
78.  
79.     int k = (n / 2);
80.     char* fk = fibois(k);
81.     char* fk1 = fibois(k + 1);
82.     char* a;
83.     char* b;
84.     if ((n % 2) == 0) {
85.         a = addis(fk1, fk1);
86.         b = subis(a, fk);
87.         res = mulis(fk, b);
88.     } else {
89.         a = mulis(fk, fk);
90.         b = mulis(fk1, fk1);
91.         res = addis(a, b);
92.     }
93.     free(a);
94.     free(b);
95.     free(fk);
96.     free(fk1);
97.     return res;
98. }
99.  
100. int main(int argc, char* argv[]) {
101.     char* f;
102.     int n = 4784969;               // The first Fibonacci number with a million digits
103.  
104.     if (argc >= 2) {
105.         n = atoi(argv[1]);
106.     }
107.  
108.     mpz_init(op1);
109.     mpz_init(op2);
110.     mpz_init(res);
111.  
112.     init_memo();
113.  
114.     f = fibois(n);
115.     writeis(f);
116.     free(f);
117.  
118.     clean_memo();
119.  
120.     mpz_clear(op1);
121.     mpz_clear(op2);
122.     mpz_clear(res);
123.  
124.     return (0);
125. }
126.  

interface.c

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 <memory.h>
10. #include <gmp.h>
11. #include "../../basext.h"
12.  
13.  
14. int BASE = 32;
15. mpz_t op1;
16. mpz_t op2;
17. mpz_t res;
18. char* memo[3];
19.  
20.  
21. static char* let(const char* s) {
22.     return strdup(s);
23. }
24.  
25.  
26. /**************************
27.  Extension Module Functions
28. **************************/
29.  
30. typedef struct _ModuleObject {
31.   void *HandleArray;
32. }ModuleObject,*pModuleObject;
33.  
34.  
35. besVERSION_NEGOTIATE
36.   return (int)INTERFACE_VERSION;
37. besEND
38.  
39.  
40. besSUB_START
41.   pModuleObject p;
42.  
43.   besMODULEPOINTER = besALLOC(sizeof(ModuleObject));
44.   if( besMODULEPOINTER == NULL )return 0;
45.  
46.   p = (pModuleObject)besMODULEPOINTER;
47.   return 0;
48. besEND
49.  
50.  
51. besSUB_FINISH
52.   pModuleObject p;
53.  
54.   p = (pModuleObject)besMODULEPOINTER;
55.   if( p == NULL )return 0;
56.   return 0;
57. besEND
58.  
59.  
60. /*************
61.  GMP Functions
62. *************/
63.  
64.  
65. besFUNCTION(fibo)
66.   int fval;
67.  
68.   besARGUMENTS("i")
69.     &fval
70.   besARGEND
71.  
72.   char buf[fval+1];
73.   memset(buf,0,1);
74.   mpz_t res;
75.   mpz_init(res);
76.  
77.   mpz_fib_ui(res, fval);
78.  
79.   gmp_snprintf( buf,sizeof(buf),"%Zd", res );
80.  
81.   besRETURN_STRING(buf);
82. besEND
83.  
84.  
85. besFUNCTION(writeis)
86.   const char *s;
87.  
88.   besARGUMENTS("z")
89.     &s
90.   besARGEND
91.  
92.   mpz_set_str (op1, s, BASE);
93.   char *buf=mpz_get_str (0, 10, op1);
94.   puts(buf);
95.   free(buf);
96.   besRETURNVALUE = NULL;
97.  
98. besEND
99.  
100.  
101. besFUNCTION(letis)
102.   const char* s;
103.  
104.   besARGUMENTS("z")
105.     &s
106.   besARGEND
107.  
108.   besRETURN_STRING(strdup(s));
109.  
110. besEND
111.  
112.  
113. besFUNCTION(addis)
114.   const char* s1;
115.   const char* s2;
116.  
117.   besARGUMENTS("zz")
118.     &s1, &s2
119.   besARGEND
120.  
121.   mpz_set_str (op1, s1, BASE);
122.   mpz_set_str (op2, s2, BASE);
123.   mpz_add (res, op1, op2);
124.   char* res_string = mpz_get_str (0, BASE, res);
125.   besRETURN_STRING(res_string);
126.  
127. besEND
128.  
129.  
130. besFUNCTION(subis)
131.   const char* s1;
132.   const char* s2;
133.  
134.   besARGUMENTS("zz")
135.     &s1, &s2
136.   besARGEND
137.  
138.   mpz_set_str (op1, s1, BASE);
139.   mpz_set_str (op2, s2, BASE);
140.   mpz_sub (res, op1, op2);
141.   char* res_string = mpz_get_str (0, BASE, res);
142.   besRETURN_STRING(res_string);
143.  
144. besEND
145.  
146.  
147. besFUNCTION(mulis)
148.   const char* s1;
149.   const char* s2;
150.  
151.   besARGUMENTS("zz")
152.     &s1, &s2
153.   besARGEND
154.  
155.   mpz_set_str (op1, s1, BASE);
156.   mpz_set_str (op2, s2, BASE);
157.   mpz_mul (res, op1, op2);
158.   char* res_string = mpz_get_str (0, BASE, res);
159.   besRETURN_STRING(res_string);
160.  
161. besEND
162.  
163.  
164. besFUNCTION(init_memo)
165.   memo[0] = let("0");
166.   memo[1] = let("1");
167.   memo[2] = let("1");
168.   besRETURNVALUE = NULL;
169.  
170. besEND
171.  
172. besFUNCTION(clean_memo)
173.   free(memo[0]);
174.   free(memo[1]);
175.   free(memo[2]);
176.   besRETURNVALUE = NULL;
177.  
178. besEND
179.  
180.  
181. besFUNCTION(init_globals)
182.   mpz_init(op1);
183.   mpz_init(op2);
184.   mpz_init(res);
185.   besRETURNVALUE = NULL;
186.  
187. besEND
188.  
189.  
190. besFUNCTION(clear_globals)
191.   mpz_clear(op1);
192.   mpz_clear(op2);
193.   mpz_clear(res);
194.   besRETURNVALUE = NULL;
195.  
196. besEND
197.  

gmp.bas

Code: ScriptBasic

1.  
2. MODULE GMP
3.  
4. DECLARE SUB ::FIBO              ALIAS "fibo"            LIB "gmp"
5. DECLARE SUB ::BI_WRITE          ALIAS "writeis"         LIB "gmp"
6. DECLARE SUB ::BI_LET            ALIAS "letis"           LIB "gmp"
7. DECLARE SUB ::BI_ADD            ALIAS "addis"           LIB "gmp"
8. DECLARE SUB ::BI_SUB            ALIAS "subis"           LIB "gmp"
9. DECLARE SUB ::BI_MUL            ALIAS "mulis"           LIB "gmp"
10. DECLARE SUB ::BI_INIT_MEMO      ALIAS "init_memo"       LIB "gmp"
11. DECLARE SUB ::BI_CLEAN_MEMO     ALIAS "clean_memo"      LIB "gmp"
12. DECLARE SUB ::BI_INIT_GLOBALS   ALIAS "init_globals"    LIB "gmp"
13. DECLARE SUB ::BI_CLEAR_GLOBALS  ALIAS "clear_globals"   LIB "gmp"
14.  
15. END MODULE
16.  

gmpfibo.sb

Code: ScriptBasic

1. ' Return the n'th Fibonacci number as a decimal string for integer n
2.  
3. IMPORT gmp.bas
4.  
5. SPLITA STRING(4,"0") BY "" To memo
6.  
7. FUNCTION fibois(n)
8.   IF n <= 2 THEN
9.     fibis = GMP::BI_LET(memo[n])
10.   END IF
11.  
12.   k = (n / 2)
13.   fk = fibois(k)
14.   fk1 = fibois(k + 1)
15.   a = ""
16.   b = ""
17.   IF n % 2 = 0 THEN
18.     a = GMP::BI_ADD(fk1, fk1)
19.     b = GMP::BI_SUB(a, fk)
20.     res = GMP::BI_MUL(fk, b)
21.   ELSE
22.     a = GMP::BI_MUL(fk, fk)
23.     b = GMP::BI_MUL(fk1, fk1)
24.     res = GMP::BI_ADD(a, b)
25.   END IF
26.   fibois = res
27. END FUNCTION
28.  
29. ' MAIN
30.  
31. n = 4784969
32. GMP::BI_INIT_GLOBALS()
33. GMP::BI_INIT_MEMO()
34. f = fibois(n)
35. GMP::BI_WRITE(f)
36. GMP::BI_CLEAN_MEMO()
37. GMP::BI_CLEAR_GLOBALS()
38.  
39. PRINT LEN(f),"\n"
40.