FBMath update
FBMath update
The FreeBASIC math library (FBMath) has been updated. The modifications are as follows:
* Added: Complex numbers and functions
* Added: Expression evaluation
* Added: HSV / RGB conversion
* The LargeInt library (contributed by S. J. Schaper) is now fully
integrated into FBMath
* The library has its own version of the "Mersenne Twister" random
number generator and does not call the built-in functions RND or
RANDOMIZE anymore (this was done for solving a problem with FB 0.24)
* The compiler directive USE_PREFIX modifies the name of some library
functions (already declared in windows.bi or math.bi) by prefixing
them with 'fb'. The following functions are concerned:
Min, Max, Exp10, Exp2, Log10, Log2, Logf, Round, Ceil, Floor,
Sinh, Cosh, Tanh, ASinh, ACosh, ATanh, Erf, Erfc
The directive renames these functions as fbMin, fbMax, etc.
It is therefore possible to include windows.bi, math.bi or crt.bi
in addition to fbmath.bi in the same program.
* The compilation scripts compil_prefix.bat and compil_prefix.sh use
this compiler directive.
* The option "-w pedantic" is always used in the compilation scripts.
All BYVAL or BYREF attributes are explicitly set.
Download link:
http://sourceforge.net/projects/fbmath/ ... mat060.zip
* Added: Complex numbers and functions
* Added: Expression evaluation
* Added: HSV / RGB conversion
* The LargeInt library (contributed by S. J. Schaper) is now fully
integrated into FBMath
* The library has its own version of the "Mersenne Twister" random
number generator and does not call the built-in functions RND or
RANDOMIZE anymore (this was done for solving a problem with FB 0.24)
* The compiler directive USE_PREFIX modifies the name of some library
functions (already declared in windows.bi or math.bi) by prefixing
them with 'fb'. The following functions are concerned:
Min, Max, Exp10, Exp2, Log10, Log2, Logf, Round, Ceil, Floor,
Sinh, Cosh, Tanh, ASinh, ACosh, ATanh, Erf, Erfc
The directive renames these functions as fbMin, fbMax, etc.
It is therefore possible to include windows.bi, math.bi or crt.bi
in addition to fbmath.bi in the same program.
* The compilation scripts compil_prefix.bat and compil_prefix.sh use
this compiler directive.
* The option "-w pedantic" is always used in the compilation scripts.
All BYVAL or BYREF attributes are explicitly set.
Download link:
http://sourceforge.net/projects/fbmath/ ... mat060.zip
Re: FBMath update
jdebord
While trying out the new fbc version I had to take a closer look into the working of compil.bat (fbmath 060), I came across this error messages.
Looks like that
DIM R() AS DOUBLE
DIM C() AS Complex
is missing.
While trying out the new fbc version I had to take a closer look into the working of compil.bat (fbmath 060), I came across this error messages.
Code: Select all
for %%a in (polynom\*.bas) do %Cmdlin%
polynom\polyutil.bas(71) error 41: Variable not declared, R in 'R(J) = Z(I).X'
polynom\polyutil.bas(74) error 41: Variable not declared, C in 'C(K).X = Z(I).X'
polynom\polyutil.bas(74) error 3: Expected End-of-Line, found 'X' in 'C(K).X = Z(I).X'
polynom\polyutil.bas(75) error 70: Array not dimensioned, before '(' in 'C(K).Y = Z(I).Y'
polynom\polyutil.bas(81) error 70: Array not dimensioned, before '(' in 'K = I: A = R(I)'
polynom\polyutil.bas(83) error 70: Array not dimensioned, before '(' in 'IF R(J) < A THEN K = J: A = R(J)'
polynom\polyutil.bas(83) error 70: Array not dimensioned, before '(' in 'IF R(J) < A THEN K = J: A = R(J)'
polynom\polyutil.bas(85) error 70: Array not dimensioned, before '(' in 'SWAP R(I), R(K)'
polynom\polyutil.bas(90) error 70: Array not dimensioned, before '(' in 'Z(I).X = R(I)'
polynom\polyutil.bas(97) error 70: Array not dimensioned, before '(' in 'Z(J).X = C(I).X'
polynom\polyutil.bas(97) error 125: Too many errors, exiting
D:\FreeBASIC-0.90.0rc1-win32\bin\win32\ar.exe r libfbmath.a polynom\*.o
DIM R() AS DOUBLE
DIM C() AS Complex
is missing.
Re: FBMath update
These arrays are dimensioned by REDIM inside the subroutine SortRoots. This does not seem to give problems with FB 0.24
For FB 0.90 you can replace the two 'offending' lines with:
Code: Select all
' ******************************************************************
' Utility functions to handle roots of polynomials
' ******************************************************************
' ------------------------------------------------------------------
' Type definition
' ------------------------------------------------------------------
TYPE Complex
X AS DOUBLE
Y AS DOUBLE
END TYPE
' ******************************************************************
FUNCTION SetRealRoots(BYVAL Deg AS INTEGER,_
Z() AS Complex,_
BYVAL Tol AS DOUBLE) AS INTEGER
' ------------------------------------------------------------------
' Set the imaginary part of a root to zero if it is less than a
' fraction Tol of its real part. This root is therefore considered
' real. The function returns the total number of real roots.
' ------------------------------------------------------------------
DIM AS INTEGER I, N
FOR I = 1 TO Deg
IF Z(I).Y <> 0 AND ABS(Z(I).Y) < Tol * ABS(Z(I).X) THEN
Z(I).Y = 0
END IF
NEXT I
' Count real roots
N = 0
FOR I = 1 TO Deg
IF Z(I).Y = 0 THEN N = N + 1
NEXT I
RETURN N
END FUNCTION
SUB SortRoots(BYVAL Deg AS INTEGER, Z() AS Complex)
' ------------------------------------------------------------------
' Sorts roots so that:
'
' (1) The N real roots are stored in elements 1..N of vector Z,
' in increasing order.
'
' (2) The complex roots are stored in elements (N + 1)..Deg of
' vector Z and are unordered.
' ------------------------------------------------------------------
DIM AS DOUBLE A
DIM AS INTEGER I, J, K, Nr, Nc
' Count real and complex roots
Nr = 0: Nc = 0
FOR I = 1 TO Deg
IF Z(I).Y = 0 THEN Nr = Nr + 1 ELSE Nc = Nc + 1
NEXT I
IF Nr > 0 THEN REDIM R(1 TO Nr) AS DOUBLE
IF Nc > 0 THEN REDIM C(1 TO Nc) AS Complex
' Store real roots in R and complex roots in (X,Y)
J = 0: K = 0
FOR I = 1 TO Deg
IF Z(I).Y = 0 THEN
J = J + 1
R(J) = Z(I).X
ELSE
K = K + 1
C(K).X = Z(I).X
C(K).Y = Z(I).Y
END IF
NEXT I
' Sort real roots
FOR I = 1 TO Nr - 1
K = I: A = R(I)
FOR J = I + 1 TO Nr
IF R(J) < A THEN K = J: A = R(J)
NEXT J
SWAP R(I), R(K)
NEXT I
' Transfer real roots into elements 1..Nr
FOR I = 1 TO Nr
Z(I).X = R(I)
Z(I).Y = 0
NEXT I
' Transfer complex roots into elements (Nr+1)..Deg
FOR I = 1 TO Nc
J = I + Nr
Z(J).X = C(I).X
Z(J).Y = C(I).Y
NEXT I
END SUB
Code: Select all
DIM R(1 TO Nr) AS DOUBLE
DIM C(1 TO Nc) AS Complex
Re: FBMath update
Your right, I messed up, did not check if older FB versions worked correctly.jdebord wrote:These arrays are dimensioned by REDIM inside the subroutine SortRoots. This does not seem to give problems with FB 0.24
The problem is with the new FB version.
Re: FBMath update
that looks like it was due to this bug: #654
-
- Site Admin
- Posts: 6323
- Joined: Jul 05, 2005 17:32
- Location: Manchester, Lancs
Re: FBMath update
In my understanding of Redim, this is what is going on:jdebord wrote:For FB 0.90 you can replace the two 'offending' lines with:
Code: Select all
DIM R(1 TO Nr) AS DOUBLE DIM C(1 TO Nc) AS Complex
In 0.24, if R() and C() exist, redim would resize them, while dim would re-define them (unless their previous definition was in the same scope). The defined arrays would AFAIK be uninitialised if the If block was skipped.
If R() and C() don't already exist, then dim/redim as type would define them. Redim (or dim) without a type would give a compile error. Again, the arrays would AFAIK be uninitialised if the If block was skipped.
In 0.90, if R() and C() exist, again redim would resize, while dim would re-define them for the (short) scope of the If block.
If they don't exist, again dim/redim as type would define them, requiring the type name to avoid a compile error. Their scope would again be limited to the single line, and they would be unusable after that.
I can't see any case where changing the redim to dim would fix the code in 0.90. But I do know that the Dim will be effectively useless.
All I can suggest is to ensure the arrays are defined somewhere before the If block, and remove the 'as type' on the redims so they can't be misconstrued as definitions.
Re: FBMath update
jdebord
Here some things that I came across and some random thoughts about FBMath. How you want to do things is your business, I's not my intention to criticize your project, merely to improve things (I know that I can be wrong). Most things will save some memory and/or computing time, most of time not much but every helps. I haved looked at all of the code.
In “complex.bas”
FB swap works also on type's so, swap x,y gives the same result.
This part does nothing because in the begin there is a test to see if Z.Y = 0.
ELSEIF AbsY = 0 THEN
RETURN 0
Better remove those two lines.
There is a situation that is not tested, if AbsY = AbsX.
The answer would be AbsY (or AbsX.) * SQR(2). The routine can handle this situation and a other question would be how often this situation occurs.
Just curios, why not simply calculate SQR(Z.X * Z.X + Z.Y * Z.Y), the FPU can internally handle numbers up to approx 1.19E4932. Two doubles would only go as high as approx 3.58E618, well below what the FPU can handle.(After the two first IF Then statements)
In SUB SinhCosh(BYVAL X AS DOUBLE, BYREF SinhX AS DOUBLE, BYREF CoshX AS DOUBLE)
' Compute Sinh(X) and Cosh(X) with only one exponential
You have the following calculation
SinhX = 0.5 * (ExpX - ExpMinusX)
CoshX = 0.5 * (ExpX + ExpMinusX)
The calculation for CoshX can be replaced by CoshX = SinhX + ExpMinusX this will save a FPU multiplication ( approx 35 clock ticks). In “ROOTPOL4.BAS” you have the same situation.
In “POLYUTIL.BAS”
Nr = 0: Nc = 0 You just declared them so this is not needed they are already “0”.
J = 0: K = 0 At this point J and K aren't used so the are still “0”.
In “HSVRGB.BAS” in the sub routine HSVtoRGB
IF S = 0 THEN ' achromatic (grey)
R = V * 255
G = V * 255
B = V * 255
EXIT SUB
END IF
Why not simply do
R = V * 255
G = R
B = R
This will save some memory and some clock ticks (R, G and B are integers and V is a double so the FPU calculates V * 255 and then converts the result into a integer).
Convert
DIM AS INTEGER I
DIM AS DOUBLE Z, F, P, Q, T
Z = H / 60 ' sector 0 to 5
I = INT(Z)
F = frac(Z)
P = V * (1 - S)
Q = V * (1 - S * F)
T = V * (1 - S * (1 – F)
into
DIM AS DOUBLE Z = H / 60 ' sector 0 to 5
DIM AS INTEGER I = INT(Z)
DIM AS DOUBLE F = frac(Z)
DIM AS DOUBLE P = V * (1 - S)
DIM AS DOUBLE Q = V * (1 - S * F)
DIM AS DOUBLE T = V * (1 - S * (1 – F))
I have seen code where you did this.
And for the FUNCTION Hex2(BYVAL i AS UBYTE) AS STRING read the wiki about HEX.
Regards.
Here some things that I came across and some random thoughts about FBMath. How you want to do things is your business, I's not my intention to criticize your project, merely to improve things (I know that I can be wrong). Most things will save some memory and/or computing time, most of time not much but every helps. I haved looked at all of the code.
In “complex.bas”
Code: Select all
SUB CSwap(BYREF X AS Complex, BYREF Y AS Complex)
' Exchanges two Complexes
DIM AS Complex Temp
ErrCode = FOk
Temp = X
X = Y
Y = Temp
END SUB
Code: Select all
FUNCTION CAbs(BYREF Z AS Complex) AS DOUBLE
' Modulus of Z
ErrCode = FOk
IF Z.X = 0 THEN RETURN ABS(Z.Y)
IF Z.Y = 0 THEN RETURN ABS(Z.X)
DIM AS DOUBLE AbsX, AbsY, R
AbsX = ABS(Z.X)
AbsY = ABS(Z.Y)
IF AbsX > AbsY THEN
R = AbsY / AbsX
RETURN AbsX * SQR(1 + R * R)
ELSEIF AbsY = 0 THEN
RETURN 0
ELSE
R = AbsX / AbsY
RETURN AbsY * SQR(1 + R * R)
END IF
END FUNCTION
ELSEIF AbsY = 0 THEN
RETURN 0
Better remove those two lines.
There is a situation that is not tested, if AbsY = AbsX.
The answer would be AbsY (or AbsX.) * SQR(2). The routine can handle this situation and a other question would be how often this situation occurs.
Just curios, why not simply calculate SQR(Z.X * Z.X + Z.Y * Z.Y), the FPU can internally handle numbers up to approx 1.19E4932. Two doubles would only go as high as approx 3.58E618, well below what the FPU can handle.(After the two first IF Then statements)
In SUB SinhCosh(BYVAL X AS DOUBLE, BYREF SinhX AS DOUBLE, BYREF CoshX AS DOUBLE)
' Compute Sinh(X) and Cosh(X) with only one exponential
You have the following calculation
SinhX = 0.5 * (ExpX - ExpMinusX)
CoshX = 0.5 * (ExpX + ExpMinusX)
The calculation for CoshX can be replaced by CoshX = SinhX + ExpMinusX this will save a FPU multiplication ( approx 35 clock ticks). In “ROOTPOL4.BAS” you have the same situation.
In “POLYUTIL.BAS”
Code: Select all
SUB SortRoots(BYVAL Deg AS INTEGER, Z() AS Complex)
<left some code out>
DIM AS INTEGER I, J, K, Nr, Nc
' Count real and complex roots
Nr = 0: Nc = 0
FOR I = 1 TO Deg
IF Z(I).Y = 0 THEN Nr = Nr + 1 ELSE Nc = Nc + 1
NEXT I
IF Nr > 0 THEN REDIM R(1 TO Nr) AS DOUBLE
IF Nc > 0 THEN REDIM C(1 TO Nc) AS Complex
' Store real roots in R and complex roots in (X,Y)
J = 0: K = 0
< removed the rest >
J = 0: K = 0 At this point J and K aren't used so the are still “0”.
In “HSVRGB.BAS” in the sub routine HSVtoRGB
IF S = 0 THEN ' achromatic (grey)
R = V * 255
G = V * 255
B = V * 255
EXIT SUB
END IF
Why not simply do
R = V * 255
G = R
B = R
This will save some memory and some clock ticks (R, G and B are integers and V is a double so the FPU calculates V * 255 and then converts the result into a integer).
Convert
DIM AS INTEGER I
DIM AS DOUBLE Z, F, P, Q, T
Z = H / 60 ' sector 0 to 5
I = INT(Z)
F = frac(Z)
P = V * (1 - S)
Q = V * (1 - S * F)
T = V * (1 - S * (1 – F)
into
DIM AS DOUBLE Z = H / 60 ' sector 0 to 5
DIM AS INTEGER I = INT(Z)
DIM AS DOUBLE F = frac(Z)
DIM AS DOUBLE P = V * (1 - S)
DIM AS DOUBLE Q = V * (1 - S * F)
DIM AS DOUBLE T = V * (1 - S * (1 – F))
I have seen code where you did this.
And for the FUNCTION Hex2(BYVAL i AS UBYTE) AS STRING read the wiki about HEX.
Regards.
Re: FBMath update
Hello frisian. Thank you for the constructive criticism :)
I have modified COMPLEX.BAS and HSVRGB.BAS according to your suggestions. I have also modified POLYUTIL.BAS so that it is now compatible with FB 0.90. The modified files are in the SVN.
Some remarks:
1) I like explicit initializations such as I = 0 because they make the algorithms more obvious, at least for me! Also, I have several versions of this library, for several compilers, and some of them do not initialize variables to zero.
2) Compute SQR(Z.X * Z.X + Z.Y * Z.Y) : an overflow occurs if the argument of SQR exceeds the maximum value of a 64-bit real.
3) In SinhCosh, your modification will probably change the error pattern so it deserves further study. Also, how does the FPU handle a multiplication by 0.5 = 2^(-1)? Right shift? This should be very fast.
I have modified COMPLEX.BAS and HSVRGB.BAS according to your suggestions. I have also modified POLYUTIL.BAS so that it is now compatible with FB 0.90. The modified files are in the SVN.
Some remarks:
1) I like explicit initializations such as I = 0 because they make the algorithms more obvious, at least for me! Also, I have several versions of this library, for several compilers, and some of them do not initialize variables to zero.
2) Compute SQR(Z.X * Z.X + Z.Y * Z.Y) : an overflow occurs if the argument of SQR exceeds the maximum value of a 64-bit real.
3) In SinhCosh, your modification will probably change the error pattern so it deserves further study. Also, how does the FPU handle a multiplication by 0.5 = 2^(-1)? Right shift? This should be very fast.
Re: FBMath update
If Z.X and Z.Y are both MaxNum (1.79e308) then the answer would be MaxNum*SQR(2). This would be greater than the maximum value. But every routine to calculate the answer will exceed the maximum value of a 64-bit real.jdebord wrote:2) Compute SQR(Z.X * Z.X + Z.Y * Z.Y) : an overflow occurs if the argument of SQR exceeds the maximum value of a 64-bit real.
jdebord wrote:3) In SinhCosh, your modification will probably change the error pattern so it deserves further study. Also, how does the FPU handle a multiplication by 0.5 = 2^(-1)? Right shift? This should be very fast.
Indeed the modification could give a slighty different result but that would be a round off error by the FPU, but that is something every floating point calculation faces, when a 80-bit real has to be stored in a 64-bit real. That is why it is better to avoid storing a result and then read it back for further calculation.
The FPU does not use left or right shift, due to way the exponent is stored and is used to signal certain situations like -INF or +INF or NAN (not a number). Also it can't enter number direct into the registers, like the CPU does with things like mov eax, [name], it has to read it from memory or by using the stack. There for it wise to avoid situations where data is moved between the CPU and FPU during calculations. The FPU can only load or use signed integers directly, unsigned integers have to be handled in a different way. For some info about FPU data types and addressing modes look here.
It seems that the speed for FPU division and multiplication can vary if the numbers are simple (single) and/or rounded (like 2 instead of 2.1).
Re: FBMath update
Another slight update, to cope with the release of FB 0.90
http://sourceforge.net/projects/fbmath/ ... mat062.zip
- Precompiled library for FB 0.90 (Windows)
- New compilation scripts for Windows, by S. J. Schaper
- Program MANDELMIX for plotting mixtures of Mandelbrot sets (already published in "Tips and Tricks")
http://sourceforge.net/projects/fbmath/ ... mat062.zip
- Precompiled library for FB 0.90 (Windows)
- New compilation scripts for Windows, by S. J. Schaper
- Program MANDELMIX for plotting mixtures of Mandelbrot sets (already published in "Tips and Tricks")
Re: FBMath update
Thanks jdebord.
I have fb.9 set up now .
I copied your libfbmath.a into the lib folder, and both .bi files into the inc folder.
I tested your mandelmix straight from your demo folder, it works fine.
By the way, I prefer Numerical integration by the Simpson method.
When I was at sea, the old Board of Trade would not accept areas by trapezium method, they were deemed too inaccurate, a no no.
If you used them at exam time you would be sure of fail.
Old habits die hard, here's a Simpson area.
I have fb.9 set up now .
I copied your libfbmath.a into the lib folder, and both .bi files into the inc folder.
I tested your mandelmix straight from your demo folder, it works fine.
By the way, I prefer Numerical integration by the Simpson method.
When I was at sea, the old Board of Trade would not accept areas by trapezium method, they were deemed too inaccurate, a no no.
If you used them at exam time you would be sure of fail.
Old habits die hard, here's a Simpson area.
Code: Select all
sub drawline(x as integer,y as integer,angle as double,length as double,col as uinteger)
angle=angle*atn(1)/45
var x2=x+length*cos(angle)
var y2=y-length*sin(angle)
line(x,y)-(x2,y2),col
end sub
'numerical integrate
'fn=function to integrate
'between a and b
'num=number of ordinates
'AR the returned area
#macro GetArea(fn,a,b,num,AR)
scope
dim as integer n=(num)
if n mod 2=0 then n=n+1 'must be an odd number of ordinates for this Simpson.
Var h=((b)-(a))/n
Var Part1=fn((a)+.5*h)
Var Part2=0.
for i as integer=1 to n-1
Part1 = Part1+fn((a)+h*i+h/2)
Part2 = Part2+fn((a)+h*i)
next i
(AR)= (h/6)*(fn((a))+fn((b))+Part1*4+Part2*2)
end scope
#endmacro
'Actual intregal of 50*sin(x/40)+x/2+50-x^2/2000+200 between L and U
#macro integrate(L,U)
(50*-40*cos((U)/40)+(1/2)*(U)^2/2+50*(U)-(1/2000)*(U)^3/3+200*(U))-_
(50*-40*cos((L)/40)+(1/2)*(L)^2/2+50*(L)-(1/2000)*(L)^3/3+200*(L))
#endmacro
'the function under scrutiny
function wavysine(x as double) as double
return 50*sin(x/40)+x/2+50-x^2/2000+200
end function
dim as double area
screen 20,32
dim as integer mx,my
do
getmouse mx,my
if mx=-1 then mx=0
screenlock
cls
'just sketch the wavysine
for x as single=0 to 1024 step .5
pset(x,wavysine(x)),rgb(255,255,255)
next x
'line from mouse UP and DOWN vertically.
drawline(mx,my,90,800,rgb(255,255,255))
drawline(mx,my,-90,800,rgb(255,255,255))
paint(0,0),rgb(200,0,0),rgb(255,255,255)
GetArea(wavysine,0,mx,40,area) 'by Simpsons Rule using ~40 ordinates
var real=integrate(0,mx) 'by direct calculation
draw string(100,380+50)," Total Screen Area " & 1024*768
draw string(100,400+50),"Simpson's (Red) Area " & (area)
draw string(100,420+50)," Real (Red) Area " & (real)
screenunlock
sleep 1,1
loop until len(inkey)
sleep
Re: FBMath update
Thanks for the Simpson method, dodicat. I will add it.
In my field (pharmacokinetics) the trapezoidal rule is the reference :) But different specialties have different references...
In my field (pharmacokinetics) the trapezoidal rule is the reference :) But different specialties have different references...
Re: FBMath update
I got o:fake errors putting the .a file in the lib dir. I usually have problems with libs.
Runs fine if I put the file with fbmath.bi
I made the mandel mix demo do small tiles here.
Runs fine if I put the file with fbmath.bi
I made the mandel mix demo do small tiles here.
Code: Select all
' =======================================================================
' This program generates fractal images by combining two Mandelbrot sets
' with different exponents, according to the formula:
'
' z(n+1) = (1 - t) * z(n)^p + t * z(n)^q + c
'
' where p and q can be integer or real (>= 2) and t varies from 0 to 1
'
' The method is demonstrated with a series of 28 pictures first published
' by Rollo Silver in his magazine "Amygdala" (#33-34, 1994). It shows the
' combination of the classical Mandelbrot set (p = 2) with the "Multibrot"
' (q = 3) through the progressive increments of t. The pictures are stored
' as BMP files, for a total of about 30 megabytes.
'
' The pictures are colored according to the distance estimator method
' (http://mrob.com/pub/muency/color.html): the value V is computed from
' the distance estimator, while the hue H and saturation S are computed
' from the iteration number.
' =======================================================================
dim shared as string kstr '' new
#include "fbmath.bi"
CONST ScrWidth = 800
CONST ScrHeight = 600
CONST HalfScrWidth = ScrWidth \ 2
CONST HalfScrHeight = ScrHeight \ 2
'CONST Scale = ScrHeight
'' dafhi
const WM=ScrWidth-1
const HM=ScrHeight-1
CONST Esc = 1.0E+10 ' Escape radius (high value required for dist. estimator)
DIM SHARED AS DOUBLE p = 2 ' Exponent (first set)
DIM SHARED AS DOUBLE q = 3 ' Exponent (second set)
DIM SHARED AS DOUBLE p1 ' p - 1
DIM SHARED AS DOUBLE q1 ' q - 1
DIM SHARED AS DOUBLE DistFact = -2 ' Factor for distance estimator
DIM SHARED AS DOUBLE ColorFact = -3 ' Color factor
DIM SHARED AS DOUBLE t ' Proportion of 2nd set
DIM SHARED AS DOUBLE x0 ' Center at (x0, y0)
DIM SHARED AS DOUBLE y0
DIM SHARED AS INTEGER MaxIter ' Max. number of iterations
DIM SHARED AS DOUBLE ZoomFact ' Zoom factor
DIM SHARED AS DOUBLE AbsColor ' Abs(ColorFact)
DIM SHARED AS DOUBLE ScaleFact ' Distance between 2 pixels
DIM SHARED AS DOUBLE Lnp ' ln(p)
DIM SHARED AS DOUBLE LLE ' ln(ln(Esc)) / ln(p)
DIM SHARED AS INTEGER Iter ' Iteration count
DIM SHARED AS DOUBLE Dist ' Distance estimator
DIM SHARED AS DOUBLE Dwell ' Continuous dwell
DIM SHARED AS Complex C_zero = (0, 0)
DIM SHARED AS Complex C_one = (1, 0)
' ********************************************************************
' Subroutines
' ********************************************************************
FUNCTION MdbCol() AS INTEGER
DIM AS DOUBLE DScale, Angle, Radius, P, H, S, V
DIM AS INTEGER R, G, B
DScale = LOG(Dist / ScaleFact) / Lnp + DistFact
IF DScale > 0 THEN
V = 1
ELSEIF DScale > -8 THEN
V = (8 + DScale) / 8
ELSE
V = 0
END IF
P = LOG(Dwell) * AbsColor
IF P < 0.5 THEN
P = 1 - 1.5 * P
Angle = 1 - P
ELSE
P = 1.5 * P - 0.5
Angle = P
END IF
Radius = SQR(P)
IF (Iter MOD 2 = 1) AND (ColorFact > 0) THEN
V = 0.85 * V
Radius = 0.667 * Radius
END IF
H = Angle * 10
H = H - INT(H)
S = Radius - INT(Radius)
HSVtoRGB H * 360, S, V, R, G, B
RETURN RGB(R, G, B)
END FUNCTION
FUNCTION MandelMix (BYVAL xt AS DOUBLE, BYVAL yt AS DOUBLE) AS INTEGER
DIM AS Complex c, z, zp1, zq1, zn, dz, dzn
DIM AS DOUBLE r, Module, LnMod
c = Cmplx(xt, yt)
z = C_zero
dz = C_zero
r = 1 - t
Iter = 0
Module = CAbs(z)
DO WHILE Iter < MaxIter AND Module < Esc
zp1 = r * z ^ p1 ' z^(p-1)
zq1 = t * z ^ q1 ' z^(q-1)
' z(n+1) = (1 - t) * z(n)^p + t * z(n)^q + c
zn = (zp1 + zq1) * z + c
' Derivative : dz(n+1)/dc (used by the distance estimator method)
dzn = (p * zp1 + q * zq1) * dz
Module = CAbs(zn)
z = zn
dz = dzn
Iter = Iter + 1
LOOP
IF Iter = MaxIter THEN RETURN &HFFFFFF
LnMod = LOG(Module)
Dist = p * Module * LnMod / CAbs(dzn)
Dwell = Iter - LOG(LnMod) / Lnp + LLE
RETURN MdbCol()
END FUNCTION
SUB MandelGraph (byval cenx as double=0, byval ceny as double=0, byval _x1 as integer=0,byval _y1 as integer=0, byval _wid as integer=160, byval _hgt as integer=120)
'' sub modified by dafhi for output to screen
if _wid < 1 then exit sub '' comment 4 noobs .. I put _ on names
if _hgt < 1 then exit sub '' for different reasons:
'' 1. usually to prevent name conflict
'' 2. temporary variables
dim as integer bypp,pitch, _wm=_wid-1, _hm=_hgt-1
dim as integer _x2=_x1+_wm, _y2=_y1+_hm, w,h
dim as double _clipx = -_wm/2, _clipy= -_hm/2
if _x1<0 then _clipx += -_x1: _x1=0
if _y1<0 then _clipy += -_y1: _y1=0
if _x2>WM then _x2=WM
if _y2>HM then _y2=HM
_wm=_x2-_x1: _hm=_y2-_y1
_clipx *= ScaleFact: _clipx += cenx
_clipy *= ScaleFact: _clipy += ceny
ScreenInfo w, h, , bypp, pitch: pitch\=bypp
dim as uinteger ptr pixel = screenptr: pixel += _y1 * pitch + _x1
dim as double Ny = _clipy
for pixY as uinteger ptr = pixel to pixel + _HM * pitch step pitch
dim as double Nx = _clipx
for pixX as uinteger ptr = pixY to pixY + _wm
*pixX = MandelMix(Nx, Ny)
Nx += ScaleFact
NEXT: screenlock: screenunlock
kstr=inkey$: if kstr=chr(27) then exit for
Ny+=ScaleFact
NEXT
END SUB
' ********************************************************************
' Main program
' ********************************************************************
CONST Npics = 28
' List of pictures. The descriptions are by Rollo Silver :)
' --------------------------------------------------------------------
' t ZoomFact x0 y0 MaxIter Description
' --------------------------------------------------------------------
DATA 0.10 , 0.5 , -4 , 0 , 100
DATA 0.10 , 2.5 , -8 , 0 , 200 ' The intruder
DATA 0.12 , 0.6 , -3 , 0 , 200 ' The intruder approaches the set
DATA 0.12 , 5 , -6 , 0 , 1000 ' The intruder close-up
DATA 0.12 , 50 , -6.3 , 0.03 , 3000 ' The ugliest fractal
DATA 0.12 , 750 , -6.28735 , 0.0064 , 3000 ' The jet
DATA 0.12 , 3000 , -6.28735 , 0.0064 , 4000 ' The jewel
DATA 0.1258 , 0.65 , -3 , 0 , 500 ' The intruder menaces the set
DATA 0.1258 , 4 , -3.753 , 0 , 2500 ' Imminent collision
DATA 0.12595 , 4 , -3.753 , 0 , 3000 ' The collision begins
DATA 0.126 , 4 , -3.753 , 0 , 3000 ' Compression and inflation
DATA 0.1262 , 4 , -3.753 , 0 , 3000 ' A common head
DATA 0.1266 , 4 , -3.753 , 0 , 3000 ' Coalescence
DATA 0.127 , 4 , -3.753 , 0 , 3000 ' Separation
DATA 0.1285 , 3 , -3.66 , 0.4 , 3000 ' The fountain
DATA 0.1285 , 0.65 , -3 , 0 , 1000 ' The fountain in perspective
DATA 0.14 , 0.65 , -3 , 0 , 1000 ' The major collision impends
DATA 0.142 , 0.65 , -3 , 0 , 1000 ' The major collision
DATA 0.17 , 0.65 , -3 , 0 , 1000 ' Coalescence
DATA 0.22 , 0.65 , -3 , 0 , 1000 ' Buds emerge
DATA 0.28 , 0.65 , -3 , 0 , 1000 ' Ovoid with buds
DATA 0.28 , 4 , -1.093 , -1.591 , 500 ' Rod-like appendage
DATA 0.30 , 4 , -1.011 , -1.591 , 500 ' An ithyphallic appendage
DATA 0.32 , 4 , -0.95 , -1.575 , 500 ' Ithyphallic hypertrophy
DATA 0.40 , 4 , -0.663 , -1.539 , 500 ' The enantiomorphic urge
DATA 0.50 , 4 , -0.37 , -1.446 , 2000 ' Almost complete
DATA 0.50 , 1 , -0.75 , 0 , 1000 ' Almost complete, the big picture
DATA 1.00 , 1.5 , 0 , 0 , 500 ' The cubic Mandelbrot set, complete
'--------------------------------------------------------------------
p1 = p - 1
q1 = q - 1
Lnp = LOG(p)
LLE = LOG(LOG(Esc)) / Lnp
ColorFact = 0.01 * ColorFact
AbsColor = ABS(ColorFact)
SCREENRES ScrWidth, ScrHeight, 32
dim as integer x,y,wid=160,hgt=120, wBy2
FOR I as integer = 1 TO Npics
READ t, Zoomfact, x0, y0, MaxIter
ScaleFact = 6.5 / (ZoomFact*sqr( (wid-1)^2 + (hgt-1)^2 ))
MandelGraph x0,y0, x,y, wid-1,hgt-1 '' -1 for border
x+=wid
if x > WM then
wBy2 = wid/2 - wBy2
x=-wBy2
y+=hgt
EndIf
if kstr=chr(27) then exit for
if y > ScrHeight then exit for
NEXT I
if kstr <> chr(27) then sleep
Re: FBMath update
The lib file is ok here Dafhi, in it's lib folder.
(Win XP, fb .9)
I might as well test out jdebord's library while I'm here.
The only extra functions I had to create were differentiate a polynomial and a Newton method to get the roots of a complex polynomial.
I've not used the root value as such, just the number of iterations to get one.
Sometimes a root is not found within the range of iterations here, but I've noted each case.
TEST FBMATHS ON FB .9:
(Win XP, fb .9)
I might as well test out jdebord's library while I'm here.
The only extra functions I had to create were differentiate a polynomial and a Newton method to get the roots of a complex polynomial.
I've not used the root value as such, just the number of iterations to get one.
Sometimes a root is not found within the range of iterations here, but I've noted each case.
TEST FBMATHS ON FB .9:
Code: Select all
#include "fbmath.bi"
#define Intrange(f,l) int(Rnd*((l+1)-(f))+(f))
'differentiate a complex polynomial
Sub differentiate(pol() As complex,Dpol() As complex)
Redim Dpol(Ubound(pol)-1)
For count As Integer=0 To Ubound(Dpol)
Dpol(count)=(count+1)*pol(count+1)
Next count
End Sub
' A Root of a complex polynomial
Sub GetARoot(f0() As complex,f1() As complex, start As complex,Byref iterations As Integer)
Dim As complex x1=start,x2
For z As Integer=0 To iterations
x2=x1-Cpoly(x1,f0(),Ubound(f0))/Cpoly(x1,f1(),Ubound(f1)) 'Newtons method --> x1=x0-(f(x0)/f'(x0))
If x1=x2 Then iterations=z: Exit For
x1=x2
Next z
End Sub
'======================= GO ========================================
Screenres 1000,700,32
Windowtitle "Press <Esc> to end"
Dim As Integer iterate,startx,finishx,starty,finishy
Dim As Integer r1,r2,r3,r4,i,flag
begin:
startx=0:finishx=100
starty=0:finishy=70
Do
iterate=IntRange(50,80)
'dimension of a polynomial
Dim As Integer dimension=IntRange(3,8)
Do
r1=IntRange(-2,2)
r2=IntRange(-4,4)
r3=IntRange(-12,12)
r4=IntRange(-3,3)
Loop Until r1+r2<>0 And r3+r4<>0
Dim As complex pol(dimension)
'make a polynomial
For z As Integer=0 To dimension
pol(z)=Type<complex>(Rnd*(r1+r2),Rnd*(r3+r4))
Next z
Redim As complex grad(0)
differentiate(pol(),grad())
Var xx=(finishx+startx)\2,yy=(finishy+starty)\2
flag=1
For y As Integer=starty To finishy
For x As Integer=startx To finishx
i=iterate
GetARoot(pol(),grad(),Type<complex>(xx-x,yy-y),i) 'iterate for a root of polynomial
If i<>iterate Then flag=0
Pset(x,y),Rgb(i*128 Mod 255,i*64 Mod 255,i*32 Mod 255)
Next x
Next y
'Case no root found:
If flag Then Draw String(xx-4*12,yy),"GOT NO ROOT.",Rgb(255,0,0)
start:
startx=startx+100:finishx=startx+100
If finishx=1100 And finishy=700 Then
'Screen is Full
Goto begin:
End If
If finishx=1100 Then
startx=-100
starty=starty+70:finishy=starty+70
Goto start
End If
Loop Until Len(Inkey)
Sleep
Re: FBMath update
Decided I didn't like that one very much, so I've adjusted it.
Code: Select all
#include "fbmath.bi"
#define Intrange(f,l) int(Rnd*((l+1)-(f))+(f))
'differentiate a complex polynomial
Sub differentiate(pol() As complex,Dpol() As complex)
Redim Dpol(Ubound(pol)-1)
For count As Integer=0 To Ubound(Dpol)
Dpol(count)=(count+1)*pol(count+1)
Next count
End Sub
' A Root of a complex polynomial
Sub GetARoot(f0() As complex,f1() As complex, start As complex,Byref iterations As Integer)
Dim As complex x1=start,x2
For z As Integer=0 To iterations
x2=x1-Cpoly(x1,f0(),Ubound(f0))/Cpoly(x1,f1(),Ubound(f1)) 'Newtons method --> x1=x0-(f(x0)/f'(x0))
If x1=x2 Then iterations=z: Exit For
x1=x2
Next z
End Sub
'Filter an image
Function Blur(Byref tim As Uinteger Pointer,rad As Single=2,destroy as integer=1) As Uinteger Pointer
Type p2
As Integer x,y
As Uinteger col
End Type
#macro ppoint(_x,_y,colour)
pixel=row+pitch*(_y)+4*(_x)
(colour)=*pixel
#endmacro
#macro ppset(_x,_y,colour)
pixel=row+pitch*(_y)+4*(_x)
*pixel=(colour)
#endmacro
#macro average()
ar=0:ag=0:ab=0:increment=0
xmin=x:If xmin>rad Then xmin=rad
xmax=rad:If x>=(_x-1-rad) Then xmax=_x-1-x
ymin=y:If ymin>rad Then ymin=rad
ymax=rad:If y>=(_y-1-rad) Then ymax=_y-1-y
For y1 As Integer=-ymin To ymax
For x1 As Integer=-xmin To xmax
increment=increment+1
ar=ar+(NewPoints(x+x1,y+y1).col Shr 16 And 255)
ag=ag+(NewPoints(x+x1,y+y1).col Shr 8 And 255)
ab=ab+(NewPoints(x+x1,y+y1).col And 255)
Next x1
Next y1
averagecolour=Rgb(ar/(increment),ag/(increment),ab/(increment))
#endmacro
Dim As Integer _x,_y
Imageinfo tim,_x,_y
Dim As Uinteger Pointer im=Imagecreate(_x,_y)
Dim As Integer pitch
Dim As Any Pointer row
Dim As Uinteger Pointer pixel
Dim As Uinteger col
Imageinfo tim,,,,pitch,row
Dim As p2 NewPoints(_x-1,_y-1)
For y As Integer=0 To (_y)-1
For x As Integer=0 To (_x)-1
ppoint(x,y,col)
NewPoints(x,y)=type<p2>(x,y,col)
Next x
Next y
Dim As Uinteger averagecolour
Dim As Integer ar,ag,ab
Dim As Integer xmin,xmax,ymin,ymax,increment
Imageinfo im,,,,pitch,row
For y As Integer=0 To _y-1
For x As Integer=0 To _x-1
average()
ppset((NewPoints(x,y).x),(NewPoints(x,y).y),averagecolour)
Next x
Next y
if destroy then ImageDestroy tim: tim = 0
Function= im
End Function
sub sorry(im as uinteger ptr)
circle im,(50,70),50,rgb(200,50,0),,,1.5
circle im,(35,60),8,rgb(0,100,0),,,.3,f
circle im,(65,60),8,rgb(0,100,0),,,.3,f
circle im,(35,60),2,rgb(0,0,0),,,.3,f
circle im,(65,60),2,rgb(0,0,0),,,.3,f
circle im,(50,160),70,rgb(255,255,255),1.3,1.8
Draw String im,(50-4*12,130),"GOT NO ROOT!",Rgb(255,0,0)
end sub
'======================= GO ========================================
Screenres 1000,700,32
dim as uinteger ptr im=imagecreate(100,70*2)
Windowtitle "Press <Esc> to end"
Dim As Integer iterate,startx,finishx,starty,finishy
Dim As Integer r1,r2,r3,r4,i,flag
begin:
startx=0:finishx=100
starty=0:finishy=70
Do
iterate=IntRange(50,80)
'dimension of a polynomial
Dim As Integer dimension=IntRange(3,8)
Do
r1=IntRange(-2,2)
r2=IntRange(-4,4)
r3=IntRange(-12,12)
r4=IntRange(-3,3)
Loop Until r1+r2<>0 And r3+r4<>0
Dim As complex pol(dimension)
'make a polynomial
For z As Integer=0 To dimension
pol(z)=Type<complex>(Rnd*(r1+r2),Rnd*(r3+r4))
Next z
Redim As complex grad(0)
differentiate(pol(),grad())
flag=1
For y As Integer=0 To 70*2
For x As Integer=0 To 100
i=iterate
GetARoot(pol(),grad(),Type<complex>(50-x,70-y),i) 'iterate for a root of polynomial
If i<>iterate Then flag=0
Pset im,(x,y),Rgb(i*128 and 255,i*64 and 255,i*32 and 255)
Next x
Next y
'smooth the image
im=blur(im,1)
line im,(0,0)-(100,70*2),rgb(0,0,0),b
'Case no root found:
If flag Then sorry(im)
put(startx,starty),im,pset
start:
startx=startx+100:finishx=startx+100
If finishx=1100 And finishy=700 Then
'Screen is Full
Goto begin:
End If
If finishx=1100 Then
startx=-100
starty=starty+70*2:finishy=starty+70*2
Goto start
End If
Loop Until Inkey=chr(27)
Sleep
imagedestroy im