(정수론/C/GMP) 소수 구하는 함수 비교

C언어에서의 소수 관련 함수를 어떻게 gmp 라이브러리 함수로 바꿀 수 있는지 한번 살펴보자.

isPrime과 getPrimes 함수를 위주로 살펴보자.

 

#include "gmp.h"
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <math.h>
#include <limits.h>

int isPrime(unsigned long num) {
  unsigned long i, root_num;
  root_num = (int)sqrt(num) + 1;
  if (num <= 1) return 0;  // 1 이하라면 소수 아님
  if (num % 2 == 0) return 0; // 2의 배수라면 소수 아님

  for (i = 3; i < root_num; i+=2)
    if (num % i == 0)
      return 0;
    //printf("i = %lu", i);
  return 1;
}


int isPrime2(mpz_t num){
  mpz_t root_val, i, r;
  mpz_inits(root_val, i, r, NULL);


  if (mpz_cmp_ui(num, 1) != 1) // 1 이하면 소수 아님
    return 0;
  // 2로 나눠진다면 소수 아님, 리턴 0
  mpz_fdiv_r_ui(r, num, 2);      // 2로 나눠진다면 소수 아님, 리턴 0
  if (mpz_cmp_ui(r, 0) == 0) return 0;

  mpz_sqrt (root_val, num);
  mpz_add_ui(root_val, root_val, 1);
  mpz_set_ui(i, 3);    // i = 3
  for (; mpz_cmp(root_val, i);) {  // root_val > i 인 동안
    mpz_fdiv_r(r, num, i);  // r = num % i
    if (mpz_cmp_ui(r, 0) == 0)  // num % i == 0
      return 0;
    // num이 2의 배수인 경우는 제외되었으므로 홀수만으로 나눔
    mpz_add_ui(i, i, 2);  // i = i + 2
  }
  return 1;
}

void getDivisor(mpz_t num) {
  mpz_t root_val, i, q, r;
  mpz_inits(root_val, i, q, r, NULL);
  if (mpz_cmp_ui(num, 1) != 1) { // 1 이하면
    printf("You must input a number lager than 1!");
    return;
  }
  mpz_sqrt(root_val, num);
  mpz_add_ui(root_val, root_val, 1);
  mpz_set_ui(i, 2);
  
  while (mpz_cmp(root_val, i) >= 0) { // root_val >= i 인 동안
    mpz_fdiv_r(r, num, i);
    if (mpz_cmp_ui(r, 0) == 0) { // r == 0
      mpz_fdiv_q(num, num, i);
      //재설정된 num의 sqrt값을 root_val로 설정
      mpz_sqrt(root_val, num);
      gmp_printf("%Zd, %Zd\n\n", i, num);
    }
    else
      if (mpz_cmp_ui(i, 2) != 0)  // i가 2가 아니라면
        mpz_add_ui(i, i, 2);      // i가 3 이상의 홀수이므로 2를 더하여 홀수로만 검사함
      else
        mpz_add_ui(i, i, 1);    // i가 2라면 3으로 놓음
  }
}

unsigned long getPrimes(unsigned long num) {
    // 속도를 위해 소수 2는 제외하고 3부터 홀수만 검사함. 2외에는 짝수에 소수 없음.
    for (unsigned long i = 3; i <= num; i+=2) {
        if(isPrime(i))   //소수면 출력하거나 개수를 하나 올린다.
            printf("%lu\n", i);
    }
}


void getPrimes2(mpz_t num) {
  mpz_t i, q, r;
  mpz_inits(i, q, r, NULL);
  mpz_set_ui(i, 3);   // 소수 2는 나중에 추가. 홀수만 검사하기 위해서

  while (mpz_cmp(num, i) >= 0) {  // num >= i 인 동안
    if(isPrime2(i)) 
      gmp_printf("%Zd\t", i);
    mpz_add_ui(i, i, 2); // i += 2
  }
}


int main (int argc, char *argv[]) {
  mpz_t n, m;
  mpz_inits(n, m, NULL);

  int MAX_NUM  = 128;   // 128자리수 아래는 각기 10^8, 10^16...
  char *big_8 = "100000000";
  char *big_16 = "10000000000000000";
  char *big_32 = "100000000000000000000000000000000";
  char *big_64 = "10000000000000000000000000000000000000000000000000000000000000000";
  char *big_128 = "100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000";
  char *big_256 = "10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000";
  char *big_512 = "100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000";
  char *big_1024 = "10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000";
  char *test_num = "14982300945943290918090909454983788547895985895893455426";
  
  unsigned long num = ULONG_MAX;      // 18446744073709551615

  //mpz_set_str(n, big_8, 10);
  //printf("%lu\n", ULONG_MAX);
  //getPrimes2(n);
  mpz_set_str(m, test_num, 10);
  getDivisor(m);

  //getPrimes(num);
  //printf("%u\n", UINT_MAX);
  
  return 1;
}

출력:

2, 7491150472971645459045454727491894273947992947946727713

1431690181, 5232382377407424196789580920854198603984134573