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

For other topics related to the FreeBASIC project or its community.
dodicat
Posts: 6645
Joined: Jan 10, 2006 20:30
Location: Scotland

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

Tourist Trap.
I'll have a look at your rotate code later.
Here is my final 5 million unit vectors.
(Getting fed up with it now).
64 bit and 32 bit gcc is crap at graphics by direct screen, but fast at the arithmetic.

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 fill(a() as pt)    dim as long X=ubound(a)    #macro subfill(L,U,d1,d2,d3)    for n as long=L to U     a(n)=type(d1 rnd, _                   d2 rnd*sqr(1-a(n).x*a(n).x), _               d3 sqr((1-(a(n).x*a(n).x+a(n).y*a(n).y))))     next n     #endmacro         subfill(1,X/8,,,)     subfill(X/8,2*X/8,-,,)   subfill(2*X/8,3*X/8,-,-,)   subfill(3*X/8,4*X/8,,-,)   subfill(4*X/8,5*X/8,,,-)   subfill(5*X/8,6*X/8,-,,-)   subfill(6*X/8,7*X/8,-,-,-)       subfill(7*X/8,X,,-,-)    end sub'=============================================================================#macro ppset(_x,_y,colour)pixel=row+pitch*(_y)+(_x) shl 2*pixel=(colour)#endmacro#define map(a,b,x,c,d) ((d)-(c))*((x)-(a))\((b)-(a))+(c)dim as long max=5000000redim as pt p(1 to max)'===============================================do     dim as double t=timer 'Fill with random unit vectors fill(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(250,240),"Shade by .z"draw string(350+250,240),"Shade by .x"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(200+150*p(z).x),clng(400+150*p(z).y),RGB(Abs(p(z).z)*255,Abs(p(z).z)*255,Abs(p(z).z)*255))   dim as long c=map(-1,1,p(z).x,10,255)   PpSet(clng(200+350+150*p(z).x),clng(400+150*p(z).y),RGB(c,c,c))Nextscreenunlockprint "press a key"sleepclsloop until inkey=chr(27) `
dafhi
Posts: 1357
Joined: Jun 04, 2005 9:51

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

@dodicat - I see about 0.6 seconds to generate a frame from my example. I think I saw you mention 0.2 with your method on your cpu

@Tourist Trap - I tend to throw stuff together. It makes me happy to know others find my code useful. Add or remove as you see fit!
Tourist Trap
Posts: 2901
Joined: Jun 02, 2015 16:24

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

dafhi wrote: Add or remove as you see fit!

Thanks a lot. I've started playing with it. Was hard to figure out why I had to keep things into untit cube (coordinates<1), but since I got it, I was able to render my Rodrigues rotation delightfully ;-)

@dodicat, I was just looking on a way to render this in 3D perspective. I really don't know how to make it, but the dahfi's renderer is absolutly neat:

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 P3D   as V3type 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 subfunction 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 functiondeclare sub ConePrecalc(A() as V3, _                        byval NumPoints as integer=25, _                        byval ConeHalfAngle as single=_pi/20)sub ConePrecalc(D() as V3, _                byval NumPoints as integer=25, _                byval ConeHalfAngle as single=_pi/20)    redim D(NumPoints - 1)    '   '_____________________________init____   'origin h of the rotation axis (U)   dim as single   hx = 0   dim as single   hy = 0   dim as single   hz = 0   dim as P3D h       h._x = hx       h._y = hy       h._z = hz   'direction unit vector of the axis (U)   dim as single   a   =   0   dim as single   b   =   1   dim as single   c   =   1   dim as single   norm    = sqr(a^2 + b^2 + c^2)   dim as P3D abc       abc._x = a       abc._y = b       abc._z = c   dim as P3D u       u._x = a/norm       u._y = b/norm       u._z = c/norm   'point M to apply rotation on   dim as single   mx   = 0.2   dim as single   my   = 0   dim as single   mz   = 0.6   dim as P3D M       m._x = mx       m._y = my       m._z = mz   'angle of the rotation   dim as double   axisRotationAngle   '____________________________________   'feed    for i as integer = 0 to NumPoints - 1        '        if i<100 then           D(i)._x = hx + i*a/100         D(i)._y = hy + i*b/100         D(i)._z = hz + i*c/100         D(i)._rad = 0.6         D(i)._c = rgb(250,110,050)      else           D(i) = RodriguesRotation(h, u, m, i/10)           with D(i)               ._rad   = 1               ._c     = rgb(255,255,255)           end with       end if    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    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(), 500, _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 = .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`

Just I guess that applying some z_order before showing wouldn't hurt.
dodicat
Posts: 6645
Joined: Jan 10, 2006 20:30
Location: Scotland

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

Tourist Trap.
My perspective is simply +z into the screen.
The larger the z the further away, so any 3D vector shrinks linearly.
if z=0 then any vector remains the same.
if z is -ve then vectors are enlarged.
By a vector I mean a line from one point to another.
If both these points are far away (big +ve z) then the line is shortened .

Code: Select all

`screen 19type pt    as single x,y,z    as ulong c'colour if needed    as long flag 'if needed    declare operator cast() as string '(printout)    #define length(Pt) sqr(Pt.x^2+Pt.y^2+Pt.z^2)    #define distance(Pt1,Pt2)  sqr((Pt1.x-Pt2.x)^2 + (Pt1.y-Pt2.y)^2 + (Pt1.z-Pt2.z)^2)end typeOperator pt.cast() As String#define gap(z) string(25-len(str(z))," ")Return Str(x)+gap(x)+Str(y)+gap(y)+Str(z)End OperatorFunction perspective(p As Pt,eyepoint As Pt) As Pt    Dim As Single   w=1+(p.z/eyepoint.z)    Return Type<Pt>((p.x-eyepoint.x)/w+eyepoint.x,(p.y-eyepoint.y)/w+eyepoint.y,(p.z-eyepoint.z)/w+eyepoint.z,p.c,p.flag)End Functiondim as pt eyepoint=type(350,100,500)'=========  Original point ==============dim as pt p1=(200,200,100)dim as pt p2=(500,400,100)line(P1.x,P1.y)-(P2.x,P2.y),5print p1print p2print distance(p1,p2)'========== New point =====================P1=perspective(p1,eyepoint)P2=perspective(P2,eyepoint)line(P1.x,P1.y)-(P2.x,P2.y)printprintprint p1print p2print distance(p1,p2)circle(350,100),8,15,,,,f 'circle at eyepointsleep `

I wasn't going to bother with any more unit stuff, but I found a little discrepancy.
There wasn't an even spread of vectors around a sphere, there was a kinda corridor of sparse points.
I fixed it, more or less.
Here are the Octants and a sample from each octant:
8 million points, one million in each Octant.

Code: Select all

`Type pt    As Double x,y,z    Declare Operator Cast() As String '(printout)    #define length(Pt) Sqr(Pt.x^2+Pt.y^2+Pt.z^2)End TypeOperator pt.cast() As String#define gap(z) string(25-len(str(z))," ")Return Str(x)+gap(x)+Str(y)+gap(y)+Str(z)End OperatorSub FillWithUnitVectors(a() As pt)    Dim As Long X=Ubound(a),lb=Lbound(a)    #macro subfill(L,U,d1,d2,d3)    For n As Long=L To U        a(n).x= d1 Rnd            a(n).y= d2 Rnd*Sqr(1-(a(n).x*a(n).x))        a(n).z= d3 Sqr((1-(a(n).x*a(n).x+a(n).y*a(n).y)))        If n And 1 Then Swap a(n).x,a(n).z    Next n    #endmacro    subfill(lb,X/8,,,)    subfill(X/8,2*X/8,-,,)    subfill(2*X/8,3*X/8,-,-,)    subfill(3*X/8,4*X/8,,-,)    subfill(4*X/8,5*X/8,,,-)    subfill(5*X/8,6*X/8,-,,-)    subfill(6*X/8,7*X/8,-,-,-)    subfill(7*X/8,X,,-,-)End Sub '++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++Randomize ,2:Redim As pt q(1 To 8000000) Dim As Double t=Timer   FillWithUnitVectors(q()) Print "Time taken",Timer-t;"  seconds  ";"for ";str(ubound(q));"  elements"Print     dim as double residue 'get accumulated error   for n as long=lbound(q) to ubound(q)       residue+=abs(1-sqr(q(n).x^2+q(n).y^2+q(n).z^2))   next n      print "error: ";residuePrint "index",".x",,".y",,".z",,"length"PrintFor n As Ulong=500000 To 8000000 Step 1000000    Print "Q(";n;")=", q(n),,length(q(n))NextPrintprint "done"Sleep `

Code: Select all

` Type pt    As Double x,y,z    Declare Operator Cast() As String '(printout)    #define length(Pt) Sqr(Pt.x^2+Pt.y^2+Pt.z^2)End TypeOperator pt.cast() As String#define gap(z) string(25-len(str(z))," ")Return Str(x)+gap(x)+Str(y)+gap(y)+Str(z)End OperatorSub FillWithUnitVectors(a() As pt)    Dim As Long X=Ubound(a),lb=Lbound(a)    #macro subfill(L,U,d1,d2,d3)    For n As Long=L To U        a(n).x= d1 Rnd            a(n).y= d2 Rnd*Sqr(1-(a(n).x*a(n).x))        a(n).z= d3 Sqr((1-(a(n).x*a(n).x+a(n).y*a(n).y)))        If n And 1 Then Swap a(n).x,a(n).z    Next n    #endmacro    subfill(lb,X/8,,,)    subfill(X/8,2*X/8,-,,)    subfill(2*X/8,3*X/8,-,-,)    subfill(3*X/8,4*X/8,,-,)    subfill(4*X/8,5*X/8,,,-)    subfill(5*X/8,6*X/8,-,,-)    subfill(6*X/8,7*X/8,-,-,-)    subfill(7*X/8,X,,-,-)End Sub '++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++Randomize ,2:Redim As pt q(1 To 8000000) Dim As Double t=Timer   FillWithUnitVectors(q()) Print "Time taken",Timer-t;"  seconds  ";"for ";str(ubound(q));"  elements"Print     dim as double residue 'get accumulated error   for n as long=lbound(q) to ubound(q)       residue+=abs(1-sqr(q(n).x^2+q(n).y^2+q(n).z^2))   next n      print "error: ";residuePrint "index",".x",,".y",,".z",,"length"PrintFor n As Ulong=500000 To 8000000 Step 1000000    Print "Q(";n;")=", q(n),,length(q(n))NextPrintprint "done"Sleep `

Here is a visual look at a sphere.
(This is where I found the sparsity of points)

Code: Select all

`screen 20,32dim as integer xres,yresscreeninfo xres,yresrandomize ,2type pt    as single x,y,z    as ulong c    as long flag 'unused    #define length(Pt) sqr(Pt.x^2+Pt.y^2+Pt.z^2)    #define map(a,b,x,c,d) ((d)-(c))*((x)-(a))\((b)-(a))+(c)end typetype Angle 'FLOATS for angles    as single sx,sy,sz    as single cx,cy,cz    declare static function construct(as single,as single,as single) as Angleend typefunction Angle.construct(x as single,y as single,z as single) as Angle    return    type <Angle>(sin(x),sin(y),sin(z), _                           cos(x),cos(y),cos(z))   end function   Sub QuickSort(array() As pt,begin As Integer,Finish As Uinteger)    Dim As long i=begin,j=finish     Dim As pt x =array(((I+J)\2))    While  I <= J        While array(I).z > X.z: I+=1:wend ' reversable <        While array(J).z < X.z: J-=1:wend ' reversable >If I<=J Then Swap array(I),array(J): I+=1:J-=1    Wend    If J > begin Then QuickSort(array(),begin,J)    If I < Finish Then QuickSort(array(),I,Finish)End Subtype point as ptFunction RotatePoint(c As Point,p As Point,a as Angle,scale As Angle=Type<Angle>(1,1,1)) As Point    Dim As Single dx=p.x-c.x,dy=p.y-c.y,dz=p.z-c.z    Return Type<Point>((scale.sx)*((a.cy*a.cz)*dx+(-a.cx*a.sz+a.sx*a.sy*a.cz)*dy+(a.sx*a.sz+a.cx*a.sy*a.cz)*dz)+c.x,_    (scale.sy)*((a.cy*a.sz)*dx+(a.cx*a.cz+a.sx*a.sy*a.sz)*dy+(-a.sx*a.cz+a.cx*a.sy*a.sz)*dz)+c.y,_    (scale.sz)*((-a.sy)*dx+(a.sx*a.cy)*dy+(a.cx*a.cy)*dz)+c.z,p.c,p.flag)End Function Function perspective(p As Point,eyepoint As Point) As Point    Dim As Single   w=1+(p.z/eyepoint.z)    Return Type<Point>((p.x-eyepoint.x)/w+eyepoint.x,(p.y-eyepoint.y)/w+eyepoint.y,(p.z-eyepoint.z)/w+eyepoint.z,p.c,p.flag)End FunctionFunction Regulate(Byval MyFps As Long,Byref fps As Long=0) As Long            Static As Double timervalue,_lastsleeptime,t3,frames            var t=Timer            frames+=1            If (t-t3)>=1 Then t3=t:fps=frames:frames=0            Var sleeptime=_lastsleeptime+((1/myfps)-T+timervalue)*1000            If sleeptime<1 Then sleeptime=1            _lastsleeptime=sleeptime            timervalue=T            Return sleeptime        End Function            sub FillWithUnitVectors(a() as pt)    dim as long X=ubound(a),lb=lbound(a)    #macro subfill(L,U,d1,d2,d3)    for n as long=L to U               a(n).x= d1 rnd                   a(n).y= d2 rnd*sqr(1-(a(n).x*a(n).x))               a(n).z= d3 sqr((1-(a(n).x*a(n).x+a(n).y*a(n).y)))                if n and 1 then swap a(n).x,a(n).z     next n     #endmacro          subfill(lb,X/8,,,)      subfill(X/8,2*X/8,-,,)   subfill(2*X/8,3*X/8,-,-,)    subfill(3*X/8,4*X/8,,-,)    subfill(4*X/8,5*X/8,,,-)   subfill(5*X/8,6*X/8,-,,-)  subfill(6*X/8,7*X/8,-,-,-)       subfill(7*X/8,X,,-,-)    end sub         redim as pt Q(1 to 30000):randomize ,2       redim as pt rot(lbound(q) to ubound(q))        FillWithUnitVectors(q())    'set colours in q()    for n as long=1 to ubound(q)        q(n).c=rgb(rnd*255,rnd*255,rnd*255)    next        #macro drawpoints(q)    for n as long=1 to ubound(q)        var rad=map(-1,1,q(n).z,3,1)        circle(xres/2+.3*xres*q(n).x,yres/2+.3*xres*q(n).y),rad,q(n).c,,,,f    next    #endmacro        #macro rotate(angle)    for n as long=1 to ubound(q)        rot(n)=rotatepoint(Type<pt>(0,0,0),q(n),angle)        rot(n)=perspective(rot(n),type<pt>(0,0,3))    next    quicksort(rot(),lbound(rot),ubound(rot))    drawpoints(rot)    #endmacro '+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++       dim as single da     dim as string key    dim as long fps    do        key=inkey        if key=" " then FillWithUnitVectors(q())        da+=.005        screenlock        cls        draw string(20,20),"Framerate = " +str(fps) +" --  Press <space> to refresh"        rotate( (Angle.construct(0,da,0)) )        screenunlock        sleep regulate(35,fps),1        loop until key=chr(27)    sleep  `
Tourist Trap
Posts: 2901
Joined: Jun 02, 2015 16:24

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

dodicat wrote:My perspective is simply +z into the screen.

Thanks. So in conclusion, the best is to change the axis as dafhi does. Otherwise the scene is very very hard to control.
dafhi
Posts: 1357
Joined: Jun 04, 2005 9:51

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

@dodicat - my point distribution breakthrough happened when I decided to try representing a sphere with dot rings
dodicat
Posts: 6645
Joined: Jan 10, 2006 20:30
Location: Scotland

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

yea, Dafhi.
I suppose this problem is really how to fill the surface of a sphere randomly.

I have speeded my method up a bit (on -gen GAS anyway).
I populate the whole Northern hemisphere.
The southern Hemisphere is then a copy of the Northern, but with Latitude South.
I.e. .z becomes -.z
It looks o.k. on my sphere code.
Here is the timer code, starting with 5 million, and continuing in excess of 5 million.

Code: Select all

`Type pt    As Double x,y,z    Declare Operator Cast() As String '(printout)    #define length(Pt) Sqr(Pt.x^2+Pt.y^2+Pt.z^2)End TypeOperator pt.cast() As String#define gap(z) string(25-len(str(z))," ")Return Str(x)+gap(x)+Str(y)+gap(y)+Str(z)End OperatorSub FillWithUnitVectors(a() As pt)    Dim As Long X=Ubound(a),lb=Lbound(a),half=X\2-1    #macro subfill(L,U,d1,d2,d3)    For n As Long=L To U        a(n).x= d1 Rnd            a(n).y= d2 Rnd*Sqr(1-(a(n).x*a(n).x))        a(n).z= d3 Sqr((1-(a(n).x*a(n).x+a(n).y*a(n).y)))    Next n    #endmacro    subfill(lb,X\8,,,)'First hemisphere, 4 quadrants.    subfill(X\8,2*X\8,-,,)    subfill(2*X\8,3*X\8,-,-,)    subfill(3*X\8,4*X\8,,-,)    For n As Long=X\2 To X 'get the other hemisphere        a(n).x=a(n-half).x        a(n).y=a(n-half).y        a(n).z=-a(n-half).z        If n And 1 Then Swap a(n).x,a(n).z:Swap a(n-half).x,a(n-half).z 'plug any gaps    Next nEnd Sub'++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#define Intrange(f,l) int(Rnd*((l+1)-(f))+(f))dim as long p=5000000doRandomize ,2:Redim  As pt q(1 To p)  Dim As Double t=TimerFillWithUnitVectors(q())Print "Time taken",Timer-t;"  seconds  ";"for ";Str(Ubound(q));"  elements"Print  Dim As Double residue 'get accumulated errorFor n As Long=Lbound(q) To Ubound(q)    residue+=Abs(1-Sqr(q(n).x^2+q(n).y^2+q(n).z^2))Next nPrint "error: ";residuePrint "index",".x",,".y",,".z",,"length"PrintFor n As Ulong=.5*p\8 To p Step p\8    Print "Q(";n;")=", q(n),,length(q(n))    NextPrintPrint "Press a key" p=5000000+ IntRange(0,3000000)sleepclsloop until inkey=chr(27)Sleep `
dodicat
Posts: 6645
Joined: Jan 10, 2006 20:30
Location: Scotland

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

Last one.
Just one octant filled randomly, the rest of the sphere surface is filled with reflections.
Timer code

Code: Select all

`Type pt    As Double x,y,z    Declare Operator Cast() As String '(printout)    #define length(Pt) Sqr(Pt.x^2+Pt.y^2+Pt.z^2)End TypeOperator pt.cast() As String#define gap(z) string(25-len(str(z))," ")Return Str(x)+gap(x)+Str(y)+gap(y)+Str(z)End OperatorSub FillWithUnitVectors(a() As pt)    Dim As Long X=Ubound(a),lb=Lbound(a),half=X\2-1,octt=X\8-1,qtr=X\4-1,trq=3*X\4-1    #macro subfill(L,U,d1,d2,d3)    For n As Long=L To U        a(n).x= d1 Rnd            a(n).y= d2 Rnd*Sqr(1-(a(n).x*a(n).x))        a(n).z= d3 Sqr((1-(a(n).x*a(n).x+a(n).y*a(n).y)))        If n And 1 Then Swap a(n).y,a(n).z    Next n    #endmacro    #macro cpy(L,U,d1,d2,d3,r)    For n As Long=L To U        a(n).x= d1 a(n-r).x           a(n).y= d2 a(n-r).y           a(n).z= d3 a(n-r).z    Next n    #endmacro    subfill(lb,X\8,,,) 'random    cpy(X\8,2*X\8,-,,,octt)   ' copies of octant 1    cpy(2*X\8,3*X\8,-,-,,octt)    cpy(3*X\8,4*X\8,,-,,qtr)        cpy(4*X\8,5*X\8,,,-,half)    cpy(5*X\8,6*X\8,-,,-,qtr)    cpy(6*X\8,7*X\8,-,-,-,half)    cpy(7*X\8,X,,-,-,trq)End Sub'++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#define Intrange(f,l) int(Rnd*((l+1)-(f))+(f))Dim As Long p=5000000Do    Randomize ,2:Redim  As pt q(1 To p)             Dim As Double t=Timer    FillWithUnitVectors(q())    Print "Time taken",Timer-t;"  seconds  ";"for ";Str(Ubound(q));"  elements"    Print      Dim As Double residue 'get accumulated error    For n As Long=Lbound(q) To Ubound(q)        residue+=Abs(1-Sqr(q(n).x^2+q(n).y^2+q(n).z^2))    Next n        Print "error: ";residue        Print "index",".x",,".y",,".z",,"length"    Print    For n As Ulong=.5*p\8 To p Step p\8        Print "Q(";n;")=", q(n),,length(q(n))            Next    Print    Print "Press a key"    p=5000000+ IntRange(0,3000000)    Sleep    ClsLoop Until Inkey=Chr(27)Sleep `

About one tenth of second for five million points on this box.
That's me done here now.
D.J.Peters
Posts: 8125
Joined: May 28, 2005 3:28
Contact:

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

no comment :-)

Joshy

Code: Select all

`type tReal as singletype v3_1  as tReal x,y,zend type#define v3Add(l,r) type<v3_1>(l.x+r.x, l.y+r.y, l.z+r.z)type v3_2  as tReal x,y,zend typeoperator +(l as v3_2,r as v3_2) as v3_2  operator = type<v3_2>(l.x+r.x, l.y+r.y, l.z+r.z)end operatortype v3_3  declare constructor(x as tReal=0,y as tReal=0,z as tReal=0)  as tReal x,y,zend typeconstructor v3_3(x as tReal=0,y as tReal=0,z as tReal=0)  this.x=x : this.y=y : this.z=zend constructoroperator +(l as v3_3,r as v3_3) as v3_3  operator = v3_3(l.x+r.x, l.y+r.y, l.z+r.z)end operatorconst as integer count = 100000000dim as v3_1 a1,b1,c1dim as v3_2 a2,b2,c2dim as v3_3 a3,b3,c3print "wait ..."var tStart = Timer()for i as integer= 1 to count  c1=v3Add(a1,b1)nextvar tResult = Timer()-tStartprint "v3_1 " & tResulttStart = Timer()for i as integer= 1 to count  c2=a2+b2nexttResult = Timer()-tStartprint "v3_2 " & tResulttStart = Timer()for i as integer= 1 to count  c3=a3+b3nexttResult = Timer()-tStartprint "v3_3 " & tResultprint "ok ..."sleep`
fxm
Posts: 9825
Joined: Apr 22, 2009 12:46
Location: Paris suburbs, FRANCE

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

return v3_3(l.x+r.x, l.y+r.y, l.z+r.z)
is faster than
operator = v3_3(l.x+r.x, l.y+r.y, l.z+r.z)
because the first syntax calls once the constructor and once operator let, instead of twice each for the second syntax.
dafhi
Posts: 1357
Joined: Jun 04, 2005 9:51

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

The examples I gave in this thread are incorrect.

I really want to write a path tracer so I will revisit dodicat's and gothon's work
dafhi
Posts: 1357
Joined: Jun 04, 2005 9:51

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

so easy now

Code: Select all

`...`
Last edited by dafhi on Dec 02, 2016 5:34, edited 3 times in total.
Tourist Trap
Posts: 2901
Joined: Jun 02, 2015 16:24

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

fxm wrote:return v3_3(l.x+r.x, l.y+r.y, l.z+r.z)
is faster than
operator = v3_3(l.x+r.x, l.y+r.y, l.z+r.z)
because the first syntax calls once the constructor and once operator let, instead of twice each for the second syntax.

Do you mean that the = sign makes the things slower (double work)?
fxm
Posts: 9825
Joined: Apr 22, 2009 12:46
Location: Paris suburbs, FRANCE

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

Tourist Trap wrote:Do you mean that the = sign makes the things slower (double work)?
Yes, but that occurs only in a particular case.

- When an object is returned (by value) from a function (or an operator), by using 'Function = x' (or 'Operator = x') assignment, the default constructor is called once at first, and then the let operator at each assignment.
- When an object is returned (by value) from a function (or an operator), by using 'Return x' statement, the copy constructor is called once at end.

But compiler optimizes 'Return x' when 'x' is already the direct call to any constructor of the UDT: in that case, the standard use of copy constructor is suppressed and only the 'x' constructor call is kept.

Test code for enthusiasts of the subject:

Code: Select all

`Type UDT  Declare Constructor()  Declare Constructor(Byval I0 As Integer)  Declare Constructor(Byref u As UDT)  Declare Operator Let(Byref u As UDT)  Declare Destructor()  Declare Function ff() As UDT  Declare Function fr() As UDT  Declare Function ffc() As UDT  Declare Function frc() As UDT  Declare Function ffci() As UDT  Declare Function frci() As UDT  Declare Function ffcc() As UDT  Declare Function frcc() As UDT  Dim As Integer IEnd TypeConstructor UDT()  Print "constructor UDT()",,, @ThisEnd ConstructorConstructor UDT(Byval I0 As Integer)  Print "constructor UDT(As Integer)",,, @This  This.I = I0End ConstructorConstructor UDT(Byref u As UDT)  Print "constructor UDT(As UDT)",,, @This, @u  This.I = u.IEnd ConstructorOperator UDT.Let(Byref u As UDT)  Print "operator UDT.Let(As UDT)",,, @This, @u  This.I = u.IEnd OperatorDestructor UDT()  Print "destructor UDT()",,, @ThisEnd DestructorFunction UDT.ff() As UDT  Function = ThisEnd FunctionFunction UDT.fr() As UDT  Return ThisEnd FunctionFunction UDT.ffc() As UDT  Function = UDT()End FunctionFunction UDT.frc() As UDT  Return UDT()End FunctionFunction UDT.ffci() As UDT  Function = UDT(123)End FunctionFunction UDT.frci() As UDT  Return UDT(123)End FunctionFunction UDT.ffcc() As UDT  Function = UDT(This)End FunctionFunction UDT.frcc() As UDT  Return UDT(This)End FunctionWidth , 47Scope  Dim u As UDT = UDT(123)  Print "------------------------ 'Function = This' ----------------------------------"  Print u.ff().I  Print "------------------------ 'Return This' --------------------------------------"  Print u.fr().I  Print "------------------------ 'Function = UDT()' ---------------------------------"  Print u.ffc().I  Print "------------------------ 'Return UDT()' -------------------------------------"  Print u.frc().I  Print "------------------------ 'Function = UDT(123)' ------------------------------"  Print u.ffci().I  Print "------------------------ 'Return UDT(123)' ----------------------------------"  Print u.frci().I  Print "------------------------ 'Function = UDT(This)' -----------------------------"  Print u.ffcc().I  Print "------------------------ 'Return UDT(This)' ---------------------------------"  Print u.frcc().I  Print "-----------------------------------------------------------------------------"End ScopeSleep`

Code: Select all

`constructor UDT(As Integer)                             1638072------------------------ 'Function = This' ----------------------------------constructor UDT()                                       1638068operator UDT.Let(As UDT)                                1638068       1638072 123destructor UDT()                                        1638068------------------------ 'Return This' --------------------------------------constructor UDT(As UDT)                                 1638064       1638072 123destructor UDT()                                        1638064------------------------ 'Function = UDT()' ---------------------------------constructor UDT()                                       1638060constructor UDT()                                       1638000operator UDT.Let(As UDT)                                1638060       1638000destructor UDT()                                        1638000 0destructor UDT()                                        1638060------------------------ 'Return UDT()' -------------------------------------constructor UDT()                                       1638056 0destructor UDT()                                        1638056------------------------ 'Function = UDT(123)' ------------------------------constructor UDT()                                       1638052constructor UDT(As Integer)                             1638000operator UDT.Let(As UDT)                                1638052       1638000destructor UDT()                                        1638000 123destructor UDT()                                        1638052------------------------ 'Return UDT(123)' ----------------------------------constructor UDT(As Integer)                             1638048 123destructor UDT()                                        1638048------------------------ 'Function = UDT(This)' -----------------------------constructor UDT()                                       1638044constructor UDT(As UDT)                                 1638000       1638072operator UDT.Let(As UDT)                                1638044       1638000destructor UDT()                                        1638000 123destructor UDT()                                        1638044------------------------ 'Return UDT(This)' ---------------------------------constructor UDT(As UDT)                                 1638040       1638072 123destructor UDT()                                        1638040-----------------------------------------------------------------------------destructor UDT()                                        1638072`

Same kind of compiler optimization (only one constructor called) than for equivalent expressions to 'Return x':
Dim As UDT u0 = UDT()
or
Dim As UDT u1 = UDT(123)
or
Dim As UDT u2 = UDT(u)

This above is more performant than that follows (one constructor and one let operator and one destructor called in addition) equivalent to 'Function = x' ('Operator = x'):
Dim As UDT u0
u0 = UDT()

or
Dim As UDT u1
u1 = UDT(123)

or
Dim As UDT u2
u2 = UDT(u)

This kind of compiler optimization in not a recent improvement but exists since the origin when constructor, destructor, method, properties and operator overloading have been added in the version fbc 0.17b.
Last edited by fxm on Dec 02, 2016 15:33, edited 7 times in total.
Stonemonkey
Posts: 646
Joined: Jun 09, 2005 0:08

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

Hi Joshy, you don't say what you're doing with the normals after, depending on what you're doing there can be ways of reducing the number of reciprocal square roots required.

Code: Select all

`type v3d    as single x,y,zend type#define dot(v0,v1) (v0.x*v1.x+v0.y*v1.y+v0.z*v1.z)#define random3d(v,s) v.x=rnd*s-.5*s:v.y=rnd*s-.5*s:v.z=rnd*s-.5*s#define normalise(v) scope:dim as single d=1.0/sqr(dot(v,v)):v.x*=d:v.y*=d:v.z*=d:end scoperandomize timerdim as v3d v,nrandom3d(v,10)random3d(n,10)print dot(v,n)/sqr(dot(v,v)*dot(n,n))normalise(v)normalise(n)print dot(v,n)sleepend`