square root a la lanczos

General FreeBASIC programming questions.
Post Reply
srvaldez
Posts: 3373
Joined: Sep 25, 2005 21:54

square root a la lanczos

Post by srvaldez »

just for fun to see if it was possible to approximate functions other than the gamma function using the Lanczos scheme
to illustrate I will use a system of linear equations of degree 5, this will not give us much precision but it will show the algorithm
lets approximate sqr between 1 and 2, here I will divide the interval between 0 and 1
matrix A

Code: Select all

[[1, 1/(0+1), 1/(0+2), 1/(0+3), 1/(0+4)],
 [1, 1/(1/4+1), 1/(1/4+2), 1/(1/4+3), 1/(1/4+4)],
 [1, 1/(2/4+1), 1/(2/4+2), 1/(2/4+3), 1/(2/4+4)],
 [1, 1/(3/4+1), 1/(3/4+2), 1/(3/4+3), 1/(3/4+4)],
 [1, 1/(4/4+1), 1/(4/4+2), 1/(4/4+3), 1/(4/4+4)]]
 
matrix B

Code: Select all

[[sqr(1+0)],
 [sqr(1+1/4)],
 [sqr(1+2/4)],
 [sqr(1+3/4)],
 [sqr(1+4/4)]]
 
solve A=B

Code: Select all

Function Sqr1(Byval z As Double) As Double
	Return 4.1699272354475147801+_
			.14280730704570221521/(z+1.)-_
			7.8444468921392051374/(z+2.)+_
			41.993569124718655109/(z+3.)-_
			53.553469885319331185/(z+4.)
End Function

Dim As Double y
Print "  x           sqrt(1+x)                error"
For x As Double=0 To 1.05 Step 0.1
	y=Sqr1(x)
	Print Using "##.###";x;
	Print "  ";y, y-sqr(1+x)
Next
sleep
output

Code: Select all

  x           sqrt(1+x)                error
 0.000   1.000000000000002   1.776356839400251e-015
 0.100   1.048791856644213  -1.699152593870323e-005
 0.200   1.095439834710485  -5.2802998475876e-006
 0.300   1.140178581579455   3.156480317345611e-006
 0.400   1.183219503184455   3.546564532053154e-006
 0.500   1.224744871391589   4.440892098500626e-016
 0.600   1.264908485536255  -2.578531096864012e-006
 0.700   1.303838815119589  -1.66592094119622e-006
 0.800   1.341642803058836   2.016558961814852e-006
 0.900   1.378409547471753   4.67226273070942e-006
 1.000   1.414213562373096   6.661338147750939e-016
 
the coefficients can be transformed and the function rewritten as

Code: Select all

Function Sqr1(Byval z As Double) As Double
	dim as double n, m
	n=24+(61.986659812467627+_
		(57.111308428075295+_
		(22.437732008780969+_
		4.1699272354475148*z)*z)*z)*z
	m=24+(50+(35+(10+z)*z)*z)*z
	return n/m
End Function
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: square root a la lanczos

Post by dodicat »

From the LolRemez thread
viewtopic.php?f=7&t=25897&p=275050&hili ... ly#p275050

Code: Select all




function sqrt(x as double) as double
' approximating sqr(x)  [1 To 2]
Return _ 
((((((+0.00488320646279881)*x-0.04697287955202312)*x+0.1946001357427462)*x-0.4743301141755174)*x _ 
+1.028216061109412)*x _ 
+0.2936052298143235)
end function

Dim As Double y
Print "  x           sqrt(1+x)                error"
For x As Double=0 To 1.05 Step 0.1
   y=sqrt(1+x)
   Print Using "##.###";x;
   Print "  ";y, y-sqr(1+x)
Next
sleep

 
srvaldez
Posts: 3373
Joined: Sep 25, 2005 21:54

Re: square root a la lanczos

Post by srvaldez »

Hi dodicat
minimax approximations are better, I presented my scheme as a novelty, I don't think anyone else would be crazy enough to come up with it
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: square root a la lanczos

Post by dodicat »

My solution was minimax.
Here is another expanded version, minimax via Chebyshev.

Code: Select all



Sub GaussJordan(matrix() As Double,rhs() As Double,ans() As Double)
      Dim As Long n=Ubound(matrix,1)
      Redim ans(0):Redim ans(1 To n)
      Dim As Double b(1 To n,1 To n),r(1 To n)
      For c As Long=1 To n 'take copies
            r(c)=rhs(c)
            For d As Long=1 To n
                  b(c,d)=matrix(c,d)
            Next d
      Next c
      #macro pivot(num)
      For p1 As Long  = num To n - 1
            For p2 As Long  = p1 + 1 To n  
                  If Abs(b(p1,num))<Abs(b(p2,num)) Then
                        Swap r(p1),r(p2)
                        For g As Long=1 To n
                              Swap b(p1,g),b(p2,g)
                        Next g
                  End If
            Next p2
      Next p1
      #endmacro
      For k As Long=1 To n-1
            pivot(k)              'full pivoting each run
            For row As Long =k To n-1
                  If b(row+1,k)=0 Then Exit For
                  Var f=b(k,k)/b(row+1,k)
                  r(row+1)=r(row+1)*f-r(k)
                  For g As Long=1 To n
                        b((row+1),g)=b((row+1),g)*f-b(k,g)
                  Next g
            Next row
      Next k
      'back substitute 
      For z As Long=n To 1 Step -1
            ans(z)=r(z)/b(z,z)
            For j As Long = n To z+1 Step -1
                  ans(z)=ans(z)-(b(z,j)*ans(j)/b(z,z))
            Next j
      Next    z
End Sub

'Interpolate through point via a polynomial (spline)
Sub Interpolate(x_values() As Double,y_values() As Double,p() As Double)
      Var n=Ubound(x_values)
      Redim p(0):Redim p(1 To n)
      Dim As Double matrix(1 To n,1 To n),rhs(1 To n)
      For a As Long=1 To n
            rhs(a)=y_values(a)
            For b As Long=1 To n
                  matrix(a,b)=x_values(a)^(b-1)
            Next b
      Next a
      'Solve the linear equations
      GaussJordan(matrix(),rhs(),p())
End Sub

Function polyval(Coefficients() As Double,Byval x As Double) As Double
      Dim As Double acc
      For i As Long=Ubound(Coefficients) To Lbound(Coefficients) Step -1
            acc=acc*x+Coefficients(i)
      Next i
      Return acc
End Function

Sub chebypoints(min As Double,max As Double,order As Integer,c()As Double)
      Redim c(1 To order)
      Dim pi As Double=4*Atn(1)
      Dim count As Integer
      For k As Integer=order To 1 Step -1
            count=count+1
            c(count)= Cos((pi/2)*(2*k-1)/order)
            c(count)=(max-min)*(c(count)+1)/2+min
      Next k
      Print "chebychev points in range"
      For n As Long=Lbound(c) To Ubound(c)
            Print n,c(n)
      Next
      Print
End Sub

Function nest(pol() As Double) As String
      Dim As String xx,tmp,ch=""
      For a As Long=Ubound(pol) To Lbound(pol) Step -1
            Dim As String op
            If a=Lbound(pol) Then xx="" Else xx="*x"
            Dim As Double p=pol(a)
            If Sgn(p)>=0 Then op="+" Else op=""
            tmp="("+ch+op+Str(p)+")"+xx
            ch="("+ch+op+Str(p)+")"+xx
            If Len(tmp)>50 Then ch=ch+" _ " +Chr(10) '80
      Next a
      Return ch
End Function

Function func(x As Double) As Double
 return (((((-0.01019205450473583)*x+0.08530805745877951)*x-0.3142876923553593)*x _
+0.9127809568362455)*x _
+0.3264033259436233)
End Function

Dim As Long n=5 '<----------------------  number of terms
Dim As Double lowerx=1 '<-----------------  LIMITS 
Dim As Double upperX=2 '<------------------ LIMITS

Dim As Double x(1 To 5)
Dim As Double y(1 To 5)

Redim As Double points()
lowerx=1
upperx=2

chebypoints(lowerX,upperX,n,points())
For z As Double=1 To 5
      x(z)=points(z)
      y(z)=sqr(points(z))
Next
Redim As Double Poly(0)
interpolate(x(),y(),poly())

Print "Polynomial Coefficients:"
Print
For z As Long=1 To Ubound(Poly)
      If z=1 Then
            Print "constant term  ";Tab(20);Poly(z)
      Else
            Print Tab(8); "x^";z-1;"   ";Tab(20);Poly(z)
      End If
Next z
Print
Print "nested form:"
Print nest(poly())
Print


Print "  x           sqrt(1+x)                error"
For x As Double=0 To 1.05 Step 0.1
      Var y=func(1+x)
      Print Using "##.###";x;
      Print "  ";y, y-Sqr(1+x)
Next

Sleep



 
srvaldez
Posts: 3373
Joined: Sep 25, 2005 21:54

Re: square root a la lanczos

Post by srvaldez »

nice :-)
jj2007
Posts: 2326
Joined: Oct 23, 2016 15:28
Location: Roma, Italia
Contact:

Re: square root a la lanczos

Post by jj2007 »

dodicat wrote:From the LolRemez thread
Little speed test:

Code: Select all

... see below, 18:32 ...
Last edited by jj2007 on Mar 03, 2021 17:33, edited 2 times in total.
srvaldez
Posts: 3373
Joined: Sep 25, 2005 21:54

Re: square root a la lanczos

Post by srvaldez »

@jj2007
if you compile with optimization using gcc the non-asm version will be completely eliminated because the calculations are not used outside the timing loop, that's why I include some kind of sum in the loop and print it out after
jj2007
Posts: 2326
Joined: Oct 23, 2016 15:28
Location: Roma, Italia
Contact:

Re: square root a la lanczos

Post by jj2007 »

Good idea, srvaldez:

Code: Select all

function sqrt(x as double) as double
' approximating sqr(x)  [1 To 2]
Return _
((((((+0.00488320646279881)*x-0.04697287955202312)*x+0.1946001357427462)*x-0.4743301141755174)*x _
+1.028216061109412)*x _
+0.2936052298143235)
end function
Dim As Double y, sum
Print " x       sqrt(1+x)           error"
For x As Double=0 To 1.05 Step 0.1
  #if 1
  	asm
		fld1
		fadd QWORD ptr [x]
		fsqrt
		fstp QWORD ptr [y]
  	end asm
  #else
	y=sqrt(1+x)
  #endif
   Print Using "##.###";x;
   Print "  ";y, y-sqr(1+x)
Next
Print

Print " x       sqrt(1+x)           error"
For x As Double=0 To 1.05 Step 0.1
  #if 0
  	asm
		fld1
		fadd QWORD ptr [x]
		fsqrt
		fstp QWORD ptr [y]
  	end asm
  #else
	y=sqrt(1+x)
  #endif
   Print Using "##.###";x;
   Print "  ";y, y-sqr(1+x)
Next
Print

Dim t as double
sum=0
t=Timer()
for outerloop as long=1 to 5000000
  For x As Double=0 To 1.05 Step 0.1
	asm
		fld1
		fadd QWORD ptr [x]
		fsqrt
		fstp QWORD ptr [y]
	end asm
	sum+=y
  Next
Next

t=timer()-t
Print cast(integer, t*1000); " ms for FPU sqrt,  sum=";sum

sum=0
t=Timer()
for outerloop as long=1 to 5000000
  For x As Double=0 To 1.05 Step 0.1
	y=sqrt(1+x)
	sum+=y
  Next
Next
t=timer()-t
Print cast(integer, t*1000); " ms for the proxy, sum=";sum

sleep

Code: Select all

GAS:
 516 ms for FPU sqrt,  sum= 66977004.90197226
 688 ms for the proxy, sum= 66977016.92774981
Gcc:
 439 ms for FPU sqrt,  sum= 66977004.90197226
 307 ms for the proxy, sum= 66977016.91725335
FB64:
 436 ms for FPU sqrt,  sum= 66977004.90197226
 251 ms for the proxy, sum= 66977016.92774981
SARG:
 434 ms for FPU sqrt,  sum= 66977004.90197226
 667 ms for the proxy, sum= 66977016.92774981
srvaldez
Posts: 3373
Joined: Sep 25, 2005 21:54

Re: square root a la lanczos

Post by srvaldez »

Code: Select all

win32-gas
 175 ms for FPU sqrt,  sum= 66977004.90197226
 332 ms for the proxy, sum= 66977016.92774981
win32-gcc -O 2
 105 ms for FPU sqrt,  sum= 66977004.90197226
 135 ms for the proxy, sum= 66977016.91725335
win64-gas64
 102 ms for FPU sqrt,  sum= 66977004.90197226
 280 ms for the proxy, sum= 66977016.92774981
win64-gcc -O 2
 105 ms for FPU sqrt,  sum= 66977004.90197226
 63 ms for the proxy, sum= 66977016.92774981
 
Post Reply