Algorithm for calculating trigonometry, logarithms or something like that. ONLY addition-subtraction

Question

I am restoring the Ascota 170 antique mechanical programmable computer. It is already working. Now I’m looking for an algorithm to demonstrate its capabilities — like calculating trigonometric or logarithmic tables. Or something like that. Unfortunately, from mathematical operations, a computer is only capable of adding and subtracting integers (55 registers from -1E12 to 1E12). There is not even a shift-to-digit operation — so that it can be programmatically implemented to multiply only by very small numbers. But its logical operations are very well developed.

Could you advise me any suitable algorithm?

Answer 1

So what you're doing is really kinda awesome. And as it happens, I can explain quite a bit about how to implement fractional logarithms using only integer addition and subtraction! This post is going to be long, but there's lots of detail included, and a working implementation at the end, and it should be enough for you to do some fun things with your weird mechanical computer.

---

Implementing Comparisons

You're going to need to be able to compare numbers. While you said you can perform comparisons == 0 and > 0, that's not really quite enough for most of the interesting algorithms you'll want to implement. You need relative comparisons, which can be determined via subtraction:

isLessThan(a, b):
  diff = b - a
  if diff > 0 then return true
  else return false

isGreaterThan(a, b):
  diff = a - b
  if diff > 0 then return true
  else return false

isLessThanOrEqual(a, b):
  diff = a - b
  if diff > 0 then return false
  else return true

isGreaterThanOrEqual(a, b):
  diff = b - a
  if diff > 0 then return false
  else return true

For the rest of this post, I'm just going to write the simpler form of a > b, but if you can't do that directly, you can substitute in one of the operations above.

---

Implementing Shifts

Now, since you don't have digit-shifting hardware, you'll have to create "routines" to implement it. A left-shift is easy: Add a number to itself, and again, and again, and then add the original number, and then add it one more time; and that's the equivalent of shifting left by 1 digit.

So shift left by one digit, or multiply-by-ten:

shiftLeft(value):
    value2 = value + value
    value4 = value2 + value2
    value5 = value4 + value
    return value5 + value5

Shifting by many digits is just repeated invocation of shiftLeft():

shl(value, count):
  repeat:
    if count <= 0 then goto done
    value = shiftLeft(value)
    count = count - 1
  done:
    return value

Shifting right by one digit is a little harder: We need to do this with repeated subtraction and addition, as in the pseudocode below:

shr(value, count):
    if count == 0 then return value

    index = 11
    shifted = 0
  repeat1:
    if index < 0 then goto done
    adder = shl(1, index - count)
    subtractor = shl(adder, count)
  repeat2:
    if value <= subtractor then goto next
    value = value - subtractor
    shifted = shifted + adder
    goto repeat2
  next:
    index = index - 1
    goto repeat1

  done:
    return count

Conveniently, since it's hard to shift right in the first place, the algorithm lets us directly choose how many digits to shift by.

---

Multiplication

It looks like your hardware might have multiplication? But if it doesn't, you can implement multiplication using repeated addition and shifting. Binary multiplication is the easiest form to implement that's actually efficient, and that requires us to first implement multiplyByTwo() and divideByTwo(), using the same basic techniques that we used to implement shiftLeft() and shr().

Once you have those implemented, multiplication involves repeatedly slicing off the last bit of one of the numbers, and if that bit is a 1, then adding a growing version of the other number to the running total:

multiply(a, b):
    product = 0
  repeat:
    if b <= 0 then goto done
    nextB = divideByTwo(b)
    bit = b - multiplyByTwo(nextB)
    if bit == 0 then goto skip
    product = product + a
  skip:
    a = a + a
    b = nextB
    goto repeat
  done:
    return product

A full implementation of this is included below, if you need it.

---

Integer Logarithms

We can use our ability to shift right by a digit to calculate the integer part of the base-10 logarithm of a number — this is really just how many times you can shift the number right before you reach a number too small to shift.

integerLogarithm(value):

    count = 0
  repeat:
    if value <= 9 then goto done
    value = shiftRight(value)
    count = count + 1
    goto repeat
  done:
    return count

So for 0-9, this returns 0; for 10-99, this returns 1; for 100-999 this returns 2, and so on.

---

Integer Exponents

The opposite of the above algorithm is pretty trivial: To calculate 10 raised to an integer power, we just shift the digits left by the power.

integerExponent(count):

    value = shl(1, count)
    return value

So for 0, this returns 1; for 1, this return 10; for 2, this returns 100; for 3, this returns 1000; and so on.

---

Splitting the Integer and Fraction

Now that we can handle integer powers and logarithms, we're almost ready to handle the fractional part. But before we can really talk about how to compute the fractional part of the logarithm, we have to talk about how to divide up the problem so we can compute the fractional part separately from the integer part. Ideally, we only want to deal with computing logarithms for numbers in a fixed range — say, from 1 to 10, rather than from 1 to infinity.

We can use our integer logarithm and exponent routines to slice up the full logarithm problem so that we're always dealing with a value in the range of [1, 10), no matter what the input number was.

First, we calculate the integer logarithm, and then the integer exponent, and then we subtract that from the original number. Whatever is left over is the fractional part that we need to calculate: And then the only remaining exercise is to shift that fractional part so that it's always in a consistent range.

normalize(value):

    intLog = integerLogarithm(value)    // From 0 to 12 (meaningful digits)
    if intLog <= 5 then goto lessThan
    value = shr(value, intLog - 5)
    goto done
  lessThan:
    value = shl(value, 5 - intLog)
  done:
    return value

You can convince yourself with relatively little effort that no matter what the original value was, its highest nonzero digit will be moved to column 7: So "12345" will become "000000123450" (i.e., "0000001.23450"). This allows us to pretend that there's always an invisible decimal point a little more than halfway down the number, so that now we only need to solve the problem of calculating logarithms of values in the range of [1, 10).

(Why "more than halfway"? We will need the upper half of the value to always be zero, and you'll see why in a moment.)

---

Fractional Logarithms

Knuth explains how to do this in The Art of Computer Programming, section 1.2.2. Our goal will be to calculate log10(x) so that for some values of b1, b2, b3 ... , where n is already 0 (because we split out the integer portion above):

log10(x) = n + b1/2 + b2/4 + b3/8 + b4/16 + ...

Knuth says that we can obtain b1, b2, b3 ... like this:

To obtain b1, b2, ..., we now set x0 = x / 10^n and, for k >= 1,

b[k] = 0, x[k] = x[k-1] ^ 2, if x[k-1] ^ 2 < 10;

b[k] = 1, x[k] = x[k-1] ^ 2 / 10, if x[k-1] ^ 2 >= 10.

That is to say, each step uses pseudocode loop something like this:

fractionalLogarithm(x):
  for i = 1 to numberOfBinaryDigitsOfPrecision:
    nextX = x * x
    if nextX < 10 then:
      b[i] = 0
    else:
      b[i] = 1
      nextX = nextX / 10

In order for this to work using the fixed-point numbers we have above, we have to implement x * x using a shift to move the decimal point back into place, which will lose some digits. This will cause error to propagate, as Knuth says, but it will give enough accuracy that it's good enough for demonstration purposes.

So given a fractional value generated by normalize(value), we can compute its fractional binary logarithm like this:

fractionalLogarithm(value):
  for i = 1 to 20:
    value = shr(value * value, 6)
    if value < 1000000 then:
      b[i] = 0
    else:
      b[i] = 1
      value = shr(value, 1)

But a binary fractional logarithm — individual bits! — isn't especially useful, especially since we computed an decimal version of the integer part of the logarithm in the earlier step. So we'll modify this one more time, to calculate a decimal fractional logarithm, to five places, instead of calculating an array of bits; for that, we'll need a table of 20 values that represent the conversions of each of those bits to decimal, and we'll store them as fixed-point as well:

table[1] = 1/(2^1) = 1/2  = 500000
table[2] = 1/(2^2) = 1/4  = 250000
table[3] = 1/(2^3) = 1/8  = 125000
table[4] = 1/(2^4) = 1/16 = 062500
table[5] = 1/(2^5) = 1/32 = 031250
table[6] = 1/(2^6) = 1/64 = 015625
...
table[17] = 1/(2^17) = 1/131072 = 000008
table[18] = 1/(2^18) = 1/262144 = 000004
table[19] = 1/(2^19) = 1/514288 = 000002
table[20] = 1/(2^20) = 1/1048576 = 000001

So now with that table, we can produce the whole fractional logarithm, using pure integer math:

fractionalLogarithm(value):
  log = 0
  for i = 1 to 20:
    value = shr(value * value, 6)
    if value >= 1000000 then:
      log = log + table[i]
      value = shr(value, 1)
  return log

---

Putting It All Together

Finally, for a complete logarithm of any integer your machine can represent, this is the whole thing, which will compute the logarithm with six digits of precision, in the form "0000XX.XXXXXX":

log(value):
  intPart = integerLogarithm(value)
  value = normalize(value)
  fracPart = fractionalLogarithm(value)
  result = shl(intPart, 6) + fracPart
  return result

---

Demonstration

To show that the math works — and that it works pretty well! — below is a JavaScript implementation of the above algorithm. It uses pure integer math: Only addition, subtraction, and relative comparison. Functions are used to organize the code, but they behave like subroutines: They're not recursive, and don't nest very deeply.

You can try it out live (click the 'Run' button and type 12345 in the input field). Compare the result to the standard Math.log() function, and you'll see how close the pure-integer version gets:

function shiftLeft(value) {
  var value2 = value + value;
  var value4 = value2 + value2;
  var value5 = value4 + value;
  return value5 + value5;
}

function shl(value, count) {
  while (count > 0) {
    value = shiftLeft(value);
    count = count - 1;
  }
  return value;
}

function shr(value, count) {
  if (count == 0) return value;

  var index = 11;
  var shifted = 0;
  while (index >= 0) {
    var adder = shl(1, index - count);
    var subtractor = shl(adder, count);
    while (value > subtractor) {
      value = value - subtractor;
      shifted = shifted + adder;
    }
    index = index - 1;
  }
  return shifted;
}

//-----------------------------------

function multiplyByTwo(value) {
  return value + value;
}

function multiplyByPowerOfTwo(value, count) {
  while (count > 0) {
    value = value + value;
	count = count - 1;
  }
  return value;
}

function divideByPowerOfTwo(value, count) {
  if (count == 0) return value;

  var index = 39;	// lg(floor(pow(10, 12)))
  var shifted = 0;
  while (index >= 0) {
    var adder = multiplyByPowerOfTwo(1, index - count);
    var subtractor = multiplyByPowerOfTwo(adder, count);
    while (value >= subtractor) {
      value = value - subtractor;
      shifted = shifted + adder;
    }
    index = index - 1;
  }
  return shifted;
}

function divideByTwo(value) {
  return divideByPowerOfTwo(value, 1);
}

function multiply(a, b) {
  var product = 0;
  while (b > 0) {
    nextB = divideByTwo(b);
    bit = b - multiplyByTwo(nextB);
    if (bit != 0) {
      product += a;
    }
    a = a + a;
	b = nextB;
  }
  return product;
}

//-----------------------------------

var logTable = {
   "1": 500000,
   "2": 250000,
   "3": 125000,
   "4":  62500,
   "5":  31250,
   "6":  15625,
   "7":   7813,
   "8":   3906,
   "9":   1953,
  "10":    977,
  "11":    488,
  "12":    244,
  "13":    122,
  "14":     61,
  "15":     31,
  "16":     15,
  "17":      8,
  "18":      4,
  "19":      2,
  "20":      1,
};

//-----------------------------------


function integerLogarithm(value) {
  var count = 0;
  while (value > 9) {
    value = shr(value, 1);
    count = count + 1;
  }
  return count;
}

function normalize(value) {
  var intLog = integerLogarithm(value);
  if (intLog > 5)
    value = shr(value, intLog - 5);
  else
    value = shl(value, 5 - intLog);
  return value;
}

function fractionalLogarithm(value) {
  var log = 0;
  for (i = 1; i < 20; i++) {
    var squaredValue = multiply(value, value);
    value = shr(squaredValue, 5);
    if (value >= 1000000) {
      log = log + logTable[i];
      value = shr(value, 1);
    }
  }
  return log;
}

function log(value) {
  var intPart = integerLogarithm(value);
  value = normalize(value);
  var fracPart = fractionalLogarithm(value);
  var result = shl(intPart, 6) + fracPart;
  return result;
}

//-----------------------------------

// Just a little jQuery event handling to wrap a UI around the above functions.
$("#InputValue").on("keydown keyup keypress focus blur", function(e) {
  var inputValue = Number(this.value.replace(/[^0-9]+/g, ''));
  var outputValue = log(inputValue);
  $("#OutputValue").text(outputValue / 1000000);
  var trueResult = Math.floor((Math.log(inputValue) / Math.log(10)) * 1000000 + 0.5) / 1000000
  $("#TrueResult").text(trueResult);
});

<script src="https://cdnjs.cloudflare.com/ajax/libs/jquery/3.3.1/jquery.min.js"></script>

Input integer: <input type="text" id="InputValue" /><br /><br />
Result using integer algorithm: <span id="OutputValue"></span><br /><br />
True logarithm: <span id="TrueResult"></span><br />

Answer 2

As I mentioned in your Original question on SE/RC for pow,sqrt,n-root,log,exp see:

• Power by squaring for negative exponents

and all the sub-links in there.

Once you got working *,/,<<,>> (which the other answer covers well) and may fixed point instead of floating you can also start computing goniometrics. For that the best is use Chebyshev series but as I lack the math behind them I can use only already precomputed ones ... Taylor is a common knowledge so computing that should be easy here what I code for my arithmetics template to cover math for arbitrary math data types (bignums):

// Taylor goniometric   https://en.wikipedia.org/wiki/Taylor_series
friend T sin     (const T &x)   // = sin(x)
    {
    int i; T z,dz,x2,a,b;
    x2=x/(pi+pi); x2-=::integer(x2); x2*=pi+pi;
    for (z=x2,a=x2,b=1,x2*=x2,i=2;;)
        {
        a*=x2; b*=i; i++; b*=i; i++; dz=a/b; z-=dz;
        a*=x2; b*=i; i++; b*=i; i++; dz=a/b; z+=dz;
        if (::abs(dz)<zero) break;
        }
    return z;
    }
friend T cos     (const T &x)   // = cos(x)
    {
    int i; T z,dz,x2,a,b;
    x2=x/(pi+pi); x2-=::integer(x2); x2*=pi+pi;
    for (z=1,a=1,b=1,x2*=x2,i=1;;)
        {
        a*=x2; b*=i; i++; b*=i; i++; dz=a/b; z-=dz;
        a*=x2; b*=i; i++; b*=i; i++; dz=a/b; z+=dz;
        if (::abs(dz)<zero) break;
        }
    return z;
    }
friend T tan     (const T &x)                                               // = tan(x)
    {
    int i; T z0,z1,dz,x1,x2,a,b;
    x1=x/pi; x1-=::integer(x1); x1*=pi; x2=x1*x1;
    for (z0=1,z1=1,a=1,b=1,i=2;;)
        {
        a*=x2; b*=i; i++; dz=a/b; z0-=dz;   // z0=cos(x)
               b*=i; i++; dz=a/b; z1-=dz;   // z1=sin(x)/x
        a*=x2; b*=i; i++; dz=a/b; z0+=dz;
               b*=i; i++; dz=a/b; z1+=dz;
        if (::abs(dz)<zero) break;
        }
    return (x1*z1)/z0;
    }
friend T ctg     (const T &x)                                               // = cotan(x)
    {
    int i; T z0,z1,dz,x1,x2,a,b;
    x1=x/pi; x1-=::integer(x1); x1*=pi; x2=x1*x1;
    for (z0=1,z1=1,a=1,b=1,i=2;;)
        {
        a*=x2; b*=i; i++; dz=a/b; z0-=dz;   // z0=cos(x)
               b*=i; i++; dz=a/b; z1-=dz;   // z1=sin(x)/x
        a*=x2; b*=i; i++; dz=a/b; z0+=dz;
               b*=i; i++; dz=a/b; z1+=dz;
        if (::abs(dz)<zero) break;
        }
    return z0/(x1*z1);
    }
friend T asin    (const T &x)                                               // = asin(x)
    {
    if (x<=-1.0) return -0.5*pi;
    if (x>=+1.0) return +0.5*pi;
    return ::atan(x/::sqrt(1.0-(x*x)));
    }
friend T acos    (const T &x){ T z; z=0.5*pi-::asin(x); return z; }         // = acos(x)
friend T atan    (const T &x)                                               // = atan(x)
    {
    bool _shift=false;
    bool _invert=false;
    bool _negative=false;
    T z,dz,x1,x2,a,b; x1=x;
    if (x1<0.0) { _negative=true; x1=-x1; }
    if (x1>1.0) { _invert=true; x1=1.0/x1; }
    if (x1>0.7) { _shift=true; b=::sqrt(3.0)/3.0; x1=(x1-b)/(1.0+(x1*b)); }
    for (x2=x1*x1,z=x1,a=x1,b=1;;)  // if x1>0.8 convergence is slow
        {
        a*=x2; b+=2; dz=a/b; z-=dz;
        a*=x2; b+=2; dz=a/b; z+=dz;
        if (::abs(dz)<zero) break;
        }
    if (_shift) z+=pi/6.0;
    if (_invert) z=0.5*pi-z;
    if (_negative) z=-z;
    return z;
    }
friend T actg    (const T &x){ T z; z=::atan(1.0/x); return z; }            // = acotan(x)
friend T atan2   (const T &y,const T &x){ return atanxy(x,y); }             // = atan(y/x)
friend T atanxy  (const T &x,const T &y)                                    // = atan(y/x)
    {
    int sx,sy; T a;
    T _zero=1.0e-30;
    sx=0; if (x<-_zero) sx=-1; if (x>+_zero) sx=+1;
    sy=0; if (y<-_zero) sy=-1; if (y>+_zero) sy=+1;
    if ((sy==0)&&(sx==0)) return 0.0;
    if ((sx==0)&&(sy> 0)) return 0.5*x.pi;
    if ((sx==0)&&(sy< 0)) return 1.5*x.pi;
    if ((sy==0)&&(sx> 0)) return 0.0;
    if ((sy==0)&&(sx< 0)) return x.pi;
    a=y/x; if (a<0) a=-a;
    a=::atan(a);
    if ((sx>0)&&(sy>0)) a=a;
    if ((sx<0)&&(sy>0)) a=x.pi-a;
    if ((sx<0)&&(sy<0)) a=x.pi+a;
    if ((sx>0)&&(sy<0)) a=x.pi+x.pi-a;
    return a;
    }

As I mentioned you need to use floating or fixed point for this as the results are not integers !!!

But as I mentioned before CORDIC is better suited for computing on integers (if you search there where some QAs here on SE/SO with C++ code for this).

IIRC it exploit some (arc)tan angle summation identity that leads to a nicely computable on integers delta angle something like sqrt(1+x*x) which is easily computable on integers. With binary search or approximation/iteration you can compute the tan of any angle and using goniometric identities you can compute any cotan sin and cos ... But I might be wrong as I do not use CORDIC and read about it a long time ago

Anyway once you got some function its inverse can be usually computed with binary search.