long double for FB-Win64

Windows specific questions.
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: long double for FB-Win64

Post by dodicat »

The Chebychev points are supposed to minimise maximum errors (minimax theorems e.t.c.)
Splitting the range equally would hardly make any difference, in fact as long as you have X amounts of points reasonably spaced the resulting polynomial coefficients should be OK.
srvaldez
Posts: 3373
Joined: Sep 25, 2005 21:54

Re: long double for FB-Win64

Post by srvaldez »

thank you dodicat for participating :-)
I wanted to know whether my latest version of clongdouble would work for you, in some cases, it can be significantly faster than double, though in most cases it's slower.
what I found interesting is how efficiently gcc implements the math functions in x64, the parameters are passed by reference with the first parameter for the result, a naked sub is ideal for this and very efficient, needing a simple jump.
on the other hand, I need to read and study about Chebyshev series.
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: long double for FB-Win64

Post by dodicat »

The Chebyshev points are the roots of functions of the form
cos(N*acos(x))
So if you want 5 points (for the first 5 coefficients) the points would be the roots of cos(5*acos(x))
The function cos(N*acos(x)) only exists between -1 and 1 (because of the acos bit)
and it ranges -1 and 1 vertically (because of the cos bit)
Here is the slow way of getting the roots for Chebyshev 3, 4 and 5.

Code: Select all



#include "crt.bi"  

Namespace globals
Dim Shared As Integer xres,yres
Dim Shared As Double minx,maxx,miny,maxy,PLOT_GRADE=5000
Dim Shared As Double MinimumY,MaximumY
Dim Shared As Double MinimumX,MaximumX
Type fun  As Function(x As Double) As Double
Dim Shared f As fun
End Namespace

Sub sketch(fn As Any Ptr,colour As Ulong,axiscolour As Ulong=Rgb(150,150,150))
    Using globals
    f=fn
    Dim As Double last=f(minx)
    For x As Double=minx To maxx Step (maxx-minx)/PLOT_GRADE
        Dim As Double x1=(xres)*(x-minx)/(maxx-minx)
        Dim As Double d=f(x)
        Dim As Double y1=(yres)*(d-maxy)/(miny-maxy)
        If Sgn(last)<> Sgn(d)  Then Circle(x1,y1),2,0,,,,f
        If x=minx Then Pset(x1,y1),colour Else Line -(x1,y1),colour
        last=d
    Next x
    'axis
    Dim As Long f1,f2
    If Sgn(minx)<>Sgn(maxx) Then
        Line(((minx/(minx-maxx))*xres),0)-(((minx/(minx-maxx))*xres),yres),(axiscolour) 'y axis
        f1=1
        If Sgn(minx)=0 Or Sgn(maxx)=0 Then f1=0
    End If
    If Sgn(miny)<>Sgn(maxy) Then
        Line(0,(yres-(miny/(miny-maxy))*yres))-(xres,(yres-(miny/(miny-maxy))*yres)),(axiscolour) 'x axis
        f2=1
        If Sgn(miny)=0 Or Sgn(maxy)=0 Then f2=0
    End If
    
    If f2 Then
        Draw String(0,(yres-(miny/(miny-maxy))*yres)),Str(minx),(axiscolour)
        Draw String(xres-8-8*(Len(Str(maxx))),(yres-(miny/(miny-maxy))*yres)),Str(maxx),(axiscolour)
    Else
        Draw String(0,yres/2),Str(minx),(axiscolour)
        Draw String(xres-8-8*(Len(Str(maxx))),yres/2),Str(maxx),(axiscolour)
    End If
    
    If f1 Then
        Draw String(((minx/(minx-maxx))*xres),0),Str(maxy),(axiscolour)
        Draw String(((minx/(minx-maxx))*xres),yres-16),Str(miny),(axiscolour)
    Else
        Draw String(xres/2,0),Str(maxy),(axiscolour)
        Draw String(xres/2,yres-16),Str(miny),(axiscolour)
    End If
End Sub

Sub getyrange(fn As Any Ptr,sx As Double,lx As Double,Byref by As Double,Byref sy As Double)
    Using globals
    f=fn
    #macro _window(topleftX,topleftY,bottomrightX,bottomrightY) 
    minx=(topleftX)
    maxx=(bottomrightX)
    miny=(bottomrightY)
    maxy=(topleftY)
    #endmacro
    MinimumY=1e50:MaximumY=-1e50
    For n As Double=MinimumX To lx Step(lx-MinimumX)/10000
        Dim As Double v=f(n)
        If MinimumY>V Then MinimumY=v
        If MaximumY<V Then MaximumY=V
    Next
    _window(MinimumX,MaximumY,MaximumX,MinimumY) 
End Sub

Sub bisect(fn As Any Ptr,min As Double,max As Double,Byref O As Double)
    Using globals
    f=fn
    Dim As Double last,st=(max-min)/100000,v
    For n As Double=min To max Step st
        v=f(n)
        If Sgn(v)<>Sgn(last) Then 
            Puts "Root "+str(n)+string(60-len(str(n))," ")+"Error =  "+str(f(n))
            O=n+st:Exit Sub
        End If
        last=v
    Next
End Sub

Sub roots(fn As Any Ptr,min As Double,max As Double)
    Using globals
    f=fn
    MinimumX=min
    MaximumX=max
    Dim As Double last,O,v,st=(max-min)/10000000
    For n As Double=min To max Step st
        v=f(n)
        If Sgn(v)<>Sgn(last) And n>min Then bisect(f,n-st,n,O):n=O
        last=v
    Next
    ''  screen plot  -- get fn moving
    getyrange(f,MinimumX,MaximumX,MinimumY,MaximumY)
    cls
    Color ,Rgb(236,233,216)
    Screeninfo globals.xres,globals.yres
    Screencontrol 100,.4*globals.xres,.4*globals.yres
    Cls
     line(0,0)-(globals.xres-1,globals.yres-1),rgb(0,0,0),b
    sketch(f,Rgb(0,100,255))
End Sub

'======================================  

#macro InputFunction(n,fn,minx,maxx)
Puts #fn +"     " +str(minx) + "   to   "+str(maxx)
puts " "
windowtitle  #fn
Function f##n(x As Double) As Double
    Return fn
End Function
roots(Procptr(f##n),minx,maxx)     ' Please note: catches roots AND discontinuaties in the given x range
Puts "---------------------"
#endmacro

'================   USE ======================
'Note InputFunction parameter 1 must be unique at each run eg 1  2  3  ...
screen 19,32
InputFunction(1,cos(3*acos(x)),-1,1)  'chebyshev functions
puts "Press a key"
sleep

InputFunction(2,cos(4*acos(x)),-1,1) 
puts "Press a key"

sleep
InputFunction(3,cos(5*acos(x)),-1,1)
puts "Done"
sleep


 
Chebyshev did all the maths of course, he reckoned if, say for four points, instead of using the lowest, then 1/3 up then 2/3 up then the last point, you should choose the roots of cos(4*acos(x)) to get the best set of coefficients.
If the range is not -1 to 1 then his -1 to 1 points are mapped to the new range.
His actual proof is complicated to follow, (it was when I was young, I am scared to look at it now)
But anyway that is the gist.
I think it is used by Intel in their cpu sin, cos, tan approximations.
srvaldez
Posts: 3373
Joined: Sep 25, 2005 21:54

Re: long double for FB-Win64

Post by srvaldez »

thank you dodicat, yes Chebyshev series often give near minimax approximations and are relatively easy to calculate.
[edit]
in case anyone reads this post, I have completely redone the clongdouble bi files using intel syntax for the inline asm, also, instead of using the ioldouble routines for input/output I use the mingw sscanf and sprintf functions, it gives you all the flexibility of the C float formatting options.
also, there's a win32 version, all in the first 3 posts.
Post Reply