手算方法的思想是利用之前已求得的所有位,求出下一个小数位。这种思想在算法领域有着广泛的应用,比如动态规划思想。
设之前已求出的整数部分和小数部分为a
,对a
不断乘以10,直到它是一个整数,仍记为a
。这一过程即移除小数点的过程。比如:设之前已经求出a = 11.45
,那么小数点移除后得a = 1145
。a
的初值为平方根的整数部分intPart
,用二分法确定。
设v
的初值为待求平方根的数,每次a
乘以10,它也要乘以100。
设要求的下一小数位为b
。于是b
是满足v >= (a + 0.1*b) ^ 2
的最大自然数。
于是100 * v >= (10a+b)^2 = 100*a*a+20*a*b+b*b
。设每一轮的差值delta = v – a*a
,于是100*delta >= 20*a*b+b*b
。
因此我们只需要找到最大自然数b
,满足100*delta >= 20*a*b+b*b
。为了实现的简洁,确定b
时我们可以直接遍历每个候选并计算。
最后确定delta
的更新式。next_v = 100 * v
,next_a = 10 * a + b
,故next_delta = next_v - next_a * next_a
。但我们也可以把它变换成一个没有“next_”的更新式。next_delta = 100*v-100*a*a-(20*a*b+b*b) = 100*delta - (20*a*b+b*b)
。
它相比于牛顿迭代等没有计算速度上的优势,牛顿迭代在此问题是平方收敛,它的速度就慢太多了。但如果只使用通常的double,牛顿迭代无法获得极其高的精度(比如小数点后114514位)。而手算方法在使用BigInt的情况下可以达到任意精度。
伪代码
function mySqrt(v,want){
const intPart = 二分法求得平方根的整数部分
let delta = v - intPart * intPart
if (!delta) return intPart// 是完全平方数
let a = intPart
let ans = []
for (let i = 1; i <= want; ++i) {
let val = 最大的b,使100 * delta >= 20 * a * b + b * b
ans.push(val)
delta = 100 * delta - (20 * a * val + val * val)
a = 10 * a + val
}
return `${intPart}.${ans.join('')}`
}
代码
"use strict";
//node 开根号.js
//若v是完全平方数,则返回整数,否则返回want位小数
function mySqrt(v,want) {
let l = 1, r = v
while (l < r - 1) {
let mid = (l + r) >> 1
if (mid * mid > v) r = mid - 1
else l = mid
}
const intPart = r * r <= v ? r : l
let delta = v - intPart * intPart
if (!delta) return intPart
let tmp = intPart
let ans = []
for (let i = 1; i <= want; ++i) {
let val = 9
for (let j = 0; j < 10; ++j) {
if (20 * tmp * j + j * j > 100 * delta) {
val = j - 1
break
}
}
ans.push(val)
delta = 100 * delta - (20 * tmp * val + val * val)
tmp = 10 * tmp + val
// console.log(delta, tmp)//因为tmp会越来越大,所以delta大体上也会越来越大
}
return `${intPart}.${ans.join('')}`
}
//1.4142135623730 41 1.7320508075688 5.0990195135927
//2.2360679774997 2.4494897427831 2.6457513110645 338.3991725758
console.log(mySqrt(2,13))
console.log(mySqrt(1681,13))
console.log(mySqrt(3,13))
console.log(mySqrt(26,13))
console.log(mySqrt(5,13))
console.log(mySqrt(6,13))
console.log(mySqrt(7,13))
console.log(mySqrt(114514,10))
接下来推广至任意次方根。对于立方根求解问题,模仿以上做法,1000*v >= 1000*a^3 + 300*a^2*b + 30*a*b^2 + b^3
,找出1000*delta >= 300*a^2*b + 30*a*b^2 + b^3
的最大b
即可。delta
更新式同理。
求其他次方根,可继续使用二项式定理展开,得到式子,但意义不大。
代码
"use strict";
// 若v是完全立方数,则返回整数,否则返回want位小数
function cubeRt(v, want) {
let p3 = x => x * x * x
let l = 1, r = v
while (l < r - 1) {
let mid = (l + r) >> 1
if (p3(mid) > v) r = mid - 1
else l = mid
}
const intPart = p3(r) <= v ? r : l
let delta = BigInt(v - p3(intPart))
if (!delta) return intPart
let a = BigInt(intPart)
let ans = []
for (let i = 1; i <= want; ++i) {
let val = 9n
const d1000 = 1000n * delta
const c1 = 300n * a, c2 = 30n * a
for (let b = 1n; b < 10n; ++b) {
if (c1 * a * b + c2 * b * b + p3(b) > d1000) {
val = --b
break
}
}
ans.push(val.toString())
delta = d1000 - (c1 * a * val + c2 * val * val + p3(val))
a = 10n * a + val
}
return `${intPart}.${ans.join('')}`
}
// 1
// 1.25992104989487316476721060727822835057025146470150
// 1.44224957030740838232163831078010958839186925349935
// 1.58740105196819947475170563927230826039149332789985
// 1.70997594667669698935310887254386010986805511054305
// 1.81712059283213965889121175632726050242821046314121
// 1.91293118277238910119911683954876028286243905034587
// 2
// 2.08008382305190411453005682435788538633780534037326
// 2.15443469003188372175929356651935049525934494219210
// 2.96249606840737050867306218934183853756635742231886
// 3
// 3.03658897187566251942080957850566963558145397724811
// 11.89020213687269261757960016357431082984522211204992
// 48.560840499788364316317086175904910280469734732
console.log(cubeRt(1, 50))
console.log(cubeRt(2, 50))
console.log(cubeRt(3, 50))
console.log(cubeRt(4, 50))
console.log(cubeRt(5, 50))
console.log(cubeRt(6, 50))
console.log(cubeRt(7, 50))
console.log(cubeRt(8, 50))
console.log(cubeRt(9, 50))
console.log(cubeRt(10, 50))
console.log(cubeRt(26, 50))
console.log(cubeRt(27, 50))
console.log(cubeRt(28, 50))
console.log(cubeRt(1681, 50))
console.log(cubeRt(114514, 45))