Optimizing sketch size with 64 bit integer math

[RyanPierce]

My project to turn a Solder:Time II into a "Mars Watch" that tells local time, UTC time, and Mars times (MTC and mission clocks for Curiosity and Opportunity) so far has been successful. But I'm rapidly running out of sketch space, and I still have code I need to add. I'm using the Arduino 1.0.1 environment.

The big offender is the function I wrote that does the Mars time conversions. My code compiles to 29,872 bytes. Commenting out the function reduces it to 16,020 bytes. So this one function is taking up more than 13k! Here's the code:

Code: Select all

uint8_t MinTens = 0;
uint8_t MinOnes = 0;
uint8_t SecTens = 0;
uint8_t SecOnes = 0;
uint8_t UTCHourTens = 1;
uint8_t UTCHourOnes = 2;
uint8_t UTCDateOnes = 1;
uint8_t UTCDateTens = 0;
uint8_t UTCMonthCode = 1;
uint8_t UTCYearsOnes = 2;
uint8_t UTCYearsTens = 1;
uint32_t MarsSol;
uint8_t MarsHourTens;
uint8_t MarsHourOnes;
uint8_t MarsMinTens;
uint8_t MarsMinOnes;

void ComputeMarsTime(uint64_t Epoch)
{
  uint32_t UnixTime = 946684800L; // Jan 1, 2000 epoch
  uint8_t temp;
  temp = UTCYearsTens * 10 + UTCYearsOnes;
  UnixTime += (uint64_t) temp * 31536000L;
  // Now correct for leap years. Note that 2000 is NOT a leap year. We are also subtracting 1 because the current
  // year is dealt with later.
  if ( temp )
    temp--;
  temp -= temp % 4;
  temp /= 4;
  UnixTime += (uint64_t) temp * 86400L;
  if ( UTCMonthCode >= 2 )
    UnixTime += 2678400L;
  if ( UTCMonthCode >= 3 ) {
    UnixTime += 2419200L;
    // Correct for current year leap year. Note that 2000 is NOT a leap year.
    if ( UTCYearsTens != 0 || UTCYearsOnes != 0 ) {
      if ( ( UTCYearsTens * 10 + UTCYearsOnes ) % 4 == 0 )
        UnixTime += 86400L;
    }
  }
  if ( UTCMonthCode >= 4 )
    UnixTime += 2678400L;
  if ( UTCMonthCode >= 5 )
    UnixTime += 2592000L; 
  if ( UTCMonthCode >= 6 )
    UnixTime += 2678400L;
  if ( UTCMonthCode >= 7 )
    UnixTime += 2592000L;   
  if ( UTCMonthCode >= 8 )
    UnixTime += 2678400L;   
  if ( UTCMonthCode >= 9 )
    UnixTime += 2678400L;
  if ( UTCMonthCode >= 10 )
    UnixTime += 2592000L; 
  if ( UTCMonthCode >= 11 )
    UnixTime += 2678400L;
  if ( UTCMonthCode >= 12 )
    UnixTime += 2592000L;
  temp = UTCDateTens * 10 + UTCDateOnes;   
  UnixTime += (uint32_t) temp * 86400L;
  temp = UTCHourTens * 10 + UTCHourOnes;
  UnixTime += (uint32_t) temp * 3600L;
  temp = MinTens * 10 + MinOnes;
  UnixTime += (uint32_t) temp * 60L;
  temp = SecTens * 10 + SecOnes;
  UnixTime += (uint32_t) temp;

  uint64_t JD_UTC = (UnixTime * 1000000LL) / 86400LL + 2440587500000LL;
  uint64_t JD_TT = JD_UTC;
  JD_TT += (35LL*1000000LL + 32184000LL) / 86400LL; // Note: TODO MAKE LEAP SECOND CONFIGURABLE
  uint64_t J2000 = JD_TT - 2451545000000LL;
  uint64_t MSD = J2000 - 4500000LL;
  MSD *= 1000000000LL;
  MSD /= 1027491252LL;
  MSD += 44796000000LL - 960LL;
  MSD -= Epoch; // Epoch of mission time. Set to 0 for MTC.
  MarsSol = MSD / 1000000LL;
  uint32_t MTC = (MSD % 1000000L) * 24L; 
  byte MTCHour = MTC / 1000000L;
  MarsHourOnes = MTCHour % 10;
  MarsHourTens = ( MTCHour - MarsHourOnes ) / 10;
  MTC = (MTC % 1000000L)*60L;
  byte MTCMin = (MTC - MTC % 1000000L)/1000000L;
  MarsMinOnes = MTCMin % 10;
  MarsMinTens = ( MTCMin - MarsMinOnes ) / 10;
}


Now I can clearly reduce the number of math operations. But given that this is taking up >13k, I'd like to get a conceptual idea of where the big offender(s) are before I dive in, so I can optimize where it counts. My fear is that the bulk of this is assembly language subroutines for the 64 bit math operations I'm using. In that case, each addition, subtraction, modulo, etc. is only taking up a handful of bytes, and reducing the number of operations won't buy me much. I don't really need to squeeze every last byte out of it; I'm just looking for low hanging fruit that could free up a little breathing room.

Any advice would be appreciated. Thanks!

[Corwyn]

Why use 64 bit math at all?

If you don't need finer than 1 second resolution 32 bits gets you 136 years (earth years).

You also manage to use all the basic arithmetic operations, if you can figure out a way to avoid one that will avoid a bunch of imported code.

[RyanPierce]

Thanks! That's an interesting point.... If I could keep everything in terms of seconds, then yes, 32 bit would suffice.

But I think I still have a problem. The conversion requires that I take Earth time (in whatever time units since an epoch that I want) and either divide it by 1.027491252 or multiply it by the reciprocal of that number. Everything else is just linear shifting, either before or after the multiplication, which could just as easily be done in 32 bit as 64 bit. But that multiplication or division is the single critical point where any loss of precision can greatly shift the result. The only way I know to do that operation without resorting to 64 bit floating point is:

1. Cast the 32 bit time in Earth seconds to a uint64_t.
2. Either:
A. Multiply the 64 bit number by 1000000000 and divide it by 1027491252, or
B. Multiply the 64 bit number by 973244295, and divide it by 1000000000
3. Cast the resulting 64 bit time in Mars seconds back to a 32 bit number.

If I do all linear shifting to the 32 bit numbers, and not the 64 bit number, I can avoid including code for 64 bit addition and subtraction. But I imagine that code is small and trivial compared to 64 bit multiplication and division code, both of which I must use either way.

[Corwyn]

The problem being that the multiplication is overflowing, yes? So split the multiplication up by digit. Instead of multiplying by 973244295 and divviding by 100000000, (multiply by 9, divided by 10) + (multiply by 7, divided by 100) + ...

Now instead of doing it in base ten, do it in base 16, so those divisions are actually just shifts.

Does that work?

[westfw]

You're using an algorithm designed for quick and short implementation on a modern computer with floating point hardware. You should be able to find or design an integer-based algorithm, even if that means taking your "epoch" and subtracting off one martian year worth of seconds, one at a time...

[odometer]

I have suggested this for another, similar, project, and now I will suggest it for this.

The way you multiply by a funny decimal number is you don't. You use an approximation.
For example, let's say I want to multiply something by pi. If I don't need too much accuracy, I will approximate pi as 3+1/7, and thus I will approximate (pi*x) as (3*x)+(x/7). If I need some more accuracy, I will approximate pi as 3+16/113, and thus approximate (pi*x) as (3*x)+((16*x)/113). If I am worried that multiplying x by 16 will give me too big a result, then I will compute ((16*x)/113) as ((x DIV 113)*16)+(((x MOD 113)*16)/113).

How did I get this magic fraction 16/113 to approximate (pi-3)? I used this page:
http://mysite.verizon.net/res148h4j/javascript/script_fractions.html
(Note: the answer will appear in a pop-up window.)

For your number, 1.027491252, I get 44925/43723 (which equals 1+1202/43723) as my approximation. So you need x+((1202*x)/43723). Since 1202*x will overflow a 32-bit integer, you can calculate (1202*x)/43723 as ((x DIV 43723)*1202) + ((x MOD 43723)*1202)/43723.

*** OOPS!! *** I just realized: that's if you want to multiply by that number. But you don't want to do that. You want to divide by it. Please see my next post in this thread.

Your expression in code (for the multiplication) will be

Code: Select all

x + ((x / 43723) * 1202) + (((x % 43723) * 1202) / 43723)

And you won't need a single 64-bit integer. Anywhere.

[adafruit_support_bill]

Thanks odometer. That is a very useful technique.

[RyanPierce]

Thanks so much! I hadn't thought of that.... I'm going to give it a try.

[odometer]

I'm sorry, I just realized something.

For your number, 1.027491252, I get 44925/43723 (which equals 1+1202/43723) as my approximation.

I made a mistake. You did not want to multiply by this number; you wanted to divide by it.

Easy fix. Flip the fraction over to get 43723/44925. We could deal with it as such, or we could treat it as 1-1202/44925 (that's supposed to be a minus sign, not a hyphen).

First method: use the fraction 43723/44925 as such.
Just to make sure we won't overflow anything, let's multiply 43723 by 44925. The result fits -- just barely -- inside 31 bits, so we're safe using a 32-bit integer. So our code is:

Code: Select all

((x / 44925) * 43723) + (((x % 44925) * 43723) / 44925)


Second method: use 1 minus 1202/44925. The code will be:

Code: Select all

x - ((x / 44925) * 1202) - (((x % 44925) * 1202) / 44925)


The two methods will sometimes produce slightly different (by 1 unit) results, having to do with the way the rounding works. Because 1 unit represents 1 second, the difference isn't important, assuming that you only care about the time to within a second or so.

[odometer]

I noticed a couple of other things.

First (and most serious), your code handles leap years incorrectly. Really, the year 2000 was a leap year. Please see here: http://www.tondering.dk/claus/cal/gregorian.php#leapyears
The simple 4-year rule for leap years works from 1901 to 2099, but fails outside this range.

Second, it seems that you are trying to calculate in millionths of a Mars day. Since you want your output in units of 1, 1/24, 1/1440, and 1/86400 of a Mars day, the number 1000000 or 1/1000000 doesn't seem to have any useful place here. I have told you how to convert the count of Earth seconds to the count of Mars seconds. Now I will show you how to transform raw seconds (Earth or Mars, the method is the same) to days, hours, minutes, and seconds.

Method A:
Divide raw_seconds by 86400. Quotient is days, remainder is used as follows:
Divide by 3600. Quotient is hours, remainder is used as follows:
Divide by 60. Quotient is minutes, remainder is seconds.

Method B:
Divide raw_seconds by 60. Remainder is seconds, quotient is used as follows:
Divide by 60. Remainder is minutes, quotient is used as follows:
Divide by 24. Remainder is hours; quotient is days.

Because divisions on the Arduino are a bit slow, I did some Googling and some tweaking to find other methods. Some are very weird indeed. For instance, to divide a 32-bit integer by 60, giving quotient and remainder, you begin by calculating the remainder, as follows:

Code: Select all

byte mod60( unsigned int n ) {
    byte r = ((byte)(n)) & 3;
    unsigned int t = (n >> 2);
    t = (t >> 16) + (t & 0xFFFF);
    t = (t >>  8) + (t & 0xFF);
    t = (t >>  4) + (t & 0xF);
    r += (4 * (byte)t);
    while (r >= 60) r -= 60;
    return r;
}


Then, the quotient is:

Code: Select all

((n - (unsigned int)(mod60(n))) >> 2) * 0xEEEEEEEF


[RyanPierce]

@odometer you got me on the 2000 leap year issue. I was aware of the century leap year rule, but for some reason I thought, incorrectly, that 2000 wasn't divisible by 400.... Ooops.

That said, I'm completely perplexed because the code that I originally posted actually works and gives the right result. It agrees with NASA's Mars24 Sunclock program, as well as the illustrated examples on James Tauber's github page. Looking at my code, it appears I'm failing to add an extra day's worth of seconds for the year 2000 for any year after 2000. (I'm also not adding a day if the current year is 2000 and the month is March or later, but that bug shouldn't really matter.) So I'd expect my Mars time to be on the order of 40 minutes off.

I'm going to have to dig into this to do more debugging. I'd like to compare my constructed Unix time with actual Unix time, and verify that I'm not compensating for being off by a day elsewhere.

The full code is here: https://github.com/RyanPierce/Mars_Watch Note that I haven't made the changes you suggested yet to keep everything 32 bit; once it worked, I figured I'd stop and optimize it after the holidays; even with the 64 bit math bloating the code, I can just barely fit in most of the functionality I need.

In other news, my rover-hugging girlfriend absolutely loves the watch! It's been quite a conversation starter.

Thanks!

[odometer]

You're making the same error I did at least once: you forgot that day 0 of your count corresponds to day 1 of the month!

But then again, it seems that the main functionality (i.e. rover-hugging girlfriend + Mars watch = Ryan-hugging girlfriend) is A-OK, and that's what matters most.

Next stage: implement Darian calendar.