※上記の広告は60日以上更新のないWIKIに表示されています。更新することで広告が下部へ移動します。

積分

/* 
 * int_seq.c
 *
 * compile:
 *   gcc int_seq.c -lm
 *
 */
 
#include <pthread.h>
#include <unistd.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>
 
typedef struct f_args
{
  double (*f)(double);
  double dest;
  double a;
  double b;
  int N;
} f_args, * f_args_t;
 
void* integrate(void * _arg)
{
  f_args_t arg = _arg;
  double (*f)(double) = arg->f;
  double a = arg->a;
  double b = arg->b;
  int N = ((f_args_t)_arg)->N;
  int i;
  double s = 0.0;
  for (i = 0; i < N; i++) {
    double x      = a + (b - a) *  i      / N;
    double x_next = a + (b - a) * (i + 1) / N;
    double dx = x_next - x;
    s += f(x) * dx;
  }
  arg->dest = s;
  return _arg;
}
 
double cur_time()
{
  struct timeval tp[1];
  gettimeofday(tp, NULL);
  return tp->tv_usec * 1.0E-6 + tp->tv_sec;
}
 
int main()
{
  //int n_threads = 2;
  //pthread_t *tids = (pthread_t *)malloc(sizeof(pthread_t *) * n_threads);
  //f_args *args = (f_args *)malloc(sizeof(f_args_t *) * n_threads);
  pthread_t tids[2];
  double a = 0.0;
  double b = M_PI / 2.0;
  double t0 = cur_time();
  //int N = 100000000;
  f_args args1 = {sin, 0.0, a, b/2, 50000000};
  f_args args2 = {sin, 0.0, b/2, b, 50000000};
  double s;
  pthread_create(&tids[0], NULL, integrate, &args1);
  pthread_create(&tids[1], NULL, integrate, &args2);
  pthread_join(tids[0], NULL);
  pthread_join(tids[1], NULL);
  s = args1.dest + args2.dest;
  double t1 = cur_time();
  printf("\\int_[%f,%f] sin(x) dx = %f\n", a, b, s);
  printf("%.6f sec\n", t1 - t0);
  return 0;
}

素数

C言語版

/* 
 * prime_seq.c
 *
 *
 */
 
#include <stdio.h>
#include <sys/time.h>
#include <stdlib.h>
#include <pthread.h>
#include <unistd.h>
 
typedef struct prime_task
{
  int a;
  int b;
} *prime_task_t;
 
typedef struct f_args
{
  prime_task_t task;
  int thread_id;
} f_args, *f_args_t;
 
volatile int task_counter = 0;
volatile int primes_counter = 0;
int n_prime_task = 100;
pthread_mutex_t m;
 
#undef DEBUG
 
int check_prime(int n)
{
  int d;
  for (d = 2; d * d <= n; d++) {
    if (n % d == 0) return 0;
  }
#ifdef DEBUG
  printf("%d is prime\n", n);
#endif
  return 1;
}
 
int check_primes(int a, int b)
{
  int n;
  int s = 0;
  for (n = a; n < b; n++) {
    s += check_prime(n);
  }
  return s;
}
 
void* check_primes_th(void * _task)
{
  prime_task_t task = ( (f_args_t)_task)->task;
  int a,b;
  while(1)
    {
      pthread_mutex_lock(&m);
      if(task_counter >= n_prime_task)
	break;
      a = task[task_counter].a;
      b = task[task_counter].b;
      ++task_counter;
      pthread_mutex_unlock(&m);
#ifdef DEBUG
      printf("task : %d, thread : %d\n",task_counter-1, ( (f_args_t)_task)->thread_id);
#endif
      primes_counter += check_primes(a,b);
    }
  pthread_mutex_unlock(&m);
  return NULL;
}
 
double cur_time()
{
  struct timeval tp[1];
  gettimeofday(tp, NULL);
  return tp->tv_usec * 1.0E-6 + tp->tv_sec;
}
 
int main()
{
  int N = 10000000;
  pthread_t tids[2];
  f_args args[2];
  prime_task_t pt = (prime_task_t)malloc(sizeof(struct prime_task) * n_prime_task);
  int t;
  for(t=0; t < n_prime_task; ++t)
    {
      pt[t].a = (N-1)*  t  /n_prime_task + 2;
      pt[t].b = (N-1)*(t+1)/n_prime_task + 1;
    }
  pthread_mutex_init(&m, NULL);
  args[0].task = args[1].task = pt;
  args[0].thread_id = 1;
  args[1].thread_id = 2;
  double t0 = cur_time();
  pthread_create(&tids[0], NULL, check_primes_th, &args[0]);
  pthread_create(&tids[1], NULL, check_primes_th, &args[1]);
  pthread_join(tids[0], NULL);
  pthread_join(tids[1], NULL);
  double t1 = cur_time();
  printf("%d primes in %f sec\n", primes_counter, t1 - t0);
  return 0;
}

C++に書き換え

//#define DEBUG // デバッグ用画面出力の制御
 
#include <cstdio>
#include <sys/time.h>
#include <deque>
#include <vector>
#include <pthread.h>
#include <unistd.h>
 
inline double cur_time()
{
  struct timeval tp[1];
  gettimeofday(tp, NULL);
  return tp->tv_usec * 1.0E-6 + tp->tv_sec;
}
 
struct PrimeTask
{
  int a, b;
  PrimeTask(int a, int b) : a(a), b(b){};
};
 
 
extern "C"
{
  typedef void* (*thread_run_t)(void*);
};
 
class CheckPrimes
{
private:
  unsigned int n_thread;
  unsigned int n_task;
  pthread_mutex_t m;
  volatile int primes_counter;
  std::deque<PrimeTask> task;
 
  static int check_prime(int n)
  {
    for (int d = 2; d * d <= n; d++) {
      if (n % d == 0) return 0;
    }
    return 1;
  };
 
  static int check_primes(int a, int b)
  {
    int n;
    int s = 0;
    for (n = a; n < b; n++) {
      s += check_prime(n);
    }
    return s;
  };
 
  void check_primes_th()
  {
    int a,b;
    while(1)
      {
	pthread_mutex_lock(&m);
	if(task.empty())
	  break;
	a = task.back().a;
	b = task.back().b;
	task.pop_back();
	pthread_mutex_unlock(&m);
	int x = check_primes(a,b);
	primes_counter += x;
#ifdef DEBUG
	pthread_t xtid = pthread_self();
	printf("task left : %d, thread : %#x, [%d,%d], detected : %6d, total : %6d\n",
	       task.size(), (unsigned int)xtid, a, b, x, primes_counter);
#endif
      }
    pthread_mutex_unlock(&m);
  };
 
  static void run(void * pParam)
  {
    CheckPrimes *pCheckPrimes = (CheckPrimes *)pParam;
    pCheckPrimes->check_primes_th();
  };
 
public:
  CheckPrimes(unsigned int threads, unsigned int tasks)
    : n_thread(threads), n_task(tasks), primes_counter(0)
  {
    pthread_mutex_init(&m, NULL);
  };
  int operator()(int Nmax)
  {
    std::vector<pthread_t> tid;
    tid.resize(n_thread);
    task.clear();
    primes_counter = 0;
 
    for(unsigned int t=0; t < n_task; ++t)
      task.push_back(PrimeTask((Nmax-1)*  t  /n_task + 2,
			       (Nmax-1)*(t+1)/n_task + 1)
		     );
 
    for(std::vector<pthread_t>::iterator it=tid.begin();
	it != tid.end(); ++it)
      pthread_create(&(*it), NULL, (thread_run_t)CheckPrimes::run, this);
    for(std::vector<pthread_t>::iterator it=tid.begin();
	it != tid.end(); ++it)
      pthread_join(*it, NULL);
    return primes_counter;
  };
};
 
int main()
{
  int N = 10000000;
  double t0 = cur_time();
  int primes = CheckPrimes(2, 100)(N);
  double t1 = cur_time();
  printf("%d primes in %f sec\n", primes, t1 - t0);
  return 0;
}
|新しいページ|検索|ページ一覧|RSS|@ウィキご利用ガイド | 管理者にお問合せ
|ログイン|