Dos Batch Math Library

[Aacini]

. . .

Note:
I have not tested, if this algorithm computes always the correct result.
Some values within [0:2^31-1] may need more iterations.
Because of the starting value the number of iterations is not monotonically increasing anymore
(if i remember right then the original algorithm (initial x=N) needs a maximum of 19 iterations).

But you should check this using c++/java/similar... batch is much too slow to check all numbers within [0:2^32-1];
my pc would take >300 days, but it isn't the fastest... .

penpen

Well, that is not really neccessary... The problems with these iterative methods appears at the limits of a range of results. For example, the Sqrt from 25 to 35 is 5, so it is enough to test the method with 25 and 35 to know that all values in this range will return the same result. This way, the test should be done on the square of a number and this value minus one for all numbers from 2 up to Sqrt(2^31-1), that is 46340. I wrote the program below to perform this test over both methods:

Code: Select all

@echo off
setlocal EnableDelayedExpansion

for /F %%a in ('copy /Z "%~F0" NUL') do set "CR=%%a"


rem Aacini's method:
set "Sqrt1(Num)=set /A "M=(Num),sqrt=M/(11*1024)+40,sqrt=(M/sqrt+sqrt)/2"&(for /L %%i in (1,1,5) do set /A "x2=(M/sqrt+sqrt)/2"^&if ^!x2^! leq ^!sqrt^! set /A sqrt=x2)"
echo Aacini start @ %time%

< NUL (for /L %%n in (2,1,46340) do (

   set /P "=%%n!CR!"

   set /A "Num2=%%n*%%n, Num2M1=Num2-1, NM1=%%n-1"

   %Sqrt1(Num):Num=Num2M1%
   if !sqrt! neq !NM1! echo Aacini failed on %%n- sqrt(!Num2M1!^) = !sqrt!
   %Sqrt1(Num):Num=Num2%
   if !sqrt! neq %%n echo Aacini failed on %%n- sqrt(!Num2!^) = !sqrt!

))

echo Aacini end   @ %time%


echo =====================


rem penpen's method:
set "Sqrt2(N)=( x=(N)/(11*1024)+40-^!(N)*9, x=((N)/x+x)/2, x=((N)/x+x)/2, x=((N)/x+x)/2, x=((N)/x+x)/2, x=((N)/x+x)/2, x+=(((N)-((x-1)*(x-1)+(x+x-1)))>>31))/(1+((N)>>31))"
echo penpen start @ %time%

< NUL (for /L %%n in (2,1,46340) do (

   set /P "=%%n!CR!"

   set /A "Num2=%%n*%%n, Num2M1=Num2-1, NM1=%%n-1"

   set /A "sqrt=!Sqrt2(N):N=Num2M1!"
   if !sqrt! neq !NM1! echo penpen failed on %%n- sqrt(!Num2M1!^) = !sqrt!
   set /A "sqrt=!Sqrt2(N):N=Num2!"
   if !sqrt! neq %%n echo penpen failed on %%n- sqrt(!Num2!^) = !sqrt!

))

echo penpen end   @ %time%

Each test takes less than 5 minutes; this is the output:

Code: Select all

Aacini start @ 23:50:44.84
Aacini end   @ 23:55:04.78
=====================
penpen start @ 23:55:04.79
penpen end   @ 23:57:23.93

My method is slower because it have a FOR command; I will try to eliminate it in order to speed up the method.

Antonio

[Aacini]

Ok. I converted the FOR loop into a long expression; after that, I realized that the IF test may be performed just in the last iteration and the result is the same. I also slightly optimized the expression changing the division by 2 for a right shift.

Code: Select all

@echo off
setlocal EnableDelayedExpansion

for /F %%a in ('copy /Z "%~F0" NUL') do set "CR=%%a"

rem Aacini's method, converted into a long expression:
set "Sqrt(N)=( M=(N),x=M/(11*1024)+40, x=(M/x+x)>>1, x=(M/x+x)>>1, x=(M/x+x)>>1, x=(M/x+x)>>1, x=(M/x+x)>>1, x2=(M/x+x)>>1,x+=-(x2-x-1>>31)*(x2-x) )"

echo Aacini2 start @ %time%

< NUL (for /L %%n in (2,1,46340) do (

   set /P "=%%n!CR!"

   set /A "Num2=%%n*%%n, Num2M1=Num2-1, NM1=%%n-1"

   set /A "sqrt=%Sqrt(N):N=Num2M1%"
   if !sqrt! neq !NM1! echo Aacini failed on %%n- sqrt(!Num2M1!^) = !sqrt!
   set /A "sqrt=%Sqrt(N):N=Num2%"
   if !sqrt! neq %%n echo Aacini failed on %%n- sqrt(!Num2!^) = !sqrt!

))

echo Aacini2 end   @ %time%

This is the output:

Code: Select all

Aacini2 start @ 11:50:33.11
Aacini2 end   @ 11:52:53.84

The timing in the original FOR-loop method was 4:19.94 minutes but this one takes just 2:20.73 minutes, practically the same than penpen's version that was 2:19.14 minutes.

Antonio

[Squashman]

56 seconds on my 7 year old Duo Core. My gosh. How old are the computers you guys are using?

[einstein1969]

@Aacini

Yes. Your method fails with numbers as less as 35 or 48. The method iterates on a series of values and in each step it must produce a value less or equal than the previous one, so the method may fail when a next value is greater. I inserted the required IF command, but then realized that the test must not be applied to the first value. I inserted the test from the second value on and cut the loop to 4 cycles, but then realized that in certain cases the result was wrong and it requires an additional step. This means that it needs a total of 6 cycles.

Thanks, now it clear. Don't worry for the 6 iterations... We optimize later when all work.

@penpen

So in your algorithm (using int calculus) you could test if (N-x*x) is greater than or equal to 0;
if this is true then sqrt(N)=x, else sqrt(N)=x-1.

It's the same solution proposed by AAcini?

set "Sqrt(N)=( x=(N)/(11*1024)+40-^!(N)*9, x=((N)/x+x)/2, x=((N)/x+x)/2, x=((N)/x+x)/2, x=((N)/x+x)/2, x=((N)/x+x)/2, x+=(((N)-((x-1)*(x-1)+(x+x-1)))>>31))/(1+((N)>>31))"

where end the ) of the block of separator "," ?

Thanks for this solution!

@Aacini
Very nice the Test. I like it!!!
This is necessary.

set "Sqrt(N)=( M=(N),x=M/(11*1024)+40, x=(M/x+x)>>1, x=(M/x+x)>>1, x=(M/x+x)>>1, x=(M/x+x)>>1, x=(M/x+x)>>1, x2=(M/x+x)>>1,x+=-(x2-x-1>>31)*(x2-x) )"

In this moment I don't undestand the difference with penpen method...
I will study better.

@Squashman
For better test empty the environment.

Code: Select all

for /f "delims==" %%v in ('set') do set "%%v="

If there are open applications the time is very different. I have closed and paused the applications that used more CPU.
On my celeron @2.00Ghz (The only pc live and this is a present... The money are finished ) take 19 second with empty environment e ReM the set/P

EDIT:
Not bad for about 90000 sqrt's

Code: Select all

Aacini2 start @ 22:48:12,25
Aacini2 end   @ 22:48:31,13

Einstein1969

[penpen]

@penpen

So in your algorithm (using int calculus) you could test if (N-x*x) is greater than or equal to 0;
if this is true then sqrt(N)=x, else sqrt(N)=x-1.

It's the same solution proposed by AAcini?

I think this is not the case, but i might have missed it:
As far as i understand Aacini's algorithm, it performs 6 iterations and chooses the minimal one of the last 2 iterations.
My algorithm performs 5 iterations, and then checks if (N > x*x) and then returns x, else it returns (x-1).
(The last step equals one iteration if needed, because in this case (signed int 32 bit) ||x-N/x||<=2.)

set "Sqrt(N)=( x=(N)/(11*1024)+40-^!(N)*9, x=((N)/x+x)/2, x=((N)/x+x)/2, x=((N)/x+x)/2, x=((N)/x+x)/2, x=((N)/x+x)/2, x+=(((N)-((x-1)*(x-1)+(x+x-1)))>>31))/(1+((N)>>31))"

where end the ) of the block of separator "," ?

It is hard to see, and it simply was a flaw (it was late when i wrote it):
set "Sqrt(N)=( x=(N)/(11*1024)+40-^!(N)*9, x=((N)/x+x)/2, x=((N)/x+x)/2, x=((N)/x+x)/2, x=((N)/x+x)/2, x=((N)/x+x)/2, x+=(((N)-((x-1)*(x-1)+(x+x-1)))>>31))/(1+((N)>>31))"

If you want, then you could place this ) to the end of the formula.


penpen

[Squashman]

56 seconds on my 7 year old Duo Core. My gosh. How old are the computers you guys are using?

27 seconds on my Quad Core at home. My Duo Core is the piece of crap that work gave me when I needed a laptop for a new job. I went from a quad core I5 desktop at work to a old shitty Duo Core. I told my boss that I would not take the job unless he gave me a laptop the same speed. Should of had that written into my contract.

[Aacini]

This is an excellent example of the reasons because I always give great importance to clear descriptions.

First penpen's explanation was a mathematical one. When I saw it I have not much spare time, so I postponed a more detailed study of it for later. However, the phrase in his last explanation is crystal clear: "My algorithm performs 5 iterations, and then checks if (N > x*x) and then returns x, else it returns (x-1)". I modified my last code using this description; the result is a method 11% faster:

Code: Select all

@echo off
setlocal EnableDelayedExpansion

for /F %%a in ('copy /Z "%~F0" NUL') do set "CR=%%a"

rem Aacini's method, optimized using this penpen's description:
rem "Perform 5 iterations; if x*x > N decrement x by 1"
set "Sqrt(N)=( M=(N),x=M/(11*1024)+40, x=(M/x+x)>>1, x=(M/x+x)>>1, x=(M/x+x)>>1, x=(M/x+x)>>1, x=(M/x+x)>>1, x+=(M-x*x)>>31 )"

echo Aacini3 start @ %time%

< NUL (for /L %%n in (2,1,46340) do (

   set /P "=%%n!CR!"

   set /A "Num2=%%n*%%n, Num2M1=Num2-1, NM1=%%n-1"

   set /A "sqrt=%Sqrt(N):N=Num2M1%"
   if !sqrt! neq !NM1! echo Aacini failed on %%n- sqrt(!Num2M1!^) = !sqrt!
   set /A "sqrt=%Sqrt(N):N=Num2%"
   if !sqrt! neq %%n echo Aacini failed on %%n- sqrt(!Num2!^) = !sqrt!

))

echo Aacini3 end   @ %time%

Output:

Code: Select all

Aacini3 start @ 21:21:31.41
Aacini3 end   @ 21:23:36.12

= 2:04.71 minutes

@Squashman: My computer is not an old one, but a very cheap one!

Antonio

[Squashman]

the result is a method 11% faster:

Yep! Went down to 46 seconds on my Duo Core. Ten seconds faster!

[Ed Dyreen]

56 seconds on my 7 year old Duo Core. My gosh. How old are the computers you guys are using?

My virtual Machine is currently running on a dual socket hyper threading Xeon server of about 10 years old that I bought for 30€.

I am planning to upgrade it to a quad socket single core Xeon Server with 3PSU's of 400 watts and room for 16 DRR1 slots and room for 16 scsi drives of about 12 years old that I bought for 50€.

That's older than my nefue cuz he's only 7

I am not a fan of single sockets no matter how many cores as their speed is always limited to the bandwidth of a single socket and that means sharing the limited bandwidth of a single socket between cores, no thanks. Due to the separate pipelines a multi socket can access the memory simultaneously as long as they are not both trying to write to the same address, can't do that with a single socket.

[penpen]

I didn't use "N-x*x", but "N-(x-1)*(x-1)-2*x+1" because i suspected the first term to fail for x>46340:
I obviously was wrong on that part.

Note that i have added the term "/(1+((N)>>31))" to guard against negative N values.
The most speed you gained would be lost if you add it (again; which i'm in favour of).


penpen

[einstein1969]

Another problem on percent expand. I probed to undenstand but i don't have good basic
Why?

Code: Select all

@echo off
setlocal EnableDelayedExpansion

rem Aacini
rem As far as i understand Aacini's algorithm, it performs 6 iterations and chooses the minimal one of the last 2 iterations.
set "Sqrt1(N)=( M=(N),x=M/(11*1024)+40, x=(M/x+x)>>1, x=(M/x+x)>>1, x=(M/x+x)>>1, x=(M/x+x)>>1, x=(M/x+x)>>1, x2=(M/x+x)>>1,x+=-(x2-x-1>>31)*(x2-x) )"


rem penpen
rem My algorithm performs 5 iterations, and then checks if (N > x*x) and then returns x, else it returns (x-1).
Rem (The last step equals one iteration if needed, because in this case (signed int 32 bit) ||x-N/x||<=2.)

set "Sqrt2(N)=( x=(N)/(11*1024)+40-^!(N)*9, x=((N)/x+x)/2, x=((N)/x+x)/2, x=((N)/x+x)/2, x=((N)/x+x)/2, x=((N)/x+x)/2, x+=(((N)-((x-1)*(x-1)+(x+x-1)))>>31))/(1+((N)>>31))"


(
  set /A "sqrt1=%Sqrt1(N):N=36%"
)
(
  set /A "sqrt2=%Sqrt2(N):N=36%"
)
echo %sqrt1% %sqrt2%
pause

result:

Code: Select all

6 -11

einstein1969

[penpen]

The delayed expansion is enabled, so this is the expansion:

Code: Select all

:: source %Sqrt2(N):N=36%
:: initial term
x=(N)/(11*1024)+40-^!(N)*9
:: %-expansion
x=(36)/(11*1024)+40-!(36)*9
:: !-expansion
x=(36)/(11*1024)+40-(36)*9
:: computes to
x=-284

The exclamation mark should have been used as the not operator:

Code: Select all

:: source !Sqrt2(N):N=36!
:: initial term
x=(N)/(11*1024)+40-^!(N)*9
:: no %-expansion
:: !-expansion
x=(36)/(11*1024)+40-!(36)*9
:: computes to
x=40

Note: Aacini has demonstrated that the "-^!(N)*9" part isn't neccessary anymore.


penpen

[einstein1969]

Ok. This is second step for optimizing the sqrt.

I developed few days ago and there is the precedent problem of oscillating series of Newton. I have not applicated the solution of penpen and Aacini for the moment.

What is the question that i raise?

When perform a sqrt the integer/integral part is performed by the precedent developed code Integer 32 BIT SQRT. But for normal application is necessary some time the fractional part.

In the previus developed code I, dbenham and trebor68 have used the classic schoolar algorithm with good result in term of precision and velocity, but I go down for a very fast solution. In this mode the 32bit of the result can be full of the fractional part.

In this i need your help.

I introduce a my understud of FIXED POINT mathematic.

For example the SQRT result can be of the type 5.4. 5 integral digit plus 4 fractional digit.
And the result can be stored in TWO variables or in the same variable at BIT/DECIMAL position different.

This is an example of start to compute a sqrt with fractional part using the INTEGER 32bit Sqrt

Code: Select all

@echo off
setlocal EnableDelayedExpansion

Rem Fixed Point
Rem If you have a fractional number, than precalc multiply by EVEN power of 10.
Rem es. SQRT(1.2)=SQRT(1.2*100)/10=SQRT(120)/10

rem define sqrt
set "Sqrt(N)=( x=(N)/(11*1024)+40-^!N*9, x=((N)/x+x)/2, x=((N)/x+x)/2, x=((N)/x+x)/2, x=((N)/x+x)/2, x=((N)/x+x)/2, x/=(1+((N)>>31)))"

rem es. sqrt 25 with 1 decimal
set /A "num=!Sqrt(N):N=25*100!"
echo Sqrt(25.0)=%num:~0,-1%.%num:~-1% & echo(

rem es. sqrt 28.0
set /A "num=!Sqrt(N):N=28*100!"
echo Sqrt(28.0)=%num:~0,-1%.%num:~-1% & echo(

rem es sqrt 2.0000 (with 4 decimal)
set /A "num=!Sqrt(N):N=2*100*100*100*100!"
echo Sqrt(2.0000)=%num:~0,-4%.%num:~-4% & echo(

result:

Code: Select all

Sqrt(25.0)=5.0

Sqrt(28.0)=5.2

Sqrt(2.0000)=1.4142

This is limitated because the 32BIT. And the max number of digit rappresentable is 123.4567890 or 12345678.90 and all other combination. EDIT: (my flame:THIS is not TRUE)

The best method is get the REMAINDER and continue from there.

I have developed in this direction. I get the integer result of sqrt and the Remainder
and I have probed to applicate the Newton iteration!

I introduce the simple mathematic (I'm not an expert...)
You can see this example on StackOverflow But there isn't the Newton implementation. I done

Call the integer sqrt "ISQRT" and "SQRT" is the real/floating point.

ISQRT(140)=11

140=11^2+19 (19 is the REMAINDER = N-ISQRT(N)*ISQRT(N)) [1]

The SQRT(140) = 11,832 that is (11+0.832)

that is (11+0.832)^2=140 [2]

we need found 0.832 and then rewrite the [2]:

(11+x)^2=140

expand this : 11^2+2*11*x+x^2=140

but if we get the [1] we see that 2*11*x+x^2=19

With a guess of 0.5 (0<x<1) iterating via newton we get the fractional part in about 3 iteration.

That is:

Code: Select all

set "Sqrt(N)=( x=(N)/(11*1024)+40-^!N*9, x=((N)/x+x)/2, x=((N)/x+x)/2, x=((N)/x+x)/2, x=((N)/x+x)/2, x=((N)/x+x)/2, Sqrt_Remai_der=(N)-x*x, x=x )"


rem approx via newton
rem es. sqrt of 60000
set /A "num=!Sqrt(N):N=60000!"
echo Sqrt(60000)=%num% Remainder=%Sqrt_Remai_der% & echo(

echo 2*%num%*x+x*x-%Sqrt_Remai_der%=0
rem set /A approx=0.5-(0.25+2*num/2-Sqrt_Remai_der)/(2*0.5+2*num) & rem guess=0.5
rem set /A approx=0.5*1000-(0.25*1000+2*num/2*1000-Sqrt_Remai_der*1000)/(2*0.5+2*num) & rem guess=0.5
set /A approx=500-(250+2*num/2*1000-Sqrt_Remai_der*1000)/(1+2*num)
For /L %%I in (1,1,2) do set /A approx=approx-(approx*approx/1000+2*num*approx-Sqrt_Remai_der*1000)*1/(2*approx/1000+2*num)
echo Sqrt(60000)=%num%.%approx%  approx.  (244.948...) & echo(

rem es. sqrt of 2
set /A "num=!Sqrt(N):N=2!"
echo Sqrt(2)=%num% Remainder=%Sqrt_Remai_der% & echo(

echo 2*%num%*x+x*x-%Sqrt_Remai_der%=0
set /A approx=500-(250+2*num/2*1000-Sqrt_Remai_der*1000)/(1+2*num)
For /L %%I in (1,1,2) do set /A approx=approx-(approx*approx/1000+2*num*approx-Sqrt_Remai_der*1000)*1/(2*approx/1000+2*num)
echo Sqrt(2)=%num%.%approx%    (1.414...) & echo(

rem es. sqrt of 3
set /A "num=!Sqrt(N):N=3!"
echo Sqrt(3)=%num% Remainder=%Sqrt_Remai_der% & echo(

echo 2*%num%*x+x*x-%Sqrt_Remai_der%=0
set /A approx=500-(250+2*num/2*1000-Sqrt_Remai_der*1000)/(1+2*num)
For /L %%I in (1,1,2) do set /A approx=approx-(approx*approx/1000+2*num*approx-Sqrt_Remai_der*1000)*1/(2*approx/1000+2*num)
echo Sqrt(3)=%num%.%approx%    (1.732...) & echo(

rem es. sqrt of 2147483647
set /A "num=!Sqrt(N):N=2147483647!"
echo Sqrt(2147483647)=%num% Remainder=%Sqrt_Remai_der% & echo(

echo 2*%num%*x+x*x-%Sqrt_Remai_der%=0
set /A approx=500-(250+2*num/2*1000-Sqrt_Remai_der*1000)/(1+2*num)
For /L %%I in (1,1,2) do set /A approx=approx-(approx*approx/1000+2*num*approx-Sqrt_Remai_der*1000)*1/(2*approx/1000+2*num)
echo Sqrt(2147483647)=%num%.%approx%    (46340.950...) & echo(

rem es. sqrt of 65535
set /A "num=!Sqrt(N):N=65535!"
echo Sqrt(65535)=%num% Remainder=%Sqrt_Remai_der% & echo(

echo 2*%num%*x+x*x-%Sqrt_Remai_der%=0
set /A approx=500-(250+2*num/2*1000-Sqrt_Remai_der*1000)/(1+2*num)
For /L %%I in (1,1,2) do set /A approx=approx-(approx*approx/1000+2*num*approx-Sqrt_Remai_der*1000)*1/(2*approx/1000+2*num)
echo Sqrt(65535)=%num%.%approx%    (255.998...) & echo(

Rem Patching

And we get the error at 65535 with the REMAINDER=-1

Code: Select all

Sqrt(2)=1 Remainder=1

2*1*x+x*x-1=0
Sqrt(2)=1.414    (1.414...)

Sqrt(3)=1 Remainder=2

2*1*x+x*x-2=0
Sqrt(3)=1.732    (1.732...)

Sqrt(2147483647)=46340 Remainder=88047

2*46340*x+x*x-88047=0
Sqrt(2147483647)=46340.950    (46340.950...)

Sqrt(65535)=256 Remainder=-1

2*256*x+x*x--1=0
Sqrt(65535)=256.-1    (255.998...)

I have not patched at the moment. I need a confirm by you if the error is the same

@penpen
I can't understund how to patch the ! problem in the percent expansion... I will go down the iteration from 5 to 4 or 3

EDIT: added the definition of ISQRT with remainder. Not optimized.
EDIT2: correct explain on ISQRT(140)

Einstein1969

[Ed Dyreen]

The exclamation mark should have been used as the not operator:

Bizarre, it is missing from the set /? help. from a simple test I assume that ! only executes if the operand NEQ the operated

Code: Select all

C:\profSys\ADMIN>set /a x=(36)-(35)
1
C:\profSys\ADMIN>set /a x=(36)-!(35)
36

Is this a correct interpretation and why is it not listed under set /? while the operator & clearly is ??

[penpen]

@penpen
I can't understund how to patch the ! problem in the percent expansion... I will go down the iteration from 5 to 4 or 3

If you want to create a function, that works with both "%"-expansion and "!"-expansion, then i think:
This is impossible.

You would have to escape the exclamation mark "^^^!" if you want the correct result using "%sqrt2(N):N=36%";
but then you lose the correct function of "!sqrt2(N):N=36!" because the tilde character is interpreted as XOR.


Sidenote:
If you want an sqrt with as many decimal digits as possible, you should use another algorithm.
I cannot remember the name, but maybe i could find a link.
(This algorithm probably is slower, but computes a digit every iteration.)

Bizarre, it is missing from the set /? help.

It is present (but not explained) in the technet page listed as unary operator:
If N equals 0 then !N equals 1, else 0.


penpen