Investigating floats

float variables are, on most current hardware, instrumented as as single precision IEEE 754 value. We know that floating point values only approximate real values, but it is often hard to visualise. When we have the binary decimal 0.12345 we are really saying something like 1*10-1 + 2*10-2 + 3*10-3 .... Thus a binary decimal 0.10110 describes 1*2-1 + 0*2-2 + 1*2-3 ...

So clearly, with one decimal digit we can describe the range 0.0 - 0.9, but with one binary digit we can describe only 0 and 1/2 = 0.5. The fractions we can represent exactly are the dyadic rationals. Other rationals repeat (just like in base 10) and others are irrational. So for example, 0.1 can not be repesented exactly without infinite precision.

The program below illustrates how a single precision floating point value works by showing what bits are set and adding them up.

For example

$ ./float .5
0.500000 = 1 * (1/2^0) * 2^-1
0.500000 = 1 * (1/1) * 2^-1
0.500000 = 1 * 1 * 0.500000
0.500000 = 0.5
$ ./float .1
0.100000 = 1 * (1/2^0 + 1/2^1 + 1/2^4 + 1/2^5 + 1/2^8 + 1/2^9 + 1/2^12 + 1/2^13 + 1/2^16 + 1/2^17 + 1/2^20 + 1/2^21 + 1/2^23) * 2^-4
0.100000 = 1 * (1/1 + 1/2 + 1/16 + 1/32 + 1/256 + 1/512 + 1/4096 + 1/8192 + 1/65536 + 1/131072 + 1/1048576 + 1/2097152 + 1/8388608) * 2^-4
0.100000 = 1 * 1.60000002384 * 0.062500
0.100000 = 0.10000000149

I use doubles to show how the errors creep in.

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

/* return 2^n */
int two_to_pos(int n)
{
    if (n == 0)
        return 1;
    return 2 * two_to_pos(n - 1);
}

double two_to_neg(int n)
{
    if (n == 0)
        return 1;
    return 1.0 / (two_to_pos(abs(n)));
}

double two_to(int n)
{
    if (n >= 0)
        return two_to_pos(n);
    if (n < 0)
        return two_to_neg(n);
    return 0;
}

/* Go through some memory "m" which is the 24 bit significand of a
   floating point number.  We work "backwards" from the bits
   furthest on the right, for no particular reason. */
double calc_float(int m, int bit)
{
    /* 23 bits; this terminates recursion */
    if (bit > 23)
        return 0;

    /* if the bit is set, it represents the value 1/2^bit */
    if ((m >> bit) & 1)
        return 1.0L/two_to(23 - bit) + calc_float(m, bit + 1);

    /* otherwise go to the next bit */
    return calc_float(m, bit + 1);
}

int main(int argc, char *argv[])
{
    float f;
    int m,i,sign,exponent,significand;

    if (argc != 2)
    {
        printf("usage: float 123.456\n");
        exit(1);
    }

    if (sscanf(argv[1], "%f", &f) != 1)
    {
        printf("invalid input\n");
        exit(1);
    }

    /* We need to "fool" the compiler, as if we start to use casts
       (e.g. (int)f) it will actually do a conversion for us.  We
       want access to the raw bits, so we just copy it into a same
       sized variable. */
    memcpy(&m, &f, 4);

    /* The sign bit is the first bit */
    sign = (m >> 31) & 0x1;

    /* Exponent is 8 bits following the sign bit */
    exponent = ((m >> 23) & 0xFF) - 127;

    /* Significand fills out the float, the first bit is implied
       to be 1, hence the 24 bit OR value below. */
    significand = (m & 0x7FFFFF) | 0x800000;

    /* print out a power representation */
    printf("%f = %d * (", f, sign ? -1 : 1);
    for(i = 23 ; i >= 0 ; i--)
    {
        if ((significand >> i) & 1)
            printf("%s1/2^%d", (i == 23) ? "" : " + ",
                   23-i);
    }
    printf(") * 2^%d\n", exponent);

    /* print out a fractional representation */
    printf("%f = %d * (", f, sign ? -1 : 1);
    for(i = 23 ; i >= 0 ; i--)
    {
        if ((significand >> i) & 1)
            printf("%s1/%d", (i == 23) ? "" : " + ",
                   (int)two_to(23-i));
    }
    printf(") * 2^%d\n", exponent);

    /* convert this into decimal and print it out */
    printf("%f = %d * %.12g * %f\n",
           f,
           (sign ? -1 : 1),
           calc_float(significand, 0),
           two_to(exponent));

    /* do the math this time */
    printf("%f = %.12g\n",
           f,
           (sign ? -1 : 1) *
           calc_float(significand, 0) *
           two_to(exponent)
        );

    return 0;
}

IEC 60027-2 Units

It has been pointed out to me that in a previous post I incorrectly suggested the Ki/Mi/Gi style of units for describing binary multiples has something to do with SI units. These style of units are officially described in IEC 60027-2: Letter symbols to be used in electrical technology - Part 2: Telecommunications and electronics, thus you might like to refer to them as "IEC Units" or similar.

The NIST reference has this to say:

It is important to recognise that the new prefixes for binary multiples are not part of the International System of Units (SI), the modern metric system. However, for ease of understanding and recall, they were derived from the SI prefixes for positive powers of ten.

It is a considerable faux pas in many circles to get these things mixed up.

printf tricks

A sometimes handy trick is registering your own specifiers for printf. Although it is GNU only, it is still a handy trick. For example, glibc comes with a built-in extension to reduce a float for you, e.g.

$ cat test.c
#include
#include

int main(void)
{
        float bytes = 1024.0 * 1024.0 * 1024.0;

        register_printf_function ('b', printf_size, printf_size_info);

        printf("%.0bB\n", bytes);

}

$ gcc -o test test.c
$ ./test
1gB

If we register the handler with a capital letter, the base will be 1000, otherwise it is 1024. It is unfortunate that this doesn't appear to support SI units, but it is still handy (and may one day).

The glibc manual has full details on customising printf further. You can easily add a specifier for your own structures, etc.

ERTOS Systems Administrator/Web Guru Job

If you're into Linux check out the ERTOS sysadmin job. Based at UNSW I can guarantee that if you're into Linux and open source you will not find a better work environment. I'm still constantly surprised by the bright people around me, and there's heaps of unique opportunities here. This isn't supporting Windows users ... you'll be hacking together solutions to problems you didn't know existed!

Feel free to contact me if you have any questions.

Google Maps - Lars Rasmussen

Google Maps has a new branch in Darling Harbour, Sydney.

GM has about 500 million squares. Within weeks they hope to have high res of Sydney; once or twice more zoom over the current even. Zachary's pizza, Berkeley is Lars' tip for best Pizza in the world.

Originally was a C++ browser with Where2. They took it to Google who said they liked the web. They tried plugins/flash. With AJAX it was 4 weeks for prototype; as opposed to years for C++ app. Javascript has inconsistencies across browsers, but compared to C++ across OS it is easy. However, compared to Google Earth it is nowhere near as rich in features.

On the first day hit top load they predicted for 6 months. Satellite imaging caused traffic to grow exponentially in one day.

Their design philosophy

  • end-user is royalty
  • go beyond browser lowest-common-denominator : because they're Google browser people ask them what they want.
  • dynamic as native apps; simple as google website
  • do what it takes
  • launch early + often : with 4 months dev, figure 3 months to launch. Larry Page played with demo and said "just go". labs.google.com playground.

Now under Google local. They listen to feedback; they discussed frame layout endlessly for the merge; launched and people hated it, changed in a week.

DHTML overview : javascript handles events, which can modify the structure of the page. browser does nothing until event handler finishes, then re-renders the page off screen and displays it.

Back-button was solved with virtual pages. What does it mean to push the back button? Back should go back through search history, not random junk. Normally the browser remembers clicks in iframes as history. To get around this you use a "virtual" iframe which copies what is happening into a main page. this way you can control the pages that are shown and consequently the back bar.

Backend map providers got nervous about hacks. housingmaps.com had inadvertently broken a link back to the data provider so Google weren't weren't being charged. API wasn't a technical challenge, but a contractual one. Google loved the maps hacks but couldn't tell anyone.

Google would never build a bus map route for Seattle, but they are inspired by it. Google watches the api key to see what cool things are happening. Katrina imagery request was from the government. If you want a job write a good maps hack.

Growth areas is mostly about integrating all of Googles data sources.

Cafe in san-fran? 20 billion web pages that probably can tell you where "Silicon Valley" is, but Google maps doesn't know. User logs show that users think google knows that stuff -- maybe comes from the "know it all" search box which people think can do anything.

Processing : different data and different confidence; can you join it? API is telatlis(?) data ... navtech data is on the site. can they merge? can they merge other databases? adding spacial infrastructure to google search infrastructure. "cheapest gas via work"

Routing: shortest path algorithm with weighted graph. If underlying data is correct works fine. However traffic not taken into account; taxi driver berated him that maps takes him straight into traffic -- some drivers love it because tourists bring it and they make money. Ridefinder -- nobody uses it but you get gps feeds from taxi drivers which you can use to re-weight your maps.

Drawing: pretty anti-aliased is nice. Others have had to catch up

Questions: When is Sydney maps coming? They are very tied to data format ... and Sydney maps don't' come in that format. working to remove that, early next year. where can you get Aussie map data? sensis, map data sciences. all comes from the government in some form.

Different grid systems around the world? No cartographers on the team really. In japan, use different datum; lat/long use different views of ovalness of the earth. datums need to be normalised.

How make money? Haven't settled on that just yet. Google figures users == money; not sure how but something will come up. Advertisements is obvious for "local search" stuff.

API via c++? no, web is focus.

Satellite tile age? higher == fresher. Google Earth looks after imagery; new Sydney about 6 months. Imagery has a "novelty" feel; big hits when first released but users don't always come back.

Censorship of photos: have to respect.

Obfuscated javascript? makes code smaller. time to download is now in seconds, needs to be milliseconds.

Government fudging? Not that he knows about ...

Projections of tiles? New different ones on the way, maybe. Night satellite, etc. winter/summer etc. geological features, historical maps. Have the data any maybe one day.

Disputed borders: no easy answers. Chinese laws say you can't move maps out of the country; but they do have maps on servers in China. Might have to present different things to different users.

Printing files side-by-side

I'm really not sure if there is an eaiser way to do this, but here is my newly most-used utility. It puts two files beside each other; but as opposed to sdiff/diff -y doesn't analyse it, and as opposed to paste keeps to fixed widths. Here is a screen shot.

$  python ./side-by-side.py --width 40 /tmp/afile.txt /tmp/afile2.txt
this is a file                           |  i have here another file
that has lines of                        |  that also has some text and
text.  to read                           |  some lines.  although it
                                         |  is slightly longer than the other
this is a really really really really re *  file with all these words
                                         |  in
                                         |  it

I'd love to hear from all the Python freaks how you could get the LOC even lower; every time I do something like this I find out about a new, quicker, way to recurse a list :)

#!/usr/bin/python

import sys, os
from optparse import OptionParser

class InFile:
    def __init__(self, filename):
        try:
            self.lines=[]
            self.maxlen = 0
            for l in open(filename).readlines():
                self.lines.append(l.rstrip())
        except IOError, (error, message):
            print "Can't read input %s : %s" % (filename, message)
            sys.exit(1)

        self.nlines = len(self.lines)
        if self.nlines == 0:
            self.lines.append("")
        self.maxlen = max(map(len, self.lines))

    # pad to the max len, with a extra space then the deliminator
    def pad_lines(self, nlines, width=0, nodiv=False, notrunc=False):
        if width == 0:
            width = self.maxlen
        pad = []
        for i in range(0, width):
            pad += " "
        # add on some extra for the divider and spaces
        pad += "   "
        padlen = len(pad)
        for i in range(0, nlines):
            try:
                linelen = len(self.lines[i])
            except IndexError:
                self.lines.append("")
                linelen = 0
            if (linelen > width):
                linelen = width
                if not notrunc:
                    pad[-2] = "*"
            elif nodiv:
                pad[-2] = " "
            else:
                pad[-2] = "|"
            self.lines[i] = self.lines[i][:linelen] +  "".join(pad[linelen - padlen:])

usage= "side-by-side [-w width] file1 file2 ... filen"

parser = OptionParser(usage, version=".1")

parser.add_option("-w", "--width", dest="width", action="store", type="int",
                      help="Set fixed width for each file", default=0)
parser.add_option("--last-div", dest="lastdiv", action="store_true",
                  help="Print divider after last file", default=False)
parser.add_option("--no-div", dest="nodiv", action="store_true",
                  help="Don't print any divider characters", default=False)
parser.add_option("--no-trunc", dest="notrunc", action="store_true",
                  help="Don't show truncation with a '*'", default=False)

(options, args) = parser.parse_args()

flist = []
if (len(args) == 0):
    print usage
    sys.exit(1)
for f in args:
    flist.append(InFile(f))

max_lines = max(map(lambda f: f.nlines, flist))

for i in range(0,len(flist)):
    if (len(flist)-1 == i):
        options.nodiv = not options.lastdiv
    flist[i].pad_lines(max_lines, options.width, options.nodiv, options.notrunc)

for l in range(0, max_lines):
    for f in flist:
        print f.lines[l],
    print

update: Leon suggests

pr -Tm file1 file2

Which is pretty close, but doesn't seem to put any divider between the files. Still might be a handy tool for your toolbox. It seems, from the pr man page, the word for doing this sort of thing is columnate.

update 2: Told you I'd learn new ways to iterate! Stephen Thorne came up with a neat solution. Some extra tricks he used

  • ' ' * n gives you n spaces. Pretty obvious when you think of it!
  • Usage of generators with the yield statement.
  • The itertools package with it's modified zip like izip function.

All very handy tips for your Python toolbox. If you're learning Python I'd reccommend solving this problem as you can really put to use some of the niceties of the language.

Speeding up a very slow apt-get update

If your laptop hard disk is so slow it sometimes seems like it would be faster to write the ones and zeros with a pencil on paper, a large apt-get update can take hours.

I have found most of the time is actually spent in /sbin/ldconfig. ldconfig is extremely disk intensive and a real bottle neck, and really only needs to be run once once you've finished, rather than every time a library is installed. Replacing ldconfig with a stub makes update times much more tolerable.

#!/bin/bash

exit 0

I think there has been some talk of apt being smart enough to only do it once for you, but current versions don't seem to implement this. When you're finished restore the old (real) binary and run it to update your shared library cache.

Upgrading the Opentel ODT4200 Firmware via Linux

The short version is : pull power from the box, setup Z-modem to send the new firmware file, turn the box on. That's it -- no Tridge style hacking required.

For those of you who missed the BBS experience and skipped straight to the world of TCP/IP, in more detail I did

$ wget http://www.kristal.com.au/opentel/downloads/ODT4200PVR_A_KRISTAL_V202.zip
$ unzip ODT4200PVR_A_KRISTAL_V202.zip
$ stty -F /dev/ttyS0 ispeed 115200 ospeed 115200
$ sz -vvv -b ODT4200PVR_A_KRISTAL_V202.pgm < /dev/ttyS0 > /dev/ttyS0
[walk over to box, pull out of power point and put back in, it downloads software and reboots]

sz is in the lrzsz package in Debian.

Running strings over the firmware doesn't uncover anything extraordinarily interesting; anyone else pulled it apart any further?

Update : the latest release appears to be a semi-public 2.04 version which can be downloaded directly from Oznetics.