{ ---------------------------------------------------------------------

Title: KSP Math Library

Author: R.D.Villwock aka 'Big Bob'

First Written: October 2, 2006

Current Version: 2.05

Last Modified: April 19, 2008

----------------------------------------------------------------------}

{ This module provides a number of useful math functions not included with the basic
  KSP scripting language. These functions can support scripts that perform various
  forms of equal-power crossfading and other forms of sophisticated volume control.
  The Log/Exp routines can also be useful for various tempo/pitch calculations.
  
  NOTE: This Library is written as an 'import' module to be compiled with your main,
  host script using Nils Liberg's KScript Editor V1.25.6 or higher. Importing this 
  Library will not 'bloat your code'. Nothing will be added to your code except
  what is needed to support the functions you actually use. If you don't use any
  of the Library's functions, no code or data declarations will be added.
  
  With the exception of the CVHex routine, all functions automatically declare any
  needed data structures so all you need do to use any Library function is to
  simply reference it in your code. To use CVHex however, you need to include a
  call to UseHex in your ICB before the first use of CVHex.
  
  ----------------------------- The Standard Functions ----------------------------- 
  ATFade(cv,rng,atn)    This function maps a MIDI CC value, 'cv', to control overall

                        attenuation, 'atn', with an audio-taper and range adjustment

  Ang90                 This constant (currently = 1000 dg) is available if your

                        script references any Trig functions. It is useful for angle

                        conversion from arbitrary right-angle units. For example,

                        if you want to use a MIDI CC to cover tha range from 0..90

                        degrees, you can use:  angle = CC*Ang90/127.

  exp(lgx,X)            This function returns an integer 'X' which is the natural
                        antilog of the scaled input 'lgx', ie  X = e^(lgx/1000000)
                        where 0.0 < lgx < 21.4875 * 10^6

  Exp2(lgx,X)           This function returns an integer 'X' which is the binary
                        antilog of the scaled input 'lgx', ie X = 2^(lgx/1000000)
                        where: 0.0 < lgx < 30.9999 * 10^6


  Exp10(lgx,X)          This function returns an integer 'X' which is the common

                        antilog of the scaled input 'lgx', ie X = 10^(X/1000000)
                        where: 0.0 < X < 9.3319 * 10^6

                         
  Ln(X,lgx)             This function accepts an input of 0 < 'X' < 2^16 and returns

                        the base e logarithm, 'lgx', scaled by 10,000. This routine

                        is less accurate and accepts a narrower input range than XLn

                        but executes faster.

                        

  Log2(X,lgx)           This function accepts an input of 0 < 'X' < 2^16 and returns

                        the base 2 logarithm, 'lgx', scaled by 10,000. This routine

                        is less accurate and accepts a narrower input range than XLog2

                        but Log2 executes faster.

                        

  Log10(X,lgx)          This function accepts an input of 0 < 'X' < 2^16 and returns

                        the base 10 logarithm, 'lgx', scaled by 10,000. This routine

                        is less accurate and accepts a narrower input range than XLog10

                        but Log10 executes faster.
                        
  Muted                 This constant (currently -200,000 mdb) is available if your
                        script references ep_to_mdb, Get_db, VR_to_mdb, or ATFade.
                        
  SinCos(ang,sin,cos)   This function accepts an input angle in the first quadrant
                        (0..1000 deci-grads) and returns both sine and cosine scaled
                        by 10,000. For the full integer range of angles, use FSinCos.
                        
  Tangent(ang,tan)      This function accepts an input angle in the first quadrant

                        (0..1000 deci-grads) and returns the tangent scaled by 10,000

                        For the full integer range of angles, use FTangent.
                        
----------------------------------- Extended Functions -------------------------------
  XLn(X,lgx)            This extended precision function also accepts any positive
                        integer (0 < 'X' < 2^31) as input. XLn returns the natural
                        (base e) logarithm, 'lgx', scaled by 1,000,000. However,
                        this routine is not as fast as Ln.

                        

  XLog2(X,lgx)          This extended precision function also accepts any positive

                        integer (0 < 'X' < 2^31) as input. XLog2 returns the binary

                        (base 2) logarithm, 'lgx', scaled by 1,000,000. However,

                        this routine is not as fast as Log2.
                        

  XLog10(X,lgx)         This extended precision function also accepts any positive

                        integer (0 < 'X' < 2^31) as input. XLog10 returns the common

                        (base 10) logarithm, 'lgx', scaled by 1,000,000. However,

                        this routine is not as fast as Log10.

  XSinCos(ang,sin,cos)  This extended function accepts any input angle in deci-grads

                        and returns both the sine and cosine scaled by 10,000

                        

  XTangent(ang,tan)     This extended function accepts any input angle in deci-grads

                        and returns the tangent scaled by 10,000, tan(1000) = 100M

------------------------------ Utility Conversion Routines -------------------------
  CVHex(num,str)        Converts any integer 'num' to its equivalent 8-character
                        Hexadecimal String, placing the result in 'str'. To use this
                        function you must first invoke the UseHex data function
                        somewhere in your ICB. This function is provided mostly for
                        use in testing rather than for final script support.
                        
  ep_to_mdb(N,Vol)      This function converts an Engine Parameter 'N' to equivalent

                        Volume using the formula: 'Vol' = 20000*log(N/Nm)^3 + 12000

                        where 'Vol' is in mdb and Nm = 1,000,000

  Get_db(VR,Vol)        This function takes an input scaled-ratio, 'VR' = 10000(V/V0)

                        and returns the equivalent volume, 'Vol', in mdb using the

                        formula Vol = 20000*log(V/V0). For continuity, this legacy

                        routine limits the input ratio V/V0 to 1.0000 but otherwise

                        performs the same function as VR_mdb. While Get_db is not as

                        accurate as VR_mdb, it executes faster since it is based on

                        DeVos' logarithm approximation.
                        
  mdb_to_ep(Vol,N)      This function converts Volume in mdb to its equivalent
                        ep value N using the formula: N = Nm*10^(Vol - 12000)/60000 
                        where 'Vol' is in mdb and Nm = 1,000,000.

  VR_to_ep(VR,N)        This function takes an input scaled-ratio, 'VR' = 10000(V/V0)

                        and returns the corresponding Engine Parameter 'N' (0..1000000),

                        where N = 1,000,000 when V/V0 = 4.0

  VR_to_mdb(VR,Vol)     This function takes an input scaled-ratio, 'VR' = 10000(V/V0)

                        and returns the equivalent volume, 'Vol', in mdb using the

                        formula Vol = 20000*log|V/V0|. This function is the same as

                        the legacy function named Get_db except that V/V0 can be as

                        high as 4.0000 and the routine is more accurate. However, the

                        legacy routine is faster since it uses DeVos' approximation.

  ---------------------------------------------------------------------------------

                           Library Functions Source Code

  ---------------------------------------------------------------------------------
  ATFade(cv,rng,atn)     Audio-Taper MIDI Fader Function 

  

  'cv'  = control value input 0..127

  'rng' = control range input 0..100%

  'atn' = attenuation output -200,000..0 mdb

  

  Accepts a linear control value, and returns an audio-taper

  attenuation value. At 'cv' = 127, 'atn' = 0 db. As 'cv' swings

  from 127..0, 'atn' increases from 0 db toward some maximum

  negative value determined by the 'rng' input. At 'rng' = 100%,

  'atn' swings through the full range from 0..-200 db (muted).

  For 'rng' values of 100..0%, the max 'atn' value (at 'cv' = 0)

  changes from -200..0 db. Thus, with 'rng' at 0%, 'atn' = 0 db

  for any 'cv' value (ie 'cv' has NO control range).

    

  For the full range of attenuation, 'rng' = 100%, a 3-segment

  linear db relationship is used. The 'working' zone from 0..-25 db

  is assigned to 'cv' = 127..50. The 'soft' zone from -25..-60 db

  is assigned to 'cv' = 50..17, and the 'fade-out' zone from

  -60..-200 db is assigned to 'cv' = 17..0. This provides smooth

  and accurate volume control throughout the entire audio range. }



function ATFade(cv,rng,atn)
  DefineMuted      { Declare Muted constant } 

  declare const 30K := 30000  

  declare const 50K := 50000

  declare const 60K := 60000

  declare const 70K := 70000

  declare const 280K := 280000  

  declare cx  { Range-mapped 'cv' }

     { node 0: c0 = 0,  a0 = -200000 mdb

       node 1: c1 = 17, a1 = -60000 mdb

       node 2: c2 = 50, a2 = -25000 mdb }

         

  cx := 127 + rng*(cv - 127)/100  { Set range for cx }

  select cx

    case 50 to 127      { 'Working' zone, -25..0 db }

      atn := (50K*(cx - 127) - 72)/154

      { atn := a2*(127 - cx)/(127 - c2) rounded }

    case 17 to 49       { Soft zone, -60..-25 db }

      atn := -60K + (70K*(cx - 17) + 33)/66

      {  atn := a1 + (a2 - a1)*(cx - c1)/(c2 - c1) rounded }

    case 0 to 16        { Fade-out zone, -200..60 db }

      atn := Muted + (280K*cx + 17)/34

      { atn := a0 + (a1 - a0)*cx/c1 rounded }

  end select

end function { ATFade }

{ -------------------------------------------------------------------------

    exp(lgx,X))  Natural AntiLog Function (base e) 

  

            'lgx' = Input value (scaled by 1000000)

              'X' = e^(lgx/1000000) ie Base e Antilog of input when
                    input is in the range 0.0 <= lgx <= 21.487562

            

    NOTE: For lgx > 21.487562, result X is clamped to MaxInt = 2,147,483,647

          For lgx < 0.0, result X is set to zero.  }



function exp(lgx,X)

  declare const MaxInt := 0x7FFFFFFF  { Maximum positive integer, 2^31 - 1 }
  declare rbx  { Re-based value of lgx }

  

  if lgx > 21487562

    X := MaxInt

  else  { Rebase rbx = lgx/Ln(2) }
    rbx := ((lgx/1000000)*70692057 + 25)/49 + ((lgx mod 1000000)*1649 + 571)/1143

    XP2.Core(rbx,X)  { Compute antilog of re-based lgx: 0.0 <= rbx < 31.0 }

  end if

end function { exp}


{ -------------------------------------------------------------------------

    Exp2(lgx,X))  Binary AntiLog Function (base 2) 

  

            'lgx' = Input value (scaled by 1000000)

              'X' = 2^(lgx/1000000) ie Base 2 Antilog of input when
                    input is in the range 0.0 <= lgx <= 30.999999              

            

    NOTE: For lgx > 30.999999, result X is clamped to MaxInt = 2,147,483,647
          For lgx < 0.0, result X is set to zero. }
          
function Exp2(lgx,X)
  declare const MaxInt := 0x7FFFFFFF  { Maximum positive integer, 2^31 - 1 }
  
  if lgx > 30999999
    X := MaxInt
  else
    XP2.Core(lgx,X)  { Compute the antilog of 0.0 <= lgx <= 30.999999 }
  end if
end function { Exp2 }

{ -------------------------------------------------------------------------

    Exp10(lgx,X))  Common AntiLog Function (base 10) 

  

            'lgx' = Input value (scaled by 1000000)

              'X' = 10^(lgx/1000000) ie Base 10 Antilog of input when

                    input is in the range 0.0 <= lgx <= 9.331929

            

    NOTE: For lgx > 9.331929, result X is clamped to MaxInt = 2,147,483,647

          For lgx < 0.0, result X is set to zero.  }



function Exp10(lgx,X)

  declare const MaxInt := 0x7FFFFFFF  { Maximum positive integer, 2^31 - 1 }
  declare rbx  { Re-based value of lgx }

  

  if lgx > 9331929

    X := MaxInt

  else  { Rebase rbx = lgx/log(2) }
    rbx := ((lgx/1000000)*106301699 + 16)/32 + ((lgx mod 1000000)*2136 + 321)/643 
    XP2.Core(rbx,X)  { Compute antilog of re-based lgx: 0.0 <= rbx < 31.0 }

  end if

end function { Exp10 }

{ -------------------------------------------------------------------------

  Ln(X,lgx))    Natural Logarithm Function (base e)

  

    'X' = Input positive integer value 1..65,536

  'lgx' = Ln(x) scaled by 10,000  ie Ln(x) = 'lnx'/10,000



          NOTE: For X <= 0, result lgx is set to NAL = -100.0000

                For X > 65,535, result lgx is clamped to Ln(65536) = 11.0904 }

                

function Ln(X,lgx)

  Log2(X,lgx)   { Get binary log of X }

  if lgx > 0    { Convert Base if necessary }

    lgx := (lgx*7050 + 5085)/10171  { Ln(X) = Ln(2)*Log2(X) }   

  end if

end function { Ln }   



{ -------------------------------------------------------------------------

  Log2(X,lgx))   Binary Logarithm Function (base 2) 

  

              'X' = Input integer value 1..65,536

            'lgx' = Log2(X) scaled by 10,000  ie Log2(X) = 'lgx'/10,000

            

        NOTE: For X <= 0, result lgx is set to NAL = -100.0000

              For X > 65,536, result lgx is clamped to Log2(65536)=16.0000 

                          

  This routine uses a fast approximation devised by John DeVos (circa 2001). }

                     

function Log2(X,lgx)

  declare const LogScale := 10000      { Scaling for output logarithm }

  declare const MaxLog := 16*LogScale  { Log(Xmax), Xmax = 65536 }

  declare const NAL := -1000000   { Sentinel value used for undefined log, 'Not a Log' }

  declare const Lm3 := 6*LogScale { 3*LmScale }  

  { Pe = 0.014830 = 16,313/1.1M   1.1M = 110*LogScale }

  declare const SPe := 32626      { Scaled-Pe = Pe*220*LogScale }

  declare const SPe3 := 97584     { Scaled-Pe3 = 2.991*SPe }

  declare const PeScale := 220    { Extra Pe scale factor (above LogScale) }

  declare const PeRound := PeScale/2 { To half-adjust Pe error correction }   

  declare err { Error Correction term for mantissa }

  declare Lc  { Log2(x) characteristic }

  declare Lm  { Log2(x) mantissa approximation }

  declare X0  { X0 = 2^Lc, Lower end of mantissa 'interval' for input X } 

  

  if X <= 0         { First, filter out the extremes in the input X }

    lgx := NAL      { Log not defined for zero or negative X }

  else if X > 65535 { Equal or greater than maximum allowed X = 65536 }

    lgx := MaxLog   { Log2(65536)*10000 = maximum log output }

  else { Valid x range }  

    X0 := 0x8000    { Highest mantissa interval base }

    Lc := 15        { Max characteristic }

    while X .and. X0 = 0   { Loop to find highest 'on' bit }

      X0 := sh_right(X0,1) { Keep track of high bit value }

      dec (Lc)      { Keep track of high bit position from 15..0 }

    end while       { On loop exit, Lc = unscaled characteristic, X0 = 2^Lc }

    if X = X0

      lgx := Lc*LogScale { For exact powers of 2, mantissa = 0 so we're done }

    else { mantissa > 0 }

      Lm := Lm3*(X - X0)/(X + X0)           { Mantissa scaled by 2*LogScale }

      err := sh_right(X*SPe,Lc - 1) - SPe3  { SPe*2*x/X0 - SPe3 }

      err := SPe - err*err/SPe              { Error scaled by PeScale*LogScale }

      err := (sh_left(err,1) + PeRound)/PeScale     { Reduce scaling to 2*LogScale }

      lgx := Lc*LogScale + sh_right(Lm - err + 1,1) { Remove 2* scaling from mantissa }

    end if  { Final Result scaled only by LogScale }

  end if

end function { Log2 }



{ -------------------------------------------------------------------------

  Log10(X,lgx))    Common Logarithm Function (base 10)

  

         'X' = Input integer value 1..65,536

      'logx' = Log10(X) scaled by 10,000  ie Log10(X) = 'lgx'/10,000

  

           NOTE: For X <= 0, result lgx is set to NAL = -100.0000

                 For X > 65,536, result lgx is clamped to Log10(65536) = 4.8165 }

                 

function Log10(X,lgx)

  Log2(X,lgx)     { Get binary log of input X }

  if lgx > 0      { Convert base if necessary }

    lgx := (lgx*12655 + 21019)/42039  { lgx = Log10(2)*Log2(X) } 

  end if  

end function { Log10 }

{ ------------------------------------------------------------------  
  SinCos(ang,sin,cos)   First Quadrant Sine/Cosine Function 

  

    'ang' = 1st quadrant input angle in deci-grads (0..1000)

          NOTE: 1000 deci-grads = 90 degrees = PI/2 radians

    'sin' = Scaled Sine of input angle: 'sin'/10000 = sin(ang)

    'cos' = Scaled Cosine of input angle: 'cos'/10000 = cos(ang) 

  

  NORE: No input angle checking is performed. Caller should confine
        'ang' to range 0 <= 'ang' <= 1000 dg. Note also that this is
        the standard Cordic version for speed. For a smaller but slower
        version, see header comments of XSinCos at the end of this Library }
        
function SinCos(ang,sin,cos)
  BuildAngleTable  { Declare Ang90, AngScale, and scaled Atan table }
  declare const Ang45 := Ang90/2 { Initial 45 degree rotation }
  declare const cK := 607252935  { Cordic K*10^9 (15 or more iterations) }
  declare const DownScale := 100000 { Reduce scaling from 10^9 to 10^4 }
  declare const DownRound := DownScale/2  { Half adjust scale down }
  declare dx { Change in X resulting from current rotation }

  declare dy { Change in Y resulting from current rotation }  

  declare n  { Iteration index }

  declare x  { X-component of Cordic vector scaled by 10^9 }

  declare y  { Y-component of Cordic vector scaled by 10^9 }

  declare z  { Residual Cordic Angle }

  

  x := cK  { Pre-rotate reference vector K by 45 degrees so  }

  y := cK  {  that iteration for n = 0 need not be performed }

  z := (ang - Ang45)*AngScale { Scaled input angle referenced to 500 dg } 

  for n := 1 to TrigBits { Use 15 iterations ( about 4-digit precision }

    dx := sh_right(y,n)  { Rotated change in current x }

    dy := sh_right(x,n)  { Rotated change in current y }

    if z < 0  { Rotate vector clockwise }

      x := x + dx   { Compute new rotated X }

      y := y - dy   { Compute new rotated Y }

      z := z + AngTbl[n]  { Update residual angle }

    else      { Rotate vector counter-clockwise }  

      x := x - dx   { Compute new rotated X }

      y := y + dy   { Compute new rotated Y }

      z := z - AngTbl[n]  { Update residual angle } 

    end if

  end for

  sin := (y + DownRound)/DownScale  { Reduce scale to 10^4 and round }

  cos := (x + DownRound)/DownScale  { Both are positive in 1st quadrant }

end function { SinCos }         



{  -------------------------------------------------------------------

  Tangent(ang,tan)   First-Quadrant Tangent Function 

  

    'ang' = 1st quadrant input angle in deci-grads (0..1000)

          NOTE: 1000 deci-grads = 90 degrees = PI/2 radians

    'tan' = Scaled Tangent of input angle: 'tan'/10000 = tan(ang)
    
    NOTE: No input angle checking is performed. Caller should confine

          'ang' to range 0 <= 'ang' <= 1000 dg. Also note that the
          tangent function is discontinuous at +/-90 degrees, that is:
          tan(1000) -> infinity.    Since tan(999) is 6,250,000

          (625 scaled by 10K), for tan(1000) this routine returns an

          even larger value MaxTan = 100M as a sentinel for tan(1000) } 


function Tangent(ang,tan)
  declare const TrigScale := 10000         { Trig output scaling }
  declare const MaxTan := 10000*TrigScale  { Scaled sentinel for tan(1000) }
  declare sin
  declare cos
  
  if ang = Ang90   { Tangent of 1000 dg is infinite, use 100M }
    tan := MaxTan  {  as a sentinel  ie tan(1000) = 10K*Scale }
  else { ang < 90 degrees }
    SinCos(ang,sin,cos)      { Compute sine and cosine }
    tan := sin*TrigScale/cos { Output is scaled ratio of sin/cos }
  end if
end function { Tangent }

{ ---------------------------- Extended Functions ------------------------------
   These routines extend the input argument range and/or extend overall accuracy
  ------------------------------------------------------------------------------}

{   XLn(X,lgx)    Extended Natural Logarithm Function (base e)

  

             'X' = Input can be any integer value

           'lgx' = Ln(X) scaled by 1,000,000  ie Ln(x) = 'lgx'/1000000



          NOTE: For X <= 0, result lgx is set to NAL = -100.000000 }

                

function XLn(X,lgx)
  declare const NAL := -100000000 { Token for Ln(0) or negative X, -100.000000 } 
  declare c  { Log characteristic }
  declare m  { Log mantissa }

  

  if X <= 0
    lgx := NAL             { Log is not defined, return NAL }
  else
    XLg.Core(X,c,m)        { Get binary log of 1 <= X <= 2,147,483,647 }

    lgx := (c*49906597+36)/72 + (m*1588+1145)/2291  { Ln(X) = Ln(2)*Log2(X) }
  end if   

end function { XLn } 



{ -------------------------------------------------------------------------

    XLog2(X,lgx)   Extended Binary Logarithm Function (base 2) 

  

              'X' = Input can be any integer value

            'lgx' = Log2(X) scaled by 1,000,000  ie Log2(X) = 'lgx'/1000000



          NOTE: For X <= 0, result lgx is set to NAL = -100.000000 }
          

function XLog2(X,lgx)

  declare const NAL := -100000000  { Token for Log2(0) or negative X, -100.000000 } 

  declare c  { Log characteristic }

  declare m  { Log mantissa }

  

  if X <= 0

    lgx := NAL            { Log is not defined, return NAL }

  else

    XLg.Core(X,c,m)       { Get binary log of 1 <= X <= 2,147,483,647 }

    lgx := c*1000000 + m  { Log2(X) }
  end if   

end function { XLog2 }

{   XLog10(X,lgx)   Extended Common Logarithm Function (base 10)

  

             'X' = Input can be any integer value

           'lgx' = Log10(X) scaled by 1,000,000  ie Log10(X) = 'lgx'/1000000



          NOTE: For X <= 0, result lgx is set to NAL = -100.000000 }

                

function XLog10(X,lgx)

  declare const NAL := -100000000 { Token for Log10(0) or negative X, -100.000000 } 

  declare c  { Log characteristic }

  declare m  { Log mantissa }

  

  if X <= 0

    lgx := NAL             { Log is not defined, return NAL }

  else

    XLg.Core(X,c,m)        { Get binary log of 1 <= X <= 2,147,483,647 }

    lgx := (c*70139989+116)/233 + (m*1268/407*83+400)/859  { Log10(X) = Log10(2)*Log2(X) }
  end if   

end function { XLog10 } 

{ -------------------------------------------------------------------

   XSinCos(ang,sin,cos)    Extended-Angle Sine/Cosine Function 

  

    'ang' = Input angle in deci-grads (1000 deci-grads = 90 degrees)

            Angle may be any value of a 32-bit signed-integer

    'sin' = Scaled Sine of input angle: 'sin'/10000 = sin(ang)

    'cos' = Scaled Cosine of input angle: 'cos'/10000 = cos(ang) }



function XSinCos(ang,sin,cos)

  declare ra    { Reduced and reflected angle }

  declare quad  { Geometric Quadrant, I..IV }

  

  ra := ang     { Copy of caller's angle }

  ReduceAngle(ra,quad) { Reduce and reflect angle as needed }

  SinCos(ra,sin,cos)   { Compute sin & cos of reduced, Q1 angle, 'ra' }

  select quad          { Attach appropriate signs per quadrant }

    case 2

      cos := -cos

    case 3

      cos := -cos

      sin := -sin

    case 4

      sin := -sin

  end select  

end function { XSinCos }



{  -------------------------------------------------------------------

    XTangent(ang,tan)   Extended-Angle Tangent Function 

  

     'ang' = Input angle in deci-grads (1000 deci-grads = 90 degrees)

             Angle may be any value of a 32-bit signed-integer

     'tan' = Scaled Tangent of input angle: 'tan'/10000 = tan(ang)

    

    NOTE: The tangent function is discontinuous at +/-90 degrees, 

          that is: tan(1000) -> infinity. Since tan(999) is 6,250,000

          (625 scaled by 10K), for tan(1000) this routine returns an

          even larger value MaxTan = 100M as a sentinel for tan(1000) } 



function XTangent(ang,tan)

  declare ra   { Reduced angle }

  declare quad

  

  ra := ang            { Copy of caller's angle }

  ReduceAngle(ra,quad) { Reduce and reflect angle as needed }

  Tangent(ra,tan)      { Compute Tangent of reduced, Q1 angle, 'ra' }

  if quad = 2 or quad = 4  { Now attach correct sign for quadrant }

    tan := -tan            { Tangent is negative in quadrants 2 & 4 }

  end if

end function { XTangent }

{ ------------------------------------------------------------------------------------
                           Utility Conversion Routines
  ------------------------------------------------------------------------------------
  CVHex(num,str)   Convert integer number 'num' to an 8-digit Hexadecimal string, 'str' 



      'num' = Input 32-bit integer

      'str' = 8-digit Hex representation of input number 
      
      NOTE: In order to use this function in your script, you must invoke the routine
            named UseHex in your ICB before the first usage of CVHex. }

function CVHex(num,str)

  declare n

  declare x

  x := num    { Copy of num }

  str := ''   { Reset the output string to null }

  for n := 1 to 8   { Derive 8 digits from low to high }

    if x .and. 0xF < 10

      str := x .and. 0xF & str  { Add new digit 0..9 on left of string }

    else                        { Fetch hex digit A..F and add to left of string }

      str := HexDgt[(x .and. 0xF) - 10] & str

    end if

    x := sh_right(x,4)          { Shift next 4-bits into view }

  end for

end function { CVHex }

function UseHex { Call this function in the ICB if you want to use CVHex in your script }
  declare global !HexDgt[6]

  HexDgt[0] := 'A'

  HexDgt[1] := 'B'

  HexDgt[2] := 'C'

  HexDgt[3] := 'D'

  HexDgt[4] := 'E'

  HexDgt[5] := 'F'  
end function { UseHex }

{----------------------------------------------------------------------------------  
  ep_to_mdb(N,Vol)   Convert Engine Par 'N' to equivalent volume in mdb 



      'N' = Input engine parameter, 0..1,000,000

    'Vol' = Corresponding Volume in mdb



     This function uses Kontakt's 18db/octave relationship as follows:

     vol = 20000*log(N/Nm)^3 + 12000   where 'Vol' is in mdb and Nm = 1,000,000
     ie Vol = 60000*log(N) - 60000*log(Nm) + 12000 = 60000*log(N) - 348000
     and thus: Vol = 0.06*Log10(N) - 348000
     
     NOTE: If N = 0, Vol is set to Muted as a token for -infinity 
           If N > 1000000, Vol is clamped to +12,000 mdb }



function ep_to_mdb(N,Vol)
  DefineMuted              { Declare Muted constant }

  XLog10(N,Vol)            { Vol = 10000*log(N) }
  if Vol < 0
    Vol := Muted           { -Infinity = Muted}
  else if N >= 1000000
    Vol := 12000           { Max engine par volume is 12000 mdb }
  else

    Vol := (3*Vol + 25)/50 - 348000  { Vol = 0.06*Log10(N) - 348000 }
  end if
end function { ep_to_mdb }

{ ------------------------------------------------------------------------------  
  Get_db(VR,Vol)   Outputs volume in mdb equivalent to input V/V0 = VR/10000

  

   'VR' = Input volume ratio 0.0 < V/V0 <= 1.0000 as an integer scaled by 10,000

  'Vol' = Vol = 20000*log(V/V0). Thus for V/V0 = 0.0001, Vol = -80,000 mdb and for

          V/V0 >= 1.0000, Vol = 0 mdb. For V/V0 <= 0, log(V/V0) is not defined, so

          Vol is set to Muted. Since VR = 10000*(V/V0), then V/V0 = VR/10000 and:

          Vol = 20000[log(VR/10000)] = 20000[log(VR) - 4] = 20000*log(VR) - 80000

      

      NOTE: This routine uses the fast Binary Log2 function (DeVos Approximation) }

          

function Get_db(VR,Vol)

  DefineMuted       { Declare Muted constant }
  

  if VR >= 10000    { If 'VR' is too big }

    Vol := 0        { Clamp output to 0db reference }

  else

    Log10(VR,Vol)   { vol = Log10(VR) = 10000*log(VR) } 

    if Vol < 0

      Vol := Muted  { 'VR' <= 0, output Muted }

    else { 0.0001 < V/V0 < 1.0000 }

      Vol := 2*Vol - 80000 { Vol = 20000*log(VR) - 80000 }

    end if

  end if

end function { Get_db }

{ ------------------------------------------------------------------------------------

  mdb_to_ep(Vl,N)   Convert Volume in mdb to its equivalent ep value N



    'Vol' = Volume in mdb

      'N' = Equivalent engine parameter, 0..1,000,000

      

     This function uses Kontakt's 18db/octave relationship as follows:

     Vol = 20000*log(N/Nm)^3 + 12000   where 'Vol' is in mdb and Nm = 1,000,000

     Therefore: log(N/Nm) = (Vol - 12000)/60000 or N/Nm = 10^(Vol - 12000)/60000

     Thus: N = Nm*10^(Vol - 12000)/60000 = Nm*10^(Vol + 200000)/60000*10^K

           where K = -212000/60000

     and therefore:  N = 292.864*10^(Vol + 200000)/60000, or
       N = 292.864*Exp10[(1000000/60000)*(Vol + 200000)] 
         = 292.864*Exp10[(100/6)*(Vol + 200000)]

     

     NOTE: If Vol <= Muted, it is treated as -infinity and N is set to 0

           If Vol > 12,000 mdb, N is clamped to 1,000,000 }



function mdb_to_ep(Vol,N)
  DefineMuted

  declare vx

  

  if Vol <= Muted

    N := 0             { Treat as -infinity }

  else if Vol >= 12000

    N := 1000000       { Max engine par volume is 12000 mdb }

  else

    vx := 100*(Vol + 200003)/6 { Rounded Positive volume exponent, scaled by 10^6 } 

    Exp10(vx,N)                { N = 10^(vx/1000000) }

    N := (N*292864 + 500)/1000 { Multiply by 292.864 }

  end if

end function { mdb_to_ep } 

{ ------------------------------------------------------------------------------------

     VR_to_ep(VR,N)  Compute Engine Parameter for Input Volume Ratio 



     'VR' = Volume ratio scaled by 10,000   ie (V/V0)*10000

      'N' = Corresponding Engine Parameter (0..1,000,000) for 0.0 <= V/V0 <= 4.0000

     

     Using Kontakt's 18db/octave control taper, this routine computes the Engine

     Parameter, ep, as N = 10^6*(0.25*V/V0)^1/3, or N = 10^4*(25*VR)^1/3, which sets
     N = 629,961 (0db) when VR = 1.0000 and to N = 1,000,000 (+12db) when VR = 4.0000. 

     The most significant digits are calculated using Newton's Method and the least

     3 digits are then derived from the remainder using a novel approximation

     technique.

     

     NOTE: For VR > 4.0000, N is clamped at 1,000,000 (+12db) 
           For VR < 0.0, N is clamped at 0 }



function VR_to_ep(VR,N)

  declare SR   { Scaled Volume Ratio }

  declare n1   { Last iteration estimate }

  declare n2   { Prior to Last iteration estimate }



  if VR <= 0

    N := 0     { Cube root of zero is zero, we're done }

  else if VR > 40000

    N := 1000000   { VR > 4.0000, clamp N to 1,000,000 }

  else { 0.0 < VR < 4.0000 }

    SR := 25000*VR { SR = 0.25*(V/V0)*10^9 }

    n1 := 30       { Use n1 to find top bit position of SR }

    while sh_left(SR,31 - n1) >= 0

      dec(n1)

    end while

    N := sh_left(1,(n1+2)/3) { Set N to first estimate for cube root of SR }

    n1 := -1         { Impossible prior estimates to force the }

    n2 := -2         {  iteration loop to run at least once    }

    while N # n1 and N # n2  { Loop until N matches either last or prior estimate }

      n2 := n1       { Update prior estimate }

      n1 := N        { Update last estimate }

      N := (2*n1 + SR/(n1*n1))/3 { Compute new estimate from last estimate }

    end while

    if n1 < N        { If last two iterations differ, return the smaller, }

      N := n1        {  this only happens when iterations oscillate }

    end if           { We now have the integer part of the root in N }

    n1 := SR - N*N*N { Set n1 to remainder and use it to approximate the root fraction }

    N := 1000*N + 333*n1/(N*N) { Add fractional part of root and scale by 1000 }

  end if

end function { VR_to_ep } 

{ ------------------------------------------------------------------------------

  VR_to_mdb(VR,Vol)   Converts scaled input volume ratio to mdb equivalent

  

   'VR' = Input volume ratio 0.0 < V/V0 <= 4.0000 as an integer scaled by 10,000

  'Vol' = Vol = 20000*log(V/V0). Thus for V/V0 = 0.0001, Vol = -80,000 mdb and for

          V/V0 = 1.0000, Vol = 0 mdb. For V/V0 = 4.0000, Vol = +12,000 mdb. And,
          For V/V0 <= 0, log(V/V0) is not defined, so Vol is set to Muted. Since
          VR = 10000*(V/V0), then V/V0 = VR/10000 and:

          Vol = 20000[log(VR/10000)] = 20000[log(VR) - 4] = 20000*log(VR) - 80000 
          
          NOTE: This routine performs the same function as Get_db but is more
                accurate and can accept volume ratios up to 4.0000. However, the
                legacy Get_db routine executes faster since it is based on the
                DeVos approximation for Log2. }

          

function VR_to_mdb(VR,Vol)

  DefineMuted       { Declare Muted constant }

  

  if VR >= 40000    { If 'VR' is too big }

    Vol := 12000    { Clamp output to +12db above zero reference }

  else

    XLog10(VR,Vol)  { Vol = Log10(VR) = 1000000*log(VR) } 

    if Vol < 0

      Vol := Muted  { 'VR' <= 0, output Muted }

    else { 0.0001 < V/Vm < 4.0000 }

      Vol := (Vol + 25)/50 - 80000  { Vol = 20000*log(VR) - 80000 }

    end if

  end if

end function { VR_to_mdb }

{ ------------------------------------------------------------------
                         Library Support Routines 
  ------------------------------------------------------------------ 
  BuildAngleTable    Data Function to declare the arctangent table, 
                     AngTbl, its scaling constant AngScale, and the
                     right-angle constant Ang90 for Trig functions.
                     As coded angles are assumed to be in deci-grads
                     and are scaled by 10^6.

  NOTE: To change to a different angular unit or resolution,
        only this Data Function need be edited appropriately,
        see Technical Guide }



function BuildAngleTable
  declare global const TrigBits := 15   { Precision in bits (#of iterations ) }
  declare global const AngScale := 1000000  { Input Angle scaling factor }
  declare global AngTbl[TrigBits+1] :=  { Angle Table (in dg scaled by 10^6) } ...

    { AngTbl[n] = arctan(2^-n) and AngTbl[0] is not used } (500000000, ...

     295167235, 155958261, 79166848,  39737049, 19887896, ...

       9946375,   4973491,   2486783,  1243396,   621699, ...

        310849,    155425,     77712,    38856,    19428)  

  declare global const Ang90 := 1000 { Right angle (90 degrees) = 1000 dg } 

end function { BuildAngleTable }

{ ------------------------------------------------------------------

  BuildLogTable      Data Function to declare the Logs table used by the
                     XP2.Core and XLg.Core functions. Currently the Logs
                     array elements are scaled by 10^9. 



  NOTE: The Logs array elements contain the value Logs[n] = Lg(1+2^-n) for each
        n from 1 to LogBits, Logs[0] is not used. Both XLg.Core and XP2.Core
        use a rounding 'constant' named LSB2 which is set to Logs[LogBits]/2.
        If LogBits is changed, the value of this rounding constant must also be
        changed. Therefore, it is defined here to keep it with he Logs table. }

function BuildLogTable

  declare global const LogBits := 20  { Precision bits (or #of loop iterations) }

  declare global Logs[LogBits+1] :=   { Logs Table Lg(1+2^-n), 0 < Logs[n] < 1 } ... 
       (1000000000, 584962501, 321928095, 169925001, 87462841, 44394119, ... 
        22367813, 11227256, 5624549, 2815016, 1408195, 704269, 352178, ...
        176099, 88052, 44028, 22014, 11006, 5504, 2751, 1376)
  declare global const LSB2 := 688    { Rounding Bias, Logs[LogBits]/2 }
end function { BuildLogTable } 

{ ------------------------------------------------------------------ 

  DefineMuted     Data function to declare 'Muted' constant

                  needed by ATFade and Get_db }



function DefineMuted

  declare global const Muted := -200000 { -200 db in mdb }

end function { DefineMuted }

{ ------------------------------------------------------------------  

  ReduceAngle(ang,quad)    Reduces and reflects the input angle

  

        'ang' = Input angle over full signed integer range. Returns

                a 1st quadrant angle from 0..1000 dg (0..90 deg)

       'quad' = Returns the geometric quadrant number from 1..4 

       

       NOTE: Quadrature angles are still considered to be in the
             'approach' quadrant. For example (in degrees), +90d is
             in Q1 while -270d is in Q2. Similarly, +180d is in Q2
             while -180 is in Q3, etc }



function ReduceAngle(ang,quad)

  declare const Ang360 := 4*Ang90 { Full circle, 360 degrees }

  declare ra                      { Reduced and reflected angle }

  

  ra := abs(ang) mod Ang360   { ra in range 0 <= ra < 4000 dg }

  quad := (ra - 1)/Ang90 + 1  { Quadrant of reduced, positive angle }

  ra := ra - Ang90*(quad - 1) { Quadrant-excess angle }

  if quad = 2 or quad = 4

    ra := Ang90 - ra   { Use complementary angle in Q2 & Q4 }

  end if

  if ang < 0

    quad := 5 - quad   { Complement quadrant for negative 'ang' }

  end if

  ang := ra  { Return reduced angle }  

end function { ReduceAngle } 

{ -------------------------------------------------------------------------

   XLg.Core(X,c,m)  Support Function for Extended Binary Logs  

  

            'X' = Input, Positive Integer in the range 0 < X < 2^31

            'c' = Output, Charcteristic of Log2(X)

            'm' = Output, Manitissa of Log2(X) scaled by 10^6

                    

          This routine uses a CORDIC-Like algorithm that utilizes the same

          Logs array as that used by the exponential support routine XP2.Core. }

          

function XLg.Core(X,c,m)

  BuildLogTable    { Declare LogBits (precision) and Logs array }

  declare const Bit30 := 0x40000000  { Value of highest positive bit position }   

  declare n    { Iteration index }

  declare nx   { Normalized X, with binary point between bits 29 & 30 }

  declare sx   { Shifted nx }

  

  nx := X      { Copy of input X to be normalized }

  c := 30      { X = nx*2^30 }

  while nx < Bit30       { Normalize 1.0 <= nx < 2.0      }

    nx := sh_left(nx,1)  {  and adjust exponent as needed }

    dec(c)

  end while

  if nx = Bit30   { nx = 1.0 }

    m := 0        { X is an integral power of 2 }

  else { Mantissa must be calculated for 1.0 < nx < 2.0 }

    m := 1000000000 - LSB2 { m = 1.0 - LSB/2 for rounding, scaled by 10^9 }

    n := 1

    while n <= LogBits and nx >= 0  { More bits to derive and nx not exactly 2.0 yet }

      { Cordic loop to drive nx -> 2.0 as m -> log2(nx) }

      sx := sh_right(nx,n) { Next increment to make nx -> 2.0 } 
      if (-nx - sx) < 0    { Is sx + nx <= 2.0, ie can we add sx without exceeding 2.0? }

        nx := nx + sx      { Yes, increase nx by sx }

        m := m - Logs[n]   { Continue to make m -> Log2(nx) }

      end if

      inc(n)

    end while

    m := (m + 500)/1000    { Downscale m from 10^9 to 10^6 } 

  end if

end function { XLg.Core }

{ -------------------------------------------------------------------------

   XP2.Core(Z,X)  Support Function for Binary AntiLogs  

  

              'Z' = Input, essentially a base 2 logarithm (scaled by 10^6)

              'X' = Output, 2^(Z/1000000) ie the Base 2 Antilog of 'Z'

                    where 'Z' is in the range from 0.0 <= Z <= 30.9999

                    

          This routine uses a CORDIC-Like algorithm that utilizes the same

          Logs array as that used by the log support routine XLg.Core. }


function XP2.Core(Z,X)

  BuildLogTable { Declare LogBits (precision) and Logs array }

  declare const Bit30 := 0x40000000  { Value of highest positive bit position }

  declare C     { Integer part of Z, ie the characteristic }

  declare M     { Fractonal part of Z, ie the mantissa }

  declare n     { Iteration index }

  

  if Z < 0
    X := 0      { 2^-Z = 0 in integer domain }
  else  { 0.0 <= Z < 31.000000 }
    C := Z/1000000        { Characteristic of Z }

    M := Z mod 1000000    { Mantissa of Z, scaled by 10^6 }

    X := Bit30            { Use X as temp for SX = 1.0 scaled by 2^30 }

    if M # 0              { Skip mantissa iteration loop for even powers of two }

      M := M*1000 + LSB2  { Rescale M to match Logs table (10^9) and bias by LSB/2 }

      n := 1

      while n <= LogBits  { Execute Cordic Loop to drive M -> 0.0 and SX -> 2^M }

        if M >= Logs[n]

          M := M - Logs[n]

          X := X + sh_right(X,n)  { Accumulating scaled product for SX = 2^M }

        end if

        inc(n)

      end while  { SX now = 2^M * 2^30 }

    end if 

    { If C = 30, X = SX; else X = SX*2^(C-30) }
    if C < 30 { Right shift SX by (30-C) with rounding }

      X := sh_right(sh_right(X,29-C)+1,1) .and. 0x7FFFFFFF  

    end if
  end if

end function { XP2.Core } 
