algorithm - 用于计算三角函数，对数等的算法。仅加减

Question

我正在还原Ascota 170古董机械可编程计算机。它已经在工作了。
现在，我正在寻找一种算法来证明其功能-例如计算三角表或对数表。或类似的东西。
不幸的是，从数学运算来看，计算机只能加和减整数（-1E12至1E12中的55个寄存器）。甚至没有移位到数字的运算-因此可以通过编程将其实现为仅乘以非常小的数字。
但是它的逻辑运算已经非常完善。

您能建议我任何合适的算法吗？

Answer 1

因此，您正在做的事情确实很棒。碰巧的是，我可以解释很多关于如何仅使用整数加减运算来实现分数对数的方法！这篇文章将很长，但是其中包含很多细节，最后还有一个可行的实现，对于您用怪异的机械计算机做一些有趣的事情应该足够了。



实施比较

您将需要能够比较数字。虽然您说可以执行== 0和> 0的比较，但是对于您要实现的大多数有趣算法来说，这还远远不够。您需要相对比较，可以通过减法确定：

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

在本文的其余部分中，我将只编写更简单的 a > b形式，但是如果您不能直接做到这一点，则可以替代上面的一种操作。



实施轮班

现在，由于您没有数字移位硬件，因此必须创建“例程”来实现它。左移很容易：给自己加上一个数字，一次又一次，然后加上原始数字，再加上一次。这相当于向左移一位数。

因此，左移一位或乘以十：

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

多次移位只是对 shiftLeft()的重复调用：

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

向右移动一位数字会有点困难：我们需要通过反复的减法和加法来做到这一点，如下面的伪代码所示：

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

方便地，由于一开始很难向右移位，因此该算法使我们可以直接选择要移位的位数。



乘法

看来您的硬件可能有乘法？但是，如果没有，您可以使用重复的加法和移位来实现乘法。二进制乘法是最有效的最简单实现形式，它要求我们首先使用与实现 multiplyByTwo()和 divideByTwo()相同的基本技术来实现 shiftLeft()和 shr()。

一旦实现了这些，乘法就涉及重复切出其中一个数字的最后一位，如果该位是 1，则将另一个数字的递增版本添加到运行总计中：

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

如果需要的话，下面将提供完整的实现。



整数对数

我们可以使用向右移动一位数字的能力来计算数字的以10为底的对数的整数部分–实际上，这是您可以将数字右移多少次，直到数字太小而无法移位。

integerLogarithm(value):

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

因此，对于0-9，这将返回0；对于0-9，则返回0。对于10-99，返回1；对于100-999，则返回2，依此类推。



整数指数

上面算法的反面是微不足道的：要计算10的整数次幂，我们只需要将幂左移几位即可。

integerExponent(count):

    value = shl(1, count)
    return value

因此，对于0，它将返回1；对于0，将返回1。为1，则返回10；对于2，返回100；对于3，返回1000；等等。



分割整数和分数

现在我们可以处理整数幂和对数了，我们几乎可以处理小数部分了。但是，在我们真正谈论如何计算对数的小数部分之前，我们必须谈论如何对问题进行除法，以便我们可以将整数部分与小数部分分开计算。理想情况下，我们只想处理固定范围内数字的计算对数-例如，从1到10，而不是从1到无穷大。

我们可以使用整数对数和指数例程对完整的对数问题进行切分，以便无论输入数字是多少，我们始终要处理[1，10）范围内的值。

首先，我们计算整数对数，然后计算整数指数，然后从原始数中减去它。剩下的就是我们需要计算的小数部分：然后剩下的唯一工作就是移动该小数部分，以使其始终处于一致的范围内。

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

您可以毫不费力地说服自己，无论原始值是多少，其最高非零数字都将移至第7列：因此，“ 12345”将变为“ 000000123450”（即“ 0000001.23450”）。这使我们能够假装总是存在一个不可见的小数点，该数字比数字的一半还小，因此现在我们只需要解决计算[1，10）范围内的对数的问题。

（为什么“超过一半”？我们将需要值的上半部分始终为零，一会儿您就会明白为什么。）



小数对数

克努斯（Knuth）在《计算机编程的艺术》第1.2.2节中说明了如何执行此操作。我们的目标是计算log10（x），以便对于 b1， b2， b3 ...的某些值，其中 n已经是 0（因为我们将上面的整数部分除了） ：


log10（x）= n + b1 / 2 + b2 / 4 + b3 / 8 + b4 / 16 + ...


Knuth说我们可以像这样获得 b1， b2， b3 ...：


为了获得b1，b2，...，我们现在设置x0 = x / 10 ^ n，并且对于k> = 1，

如果x [k-1] ^ 2 <10，则b [k] = 0，x [k] = x [k-1] ^ 2；

b [k] = 1，如果x [k-1] ^ 2> = 10，则x [k] = x [k-1] ^ 2/10。


也就是说，每个步骤都使用伪代码循环，如下所示：

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

为了使用上面的定点数字来实现此目的，我们必须使用移位将小数点移回原位的方式来实现 x * x，这将丢失一些数字。如Knuth所说，这将导致错误传播，但是它将提供足够的准确性，足以用于演示目的。

因此，给定由 normalize(value)生成的分数值，我们可以像这样计算其分数二进制对数：

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)

但是，二进制小数对数-单个位！ —并不是特别有用，尤其是因为我们在前面的步骤中计算了对数整数部分的十进制形式。因此，我们将再次修改此时间，以计算十进制小数对数到五个位，而不是计算位数组；为此，我们需要一个由20个值组成的表格，这些表格表示每个位到十进制的转换，并将它们存储为定点数：

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

因此，现在有了该表，我们可以使用纯整数数学运算来生成整个分数对数：

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

放在一起

最后，对于您的机器可以代表的任何整数的完整对数，这就是全部，它将以六位数的精度计算对数，形式为“ 0000XX.XXXXXX”：

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

示范

证明数学是有效的，并且效果很好！ —以下是上述算法的JavaScript实现。它使用纯整数数学：仅加法，减法和相对比较。函数用于组织代码，但它们的行为类似于子例程：它们不是递归的，并且嵌套不深。

您可以现场试用（单击“运行”按钮，然后在输入字段中键入 12345）。将结果与标准 Math.log()函数进行比较，您将看到纯整数版本的接近程度：

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 />