recreational function approximations

[srvaldez]

this thread is not for serious applications, it's only for entertainment, in the thread viewtopic.php?f=5&t=27184 I posted a tweaked version of a very simple minimax program, I post it here with some minor modifications

Code: Select all

/'
dr. hobbs journal vol 9, number 7
"Steven A. Ruzinsky"
"A Simple Minimax Algorithm"
Jul 1984
vol 9 number 7
'/

'' uncomment one of the following
#define mytype double
'#define mytype mpfr

#If mytype=mpfr
	''=============================================================================
	#Include Once "/usr/local/include/freebasic/gmp.bi"
	#Include Once "/usr/local/include/freebasic/mpfr.bi"

	'we need to set mpfr_digits before including the overloded operators
	Dim Shared As Long mpfr_digits_precision = 62 '<- set mpfr digits-precision
	Dim Shared As mpfr_rnd_t mpfr_rounding_mode = MPFR_RNDN '<- set the rounding mode
	' rounding options are
	'  MPFR_RNDN    round to nearest, with ties to even
	'  MPFR_RNDZ    round toward zero
	'  MPFR_RNDU    round toward +Inf
	'  MPFR_RNDD    round toward -Inf
	'  MPFR_RNDA    round away from zero
	'  MPFR_RNDF    faithful rounding (not implemented yet)
	'  MPFR_RNDNA    round to nearest, with ties away from zero (mpfr_round)

	#Include Once "/usr/local/include/freebasic/mpfr-fb.bi" 'mpfr overloaded operators
	#define mpfrt 1
#endif

#define n 8
#define m 1000
#define iterations 10000'00

#define pivot 100000000000
#if mpfrt
	#define k0 mpfr(-1) '' lower limit of the range of the function
	#define kp1 mpfr(1) '' upper limit of the range of the function
	#define kpi pi_const()
	#define kpio2 kpi/2
#else
	#define k0 -1 '' lower limit of the range of the function
	#define kp1 1 '' upper limit of the range of the function
	#define kpi 3.1415926535897932384626433832795
	#define kpio2 1.5707963267948966192313216916398
#endif
dim as double elapsed=timer
common shared i as integer
common shared j as integer
common shared jmax as integer
common shared k as integer

dim shared d as mytype
dim shared e as mytype
dim shared ebest as mytype
dim shared emax as mytype
dim shared pi_half as mytype
dim shared pi as mytype
dim shared t as mytype

dim shared x(m,n) as mytype
dim shared y(m) as mytype
dim shared q(n,n) as mytype
dim shared xk(n) as mytype
dim shared qx(n) as mytype
dim shared ak(n) as mytype
dim shared a(n) as mytype 
dim shared as string fun_name

'' define fun between 1 and 3 to select the function to be approximated

#define fun 4

#if fun=0

	function func(byval x as mytype) as mytype
		return atn(x)/x
	end function

	fun_name="atn(x)/x"

	#macro powers(n)
		(2*(n)-2) 'even powers, namely 0, 2, 4 etc
	#endmacro
	
#elseif fun=1

	function func(byval x as mytype) as mytype
		return atn(x)
	end function

	fun_name="atn(x)"

	#macro powers(n)
		(2*(n)-1) 'odd powers for atn, namely 1, 3, 5 etc
	#endmacro

#elseif fun=2

	function func(byval x as mytype) as mytype
		return cos(x)
	end function

	fun_name="cos(x)"

	#macro powers(n)
		(2*(n)-2) 'even powers for cos, namely 0, 2, 4 etc
	#endmacro

#elseif fun=3

	function func(byval x as mytype) as mytype
		return sqr(1+x)
	end function

	fun_name="sqrt(1+x)"

	#macro powers(n)
		((n)-1) 'powers for sqr(1+x), namely 0, 1, 2 etc
	#endmacro

#elseif fun=4

	function func(byval x as mytype) as mytype
		return 1/tan(x)
	end function

	fun_name="1/tan(x)"

	#macro powers(n)
		(2*(n)-3) 'powers for 1/tan(x), namely -1, 1, 3 etc
	#endmacro
	
#else 'if fun is not in the range of predefined functions

	function func(byval x as mytype) as mytype
		return sin(x)
	end function

	fun_name="sin(x)"

	#macro powers(n)
		(2*(n)-1) 'odd powers for sin, namely 1, 3, 5 etc
	#endmacro

#endif

sub inc_dat
	e=y(k)
	for i=1 to n
		t=x(k,i)
		xk(i)=t
		e=e-ak(i)*t
	next
end sub

sub save_best
	for i=1 to n
		a(i)=ak(i)
	next
end sub

sub find_max
	emax=k0
	jmax=1
	for j=1 to m
		e=y(j)
		for i=1 to n
			e=e-x(j,i)*ak(i)
		next
		e=abs(e)
		if (e>emax) then
			emax=e
			jmax=j
		end if
	next

	e=y(jmax)
	for i=1 to n
		t=x(jmax,i)
		xk(i)=t
		e=e-ak(i)*t
	next
	if (emax<ebest) then
		ebest=emax
		save_best
	end if
end sub

sub printout
	dim as long p, s
	dim as  string sign, sp
	print "Coefficiennts for the function ";fun_name
	print fun_name;" ="

	for i=1 to n
		p=powers(i)
		sign=""
		sp=""
		if p<10 then sp=" "
		ak(i)=a(i)
		if i<n then s=sgn(a(i+1))
		if s<0 then
			sign="-"
		elseif s>0 then
			sign="+"
		end if
		if i=n then sign=""
		if p<>0 then
			if p=1 then
				print "x    *";abs(ak(i));" ";sign
			else
				print "x^";ltrim(str(p));sp;
				if p<0 then
					print "*";abs(ak(i));" ";sign
				else
					print " *";abs(ak(i));" ";sign
				end if
			endif
		else
			print "      ";ak(i);" ";sign
		endif
	next
	find_max
end sub

cls
pi=kpi

pi_half=kpio2
for i=0 to m
	for j=0 to n
		x(i,j)=k0
	next
next
for i=0 to n
	for j=0 to n
		q(i,j)=k0
	next
next
for i=0 to m
	y(i)=k0
next
for i=0 to n
	xk(i)=k0
	qx(i)=k0
	ak(i)=k0
	a(i)=k0
next
t=kp1
for i=1 to m
	y(i)=kp1
	e=t/m
	d=1/func(e)      ' FUNCTION IS func
	for j=1 to n
		x(i,j)=d*e^powers(j)   ' powers
	next
	t=t+1
next
for i=1 to n
	q(i,i)=k0+pivot
next
ebest=kp1
for k=1 to iterations+m
	if (k>m) then find_max else inc_dat
	d=kp1
	for j=1 to n
		t=k0
		for i=1 to n
			t=t+xk(i)*q(j,i)
		next
		qx(j)=t
		d=d+xk(j)*t
	next
	for j=1 to n
		t=qx(j)/d
		for i=1 to n
			q(i,j)=q(i,j)-qx(i)*t
		next
	next
	for j=1 to n
		t=k0
		for i=1 to n
			t=t+xk(i)*q(j,i)
		next
		ak(j)=ak(j)+t*e
	next
next
elapsed=timer-elapsed

printout
print
print "error = ";emax

'' evaluate for x=1
dim as mytype z
z=1
t=0
for i=1 to n
	t=t+a(i)
next
? "relative error for x=  1 -> ";1-t/func(z)

'' evaluate for lower limit
t=0
z = -1
for i=1 to n
	t=t+z^powers(i)*a(i)
next
? "relative error for x= -1 -> ";1-t/func(z)

'' evaluate for x = 0.5
t=0
#if mpfrt
	z = mpfr("0.5") '' z = 0.5 should be OK because 0.5 is exact in double precision
#else               '' but z = 0.1 would not be OK, so you would use: z = mpfr("0.1")
	z = 0.5
#endif
for i=1 to n
	t=t+z^powers(i)*a(i)
next
? "relative error for x=0.5 -> ";1-t/func(z)
print "elapsed time = ";elapsed;" seconds"
print "press return to end ";
sleep

[srvaldez]

just in case anyone wants this
mpfr-fb.bi

Code: Select all

Type mpfr
    Declare Constructor ( )
    Declare Constructor ( Byval rhs As Long )
    Declare Constructor ( Byval rhs As LongInt )
    Declare Constructor ( Byval rhs As Integer )
    Declare Constructor ( Byval rhs As Double )
    Declare Constructor ( Byval rhs As Single )
    Declare Constructor ( Byref rhs As String )
    Declare Constructor ( Byref rhs As mpfr )
    Declare Destructor ( )
    Declare Operator Let ( Byref rhs As mpfr )
    Declare Operator Let ( Byval rhs As Long )
    Declare Operator Let ( Byval rhs As LongInt )
    Declare Operator Let ( Byval rhs As Integer )
    Declare Operator Let ( Byval rhs As Double )
    Declare Operator Let ( Byval rhs As Single )
    Declare Operator Let ( Byref rhs As String )
    Declare Operator Cast ( ) As String
    Declare Operator Cast ( ) As Long
    Declare Operator Cast ( ) As Double
    Declare Operator Cast ( ) As Single
    
    '----------------------------------------------
    Declare Operator += (Byref rhs as mpfr)
    Declare Operator += (Byval rhs as long)
    Declare Operator += (byval rhs as double)
    Declare Operator += (byval rhs as single)
    Declare Operator += (byval rhs as longint)
    Declare Operator += (byval rhs as integer)
    Declare Operator += (byref rhs as String)
    Declare Operator -= (Byref rhs as mpfr)
    Declare Operator -= (byval rhs as double)
    Declare Operator -= (byval rhs as single)
    Declare Operator -= (byval rhs as longint)
    Declare Operator -= (byval rhs as integer)
    Declare Operator -= (byval rhs as long)
    Declare Operator -= (byref rhs as String)
    Declare Operator *= (Byref rhs as mpfr)
    Declare Operator *= (byval rhs as double)
    Declare Operator *= (byval rhs as single)
    Declare Operator *= (byval rhs as longint)
    Declare Operator *= (byval rhs as integer)
    Declare Operator *= (byval rhs as long)
    Declare Operator *= (byref rhs as String)
    Declare Operator /= (Byref rhs as mpfr)
    Declare Operator /= (byval rhs as double)
    Declare Operator /= (byval rhs as single)
    Declare Operator /= (byval rhs as longint)
    Declare Operator /= (byval rhs as integer)
    Declare Operator /= (byval rhs as long)
    Declare Operator /= (byref rhs as String)
    Declare Operator ^= (Byref rhs as mpfr)
    Declare Operator ^= (byval rhs as double)
    Declare Operator ^= (byval rhs as single)
    Declare Operator ^= (byval rhs as longint)
    Declare Operator ^= (byval rhs as integer)
    Declare Operator ^= (byval rhs as long)
    Declare Operator ^= (byref rhs as String)
	
    ' For Next Implicit step = +1
    Declare Operator For ( )
    Declare Operator Step( )
    Declare Operator Next( Byref end_cond As mpfr ) As Integer
    ' For Next Exlicit step
    Declare Operator For ( Byref stp As mpfr )
    Declare Operator Step( Byref stp As mpfr )
    Declare Operator Next( Byref end_cond As mpfr, Byref step_var As mpfr ) As Integer
    '----------------------------------------------
    declare function toString( byval frmt as string = "%32.28g") as string
    declare Function toLong ( ) As Long
    declare Function toDouble ( ) As Double
    declare Function toSingle ( ) As Single
    num As mpfr_t
End Type

Declare Operator + ( Byref lhs As mpfr, Byref rhs As mpfr ) As mpfr
Declare Operator - ( Byref lhs As mpfr, Byref rhs As mpfr ) As mpfr
Declare Operator * ( Byref lhs As mpfr, Byref rhs As mpfr ) As mpfr
Declare Operator \ ( Byref lhs As mpfr, Byref rhs As mpfr ) As mpfr
Declare Operator / ( Byref lhs As mpfr, Byref rhs As mpfr ) As mpfr
Declare Operator ^ ( Byref lhs As mpfr, Byval rhs As mpfr ) As mpfr

Declare Operator - ( Byref rhs As mpfr ) As mpfr

Declare Operator = ( Byref lhs As mpfr, Byref rhs As mpfr ) As Integer
Declare Operator < ( Byref lhs As mpfr, Byref rhs As mpfr ) As Integer
Declare Operator > ( Byref lhs As mpfr, Byref rhs As mpfr ) As Integer
Declare Operator <= ( Byref lhs As mpfr, Byref rhs As mpfr ) As Integer
Declare Operator >= ( Byref lhs As mpfr, Byref rhs As mpfr ) As Integer
Declare Operator <> ( Byref lhs As mpfr, Byref rhs As mpfr ) As Integer

function mpfr.toString( byval frmt as string = "%32.28g") as string
	dim as string sp =" "
	frmt = left(frmt,len(frmt)-1)+"R*"+right(frmt,1)
	if instr(frmt , "%")=0 then frmt = "%"+frmt
    dim as long length
	length = mpfr_snprintf(0, 0, frmt, mpfr_rounding_mode, @num)
	dim as zstring ptr rstring=allocate(length+20)
	mpfr_sprintf(rstring, frmt, mpfr_rounding_mode, @num)
	if mpfr_sgn(@num)<0 then sp = ""
	function = sp + ltrim(*rstring)
	deallocate(rstring)
end function

Function mpfr.toLong ( ) As Long
	Function = mpfr_get_si(@num, mpfr_rounding_mode)
End Function

Function mpfr.toDouble ( ) As Double
	Function = mpfr_get_d(@num, mpfr_rounding_mode)
End Function

Function mpfr.toSingle ( ) As Single
	Function = mpfr_get_flt(@num, mpfr_rounding_mode)
End Function

Constructor mpfr ( )
	mpfr_init2(@num, mpfr_digits_precision * 3.33)
	mpfr_set_si(@num, 0, mpfr_rounding_mode)
End Constructor

Constructor mpfr ( Byval rhs As Long )
	mpfr_init2(@num, mpfr_digits_precision * 3.33)
	mpfr_set_si(@num, rhs, mpfr_rounding_mode)
End Constructor

Constructor mpfr ( Byref rhs As String )
	mpfr_init2(@num, mpfr_digits_precision * 3.33)
	mpfr_set_str(@num, rhs, 10, mpfr_rounding_mode)
End Constructor

Constructor mpfr ( Byref rhs As mpfr )
	mpfr_init2(@num, mpfr_digits_precision * 3.33)
	mpfr_set(@num, @rhs.num, mpfr_rounding_mode)
End Constructor

Constructor mpfr ( Byval rhs As LongInt )
	Dim As String s=str(rhs)
	mpfr_init2(@num, mpfr_digits_precision * 3.33)
	mpfr_set_str(@num, s, 10, mpfr_rounding_mode)
End Constructor

Constructor mpfr ( Byval rhs As Integer )
	Dim As String s=str(rhs)
	mpfr_init2(@num, mpfr_digits_precision * 3.33)
	mpfr_set_str(@num, s, 10, mpfr_rounding_mode)
End Constructor

Constructor mpfr ( Byval rhs As Double )
	mpfr_init2(@num, mpfr_digits_precision * 3.33)
	mpfr_set_d(@num, rhs, mpfr_rounding_mode)
End Constructor

Constructor mpfr ( Byval rhs As Single )
	mpfr_init2(@num, mpfr_digits_precision * 3.33)
	mpfr_set_flt(@num, rhs, mpfr_rounding_mode)
End Constructor

Destructor mpfr ( )
	mpfr_clear(@num)
End Destructor

Operator mpfr.let ( Byref rhs As mpfr )
	mpfr_set(@num, @rhs.num, mpfr_rounding_mode)
End Operator

Operator mpfr.let ( Byval rhs As Long )
	mpfr_set_si(@num, rhs, mpfr_rounding_mode)
End Operator

Operator mpfr.let ( Byref rhs As String )
	mpfr_set_str(@num, rhs, 10, mpfr_rounding_mode)
End Operator

Operator mpfr.Let ( Byval rhs As LongInt )
	Dim As String s=str(rhs)
	mpfr_set_str(@num, s, 10, mpfr_rounding_mode)
End Operator

Operator mpfr.Let ( Byval rhs As Integer )
	Dim As String s=str(rhs)
	mpfr_set_str(@num, s, 10, mpfr_rounding_mode)
End Operator

Operator mpfr.Let ( Byval rhs As Double )
	mpfr_set_d(@num, rhs, mpfr_rounding_mode)
End Operator

Operator mpfr.Let ( Byval rhs As Single )
	mpfr_set_flt(@num, rhs, mpfr_rounding_mode)
End Operator

Operator mpfr.cast ( ) As String
	dim as long length = mpfr_digits_precision
	dim as string sp =" ", s = "%"+str(length+4)+"."+str(length-1)+"R*g"
	length = mpfr_snprintf(0, 0, s, mpfr_rounding_mode, @num)
	dim as zstring ptr rstring=allocate(length+20)
	mpfr_sprintf(rstring, s, mpfr_rounding_mode, @num)
	if mpfr_sgn(@num)<0 then sp = ""
	Operator = sp + ltrim(*rstring)
	deallocate(rstring)
End Operator

Operator mpfr.cast ( ) As Long
	Operator = mpfr_get_si(@num, mpfr_rounding_mode)
End Operator

Operator mpfr.cast ( ) As Double
	Operator = mpfr_get_d(@num, mpfr_rounding_mode)
End Operator

Operator mpfr.cast ( ) As Single
	Operator = mpfr_get_flt(@num, mpfr_rounding_mode)
End Operator

'============================================================================
'' For Next for mpfr type
''
'' implicit step versions
'' 
'' In this example, we interpret implicit step
'' to mean 1
Operator mpfr.for ( )
End Operator
 
Operator mpfr.step ( )
        this += 1 'this = this+1 '
End Operator 
 
Operator mpfr.next ( Byref end_cond As mpfr ) As Integer
        Return this <= end_cond
End Operator
 
 
'' explicit step versions
'' 
Operator mpfr.for ( Byref step_var As mpfr )
End Operator
 
Operator mpfr.step ( Byref step_var As mpfr )
        this += step_var 'this = this + step_var '    
End Operator 
 
Operator mpfr.next ( Byref end_cond As mpfr, Byref step_var As mpfr ) As Integer
        If step_var < 0 Then
                Return this >= end_cond
        Else
                Return this <= end_cond
        End If
End Operator

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

Operator + ( Byref lhs As mpfr, Byref rhs As mpfr ) As mpfr
	Dim As mpfr result
	mpfr_add(@result.num, @lhs.num, @rhs.num, mpfr_rounding_mode)
	Operator = result
End Operator

Operator - ( Byref rhs As mpfr ) As mpfr
	Dim As mpfr result
	mpfr_neg(@result.num, @rhs.num, mpfr_rounding_mode)
	Operator = result
End Operator

Operator - ( Byref lhs As mpfr, Byref rhs As mpfr ) As mpfr
	Dim As mpfr result
	mpfr_sub(@result.num, @lhs.num, @rhs.num, mpfr_rounding_mode)
	Operator = result
End Operator

Operator * ( Byref lhs As mpfr, Byref rhs As mpfr ) As mpfr
	Dim As mpfr result
	mpfr_mul(@result.num, @lhs.num, @rhs.num, mpfr_rounding_mode)
	Operator = result
End Operator

Operator / ( Byref lhs As mpfr, Byref rhs As mpfr ) As mpfr
	Dim As mpfr result
	mpfr_div(@result.num, @lhs.num, @rhs.num, mpfr_rounding_mode)
	Operator = result
End Operator

Operator ^ ( Byref lhs As mpfr, Byval rhs As mpfr ) As mpfr
	Dim As mpfr result
	mpfr_pow(@result.num, @lhs.num, @rhs.num, mpfr_rounding_mode)
	Operator = result
End Operator
'=================
Operator mpfr.+= (byref rhs as mpfr)
	mpfr_add(@this.num, @this.num, @rhs.num, mpfr_rounding_mode)
End Operator

Operator mpfr.+= (Byval rhs as long)
	mpfr_add_si(@this.num, @this.num, rhs, mpfr_rounding_mode)
End Operator

Operator mpfr.+= (byval rhs as double)
	mpfr_add_d(@this.num, @this.num, rhs, mpfr_rounding_mode)
End Operator

Operator mpfr.+= (byval rhs as single)
	mpfr_add_d(@this.num, @this.num, rhs, mpfr_rounding_mode)
End Operator

Operator mpfr.+= (byval rhs as longint)
	Dim As mpfr tmp
	tmp=str(rhs)
	mpfr_add(@this.num, @this.num, @tmp.num, mpfr_rounding_mode)
End Operator

Operator mpfr.+= (byval rhs as integer)
	Dim As mpfr tmp
	tmp=str(rhs)
	mpfr_add(@this.num, @this.num, @tmp.num, mpfr_rounding_mode)
End Operator

Operator mpfr.+= (byref rhs as String)
	Dim As mpfr tmp
	tmp=rhs
	mpfr_add(@this.num, @this.num, @tmp.num, mpfr_rounding_mode)
End Operator

Operator mpfr.-= (byref rhs as mpfr)
	mpfr_sub(@this.num, @this.num, @rhs.num, mpfr_rounding_mode)
End Operator

Operator mpfr.-= (Byval rhs as long)
	mpfr_sub_si(@this.num, @this.num, rhs, mpfr_rounding_mode)
End Operator

Operator mpfr.-= (byval rhs as double)
	mpfr_sub_d(@this.num, @this.num, rhs, mpfr_rounding_mode)
End Operator

Operator mpfr.-= (byval rhs as single)
	mpfr_sub_d(@this.num, @this.num, rhs, mpfr_rounding_mode)
End Operator

Operator mpfr.-= (byval rhs as longint)
	Dim As mpfr tmp
	tmp=str(rhs)
	mpfr_sub(@this.num, @this.num, @tmp.num, mpfr_rounding_mode)
End Operator

Operator mpfr.-= (byval rhs as integer)
	Dim As mpfr tmp
	tmp=str(rhs)
	mpfr_sub(@this.num, @this.num, @tmp.num, mpfr_rounding_mode)
End Operator

Operator mpfr.-= (byref rhs as String)
	Dim As mpfr tmp
	tmp=rhs
	mpfr_sub(@this.num, @this.num, @tmp.num, mpfr_rounding_mode)
End Operator

Operator mpfr.*= (byref rhs as mpfr)
	mpfr_mul(@this.num, @this.num, @rhs.num, mpfr_rounding_mode)
End Operator

Operator mpfr.*= (Byval rhs as long)
	mpfr_mul_si(@this.num, @this.num, rhs, mpfr_rounding_mode)
End Operator

Operator mpfr.*= (byval rhs as double)
	Dim As mpfr tmp
	tmp=str(rhs)
	mpfr_mul_d(@this.num, @this.num, rhs, mpfr_rounding_mode)
End Operator

Operator mpfr.*= (byval rhs as single)
	mpfr_mul_d(@this.num, @this.num, rhs, mpfr_rounding_mode)
End Operator

Operator mpfr.*= (byval rhs as longint)
	Dim As mpfr tmp
	tmp=str(rhs)
	mpfr_mul(@this.num, @this.num, @tmp.num, mpfr_rounding_mode)
End Operator

Operator mpfr.*= (byval rhs as integer)
	Dim As mpfr tmp
	tmp=str(rhs)
	mpfr_mul(@this.num, @this.num, @tmp.num, mpfr_rounding_mode)
End Operator

Operator mpfr.*= (byref rhs as String)
	Dim As mpfr tmp
	tmp=rhs
	mpfr_mul(@this.num, @this.num, @tmp.num, mpfr_rounding_mode)
End Operator

Operator mpfr./= (byref rhs as mpfr)
	mpfr_div(@this.num, @this.num, @rhs.num, mpfr_rounding_mode)
End Operator

Operator mpfr./= (Byval rhs as long)
	mpfr_div_si(@this.num, @this.num, rhs, mpfr_rounding_mode)
End Operator

Operator mpfr./= (byval rhs as double)
	mpfr_div_d(@this.num, @this.num, rhs, mpfr_rounding_mode)
End Operator

Operator mpfr./= (byval rhs as single)
	mpfr_div_d(@this.num, @this.num, rhs, mpfr_rounding_mode)
End Operator

Operator mpfr./= (byval rhs as longint)
	Dim As mpfr tmp
	tmp=str(rhs)
	mpfr_div(@this.num, @this.num, @tmp.num, mpfr_rounding_mode)
End Operator

Operator mpfr./= (byval rhs as integer)
	Dim As mpfr tmp
	tmp=str(rhs)
	mpfr_div(@this.num, @this.num, @tmp.num, mpfr_rounding_mode)
End Operator

Operator mpfr./= (byref rhs as String)
	Dim As mpfr tmp
	tmp=rhs
	mpfr_div(@this.num, @this.num, @tmp.num, mpfr_rounding_mode)
End Operator

Operator mpfr.^= (byref rhs as mpfr)
	mpfr_pow(@this.num, @this.num, @rhs.num, mpfr_rounding_mode)
End Operator

Operator mpfr.^= (Byval rhs as long)
	mpfr_pow_si(@this.num, @this.num, rhs, mpfr_rounding_mode)
End Operator

Operator mpfr.^= (byval rhs as double)
	Dim As mpfr tmp
	tmp=str(rhs)
	mpfr_pow(@this.num, @this.num, @tmp.num, mpfr_rounding_mode)
End Operator

Operator mpfr.^= (byval rhs as single)
	Dim As mpfr tmp
	tmp=str(rhs)
	mpfr_pow(@this.num, @this.num, @tmp.num, mpfr_rounding_mode)
End Operator

Operator mpfr.^= (byval rhs as longint)
	Dim As mpfr tmp
	tmp=str(rhs)
	mpfr_pow(@this.num, @this.num, @tmp.num, mpfr_rounding_mode)
End Operator

Operator mpfr.^= (byval rhs as integer)
	Dim As mpfr tmp
	tmp=str(rhs)
	mpfr_pow(@this.num, @this.num, @tmp.num, mpfr_rounding_mode)
End Operator

Operator mpfr.^= (byref rhs as String)
	Dim As mpfr tmp
	tmp=rhs
	mpfr_pow(@this.num, @this.num, @tmp.num, mpfr_rounding_mode)
End Operator
    
Operator = ( Byref lhs As mpfr, Byref rhs As mpfr ) As Integer
	Operator = (mpfr_cmp(@lhs.num, @rhs.num) = 0)
End Operator

Operator < ( Byref lhs As mpfr, Byref rhs As mpfr ) As Integer
	Operator = (mpfr_cmp(@lhs.num, @rhs.num) < 0)
End Operator

Operator > ( Byref lhs As mpfr, Byref rhs As mpfr ) As Integer
	Operator = (mpfr_cmp(@lhs.num, @rhs.num) > 0)
End Operator

Operator <= ( Byref lhs As mpfr, Byref rhs As mpfr ) As Integer
	Operator = (mpfr_cmp(@lhs.num, @rhs.num) <= 0)
End Operator

Operator >= ( Byref lhs As mpfr, Byref rhs As mpfr ) As Integer
	Operator = (mpfr_cmp(@lhs.num, @rhs.num) >= 0)
End Operator

Operator <> ( Byref lhs As mpfr, Byref rhs As mpfr ) As Integer
	Operator = (mpfr_cmp(@lhs.num, @rhs.num) <> 0)
End Operator

operator Abs(byref rhs as mpfr) as mpfr
	Dim As mpfr result
	mpfr_abs(@result.num, @rhs.num, mpfr_rounding_mode)
	operator = result
end operator

operator Acos(byref rhs as mpfr) as mpfr
	Dim As mpfr result
	mpfr_acos(@result.num, @rhs.num, mpfr_rounding_mode)
	operator = result
end operator

function Acosh overload(byref rhs as mpfr) as mpfr
	Dim As mpfr result
	mpfr_acosh(@result.num, @rhs.num, mpfr_rounding_mode)
	Return result
end function

function Ai(byref rhs as mpfr) as mpfr
	Dim As mpfr result
	mpfr_ai(@result.num, @rhs.num, mpfr_rounding_mode)
	Return result
end function

operator Asin(byref rhs as mpfr) as mpfr
	Dim As mpfr result
	mpfr_asin(@result.num, @rhs.num, mpfr_rounding_mode)
	operator = result
end operator

function Asinh overload(byref rhs as mpfr) as mpfr
	Dim As mpfr result
	mpfr_asinh(@result.num, @rhs.num, mpfr_rounding_mode)
	Return result
end function

function Atan(byref rhs as mpfr) as mpfr
	Dim As mpfr result
	mpfr_atan(@result.num, @rhs.num, mpfr_rounding_mode)
	Return result
end function

operator Atn(byref rhs as mpfr) as mpfr
	Dim As mpfr result
	mpfr_atan(@result.num, @rhs.num, mpfr_rounding_mode)
	operator = result
end operator

Function Atan2_ ( ByVal y As Double, ByVal x As Double ) As Double
	Function = Atan2 ( y, x )
End Function

#undef Atan2

Declare Function Atan2 Overload ( ByVal y As Double, ByVal x As Double ) As Double

Function Atan2 ( ByVal y As Double, ByVal x As Double ) As Double
	Function = Atan2_ ( y, x )
End Function

function Atan2(byref lhs as mpfr, byref rhs as mpfr) as mpfr
	Dim As mpfr result
	mpfr_atan2 (@result.num, @lhs.num, @rhs.num, mpfr_rounding_mode)
	Return result
end function

function Atanh overload(byref rhs as mpfr) as mpfr
	Dim As mpfr result
	mpfr_atanh(@result.num, @rhs.num, mpfr_rounding_mode)
	Return result
end function

function cbrt(byref x as mpfr) as mpfr
	Dim As mpfr result
	mpfr_cbrt(@result.num, @x.num, mpfr_rounding_mode)
	return result
end function

function ceil(byref x as mpfr) as mpfr
	Dim As mpfr result
	mpfr_ceil(@result.num, @x.num)
	function = result
end function

function floor(byref x as mpfr) as mpfr
	Dim As mpfr result
	mpfr_floor(@result.num, @x.num)
	return result
end function

operator int(byref x as mpfr) as mpfr
	Dim As mpfr result
	mpfr_floor(@result.num, @x.num)
	operator = result
end operator

operator frac(byref x as mpfr) as mpfr
	Dim As mpfr result
	mpfr_frac(@result.num, @x.num, mpfr_rounding_mode)
	operator = result
end operator

function catalan_const() as mpfr
	Dim As mpfr result
	mpfr_const_catalan(@result.num, mpfr_rounding_mode)
	return result
end function

function euler_const() as mpfr
	Dim As mpfr result
	mpfr_const_euler(@result.num, mpfr_rounding_mode)
	return result
end function

function log2_const() as mpfr
	Dim As mpfr result
	mpfr_const_log2(@result.num, mpfr_rounding_mode)
	return result
end function

function pi_const() as mpfr
	Dim As mpfr result
	mpfr_const_pi(@result.num, mpfr_rounding_mode)
	return result
end function

operator Cos(byref rhs as mpfr) as mpfr
	Dim As mpfr result
	mpfr_cos(@result.num, @rhs.num, mpfr_rounding_mode)
	operator = result
end operator

function Cosh overload(byref rhs as mpfr) as mpfr
	Dim As mpfr result
	mpfr_cosh(@result.num, @rhs.num, mpfr_rounding_mode)
	Return result
end function

function Cot(byref rhs as mpfr) as mpfr
	Dim As mpfr result
	mpfr_cot(@result.num, @rhs.num, mpfr_rounding_mode)
	Return result
end function

function Coth(byref rhs as mpfr) as mpfr
	Dim As mpfr result
	mpfr_coth(@result.num, @rhs.num, mpfr_rounding_mode)
	Return result
end function

function Csc(byref rhs as mpfr) as mpfr
	Dim As mpfr result
	mpfr_csc(@result.num, @rhs.num, mpfr_rounding_mode)
	Return result
end function

function Csch(byref rhs as mpfr) as mpfr
	Dim As mpfr result
	mpfr_csch(@result.num, @rhs.num, mpfr_rounding_mode)
	Return result
end function

function Digamma(byref rhs as mpfr) as mpfr
	Dim As mpfr result
	mpfr_digamma(@result.num, @rhs.num, mpfr_rounding_mode)
	Return result
end function

function Eint(byref rhs as mpfr) as mpfr
	Dim As mpfr result
	mpfr_eint(@result.num, @rhs.num, mpfr_rounding_mode)
	Return result
end function

function Erf(byref rhs as mpfr) as mpfr
	Dim As mpfr result
	mpfr_erf(@result.num, @rhs.num, mpfr_rounding_mode)
	Return result
end function

function Erfc(byref rhs as mpfr) as mpfr
	Dim As mpfr result
	mpfr_erfc(@result.num, @rhs.num, mpfr_rounding_mode)
	Return result
end function

operator Exp(byref x as mpfr) as mpfr
	dim as mpfr result
	mpfr_exp(@result.num, @x.num, mpfr_rounding_mode)
	operator = result
end operator

function Exp10 overload(byref x as mpfr) as mpfr
	dim as mpfr result
	mpfr_exp10(@result.num, @x.num, mpfr_rounding_mode)
	return result
end function

function Exp2(byref x as mpfr) as mpfr
	dim as mpfr result
	mpfr_exp2(@result.num, @x.num, mpfr_rounding_mode)
	return result
end function

function Expm1(byref x as mpfr) as mpfr
	dim as mpfr result
	mpfr_expm1(@result.num, @x.num, mpfr_rounding_mode)
	return result
end function

function facui(byval x as culong) as mpfr
	Dim As mpfr result
	mpfr_fac_ui(@result.num, x, mpfr_rounding_mode)
	return result
end function

function fma(byref a as mpfr, byref b as mpfr, byref c as mpfr) as mpfr
	Dim As mpfr result
	mpfr_fma(@result.num, @a.num, @b.num, @c.num, mpfr_rounding_mode)
	return result
end function

function fms(byref a as mpfr, byref b as mpfr, byref c as mpfr) as mpfr
	Dim As mpfr result
	mpfr_fms(@result.num, @a.num, @b.num, @c.num, mpfr_rounding_mode)
	return result
end function

function gamma(byref x as mpfr) as mpfr
	Dim As mpfr result
	mpfr_gamma(@result.num, @x.num, mpfr_rounding_mode)
	return result
end function

function hypot(byref x as mpfr, byref y as mpfr) as mpfr
	Dim As mpfr result
	mpfr_hypot(@result.num, @x.num, @y.num, mpfr_rounding_mode)
	return result
end function

function j0(byref x as mpfr) as mpfr
	Dim As mpfr result
	mpfr_j0(@result.num, @x.num, mpfr_rounding_mode)
	return result
end function

function j1(byref x as mpfr) as mpfr
	Dim As mpfr result
	mpfr_j1(@result.num, @x.num, mpfr_rounding_mode)
	return result
end function

function jn(byref x as mpfr, byval n as clong) as mpfr
	Dim As mpfr result
	mpfr_jn(@result.num, n, @x.num, mpfr_rounding_mode)
	return result
end function

function y0(byref x as mpfr) as mpfr
	Dim As mpfr result
	mpfr_y0(@result.num, @x.num, mpfr_rounding_mode)
	return result
end function

function y1(byref x as mpfr) as mpfr
	Dim As mpfr result
	mpfr_y1(@result.num, @x.num, mpfr_rounding_mode)
	return result
end function

function yn(byref x as mpfr, byval n as clong) as mpfr
	Dim As mpfr result
	mpfr_yn(@result.num, n, @x.num, mpfr_rounding_mode)
	return result
end function

function Zeta(byref x as mpfr) as mpfr
	Dim As mpfr result
	mpfr_zeta(@result.num, @x.num, mpfr_rounding_mode)
	return result
end function

function ZetaUi(byval n as clong) as mpfr
	Dim As mpfr result
	mpfr_zeta_ui(@result.num, n, mpfr_rounding_mode)
	return result
end function

function lgamma(byval signp as long ptr, byref x as mpfr) as mpfr
	Dim As mpfr result
	mpfr_lgamma(@result.num, signp, @x.num, mpfr_rounding_mode)
	return result
end function

function Li2(byref x as mpfr) as mpfr
	Dim As mpfr result
	mpfr_li2(@result.num, @x.num, mpfr_rounding_mode)
	return result
end function

function lngamma(byref x as mpfr) as mpfr
	Dim As mpfr result
	mpfr_lngamma(@result.num, @x.num, mpfr_rounding_mode)
	return result
end function

operator Log(byref x as mpfr) as mpfr
	Dim As mpfr result
	mpfr_log(@result.num, @x.num, mpfr_rounding_mode)
	operator = result
end operator

function Log10 overload(byref x as mpfr) as mpfr
	Dim As mpfr result
	mpfr_log10(@result.num, @x.num, mpfr_rounding_mode)
	return result
end function

function Log1p(byref x as mpfr) as mpfr
	Dim As mpfr result
	mpfr_log1p(@result.num, @x.num, mpfr_rounding_mode)
	return result
end function

function Log2(byref x as mpfr) as mpfr
	Dim As mpfr result
	mpfr_log2(@result.num, @x.num, mpfr_rounding_mode)
	return result
end function

function Sec(byref x as mpfr) as mpfr
	Dim As mpfr result
	mpfr_sec(@result.num, @x.num, mpfr_rounding_mode)
	return result
end function

function Sech(byref x as mpfr) as mpfr
	Dim As mpfr result
	mpfr_sech(@result.num, @x.num, mpfr_rounding_mode)
	return result
end function

operator Sin(byref x as mpfr) as mpfr
	Dim As mpfr result
	mpfr_sin(@result.num, @x.num, mpfr_rounding_mode)
	operator = result
end operator

function Sinh overload(byref x as mpfr) as mpfr
	Dim As mpfr result
	mpfr_sinh(@result.num, @x.num, mpfr_rounding_mode)
	return result
end function

Function Sqr_ ( ByVal number As Double ) As Double
	Function = Sqr( number )
End Function

#undef Sqr

Declare Function Sqr Overload ( ByVal number As Double ) As Double

Function Sqr ( ByVal number As Double ) As Double
	Function = Sqr_( number )
End Function

function Sqr(byref x as mpfr) as mpfr
	Dim As mpfr result
	mpfr_sqrt(@result.num, @x.num, mpfr_rounding_mode)
	return result
end function

function Square(byref x as mpfr) as mpfr
	Dim As mpfr result
	mpfr_sqr(@result.num, @x.num, mpfr_rounding_mode)
	return result
end function

operator Tan(byref x as mpfr) as mpfr
	Dim As mpfr result
	mpfr_tan(@result.num, @x.num, mpfr_rounding_mode)
	operator = result
end operator

function Tanh overload(byref x as mpfr) as mpfr
	Dim As mpfr result
	mpfr_tanh(@result.num, @x.num, mpfr_rounding_mode)
	return result
end function

function Trunc(byref x as mpfr) as mpfr
	Dim As mpfr result
	mpfr_trunc(@result.num, @x.num)
	return result
end function
/'
function str_(byval x as double) as string
	return str(x)
end function

#undef str

declare function str overload (byval x as double) as string

function str(byval x as double) as string
	return str_(x)
end function

function str (byref x as mpfr ) As String
	dim as long length = mpfr_digits_precision
	dim as string sp =" ", s = "%"+str(length+4)+"."+str(length-1)+"R*g"
	length = mpfr_snprintf(0, 0, s, mpfr_rounding_mode, @x.num)
	dim as zstring ptr rstring=allocate(length+20)
	mpfr_sprintf(rstring, s, mpfr_rounding_mode, @x.num)
	if mpfr_sgn(@x.num)<0 then sp = ""
	return sp + ltrim(*rstring)
	deallocate(rstring)
End function
'/

[srvaldez]

the result of the first post using double precision

Code: Select all

Coefficiennts for the function 1/tan(x)
1/tan(x) =
x^-1 * 0.9999998020688422 -
x    * 0.333323900070707 -
x^3  * 0.02234139534929878 -
x^5  * 0.001480866524972864 -
x^7  * 0.001904747469230853 +
x^9  * 0.002349009984637363 -
x^11 * 0.001668966115876747 +
x^13 * 0.0004636904788444083 

error =  1.979217905467521e-07
relative error for x=  1 -> -1.723724407654004e-08
relative error for x= -1 -> -1.723724407654004e-08
relative error for x=0.5 ->  3.843739826425718e-08
elapsed time =  0.07604098320007324 seconds

result using mpfr

Code: Select all

Coefficiennts for the function 1/tan(x)
1/tan(x) =
x^-1 * 1.000000000067919750594788500783910857240322905811022265102052 -
x    * 0.333333333235636776694714738097455546905902779599336266954574 -
x^3  * 0.0222222227162403178330842260590083360568861492842891724008522 -
x^5  * 0.002116406633270569954632617968645639855173204749875414475236964 -
x^7  * 0.0002116020744501638786828692328701762983705541051988735986960567 -
x^9  * 2.148910404795801007184487179563973277139470187989022801036064e-05 -
x^11 * 2.00927224427890133199538817851653537394562391911124783341476e-06 -
x^13 * 3.211019048188598209448689971257997031495795870003330668035861e-07 

error =  7.356747960677978820219615799439810124280595131186035166692092e-11
relative error for x=  1 ->  6.550202322215880269477162625665154966466395887223503124289759e-12
relative error for x= -1 ->  6.550202322215880269477162625665154966466395887223503124289759e-12
relative error for x=0.5 -> -6.902504101480611696226961956416046911356492084357866302071388e-11
elapsed time =  70.3135039806366 seconds

[srvaldez]

I found this at https://stackoverflow.com/a/26825029

Code: Select all

'' by njuffa
'' https://stackoverflow.com/a/26825029

function atanf_poly(byval a as single) as single
	dim as single r, s

	s = a * a
	'r = (1+(((((((0.278569828e-2*s-0.158660226e-1)*s+0.424722321e-1)*s-0.749753043e-1)*s+.106448799)*s-.142070308)*s+.199934542)*s-.333331466)*s)*a
	r = a*s*(((((((0.278569828e-2*s-0.158660226e-1)*s+0.424722321e-1)*s-0.749753043e-1)*s+.106448799)*s-.142070308)*s+.199934542)*s-.333331466)+a
	/'
	r = 2.78569828e-3
	r = (r * s) - 1.58660226e-2
	r = (r * s) + 4.24722321e-2
	r = (r * s) - 7.49753043e-2
	r = (r * s) + 1.06448799e-1
	r = (r * s) - 1.42070308e-1
	r = (r * s) + 1.99934542e-1
	r = (r * s) - 3.33331466e-1
	r = r * s
	r = (r * a) + a'/
	return r
end function

? atanf_poly(-1), atanf_poly(-1)-atn(-1)
? atanf_poly(1), atanf_poly(1)-atn(1)
? atanf_poly(0.5), atanf_poly(0.5)-atn(0.5)
? atanf_poly(0.1), atanf_poly(0.1)-atn(0.1)

[dodicat]

Here is a double [-2,2] approximation.

Code: Select all

  
screen 19,32
 #macro sketch(_function,minx,maxx,miny,maxy,col)
For x As Double=minx To maxx Step (maxx-minx)/500
 Dim As Double x1=(800)*(x-minx)/(maxx-minx)
 Dim As Double y1=(600)*(_function-maxy)/(miny-maxy)
 If x=minx Then Pset(x1,y1),col Else Line -(x1,y1),col
Next x
#endmacro
 

function f(x as double) as double
' approximating atn(x)  [-2 To 2]
Return _ 
((((((((((+0.001918512641699786)*x+5.801349831568384e-017)*x-0.02233443494461894)*x-5.33883172588491e-016)*x _ 
+0.1053572142797391)*x _ 
+1.601515686559807e-015)*x _ 
-0.297623118874788)*x _ 
-1.692981692336115e-015)*x _ 
+0.9973075167764144)*x _ 
+5.282232984349378e-016) _ 

end function
 
sketch(atn(x),-2,2,-1.10714871779409,1.099018815614147,Rgb(200,0,0))
sketch(f(x),-2,2,-1.10714871779409,1.099018815614147,Rgb(0,0,200))

? f(-1), f(-1)-atn(-1)
? f(1), f(1)-atn(1)
? f(0.5), f(0.5)-atn(0.5)
? f(0.1), f(0.1)-atn(0.1)
print
print "RED atn, BLUE approximation"
sleep
 

[srvaldez]

hi dodicat
you can get an additional digit of accuracy by using this

Code: Select all

(81.54559687151539*(4.531540274801958e-14+(.4118594774380018+(9.654405239783765e-11+(.1235027084810615+1.063666015359081e-11*x)*x)*x)*x))/((19.20815561560766+(2.979586874916466e-9+x)*x)*(1.750707902811483+(1.447504036147797e-10+x)*x))

btw, nice plot

[srvaldez]

hello dodicat
would you have a look at this pdf and let me know what you think? https://pdfs.semanticscholar.org/43b1/b ... 82c447.pdf
it has source code in TurboBasic, so it should be relatively easy to port to FB if you think it's worth the trouble

[dodicat]

Thanks srvaldez.
with their first function
2*exp(-x)*sin(x) +exp(-2*x) -x^2 +4
although they have it
2*exp(-ɯ)*sin(ɯ) +exp(-2*ɯ) -x^2 +4
(mixing omega and x for some reason)

I tried with chebyshev, got slightly different coefficients.
Their basic code is very long, and very hard to copy and paste from a .pdf file.

[marcov]

Code: Select all

function rolldice () as double
return 4   ' I rolled a dice, and this was the outcome. Really random,  no cheating!
end function

I tried to make an AI using Tensorflow, but decided it was useless since it couldn't even predict the roll of simple dice.

[dodicat]

A great function marcov, for chucking a die.

Code: Select all

#undef return 
#define Intrange(f,l) int(Rnd*(((l)+1)-(f)))+(f)
#define return  function = intRange(5,10)- 


function rolldice () as double
return 4   ' I rolled a dice, and this was the outcome. Really random,  no cheating!
end function


dim as integer a(1 to 6)
for n as long=1 to 1000000
    a(rolldice)+=1
next

for n as long=1 to 6
    print a(n)
next
sleep

  

Thank you very much, I shall always use it.