Hi Roger,
I'm not sure this algo is really efficient on large number...
Anyway, here a little port of your code using GMP ( i removed the not used f1 and f2 functions to fit here)
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <pthread.h>
#include <gmp.h>
typedef struct arg {
mpz_t e;
mpz_t end;
mpz_t n;
arg(mpz_t _e,
mpz_t _end,
mpz_t _n)
{
mpz_init_set(e, _e);
mpz_init_set(end, _end);
mpz_init_set(n, _n);
}
} arg_t;
bool is_prime(mpz_t i)
{
return (mpz_probab_prime_p(i, 5) > 0);
}
void maxGB(mpz_t i, mpz_t* retValue)
{
mpz_t t, s, tmp;
mpz_init(t);
mpz_init(s);
mpz_init(tmp);
mpz_set(t, i);
mpz_tdiv_q_ui(t, t, 2);
mpz_set(s, t);
while(mpz_cmp_ui (t, 0) != 0)
{
mpz_sub_ui(t, t, 1);
mpz_add_ui(s, s, 1);
if(is_prime(t) && is_prime(s))
{
mpz_add(tmp, t, s);
if(mpz_cmp(tmp, i) == 0)
{
mpz_mul(tmp, t, s);
mpz_set(*retValue, tmp);
return;
}
}
}
mpz_set_ui(*retValue, 0);
}
bool Q(mpz_t n, mpz_t e, mpz_t *p, mpz_t *q)
{
bool bReturn = false;
mpz_t p_q, tmp1, tmp2;
mpz_init(p_q);
mpz_init(tmp1);
mpz_init(tmp2);
mpz_set(tmp1, e);
mpz_mul(tmp1, e, e);
mpz_mul_ui(tmp2, n, 4);
mpz_sub(tmp1, tmp1, tmp2);
mpz_abs(tmp1, tmp1);
mpz_sqrt(p_q, tmp1);
mpz_sub(*p, e, p_q);
mpz_tdiv_q_ui(*p,*p, 2);
mpz_add(*q, e, p_q);
mpz_tdiv_q_ui(*q, *q, 2);
mpz_mul(tmp1, *p, *q);
if(mpz_cmp(tmp1, n) == 0)
bReturn = true;
exit:
mpz_clear(p_q);
mpz_clear(tmp1);
mpz_clear(tmp2);
return bReturn;
}
void f(mpz_t n, mpz_t min, mpz_t max, mpz_t* retValue)
{
mpz_t e, p, q, p_e, tmp, tmp1, tmp2;
mpz_init_set_ui(e, 0);
mpz_init_set_ui(p, 0);
mpz_init_set_ui(q, 0);
mpz_init_set_ui(p_e, 0);
mpz_init_set_ui(tmp, 0);
mpz_init(tmp1);
mpz_init(tmp2);
mpz_mod_ui(tmp1, min, 2);
if(mpz_cmp_ui(tmp1, 0) != 0)
mpz_add_ui(min, min, 1);
mpz_mod_ui(tmp1, max, 2);
if(mpz_cmp_ui(tmp1, 0) != 0)
mpz_add_ui(max, max, 1);
while(mpz_cmp(max, min) > 0)
{
mpz_add(tmp1, min, max);
mpz_tdiv_q_ui(e, tmp1, 2);
if(mpz_fdiv_ui(e, 2) == 1)
mpz_add_ui(e, e, 1);
if(mpz_cmp(p_e, e) == 0)
goto exit;
maxGB(e, &tmp);
if(mpz_cmp(tmp, n) > 0)
{
mpz_sub_ui(max, e, 1);
}
else if(mpz_cmp(tmp, n) < 0)
{
mpz_add_ui(min, e, 1);
}
else
{
Q(n, e, &p, &q);
goto exit;
}
if(Q(n, e, &p, &q))
goto exit;
mpz_set(p_e, e);
}
exit:
mpz_set(*retValue, e);
mpz_clear(e);
mpz_clear(p);
mpz_clear(q);
mpz_clear(p_e);
mpz_clear(tmp);
mpz_clear(tmp1);
mpz_clear(tmp2);
}
void *worker_thread(void *args)
{
arg_t *a = (arg_t *)args;
mpz_t trials, p, q, i;
mpz_init_set_ui(trials, 0);
mpz_init_set_ui(p, 0);
mpz_init_set_ui(q, 0);
mpz_init_set(i, a->e);
while(mpz_cmp(i, a->end) < 0)
{
if(Q(a->n, i, &p, &q))
{
gmp_printf("factored: e=%Zd\tp=%Zd\tq=%Zd\ttrials=%Zd\n", i,p,q,trials);
exit(0);
}
mpz_add_ui(i, i, 2);
mpz_add_ui(trials, trials, 1);
}
return NULL;
}
#define NR_THREADS 16
void f3(mpz_t n)
{
mpz_t e, p, q, i, trials, quantum, tmp1, tmp2;
mpz_init(e);
mpz_init(quantum);
mpz_init_set_ui(p, 0);
mpz_init_set_ui(q, 0);
mpz_init_set_ui(i, 0);
mpz_init_set_ui(trials, 0);
mpz_init(tmp1);
mpz_init(tmp2);
mpz_set_ui(tmp1, 2);
mpz_add_ui(tmp2, n, 1);
f(n, tmp1, tmp2, &e);
gmp_printf("e: %Zd\n", e);
if(Q(n, e, &p, &q))
{
gmp_printf("factored: e=%Zd\tp=%Zd\tq=%Zd\ttrials=%Zd\n", e,p,q,trials);
return;
}
else
{
mpz_sub(tmp1, n, e);
if(NR_THREADS > 1)
mpz_tdiv_q_ui(quantum, tmp1, NR_THREADS - 1);
else
mpz_set(quantum, tmp1);
gmp_printf("quantum: %Zd\n", quantum);
pthread_t threads[NR_THREADS]={0};
for(int nr_threads = 0; nr_threads < NR_THREADS; nr_threads++)
{
mpz_mul_ui(tmp1, quantum, nr_threads);
mpz_add(tmp1, tmp1, e);
mpz_add_ui(tmp2, n, 1);
arg_t *a = new arg_t(tmp1 , tmp2, n);
pthread_create(&threads[nr_threads], NULL, worker_thread, (void *)a);
}
for(int nr_threads = 0; nr_threads < NR_THREADS; nr_threads++) {
pthread_join(threads[nr_threads], NULL);
}
}
printf("failed\n");
}
int main(int argc, char **argv)
{
mpz_t n;
mpz_init_set_str(n, argv[1], 10);
f3(n);
}
all the best,
orkblutt
|