Skip to content

Instantly share code, notes, and snippets.

@Softwave
Last active May 3, 2024 21:22
Show Gist options
  • Star 87 You must be signed in to star a gist
  • Fork 9 You must be signed in to fork a gist
  • Save Softwave/f61091aed8c8d8249014b5056447a698 to your computer and use it in GitHub Desktop.
Save Softwave/f61091aed8c8d8249014b5056447a698 to your computer and use it in GitHub Desktop.
Fibonacci Program

Fib

Simple fibonacci number calculator.

Usage: fib nth Fibonacci number

// Calculate Fibonacci Numbers
// Public Domain
// https://creativecommons.org/publicdomain/zero/1.0/
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <ctype.h>
#include <gmp.h>
#include <time.h>
long limit, i = 0;
int main(int argc, char *argv[])
{
// Get User Input
if (argc != 2)
{
printf("Improper input. Exiting.\n");
return -1;
}
limit = strtol(argv[1], NULL, 10);
// Setup GMP
mpz_t a, b, c;
mpz_init_set_ui(a, 1);
mpz_init_set_ui(b, 0);
mpz_init(c);
// Start timing
clock_t start_time = clock();
for (i = 0; i < limit; i++)
{
// Perform the Fibonacci Calculation
mpz_add(c, a, b);
mpz_set(a, b);
mpz_set(b, c);
}
// End timing
clock_t end_time = clock();
// Print the results to stdout
printf("Fibonacci Number %ld: ", i);
mpz_out_str(stdout, 10, b);
printf("\n");
// Cleanup
mpz_clear(a);
mpz_clear(b);
mpz_clear(c);
// Print time taaken
double time_taken = ((double) end_time - start_time) / CLOCKS_PER_SEC;
printf("Calculation Time: %f seconds\n", time_taken);
return 0;
}
@xBiggs
Copy link

xBiggs commented Oct 5, 2023

Saw your video and loved it!

I added some changes:

  • Migrated to Modern C >= C99
  • Added Error Handling
  • Removed Unused Imports
  • Added const Correctness

Compiled with gcc-13 -std=c99 -Ofast -march=native -Werror -Wall -Wextra -Wpedantic -Wshadow -Wconversion -pedantic-errors fib.c -o fib -lgmp

Screenshot 2023-10-05 075603

fib.c

#if __STDC_VERSION__ < 199901L
  #error "C99 Minimum Required"
#endif

// Calculate Fibonacci Numbers
// Public Domain
// https://creativecommons.org/publicdomain/zero/1.0/

#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
#include <time.h>
#include <errno.h>

int main(int argc, char* argv[argc + 1]) {
  // Get User Input
  if (argc != 2) {
    fprintf(stderr, "Usage: %s <limit>\n", argv[0]);
    return EXIT_FAILURE;
  }

  char* end = argv[1];
  errno = 0;
  const long limit = strtol(argv[1], &end, 10);
  if ((argv[1] == end) || *end) {
    fprintf(stderr, "Error Parsing %s\n", argv[1]);
    return EXIT_FAILURE;
  } else if (errno == ERANGE) {
    perror(argv[1]);
    return EXIT_FAILURE;
  }

  // Setup GMP
  mpz_t a, b, c;
  mpz_init_set_ui(a, 1);
  mpz_init_set_ui(b, 0);
  mpz_init(c);

  // Start timing
  const clock_t start_time = clock();
  if (start_time == (clock_t) {-1}) {
    fprintf(stderr, "Error start_time clock()\n");
    return EXIT_FAILURE;
  }

  long i = 0;
  while (i < limit) {
    // Perform the Fibonacci Calculation
    mpz_add(c, a, b);
    mpz_set(a, b);
    mpz_set(b, c);
    ++i;
  }

  // End timing
  const clock_t end_time = clock();
  if (end_time == (clock_t) {-1}) {
    fprintf(stderr, "Error end_time clock()\n");
    return EXIT_FAILURE;
  }

  // Print the results to stdout
  if (printf("Fibonacci Number %zu: ", i) < 0) return EXIT_FAILURE;
  if (!mpz_out_str(stdout, 10, b)) return EXIT_FAILURE;
  if (putchar('\n') == EOF) return EXIT_FAILURE;

  // Cleanup
  mpz_clear(a);
  mpz_clear(b);
  mpz_clear(c);

  // Print time taken
  const double time_taken = ((double) (end_time - start_time)) / (double) CLOCKS_PER_SEC;
  if (printf("Calculation Time: %lf seconds\n", time_taken) < 0) return EXIT_FAILURE;
  if (fflush(stdout) == EOF) return EXIT_FAILURE;
  return EXIT_SUCCESS;
}

makefile

cc = gcc
cc_standard = -std=c99
cc_optimization = -Ofast -march=native
cc_warnings = -Werror -Wall -Wextra -Wpedantic -Wshadow -Wconversion -pedantic-errors
cc_link = -lgmp

fib: fib.c
	${cc} ${cc_standard} ${cc_optimization} ${cc_warnings} $^ -o $@ ${cc_link}

.PHONY: clean
clean:
	rm -rf fib

@LeBigCronch
Copy link

optimization: got rid of variable c, changed to only 1 loop positional check and and 1 operation per loop cycle.

	mpz_t a, b;
	mpz_set_ui(a,0);
	mpz_set_ui(b,1);

	for(i = 1; i < limit; i++){
		if(i&1){ 
			mpz_add(b,a,b); 
		}
		else{ 
			mpz_add(a,a,b); 
		}
	}

this here will check for which variable to print, this is after the clock stops

	if(i%2){ 
		gmp_printf("%Zd",a);
	} 
	else { 
		gmp_printf("%Zd",b); 
	}

@LizzyFleckenstein03
Copy link

Immense speedup:

// Calculate Fibonacci Numbers
// Public Domain
// https://creativecommons.org/publicdomain/zero/1.0/
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <ctype.h>
#include <gmp.h>
#include <time.h>

int main(int argc, char *argv[])
{
	// Get User Input
	if (argc != 2) {
		printf("usage: %s NUM\n", argv[0]);
		return EXIT_FAILURE;
	}

	long count = strtol(argv[1], NULL, 10);

	// Setup GMP
	mpz_t a, b, p, q;
	mpz_init_set_ui(a, 1);
	mpz_init_set_ui(b, 0);
	mpz_init_set_ui(p, 0);
	mpz_init_set_ui(q, 1);

	mpz_t tmp;
	mpz_init(tmp);

   	while (count > 0) {
		if (count % 2 == 0) {
			mpz_mul(tmp, q, q);
			mpz_mul(q, q, p);
			mpz_mul_2exp(q, q, 1);
			mpz_add(q, q, tmp);

			mpz_mul(p, p, p);
			mpz_add(p, p, tmp);

			count /= 2;
		} else {
			mpz_mul(tmp, a, q);

			mpz_mul(a, a, p);
			mpz_addmul(a, b, q);
			mpz_add(a, a, tmp);

			mpz_mul(b, b, p);
			mpz_add(b, b, tmp);

			count -= 1;
		}
   	}

   	mpz_out_str(stdout, 10, b);
   	printf("\n");

	// Cleanup
   	mpz_clear(a);
   	mpz_clear(b);
   	mpz_clear(p);
   	mpz_clear(q);
   	mpz_clear(tmp);
}

This uses an algorithm described in SICP.
I also removed the timing code since one can simply use time ./fib ... to benchmark.

On my machine, my implementation can get fib(500000) in 0.01s, as opposed to 1.6s with the old implementation. I believe that now we can truly claim that "computers are fast".

@xBiggs
Copy link

xBiggs commented Oct 17, 2023

@LizzyFleckenstein03 The clock timing code was purposefully included so the IO is not counted.

@Francesco146
Copy link

Francesco146 commented Oct 27, 2023

The instruction (count % 2 == 0) can be optimized with count & 1, and the instruction count /= 2; with count *= 0.5. I've tested and:

with @LizzyFleckenstein03 implementation it got:
fib(500 000) in 0.006359
in mine implemention:
fib(500 000) in 0.001956, with -O3 in the Makefile, the time it's nearly the same, with the record of 0.001920

For the record:
fib(100 000 000) in 0.699801 and
fib(1 000 000 000) in 8.852119

I was thinking that division by two could be improved by replacing it with the instruction count >>= 1; but it resulted in an average performance degradation

here's the whole program (same Makefile):

// Calculate Fibonacci Numbers
// Public Domain
// https://creativecommons.org/publicdomain/zero/1.0/
#include <stdio.h>
#include <ctype.h>
#include <gmp.h>
#include <stdint.h>
#include <stdlib.h>
#include <time.h>

int main(int argc, char *argv[])
{
    // Get User Input
    if (argc != 2)
    {
        fprintf(stderr, "Usage: %s <number>\n", argv[0]);
        return EXIT_FAILURE;
    }
    long count = strtol(argv[1], NULL, 10);
    // Setup GMP
    mpz_t a, b, p, q;
    mpz_init_set_ui(a, 1);
    mpz_init_set_ui(b, 0);
    mpz_init_set_ui(p, 0);
    mpz_init_set_ui(q, 1);
    mpz_t tmp;
    mpz_init(tmp);
    // Start timing
    const clock_t start_time = clock();
    while (count > 0) 
    {
        if (count & 1)
        {
            mpz_mul(tmp, q, q);
            mpz_mul(q, q, p);
            mpz_mul_2exp(q, q, 1);
            mpz_add(q, q, tmp);
            mpz_mul(p, p, p);
            mpz_add(p, p, tmp);
            count *= 0.5;
        } 
        else 
        {
            mpz_mul(tmp, a, q);
            mpz_mul(a, a, p);
            mpz_addmul(a, b, q);
            mpz_add(a, a, tmp);
            mpz_mul(b, b, p);
            mpz_add(b, b, tmp);
            count -= 1;
        }
    }
    // End timing
    const clock_t end_time = clock();
    if (end_time == (clock_t){-1})
    {
        fprintf(stderr, "Error end_time clock()\n");
        return EXIT_FAILURE;
    }
    mpz_out_str(stdout, 10, b);
    printf("\n");
    // Cleanup
    mpz_clear(a);
    mpz_clear(b);
    mpz_clear(p);
    mpz_clear(q);
    mpz_clear(tmp);
    // Print time taken
    const double time_taken = ((double) (end_time - start_time)) / (double) CLOCKS_PER_SEC;
    if (printf("Calculation Time: %lf seconds\n", time_taken) < 0) 
        return EXIT_FAILURE;
    if (fflush(stdout) == EOF) 
        return EXIT_FAILURE;
    return EXIT_SUCCESS;
}

@lolzhunter
Copy link

lolzhunter commented Nov 13, 2023

one of the most basic things a compiler does is change the /2 to a bitshift operation like you wanted to do, because yes it is faster, not sure why manually doing it on your end made it slower but trust that it is doing that optimisation itself when you compile

@Francesco146
Copy link

not sure why manually doing it on your end made it slower

if i have time this weekend, I'll look into the assembly instructions to check

@LeBigCronch
Copy link

LeBigCronch commented Nov 13, 2023

not sure why manually doing it on your end made it slower

if I had to guess I'd say there's some other smaller optimizations that go along with it but when you optimize for it it didn't recognize it and doesn't implement them

@mandreamorim
Copy link

@Francesco146 did you check your output? it may be my compiler but it doesn't check to me

@Francesco146
Copy link

Francesco146 commented Nov 15, 2023

one of the most basic things a compiler does is change the /2 to a bitshift operation like you wanted to do, because yes it is faster, not sure why manually doing it on your end made it slower but trust that it is doing that optimization itself when you compile

so, i've done some research and did some experiments, here's what i got:

the compiler actually did something really stupid with the /2 and *0.5, making the assembly code unnecessarily long by doing some stuff with the floating pointer instruction:

vxorpd  xmm1, xmm1, xmm1
vcvtsi2sd xmm0, xmm1, r15
vmulsd  xmm2, xmm0, QWORD PTR .LC1[rip]
vcvttsd2si r15, xmm2
test    r15, r15
jg      .L7

so i manually inserted the count >>= 1; and ran 10'000 iterations of the code and calculated the average time:

  1. fib(100000) - count >>= 1; : 0,000117 seconds
  2. fib(100000) - count *= 0.5; : 0,000119 seconds
  3. fib(100000) - count /= 2; : 0,000118 seconds

The times are really close, but the assembly code compiled is shorter than the others, soo that's it. I'll leave the tester script, capable of testing the speed of the code and stress it with the random mode.

#!/bin/bash

# Default values
num_iterations=10000
max_fibonacci_number=100000
use_random=false

# Function to display script usage
display_help() {
    echo "Usage: $0 [OPTIONS]"
    echo "Options:"
    echo "  --random        Use random Fibonacci numbers."
    echo "  -n, --iterations NUM   Set the number of iterations (default: $num_iterations)."
    echo "  -m, --max-fibonacci NUM  Set the maximum Fibonacci number (default: $max_fibonacci_number)."
    echo "  -h, --help      Display this help message."
    exit 1
}

# Process command-line options
while [[ "$#" -gt 0 ]]; do
    case "$1" in
        --random)
            use_random=true
            ;;
        -n|--iterations)
            num_iterations="$2"
            shift
            ;;
        -m|--max-fibonacci)
            max_fibonacci_number="$2"
            shift
            ;;
        -h|--help)
            display_help
            ;;
        *)
            echo "Unknown option: $1"
            display_help
            ;;
    esac
    shift
done

make clean
make

# Run the program for the first time to initialize libraries and perform a warm-up
./fib $max_fibonacci_number &>/dev/null

# Initialize the variable to accumulate total time
total_time=0

# Run iterations and accumulate times
for ((i=1; i<=$num_iterations; i++)); do
    if [ "$use_random" = true ]; then
        # Generate a random Fibonacci number within the specified range
        fibonacci_number=$(( (RANDOM % $max_fibonacci_number) + 1 ))
    else
        # Use a fixed Fibonacci number
        fibonacci_number=$max_fibonacci_number
    fi

    # Run the command and capture the time
    exec_time=$(./fib $fibonacci_number | awk '/seconds/ {print $2}')
    
    total_time=$(echo "$total_time + $exec_time" | bc -l)

    echo -ne "fib($fibonacci_number):\t$exec_time seconds\n"
done

# Calculate the average time
average_time=$(awk "BEGIN {printf \"%.6f\", $total_time / $num_iterations}")

echo "----------------------------------------"
echo "Number of iterations: $num_iterations"
if [ "$use_random" = true ]; then
    echo "Fibonacci number range: 1 to $max_fibonacci_number (randomly generated each time)"
else
    echo "Fibonacci number: $max_fibonacci_number"
fi
echo "Average time: $average_time seconds"
echo "----------------------------------------"

also (not really important) i've added this, away from the time calculation, to the main fib.c file because we don't really need the actual value saved

    ...
    _Bool need_result = 0;
    // Get User Input
    if (argc < 2 || argc > 3)
    {
        fprintf(stderr, "Usage: %s <number> [--need-result]\n", argv[0]);
	return EXIT_FAILURE;
    }

    if (argc == 3 && strcmp(argv[2], "--need-result") == 0)
        need_result = 1;

    ...


    // Check if --need-result option is provided
    if (need_result) 
    {
        // open file to write the result (out.txt)
        FILE *fp;
        fp = fopen("out.txt", "w");
        if (fp == NULL)
        {
            fprintf(stderr, "Error opening file\n");
            return EXIT_FAILURE;
        }
        mpz_out_str(fp, 10, b);
    }
    ...

@Francesco146
Copy link

Francesco146 commented Nov 15, 2023

if I had to guess I'd say there's some other smaller optimizations that go along with it but when you optimize for it it didn't recognize it and doesn't implement them

yes, i think that too. might also be that i run this code on ARM, and the M1 is doing stuff with Rosetta. hard to know

@Francesco146
Copy link

@Francesco146 did you check your output? it may be my compiler but it doesn't check to me

not sure what are you referring to, have you installed the basic libraries?

@NoahTheBaker
Copy link

NoahTheBaker commented Nov 26, 2023

The instruction (count % 2 == 0) can be optimized with count & 1, and the instruction count /= 2; with count *= 0.5. I've tested and:

. . .
while (count > 0)
{
if (count & 1)
{
mpz_mul(tmp, q, q);
mpz_mul(q, q, p);
mpz_mul_2exp(q, q, 1);
mpz_add(q, q, tmp);
mpz_mul(p, p, p);
mpz_add(p, p, tmp);
count *= 0.5;
}
else
{
mpz_mul(tmp, a, q);
mpz_mul(a, a, p);
mpz_addmul(a, b, q);
mpz_add(a, a, tmp);
mpz_mul(b, b, p);
mpz_add(b, b, tmp);
count -= 1;
}
}
// End timing
//.......
return EXIT_SUCCESS;
}

@Francesco146 Testing your code out I found some weird behavior (unless I'm misunderstanding something).

Here's the output of the above code (fibtest.exe) compared to a program I made with that and some other code in the comment chain (fib.exe). I checked the output of fib.exe with WolframAlpha to make sure the results were the correct element of the fib sequence.

./fibtest.exe 1
0
Calculation Time: 0.000000 seconds
./fib.exe 1    
1
Calculation Time: 0.000000 seconds
Number of Digits: 1
./fibtest.exe 2
1
Calculation Time: 0.000000 seconds
./fib.exe 2
1
Calculation Time: 0.000000 seconds
Number of Digits: 1
./fibtest.exe 5
1
Calculation Time: 0.000000 seconds
./fib.exe 5
5
Calculation Time: 0.000000 seconds
Number of Digits: 1
./fibtest.exe 6
2
Calculation Time: 0.000000 seconds
./fib.exe 6
8
Calculation Time: 0.000000 seconds
Number of Digits: 1
./fibtest.exe 15
0
Calculation Time: 0.000000 seconds
./fib.exe 15
610
Calculation Time: 0.000000 seconds
Number of Digits: 3
./fibtest.exe 30
610
Calculation Time: 0.000000 seconds
./fib.exe 30    
832040
Calculation Time: 0.001000 seconds
Number of Digits: 6

The code I used is here, it's a mix from @Francesco146 and @LizzyFleckenstein03:

// Calculate Fibonacci Numbers
// Public Domain
// https://creativecommons.org/publicdomain/zero/1.0/
#include <stdio.h>
#include <ctype.h>
#include <gmp.h>
#include <stdint.h>
#include <stdlib.h>
#include <time.h>
#include <string.h>

int main(int argc, char *argv[])
{
    // Get User Input
    if (argc != 2)
    {
        fprintf(stderr, "Usage: %s <number>\n", argv[0]);
        return EXIT_FAILURE;
    }
    long count = strtol(argv[1], NULL, 10);
    // Setup GMP
    mpz_t a, b, p, q;
    mpz_init_set_ui(a, 1);
    mpz_init_set_ui(b, 0);
    mpz_init_set_ui(p, 0);
    mpz_init_set_ui(q, 1);
    mpz_t tmp;
    mpz_init(tmp);
    // Start timing
    const clock_t start_time = clock();
    while (count > 0) {
		if (count % 2 == 0) {
			mpz_mul(tmp, q, q);
			mpz_mul(q, q, p);
			mpz_mul_2exp(q, q, 1);
			mpz_add(q, q, tmp);

			mpz_mul(p, p, p);
			mpz_add(p, p, tmp);

			count /= 2;
		} else {
			mpz_mul(tmp, a, q);

			mpz_mul(a, a, p);
			mpz_addmul(a, b, q);
			mpz_add(a, a, tmp);

			mpz_mul(b, b, p);
			mpz_add(b, b, tmp);

			count -= 1;
		}
   	}
    // End timing
    const clock_t end_time = clock();
    if (end_time == (clock_t){-1})
    {
        fprintf(stderr, "Error end_time clock()\n");
        return EXIT_FAILURE;
    }
    mpz_out_str(stdout, 10, b);
    unsigned int digits = strlen(mpz_get_str(NULL, 10, b));
    printf("\n");
    // Cleanup
    mpz_clear(a);
    mpz_clear(b);
    mpz_clear(p);
    mpz_clear(q);
    mpz_clear(tmp);
    // Print time taken
    const double time_taken = ((double) (end_time - start_time)) / (double) CLOCKS_PER_SEC;
    if (printf("Calculation Time: %lf seconds\n", time_taken) < 0) 
        return EXIT_FAILURE;
    if (fflush(stdout) == EOF) 
        return EXIT_FAILURE;
    printf("Number of Digits: %u\n", digits);
	return EXIT_SUCCESS;
}

Which also gave me these times:

./fib.exe 100000
Calculation Time: 0.000000 seconds
Number of Digits: 20899
./fib.exe 1000000
Calculation Time: 0.008000 seconds
Number of Digits: 208988
./fib.exe 10000000
Calculation Time: 0.097000 seconds
Number of Digits: 2089877
./fib.exe 100000000
Calculation Time: 1.316000 seconds
Number of Digits: 20898764
./fib.exe 1000000000
Calculation Time: 16.641000 seconds
Number of Digits: 208987640

Fun!

@Francesco146
Copy link

@NoahTheBaker not sure what's going on here, can you test the latest code? that version is a bit outdated

@NoahTheBaker
Copy link

@Francesco146 I believe I used the latest one you posted (Oct 27th), at least that's the last one I see aside from your timing script? I could be pulling a dumb

@Francesco146
Copy link

@Francesco146 I believe I used the latest one you posted (Oct 27th), at least that's the last one I see aside from your timing script? I could be pulling a dumb

after the bit-shift analysis, I replaced the 0.5 with that, I'll send you the full version which I ask you to test right here. perhaps with some optimizations I made a simplification that changes the actual result. I would be happy if you decided to test it in order to exclude quick but incorrect calculations.

@NoahTheBaker
Copy link

NoahTheBaker commented Nov 26, 2023

@Francesco146 I believe I used the latest one you posted (Oct 27th), at least that's the last one I see aside from your timing script? I could be pulling a dumb

after the bit-shift analysis, I replaced the 0.5 with that, I'll send you the full version which I ask you to test right here. perhaps with some optimizations I made a simplification that changes the actual result. I would be happy if you decided to test it in order to exclude quick but incorrect calculations.

3rd Edit: Fixed a mistake I made
When running your code it gives the following values:

./fibtest.exe 15      
fib(15): 0.000000 seconds
0
./fibtest.exe 30
fib(30): 0.000000 seconds
610
./fibtest.exe 1000
fib(1000): 0.000000 seconds
700601527838707533318724871765523284890968696066716357168446685359555050818763462119969522887130023714
./fibtest.exe 10000
fib(10000): 0.000000 seconds
511865913905822858393252140857692303019370496328583286846703611276316160829121874839117811096377890223341760784520464441363883635008221595005409607506796552046773590457456727956781799070386310526276130521173224078490428727580439211685193830890011304948013193753096345077838097891781185513053353450228725836492609250758663599068124138189622708238726010261183052860309857905748834

But running my code gives:

./fib.exe 15    
Fibonacci Value: 610
Number of Digits: 3
Calculation Time: 0.0000 seconds
Total Time Taken: 0.0020 seconds
./fib.exe 30
Fibonacci Value: 832040
Number of Digits: 6
Calculation Time: 0.0000 seconds
Total Time Taken: 0.0020 seconds
./fib.exe 1000     
Fibonacci Value: 43466557686937456435688527675040625802564660517371780402481729089536555417949051890403879840079255169295922593080322634775209689623239873322471161642996440906533187938298969649928516003704476137795166849228875
Number of Digits: 209
./fib.exe 10000
Fibonacci Value: 33644764876431783266621612005107543310302148460680063906564769974680081442166662368155595513633734025582065332680836159373734790483865268263040892463056431887354544369559827491606602099884183933864652731300088.....
Number of Digits: 2090
Calculation Time: 0.0000 seconds
Total Time Taken: 0.0520 seconds

And checking with WolframAlpha gives:

fib(15) = 610
fib(30) = 832040
fib(1000) = 43466557686937456435688527675040625802564660517371780402481729089536555417949051890403879840079255169295922593080322634775209689623239873322471161642996440906533187938298969649928516003704476137795166849228875
fib(10000) = 33644764876431783266621612005107543310302148460680063906564769974680081442166662368155595513633734025582065332680836159373734790483865268263040892463056431887354544369559827491606602099884183933864652731300088.....

2nd Edit: Here is the current version that I am using with verified Fibonacci values up to F(100,000,000), but I cannot check further than that as WolframAlpha seems to break.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment