523066680

改好了，去掉了二分搜索，5位以内的情况下粗暴遍历。
通过0到100的数字开根测试，支持大数字。在这个基础上再去优化会舒服很多，但还是要做出浮点数开根以后再说。

CPU: i7 4GHz
SQRT 2 精度 80 耗时约 2.3s。精度 300，耗时约 20s

:: Bignum(integer) Square Root, Decimal Solution
:: https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Decimal_(base_10)
:: 523066680/vicyang
:: 2019-01

@echo off
setlocal enabledelayedexpansion
:init
    rem 创建用于计算字符串长度的模板，长度限制为 2^pow
    set "sharp=#"
    set "serial=9876543210"
    set /a pow=11, maxlen=1^<^<pow
    for /l %%a in (1,1,%pow%) do set sharp=!sharp!!sharp!

set precision=80
rem call :check_one 2
call :check_all
pause
exit /b

:: 独立测试
:check_one
    set ta=!time!
    rem call :check_first %1 !precision!
    call :decimal_solution %1
    call :time_delta !ta! !time! tu
    echo time used: !tu!
    goto :eof

:: 批量测试
:check_all
    for /l %%a in (1,1,99) do (
        echo test number: %%a
        rem call :check_first %%a !precision!
        call :decimal_solution %%a
        echo,
    )
    goto :eof

:: 使用其他工具校验/对比结果
:check_first
    perl -Mbignum=p,-%2 -le "print sqrt(%1)" 2>nul
    goto :eof

:: 手算开根方案
:decimal_solution
    setlocal
    set num=%1
    set tnum=%1
    :: 计算长度，判断需要截取的目标长度（1 or 2）
    call :length %num% len
    set /a mod=len %% 2, skip = 2 - mod
    set target=!tnum:~0,%skip%!
    set tnum=!tnum:~%skip%!
    set "base="

    :: prec 当前精度
    set /a prec = 0, base_len=0, target_len=skip

    :dec_loop
        :: 推算下一个数
        :estimate
            :: 如果目标值 小于 基数，下一个数字判定为0
            call :cmp %target% %base%0 %target_len% %base_len%+1 cmp
            if !cmp! equ -1 (
                set /a mid=0, mp=0, mplen=0
                goto :out_estimate
            )

            if %base_len% gtr 5 (
                set /a est=!target:~0,6!/!base:~0,5!
            ) else (
                :: 在set/a计算范围内的，[粗暴]遍历
                for /l %%a in (0,1,10) do (
                    set /a mp=%base%%%a*%%a
                    if !mp! gtr !target! (set /a est=%%a-1 &goto :out_est_for)
                )
            )
            :out_est_for

            :: 199999996400/1999999988 = 99.9999988
            :: but 199999/19999 = 10
            if %est% geq 10 (
                set /a tbase_len=base_len+1
                if !target_len! gtr !tbase_len! (set /a est=9)
            )

            set /a mid=!est:~0,1!
            call :bignum_mp_single !base!!mid! !mid! !base_len!+1 1 mp mplen
            call :cmp !mp! !target! !mplen! !target_len! cmp
            :: 如果mp超出目标范围
            if !cmp! equ 1 (
                set /a mid-=1
                call :bignum_mp_single !base!!mid! !mid! !base_len!+1 1 mp mplen
            )
        :out_estimate

        set /p inp="%mid%"<nul
        rem echo,&echo tg !target!, mp !mp!, base !base!, mid !mid!, est !est!
        if "%tnum%" == "" (
            :: 如果target只剩下 00，方案结束
            if "%target%" == "00" ( goto :dec_loop_out )
            if %cmp% == 0 ( goto :dec_loop_out )
        )

        :: 计算下一段target的值
        call :bignum_minus %target% %mp% %target_len% %mplen% target target_len

        :: 扩充target，如果被开根数已经截取完，直接补0，精度+1
        if %skip% geq %len% (
            set target=%target%00
            set /a prec+=1
            if !prec! equ 1 set /p inp="."<nul
        ) else (
            if "%target%" == "0" (set target=!tnum:~0,2!
                          ) else (set target=!target!!tnum:~0,2!)
            set tnum=!tnum:~2!
            set /a skip+=2
        )
        set /a target_len+=2

        :: 更新基数 - base
        rem base=base*10+mid*2
        if "%base%" == "0" (
            set /a base=mid*2, base_len=1+base/10
        ) else (
            set /a db_mid=mid*2, dbmidlen=1+db_mid/10
            call :bignum_plus !base!0 !db_mid! !base_len!+1 !dbmidlen! base base_len
        )

    if %prec% leq %precision% (goto :dec_loop)
    :dec_loop_out
    echo,
    endlocal
    goto :eof

:: 比较
:cmp
    set /a La=%3, Lb=%4
    if %La% gtr %Lb% (set /a %5=1&goto :eof)
    if %La% lss %Lb% (set /a %5=-1&goto :eof)
    :: 如果长度相同，直接按字符串对比
    if "%1" gtr "%2" (set /a %5=1&goto :eof)
    if "%1" lss "%2" (set /a %5=-1&goto :eof)
    if "%1" equ "%2" (set /a %5=0&goto :eof)

:: 大数 乘以 单位数
:bignum_mp_single
    setlocal
    set num_a=%1
    set num_b=%2
    set /a pool = 0, maxid = %3
    set "res="
    for /l %%a in ( 1, 1, %maxid% ) do (
        set /a mp = !num_a:~-%%a,1! * num_b + pool, t = mp %% 10, pool = mp / 10
        set res=!t!!res!
    )

    if %pool% neq 0 (
        set /a maxid+=1
        set res=!pool!!res!
    )
    endlocal&set %5=%res%&set %6=%maxid%
    goto :eof

::大数加法
:bignum_plus
    setlocal
    set num_a=%1
    set num_b=%2
    set /a len_a=%3, len_b=%4
    set /a max = len_a
    if %len_b% gtr %len_a% (set /a max=len_b, len_b=len_a&set num_a=%num_b%&set num_b=%num_a%)
    set /a pool=0
    set res=
    for /l %%n in ( 1, 1, %max% ) do (
        if %%n leq %len_b% (
            set /a t = !num_a:~-%%n,1! + !num_b:~-%%n,1! + pool
        ) else (
            set /a t = !num_a:~-%%n,1! + pool
        )
        set /a mod = t %% 10, pool = t / 10
        set res=!mod!!res!
    )
    if %pool% gtr 0 (set /a max+=1 &set res=1%res%)
    endlocal &set %5=%res%&set %6=%max%
    goto :eof

::大数减法
:bignum_minus
    setlocal
    set num_a=%1
    set num_b=%2
    set /a len_a=%3, len_b=%4
    set /a max = len_a
    if %len_b% gtr %len_a% (set /a max=len_b, len_b=len_a&set num_a=%num_b%&set num_b=%num_a%)

    set /a minus = 0
    set "res="
    for /l %%n in ( 1, 1, %max% ) do (
        if %%n leq %len_b% (
            set /a dt = !num_a:~-%%n,1! - !num_b:~-%%n,1! - minus
        ) else (
            set /a dt = !num_a:~-%%n,1! - minus
        )
        if !dt! lss 0 (
            set /a t = dt + 10, minus=1
        ) else (
            set /a t = dt, minus=0
        )
        set res=!t!!res!
        if !t! equ 0 (set /a zero+=1) else (set /a zero=0)
    )
    set res=!res:~%zero%!
    endlocal &set %5=%res%&set /a %6=%max%-%zero%
    goto :eof

::字符串长度计算
:length %str% %vname%
    setlocal
    set test=%~1_%sharp%
    set test=!test:~0,%maxlen%!
    set test=%test:*_=%
    set /a len=maxlen-(%test:#=1+%1)
    endlocal &set %2=%len%
    goto :eof

:: plp626的时间差函数 时间跨度在1分钟内可调用之；用于测试一般bat运行时间
:time_delta <beginTimeVar> <endTimeVar> <retVar> // code by plp626
    setlocal
    set ta=%1&set tb=%2
    set /a "c=1!tb:~-5,2!!tb:~-2!-1!ta:~-5,2!!ta:~-2!,c+=-6000*(c>>31)"
    if defined %3 set /a c+=!%3!
    endlocal&set %3=%c%
    goto:eof
复制代码

523066680

回复 33# SQYSQYSQY

是的哟，if 在某些情况下会进行字符串对比。bbbbbbbbbbbbbbb gtr aaaaaaaaaaaaaaa 是成立的，换成数字，只要长度相同就会得到对的结果

42行是我用来校验结果是否正确的，调用 Perl 显示对应的开根结果。要正确执行需要安装 Perl 环境。
http://strawberryperl.com/releases.html 这玩意儿好像都不自带环境配置，装完还要自己把perl.exe执行路径添加到PATH中

对应的校验语句是临时注释掉的：
rem call :check_first %1 !precision!

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
复制代码

523066680

回复 38# SQYSQYSQY

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

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

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

523066680

回复 40# SQYSQYSQY

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

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

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


浮点数版本补充在39楼。

523066680

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

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

回复 40# SQYSQYSQY

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

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

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

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

用了两天 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.

老刘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
复制代码

老刘1号

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

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号

回复 46# SQYSQYSQY
回复 47# 523066680


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

523066680

回复 48# 老刘1号


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

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