## [solved]Math Guru II How to make a Normal little bit random FAST?

For other topics related to the FreeBASIC project or its community.
gothon
Posts: 224
Joined: Apr 11, 2011 22:22

### Re: Math Guru II How to make a Normal little bit random FAST?

D.J.Peters wrote:@dafhi @gothon first you think it's clever but in practice it looks different.

If you get any surface normal (it points in any direction)

1) how you would select a precalculated random vector ?
(as a minimum you have to get the right qudrant from sphere to chose one)

I suggested pre-computing a random low angle rotation matrix not a vector. A simple way to do this might be to make a small random rotation about the X, Y and Z axis, then multiplying all three rotations into one 3x3 matrix.

https://en.wikipedia.org/wiki/Rotation_ ... _rotations

Though there are probably other methods as-well since you might find a method giving you a more desirable distribution of rotations.

D.J.Peters wrote:2) more important how to rotate any new random vector does it fits the cone ?
(I think a hour about it but I don't get it)

Once you have a rotation matrix, you simply multiply it by your vector to apply the rotation.

Something like:
U.X = V.X * m_00 + V.Y * m_01 + V.Z * m_02
U.Y = V.X * m_10 + V.Y * m_11 + V.Z * m_12
U.Z = V.X * m_20 + V.Y * m_21 + V.Z * m_22

D.J.Peters wrote:Show me the code if you like please

Yep, real code is best, but hold on as it will take me a bit of time to write it and debug it. :)
gothon
Posts: 224
Joined: Apr 11, 2011 22:22

### Re: Math Guru II How to make a Normal little bit random FAST?

Ok here is some code

Code: Select all

`Type Vect3D    As Single X, Y, ZEnd TypeType Mat3x3    M(2, 2) As SingleEnd Type'Multiply two 3x3 Matricies possably backwards? :PFunction MatMult3x3(A As Mat3x3, B As Mat3x3) As Mat3x3    Dim C As Mat3x3        For I As Integer = 0 To 2        For J As Integer = 0 To 2            For K As Integer = 0 To 2                C.M(I, J) += A.M(I, K) * B.M(K, J)            Next K        Next J    Next I        Return CEnd Function'Fill parameters with a 2D rotation matrixSub RotMat2D(Angle As Double, ByRef M_00 As Single, ByRef M_01 As Single, ByRef M_10 As Single, ByRef M_11 As Single)    M_00 = Cos(Angle): M_01 = -Sin(Angle) '<- Trig Functions    M_10 = -M_01: M_11 = M_00End SubDim Normals() As Vect3DReDim Normals(5000000 - 1)Dim Low_Angle_Rot(255, 9) As Mat3x3' Precompute Rotation MatriciesFor I As Integer = 0 To UBound(Low_Angle_Rot, 1)    For J As Integer = 0 To UBound(Low_Angle_Rot, 2)        Dim As Mat3x3 Rx, Ry, Rz, R_Tmp        Rx.M(0, 0) = 1        RotMat2D (J + 1) / 10 * (Rnd() - Rnd()), Rx.M(1, 1), Rx.M(1, 2), Rx.M(2, 1), Rx.M(2, 2)        Ry.M(1, 1) = 1        RotMat2D (J + 1) / 10 * (Rnd() - Rnd()), Ry.M(0, 0), Ry.M(0, 2), Ry.M(2, 0), Ry.M(2, 2)        Rz.M(2, 2) = 1        RotMat2D (J + 1) / 10 * (Rnd() - Rnd()), Rz.M(0, 0), Rz.M(0, 1), Rz.M(1, 0), Rz.M(1, 1)        R_Tmp = MatMult3x3(Ry, Rx)        Low_Angle_Rot(I, J) = MatMult3x3(Rz, R_Tmp)    Next JNext I' Initialize normal vectors to something random but length 1For I As Integer = 0 To UBound(Normals)    Normals(I).X = Rnd * 2 - 1    Normals(I).Y = Rnd * 2 - 1    Normals(I).Z = Rnd * 2 - 1    Var L = Sqr(Normals(I).X ^ 2 + Normals(I).Y ^ 2 + Normals(I).Z ^ 2)    Normals(I).X /= L    Normals(I).Y /= L    Normals(I).Z /= LNext I#Define SMALL_RND (rnd()-rnd()) * surface_roughnessVar surface_roughness = 0.1Var Start = Timer 'Randomize and RenormalizeFor I As Integer = 0 To UBound(Normals)    Normals(I).X += SMALL_RND    Normals(I).Y += SMALL_RND    Normals(I).Z += SMALL_RND    Var L = Normals(I).X*Normals(I).X + Normals(I).Y*Normals(I).Y + Normals(I).Z*Normals(I).Z    L=1/Sqr(L) : Normals(I).X*=L : Normals(I).Y*=L : Normals(I).Z*=LNext IVar T1 = Timer - StartDim As Double E = 0 'Determine how close to 1 the length isFor I As Integer = 0 To UBound(Normals)    Var L = Normals(I).X*Normals(I).X + Normals(I).Y*Normals(I).Y + Normals(I).Z*Normals(I).Z    E += Abs(1 - L)Next IPrintPrint "Time to randomize and re-normalize vectors: " & T1Print "Average Float Error: " & E / (UBound(Normals) + 1)Start = Timer 'Multiply random rotation matriciesFor I As Integer = 0 To UBound(Normals)    Dim J As Integer = Int(Rnd * (UBound(Low_Angle_Rot) + 1))    Dim N As Vect3D    With Low_Angle_Rot(J, Int(surface_roughness* 10))        N.X = Normals(I).X * .M(0,0) + Normals(I).Y * .M(0,1) + Normals(I).Z * .M(0,2)        N.Y = Normals(I).X * .M(1,0) + Normals(I).Y * .M(1,1) + Normals(I).Z * .M(1,2)        N.Z = Normals(I).X * .M(2,0) + Normals(I).Y * .M(2,1) + Normals(I).Z * .M(2,2)    End With    Normals(I).X = N.X    Normals(I).Y = N.Y    Normals(I).Z = N.ZNext IVar T2 = Timer - StartE = 0 'Determine how close to 1 the length isFor I As Integer = 0 To UBound(Normals)    Var L = Normals(I).X*Normals(I).X + Normals(I).Y*Normals(I).Y + Normals(I).Z*Normals(I).Z    E += Abs(1 - L)Next IPrintPrint "Time to rotate vectors: " & T2Print "Average Float Error: " & E / (UBound(Normals) + 1)PrintPrint "Speedup factor: " & T1 / T2Sleep`

I get a speedup factor of around 2-3 on my computer, although that may also have something to do with there being only 1 Rnd() call in the loop instead of 6. I had to discretize surface_roughness to 10 levels so I could have different levels pre-computed. Though I just let it be 0.1 since there is no source of that information in this simple demonstration anyway.
D.J.Peters
Posts: 7838
Joined: May 28, 2005 3:28

### Re: Math Guru II How to make a Normal little bit random FAST?

@gothon I know how to use matrices I do every day and thank for your speed compare

Do you sure your random rotation forms a small cone around the original normal ?
(I will test if i'm back home)

Init 4 normals with for example an up vector 0,1,0 rotate it and plot it on screen (you can ignore z for this short test)
does it looks like a cone ?

Joshy Last edited by D.J.Peters on Sep 25, 2017 21:40, edited 1 time in total.
dodicat
Posts: 5988
Joined: Jan 10, 2006 20:30
Location: Scotland

### Re: Math Guru II How to make a Normal little bit random FAST?

Try making the unit vectors in situ
(Excuse my Latin)
The time killer is getting a random + or -
I use doubles to check errors.

Code: Select all

`screen 19,32type pt    as double x,y,zend typedim as long max=5000000redim as pt p(1 to max)redim as byte sign(1 to max+3)'get random signsdo     dim as double t=timerfor n as long=1 to max+3    sign(n)=IIf(Rnd>.5,1,-1)    next'dim as double t=timerfor n as long=1 to max 'create unit vectors on the flyp(n)=type(sign(n)*rnd , _            sign(n+1)*rnd*sqr(1-p(n).x*p(n).x), _          sign(n+2)*sqr((1-(p(n).x*p(n).x+p(n).y*p(n).y))))next nprint "Time taken ",timer-t,"  Press a key"'get error residuedim as single Rfor n as long=1 to max   r+=abs( 1-sqr(p(n).x^2+p(n).y^2+p(n).z^2))     next nprint "ERROR ";rdim as long ctrx,ctry,ctrzfor z as long=lbound(p) to ubound(p)if (sqr(p(z).x^2+p(z).y^2+p(z).z^2) -1) >1e-5 then print z, sqr(p(z).x^2+p(z).y^2+p(z).z^2)'extra ckeckif p(z).x>0 then ctrx+=1if p(z).y>0 then ctry+=1if p(z).z>0 then ctrz+=1next zprint max/2, "Expected for each"print ctrx,ctry,ctrz, "Got"draw string(350,150),"Shade by .z"draw string(350+250,150),"Bare bones"Circle(400,300),100,4For z As Long=LBound(p) To UBound(p)       PSet(400+100*p(z).x,300+100*p(z).y),RGB(Abs(p(z).z)*255,Abs(p(z).z)*255,Abs(p(z).z)*25)    PSet(400+250+100*p(z).x,300+100*p(z).y)Nextprint "press a key"sleepclsloop until inkey=chr(27) `
Tourist Trap
Posts: 2768
Joined: Jun 02, 2015 16:24

### Re: Math Guru II How to make a Normal little bit random FAST?

dodicat wrote:(Excuse my Latin)

Hi, dodicat,

I'm lost with your perspective. I don't know how to zoom/unzoom, or whatever (rotate view) with the perspective. I would need it to test some rodrigues' rotation (axis rotation). I have some problem obtaining the correct rotation, and to vizualise the things would help.

Code: Select all

`'attempt of application of the Rodrigues rotation'_________________________the dodi's perspectivetype P3D    as single   x,y,zend typefunction Perspective(P as P3D, Eyepoint as P3D) as P3D    dim as single   w => 1 + ( P.z/Eyepoint.z )    '    return type<P3D>( (P.x-Eyepoint.x)/w + Eyepoint.x, _                       (P.y-Eyepoint.y)/w + Eyepoint.y, _                       (P.z-Eyepoint.z)/w + Eyepoint.z )end function'_________________________data'origin h of the rotation axis (U)dim as integer   hx = 295dim as integer   hy = 400dim as integer   hz = 160dim as P3D h    h.x = hx    h.y = hy    h.z = hz'unit direction vector of axis (U)dim as single   a   =   0dim as single   b   =   100dim as single   c   =   0dim as single   norm    = sqr(a^2 + b^2 + c^2)a = a/normb = b/normc = c/normdim as P3D u    u.x = a    u.y = b    u.z = c'point M to apply rotation ondim as integer   mx   = 225dim as integer   my   = 200dim as integer   mz   = 160dim as P3D M    M.x = mx    M.y = my    M.z = mz'transformdim as P3D MMdim as P3D MM2'angle of the rotationdim as double   axisRotationAngle'the Rodrigues transformfunction RodriguesRotation(H as P3D, U as P3D, M as P3D, angle as double) as P3D    dim as P3D    result    '    'vector position of M relative to H    dim as P3D G        G.x = M.x - H.x        G.y = M.y - H.y        G.z = M.z - H.z    '    'K = 1 - cos angle    dim as double   K = 1 - cos(angle)    '    dim as single aam = U.x*(U.x*G.x + U.y*G.y + U.z*G.z)    dim as single bbm = U.y*(U.x*G.x + U.y*G.y + U.z*G.z)    dim as single ccm = U.z*(U.x*G.x + U.y*G.y + U.z*G.z)    '    result.x = H.x + G.x*( U.x*U.x*K + cos(angle) ) + _                      G.y*( U.x*U.y*K - U.z*sin(angle) ) + _                      G.z*( U.x*U.z*K + U.y*sin(angle) )    result.y = H.y + G.x*( U.x*U.y*K + U.z*sin(angle) ) + _                      G.y*( U.y*U.y*K + cos(angle) ) + _                      G.z*( U.y*U.z*K - U.x*sin(angle) )    result.z = H.z + G.x*( U.x*U.z*K - U.y*sin(angle) ) + _                      G.y*( U.y*U.z*K + U.x*sin(angle) ) + _                      G.z*( U.z*U.z*K + cos(angle) )   '    return type<P3D>( result.x , result.y, result.z)end function'---------------------------------------_________'---------------------------------------_________screenRes 800, 550, 32dim as P3D      pts(1 to 4)dim as P3D      eyepointdim as double   angledim as integer  gmx,gmy,wheeldo    angle += 0.04    if angle>7 then angle = 0        getmouse gmx,gmy,wheel        eyepoint.x = (gmx - 100)    eyepoint.y = 1.2*(gmy) + 20    eyepoint.z = 455 - wheel*10       pts(1) = type(50,400,0)    pts(2) = type(540,400,0)    pts(3) = type(295,400,160)    pts(4) = type(295,200,160)       MM      => RodriguesRotation(H, U, M, angle)        ''****    MM2     => RodriguesRotation(H, U, M, 2*angle)      ''****    screenlock    cls        circle (eyepoint.x, eyepoint.y - 20), 1                for i as long=1 to 20             var p1 => Perspective(pts(1), eyepoint)            var p2 => Perspective(pts(2), eyepoint)            pts(1).z += 20            pts(2).z += 20            line(p1.x, p1.y)-(p2.x, p2.y)        next i                var p3  => perspective(pts(3), eyepoint)        var p4  => perspective(pts(4), eyepoint)        var pH  => perspective(H, eyepoint)        var pM => perspective(M, eyepoint)        var pMM => perspective(MM, eyepoint)        var pMM2 => perspective(MM2, eyepoint)                circle (p3.x, p3.y), 4        circle (p4.x, p4.y), 4        circle (pM.x, pM.y), 4, rgb(200,0,200), , , , f                circle (pMM.x, pMM.y), 4, rgb(100,200,100), , , , f        MM = RodriguesRotation(H, U, M, angle + 0.01)        MM2 = RodriguesRotation(H, U, M, angle + 0.01)                circle (pMM.x, pMM.y), 2, rgb(100,200,100), , , , f        circle (pMM2.x, pMM2.y), 2, rgb(200,200,100), , , , f        MM = RodriguesRotation(H, U, M, angle + 0.02)                pMM=perspective(MM,eyepoint)        circle (pMM.x, pMM.y), 2, rgb(100,200,100), , , , f                pMM=perspective(M,eyepoint)        MM = RodriguesRotation(H, U, M, angle)        circle (pH.x, pH.y), 4, rgb(100,100,200)        line(p3.x,p3.y) - (p4.x,p4.y)        line(int(pMM.x), int(pMM.y))-(int(pH.x), int(pH.y))            screenunlock    sleep 15loop until len(inkey)sleep'(eof)`
dafhi
Posts: 1258
Joined: Jun 04, 2005 9:51

### Re: Math Guru II How to make a Normal little bit random FAST?

Finally got it done :-)

[update] translate() was not correct

Code: Select all

`...`
Last edited by dafhi on Aug 23, 2016 7:36, edited 3 times in total.
Tourist Trap
Posts: 2768
Joined: Jun 02, 2015 16:24

### Re: Math Guru II How to make a Normal little bit random FAST?

dafhi wrote:Finally got it done :-)

Very nice rendering.
dafhi
Posts: 1258
Joined: Jun 04, 2005 9:51

### Re: Math Guru II How to make a Normal little bit random FAST?

Thanks man :-)
D.J.Peters
Posts: 7838
Joined: May 28, 2005 3:28

### Re: Math Guru II How to make a Normal little bit random FAST?

@dafhi
i'm back at home and tested your code looks ok
good job man

Joshy
dodicat
Posts: 5988
Joined: Jan 10, 2006 20:30
Location: Scotland

### Re: Math Guru II How to make a Normal little bit random FAST?

Yea, nice Dafhi.
I cannot fill 5 million random unit vectors in any less than about 1/5 second on 32 bit.
64 bit is a tad faster
Tried the asm sqrt, but it is slower and only 32 bit.
randomize,2 seems the fastest rnd.
A lookup table for rnd, and use only once speeds it up a bit.

Code: Select all

`randomize ,2screen 19,32type pt    as double x,y,zend typedim as long max=5000000redim as pt p(1 to max)redim as byte sign(1 to max+3)'get random signsredim as double rnds(1 to max)'create a look up tablefor n as long=1 to max    rnds(n)= rnd-rndnext'===============================================do     dim as double t=timerfor n as long=1 to max+3     sign(n)=IIf(Rnd(n)>.5,1,-1)    next''dim as double t=timer'======================  CREATE THE UNIT VECTORS ==============for n as long=1 to max ' on the flyp(n)=type(sign(n)*rnds(n) , _    'use the look up table here only          sign(n+1)*rnd*sqr(1-p(n).x*p(n).x), _          sign(n+2)*sqr((1-(p(n).x*p(n).x+p(n).y*p(n).y))))next nprint "Time taken ",timer-t'get error residuedim as double Rfor n as long=1 to max   r+=abs( 1-sqr(p(n).x^2+p(n).y^2+p(n).z^2))     next nprint "ERROR ";rdim as long ctrx,ctry,ctrzfor z as long=lbound(p) to ubound(p)if (sqr(p(z).x^2+p(z).y^2+p(z).z^2) -1) >1e-5 then print z, sqr(p(z).x^2+p(z).y^2+p(z).z^2)'extra ckeckif p(z).x>0 then ctrx+=1if p(z).y>0 then ctry+=1if p(z).z>0 then ctrz+=1next zprint max/2, "Expected for each"print ctrx,ctry,ctrz, "Got"draw string(350,150),"Shade by .z"draw string(350+250,150),"Bare bones"Circle(400,300),100,4For z As Long=LBound(p) To UBound(p)       PSet(400+100*p(z).x,300+100*p(z).y),RGB(Abs(p(z).z)*255,Abs(p(z).z)*255,Abs(p(z).z)*25)    PSet(400+250+100*p(z).x,300+100*p(z).y)Nextprint "press a key"sleepclsloop until inkey=chr(27) `
Tourist Trap
Posts: 2768
Joined: Jun 02, 2015 16:24

### Re: Math Guru II How to make a Normal little bit random FAST?

A little off-topic now, but even if can't be displayed in perspective, here is the Rodrigues rotation formula, both normal and for small angles. It seems well working. It may serve somewhere some day.

Code: Select all

`#include "fbgfx.bi"'useless here!nameSpace cst    const as double _pi => 4*atn(1)    dim shared as double    presin(0 to 255)    dim shared as double    precos(0 to 255)    dim shared as double    cycleStep => 1*256/360end nameSpacesub PrecalcultaSinCos()    dim as ubyte    currentStep    dim as ubyte    furtherStep    do        currentStep=furtherStep        cst.presin(currentStep) = sin(currentStep)        cst.precos(currentStep) = cos(currentStep)        furtherStep += cst.cycleStep    loop until furtherStep<=currentStepend sub'starts here...type P3D    declare property PerspFactor() as double    declare property PerspFactor(byval as double)    declare property X2D() as integer    declare property Y2D() as integer    declare sub Rotate(byval AngleAroundX as double, _                        byval AngleAroundY as double, _                        byval AngleAroundZ as double)    declare sub DrawPoint(byval Colour as ulong=rgb(200,175,130))        as integer  _x        as integer  _y        as integer  _z        as double   _perspFactor    static as integer   _xPersp    static as integer   _yPersp    static as integer   _zPerspend typedim as integer   P3D._xPerspdim as integer   P3D._yPerspdim as integer   P3D._zPerspproperty P3D.PerspFactor() as double    '---->    return THIS._perspFactorend propertyproperty P3D.PerspFactor(byval SetValue as double)    THIS._perspFactor = SetValueend propertyproperty P3D.X2D() as integer    '---->    return ( THIS._x )end propertyproperty P3D.Y2D() as integer    '---->    return ( THIS._y )end propertysub P3D.Rotate(byval AngleAroundXaxis as double, _                byval AngleAroundYaxis as double, _                byval AngleAroundZaxis as double)    THIS._y = THIS._y*cos(AngleAroundXaxis) - THIS._z*sin(AngleAroundXaxis)    THIS._z = THIS._y*sin(AngleAroundXaxis) + THIS._z*cos(AngleAroundXaxis)    THIS._x = THIS._x*cos(AngleAroundYaxis) - THIS._z*sin(AngleAroundYaxis)    THIS._z = THIS._x*sin(AngleAroundYaxis) + THIS._z*cos(AngleAroundYaxis)    THIS._x = THIS._x*cos(AngleAroundZaxis) - THIS._y*sin(AngleAroundZaxis)    THIS._y = THIS._x*sin(AngleAroundZaxis) + THIS._y*cos(AngleAroundZaxis)end subsub P3D.DrawPoint(byval Colour as ulong=rgb(200,175,130))    circle (THIS.X2D, THIS.Y2D), 8 - (6 - 6*THIS._z/600), Colour, , , , fend sub    function RodriguesRotation(byref H as P3D, _                            byref U as P3D, _                            byref M as P3D, _                            byval angle as double) _                            as P3D    '_the Rodrigues transform_    '    dim as P3D    result    '    'vector position of M relative to H    dim as P3D G        G._x = M._x - H._x        G._y = M._y - H._y        G._z = M._z - H._z    '    'K = 1 - cos angle    dim as double   K = 1 - cos(angle)    '    result._x = H._x + G._x*( U._x*U._x*K + cos(angle) ) + _                        G._y*( U._x*U._y*K - U._z*sin(angle) ) + _                        G._z*( U._x*U._z*K + U._y*sin(angle) )    result._y = H._y + G._x*( U._x*U._y*K + U._z*sin(angle) ) + _                        G._y*( U._y*U._y*K + cos(angle) ) + _                        G._z*( U._y*U._z*K - U._x*sin(angle) )    result._z = H._z + G._x*( U._x*U._z*K - U._y*sin(angle) ) + _                        G._y*( U._y*U._z*K + U._x*sin(angle) ) + _                        G._z*( U._z*U._z*K + cos(angle) )    '    '---->    return resultend functionfunction SmallRodriguesRotation(byref H as P3D, _                                 byref U as P3D, _                                 byref M as P3D, _                                 byval angle as double) _                                 as P3D    '_the Rodrigues transform_    '    dim as P3D    result    '    'vector position of M relative to H    dim as P3D G        G._x = M._x - H._x        G._y = M._y - H._y        G._z = M._z - H._z    '    'K = 1 - cos angle    dim as double   K = 0    '    result._x = H._x + G._x*( U._x*U._x*K + 1 ) + _                        G._y*( U._x*U._y*K - U._z*(angle) ) + _                        G._z*( U._x*U._z*K + U._y*(angle) )    result._y = H._y + G._x*( U._x*U._y*K + U._z*(angle) ) + _                        G._y*( U._y*U._y*K + 1 ) + _                        G._z*( U._y*U._z*K - U._x*(angle) )    result._z = H._z + G._x*( U._x*U._z*K - U._y*(angle) ) + _                        G._y*( U._y*U._z*K + U._x*(angle) ) + _                        G._z*( U._z*U._z*K + 1 )    '    '---->    return resultend function'- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - .'.- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  -randomize TIMERscreenRes 800, 540, 32, 2screenSet 1, 0color , rgb(45,58,72)clsPrecalcultaSinCos()'_____________________________data____'origin h of the rotation axis (U)dim as single   hx = 400dim as single   hy = 250dim as single   hz = 0dim as P3D h    h._x = hx    h._y = hy    h._z = hz'direction unit vector of the axis (U)dim as single   a   =   0dim as single   b   =   0dim as single   c   =   100dim as single   norm    = sqr(a^2 + b^2 + c^2)dim as P3D abc    abc._x = a    abc._y = b    abc._z = cdim as P3D u    u._x = a/norm    u._y = b/norm    u._z = c/norm'point M to apply rotation ondim as single shifter = 10dim as single   mx   = 400 - shifterdim as single   my   = 250dim as single   mz   = 160dim as P3D M    m._x = mx    m._y = my    m._z = mz'angle of the rotationdim as double   axisRotationAngle'_____________________________________'- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - _dim as P3D      randomSelection(0 to 260)dim as double   randomAngle(0 to 260)dim as integer  selectionIndexdim as P3D      rotationResultPointdim as integer  gmX, gmY, gmWheel, gmBtndim as double   axisRotationDecorateAngledo    getMouse gmX, gmY, gmWheel, gmBtn    line (0, 250)-(799, 250)        m.DrawPoint(rgb(100,100,80))    h.DrawPoint(rgb(140,210,180))    rotationResultPoint = RodriguesRotation(h, u, m, axisRotationDecorateAngle)    rotationResultPoint.DrawPoint()    rotationResultPoint = RodriguesRotation(h, u, m, axisRotationDecorateAngle + 0.1)    rotationResultPoint.DrawPoint()    rotationResultPoint = RodriguesRotation(h, u, m, axisRotationDecorateAngle + 0.2)    rotationResultPoint.DrawPoint()    rotationResultPoint = RodriguesRotation(h, u, m, axisRotationDecorateAngle + 0.3)    rotationResultPoint.DrawPoint()    rotationResultPoint = RodriguesRotation(h, u, m, axisRotationDecorateAngle + 0.4)    rotationResultPoint.DrawPoint(rgb(200,100,130))    if shifter<2600 then        randomAngle(selectionIndex) = (selectionIndex/200 + rnd())*cst._pi/8        randomSelection(selectionIndex) = SmallRodriguesRotation(h, u, m, randomAngle(selectionIndex))        for i as integer = 0 to selectionIndex            randomSelection(i).DrawPoint(rgb(200,100,130))          next i    else        for i as integer = 0 to 260            randomSelection(i).DrawPoint(rgb(100,200,130))          next i    end if    draw string (10,10), "current angle  --> "& left(str(axisRotationAngle), 4)    draw string (10,20), "current radius --> "& str(shifter)    draw string (10,30), "current index --> "& str(selectionIndex)        screenSync()    flip 1, 0    cls    if shifter<2600 then        shifter += 10                   ''increase orbit radius        selectionIndex += 1        axisRotationAngle += cst._pi/8        if axisRotationAngle>2*cst._pi then axisRotationAngle = 0        m._x += shifter/100*(cos(axisRotationAngle) + sin(axisRotationAngle))        m._y -= shifter/100*(sin(axisRotationAngle) - cos(axisRotationAngle))        m._z = rnd()*180 + 10    else        m._x = mx - 120        m._y = my        m._z = mz    end if        axisRotationDecorateAngle += 0.01    if axisRotationDecorateAngle>2*cst._pi then axisRotationDecorateAngle = 0        sleep 15loop until inkey()=chr(27)'.- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  -'- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - .getKey()'(eof)`
dafhi
Posts: 1258
Joined: Jun 04, 2005 9:51

### Re: Math Guru II How to make a Normal little bit random FAST?

@dodicat - I create a set of random points. After that I choose randomly 5 mill times

Code: Select all

`'  subs Cone_Precalc and translate() have all the detailstype imagevars ' 2016 Jan 31  as any ptr              im, pixels  as integer              w,h,bpp,bypp,pitch,numpages,flags,rate, is_screen, pitchBy, wm,hm  as single               midx,midy, midxm, midym, diagonal  as string               driver_name  declare sub             scr_inf  Declare Sub             screen_init(wid As single=-1, hgt As single=-1, bpp as UInteger=32, numPages as integer=1, Flags as integer=0)  Declare Function        create(wid As long = -1, hgt As long=0, color As Ulong=RGB(0,0,0)) As Any ptr  declare sub             cls(col as ulong=&HFF000000)  declare sub             checkers(pColor As ULong=RGBA(145,145,145,255), size As UInteger = 12)  Declare Destructor private:  Declare Sub             destroy  declare sub             vars_commonend typeDestructor imagevars:  Destroy:  End DestructorSub imagevars.Destroy():  If im <> 0 Then ImageDestroy im: im = 0: endif:  End SubSub imagevars.cls(col as ulong):  line (0,0)-(wm,hm),col,bf:  end subsub imagevars.vars_common  wm=w-1: midx=w/2: midxm = wm/2  hm=h-1: midy=h/2: midym = hm/2  diagonal=sqr(w*w+h*h): if bypp<>0 then pitchBy=pitch\byppend subsub imagevars.scr_inf  ScreenInfo w,h, bpp, bypp, pitch, rate, driver_name:  pixels = ScreenPtr  vars_commonend subSub imagevars.screen_init(_w As single, _h As single, _bpp as UInteger, _numpages as integer, _flags as integer)  dim as integer ww,hh  ScreenInfo ww,hh  _w=abs(_w): if _w<=1 then _w*=ww: _flags or= 8 'borderless  _h=abs(_h): if _h<=1 then _h*=hh: _flags or= 8  numpages=_numpages: flags=_flags  Destroy: ScreenRes _w,_h,_bpp,numpages,flags: scr_inf: is_screen = -1  if numPages > 1 then screenset 0,1End subFunction imagevars.create(wid As long, hgt As long, col As ULong) As Any Ptr  if hgt=0 then scr_inf: wid=w: hgt=h  Destroy: im = ImageCreate( wid, hgt, col )  ImageInfo im, w, h, bypp, pitch, pixels:  bpp = bypp * 8  vars_common: is_screen = 0: Return imEnd FunctionSub imagevars.checkers(pColor As ULong, size As UInteger)  Dim As UInteger SizeDouble=size*2,SizeM=size-1  For Y as integer = 0 To hm Step size    For X as integer = -size * ((Y/SizeDouble)=Int(Y/SizeDouble)) To wm Step SizeDouble      Line im,(X,Y)-(X+SizeM,Y+SizeM),pColor,BF    Next  NextEnd Sub#Macro Alpha256(ret,back, fore, am, a256)  ret=((_  (fore And &Hff00ff) * a256 + _  (back And &Hff00ff) * am + &H800080) And &Hff00ff00 Or (_  (fore And &H00ff00) * a256 + _  (back And &H00ff00) * am + &H008000) And &H00ff0000) Shr 8#EndMacrotype aadot  as single               rad = 0.65, alpha = 1  as ulong                outlined = 0  as imagevars ptr        p  as any ptr              pixels  declare sub             render_target(pti as imagevars ptr)  declare sub             draw(x as single=0, y as single=0, col as ulong=&HFFFFFFFF) private:  as single               slope, slope_X2end typesub aadot.render_target(pti as imagevars ptr)  if pti->is_screen then: pixels = ScreenPtr  else: pixels = pti->pixels  end if:  p = ptiend subsub aadot.draw(x as single, y as single, col as ulong)  if p->h < 1 then exit sub  slope = alpha*256: slope_X2 = slope*2  dim as single           cone_hgt = rad * slope  dim as integer          _X = x-rad, X2 = x+rad:  _X += _X * (_X < 0)  dim as integer          _Y = y-rad, Y2 = y+rad:  _Y += _Y * (_Y < 0)  dim as single           dxClip = slope * (x - _X), dy = slope * (y - _Y)  X2 += (X2 - p->wm) * (X2 > p->wm)  Y2 += (Y2 - p->hm) * (Y2 > p->hm)  dim as integer          pwidm = X2 - _X  dim as ulong ptr        pCorner = pixels:  pCorner += _Y * p->pitchBy + _X  if outlined then    for py as ulong ptr = pCorner to @pCorner[ (Y2-_Y)*p->pitchBy ] step p->pitchBy      dim as single ySq = dy*dy, dx = dxClip      for px as ulong ptr = py to py + pwidm        dim as integer  alph = cone_hgt - sqr(dx*dx + ySq)        if alph > slope then alph = slope_X2 - alph        if alph > 0 then          dim as integer  AlphM = 256 - alph          Alpha256( *px, *px, col, AlphM, alph )        end if        dx -= slope      next:  dy -= slope    next  else    for py as ulong ptr = pCorner to @pCorner[ (Y2-_Y)*p->pitchBy ] step p->pitchBy      dim as single ySq = dy*dy, dx = dxClip      for px as ulong ptr = py to py + pwidm        dim as integer  alph = cone_hgt - sqr(dx*dx + ySq)        if alph > slope then:  alph = slope        else:  alph += alph * (alph < 0)        end if:  dim as integer  AlphM = 256 - alph        Alpha256( *px, *px, col, AlphM, alph )        dx -= slope      next:  dy -= slope    next  end ifend sub#ifndef PiConst                   TwoPi = 8 * atn(1)Const                   pi    = 4 * Atn(1)#EndIfType v3  As double             x,y,z  as double             rad  as long               colEnd TypeType Axis3D  As v3                 AxisX=(1,0,0), AxisY=(0,-1,0), AxisZ=(0,0,1)  As double             x,y,zEnd Type#Macro RotC_(dsta,dstb,srca,srcb)    temp = cosa_*srca - sina_*srcb    dstb = cosa_*srcb + sina_*srca    dsta = temp#EndMacro#Macro RotC(a_,dst,src,dota,dotb)  Scope    Dim As single a = (a_)    Dim As double cosa_ = Cos(a), sina_ = Sin(a), temp    RotC_( (dst.AxisX)dota, (dst.AxisX)dotb, (src.AxisX)dota, (src.AxisX)dotb )    RotC_( (dst.AxisY)dota, (dst.AxisY)dotb, (src.AxisY)dota, (src.AxisY)dotb )    RotC_( (dst.AxisZ)dota, (dst.AxisZ)dotb, (src.AxisZ)dota, (src.AxisZ)dotb )  End Scope#EndMacro#Macro xRot(dst,a_)  RotC(a_,dst,dst,.y,.z)#EndMacro#Macro yRot(dst,a_)  RotC(a_,dst,dst,.z,.x)#EndMacro#Macro zRot(dst,a_)  RotC(a_,dst,dst,.x,.y)#EndMacrotype tMouse  as integer        x,y,b, xp,yp,bp  as single         scalar = 1  declare function  buttons as integer  declare property  dx as single  declare property  dy as singleend typeproperty tmouse.dx as single  return scalar * (x-xp)end propertyproperty tmouse.dy as single  return scalar * (y-yp)end propertyfunction tmouse.buttons as integer  bp = b: xp = x: yp = y  getmouse x,y,,b  return bend functionsub render(buf as imagevars, back_img as imagevars, points() as v3, axis as Axis3D, zoom as single)  screenlock    put (0,0), back_img.im, pset    draw string (0,0), "cone precalculated", rgb(255,255,255)    draw string (0,10), "translation", rgb(0,255,0)    dim as aadot dot    dot.render_target @buf    dot.alpha = 0.75    dim as single z_point_scale = 0.011    For p as v3 ptr = @Points(LBound(points)) to @Points(UBound(points))      dim as single rz_ = (axis.z + axis.AxisZ.z*p->z + axis.AxisX.z*p->x + axis.AxisY.z*p->y)      If rz_ > 0.1 Then        Dim As Single rz = zoom/rz_        Dim as single y = buf.midy - rz*(axis.y + axis.AxisY.y*p->y + axis.AxisZ.y*p->z + axis.AxisX.y*p->x)        dim as single x = buf.midx + rz*(axis.x + axis.AxisX.x*p->x + axis.AxisY.x*p->y + axis.AxisZ.x*p->z)        dot.rad = z_point_scale * rz * p->rad        dot.draw x,y, p->col      end if    Next  screenunlockend subsub translate(byref ret as v3, vec as v3, precalc as v3)      '.axisY.z = .axisX.y  ' Referenced from sub Cone_Precalc()      '.axisZ.x = .axisX.y           '.axisY.x = .axisX.z      '.axisZ.y = .axisX.z         '.axisY.y = .axisX.x      '.axisZ.z = .axisX.x  ret.x = vec.x * precalc.x + vec.y * precalc.z + vec.z * precalc.y  ret.y = vec.y * precalc.x + vec.z * precalc.z + vec.x * precalc.y  ret.z = vec.z * precalc.x + vec.x * precalc.z + vec.y * precalc.y'  ret.x = vec.x * axisX.x + vec.y * axisY.x + vec.z * axisZ.x'  ret.y = vec.y * axisY.y + vec.z * axisZ.y + vec.x * axisX.y'  ret.z = vec.z * axisZ.z + vec.x * axisX.z + vec.y * axisY.zEnd Subsub TranslateAll(a() as v3, surfnorm as v3 = (0, 0, -1), plots as integer = 5000000)  dim as integer ub = ubound(a)\2, num_points = ub+1  for i as integer = 0 to plots - 1    var source = int(rnd * num_points)    var dest = num_points+source    translate a(dest), surfnorm, a(source)    a(dest).rad = 1    a(dest).col = rgb(0,255,0)  nextEnd Subsub Cone_Precalc(a() as v3, num_points as integer=25, cone_half_angle as single = pi / 20)  redim a(num_points-1)  for i as integer = 0 to num_points-1    var mincos = cos(cone_half_angle)    var cosdelta = 1 - mincos    with a(i)         var cosa = mincos + rnd * cosdelta      var sina = sqr(1 - cosa*cosa)      var x_angle = rnd * twopi           .x = cosa      .y = sina * cos(x_angle)      .z = sina * sin(x_angle)           ' 1. imagined to X matrix           '.axisX.x = cosa      '.axisX.y = sina * cos(x_angle)      '.axisX.z = sina * sin(x_angle)           ' 2. perpendicular vectors         '.axisY.z = .axisX.y      '.axisZ.x = .axisX.y           '.axisY.x = .axisX.z      '.axisZ.y = .axisX.z         '.axisY.y = .axisX.x      '.axisZ.z = .axisX.x           .rad = 1      .col = rgb(255,255,255)       end with  nextend subSub Main  dim as imagevars   buf, img  dim as Axis3D   Axis  buf.screen_init 640, 480  img.create      buf.w, buf.h, rgb(48,48,48)  img.checkers rgb(92,92,92), 30  dim as single   anim_fps = 26  dim as single   phys_fps = 55   Axis.z  = 3.0  dim as single   zoom = 0.6 * buf.diagonal  dim as single   rot_plane = pi / 9, rot_plane_inc = 0.0001   dim as v3 points()  cone_precalc points()    redim preserve points( (ubound(points)+1) * 2 - 1)   ' ~~~~~~ anim / physics / input ~~~~~~~  dim as single   anim_f, phys_f, ianim = 1 / anim_fps, iphys = 1 / phys_fps  dim as double   tNow = Timer, td, tp = tNow  dim as double   tDemoExit = tNow+5, tMessageTime = 3  dim as string   kstr  dim as tMouse   mouse    do     if anim_f <= 0 then      translateall points()      render buf, img, points(), axis, zoom      tMessageTime -= ianim      anim_f += ianim    end If    sleep 15    tp = tNow:  tNow = Timer:  td = tNow - tp    anim_f -= td       if mouse.buttons > 0 then      mouse.scalar = Pi / buf.diagonal      yrot(axis, -mouse.dx)      xrot(axis, -mouse.dy)    else      dim as single rspeed = TwoPi / 8500      phys_f += td      while phys_f > 0        zRot(axis, -rot_plane)        yRot(axis, rspeed)        zRot(axis, rot_plane)        rot_plane += rot_plane_inc        phys_f -= iphys      wend    end if       kstr=left(inkey\$,1)     Loop until kstr = chr(27) 'esc End SubMain`
dafhi
Posts: 1258
Joined: Jun 04, 2005 9:51

### Re: Math Guru II How to make a Normal little bit random FAST?

Code: Select all

`' plots visualization with translationtype imagevars ' 2016 Jan 31  as any ptr              im, pixels  as integer              w,h,bpp,bypp,pitch,numpages,flags,rate, is_screen, pitchBy, wm,hm  as single               midx,midy, midxm, midym, diagonal  as string               driver_name  declare sub             scr_inf  Declare Sub             screen_init(wid As single=-1, hgt As single=-1, bpp as UInteger=32, numPages as integer=1, Flags as integer=0)  Declare Function        create(wid As long = -1, hgt As long=0, color As Ulong=RGB(0,0,0)) As Any ptr  declare sub             cls(col as ulong=&HFF000000)  declare sub             checkers(pColor As ULong=RGBA(145,145,145,255), size As UInteger = 12)  Declare Destructor private:  Declare Sub             destroy  declare sub             vars_commonend typeDestructor imagevars:  Destroy:  End DestructorSub imagevars.Destroy():  If im <> 0 Then ImageDestroy im: im = 0: endif:  End SubSub imagevars.cls(col as ulong):  line (0,0)-(wm,hm),col,bf:  end subsub imagevars.vars_common  wm=w-1: midx=w/2: midxm = wm/2  hm=h-1: midy=h/2: midym = hm/2  diagonal=sqr(w*w+h*h): if bypp<>0 then pitchBy=pitch\byppend subsub imagevars.scr_inf  ScreenInfo w,h, bpp, bypp, pitch, rate, driver_name:  pixels = ScreenPtr  vars_commonend subSub imagevars.screen_init(_w As single, _h As single, _bpp as UInteger, _numpages as integer, _flags as integer)  dim as integer ww,hh  ScreenInfo ww,hh  _w=abs(_w): if _w<=1 then _w*=ww: _flags or= 8 'borderless  _h=abs(_h): if _h<=1 then _h*=hh: _flags or= 8  numpages=_numpages: flags=_flags  Destroy: ScreenRes _w,_h,_bpp,numpages,flags: scr_inf: is_screen = -1  if numPages > 1 then screenset 0,1End subFunction imagevars.create(wid As long, hgt As long, col As ULong) As Any Ptr  if hgt=0 then scr_inf: wid=w: hgt=h  Destroy: im = ImageCreate( wid, hgt, col )  ImageInfo im, w, h, bypp, pitch, pixels:  bpp = bypp * 8  vars_common: is_screen = 0: Return imEnd FunctionSub imagevars.checkers(pColor As ULong, size As UInteger)  Dim As UInteger SizeDouble=size*2,SizeM=size-1  For Y as integer = 0 To hm Step size    For X as integer = -size * ((Y/SizeDouble)=Int(Y/SizeDouble)) To wm Step SizeDouble      Line im,(X,Y)-(X+SizeM,Y+SizeM),pColor,BF    Next  NextEnd Sub#Macro Alpha256(ret,back, fore, am, a256)  ret=((_  (fore And &Hff00ff) * a256 + _  (back And &Hff00ff) * am + &H800080) And &Hff00ff00 Or (_  (fore And &H00ff00) * a256 + _  (back And &H00ff00) * am + &H008000) And &H00ff0000) Shr 8#EndMacrotype aadot  as single               rad = 0.65, alpha = 1  as ulong                outlined = 0  as imagevars ptr        p  as any ptr              pixels  declare sub             render_target(pti as imagevars ptr)  declare sub             draw(x as single=0, y as single=0, col as ulong=&HFFFFFFFF) private:  as single               slope, slope_X2end typesub aadot.render_target(pti as imagevars ptr)  if pti->is_screen then: pixels = ScreenPtr  else: pixels = pti->pixels  end if:  p = ptiend subsub aadot.draw(x as single, y as single, col as ulong)  if p->h < 1 then exit sub  slope = alpha*256: slope_X2 = slope*2  dim as single           cone_hgt = rad * slope  dim as integer          _X = x-rad, X2 = x+rad:  _X += _X * (_X < 0)  dim as integer          _Y = y-rad, Y2 = y+rad:  _Y += _Y * (_Y < 0)  dim as single           dxClip = slope * (x - _X), dy = slope * (y - _Y)  X2 += (X2 - p->wm) * (X2 > p->wm)  Y2 += (Y2 - p->hm) * (Y2 > p->hm)  dim as integer          pwidm = X2 - _X  dim as ulong ptr        pCorner = pixels:  pCorner += _Y * p->pitchBy + _X  if outlined then    for py as ulong ptr = pCorner to @pCorner[ (Y2-_Y)*p->pitchBy ] step p->pitchBy      dim as single ySq = dy*dy, dx = dxClip      for px as ulong ptr = py to py + pwidm        dim as integer  alph = cone_hgt - sqr(dx*dx + ySq)        if alph > slope then alph = slope_X2 - alph        if alph > 0 then          dim as integer  AlphM = 256 - alph          Alpha256( *px, *px, col, AlphM, alph )        end if        dx -= slope      next:  dy -= slope    next  else    for py as ulong ptr = pCorner to @pCorner[ (Y2-_Y)*p->pitchBy ] step p->pitchBy      dim as single ySq = dy*dy, dx = dxClip      for px as ulong ptr = py to py + pwidm        dim as integer  alph = cone_hgt - sqr(dx*dx + ySq)        if alph > slope then:  alph = slope        else:  alph += alph * (alph < 0)        end if:  dim as integer  AlphM = 256 - alph        Alpha256( *px, *px, col, AlphM, alph )        dx -= slope      next:  dy -= slope    next  end ifend sub#ifndef PiConst                   TwoPi = 8 * atn(1)Const                   pi    = 4 * Atn(1)#EndIfType v3  As double             x,y,z  as double             rad  as long               colEnd TypeType Axis3D  As v3                 AxisX=(1,0,0), AxisY=(0,-1,0), AxisZ=(0,0,1)  As double             x,y,zEnd Type#Macro RotC_(dsta,dstb,srca,srcb)    temp = cosa_*srca - sina_*srcb    dstb = cosa_*srcb + sina_*srca    dsta = temp#EndMacro#Macro RotC(a_,dst,src,dota,dotb)  Scope    Dim As single a = (a_)    Dim As double cosa_ = Cos(a), sina_ = Sin(a), temp    RotC_( (dst.AxisX)dota, (dst.AxisX)dotb, (src.AxisX)dota, (src.AxisX)dotb )    RotC_( (dst.AxisY)dota, (dst.AxisY)dotb, (src.AxisY)dota, (src.AxisY)dotb )    RotC_( (dst.AxisZ)dota, (dst.AxisZ)dotb, (src.AxisZ)dota, (src.AxisZ)dotb )  End Scope#EndMacro#Macro xRot(dst,a_)  RotC(a_,dst,dst,.y,.z)#EndMacro#Macro yRot(dst,a_)  RotC(a_,dst,dst,.z,.x)#EndMacro#Macro zRot(dst,a_)  RotC(a_,dst,dst,.x,.y)#EndMacrotype tMouse  as integer        x,y,b, xp,yp,bp  as single         scalar = 1  declare function  buttons as integer  declare property  dx as single  declare property  dy as singleend typeproperty tmouse.dx as single  return scalar * (x-xp)end propertyproperty tmouse.dy as single  return scalar * (y-yp)end propertyfunction tmouse.buttons as integer  bp = b: xp = x: yp = y  getmouse x,y,,b  return bend functionsub render(buf as imagevars, back_img as imagevars, points() as v3, axis as Axis3D, zoom as single, LB as integer = -1)  screenlock    put (0,0), back_img.im, pset    draw string (0,0), "original"    draw string (0,10), "translation", rgb(0,255,0)    dim as aadot dot    dot.render_target @buf    dot.alpha = 0.75    dim as single z_point_scale = 0.011    if lb < 0 then lb = lbound(points)    For p as v3 ptr = @Points(LB) to @Points(UBound(points))      dim as single rz_ = (axis.z + axis.AxisZ.z*p->z + axis.AxisX.z*p->x + axis.AxisY.z*p->y)      If rz_ > 0.1 Then        Dim As Single rz = zoom/rz_        Dim as single y = buf.midy - rz*(axis.y + axis.AxisY.y*p->y + axis.AxisZ.y*p->z + axis.AxisX.y*p->x)        dim as single x = buf.midx + rz*(axis.x + axis.AxisX.x*p->x + axis.AxisY.x*p->y + axis.AxisZ.x*p->z)        dot.rad = z_point_scale * rz * p->rad        dot.draw x,y, p->col      end if    Next  screenunlockend subsub translate(byref ret as v3, vec as v3, precalc as v3)  ret.x = vec.x * precalc.x + vec.y * precalc.z + vec.z * precalc.y  ret.y = vec.y * precalc.x + vec.z * precalc.z + vec.x * precalc.y  ret.z = vec.z * precalc.x + vec.x * precalc.z + vec.y * precalc.yEnd Subsub TranslateAll(a() as v3, surfnorm as v3 = (0, 0, -1))  dim as integer ub = ubound(a)\2, num_points = ub+1  for i as integer = 0 to ub    var dest = i+num_points    translate a(dest), surfnorm, a(i)    a(dest).rad = a(i).rad    a(dest).col = rgb(0,255,0)  nextEnd Subsub Cone_Precalc(a() as v3, num_points as integer=25, cone_half_angle as single = pi / 20)  redim a(num_points-1)  for i as integer = 0 to num_points-1    var mincos = cos(cone_half_angle)    var cosdelta = 1 - mincos    with a(i)      var cosa = mincos + rnd * cosdelta      var sina = sqr(1 - cosa*cosa)      var x_angle = rnd * twopi      .x = cosa      .y = sina * cos(x_angle)      .z = sina * sin(x_angle)      .rad = 1      .col = rgb(255,255,255)    end with  nextend subSub Main  dim as imagevars   buf, img  dim as Axis3D   Axis  buf.screen_init 640, 480  img.create      buf.w, buf.h, rgb(48,48,48)  img.checkers rgb(92,92,92), 30  dim as single   anim_fps = 26  dim as single   phys_fps = 55   Axis.z  = 3.0  dim as single   zoom = 0.6 * buf.diagonal  dim as single   rot_plane = pi / 9, rot_plane_inc = 0.0001   dim as v3 points()  cone_precalc points(), 250, pi    redim preserve points( (ubound(points)+1) * 2 - 1)  translateall points()   ' ~~~~~~ anim / physics / input ~~~~~~~  dim as single   anim_f, phys_f, ianim = 1 / anim_fps, iphys = 1 / phys_fps  dim as double   tNow = Timer, td, tp = tNow  dim as double   tDemoExit = tNow+5, tMessageTime = 3  dim as string   kstr  dim as tMouse   mouse    do     if anim_f <= 0 then      var translations_lbound = 0'(ubound(points)+1)\2      render buf, img, points(), axis, zoom, translations_lbound      tMessageTime -= ianim      anim_f += ianim    end If    sleep 15    tp = tNow:  tNow = Timer:  td = tNow - tp    anim_f -= td       if mouse.buttons > 0 then      mouse.scalar = Pi / buf.diagonal      yrot(axis, -mouse.dx)      xrot(axis, -mouse.dy)    else      dim as single rspeed = TwoPi / 8500      phys_f += td      while phys_f > 0        zRot(axis, -rot_plane)        yRot(axis, rspeed)        zRot(axis, rot_plane)        rot_plane += rot_plane_inc        phys_f -= iphys      wend    end if       kstr=left(inkey\$,1)     Loop until kstr = chr(27) 'esc End SubMain`
dodicat
Posts: 5988
Joined: Jan 10, 2006 20:30
Location: Scotland

### Re: Math Guru II How to make a Normal little bit random FAST?

Thanks Dafhi.
I think you have about .3 secomds per 5000000 (from your code 2 up)
(within the TranslateAll)
Is that correct?
If I reduce a sphere to it's eight quadrants, I can speed mine up a bit.
By the way -gen gas is by far the fastest graphics.
Gcc is very slow.

Code: Select all

`randomize ,2screen 19,32dim as any ptr row=screenptrdim as integer pitchdim as ulong ptr pixelscreeninfo ,,,,pitchtype pt    as double x,y,zend type'=============  EIGHT SPHERICAL QUADRANTS ================sub f1(a() as pt)    for n as long=1 to ubound(a)/8      a(n)=type(rnd, _                   rnd*sqr(1-a(n).x*a(n).x), _               sqr((1-(a(n).x*a(n).x+a(n).y*a(n).y))))next nend subsub f2(a() as pt)    for n as long=ubound(a)/8 to 2*ubound(a)/8        a(n)=type(-rnd , _                     rnd*sqr(1-a(n).x*a(n).x), _                 sqr((1-(a(n).x*a(n).x+a(n).y*a(n).y))))next nend subsub f3(a() as pt)    for n as long=2*ubound(a)/8 to 3*ubound(a)/8        a(n)=type(-rnd , _                    -rnd*sqr(1-a(n).x*a(n).x), _                 sqr((1-(a(n).x*a(n).x+a(n).y*a(n).y))))next nend subsub f4(a() as pt)    for n as long=3*ubound(a)/8 to 4*ubound(a)/8        a(n)=type(rnd , _                    -rnd*sqr(1-a(n).x*a(n).x), _                 sqr((1-(a(n).x*a(n).x+a(n).y*a(n).y))))next nend subsub f5(a() as pt)    for n as long=4*ubound(a)/8 to 5*ubound(a)/8        a(n)=type(rnd , _                     rnd*sqr(1-a(n).x*a(n).x), _                 -sqr((1-(a(n).x*a(n).x+a(n).y*a(n).y))))next nend subsub f6(a() as pt)    for n as long=5*ubound(a)/8 to 6*ubound(a)/8        a(n)=type(-rnd , _                     rnd*sqr(1-a(n).x*a(n).x), _                 -sqr((1-(a(n).x*a(n).x+a(n).y*a(n).y))))next nend subsub f7(a() as pt)    for n as long=6*ubound(a)/8 to 7*ubound(a)/8        a(n)=type(-rnd , _                     -rnd*sqr(1-a(n).x*a(n).x), _                 -sqr((1-(a(n).x*a(n).x+a(n).y*a(n).y))))next nend subsub f8(a() as pt)    for n as long=7*ubound(a)/8 to ubound(a)        a(n)=type(rnd , _                     -rnd*sqr(1-a(n).x*a(n).x), _                 -sqr((1-(a(n).x*a(n).x+a(n).y*a(n).y))))next nend subdim f(1 to 8)as sub(() as pt) ={@f1,@f2,@f3,@f4,@f5,@f6,@f7,@f8}'=============================================================================#macro ppset(_x,_y,colour)pixel=row+pitch*(_y)+(_x) shl 2*pixel=(colour)#endmacrodim as long max=5000000redim as pt p(1 to max)'===============================================do     dim as double t=timer 'get the unit vectorsfor n as long=1 to 8  f(n)(p())next n'same as'f1(p()):f2(p()):f3(p()):f4(p()):f5(p()):f6(p()):f7(p()):f8(p())print "Time taken ",timer-t'get error residuedim as double Rfor n as long=1 to max   r+=abs( 1-sqr(p(n).x^2+p(n).y^2+p(n).z^2))     next nprint "ERROR ";rdim as long ctrx,ctry,ctrzfor z as long=lbound(p) to ubound(p)if p(z).x>0 then ctrx+=1if p(z).y>0 then ctry+=1if p(z).z>0 then ctrz+=1next zprint max/2, "Expected for each"print ctrx,ctry,ctrz, "Got"draw string(350,150),"Shade by .z"draw string(350+250,150),"Bare bones"draw string(0,100),"Vector p(10) = " +str(p(10).x)+" , "+str(p(10).y)+" , " +str(p(10).z)screenlockFor z As Long=LBound(p) To UBound(p)   PpSet(clng(400+100*p(z).x),clng(300+100*p(z).y),RGB(Abs(p(z).z)*255,Abs(p(z).z)*255,Abs(p(z).z)*25))   PpSet(clng(400+250+100*p(z).x),clng(300+100*p(z).y),rgb(255,255,255))Nextscreenunlockprint "press a key"sleepclsloop until inkey=chr(27) `
Tourist Trap
Posts: 2768
Joined: Jun 02, 2015 16:24

### Re: Math Guru II How to make a Normal little bit random FAST?

Hi Dahfi, don't take it bad, but I wanted to really look at your nice work, but were enable to read it well without some "clean-up". Same code, just better fit my eyes.

And added some quick wheel zooming.

Code: Select all

`'  subs Cone_Precalc and translate() have all the details'[a bit refactored for comfort (of mine..)]#ifndef _pi    const as double _twoPi => 8*atn(1)    const as double _pi    => 4*atn(1)#endIf#macro _ALPHA256(ret,back, fore, am, a256)  ret=((_          ( fore and &hFF00FF)*a256 + _          ( back and &hFF00FF)*am + &h800080 ) and &hFF00FF00 or _           (_              ( fore and &H00FF00)*a256 + _              ( back and &h00FF00)*am + &h008000 ) and &h00FF0000) shr 8#endMacro#macro _ROTC0(dsta, dstb, srca, srcb)    temp = cosa_*srca - sina_*srcb    dstb = cosa_*srcb + sina_*srca    dsta = temp#endMacro#macro _ROTC(a_,dst, src, dota, dotb)    scope        dim as single   a       = ( a_ )        dim as double   cosa_   = cos(a), sina_ = sin(a), temp        _ROTC0( (dst._axisX)dota, (dst._axisX)dotb, (src._axisX)dota, (src._axisX)dotb )        _ROTC0( (dst._axisY)dota, (dst._axisY)dotb, (src._axisY)dota, (src._axisY)dotb )        _ROTC0( (dst._axisZ)dota, (dst._axisZ)dotb, (src._axisZ)dota, (src._axisZ)dotb )    end scope#endMacro#macro _XROT(dst,a_)    _ROTC(a_, dst, dst, ._y, ._z)#endMacro#macro _YROT(dst,a_)    _ROTC(a_, dst, dst, ._z, ._x)#endMacro#macro _ZROT(dst,a_)    _ROTC(a_, dst, dst, ._x, ._y)#endMacrotype IMAGEVARS ' 2016 Jan 31    declare destructor()    declare sub ClearScreen(byval Col as ulong=&HFF000000)    private:    declare sub Destroy()    declare sub VarsCommon()    public:    declare sub ScrInf()    declare sub ScreenInit(byval Wid as single=-1, _                            byval Hgt as single=-1, _                            byval Bpp as uinteger=32, _                            byval NumPages as integer=1, _                            byval Flags as integer=0)    declare function Create(byval Wid as long=-1, _                             byval Hgt as long=0, _                             byval Colour as ulong=rgb(0,0,0)) _                             as any ptr        declare sub Checkers(byval PColor as ulong=rgba(145,145,145,255), _                          byval Size as uinteger=12)        as any ptr      _im, _pixels        as integer      _w, _h, _bpp, _                         _bypp, _pitch, _numpages, _                         _flags, _rate, _isScreen, _                         _pitchBy, _wm, _hm        as single       _midx, _midy, _                         _midxm, _midym, _diagonal        as string       _driverNameend typedestructor IMAGEVARS()    Destroy()end destructorsub IMAGEVARS.Destroy()    if THIS._im<>0 then        imageDestroy(THIS._im)        THIS._im = 0    end ifend subsub IMAGEVARS.VarsCommon()    THIS._wm    = THIS._w - 1    THIS._midx  = THIS._w/2    THIS._midxm = THIS._wm/2    THIS._hm    = THIS._h - 1    THIS._midy  = THIS._h/2    THIS._midym = THIS._hm/2    THIS._diagonal = sqr(THIS._w*THIS._w + THIS._h*THIS._h)    if THIS._bypp<>0 then THIS._pitchBy = THIS._pitch\THIS._byppend subsub IMAGEVARS.ClearScreen(byval Col as ulong)    line (0,0)-(THIS._wm, THIS._hm), Col, bfend subsub IMAGEVARS.ScrInf()    screenInfo THIS._w, THIS._h, THIS._bpp, _                THIS._bypp, THIS._pitch, THIS._rate, _                THIS._driverName    THIS._pixels = screenPtr()    THIS.VarsCommon()end subsub IMAGEVARS.ScreenInit(byval Wid as single=-1, _                          byval Hgt as single=-1, _                          byval Bpp as uinteger=32, _                          byval NumPages as integer=1, _                          byval Flags as integer=0)    dim as integer ww, hh    screenInfo ww, hh    '    Wid = abs(Wid)    if Wid<=1 then         Wid *= ww    end if    Hgt = abs(Hgt)    if Hgt<=1 then         Hgt *=  hh    end if    THIS._w         = Wid    THIS._h         = Hgt    THIS._Bpp       = Bpp    THIS._flags     or= 8    THIS._numpages  = Numpages    THIS._flags     = Flags     '    THIS.Destroy()    '    screenRes Wid, Hgt, Bpp, Numpages, Flags    '    THIS.ScrInf()    THIS._isScreen = -1    if NumPages> 1 then screenSet 0, 1end subfunction IMAGEVARS.Create(byval Wid as long=-1, _                           byval Hgt as long=0, _                           byval Colour as ulong=rgb(0,0,0)) _                           as any ptr    if Hgt=0 then        ScrInf()        Wid = THIS._w        Hgt = THIS._h    end if    '    THIS.Destroy()    '    THIS._im => imageCreate( Wid, Hgt, Colour )    imageInfo THIS._im, THIS._w, THIS._h, THIS._bypp, THIS._pitch, THIS._pixels    THIS._bpp = THIS._bypp*8    '    THIS.VarsCommon()    THIS._isScreen = 0    '    return THIS._imend functionsub IMAGEVARS.Checkers(byval PColor as ulong, byval Size as uinteger)    dim as uinteger sizeDouble  => Size*2    dim as uinteger sizeM       => Size - 1    for y as integer = 0 to THIS._hm step Size        for x as integer = -Size*((y/sizeDouble)=int(y/sizeDouble)) to _                             THIS._wm step _                             sizeDouble            line THIS._im, (x, y)-(x + sizeM, y + sizeM), PColor, bf        next x    next yEnd Subtype AADOT    declare sub             RenderTarget(byval Pti as IMAGEVARS ptr)    declare sub             DrawAadot(byval X as single=0, _                                       byval Y as single=0, _                                       byval C as ulong=&hFFFFFFFF)        as single               _rad        = 0.65        as single               _alpha      = 1        as ulong                _outlined   = 0        as IMAGEVARS ptr        _p        as any ptr              _pixels    private:        as single               _slope        as single               _slope_X2end typesub AADOT.RenderTarget(byval Pti as IMAGEVARS ptr)    if pti->_isScreen then        THIS._pixels = screenPtr()    else        THIS._pixels = Pti->_pixels    end if    '    THIS._p = Ptiend subsub AADOT.DrawAadot(byval X as single=0, _                     byval Y as single=0, _                     byval C as ulong=&hFFFFFFFF)    if THIS._p->_h<1 then exit sub    '    THIS._slope       = THIS._alpha*256    THIS._slope_X2    = THIS._slope*2    '    dim as single     coneHgt   = THIS._rad*THIS._slope    dim as integer    _x        = X - THIS._rad     dim as integer    x2        = X + THIS._rad          _x += _x*(_x<0)    dim as integer    _y        = Y - THIS._rad    dim as integer    y2        = Y + THIS._rad          _y += _y*(_y<0)    dim as single     dxClip    = THIS._slope*(X - _x)    dim as single     dy        = THIS._slope * (y - _Y)        x2 += (x2 - THIS._p->_wm)*(x2>THIS._p->_wm)        y2 += (y2 - THIS._p->_hm)*(y2>THIS._p->_hm)    dim as integer    pwidm     = x2 - _X    dim as ulong ptr  pCorner   = THIS._pixels        pCorner += _y*THIS._p->_pitchBy + _x    '    if THIS._outlined then        for py as ulong ptr = pCorner to _                               @pCorner[ (y2-_y)*THIS._p->_pitchBy ] step _                               THIS._p->_pitchBy            dim as single ySq   = dy*dy            dim as single dx    = dxClip            for px as ulong ptr = py to py + pwidm                dim as integer  alph = coneHgt - sqr(dx*dx + ySq)                if alph>THIS._slope then alph = THIS._slope_X2 - alph                if alph>0 then                dim as integer  alphM = 256 - alph                _ALPHA256( *px, *px, C, alphM, alph )                end if                '                dx -= THIS._slope            next px            dy -= THIS._slope        next py    else        for py as ulong ptr = pCorner to @pCorner[ (Y2-_Y)*THIS._p->_pitchBy ] step THIS._p->_pitchBy            dim as single ySq = dy*dy, dx = dxClip            for px as ulong ptr = py to py + pwidm                dim as integer  alph = coneHgt - sqr(dx*dx + ySq)                '                if alph> THIS._slope then                    alph = THIS._slope                else                    alph += alph*(alph<0)                end if                '                dim as integer  alphM = 256 - alph                _ALPHA256( *px, *px, C, alphM, alph )                dx -= THIS._slope            next px            dy -= THIS._slope        next py    end ifend subtype V3    as double   _x, _y, _z    as double   _rad    as long     _cend typetype AXIS3D    as V3       _axisX = (1,0,0)    as V3       _axisY = (0,-1,0)    as V3       _axisZ = (0,0,1)    as double   _x, _y, _zend typetype TMOUSE    declare property  Dx() as single    declare property  Dy() as single    declare function  Buttons() as integer        as integer      _x, _y, _b        as integer      _xp, _yp, _bp        as single       _scalar = 1    end typeproperty TMOUSE.dx as single    '    return THIS._scalar*(THIS._x - THIS._xp)end propertyproperty TMOUSE.dy as single    '    return THIS._scalar*(THIS._y - THIS._yp)end propertyfunction TMOUSE.buttons as integer    THIS._bp = THIS._b    THIS._xp = THIS._x    THIS._yp = THIS._y    '    getMouse    THIS._x, THIS._y, , THIS._b    '    return THIS._bend function'. _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ '. + + + + + + + + + + + + + + + + + + + + declare sub Render(byref Buf as IMAGEVARS, _                    byref BackImg as IMAGEVARS, _                    Points() as V3, _                    byref Axis as AXIS3D, _                    byval Zoom as single, _                    byval Lb as integer=-1)sub Render(byref Buf as IMAGEVARS, _            byref BackImg as IMAGEVARS, _            Points() as V3, _            byref Axis as AXIS3D, _            byval Zoom as single, _           byval lB as integer=-1)    screenLock        put (0,0), BackImg._im, PSET        draw string (0,0), "cone precalculated", rgb(255,255,255)        draw string (0,10), "translation", rgb(0,255,0)        '        dim as AADOT    dot                        dot.RenderTarget (@Buf)                        dot._alpha = 0.75        '        dim as single z_point_scale = 0.011        '        if Lb<0 then Lb = lBound(points)        for p as V3 ptr = @Points(Lb) to @Points( uBound(Points) )            dim as single rz_ = ( Axis._z + _                                   Axis._axisZ._z*p->_z + _                                  Axis._axisX._z*p->_x + _                                   Axis._axisY._z*p->_y )            If rz_>0.1 Then                dim as single rz    = Zoom/rz_                dim as single y     = Buf._midy - _                                       rz*(Axis._y + _                                           Axis._axisY._y*p->_y + _                                           Axis._axisZ._y*p->_z + _                                           Axis._AxisX._y*p->_x)                dim as single x     = Buf._midx + _                                       rz*(Axis._x + _                                           Axis._axisX._x*p->_x + _                                           Axis._axisY._x*p->_y + _                                          Axis._axisZ._x*p->_z)                '                dot._rad = z_point_scale*rz*p->_rad                dot.DrawAadot(x, y, p->_c)            end if        Next p    screenUnlockend subdeclare sub Translate(byref Ret as V3, Vec as V3, Precalc as V3)sub Translate(byref Ret as V3, Vec as V3, Precalc as V3)    Ret._x = Vec._x*Precalc._x + Vec._y*Precalc._z + Vec._z*Precalc._y    Ret._y = Vec._y*Precalc._x + Vec._z*Precalc._z + Vec._x*Precalc._y    Ret._z = Vec._z*Precalc._x + Vec._x*Precalc._z + Vec._y*Precalc._yend subdeclare sub TranslateAll(A() as V3, byref SurfNorm as V3=(0, 0, -1))sub TranslateAll(A() as V3, byref SurfNorm as V3=(0, 0, -1))    dim as integer ub           = uBound(A)\2    dim as integer num_points   = ub+1    for i as integer = 0 to ub        var dest    = i+num_points        '        Translate ( A(dest), SurfNorm, A(i) )        A(dest)._rad = A(i)._rad        A(dest)._c = rgb(0,255,0)    next iend subdeclare sub ConePrecalc(A() as V3, _                         byval NumPoints as integer=25, _                         byval ConeHalfAngle as single=_pi/20)sub ConePrecalc(A() as V3, _                 byval NumPoints as integer=25, _                 byval ConeHalfAngle as single=_pi/20)    redim A(NumPoints - 1)    for i as integer = 0 to NumPoints - 1        var mincos      = cos(ConeHalfAngle)        var cosdelta    = 1 - mincos        var cosa        = mincos + rnd()*cosdelta        var sina        = sqr(1 - cosa*cosa)        var x_angle     = rnd()*_twoPi        '        with A(i)            ._x     = cosa            ._y     = sina*cos(x_angle)            ._z     = sina*sin(x_angle)            ._rad   = 1            ._c     = rgb(255,255,255)        end with    next iend subdeclare sub Main()sub Main()    dim as IMAGEVARS      buf, img    dim as AXIS3D         axis    '    buf.ScreenInit(640, 480)    img.Create      buf._w, buf._h, rgb(48,48,48)    img.Checkers    rgb(92,92,92), 30    '    dim as single   anim_fps = 26    dim as single   phys_fps = 55    '    axis._z  = 3.0    dim as single   zoom            = 0.6*buf._diagonal    dim as single   rot_plane       = _pi/9, _                     rot_plane_inc   = .0001    dim as V3       points()    ConePrecalc( points(), 5000, _pi )    '    redim preserve points( (ubound(points) + 1)*2 - 1)    TranslateAll( points() )    '    '~~~~~~ anim / physics / input ~~~~~~    dim as single   anim_f, _                     phys_f, _                     ianim = 1/anim_fps, _                     iphys = 1/phys_fps    '    dim as double   tNow = TIMER, _                     td, _                     tp = tNow    dim as double   tDemoExit = tNow+5, _                     tMessageTime = 3    dim as string   kstr    dim as tMouse   mouse    '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~    do        'addition--------------------------        dim as integer gmX, gmY, gmWheel        getMouse gmX, gmY, gmWheel        zoom = 0.6*buf._diagonal + gmWheel/(40 + log(abs(gmWheel - 100)))*buf._diagonal        '----------------------------------                        if anim_f<=0 then            var translations_lbound = 0 ''(ubound(points)+1)\2            '            Render(buf, img, points(), axis, zoom, translations_lbound)            tMessageTime    -= ianim            anim_f          += ianim        end If        '        tp  =   tNow                tNow = TIMER        td  =   tNow - tp        anim_f -= td        '        if mouse.Buttons()>0 then            mouse._scalar = _pi/buf._diagonal            _YROT(axis, -mouse.Dx)            _XROT(axis, -mouse.Dy)        else            dim as single rspeed    = _twoPi/8500            phys_f += td            while phys_f>0                _ZROT(axis, -rot_plane)                _YROT(axis, rspeed)                _ZROT(axis, rot_plane)                rot_plane += rot_plane_inc                phys_f -= iphys            wend        end if        '        kstr=left(inkey\$,1)        sleep 15    loop until ( kstr=chr(27) )    'end sub'. _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ '. _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ Main()'. _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 'END_OF_SOURCE_FILE`