[数值计算] [代码讨论]批处理开平方(计算平方根的值)

SQYSQYSQY

致浏览者，想对大数字（超过1000）计算平方根，请见31楼（约18秒算到150位），想对小数计算平方根，请见34楼（是个文件）（约23秒算到150位）
注意，部分计算器对最后一位进行了四舍五入，本程序可能会与计算器的最后一位相差1。这不是本程序的错误。
致浏览者，被开方数不超过1000，可用下面高效率程序进行计算（约3.5秒算到150位）。

@echo off & SetLocal EnableDelayedExpansion
::变量c的值是被开方数，不得超过1000！
set c=2
for %%a in (16 8 4 2 1) do (
        set /a "d=(a*%%a<<1)+%%a*%%a+b"
        if !c! geq !d! set /a "a+=%%a","b=d"
)
if !b! equ !c! echo !a! & pause & exit /b
set /p "=!a!."<nul
set "d=0"
for %%a in (512 256 128 64 32 16 8 4 2 1) do (
        set /a "d+=%%a","e=((a*d<<1)+d*d/1000)/1000+b"
        if !e! geq !c! set /a "d-=%%a"
)
set "e=00!d!" & set /p "=!e:~-3!"<nul
set /a "e=(!a!!e:~-3!*!a!!e:~-3!)%%1000000","b=e/1000","e=e%%1000"
set "e=  !e!" & set "e=!e:~-3!"
set "b=000!e!!b!"
set "e=  !d!"
set "a=!e:~-3!!a!"
for /l %%a in (3 3 2147483646) do (
        set "d=0"
        for %%b in (512 256 128 64 32 16 8 4 2 1) do (
                set /a "d+=%%b"
                if 1000 gtr !d! (
                        set /a "c=d*d/1000"
                        for /l %%c in (0 3 %%a) do set /a "c=((!a:~%%c,3!*d<<1)+c+!b:~%%c,3!)/1000"
                        set /a "e=%%a+3"
                        for %%c in (!e!) do set /a "c+=!b:~%%c,3!" & if !c! geq 1000 set /a "d-=%%b"
                ) else set /a "d-=%%b"
        )
        set "e=00!d!" & set /p "=!e:~-3!"<nul
        set "c=" & set /a "f=d*d/1000"
        for /l %%b in (0 3 %%a) do (
        set /a "e=(!a:~%%b,3!*d<<1)+!b:~%%b,3!+f","f=e/1000","e%%=1000"
        set "e=  !e!" & set "c=!c!!e:~-3!"
        )
        set /a "e=(d*d)%%1000" & set "e=  !e!" & set "b=000!e:~-3!!c!"
        set "e=  !d!" & set "a=!e:~-3!!a!"
)
echo. & pause & exit /b
复制代码

老刘1号

回复 49# 523066680


    是啊，准备有空的时候把happy的代码研究一遍

523066680

回复 48# 老刘1号


      这都是牛顿的功劳呀，牛迭是N次方，1/N次方都可以算，很多其他方程的根也可以算。因为他是在一个函数曲线任意一点(x)上求切线，然后这个切线不断迭代，向函数曲线和x轴的交点（根或者函数的解）逼近。
逼近速度非常快（差不多精度翻倍）。

现在回想，happy886r 当时写了好多大工程
happy886r - 数学计算工具 i 的重构版new i

老刘1号

回复 46# SQYSQYSQY
回复 47# 523066680


    哈哈发现了，我大致扫了一遍回帖没注意到，就贴了一下
没专门研究过这个算法，不知是牛跌

523066680

回复 44# 老刘1号

      这个牛跌前面楼层有提到过呀。
我们主要在折腾大数字和浮点数。如果只在signed int范围内开整数根，早就可以结帖了。
小程说他要写一个版本的时候我以为“硬核”要来了，结果 ……
那本《MCA》非常硬核！打算用C艹实践其中一小部分算法（大部分看不懂，小部分似懂非懂）。

回复 43# SQYSQYSQY

批处理的“数组”属于伪数组，和其他语言的数组有本质的区别，因为批的“索引”不只是数字，还有字母和符号，变量长度亦不确定，这和其他语言的哈希表（hash table）相似。
在C/C++里面的数组操作不需要“搜索”，可直接通过计算偏移量取内存地址，因为定长的数组在初始化后占用一段连续内存空间，且每个单元占用相同字节，
给定一个编号，通过 编号*字节大小+起点地址 可得目标内存地址，直接存取。
C++ vector容器的 [] 操作符消耗占用高，和其独立的实现有关（可能做了各种判断和转换处理），换用更纯粹的C语言数组可以减少这种消耗。

实测，元素数量为80W的容器/数组，1000次遍历每一个元素并写入int值，

#include <iostream>
#include <vector>
#include <chrono>
using namespace std;
using namespace std::chrono;

const int SIZE = 800000;
void vec_test(void);
void c_array_test(void);
void c_pointer_test(void);
void time_used(system_clock::time_point& time_a);

vector<int> vec(SIZE);
int array1[SIZE];
int array2[SIZE];

int main(int argc, char *argv[])
{
    system_clock::time_point start = system_clock::now();
    for (int i = 0; i < 1000; i++) vec_test();
    time_used(start);

    for (int i = 0; i < 1000; i++)  c_array_test();
    time_used(start);

    for (int i = 0; i < 1000; i++)  c_pointer_test();
    time_used(start);
    return 0;
}

void vec_test(void) {
    register int it;
    for (it = 0; it < SIZE; it++) vec[it] = it;
}

void c_array_test(void) {
    register int it;
    for (it = 0; it < SIZE; it++) array1[it] = it;
}

void c_pointer_test(void) {
    register int it;
    register int *pt = array2;
    for (it = 0; it < SIZE; it++) *(pt + it) = it;
}

void time_used(system_clock::time_point& time_a) {
    duration<double> diff;
    diff = chrono::system_clock::now() - time_a;
    cout << "Time used: " << diff.count() << endl;
    time_a = chrono::system_clock::now();
}
复制代码

g++编译测试结果：
Time used: 1.24807
Time used: 0.400023
Time used: 0.393022

Visual Studio编译，差距更明显
Time used: 3.92722
Time used: 0.0410024
Time used: 0.0360021


如果一定要用 vector 又不想用它的 [] ，可以申请一个指针，int *vp = vec.data()
通过*vp指针直接算地址读写内存，速度和c_array一样快，想要更快，申请一个寄存器指针 register int *vp。
这就是自由度，可定制，可接管。

讨论算法，可以不分语言。谈论模块化思想，也可以不分语言，用批处理同样可以表达。
但要说极限压榨性能，怎么也轮不到批处理。

老刘1号

这贴的楼层有毒吧，看的好蛋疼

老刘1号

贴吧吧友写的（Q群里聊天得知，感谢小程）
可开开方后结果为整数的整数

@set /p a=输入a=
@set /a b=a
@for /l %%a in (1,1,30) do @set /a b=(b+a/b)/2
@set b
复制代码

523066680

用了两天 VS2015，不愧是宇宙最强编辑器（是我墨守成规，一直装在系统，从来不用）

s_minus 函数和 vec_minus 函数性能分析，1W位数，各调用2W次的分析结果：





可以看到涉及数组[]操作的语句消耗较高，采用高的进制(base)处理势在必行。

再聊聊《Modern Computer Arithmetic》（后面简称MCA）
对于64位平台 unsigned long long int 支持的最大数字是 2^64-1=18446744073709551615，20位，如果我们充分利用，使用2^64作为进制，或者10^19作为进制，
很容易会遇到溢出问题，MCA中给出了几种方案：
Let T be the number of different values taken by the data type representing the coefficients ai, bi. (Clearly, β ≤ T, but equality does not necessarily hold,
for example β = 10^9 and T = 2^32.) At step 3, the value of s can be as large as 2β − 1, which is not representable if β = T. Several workarounds
are possible:
either use a machine instruction that gives the possible carry of ai + bi,
or use the fact that, if a carry occurs in ai + bi, then the computed sum – if performed modulo T – equals t := ai +bi −T < ai; thus, comparing t and ai will determine if a carry occurred.
A third solution is to keep a bit in reserve, taking β ≤ T/2.

用 T 表示一个存储单元所能表示的数的量，则有 β <= T （这里 β 表示采用的进制，以及β不一定等于T，例如T=2^32，但采用的进制为10^9）。
考虑加法操作 s=a+b+d （d为1或0，是上一次加法补进的数值），s的最大可能值为 2*β-1，当 β = T 时该公式无法正确计算。考虑以下方案：
1. 内部编码实现 （水平有限，暂时这么翻译）
2. 通过-T取余数判断，t := a + b - T < a ，对比 t 和 a 断定是否进 1
3. 采用一个保守的进制数 β，令 β ≤ T/2.

523066680

回复 28# search_Sudoku

    将手算法原样撸了一个C++初始版本，未优化，1W精度（含输出）0.9秒。

返回来看gmp库的说明，这理论好深（估计happy能看懂，我只能看看概述）
Introduction
The current asymptotically fastest known method to compute the square-root of a n-word number is using Fast Fourier Transform (FFT) multiplication and Newton'smethod,
with a complexity of 5M(n)' [1, p. 155].2 The algorithm presented here is basedon Burnikel-Ziegler Karatsuba division [2].
这里面用到 快速傅里叶变换乘法、牛顿迭代法，以及 Burnikel-Ziegler Karatsuba 除法。

牛顿迭代法的精度增长很快（差不多翻倍增长），需要用到大数除法。实践起来，具体的精度不太好确定，因此迭代到一定次数后，除法的精度取舍也是个问题，
我也发现引用了牛顿迭代法的方案大多是给出一个预判的精度，
比如 gmp 中设置 mpf_set_default_prec (10000);  开根后 gmp_printf ("Square root: %.10000Ff\n\n", sq_out);
实际有效数字只有3068位，末尾由大量的0填充。

几种接口效率对比(i7 8核 4GHz)
Perl bignum 模块，1W精度 耗时4.5s
Python Decimal 模块，10W精度 耗时 1.0s
C++ 使用 GMP库，100W有效精度，含输出，1.5秒；不含输出，0.2秒
                           1000W位有效精度，不含输出， 2.5秒。

libgmp太凶悍

C++ 测试GMP速度以及判断有效精度的代码：

#include <iostream>
#include <sstream>
#include <iomanip>
#include <chrono>
#include "gmpxx.h"
#include "gmp.h"

using namespace std;
int main(int argc, char *argv[] ) 
{
    ostringstream os;
    mpf_class n(2, 10), r(1, 3600000);
    auto start = chrono::system_clock::now();

        r = sqrt(n);
        os << setprecision(10000000) << r;

    auto end = chrono::system_clock::now();
    chrono::duration<double> diff = end-start;
    cout << "Actual length of r: " << os.str().length() << endl;
    cout << "Time Used: " << setprecision(3) << diff.count() << " s" << endl;
    return 0;
}
复制代码

Actual length of r: 1083710
Time Used: 0.561 s （这里包含 os << r 产生的时间）

GNU MP 库中用到的算法和理论，完整版：
《Modern Computer Arithmetic》
http://maths.anu.edu.au/~brent/pd/mca-cup-0.5.9.pdf

523066680

回复 40# SQYSQYSQY

    先写一个测试模块，发布前测试过会比较妥。更新后的代码测试26，37，50，65，82不通过。更新前测试通过（精确到80位）。

for /l %%a in (1,1, 100%或者其他数字%) do ( call :开根函数 )
exit/b

:开根函数
setlocal
endlocal& goto :eof
复制代码

523066680

回复 37# SQYSQYSQY

      关于自由度，我不说数独，你用批处理做一做猜数字游戏的搜索树就知道了，但凡跟递归、多重镶嵌数据结构有关的东西，用批处理都不划算。
http://bbs.bathome.net/thread-49436-1-2.html
（非杠，我就是觉得你学其他语言，起跑线就会有一个质的飞跃。也不推荐Perl，现在的主流推荐就是Python，Java，C）

百度百科的描述
https://baike.baidu.com/item/%E7 ... 97/83200?fr=aladdin

举一个简单的例子：
我要生成50个随机字符，内容可以是数字和小写字母；将字符串拆割重组、输出显示整个数组：

use Modern::Perl;
# 生成 50 个随机字符，范围是数字和小写字母
my @arr = map { ('0'..'9', 'a'..'z')[rand(36)] } (1 .. 50 );
my $s = "bathomenet";
say join(",", @arr);
say join(",", split //, $s);
print length($s);
复制代码

输出

j,x,c,z,7,k,u,p,d,w,u,k,d,0,s,9,l,a,7,n,i,9,h,s,l,5,1,q,d,m,j,i,p,m,3,x,4,h,s,d,r,s,3,q,w,i,d,z,f,y
b,a,t,h,o,m,e,n,e,t
10
复制代码

这些在 python ruby perl 里面都是常规操作，字符串长度允许范围几乎就是对应内存的剩余空间。

批处理可以做到，但是会很烦，vbs同样烦，js、Powershell应该好写

回复 39# SQYSQYSQY

      我也是业余的，工作和IT不相关。可以预见的是你的路线可以是 x^2。
而我这种职业不相关的，路线基本就只能是平缓直线甚至是上下波动了

523066680

1月17日
补充 —— 是楼层穿越问题。此信息PASS

523066680

回复 40# SQYSQYSQY

      你这起跑线多好，我那会儿没零钱，严重的时候一周只有1元上一个小时网吧，要做点什么得在稿纸上写好，不然刷的一小时没了。

空间换时间/时间换空间
查表就是为了"空间换时间"，几百个上千元素的数组查询都会慢的也就是批处理了，不如趁早换用自由度高的语言。

那个例子可以用函数表达，x^2，和 5*x 赛跑，一开始 5*x 比较快，后来 x^2 …… 这俩太陡了，应该有更合适的曲线，暂时先这样


浮点数版本补充在39楼。

523066680

回复 38# SQYSQYSQY

非常短（难道是受到了27楼小朋友的刺激？他短任他短，清风拂山冈。他横任他横，明月照大江）。可能有小细节需要完善：
1.414找不到操作数。
999999999999999999

    不打算先把大数和浮点数做进去再精简吗。

浮点数版本，效率老样子，代码不打算优化了，要也用Perl，或者换C艹

523066680

回复 35# SQYSQYSQY

    查阅 ASCII 字符表，如果想学习更多关于字符编码，建议：学习 python，用python去读写各种网页文件。
关键词：Encode, Decode, Unicode, UTF8, UTF16, GBK, CP936

用其他语言显示数字符号对应的ASCII码。

>perl -e "grep {printf qq('%d' - %d\n), $_, ord($_) } (0..9)"
'0' - 48
'1' - 49
'2' - 50
'3' - 51
'4' - 52
'5' - 53
'6' - 54
'7' - 55
'8' - 56
'9' - 57
复制代码