brainwagon

As some of you may have noticed, I occasionally like to write small programs to compute odd little mathematical curiousities. Something I hadn’t done in a long while was to use the Sieve of Eratosthenes to compute a bunch of prime numbers. I suspect that I wrote such a program very early in my exploration of computers, maybe thirty years ago. The basic algorithm is pretty straightforward, but takes O(N) space to find the primes less than N. It’s not hard to reduce the storage to N bits, and with a trivial bit of work you can reduce it to N/2 bits, and with a little more, you can reduce it to 4/15 N bits. That was fun to work out.

But last night I did something a bit different: I implemented a simple “segmented sieve”. Basically, the idea is that to generate all the primes up to N, you find all the primes up to sqrt(N), and save those. Then, you process the rest of the numbers in buffers of sqrt(N) at a time, by sieving each buffer with your saved primes. It’s a really simple idea, but makes sieving larger numbers more practical. I implemented the simplest possible version of this last night, and set it going to compute all the prime numbers up to 1E12. Here’s the log from the run:

[chessboard] % time ./bigsieve 1000000000000 > output
20580.222u 18.741s 5:43:39.10 99.9%     0+0k 0+0io 0pf+0w

And here’s the output file:

78498 primes less than 1000000
148933 primes less than 2000000
216816 primes less than 3000000
283146 primes less than 4000000
348513 primes less than 5000000
412849 primes less than 6000000
476648 primes less than 7000000
539777 primes less than 8000000
602489 primes less than 9000000
664579 primes less than 10000000
...999989 lines deleted...
37607912018 primes less than 1000000000000

The program doesn’t include any of the basic optimizations, I suspect it woudn’t be difficult to make this thing a factor of 2-4 faster without adding a ton of code. I’ll probably see if I can do that over the next few days. It’s a useless but fun project.

Well, I dusted off my checkers program source code and compiled it on my Mac. It managed to actually solve this position: a rather common one that checkers novices can’t often convert into a win. Black to move and win. My program seems to find the winning line, but only using a search significantly deeper than I think should be necessary to find it (the forced capture is 15 ply out from the start position). I’ll investigate more.

In fact, at the moment, it’s the largest known prime, with over 9.8 million digits. As part of my pi day celebration yesterday, I was trying to review how I might speed up my C code which calculates pi to large numbers of digits. Most of the fast ways rely on fast multiplication, utilizing the FFT algorithm. I wasn’t sure how that really worked.

So… I decided to write a program to learn and test my knowledge.

Rather than compute pi though, I thought I might try a somewhat different but similar task that relied on big number computation. I remembered that Fabrice Bellard had an obfuscated C code entry that used an integer Fast Fourier Transform to print out the largest (then) known prime. It’s really quite nifty. I decided to try to implement a similar thing, but rather than starting with his obfuscated code, I decided that I’d try to use the FFTW library to do the same.

It’s a tiny bit tricky: after all, we are using a floating point FFT to multiply large integers, and getting the precision issues nailed isn’t entirely obvious. My original idea was to represent the large numbers in base 100000 (giving five decimal digits per place) but for reasons which aren’t entirely clear to me, when computing the current largest prime, I seemed to run out of precision. Limiting myself to base 10000 seems to have cured the problem.

My program seems to work (at least minimally) now. I can compute Bellard’s number (2**6972593-1) in about 9.7 seconds, and the largest known prime (2**32582657-1) in about 53 seconds. Ironically, such numbers are really easy to output in base two: they are just lots of ones. But to convert them to base ten for formatting in the way us humans like is much more complicated.

The final output looks like:

 12,457,502,601,536,945,540,085,550,157,479,950,312,279,598,515,115,184
284,367,047,566,259,111,523,599,739,738,055,975,960,661,684,593,910,041
988,688,211,130,870,620,428,490,430,485,634,271,939,241,796,764,631,759
...
... 181631 lines deleted for brevity
...
826,726,604,958,937,322,582,512,072,612,621,443,114,535,641,869,584,273
577,446,330,457,465,821,333,212,445,737,104,635,692,000,092,659,011,752
880,154,053,967,871

I have double checked this against the official value on mersenne.org, and it matches precisely. I know that the program hides at least one gross inefficiency still: I suspect that I can actually compute the large prime in just about twelve seconds. But my head hurts for the moment, and I think I’ll pass for now.

It now can put out ground tracks as well as more detailed tracking information. Just a few more lines of code.

ARISS will be visible from grid CM87ux starting in 01:36:47 at 23:08:11
  23:08:11  +0.0°  210.7° ?  21.8°N 132.1°W  AOS
  23:09:00  +3.4°  207.9° ?  24.2°N 130.0°W
  23:10:00  +8.5°  202.5° ?  27.0°N 127.2°W
  23:11:00 +15.8°  192.4° ?  29.8°N 124.3°W
  23:12:00 +26.4°  170.1° ?  32.5°N 121.2°W
  23:12:51 +32.2°  132.9° ?  34.8°N 118.4°W  MAX
  23:13:00 +31.8°  125.4° ?  35.2°N 117.9°W
  23:14:00 +22.3°   88.7° ?  37.7°N 114.4°W
  23:15:00 +13.0°   72.4° ?  40.1°N 110.6°W
  23:16:00  +6.7°   64.6° ?  42.3°N 106.5°W
  23:17:00  +2.0°   60.2° ?  44.4°N 102.1°W
  23:17:31  +0.0°   58.7° ?  45.3°N  99.7°W  LOS

P.S. Sigh, wordpress is apparently working very hard to get rid of the Unicode characters I nicely output. To the right of the azimuth/elevation should appear a sequence of arrows to indicate if you should look north, south, east or west. When I run this on my terminal, I get:

I’ll have to work on deciding whether this Unicode stuff is worth the trouble. Probably not.

Well, my plan13 library has been joined with a library that decodes grid squares and the like, and another which downloads orbital elements and stores them in an sqlite3 database. The combination of all these allows you to write simple programs like the one I illustrate below, which gives predictions of the named satellites from any grid (defaulting of course to my own gridsquare). Witness:

[chessboard] % ./nextpass -v ARISS AO-27 SO-50 AO-51
ARISS [+] will be visible from CM87 in 0:54:30 at 02:51:30
        AOS: 02:51:29  -0.0°  163.2°
        MAX: 02:54:29  +4.4°  124.1°
        LOS: 02:57:27  +0.0°   86.0°
AO-27 [+] will be visible from CM87 in 7:53:15 at 09:50:15
        AOS: 09:50:11  -0.0°   46.9°
        MAX: 09:55:01  +6.2°   87.2°
        LOS: 09:59:45  +0.0°  127.0°
SO-50 [+] will be visible from CM87 in 6:16:45 at 08:13:45
        AOS: 08:13:44  -0.0°  159.7°
        MAX: 08:19:09 +12.3°  106.6°
        LOS: 08:24:47  +0.0°   53.6°
AO-51 [+] will be visible from CM87 in 0:42:30 at 02:39:30
        AOS: 02:39:23  -0.0°  146.5°
        MAX: 02:46:48 +41.0°   69.8°
        LOS: 02:54:14  +0.0°  356.3°

My old program for drawing globes made some nice postscript output, but in reexamining the source code, I can only imagine that I was doped up on cough medicine while working on it. I started a bit on a revamp of it, starting by purloining the matrix and vector library that I used in my old raytracer and swiping the outline data from xearth. A couple of hours of TV watching and coding later, and I have the basics fleshed out:

I’m not 100% satisfied with it, since the coordinate system that it currently uses is the default one that xearth uses, and which is not the same as used by the internals of the plan13 code that I will eventually interface to it. It’s not
especially critical, but I think it makes a less clear presentation. I’ll tighten it up later.

It also doesn’t do any filling of the continent outlines, which I’m not really sure is a flaw: my goal is not to present a photorealistic view of the earth, but rather to show a clear schematic of the path of satellites as they orbit the earth. Still, I’ll give it some thought.

If you’ve listened to some of my satellite audio, you’ll notice that in addition to the callsigns, people are exchanging things that sound like “Delta Mike 41″ or “Charlie Mike 87″. These are Maidenhead gridsquares: a system of rapidly transmitting your rough location. The kind most commonly heard are the ones that are 4 characters long (two alpha, followed by two digits), but it’s also not uncommon to have them be length six (two uppercase alpha, two digits, two lower case). Converting back and forth between grids and latitude and longitude is fairly simple, there are existing programs like geoid and wwl that can do it. But I decided to code up a Python library to do it. As a simple test, it takes two grid descriptors, and determines the bearing and the distance between both points. For instance, when I run python maidensquare.py CM87 BL11, I get the following:

CM87 -> BL11: bearing 251.0° distance: 2315.3 miles
BL11 -> CM87: bearing  54.0° distance: 2315.3 miles

BL11 is the gridsquare of NH7WN in Hawaii, which is, as you can see, south west of my location. You should also note that the return trip isn’t 180 degrees opposite the outgoing trip direction. That’s because these are directions on the sphere along great circle paths, and normal geometric invariants (such as the angles of a triangle adding up to 180°) simply don’t apply on the surface of a sphere.

I’m going to add this to my growing body of Python code, and will release it someday soon. Drop me an email if you’d like to try to test it out a bit more (warning: it’s mostly for programmers).

I got a pretty good recording of the NOAA17 pass today, and converted it with my software. Turned out very nice. You can see Catalina and San Clemente Island off the coast of California. One of my better ones to date.

Addendum: I didn’t use my Radio Shack scanner for this, I used my little Yaesu VX-3R instead. I’m beginning to suspect it is a good deal more sensitive than the scanner, although it’s really hard to tell.

Well, I’ve made some headway on a project that I thought would be cool to write: porting G3RUH’s Plan 13 Satellite Prediction algorithm to a more palateable language than BASIC. I chose python, and it appears to be mostly working. It reads in the TLE orbital elements (same ones I use in “predict” or “gpredict”) and then allows you to create a bunch of satellite objects, and query their positions over time.

Here’s a screendump of a simple test program that I was running this morning:

Satellites which are above the horizon are marked in bold. They are sorted by elevation. The datafields displayed are elevation, azimuth, latitude and longitude of the subsatellite point, velocity, the Doppler velocity, and the frequency of a signal Doppler shifted from 435.845Mhz (just a value I did to check, since I was using PolySat CP3 at the time, which has APRS telemetry downlinked on that frequency). The code requires some additional cleanup, and once I have it all ready to go and documented, I’ll make it available. I think it will have a lot of uses.

I can never remember these formulas, so I wrote this program. I’m putting it here so I won’t lose it, and so others may benefit.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/* $Id$
 *
 * Written by Mark VandeWettering.
 *
 * Any copyright would be pretty stupid.  Use this code as you see
 * fit, for any purpose, commercial or non-commercial, with or
 * without attribution.
 *
 * I wrote this little test program mostly as documentation for some
 * interesting recurrence relationships which I keep forgetting and
 * have to dig up from time to time.
 *
 * To generate a continuous tone, you basically are repeatedly
 * evaluating sin(n * ang), for n = 0... however many you need (or,
 * alternatively, cos).  Writing the loop in this naive way is pretty
 * horrible though, it takes many cycles to evaluate a sine or cosine.
 * Another way to do so is to compute the complex number m = cos(i *
 * ang) (where i is the sqrt(-1)) and then maintain an accumulator which
 * is repeatedly multiplied by m.  Unfortunately, eventually numerical
 * roundoff makes the resulting vector diverge.
 *
 * Luckily, there is a third option: to use a recurrence relationship
 * for cos(n * ang) and sin(n * ang).
 *
 * From section 5.4 of Numerical Recipes: The Art of Scientific Computing:
 *
 * cos(n*ang) = 2 cos(ang) cos((n-1) * ang) - cos((n-2)*ang)
 * sin(n*ang) = 2 cos(ang) sin((n-1) * ang) - sin((n-2)*ang)
 *
 * Somewhat amazingly, these recurrence relationships are _incredibly_
 * stable.  The program below implements both of them (generating a
 * sin and a cosine in each iteration) using only two multiplies and
 * two adds per iteration.  On the 3.4Ghz Intel machine I have, this
 * entire program runs in about 11 seconds, and the vectors traced
 * out by the program are still unit vectors to the precision of
 * a normal printf.  You can try replacing the "doubles", with "floats",
 * and you will see the precision isn't as good, but might be adequate
 * for your needs.  It should be noted, however, that the program runs
 * in exactly the same time in either case (doubles and floats are both
 * single cycle operations) with my machine/compiler.
 */

/* This test program runs for a billion cycles, which at the sample rate
 * specified is just about a full day.  The idea here is to see if the
 * generated sines and * cosines continue to trace out unit vectors. */

#define NCYCLES         1000000000
#define SAMPLE_RATE     11025
#define FREQ            (800.0)

main()
{
    double ang = 2.0 * M_PI * FREQ / SAMPLE_RATE ;
    double c0, c1, c2 ;
    double s0, s1, s2 ;
    double cm =  2.0 * cos(ang) ;
    int i ;

    c0 = cos(-2.0 * ang) ;
    c1 = cos(-1.0 * ang) ;

    s0 = sin(-2.0 * ang) ;
    s1 = sin(-1.0 * ang) ;

    for (i=0; i<NCYCLES; i++) {
            c2 = cm * c1 - c0 ; /* c2 is cos(i * ang) ; */
            s2 = cm * s1 - s0 ; /* s2 is sin(i * ang) ; */
            c0 = c1 ;
            c1 = c2 ;
            s0 = s1 ;
            s1 = s2 ;
    }

    for (i=0; i<100; i++) {
            c2 = cm * c1 - c0 ; /* c2 is cos(i * ang) ; */
            s2 = cm * s1 - s0 ; /* s2 is sin(i * ang) ; */
            printf("%f %f %f\n", s2, c2, s2*s2+c2*c2) ;
            c0 = c1 ;
            c1 = c2 ;
            s0 = s1 ;
            s1 = s2 ;
    }
}

/* $Log$ */

Addendum: Jim Blinn has a chapter on circle drawing in his book, JIm Blinn’s Corner which is probably relevent. It’s been some time since I read it (and I can’t find it on my shelf at the moment), but I believe it also points you to Minksy’s circle drawing algorithm that appeared in HAKMEM. It’s clever, but doesn’t really draw circles, and is no faster than the method above. What’s really odd is that it uses less storage for temporaries, which makes it look as if it were buggy.

Addendum²: Well, I coded that program up, just for fun. It turns out that you can use it as well. Gosper explains how to compute the desired value of epsilon in the Minsky’s algorithm, and unlike the code above, it is stable when evaluated using floating point arithmetic. Unfortunately, it generates a fairly notable ellipse, which doesn’t change over repeated evaluation, but isn’t very round either. As the speed the points go around the curve decrease, the curve becomes more circular.

Still, pretty interesting.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/*
 * According to Gosper, Item 151 in HAKMEM,
 * 1-epsilon^2/2 = cos(theta), where theta
 * is the desired angular increment.
 * This means that...
 * 1 - cos(theta) = epsilon^2 / 2
 * epsilon = sqrt(2.0 * (1-cos(theta))) ;
 */

#define SAMPLE_RATE 11025
#define FREQ        800.0

main()
{
    float x = 1.0, y = 0.0 ;
    float ang = 2.0 * M_PI * FREQ / SAMPLE_RATE ;
    int i ;
    float epsilon = sqrt(2.0 * (1-cos(ang))) ;

    for (i=0; i<1000000000; i++) {
        x -= epsilon * y ;
        y += epsilon * x ;
    }
    for (i=0; i<1000; i++) {
        x -= epsilon * y ;
        y += epsilon * x ;
        printf("%f %f %f\n", x, y, sqrt(x*x+y*y)) ;
    }
}

Addendum³: I’m not sure what I really said about the “complex multiply” method really is true. Using the builtin gcc complex support, it runs in about the same time (even though I think it should have twice as many multiplies, perhaps it pipelines better?) and works just as well with respect to roundoff.

For completeness:

#include <stdio.h>
#include <stdlib.h>
#include <complex.h>
#include <math.h>

#define SAMPLE_RATE 11025
#define FREQ        800.0

main()
{
        int i ;
        complex m = cexp(2.0i * M_PI * FREQ / SAMPLE_RATE) ;
        complex acc = 1 ;

        for (i=0; i<1000000000; i++)
                acc *= m ;
        for (i=0; i<1000; i++) {
                printf("%f %f %f\n", creal(acc), cimag(acc), cabs(acc)) ;
                acc *= m ;
        }
}

Well, during a road trip with my wife to San Diego and back, I managed to begin to type up my notes for an upcoming tutorial article (cross fingers) on reception and decoding of weather satellite imagery. I basically reimplemented what I had, stripping it down to its barest essence, and trying to make it easy to understand, yet still capable of creating good imagery. Oddly enough though, in the process, I inadvertently seemed to introduce some aliasing artifacts which I must admit, are puzzling me mightily. Oh well.

This morning I decided to try to record a low 25 degree maximum pass of NOAA17 that occurred to the east. It starts out fairly noisy because I have a pretty high horizon to the northeast. The first is my “advanced” decoder, which is about six times longer than the simpler decoder I wrote.

The second is the same data file, processed with my “simple”, easy-to-follow decoder. The only serious feature it is lacking is the sync detection and rectification, which as soon as I can figure the simplest way to add, I’ll try to get in. The simple version is two pages of code, not including the data for the filter tables.

Yawn. Recorded another satellite pass. Decoded it with my software. Played with it in GIMP. Part of the pass spoiled by an oddly synchronous signal, which seemed to also be Doppler shifted. Problem on the satellite? I don’t know.

NOAA17 Pass, Dec 14, 2007

Well, I was awake for a decent daytime pass of NOAA17, so I wandered out into my front yard, and recorded the pass. It was a westward pass, covering from Canada all the way down to Baja California in the south, and was reasonably noise free over a great amount of it. I hauled it into gimp and did a bit of judicious image editing, and this is what I came up with:

I consider this to be pretty darned good for as ad-hoc as my approach to satellite reception actually is.

I haven’t had the time to record some of the NOAA17 passes during the day, but once night falls, I’ve managed to record a few passes. Tonight, I decided not to use my Yaesu VX-3R, but instead tried to record the pass from my old Radio Shack PRO60 scanner. While I recorded a bit less of the pass than last night, the overall quality seemed really pretty good.

A bit of image processing with gimp yielded the following on a night time pass.  It wasn’t particularly good with respect to noise:  on the night passes the contrast is pretty low , so you get difficulties  when you try to extract good images.  Still, I’m not displeased…