use multiple threads in lea406 computing

test on 4-core-4-hyperthread i5 and 2-core-4-hyperthread i5, the speed of
generating lunar calendar almost doubled.

c/Makefile

@@ -1,30 +1,41 @@
+
 CC = gcc
 CFLAGS = -Wall -O2
-LIBS = -lm
-SRCS = lunarcal.c lunarcalbase.c lunarcalbase.h \
-	astro.c astro.h vsop.c nutation.c julian.c \
-	lea406-full.c lea406-full.h 
-OBJS = $(SRCS:.c=.o)
-TESTSRCS = testastro.c \
-	astro.c astro.h vsop.c nutation.c julian.c \
-	lea406-full.c lea406-full.h 
-TESTOBJS = $(TESTSRCS:.c=.o)
+LIBS = -lm -lpthread
+
 LUNARCAL = lunarcal
 TESTASTRO = testastro
-EXECUTABLE = lunarcal
 
+# default target
+.PHONY : all
 all: $(LUNARCAL) $(TESTASTRO)
-	@echo all compiled
+	@echo all done!
 
-$(LUNARCAL): $(OBJS) 
-	$(CC) $(CFLAGS) -o $(LUNARCAL) $(OBJS) $(LIBS)
+OBJS =
+OBJS += astro.o
+OBJS += vsop.o
+OBJS += nutation.o
+OBJS += julian.o
+OBJS += lea406-full.o
 
-$(TESTASTRO): $(TESTOBJS) 
-	$(CC) $(CFLAGS) -o $(TESTASTRO) $(TESTOBJS) $(LIBS)
+LUNARCAL_OBJS = $(OBJS)
+LUNARCAL_OBJS += lunarcalbase.o
+LUNARCAL_OBJS += lunarcal.o
 
-.c.o:
-	$(CC) $(CFLAGS) -c $< -o $@
-	    
+TESTASTRO_OBJS = $(OBJS)
+TESTASTRO_OBJS += testastro.o
+
+$(LUNARCAL_OBJS) $(TESTASTRO_OBJS): astro.h
+lunarcalbase.o lunarcal.o: lunarcalbase.h
+lea406-full.o:  lea406-full.h
+
+$(LUNARCAL): $(LUNARCAL_OBJS)
+	$(CC) $(CFLAGS) -o $(LUNARCAL) $(LUNARCAL_OBJS) $(LIBS)
+
+$(TESTASTRO): $(TESTASTRO_OBJS)
+	$(CC) $(CFLAGS) -o $(TESTASTRO) $(TESTASTRO_OBJS) $(LIBS)
+
+
+.PHONY : clean
 clean:
 	rm -f *.o core a.out astro lunarcal testastro
-

c/astro.h

@@ -16,6 +16,8 @@
 #define SUN_SPEED   TWOPI / TROPICAL_YEAR  /* longitude change per day*/
 #define NMCOUNT  15    /* default search total 15 new moons */
 #define ISODTLEN 30    /* max length of ISO date string */
+#define MAX_THREADS 32  /* max number of threads for compute lea406-full */
+#define MAX_CPUINFO_LEN 1000  /* max line buf size when parse /proc/cpuinfo */
 
 typedef struct {
     int year;
@@ -23,6 +25,10 @@ typedef struct {
     double day;
 } GregorianDate;
 
+struct worker_param {
+    int tid;
+    double tc;           /* t in century */
+};
 
 /* Function prototypes */
 
@@ -47,6 +53,8 @@ double apparentmoon(double jd, int ignorenutation);
 
 double lea406(double jd, int ignorenutation);
 
+void *lea406worker(void *args);
+
 double nutation(double jd);
 
 double lightabbr_high(double jd);
@@ -68,5 +76,6 @@ void findnewmoons(double newmoons[], int nmcount, double startjd);
 
 double solarterm(int year, double angle);
 
-
 int findastro(int year);
+
+int cpucount(void);

c/lea406-full.c

@@ -11,10 +11,63 @@ Reference:
 */
 
 #include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
 #include <math.h>
+#include <pthread.h>
+#include <assert.h>
 #include "astro.h"
 #include "lea406-full.h"
 
+static double vlea406[MAX_THREADS];
+static int num_threads = 0;  /* number of threads for compute lea406-full */
+static int nelems_per_thread;   /* number of elements assigned to a thread */
+
+/* count logical CPU by parsing /proc/cpuinfo */
+int cpucount(void)
+{
+    FILE *fp;
+    int logic_cpu;
+    char *s;
+    char buf[MAX_CPUINFO_LEN];
+    fp = fopen("/proc/cpuinfo", "rb");
+    logic_cpu = 0;
+    while ((s = fgets(buf, MAX_CPUINFO_LEN, fp)) != NULL) {
+        if (s == strstr(buf, "processor"))
+            logic_cpu += 1;
+    }
+    logic_cpu = (logic_cpu > MAX_THREADS) ? MAX_THREADS : logic_cpu;
+    return logic_cpu;
+}
+
+
+/* the thread worker for lea406 */
+void *lea406worker(void *args)
+{
+    int tid, i, start, end;
+    double t, tm, tm2, V, arg;
+    tid = ((struct worker_param *) args)->tid;
+    t =   ((struct worker_param *) args)->tc;
+    tm = t / 10.0;
+    tm2 = tm * tm;
+
+    start = tid * nelems_per_thread;
+    end = start + nelems_per_thread;
+    end = (end > LEA406TERMS) ? LEA406TERMS : end;
+
+    V = 0.0;
+    for (i = start; i < end; i++) {
+        arg = (M_ARG[i][0] + t * (M_ARG[i][1] + M_ARG[i][2] * t)) * ASEC2RAD;
+        V +=    M_AP[i][0] * sin(arg + M_AP[i][3] * DEG2RAD)
+              + M_AP[i][1] * sin(arg + M_AP[i][4] * DEG2RAD) * tm
+              + M_AP[i][2] * sin(arg + M_AP[i][5] * DEG2RAD) * tm2;
+    }
+
+    vlea406[tid] = V;
+    return NULL;
+}
+
+
 /*
  * LEA-406 Moon Solution
  *
@@ -26,20 +79,37 @@ Reference:
 
 /* compute moon ecliptic longitude using lea406 */
 double lea406(double jd, int ignorenutation) {
-    double t, tm, tm2;
+    int rc, i;
+    double t, V;
     t = (jd - J2000) / 36525.0;
-    tm = t / 10.0;
-    tm2 = tm * tm;
 
-    double V, arg;
+    /* set number of threads number of logical CPU */
+    num_threads = (num_threads) ? num_threads : cpucount();
+    pthread_t threads[num_threads];
+    struct worker_param *tmp;
+    struct worker_param *thread_args[num_threads];
+    /* init thread args */
+    for (i = 0; i < num_threads; i++) {
+        tmp = (struct worker_param *) malloc(sizeof(struct worker_param));
+        tmp->tc = (jd - J2000) / 36525.0;
+        tmp->tid = i;
+        thread_args[i] = tmp;
+    }
+
+    nelems_per_thread = LEA406TERMS / num_threads;
+    for (i = 0; i < num_threads; i++) {
+        rc = pthread_create(&threads[i], NULL, lea406worker, thread_args[i]);
+        assert(0 == rc);
+    }
+
     V = FRM[0] + (((FRM[4] * t + FRM[3]) * t + FRM[2]) * t + FRM[1]) * t;
+    for (i = 0; i < num_threads; i++) {
+        rc = pthread_join(threads[i], NULL);
+        assert(0 == rc);
+    }
 
-    int i;
-    for (i = 0; i < LEA406TERMS; i++) {
-        arg = (M_ARG[i][0] + t * (M_ARG[i][1] + M_ARG[i][2] * t)) * ASEC2RAD;
-        V +=    M_AP[i][0] * sin(arg + M_AP[i][3] * DEG2RAD)
-              + M_AP[i][1] * sin(arg + M_AP[i][4] * DEG2RAD) * tm
-              + M_AP[i][2] * sin(arg + M_AP[i][5] * DEG2RAD) * tm2;
+    for (i = 0; i < num_threads; i++) {
+        V += vlea406[i];
     }
 
     V *= ASEC2RAD;
@@ -51,6 +121,7 @@ double lea406(double jd, int ignorenutation) {
     return V;
 }
 
+
 /* calculate the apparent position of the Moon, it is an alias to the
  * lea406 function
  */

c/lunarcalbase.c

@@ -34,8 +34,8 @@ static double newmoons[MAX_NEWMOONS];
 static double solarterms[MAX_SOLARTERMS];
 static int firstnm_offset;
 
-static int cached_year[CACHESIZE];  /* caches intermedia lunar cal*/
-static struct lunarcal *cached_lcs[CACHESIZE][MAX_DAYS];
+static struct lunarcal_cache *cached_lcs[CACHESIZE];
+static int cache_initialized = 0;  /* flag for initialized cache */
 static int cachep = 0;  /* next free location in cache */
 static int rewinded = 0;  /* cache rewinded? free pointers */
 
@@ -58,6 +58,23 @@ double normjd(double jd, double tz)
 }
 
 
+/* initialize cache storage */
+void init_cache(void)
+{
+    if (cache_initialized)
+        return;
+
+    int i;
+    struct lunarcal_cache *tmp;
+    for (i = 0; i < CACHESIZE; i++) {
+        tmp = (struct lunarcal_cache *) malloc(sizeof(struct lunarcal_cache));
+        tmp->year = -1;
+        tmp->len = -1;
+        cached_lcs[i] = tmp;
+    }
+    cache_initialized = 1;
+}
+
 void cn_lunarcal(int year)
 {
     int i, k, len1, len2;
@@ -65,6 +82,7 @@ void cn_lunarcal(int year)
     struct lunarcal *thisyear[MAX_DAYS];
     struct lunarcal *nextyear[MAX_DAYS];
     struct lunarcal *output[MAX_DAYS];
+    init_cache();
     len1 = get_cached_lc(thisyear, year);
     len2 = get_cached_lc(nextyear, year + 1);
 
@@ -75,8 +93,6 @@ void cn_lunarcal(int year)
     yend = g2jd(year, 12, 31.0);
     k = 0;
 
-    for (i = 0; i < MAX_DAYS; output[i++] = NULL)
-        ;
     for (i = 0; i < len1; i++) {
         if (thisyear[i]->jd >= ystart)
             output[k++] = thisyear[i];
@@ -94,7 +110,7 @@ int get_cache_index(int year)
 {
     int i;
     for (i = 0; i < CACHESIZE; i++) {
-        if (cached_year[i] == year)
+        if (cached_lcs[i]->year == year)
             return i;
     }
     return -1;
@@ -108,10 +124,8 @@ int get_cached_lc(struct lunarcal *p[], int year)
     double jd_nm, est_nm;
 
     if ((k = get_cache_index(year)) > -1) {
-        for (i = 0; i < MAX_DAYS; i++) {
-            if (cached_lcs[k][i] == NULL)
-                break;
-            p[i] = cached_lcs[k][i];
+        for (i = 0; i < cached_lcs[k]->len; i++) {
+            p[i] = cached_lcs[k]->lcs[i];
         }
         return i;
     }
@@ -136,12 +150,18 @@ int get_cached_lc(struct lunarcal *p[], int year)
         cachep = 0;
         rewinded = 1;
     }
+
+    if (rewinded)
+        for (i = 0; i < len; i++) {
+            free(cached_lcs[cachep]->lcs[i]);
+        }
+
     for (i = 0; i < len; i++) {
-        if (rewinded)
-            free(cached_lcs[cachep][i]);
-        cached_lcs[cachep][i] = p[i];
+        cached_lcs[cachep]->lcs[i] = p[i];
     }
-    cached_year[cachep++] = year;
+    cached_lcs[cachep]->year = year;
+    cached_lcs[cachep]->len = len;
+    cachep++;
 
     return len;
 }
@@ -155,9 +175,6 @@ int mark_month_day(struct lunarcal *lcs[])
     double lc_start, lc_end, jd, month_day1;
     struct lunarcal *lc;
 
-    for (k = 0; k < MAX_DAYS; lcs[k++] = NULL)
-        ;
-
     month = 11;
     leapmonth = find_leap();
     lc_start = newmoons[firstnm_offset];

c/lunarcalbase.h

@@ -11,6 +11,12 @@ struct lunarcal {
     int day;          /* day in Lunar Calendar */
 };
 
+struct lunarcal_cache {  /* the item in cache */
+    int year;
+    int len;             /* days count of this cached lunar calendar */
+    struct lunarcal *lcs[MAX_DAYS];   /* the cached lunar calendar */
+};
+
 /* Function prototypes */
 void cn_lunarcal(int year);
 
@@ -27,3 +33,5 @@ struct lunarcal *lcalloc(double jd);
 void print_lunarcal(struct lunarcal *p[], int len);
 
 int get_cache_index(int year);
+
+void init_cache(void);

c/testastro.c

@@ -170,7 +170,7 @@ void verify_apparent_sun_moon(void)
 
     lensun  = parsejplhorizon("jpl_sun.txt",  jplsun);
     lenmoon = parsejplhorizon("jpl_moon.txt", jplmoon);
-    step = 2;
+    step = 1;
     i = 0;
     count = 0;
     delta_sun_n = 0;