/****************************************************************************** Floating Point Performance, Tested on Intel i7-3632QM 2.2 Ghz Operations Flt x Flt, have the same performance as Dbl x Dbl. However Flt x Dbl (mixed precision) are slower. Casting to the same types restores full performance: Flt x Flt(Dbl), without any cost for casting. Floating Point Precision: After 1 month of running an app: Flt time=60*60*24*30; (60s*60m*24h*30d = 2592000 seconds) At this value, 'time' + average frame time delta (1.0f/60) is the same as 'time' time+=1.0f/60; // does not modify 'time' /******************************************************************************/ #include "stdafx.h" namespace EE{ /******************************************************************************/ const Half HalfZero=false, HalfOne=true; const Int PrimeNumbers[16]={2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53}; /******************************************************************************/ Bool Special(C Flt &r) {UInt *d=(UInt*)&r; return (d[0]&0x7F800000)==0x7F800000 ;} // all exponent bits are on Bool Special(C Dbl &r) {UInt *d=(UInt*)&r; return (d[1]&0x7FF00000)==0x7FF00000 ;} // all exponent bits are on Bool NaN (C Flt &r) {UInt *d=(UInt*)&r; return Special(r) && (d[0]&0x7FFFFF );} // significand!=0 Bool Inf (C Flt &r) {UInt *d=(UInt*)&r; return Special(r) && !(d[0]&0x7FFFFF );} // significand==0 Bool NaN (C Dbl &r) {UInt *d=(UInt*)&r; return Special(r) && (d[1]&0x0FFFFF || d[0]);} // significand!=0 Bool Inf (C Dbl &r) {UInt *d=(UInt*)&r; return Special(r) && !(d[1]&0x0FFFFF || d[0]);} // significand==0 /******************************************************************************/ Vec2 ScaleFactor(C Vec2 &vec) {return Vec2 (ScaleFactor(vec.x), ScaleFactor(vec.y));} VecD2 ScaleFactor(C VecD2 &vec) {return VecD2(ScaleFactor(vec.x), ScaleFactor(vec.y));} Vec ScaleFactor(C Vec &vec) {return Vec (ScaleFactor(vec.x), ScaleFactor(vec.y), ScaleFactor(vec.z));} VecD ScaleFactor(C VecD &vec) {return VecD (ScaleFactor(vec.x), ScaleFactor(vec.y), ScaleFactor(vec.z));} Vec4 ScaleFactor(C Vec4 &vec) {return Vec4 (ScaleFactor(vec.x), ScaleFactor(vec.y), ScaleFactor(vec.z), ScaleFactor(vec.w));} VecD4 ScaleFactor(C VecD4 &vec) {return VecD4(ScaleFactor(vec.x), ScaleFactor(vec.y), ScaleFactor(vec.z), ScaleFactor(vec.w));} Vec2 ScaleFactorR(C Vec2 &vec) {return Vec2 (ScaleFactorR(vec.x), ScaleFactorR(vec.y));} VecD2 ScaleFactorR(C VecD2 &vec) {return VecD2(ScaleFactorR(vec.x), ScaleFactorR(vec.y));} Vec ScaleFactorR(C Vec &vec) {return Vec (ScaleFactorR(vec.x), ScaleFactorR(vec.y), ScaleFactorR(vec.z));} VecD ScaleFactorR(C VecD &vec) {return VecD (ScaleFactorR(vec.x), ScaleFactorR(vec.y), ScaleFactorR(vec.z));} Vec4 ScaleFactorR(C Vec4 &vec) {return Vec4 (ScaleFactorR(vec.x), ScaleFactorR(vec.y), ScaleFactorR(vec.z), ScaleFactorR(vec.w));} VecD4 ScaleFactorR(C VecD4 &vec) {return VecD4(ScaleFactorR(vec.x), ScaleFactorR(vec.y), ScaleFactorR(vec.z), ScaleFactorR(vec.w));} /******************************************************************************/ void MinMax(C Flt *f, Int elms, Flt &min, Flt &max) { if(elms)for(min=max=*f++; --elms; ) { Flt r=*f++; if(rmax)max=r; } } void MinMax(C Dbl *f, Int elms, Dbl &min, Dbl &max) { if(elms)for(min=max=*f++; --elms; ) { Dbl r=*f++; if(rmax)max=r; } } void MinMaxI(C Flt *f, Int elms, Int &min, Int &max) { if(elms) { Flt _min, _max; _min=_max=f[0]; min= max= 0 ; for(Int i=1; i_max){_max=r; max=i;} } } } void MinMaxI(C Dbl *f, Int elms, Int &min, Int &max) { if(elms) { Dbl _min, _max; _min=_max=f[0]; min= max= 0 ; for(Int i=1; i_max){_max=r; max=i;} } } } /******************************************************************************/ UInt SqrtI(UInt x, Int max_steps) { if(x<=1)return x; UInt y=1<<(BitHi(x)>>1); // approximate the sqrt REP(max_steps) { UInt prev=y; y=(y + x/y + 1)/2; if(prev==y)break; // if no change was made then stop } return y; } UInt SqrtI(UInt x) { UInt res=0, one=1<<30; for(; one>x; one>>=2); for(; one; ) { if(x>=res+one) { x -=res+one; res+=one<<1; } res>>=1; one>>=2; } if(x>res)res++; // use rounding return res; } UInt SqrtI(ULong x) { ULong res=0, one=1ll<<62; for(; one>x; one>>=2); for(; one; ) { if(x>=res+one) { x -=res+one; res+=one<<1; } res>>=1; one>>=2; } if(x>res)res++; // use rounding return res; } Int SqrtI(Int x) {return (x>0) ? SqrtI(UInt (x)) : 0;} Int SqrtI(Long x) {return (x>0) ? SqrtI(ULong(x)) : 0;} Flt Log(Flt x, Flt base) {return logf(x)/logf(base);} Dbl Log(Dbl x, Dbl base) {return log (x)/log (base);} Flt Pinch(Flt x, Flt pinch) {return x*pinch/(1+x*(pinch-1));} /******************************************************************************/ // ROUNDING /******************************************************************************/ Dbl TruncD(Dbl x) { Dbl trunc; modf(x, &trunc); return trunc; } /******************************************************************************/ Flt SplitAlpha(Flt alpha, Int steps) { /* For 2 steps: apply regular alpha blending: alpha*src + (1-alpha)*dest = final apply regular alpha blending but two times in a row (now using unknown 'beta' instead of 'alpha'): beta*src + (1-beta)*(beta*src + (1-beta)*dest) = final now we want 'final' result to be the same in both cases: alpha*src + (1-alpha)*dest = beta*src + (1-beta)*(beta*src + (1-beta)*dest) alpha*src + dest - dest*alpha = beta*src + (1-beta)*(beta*src + dest - dest*beta) alpha*src + dest - dest*alpha = beta*src + beta*src + dest - dest*beta - beta*beta*src - dest*beta + dest*beta*beta alpha*src - dest*alpha = beta*src + beta*src - dest*beta - beta*beta*src - dest*beta + dest*beta*beta alpha*src - dest*alpha = beta*beta*(dest - src) + beta*2*(src - dest) alpha*src - dest*alpha = beta*beta*(dest - src) - beta*2*(dest - src) (alpha*src - dest*alpha)/(dest-src) = beta*beta - beta*2 -alpha = beta*beta - beta*2 -> assumed dest=1, src=0 beta*beta + (-2)*beta + alpha = 0 delta = -2*-2 - 4*alpha*1 return 1-Sqrt(4-4*alpha)*0.5 return 1-2*Sqrt(1-alpha)*0.5; gives final result of: return 1-Sqrt(1-alpha); it was tested that for 3 steps: return 1-Cbrt(1-alpha); gives positive results So the final formula is: */ return (steps>1) ? 1-Pow(1-alpha, 1.0f/steps) : alpha; } Flt VisibleOpacity(Flt density, Flt range) {return Pow(1-density, range);} Flt AccumulatedDensity(Flt density, Flt range) {return 1-Pow(1-density, range);} /******************************************************************************/ // ANGLES /******************************************************************************/ Vec2 Tan(C Vec2 &angle) {return Vec2(Tan(angle.x), Tan(angle.y));} void CosSin(Flt &cos, Flt &sin, Flt angle) { #if APPLE __sincosf(angle, &sin, &cos); #elif !WINDOWS sincosf(angle, &sin, &cos); #elif X64 || ARM sin=Sin(angle); cos=Cos(angle); #else _asm { mov edx, cos mov eax, sin fld angle fsincos fstp dword ptr[edx] fstp dword ptr[eax] } #endif } void CosSin(Dbl &cos, Dbl &sin, Dbl angle) { #if APPLE __sincos(angle, &sin, &cos); #elif !WINDOWS sincos(angle, &sin, &cos); #elif X64 || ARM sin=Sin(angle); cos=Cos(angle); #else _asm { mov edx, cos mov eax, sin fld angle fsincos fstp qword ptr[edx] fstp qword ptr[eax] } #endif } Flt Acos(Flt cos) { if(cos>= 1)return 0; if(cos<=-1)return PI; return acosf(cos); } Dbl Acos(Dbl cos) { if(cos>= 1)return 0; if(cos<=-1)return PID; return acos(cos); } Flt Asin(Flt sin) { if(sin>= 1)return PI_2; if(sin<=-1)return -PI_2; return asinf(sin); } Dbl Asin(Dbl sin) { if(sin>= 1)return PID_2; if(sin<=-1)return -PID_2; return asin(sin); } Flt ACosSin(Flt cos, Flt sin) // assumes "sin>=0" { // use Acos for angles 45..135 deg return (sin>=SQRT2_2) ? Acos(cos) : (cos>=0) ? Asin(sin) : PI-Asin(sin); } Dbl ACosSin(Dbl cos, Dbl sin) // assumes "sin>=0" { // use Acos for angles 45..135 deg return (sin>=SQRT2_2) ? Acos(cos) : (cos>=0) ? Asin(sin) : PID-Asin(sin); } Vec2 Atan(C Vec2 &tan) {return Vec2(Atan(tan.x), Atan(tan.y));} Flt AngleFast(Flt x, Flt y) { Flt r=Abs(x)+Abs(y); if(!r)return 0; return (y>=0) ? (1-x/r)*PI_2 : (x/r-1)*PI_2; } Flt AngleFast(C Vec2 &v) { Flt r=Abs(v.x)+Abs(v.y); if(!r)return 0; return (v.y>=0) ? (1-v.x/r)*PI_2 : (v.x/r-1)*PI_2; } Flt AngleFast(C Vec &v, C Vec &x, C Vec &y) {return AngleFast(Dot(v, x), Dot(v, y));} Flt Angle (C Vec &v, C Vec &x, C Vec &y) {return Angle (Dot(v, x), Dot(v, y));} /******************************************************************************/ Flt AngleNormalize(Flt angle) {angle=AngleFull(angle); return (angle>PI ) ? angle-PI2 : angle;} Dbl AngleNormalize(Dbl angle) {angle=AngleFull(angle); return (angle>PID) ? angle-PID2 : angle;} /******************************************************************************/ Flt CosBetween(C Vec2 &a, C Vec2 &b) {if(Flt l2=a.length2()*b.length2())return Dot(a, b)/SqrtFast(l2); return 0;} Dbl CosBetween(C VecD2 &a, C VecD2 &b) {if(Dbl l2=a.length2()*b.length2())return Dot(a, b)/SqrtFast(l2); return 0;} Flt CosBetween(C Vec &a, C Vec &b) {if(Flt l2=a.length2()*b.length2())return Dot(a, b)/SqrtFast(l2); return 0;} Dbl CosBetween(C VecD &a, C VecD &b) {if(Dbl l2=a.length2()*b.length2())return Dot(a, b)/SqrtFast(l2); return 0;} Flt SinBetween(C Vec2 &a, C Vec2 &b) {if(Flt l2=a.length2()*b.length2())return Cross(a, b)/ SqrtFast(l2); return 0;} Dbl SinBetween(C VecD2 &a, C VecD2 &b) {if(Dbl l2=a.length2()*b.length2())return Cross(a, b)/ SqrtFast(l2); return 0;} Flt AbsSinBetween(C Vec &a, C Vec &b) {if(Flt l2=a.length2()*b.length2())return SqrtFast(Cross(a, b).length2()/l2); return 0;} Dbl AbsSinBetween(C VecD &a, C VecD &b) {if(Dbl l2=a.length2()*b.length2())return SqrtFast(Cross(a, b).length2()/l2); return 0;} Flt AngleBetween (C Vec2 &a, C Vec2 &b) {return AngleDelta(Angle(a), Angle(b));} Dbl AngleBetween (C VecD2 &a, C VecD2 &b) {return AngleDelta(Angle(a), Angle(b));} Flt AbsAngleBetween (C Vec &a, C Vec &b) {return Angle (CosBetweenN(a, b), AbsSinBetweenN(a, b));} // we don't need to do any scaling because 'CosBetweenN' and 'AbsSinBetweenN' will have proportional lengths, and we use 'Angle' which doesn't require lengthts to be normalized Dbl AbsAngleBetween (C VecD &a, C VecD &b) {return Angle (CosBetweenN(a, b), AbsSinBetweenN(a, b));} // we don't need to do any scaling because 'CosBetweenN' and 'AbsSinBetweenN' will have proportional lengths, and we use 'Angle' which doesn't require lengthts to be normalized Flt AbsAngleBetweenN(C Vec &a, C Vec &b) {return ACosSin(CosBetweenN(a, b), AbsSinBetweenN(a, b));} Dbl AbsAngleBetweenN(C VecD &a, C VecD &b) {return ACosSin(CosBetweenN(a, b), AbsSinBetweenN(a, b));} Flt AngleBetween(C Vec &a, C Vec &b, C Vec &z) {Flt angle=AbsAngleBetween(a, b); if(Dot(Cross(a, b), z)<0)CHS(angle); return angle;} Dbl AngleBetween(C VecD &a, C VecD &b, C VecD &z) {Dbl angle=AbsAngleBetween(a, b); if(Dot(Cross(a, b), z)<0)CHS(angle); return angle;} /******************************************************************************/ Vec DequantizeNormal(C Vec &n) { Flt z=CalcZ(n.xy); Vec dn(n.x, n.y, Lerp((n.z>=0) ? z : -z, n.z, n.xy.length())); // error 5232, errorNoNormalize 1583174 //Vec dn(n.x, n.y, Lerp((n.z>=0) ? z : -z, n.z, CosSin(z))); // error 5232, errorNoNormalize 1583103 //Vec dn(n.x, n.y, Lerp(n.z, (n.z>=0) ? z : -z, Sqr(z))); // error 5425, errorNoNormalize 1410261 //Vec dn(n.x, n.y, Lerp((n.z>=0) ? z : -z, n.z, 1-Sqr(z))); // error 5425, errorNoNormalize 1410261 //Vec dn(n.x, n.y, Lerp((n.z>=0) ? z : -z, n.z, n.xy.length2())); // error 5426, errorNoNormalize 1410402 //Vec dn(n.x, n.y, Lerp(n.z, (n.z>=0) ? z : -z, z)); // error 7035, errorNoNormalize 1101598 //Vec dn(n.x, n.y, Lerp((n.z>=0) ? z : -z, n.z, 1-z)); // error 7035, errorNoNormalize 1101598 //Vec dn(n.x, n.y, Lerp((n.z>=0) ? z : -z, n.z, Sqr(1-z))); // error 10519, errorNoNormalize 793122 //Vec dn(n.x, n.y, (n.z>=0) ? z : -z); // error 120477, errorNoNormalize 144597 dn.normalize(); return dn; } /******************************************************************************/ // INTERPOLATION /******************************************************************************/ #define TANGENT_FULL 0.817538f #define TANGENT 0.5f /******************************************************************************/ Flt GetTangent ( Flt prev, Flt next) {return (next-prev)*TANGENT;} Vec2 GetTangent (C Vec2 &prev, C Vec2 &next) {return (next-prev)*TANGENT;} Vec GetTangent (C Vec &prev, C Vec &next) {return (next-prev)*TANGENT;} Vec4 GetTangent (C Vec4 &prev, C Vec4 &next) {return (next-prev)*TANGENT;} Vec2 GetTangentDir(C Vec2 &prev, C Vec2 &next) {return (next-prev)*Lerp(TANGENT, TANGENT_FULL, Sqr(AbsAngleBetween(prev, next)/PI));} Vec GetTangentDir(C Vec &prev, C Vec &next) {return (next-prev)*Lerp(TANGENT, TANGENT_FULL, Sqr(AbsAngleBetween(prev, next)/PI));} Flt GetTangent( Flt prev, Flt cur, Flt next) {return GetTangent (prev , next );} Vec2 GetTangent(C Vec2 &prev, C Vec2 &cur, C Vec2 &next) {return GetTangentDir(prev-cur, next-cur);} Vec GetTangent(C Vec &prev, C Vec &cur, C Vec &next) {return GetTangentDir(prev-cur, next-cur);} Vec4 GetTangent(C Vec4 &prev, C Vec4 &cur, C Vec4 &next) {return GetTangent (prev , next );} /******************************************************************************/ VecB4 Lerp(C VecB4 &from, C VecB4 &to, Flt step) { Flt step1=1-step; return VecB4(Mid(RoundPos(from.x*step1 + to.x*step), 0, 255), Mid(RoundPos(from.y*step1 + to.y*step), 0, 255), Mid(RoundPos(from.z*step1 + to.z*step), 0, 255), Mid(RoundPos(from.w*step1 + to.w*step), 0, 255)); } VecB4 Lerp(C VecB4 &a, C VecB4 &b, C VecB4 &c, C Vec &blend) { return VecB4(Mid(RoundPos(a.x*blend.x + b.x*blend.y + c.x*blend.z), 0, 255), Mid(RoundPos(a.y*blend.x + b.y*blend.y + c.y*blend.z), 0, 255), Mid(RoundPos(a.z*blend.x + b.z*blend.y + c.z*blend.z), 0, 255), Mid(RoundPos(a.w*blend.x + b.w*blend.y + c.w*blend.z), 0, 255)); } /******************************************************************************/ Flt LerpR(C Vec2 &from, C Vec2 &to, C Vec2 &value) {Int i=Abs(to-from).maxI(); return LerpR(from.c[i], to.c[i], value.c[i]);} Dbl LerpR(C VecD2 &from, C VecD2 &to, C VecD2 &value) {Int i=Abs(to-from).maxI(); return LerpR(from.c[i], to.c[i], value.c[i]);} Flt LerpR(C Vec &from, C Vec &to, C Vec &value) {Int i=Abs(to-from).maxI(); return LerpR(from.c[i], to.c[i], value.c[i]);} Dbl LerpR(C VecD &from, C VecD &to, C VecD &value) {Int i=Abs(to-from).maxI(); return LerpR(from.c[i], to.c[i], value.c[i]);} Flt LerpR(C Vec4 &from, C Vec4 &to, C Vec4 &value) {Int i=Abs(to-from).maxI(); return LerpR(from.c[i], to.c[i], value.c[i]);} Dbl LerpR(C VecD4 &from, C VecD4 &to, C VecD4 &value) {Int i=Abs(to-from).maxI(); return LerpR(from.c[i], to.c[i], value.c[i]);} /******************************************************************************/ Flt LerpAngle(Flt from, Flt to, Flt step) {return step*AngleDelta(from, to)+from;} /******************************************************************************/ void Lerp4Weights(Vec4 &weights, Flt step) { Flt s =step, s2=step*step, s3=step*step*step; weights.x=(( -TANGENT)*(s3+s ) + (2*TANGENT )*s2 ); weights.y=(( 2-TANGENT)* s3 + ( TANGENT-3)*s2 + 1); //weights.z=((-2+TANGENT)* s3 - (2*TANGENT-3)*s2 + TANGENT*s ); don't calculate it manually, use the fact that the sum is always equal to 1, using Z component was the fastest version tested weights.w=( TANGENT *(s3-s2) ); weights.z=1-weights.x-weights.y-weights.w; } /******************************************************************************/ Flt Lerp4(Flt v0, Flt v1, Flt v2, Flt v3, Flt step) { Vec4 weights; Lerp4Weights(weights, step); return v0*weights.x + v1*weights.y + v2*weights.z + v3*weights.w; } Vec2 Lerp4(C Vec2 &v0, C Vec2 &v1, C Vec2 &v2, C Vec2 &v3, Flt step) { Vec4 weights; Lerp4Weights(weights, step); return v0*weights.x + v1*weights.y + v2*weights.z + v3*weights.w; } Vec Lerp4(C Vec &v0, C Vec &v1, C Vec &v2, C Vec &v3, Flt step) { Vec4 weights; Lerp4Weights(weights, step); return v0*weights.x + v1*weights.y + v2*weights.z + v3*weights.w; } Vec4 Lerp4(C Vec4 &v0, C Vec4 &v1, C Vec4 &v2, C Vec4 &v3, Flt step) { Vec4 weights; Lerp4Weights(weights, step); return v0*weights.x + v1*weights.y + v2*weights.z + v3*weights.w; } /******************************************************************************/ Flt LerpTan(Flt from, Flt to, Flt step, Flt tan0, Flt tan1) { return step*step*step * ((from-to)* 2 + tan0 + tan1) + step*step * ((from-to)*-3 - 2*tan0 - tan1) + step * ( tan0 ) + from; } Vec2 LerpTan(C Vec2 &from, C Vec2 &to, Flt step, C Vec2 &tan0, C Vec2 &tan1) { Flt s2=step*step, s3=step*step*step; return from + (from-to) * (2*s3 - 3*s2 ) + tan0 * ( s3 - 2*s2 + step) + tan1 * ( s3 - s2 ); } Vec LerpTan(C Vec &from, C Vec &to, Flt step, C Vec &tan0, C Vec &tan1) { Flt s2=step*step, s3=step*step*step; return from + (from-to) * (2*s3 - 3*s2 ) + tan0 * ( s3 - 2*s2 + step) + tan1 * ( s3 - s2 ); } Vec4 LerpTan(C Vec4 &from, C Vec4 &to, Flt step, C Vec4 &tan0, C Vec4 &tan1) { Flt s2=step*step, s3=step*step*step; return from + (from-to) * (2*s3 - 3*s2 ) + tan0 * ( s3 - 2*s2 + step) + tan1 * ( s3 - s2 ); } /******************************************************************************/ // ADJUST VALUE /******************************************************************************/ void AdjustValDir(Flt &value, Int dir, Flt dv) { if(dir) // adjust by direction { value+=dir*dv; Clamp(value, -1, 1); }else if(Int s=Sign(value)) // move towards zero { value-=s*dv; if(Sign(value)!=s)value=0; } } void AdjustValDir(Flt &value, Int dir, Flt change, Flt reset) { if(dir) // adjust by direction { value+=dir*change; Clamp(value, -1, 1); }else if(Int s=Sign(value)) // move towards zero { value-=s*reset; if(Sign(value)!=s)value=0; } } /******************************************************************************/ void AdjustValBool(Flt &value, Bool on, Flt dv) { if(on)value+=dv; // increase else value-=dv; // decrease SAT(value); // saturate } void AdjustValBool(Flt &value, Bool on, Flt inc, Flt dec) { if(on)value+=inc; // increase else value-=dec; // decrease SAT(value); // saturate } /******************************************************************************/ void AdjustValBoolSet(Flt &value, Bool on, Bool set, Flt dv ) {if(set)value=on;else AdjustValBool(value, on, dv );} // immediately set or smoothly adjust void AdjustValBoolSet(Flt &value, Bool on, Bool set, Flt inc, Flt dec) {if(set)value=on;else AdjustValBool(value, on, inc, dec);} // immediately set or smoothly adjust /******************************************************************************/ void AdjustValTime(Flt &value, Flt target, Flt exponent, Flt dt) {value=Pow(exponent, dt)*(value-target)+target;} // smoothly adjust value to 'target' according to 'exponent' and 'time delta' void AdjustValTime(Dbl &value, Dbl target, Flt exponent, Flt dt) {value=Pow(exponent, dt)*(value-target)+target;} // smoothly adjust value to 'target' according to 'exponent' and 'time delta' void AdjustValTime(Vec2 &value, C Vec2 &target, Flt exponent, Flt dt) {value=Pow(exponent, dt)*(value-target)+target;} // smoothly adjust value to 'target' according to 'exponent' and 'time delta' void AdjustValTime(Vec &value, C Vec &target, Flt exponent, Flt dt) {value=Pow(exponent, dt)*(value-target)+target;} // smoothly adjust value to 'target' according to 'exponent' and 'time delta' void AdjustValTime(VecD &value, C VecD &target, Flt exponent, Flt dt) {value=Pow(exponent, dt)*(value-target)+target;} // smoothly adjust value to 'target' according to 'exponent' and 'time delta' /******************************************************************************/ void AdjustValTarget(Flt &value, Flt target, Flt dv) { if(value>target){value-=dv; if(valuetarget)value=target;} // move towards 'target' by decreasing 'value' by 'dv' } void AdjustValTarget(Dbl &value, Dbl target, Dbl dv) { if(value>target){value-=dv; if(valuetarget)value=target;} // move towards 'target' by decreasing 'value' by 'dv' } void AdjustValTarget(Flt &value, Flt target, Flt inc, Flt dec) { if(value>target){value-=dec; if(valuetarget)value=target;} // move towards 'target' by decreasing 'value' by 'dec' } void AdjustValTarget(Vec2 &value, C Vec2 &target, Flt dv) { AdjustValTarget(value.x, target.x, dv); AdjustValTarget(value.y, target.y, dv); } void AdjustValTarget(Vec &value, C Vec &target, Flt dv) { AdjustValTarget(value.x, target.x, dv); AdjustValTarget(value.y, target.y, dv); AdjustValTarget(value.z, target.z, dv); } void AdjustValTarget(VecD &value, C VecD &target, Dbl dv) { AdjustValTarget(value.x, target.x, dv); AdjustValTarget(value.y, target.y, dv); AdjustValTarget(value.z, target.z, dv); } /******************************************************************************/ // SMOOTH OFFSET /******************************************************************************/ Flt SmoothOffset(Flt &offset, Flt max_length) { Flt ofs=offset; if(Flt length=Abs(ofs)) { ofs=Sign(ofs); const Flt factor=1.0f/PI_2; // mul by this because we're operating on 'Sin' below, and for low values "Sin(x*PI_2)" is equal to "x*PI_2/PI_2" (tested via "Functions" tool) length*=factor; if(length>max_length){offset=ofs*(max_length/factor); length=max_length;} Flt frac=length/max_length; frac=Sin(frac*PI_2); ofs*=max_length*frac; } return ofs; } Vec2 SmoothOffset(Vec2 &offset, Flt max_length) { Vec2 ofs=offset; if(Flt length=ofs.normalize()) { const Flt factor=1.0f/PI_2; // mul by this because we're operating on 'Sin' below, and for low values "Sin(x*PI_2)" is equal to "x*PI_2/PI_2" (tested via "Functions" tool) length*=factor; if(length>max_length){offset=ofs*(max_length/factor); length=max_length;} Flt frac=length/max_length; frac=Sin(frac*PI_2); ofs*=max_length*frac; } return ofs; } Vec SmoothOffset(Vec &offset, Flt max_length) { Vec ofs=offset; if(Flt length=ofs.normalize()) { const Flt factor=1.0f/PI_2; // mul by this because we're operating on 'Sin' below, and for low values "Sin(x*PI_2)" is equal to "x*PI_2/PI_2" (tested via "Functions" tool) length*=factor; if(length>max_length){offset=ofs*(max_length/factor); length=max_length;} Flt frac=length/max_length; frac=Sin(frac*PI_2); ofs*=max_length*frac; } return ofs; } /******************************************************************************/ // SMOOTH CURVES /******************************************************************************/ Flt SmoothSqr(Flt x) { if(x<=0 )return 0; if(x>=1 )return 1; if(x<=0.5f)return 2*Sqr( x); return 1-2*Sqr(1-x); // 1-2*Sqr(1-x) -> 1-2*(1+x*x-2*x) -> 1 - 2 - 2*x*x + 4*x -> -2*x*x + 4*x - 1 (however this formula has more operations) } Flt SmoothCube(Flt x) { if(x<=0)return 0; if(x>=1)return 1; return _SmoothCube(x); } Flt SmoothCubeInv(Flt y) { if(y<=0)return 0; if(y>=1)return 1; return 0.5f-Sin(asinf(1-2*y)/3); } Flt SmoothCube2(Flt x) { if(x<=0 )return 0; if(x>=1 )return 1; if(x<=0.5f)return 4*Cube( x); return 1-4*Cube(1-x); } Flt _SmoothQuintic(Flt x) { return x*x*x*(x*(x*6-15)+10); } Flt SmoothSextic(Flt x) { if(x<=0)return 0; if(x>=1)return 1; x*=x; return (4.0f/9)*x*x*x - (17.0f/9)*x*x + (22.0f/9)*x; } Flt SmoothSin(Flt x) { if(x<=0)return 0; if(x>=1)return 1; return 0.5f-0.5f*Cos(x*PI); } Flt SmoothPow(Flt x, Flt pow) { if(x<=0 )return 0; if(x>=1 )return 1; if(x<=0.5f)return 0.5f*Pow( 2*x, pow); return 1-0.5f*Pow(2-2*x, pow); } Flt SmoothPinch(Flt x, Flt pinch) { if(x<=0 )return 0; if(x>=1 )return 1; if(x<=0.5f)return 0.5f*Pinch( 2*x, pinch); return 1-0.5f*Pinch(2-2*x, pinch); } /******************************************************************************/ // SIGMOID /******************************************************************************/ Flt gd(Flt x) {return Atan(sinhf(x));} Dbl gd(Dbl x) {return Atan(sinh (x));} Flt SigmoidExp (Flt x) {return 2/(1+expf(-x))-1;} Flt SigmoidDiv (Flt x) {return x/(1+x);} Flt SigmoidAtan (Flt x) {return Atan(PI_2*x)*2/PI;} Flt SigmoidSqrt (Flt x) {return x/SqrtFast(1+x*x);} Flt SigmoidSqrtInv(Flt y) {return y/SqrtFast(1-y*y);} Flt SigmoidGd (Flt x) {return (2/PI)*gd(PI_2*x);} Flt SigmoidTanh (Flt x) {return tanhf(x);} Flt SigmoidErf (Flt x) {return erff((SqrtFast(PI)/2)*x);} /******************************************************************************/ // BLENDING /******************************************************************************/ Flt BlendSqr(Flt x) { x*=x; return (x<1) ? 1-x : 0; } Flt BlendSmoothSqr(Flt x) { ABS(x); return (x<=0.5f) ? 1-2*Sqr( x) : (x< 1 ) ? 2*Sqr(1-x) : 0; } Flt BlendSmoothCube(Flt x) // !! avoid changing this formula, but if it's changed, then adjust 'BlendSmoothCubeSum' and 'BlendSmoothCubeSumHalf' !! { ABS(x); return (x<1) ? 1-(3-2*x)*x*x : 0; } Flt BlendSmoothSextic(Flt x) { x*=x; return (x<1) ? (-4.0f/9)*x*x*x + (17.0f/9)*x*x - (22.0f/9)*x + 1 : 0; } Flt BlendSmoothSin(Flt x) { ABS(x); return (x<1) ? Cos(x*PI)*0.5f+0.5f : 0; } /******************************************************************************/ Flt Gaussian(Flt x) { return expf(-x*x); } /******************************************************************************/ // EQUATIONS /******************************************************************************/ Int Polynominal1(Flt a, Flt b, Flt &x ) {if(!a)return Polynominal0(b); x=-b/a; return 1;} Int Polynominal1(Dbl a, Dbl b, Dbl &x ) {if(!a)return Polynominal0(b); x=-b/a; return 1;} Int Polynominal2(Flt a, Flt b, Flt c, Flt &x0, Flt &x1) { if(!a)return Polynominal1(b, c, x0); Flt delta=b*b-4*a*c; if( delta<0)return 0; a+=a; if(!delta){x0= -b /a; return 1;} delta=Sqrt(delta); x0=(-b-delta)/a; x1=(-b+delta)/a; return 2; } Int Polynominal2(Dbl a, Dbl b, Dbl c, Dbl &x0, Dbl &x1) { if(!a)return Polynominal1(b, c, x0); Dbl delta=b*b-4*a*c; if( delta<0)return 0; a+=a; if(!delta){x0= -b /a; return 1;} delta=Sqrt(delta); x0=(-b-delta)/a; x1=(-b+delta)/a; return 2; } Int Polynominal3(Flt a, Flt b, Flt c, Flt d, Flt &x0, Flt &x1, Flt &x2) { if(!a)return Polynominal2(b, c, d, x0, x1); Flt f=c/a - b*b/(a*a*3), g=2*b*b*b/(a*a*a*27) - b*c/(a*a*3) + d/a, h=g*g/4 + f*f*f/27; if(h>0) { g/=-2; h =Sqrt(h); x0=Cbrt(g+h)+Cbrt(g-h) - b/(3*a); return 1; }else if(!f && !g && !h) { x0=x1=x2=-Cbrt(d/a); return 3; }else if(h<=0) { Flt i= g*g/4 - h, p=-b/(3*a), m, n; CosSin(m, n, Acos(-0.5f*g/Sqrt(i))/3); i=Pow(i, 1.0f/6); n*=SQRT3; x0=2*i* m +p; x1= -i*(m+n)+p; x2= -i*(m-n)+p; return 3; } return 0; } Int Polynominal3(Dbl a, Dbl b, Dbl c, Dbl d, Dbl &x0, Dbl &x1, Dbl &x2) { if(!a)return Polynominal2(b, c, d, x0, x1); Dbl f=c/a - b*b/(a*a*3), g=2*b*b*b/(a*a*a*27) - b*c/(a*a*3) + d/a, h=g*g/4 + f*f*f/27; if(h>0) { g/=-2; h =Sqrt(h); x0=Cbrt(g+h)+Cbrt(g-h) - b/(3*a); return 1; }else if(!f && !g && !h) { x0=x1=x2=-Cbrt(d/a); return 3; }else if(h<=0) { Dbl i= g*g/4 - h, p=-b/(3*a), m, n; CosSin(m, n, Acos(-0.5*g/Sqrt(i))/3); i=Pow(i, 1.0/6); n*=SQRT3; x0=2*i* m +p; x1= -i*(m+n)+p; x2= -i*(m-n)+p; return 3; } return 0; } /******************************************************************************/ Int Solve(Flt a1, Flt a2, Flt b1, Flt b2, Flt c1, Flt c2, Flt &x, Flt &y) { // x*a1 + y*b1 = c1 // x*a2 + y*b2 = c2 // W = | a1 b1 | = a1*b2 - b1*a2, Wx = | c1 b1 | = c1*b2 - c2*b1, Wy = | a1 c1 | = a1*c2 - a2*c1 // | a2 b2 | | c2 b2 | | a2 c2 | if(Flt W=a1*b2 - b1*a2) { Flt Wx=c1*b2 - c2*b1, Wy=a1*c2 - a2*c1; if(!Wx && !Wy)return -1; x=Wx/W; y=Wy/W; return 1; } return 0; } Int Solve(Dbl a1, Dbl a2, Dbl b1, Dbl b2, Dbl c1, Dbl c2, Dbl &x, Dbl &y) { // x*a1 + y*b1 = c1 // x*a2 + y*b2 = c2 // W = | a1 b1 | = a1*b2 - b1*a2, Wx = | c1 b1 | = c1*b2 - c2*b1, Wy = | a1 c1 | = a1*c2 - a2*c1 // | a2 b2 | | c2 b2 | | a2 c2 | if(Dbl W=a1*b2 - b1*a2) { Dbl Wx=c1*b2 - c2*b1, Wy=a1*c2 - a2*c1; if(!Wx && !Wy)return -1; x=Wx/W; y=Wy/W; return 1; } return 0; } /****************************************************************************** 64 BIT FLOAT - S EEEEEEEEEEE FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF (S-1 bit Sign, E-11 bit Exponent, F-52 bit Fraction) 32 BIT FLOAT - S EEEEEEEE FFFFFFFFFFFFFFFFFFFFFFF (S-1 bit Sign, E- 8 bit Exponent, F-23 bit Fraction) 16 BIT FLOAT - S EEEEE FFFFFFFFFF (S-1 bit Sign, E- 5 bit Exponent, F-10 bit Fraction) 32: Exponent | Fraction Zero | Fraction Non-Zero | Equation 0 | +0 -0 | subnormal numbers | Pow(-1, sign) * Pow(2, -126) * 0.Fraction 1..254 | value | value | Pow(-1, sign) * Pow(2, exponent-127) * 1.Fraction 255 | +Inf -Inf | NaN | 16: Exponent | Fraction Zero | Fraction Non-Zero | Equation 0 | +0 -0 | subnormal numbers | Pow(-1, sign) * Pow(2, -14) * 0.Fraction 1..30 | value | value | Pow(-1, sign) * Pow(2, exponent-15) * 1.Fraction 31 | +Inf -Inf | NaN | /******************************************************************************/ Half::operator Flt()C { UInt f; Int e=(data&0x7C00); if(e== 0)f=((data&0x8000)<<16);else if(e==0x7C00)f=((data&0x8000)<<16) | (0xFF <<23) | ((data&0x03FF)<<13);else f=((data&0x8000)<<16) | ((e+0x1C000)<<13) | ((data&0x03FF)<<13); // (127-15)<<10 = 0x1C000 return (Flt&)f; } /*Half::Half(Bool b) {data=(b ? 15360 : 0);} Half::Half(Int i) : Half(Flt(i)) {} Half::Half(UInt i) : Half(Flt(i)) {}*/ Half::Half(Flt f) { UInt u=(U32&)f; // Int e=((u>>23)&0xFF)+(15-127); // if(e<= 0)data= (u>>16)&0x8000;else // if(e> 31)data=((u>>16)&0x8000) | (31<<10) | ((u>>13)&0x03FF);else // data=((u>>16)&0x8000) | ( e<<10) | ((u>>13)&0x03FF); Int e =(u&0x7F800000); if( e<= 0x38000000)data= (u>>16)&0x8000;else if( e>= 0x47800000)data=((u>>16)&0x8000) | ( 31 <<10) | ((u>>13)&0x03FF);else data=((u>>16)&0x8000) | ((e-0x38000000)>>13) | ((u>>13)&0x03FF); } /******************************************************************************/ void DecRealByBit(Flt &r) { UInt u =(U32&)r, f =(u&0x7FFFFF), exp =((u>>23)&0xFF); Bool sign=(u>>31); if( !sign) { f--; f&=0x7FFFFF; switch(f) { case 0x7FFFFF: if(exp)exp--;else{sign=true; /*exp=0;*/ f=1;} break; // -eps } }else if(exp<0xFF) { f++; f&=0x7FFFFF; if(!f)exp++; } (U32&)r=(sign<<31)|(exp<<23)|f; } void IncRealByBit(Flt &r) { UInt u =(U32&)r, f =(u&0x7FFFFF), exp =((u>>23)&0xFF); Bool sign=(u>>31); if( sign) { f--; f&=0x7FFFFF; switch(f) { case 0x7FFFFF: if(exp)exp--;else{sign=false; /*exp=0;*/ f=1;} break; // eps case 0x000000: if(!exp){sign=false; /*exp=0;*/ f=0;} break; // fix -0 -> 0 } }else if(exp<0xFF) { f++; f&=0x7FFFFF; if(!f)exp++; } (U32&)r=(sign<<31)|(exp<<23)|f; } void DecRealByBit(Dbl &r) { ULong u =(U64&)r, f =(u&0xFFFFFFFFFFFFF), exp =((u>>52)&0x7FF), sign=(u>>63); if( !sign) { f--; f&=0xFFFFFFFFFFFFF; switch(f) { case 0xFFFFFFFFFFFFF: if(exp)exp--;else{sign=true; /*exp=0;*/ f=1;} break; // -eps } }else if(exp<0x7FF) { f++; f&=0xFFFFFFFFFFFFF; if(!f)exp++; } (U64&)r=(sign<<63)|(exp<<52)|f; } void IncRealByBit(Dbl &r) { ULong u =(U64&)r, f =(u&0xFFFFFFFFFFFFF), exp =((u>>52)&0x7FF), sign=(u>>63); if( sign) { f--; f&=0xFFFFFFFFFFFFF; switch(f) { case 0xFFFFFFFFFFFFF: if(exp)exp--;else{sign=false; /*exp=0;*/ f=1;} break; // eps case 0x0000000000000: if(!exp){sign=false; /*exp=0;*/ f=0;} break; // fix -0 -> 0 } }else if(exp<0x7FF) { f++; f&=0xFFFFFFFFFFFFF; if(!f)exp++; } (U64&)r=(sign<<63)|(exp<<52)|f; } /******************************************************************************/ // Standard Deviational Ellipse // http://resources.esri.com/help/9.3/arcgisengine/java/gp_toolref/spatial_statistics_tools/how_directional_distribution_colon_standard_deviational_ellipse_spatial_statistics_works.htm void AvgDirU(Vec2 &dir, C MemPtr &points) { Flt x2=0, y2=0, xy=0; REPA(points) { C Vec2 &p=points[i]; x2+=Sqr(p.x); y2+=Sqr(p.y); xy+=p.x*p.y; } Flt a=x2-y2, c=2*xy, b=Dist(a, c); dir.set(a+b, c); // if we would normalize, then dir.x=cos, dir.y=sin } /* void AvgDirU(Vec &dir, C MemPtr &points); // get average direction from array of points, points should be located around center Vec (0,0,0), direction is not normalized void AvgDirU(Vec &dir, C MemPtr &points) { this is wrong Flt x2=0, y2=0, z2=0, xy=0, xz=0; REPA(points) { C Vec &p=points[i]; x2+=Sqr(p.x); y2+=Sqr(p.y); z2+=Sqr(p.z); xy+=p.x*p.y; xz+=p.x*p.z; } Flt a=x2-y2, c=2*xy, b=Dist(a, c); dir.x=a+b; dir.y=c; dir.z=c; }*/ /****************************************************************************** void MathAssert() { DEBUG_ASSERT(FloorSpecial(0)==-1 && FloorSpecial(0.01f)==0 && FloorSpecial(0.5f)==0 && FloorSpecial(0.99f)==0 && FloorSpecial(1.0f)==0, "FloorSpecial"); DEBUG_ASSERT( CeilSpecial(0)== 1 && CeilSpecial(0.01f)==1 && CeilSpecial(0.5f)==1 && CeilSpecial(0.99f)==1 && CeilSpecial(1.0f)==2, "CeilSpecial"); } /******************************************************************************/ } /******************************************************************************/