我正在修复 Ascota 170 古董机械可编程计算机。它已经在工作了。 现在我正在寻找一种算法来展示它的功能——比如计算三角函数表或对数表。或类似的东西。 不幸的是，从数学运算来看，计算机只能加减整数（从 -1E12 到 1E12 的 55 个寄存器）。甚至没有移位到数字的操作——因此它可以通过编程实现只乘以非常小的数字。 但是它的逻辑运算非常发达。
你能告诉我任何合适的算法吗？
【问题讨论】：
所以你的所作所为真的很棒。碰巧，我可以解释一下如何仅使用整数加减法来实现小数对数！这篇文章会很长，但包含很多细节，最后还有一个可行的实现，应该足以让你用你奇怪的机械计算机做一些有趣的事情。
实施比较
您将需要能够比较数字。虽然您说您可以执行比较 == 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；对于 10-99，返回 1；对于 100-999，这将返回 2，依此类推。
整数指数
上述算法的反面非常简单：要计算 10 的整数幂，我们只需将数字左移幂即可。
integerExponent(count):
value = shl(1, count)
return value
所以对于 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，
b[k] = 0, x[k] = x[k-1] ^ 2, 如果 x[k-1] ^ 2
b[k] = 1，x[k] = x[k-1] ^ 2 / 10，如果 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 />【讨论】：
正如我在您的Original question on SE/RC 中提到的pow,sqrt,n-root,log,exp 看到：
以及其中的所有子链接。
一旦你开始工作*,/,<<,>>（另一个答案很好地涵盖了）并且可以固定点而不是浮动，你也可以开始计算测角。为此最好的是使用切比雪夫系列，但由于我缺乏它们背后的数学，我只能使用已经预先计算的那些......Taylor 是一个常识，所以计算应该很容易在这里我为我的算术模板编写代码以涵盖数学对于任意数学数据类型（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;
}
正如我提到的，您需要为此使用浮点数或定点数，因为结果不是整数！！！
但正如我之前提到的，CORDIC 更适合于整数计算（如果你在 SE/SO 上搜索一些 QAs 并为此使用 C++ 代码）。
IIRC 它利用了一些 (arc) 正切角求和恒等式，可以很好地计算整数 delta 角，例如 sqrt(1+x*x)，它很容易在整数上计算。通过二分搜索或近似/迭代，您可以计算任何角度的tan，并使用测角恒等式计算任何cotan sin 和cos ... 但我可能错了，因为我没有使用 CORDIC 并在很久以前阅读过它
无论如何，一旦你得到一些函数，它的逆通常可以用二分搜索来计算。
【讨论】：