FBMath update

Headers, Bindings, Libraries for use with FreeBASIC, Please include example of use to help ensure they are tested and usable.
jdebord
Posts: 547
Joined: May 27, 2005 6:20
Location: Limoges, France
Contact:

FBMath update

Post by jdebord »

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
frisian
Posts: 249
Joined: Oct 08, 2009 17:25

Re: FBMath update

Post by frisian »

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.

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
Looks like that
DIM R() AS DOUBLE
DIM C() AS Complex

is missing.
jdebord
Posts: 547
Joined: May 27, 2005 6:20
Location: Limoges, France
Contact:

Re: FBMath update

Post by jdebord »

These arrays are dimensioned by REDIM inside the subroutine SortRoots. This does not seem to give problems with FB 0.24

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
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
frisian
Posts: 249
Joined: Oct 08, 2009 17:25

Re: FBMath update

Post by frisian »

jdebord wrote:These arrays are dimensioned by REDIM inside the subroutine SortRoots. This does not seem to give problems with FB 0.24
Your right, I messed up, did not check if older FB versions worked correctly.
The problem is with the new FB version.
dkl
Site Admin
Posts: 3235
Joined: Jul 28, 2005 14:45
Location: Germany

Re: FBMath update

Post by dkl »

that looks like it was due to this bug: #654
counting_pine
Site Admin
Posts: 6323
Joined: Jul 05, 2005 17:32
Location: Manchester, Lancs

Re: FBMath update

Post by counting_pine »

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 my understanding of Redim, this is what is going on:
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.
frisian
Posts: 249
Joined: Oct 08, 2009 17:25

Re: FBMath update

Post by frisian »

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”

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
FB swap works also on type's so, swap x,y gives the same result.

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
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”

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 >
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.
jdebord
Posts: 547
Joined: May 27, 2005 6:20
Location: Limoges, France
Contact:

Re: FBMath update

Post by jdebord »

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.
frisian
Posts: 249
Joined: Oct 08, 2009 17:25

Re: FBMath update

Post by frisian »

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.
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: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).
jdebord
Posts: 547
Joined: May 27, 2005 6:20
Location: Limoges, France
Contact:

Re: FBMath update

Post by jdebord »

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")
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: FBMath update

Post by dodicat »

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.

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
     
jdebord
Posts: 547
Joined: May 27, 2005 6:20
Location: Limoges, France
Contact:

Re: FBMath update

Post by jdebord »

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...
dafhi
Posts: 1641
Joined: Jun 04, 2005 9:51

Re: FBMath update

Post by dafhi »

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.

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
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: FBMath update

Post by dodicat »

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:

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
 
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: FBMath update

Post by dodicat »

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
 
 
Post Reply