Let's start with some generic theory behind LNS....
Unlike floating point, which is now largely defined by the IEEE 754 standard, LNS (Logarithmic Number System) representations are not really well standardized. There are a variety of variations, none of which has really become a universally-accepted standard... which is odd, because LNS can deliver better accuracy and higher speed than similar-precision floating point. Generally, everybody agrees to use a base of 2, but there's very little else agreed upon. Let's initially ignore the notation issues and talk about algorithms:
Here are some generic references to LNS:
Bit Positions | ||||||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
Value Type | 15 | 14 | 13 | 12 | 11 | 10 | 9 | 8 | 7 | 6 | 5 | 4 | 3 | 2 | 1 | 0 |
Normal LNS Value | Sign | Fixed-Point Positive Log representing value in [2.95469e-39 .. 3.36617e+38] | ||||||||||||||
Zero | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
Not-A-Number | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
+Infinity | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 |
-Infinity | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 |
Because 0x0000 and 0x7fff are treated specially (as 0 and infinity), the minimum representable value corresponds to 0x0001 and the maximum is 0x7ffe. These values are magnitudes only, made signed using the separate sign bit. Because the magnitude range exceeds that of a 32-bit IEEE 754 float, we cannot use float computations to encode or decode magnitudes; however, it is well within the range and accuracy of double. Thus, the encoding of the magnitudes is sepcifed by two functions mapping between double and our LNS magnitudes. The first converts a mag into the corresponding double:
double mag2double(mag m) { return(pow(2.0, ((double) m) / 128.0) / pow(2.0, 128.0)); }
Note that pow(2.0, 128.0) is already problematic for a float representation. The second performs the opposite transformation:
mag double2mag(double d) { return((mag) (log2(d * pow(2.0, 128.0)) * 128.0)); }
However, there is a slight issue here involving rounding -- these routines don't round. This means that values can be as much as one least-significant-bit off from the best value.
LNS encoding is strange enough, and we're using a slightly stranger encoding still, so it's useful to be able to play with it. For that, I've created this 16-bit LNS Conversions CGI. It's a simple WWW form that you can use to convert between decimal values and LNS, as well as to perform LNS negate, multiply, divide, add/subtract, and compare.
Ok, so you get the basic idea... the catch is that you need algorithms. Here, algorithms are presented for each of the LNS operations in the Logick instruction set. For all of these algorithms, a few definitions are needed:
#define SIGN 0x8000 #define ONE 0x4000 #define ZERO 0x0000 #define LNSNAN 0x8000 #define POSINF 0x7fff #define NEGINF 0xffff #define INFINITE(V) (((V) & POSINF) == POSINF)
Also note that lns refers to a complete LNS value as an integer, while mag refers simply to the unsigned LNS magnitude.
Aside from a variety of special cases, the key to multiplication is addition of signed log values... which, unfortunately, our aren't. The encoding we use essentially adds a bias of ONE, so that bias needs to be subtracted before the add, and reapplied after.
lns lnsmul(lns a, lns b) { if ((a == LNSNAN) || (b == LNSNAN)) return(LNSNAN); if ((a == ZERO) || (b == ZERO)) return(ZERO); if (INFINITE(a) || INFINITE(b)) { return(POSINF | ((a ^ b) & SIGN)); } lns s = (a ^ b) & SIGN; mag m = ((a & POSINF) -ONE) + ((b & POSINF) - ONE); m += ONE; if (m >= POSINF) m = POSINF; return(s | m); }
Much like multiply, divide becomes subtract after checking various special cases. The same ONE bias removal and reapplication is needed.
lns lnsdiv(lns a, lns b) { if ((a == LNSNAN) || (b == LNSNAN)) return(LNSNAN); if (b == ZERO) return(LNSNAN); if (a == ZERO) return(ZERO); if (INFINITE(b)) { if (INFINITE(a)) return(LNSNAN); return(ZERO); } if (INFINITE(a)) return(a); lns s = (a ^ b) & SIGN; mag m = ((a & POSINF) -ONE) - ((b & POSINF) - ONE); m += ONE; if (m >= POSINF) m = POSINF; return(s | m); }
Comparisons of LNS value are fairly straightforward. Although bit positions do not have simple weightings, incrementing the unsigned LNS magnitude as a binary number still produces a monotonically increasing sequence of values. Thus, comparisons are essentially identical to how sign+magnitude integer values would be compared. They are not quite the same, however, because we have four special values that require some special handling.
Rather than testing a specific condition, the compare logic here returns a 2-bit number specifying what the relationship between the values compared is:
#define ISEQ 0 // equal #define ISLT 1 // less than #define ISGT 2 // greater than #define ISUNDEF 3 // relationship not defined
That seems simple enough, but what's that ISUNDEF thing? Well, some comparisons have undefined results. It's fairly obvious that nan isn't related to any value, including nan. There is some ambiguity about relationships involving inf; this suggests that inf is treated much like any other value, but it really is a little dicey to be saying inf==inf, because they might really be representing entirely different numbers that both happen to be larger than is directly representable. In any case, the IEEE 754 standard suggests inf==inf, so we'll stick with that.
int lnsrelate(lns a, lns b) { if ((a == LNSNAN) || (b == LNSNAN)) return(ISUNDEF); if ((a & SIGN) != (b & SIGN)) return((a & SIGN) ? ISLT : ISGT); if (a == b) return(ISEQ); /* Unsigned compare magnitude, but flip result if sign was negative */ mag ma = (a & POSINF); mag mb = (b & POSINF); if (a & SIGN) { return((ma < mb) ? ISGT : ISLT); } else { return((ma < mb) ? ISLT : ISGT); } }
Negation is simply a matter of flipping the sign bit -- except if the value is 0 or nan.
lns lnsneg(lns a) { return((a & POSINF) ? (a ^ SIGN) : a); }
As always for LNS, addition is a complicated thing. In fact, it's a bit more complicated by the fact that it is adding signed values, so an addition of values with differing signs is a subtract. The complex non-linear mappings described above for add and subtract need to be implemented by table lookup from the appropriate one of two tables. The catch is that we don't want huge tables, and the obvious structure would be two tables that each contain as many entries as the range of b-a -- because either a or b could be larger, the 15-bit magnitude subtraction has a range of 16 bits, so there would be 65,536 entries in each table. In other words, the two tables would occupy a total of 262,144 bytes... which is too much.
The first thing we do to simplify the tables is to impose the constraint that b-a is non-negative, which requires sorting the a and b values, but literally halves the size of the tables. After that, the interesting thing is that the values in each of the two tables look pretty boring -- as in monotonous. In fact, the values in each of the tables are monotonic sequences, so various compression methods could easily be applied. However, even more simply, we can use the fact that the majority of the table entries have an obvious "default" value. The add table entries from 0x03c4 to 0x7ffe are each simply a copy of their index value; this means there are only 964 interesting entries, and a table of 964 entries (or 1,024, to keep indexing easy) can suffice. The subtract table has an even more obvious structure, with all entries from 0x1850 to 0x0x7ffe having the value 0. Of the 6,224 remaining entries, those from 0x03c4 to 0x1850 all have the value 0x0001. Thus, with a little extra logic, the total space needed is no more than 2*964*2=3,856 bytes (or 4,096 bytes if you make the tables 1,024) -- a very reasonable lookup table size.
The following C code works with the 964-element tables:
lns lnsadd(lns a, lns b) { if ((a == LNSNAN) || (b == LNSNAN)) return(LNSNAN); if (INFINITE(a) && INFINITE(b)) return(LNSNAN); if (INFINITE(a)) return(a); if (INFINITE(b)) return(b); if (a == ZERO) return(b); if (b == ZERO) return(a); lns sa = (a & SIGN); lns sb = (b & SIGN); mag ma = a & POSINF; mag mb = b & POSINF; mag m; if (ma > mb) { ma ^= mb; mb ^= ma; ma ^= mb; sa ^= sb; sb ^= sa; sa ^= sb; } m = mb - ma; /* Same sign? */ if (sa == sb) { /* Add */ if (m < 0x03c4) m = addtab[m]; return((ma + m) | sa); } else { /* Subtract */ if (m > 0x1850) return(mb | sb); if (m > 0x03c4) return((mb - 1) | sb); return((mb - subtab[m]) | sb); } }
What do the tables look like? Well, here they are in VMEM format. First, addtab:
//addtab[0x0000..0x03c3] @0000 0080 0080 0081 0081 0082 0082 0083 0083 0084 0084 0085 0085 0086 0086 0087 0087 0088 0088 0089 0089 008a 008a 008b 008b 008c 008c 008d 008d 008e 008f 008f 0090 0090 0091 0091 0092 0092 0093 0093 0094 0095 0095 0096 0096 0097 0097 0098 0098 0099 009a 009a 009b 009b 009c 009c 009d 009e 009e 009f 009f 00a0 00a1 00a1 00a2 00a2 00a3 00a3 00a4 00a5 00a5 00a6 00a6 00a7 00a8 00a8 00a9 00a9 00aa 00ab 00ab 00ac 00ac 00ad 00ae 00ae 00af 00af 00b0 00b1 00b1 00b2 00b3 00b3 00b4 00b4 00b5 00b6 00b6 00b7 00b8 00b8 00b9 00b9 00ba 00bb 00bb 00bc 00bd 00bd 00be 00bf 00bf 00c0 00c1 00c1 00c2 00c2 00c3 00c4 00c4 00c5 00c6 00c6 00c7 00c8 00c8 00c9 00ca 00ca 00cb 00cc 00cc 00cd 00ce 00ce 00cf 00d0 00d0 00d1 00d2 00d2 00d3 00d4 00d5 00d5 00d6 00d7 00d7 00d8 00d9 00d9 00da 00db 00db 00dc 00dd 00de 00de 00df 00e0 00e0 00e1 00e2 00e2 00e3 00e4 00e5 00e5 00e6 00e7 00e7 00e8 00e9 00ea 00ea 00eb 00ec 00ec 00ed 00ee 00ef 00ef 00f0 00f1 00f2 00f2 00f3 00f4 00f4 00f5 00f6 00f7 00f7 00f8 00f9 00fa 00fa 00fb 00fc 00fd 00fd 00fe 00ff 0100 0100 0101 0102 0103 0103 0104 0105 0106 0106 0107 0108 0109 0109 010a 010b 010c 010c 010d 010e 010f 0110 0110 0111 0112 0113 0113 0114 0115 0116 0117 0117 0118 0119 011a 011a 011b 011c 011d 011e 011e 011f 0120 0121 0122 0122 0123 0124 0125 0126 0126 0127 0128 0129 012a 012a 012b 012c 012d 012e 012e 012f 0130 0131 0132 0132 0133 0134 0135 0136 0136 0137 0138 0139 013a 013b 013b 013c 013d 013e 013f 013f 0140 0141 0142 0143 0144 0144 0145 0146 0147 0148 0149 0149 014a 014b 014c 014d 014e 014e 014f 0150 0151 0152 0153 0153 0154 0155 0156 0157 0158 0158 0159 015a 015b 015c 015d 015e 015e 015f 0160 0161 0162 0163 0164 0164 0165 0166 0167 0168 0169 016a 016a 016b 016c 016d 016e 016f 0170 0170 0171 0172 0173 0174 0175 0176 0176 0177 0178 0179 017a 017b 017c 017d 017d 017e 017f 0180 0181 0182 0183 0184 0184 0185 0186 0187 0188 0189 018a 018b 018c 018c 018d 018e 018f 0190 0191 0192 0193 0193 0194 0195 0196 0197 0198 0199 019a 019b 019b 019c 019d 019e 019f 01a0 01a1 01a2 01a3 01a4 01a4 01a5 01a6 01a7 01a8 01a9 01aa 01ab 01ac 01ad 01ad 01ae 01af 01b0 01b1 01b2 01b3 01b4 01b5 01b6 01b6 01b7 01b8 01b9 01ba 01bb 01bc 01bd 01be 01bf 01c0 01c0 01c1 01c2 01c3 01c4 01c5 01c6 01c7 01c8 01c9 01ca 01cb 01cb 01cc 01cd 01ce 01cf 01d0 01d1 01d2 01d3 01d4 01d5 01d6 01d7 01d7 01d8 01d9 01da 01db 01dc 01dd 01de 01df 01e0 01e1 01e2 01e3 01e3 01e4 01e5 01e6 01e7 01e8 01e9 01ea 01eb 01ec 01ed 01ee 01ef 01f0 01f0 01f1 01f2 01f3 01f4 01f5 01f6 01f7 01f8 01f9 01fa 01fb 01fc 01fd 01fe 01fe 01ff 0200 0201 0202 0203 0204 0205 0206 0207 0208 0209 020a 020b 020c 020d 020e 020e 020f 0210 0211 0212 0213 0214 0215 0216 0217 0218 0219 021a 021b 021c 021d 021e 021f 021f 0220 0221 0222 0223 0224 0225 0226 0227 0228 0229 022a 022b 022c 022d 022e 022f 0230 0231 0232 0232 0233 0234 0235 0236 0237 0238 0239 023a 023b 023c 023d 023e 023f 0240 0241 0242 0243 0244 0245 0246 0247 0247 0248 0249 024a 024b 024c 024d 024e 024f 0250 0251 0252 0253 0254 0255 0256 0257 0258 0259 025a 025b 025c 025d 025e 025f 025f 0260 0261 0262 0263 0264 0265 0266 0267 0268 0269 026a 026b 026c 026d 026e 026f 0270 0271 0272 0273 0274 0275 0276 0277 0278 0279 027a 027b 027b 027c 027d 027e 027f 0280 0281 0282 0283 0284 0285 0286 0287 0288 0289 028a 028b 028c 028d 028e 028f 0290 0291 0292 0293 0294 0295 0296 0297 0298 0299 029a 029b 029c 029c 029d 029e 029f 02a0 02a1 02a2 02a3 02a4 02a5 02a6 02a7 02a8 02a9 02aa 02ab 02ac 02ad 02ae 02af 02b0 02b1 02b2 02b3 02b4 02b5 02b6 02b7 02b8 02b9 02ba 02bb 02bc 02bd 02be 02bf 02c0 02c1 02c2 02c3 02c4 02c5 02c5 02c6 02c7 02c8 02c9 02ca 02cb 02cc 02cd 02ce 02cf 02d0 02d1 02d2 02d3 02d4 02d5 02d6 02d7 02d8 02d9 02da 02db 02dc 02dd 02de 02df 02e0 02e1 02e2 02e3 02e4 02e5 02e6 02e7 02e8 02e9 02ea 02eb 02ec 02ed 02ee 02ef 02f0 02f1 02f2 02f3 02f4 02f5 02f6 02f7 02f8 02f9 02fa 02fa 02fb 02fc 02fd 02fe 02ff 0300 0301 0302 0303 0304 0305 0306 0307 0308 0309 030a 030b 030c 030d 030e 030f 0310 0311 0312 0313 0314 0315 0316 0317 0318 0319 031a 031b 031c 031d 031e 031f 0320 0321 0322 0323 0324 0325 0326 0327 0328 0329 032a 032b 032c 032d 032e 032f 0330 0331 0332 0333 0334 0335 0336 0337 0338 0339 033a 033b 033c 033d 033e 033f 0340 0341 0342 0343 0344 0344 0345 0346 0347 0348 0349 034a 034b 034c 034d 034e 034f 0350 0351 0352 0353 0354 0355 0356 0357 0358 0359 035a 035b 035c 035d 035e 035f 0360 0361 0362 0363 0364 0365 0366 0367 0368 0369 036a 036b 036c 036d 036e 036f 0370 0371 0372 0373 0374 0375 0376 0377 0378 0379 037a 037b 037c 037d 037e 037f 0380 0381 0382 0383 0384 0385 0386 0387 0388 0389 038a 038b 038c 038d 038e 038f 0390 0391 0392 0393 0394 0395 0396 0397 0398 0399 039a 039b 039c 039d 039e 039f 03a0 03a1 03a2 03a3 03a4 03a5 03a6 03a7 03a8 03a9 03aa 03ab 03ac 03ad 03ae 03af 03b0 03b1 03b2 03b3 03b4 03b5 03b6 03b7 03b8 03b9 03ba 03bb 03bc 03bd 03be 03bf 03c0 03c1 03c2 03c3 03c4
Next, subtab:
//subtab[0x0000..0x03c3] @0000 0000 03c4 0344 02fa 02c5 029c 027b 025f 0247 0232 021f 020e 01fe 01f0 01e3 01d7 01cb 01c0 01b6 01ad 01a4 019b 0193 018c 0184 017d 0176 0170 016a 0164 015e 0158 0153 014e 0149 0144 013f 013b 0136 0132 012e 012a 0126 0122 011e 011a 0117 0113 0110 010c 0109 0106 0103 0100 00fd 00fa 00f7 00f4 00f2 00ef 00ec 00ea 00e7 00e5 00e2 00e0 00de 00db 00d9 00d7 00d5 00d2 00d0 00ce 00cc 00ca 00c8 00c6 00c4 00c2 00c1 00bf 00bd 00bb 00b9 00b8 00b6 00b4 00b3 00b1 00af 00ae 00ac 00ab 00a9 00a8 00a6 00a5 00a3 00a2 00a1 009f 009e 009c 009b 009a 0098 0097 0096 0095 0093 0092 0091 0090 008f 008d 008c 008b 008a 0089 0088 0087 0086 0085 0084 0083 0082 0081 0080 0080 007f 007e 007d 007c 007b 007a 0079 0078 0077 0076 0075 0074 0073 0073 0072 0071 0070 006f 006e 006e 006d 006c 006b 006a 006a 0069 0068 0067 0067 0066 0065 0065 0064 0063 0062 0062 0061 0060 0060 005f 005e 005e 005d 005c 005c 005b 005a 005a 0059 0059 0058 0057 0057 0056 0056 0055 0054 0054 0053 0053 0052 0052 0051 0051 0050 004f 004f 004e 004e 004d 004d 004c 004c 004b 004b 004a 004a 0049 0049 0048 0048 0047 0047 0047 0046 0046 0045 0045 0044 0044 0043 0043 0043 0042 0042 0041 0041 0040 0040 0040 003f 003f 003e 003e 003e 003d 003d 003c 003c 003c 003b 003b 003b 003a 003a 0039 0039 0039 0038 0038 0038 0037 0037 0037 0036 0036 0036 0035 0035 0035 0034 0034 0034 0033 0033 0033 0032 0032 0032 0031 0031 0031 0031 0030 0030 0030 002f 002f 002f 002f 002e 002e 002e 002d 002d 002d 002d 002c 002c 002c 002c 002b 002b 002b 002b 002a 002a 002a 002a 0029 0029 0029 0029 0028 0028 0028 0028 0027 0027 0027 0027 0026 0026 0026 0026 0026 0025 0025 0025 0025 0024 0024 0024 0024 0024 0023 0023 0023 0023 0023 0022 0022 0022 0022 0022 0021 0021 0021 0021 0021 0020 0020 0020 0020 0020 001f 001f 001f 001f 001f 001f 001e 001e 001e 001e 001e 001e 001d 001d 001d 001d 001d 001d 001c 001c 001c 001c 001c 001c 001b 001b 001b 001b 001b 001b 001a 001a 001a 001a 001a 001a 001a 0019 0019 0019 0019 0019 0019 0019 0018 0018 0018 0018 0018 0018 0018 0018 0017 0017 0017 0017 0017 0017 0017 0016 0016 0016 0016 0016 0016 0016 0016 0015 0015 0015 0015 0015 0015 0015 0015 0015 0014 0014 0014 0014 0014 0014 0014 0014 0014 0013 0013 0013 0013 0013 0013 0013 0013 0013 0012 0012 0012 0012 0012 0012 0012 0012 0012 0012 0011 0011 0011 0011 0011 0011 0011 0011 0011 0011 0011 0010 0010 0010 0010 0010 0010 0010 0010 0010 0010 0010 0010 000f 000f 000f 000f 000f 000f 000f 000f 000f 000f 000f 000f 000e 000e 000e 000e 000e 000e 000e 000e 000e 000e 000e 000e 000e 000d 000d 000d 000d 000d 000d 000d 000d 000d 000d 000d 000d 000d 000d 000c 000c 000c 000c 000c 000c 000c 000c 000c 000c 000c 000c 000c 000c 000c 000c 000b 000b 000b 000b 000b 000b 000b 000b 000b 000b 000b 000b 000b 000b 000b 000b 000b 000a 000a 000a 000a 000a 000a 000a 000a 000a 000a 000a 000a 000a 000a 000a 000a 000a 000a 000a 0009 0009 0009 0009 0009 0009 0009 0009 0009 0009 0009 0009 0009 0009 0009 0009 0009 0009 0009 0009 0009 0008 0008 0008 0008 0008 0008 0008 0008 0008 0008 0008 0008 0008 0008 0008 0008 0008 0008 0008 0008 0008 0008 0008 0008 0007 0007 0007 0007 0007 0007 0007 0007 0007 0007 0007 0007 0007 0007 0007 0007 0007 0007 0007 0007 0007 0007 0007 0007 0007 0007 0007 0007 0006 0006 0006 0006 0006 0006 0006 0006 0006 0006 0006 0006 0006 0006 0006 0006 0006 0006 0006 0006 0006 0006 0006 0006 0006 0006 0006 0006 0006 0006 0006 0006 0006 0005 0005 0005 0005 0005 0005 0005 0005 0005 0005 0005 0005 0005 0005 0005 0005 0005 0005 0005 0005 0005 0005 0005 0005 0005 0005 0005 0005 0005 0005 0005 0005 0005 0005 0005 0005 0005 0005 0005 0005 0005 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0004 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0003 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002 0002