# Square root



## matthiasA (Jun 29, 2012)

This is probably a very stupid question, but can somebody help me calculating a simple square root with a script?


----------



## Big Bob (Jun 29, 2012)

Hi Matthias,

There is a cube root algorithm in my math library but I never had any occasion to do a square root. If I get a little time later I'll see what I can cook up for you but in the meantime, you may want to look at the cube root implementation discussed in the technical guide.

https://dl.dropbox.com/u/80404485/Kontakt/Math%20Library/Docs/KSPMathTechGuide_V215.pdf (https://dl.dropbox.com/u/80404485/Konta ... e_V215.pdf)

If I remember correctly, I did the cube root function with a combination of a Newton iteration and a Haley. However, I would think that for a simple square root that Newton's method would converge rapidly so you might want to just try that. 

Anyway, as soon as I get some fun-stuff time I'll give it a shot for you if no one has responded in the meantime. But hopefully someone else has already done this and if so, I'm sure they will chime in :lol: 

Rejoice,

Bob

BTW What situation did you run into that requires extracting a square root?


----------



## matthiasA (Jun 29, 2012)

Hi Bob,

my problem with the Newton iteration is to place the first guess. For a cubical convergence the Halley's rational formula works perfectly, but I'm yet to find a way how to translate this to square roots.

Doing iterations starting from 0 take much too long to calculate, obviously.

By the way, I hope I am not basically mistaking when I assume that the formula for making the intensity sliders linear (instead of exponential) is
*
ep= ep^0,5 * 1000 ?*


----------



## Big Bob (Jun 29, 2012)

> Doing iterations starting from 0 take much too long to calculate, obviously.



You could always use a small lookup table to get a better initial guess.



> By the way, I hope I am not basically mistaking when I assume that the formula for making the intensity sliders linear (instead of exponential) is
> 
> ep= ep^0,5 * 1000 ?



Are you talking about the Intensity sliders associated with modulation?

BTW You could always resort to using the relationship that:


```
log[x^0.5] = 0.5*log(x)
```
 which can be computed with the existing math library routines (log and exp).


----------



## matthiasA (Jun 29, 2012)

I have found some algorithm that sounds promising, since it deals with integers. Unfortunately, I'm neither a programmer nor a math genius, so I greatly failed to translate it into KSP. Here it is:

x = 2^ceil(numbits(N)/2)
loop:
y = floor((x + floor(N/x))/2)
if y >= x
return x
x = y

Any ideas?


----------



## matthiasA (Jun 29, 2012)

Big Bob @ Fri Jun 29 said:


> Are you talking about the Intensity sliders associated with modulation?



Yes. If you dial in 0 they stay at 0. If you dial in 500000, they go to 25%. If you dial in 707000, they go to 50%. If you dial in 1000000, they go to 100%. This matches the formula y=x^0,5*1000 perfectly, doesn't it?


----------



## Big Bob (Jun 29, 2012)

Yes indeed it do :lol: 

I'm sorry to sort of drop out on you right in the middle of this but it seems my heart has decided to start acting up again so I guess I'll have to shut down for a while (hopefully just for a while :roll: )

I there some reason you couldn't cobble this together with the log/exp routines? If execution time seems like it would be a problem you can use V406 of the Math Library and use the Fast mode.

In case you don't have it already, here's a link to the V406 package.

https://dl.dropbox.com/u/80404485/Kontakt/Math%20Library/V406/MathLibraryV406.zip (https://dl.dropbox.com/u/80404485/Konta ... ryV406.zip)

Anyway, I'll be offline for a while (boy I'll sure be glad when I get my new Glorified body) o=< 

Rejoice,

Bob


----------



## matthiasA (Jun 30, 2012)

The math 406 and its NKA tables are already heavily in use  But in this case I'm too stupid to make the calculations out of log/exp. 

No problem. I already wrote a fast routine based on Newton and a lookup table as you suggested. It's not very accurate, but it it's really fast and works well in terms of sound.

I keep all my fingers crossed for you! My plan is to bring you some fresh apple pie once my instrument is released (and sells at least a little bit so I can afford to visit the US), so please recover! 

God bless you,
Matthias


----------



## matthiasA (Jun 30, 2012)

Bob,

thanks for your hints, I got it. Since {log(b) of n_root(x)} = {1/n * log(b) of x} we simply have to calculate the ep with LOG(ep;2)/2 and then perform 2^result.

You're the greatest !


----------



## Big Bob (Jul 2, 2012)

Hi Matthias,

Here's the square-root routine I promised you.

*on init*
``make_perfview
``set_script_title('SqRoot Demo')
``*declare* ui_value_edit X (0,1000000,1)
````move_control(X,1,1)
``*declare* ui_label Show (1,1)
````move_control(Show,2,1)
``*declare* root
``message('')
*end* on

*on ui_control*(X)
``root := SqRootx1000(X)
``set_text(Show,root)
*end* on

_{ Computes the integer square root of X, where X is any
positive integer. The result returned is the highest
integer below the correct root (ie the floor value) }_

*function* ISqRoot(X) -> result _{ X >= 0 }_
``*declare* n``_{ next newton iteration }_
``*if* X < 2```_{ For X = 1 or 0, result is X itself }_
````result := X
``*else* _{ X > 1 }_
````n := num_lead_zeros(X-1)
````n := 16 - n/2`````````````_{ use lowest power of 2 above }_
````result := sh_left(1,n)````_{ the root as initial guess }_
````n := sh_right(result + sh_right(X,n),1) _{ 1st iteration }_
````*while* n < result```_{ repeat Newton iterations as needed }_
``````result := n _{ loop body will execute a max of 5 times }_
``````n := (result + X/result)/2
````*end while*
``*end* *if*
*end* *function* _{ ISqRoot }_

*function* SqRootx1000(X) -> result _{ Scaled square root }_
``result := ISqRoot(X)```_{ Integer root }_
``result := 1000*result + 500*(X-result*result)/result
``_{ scale by 1000 and add approximate fractional part }_
*end* *function* _{ SqRoot }_

_{ Obtains the number of leading zero bits
in X, where X is any positive number }_
*function* num_lead_zeros(X) -> result
``result := 1 + sh_right(((0x7FFF8000.*and*. X) - 1),31)*(-16)
``result := result + sh_right(((sh_right(0xFF000000,result) .*and*. X) - 1),31)*(-8)
``result := result + sh_right(((sh_right(0xF0000000,result) .*and*. X) - 1),31)*(-4)
``result := result + sh_right(((sh_right(0xC0000000,result) .*and*. X) - 1),31)*(-2)
``result := result + sh_right(((sh_right(0x80000000,result) .*and*. X) - 1),31)*(-1) 
*end* *function* _{ num_lead_zeros }_



*ISqRoot(X)* computes the floor integer root for any positive X value. However, since you plan to mulitply the root by 1000 anyway, you might want to include a fractional part also. So, I've included *SqRootx1000(X)* which outputs the square root of X scaled by 1000. The fractional part of the root is fairly accurate as long as X > 100 or so. But in any case, it will be more accurate than just the floor integer root by itself. I'm attaching a simple demo instrument so you can play with this a little. Set the Edit box to a value between 0 and 1000000 and the label will display 1000*X^0.5.

*ISqRoot* executes in about 2.3 us and *SqRootx1000* executes in about 2.7 us. If this isn't fast enough, there are a few things can be done to speed it up but they will bloat the code. For example, the *num_lead_zeros *routine can be made to run about twice as fast as the 5-line routine I've posted but, it requires a 32-way select case construct (and therefore about 66 code lines).

Sorry I had to bail out on you last week but I have a heart condition that often acts up these days. So far, the Good Lord has always raised me back up from these episodes and thankfully this time was no exception. So, praise the Lord!

Please let me know if you need anything else to make this thing work.

Rejoice,

Bob

BTW Why do you want to linearize the Intensity slider (since when its at say 50% of its traversal it's still going to display 25%)? I'm wondering what application this might have and whether or not I should include such a converter in the next Math Library update?



> I keep all my fingers crossed for you! My plan is to bring you some fresh apple pie once my instrument is released (and sells at least a little bit so I can afford to visit the US), so please recover!



I'll be looking forward to tasting some of that fresh apple pie. :lol:


----------



## Big Bob (Jul 10, 2012)

Hey Mathias,

Can you tell me with which modules the Intensity sliders behave this way? For most things like the Source module or Amplifier module I think the Intensity sliders track a cubic function not quadratic. For example, I just finished helping Paul linearize the Intensity slider for the source module and it tracks the usual Kontakt cubic funtion:

percent = 100*[(ep - 500000)/500000]^3

Apparently you've found somewhere where the Intensity sliders track with a second-order curve? I'd like to know where that is :roll: 

Rejoice,

Bob

BTW: The reason I'm asking about this is that I'm now in the process of updating the Math Library so, if I knew it might be useful for something, I might include a Root2 function.


----------



## matthiasA (Jul 14, 2012)

Unless I'm mistaking the percent slider for constants is a square root. I needed it for keyboard tracking. Unfortunately I dont have any computer handy except my iPhone, so I cannot check that before monday evening.

But I hope (and guess) you will be right. That would be great since the LN/EXP combo is quite CPU hungry...


----------



## Big Bob (Jul 14, 2012)

> Unless I'm mistaking the percent slider for constants is a square root



Not generally. Everywhere I've assigned a constant modulator, the Intensity slider takes on a cubic curvature. For example, see this thread:

http://www.vi-control.net/forum/viewtop ... f32d9e84db

So, can you give me an example of where you assign a constant modulator so it follows a quadratic curvature?



> That would be great since the LN/EXP combo is quite CPU hungry...


That shouldn't be the case if you include F1+F2 in the SetMathMode option list.

Rejoice,

Bob


----------



## matthiasA (Jul 15, 2012)

I was having a hard time to install KONTAKT in an internet cafe here, and my wife will probably kill me for spending my time in front of a computer during holidays... but here is the result for an constant intensity slider controlling cutoff.

'square' is the calculated value=(EP/1000000)^2
'get_disp' is the measured result and 
'delta' is the deviation. 


```
EP	square	get_disp	delta
50000	0.25	0.25	0
100000	1.00	1,00	0
125000	1.56	1.60	0.03
250000	6.25	6.30	0.04
500000	25.00	25.00	0
750000	56.25	56.30	0.04
800000	64.00	64.00	0
900000	81.00	81.00	0
950000	90.25	90.30	0.04
1000000	100.00	100.00	0
```

I dunno how to achieve similar results with a cubic function, especially because the result matches the KONTAKT algorithm by 100%... right?


----------



## Big Bob (Jul 15, 2012)

I don't know how you are obtaining the data you just posted :? 

I just inserted a LP1 legacy filter in the first Group Insert Fx slot in K4. Then I assigned a Constant modulator to it which K4 gives the name CV_CUTOFF. Then I loaded this very simple script.

*on init* 
``message('') 
``*declare* ui_value_edit ep (0,1000000,1)
``*declare* ui_label Show (1,1)
*end* on 

*on ui_control*(ep)
``set_engine_par(ENGINE_PAR_INTMOD_INTENSITY,ep,0,find_mod(0,'CV_CUTOFF'),-1) 
``set_text(Show,get_engine_par_disp(ENGINE_PAR_INTMOD_INTENSITY,0,find_mod(0,'CV_CUTOFF'),-1)) 
*end* on

The data I obtained (which agrees with Paul's data posted in the SS thread) looks like this:


```
ep       disp
   0     -100%
100K    -51.2%
200K    -21.6%
300K     -6.4%
400K     -0.8%
500K        0%
600K      0.8%
700K      6.4%
800K     21.6%
900K     51.2%
  1M      100%
```

This clearly follows the cubic function:


```
disp = 100*[(ep - 5000000)/500000]^3
```

Perhaps (when you get back from your vacation, because I wouldn't want to get you into trouble with your wife), you could post the exact details along with the script you are using to obtain the data you are seeing. 

Believe me, I don't doubt that you are getting the data that you are reporting. I just don't know how you are doing it. :roll: 

Rejoice,

Bob

BTW I notice that the forum display of my script is inserting some spaces in the CV_CUTOFF string. Those will have to be removed if you try the posted script.


----------



## mk282 (Jul 15, 2012)

Bob, use the CODE tag, it shouldn't add arbitrary spaces then.


----------



## Big Bob (Jul 15, 2012)

Hi Mario
Usually, I can't get the code tags to work in conjuction with the KSE syntax highlighting, which generally I prefer. But, I guess when it messes up like it did, I should probably settle for just code tags :lol: 

BTW Mario, do you know of anywhere in Kontakt where a Constant modulator exhibits the quadratic characteristic that Matthias is seeing? So far, I can't seem to repro his numbers but I probably haven't hit the right place. :roll:


----------



## matthiasA (Jul 15, 2012)

The code is almost identical.


```
on init
declare ui_slider $value (0,1000000)
end on

on ui_control ($value)
set_engine_par($ENGINE_PAR_MOD_TARGET_INTENSITY,$value,0,18,-1)
message($value&" "&get_engine_par_disp($ENGINE_PAR_MOD_TARGET_INTENSITY,0,18,-1))
end on
```

Since I never used K4 I dunno how it has been handled there, but for K5 it seems that negative values will be clipped. Can't check it at the moment, this Kontakt engine 5.03 I installed here in the cafe is ... well... you know... and behaves quite strange (it won't accept on ui_control!). 


PS: Well, definitely back on tuesday  It's better like that, I love her.

However, to subtract constant intensities from any modulation target in K5, I need to adress the "invert" button. Not?


----------



## Big Bob (Jul 15, 2012)

Aha! :idea: At last you have disclosed the key piece of information that was missing. 8) 

You are using ENGINE_PAR_MOD_TARGET_INTENSITY instead of ENGINE_PAR_INTMOD_INTENSITY. 

It has been known for some time that the former is uni-polar while the latter is bipolar. But, what has not been exposed previously is that the non-linearty of the former is apparently quadratic while the latter is cubic. This is a very important fact that has previously been overlooked, at least by me.

As to K5 versus K4, the slider control behaves the same as in K4 (when using INTMOD) but there *is* a difference in that the get_ep_disp function doesn't display any values below 1.0% (which occurs at about ep = 607722) whereas K4 displays all the values correctly. *NOTE MK's later post because this has been corrected with K5.03*

As to the Invert button situation, when using INTMOD (vs using MOD_TARGET), for ep values below 500000, the Invert button activates automatically whereas for MOD_TARGET, you have to turn it on with a separate ep. 

There are then two ways to handle this. The first is use INTMOD and then you can use the standard cubic linearization described in following thread (of course in K5 ep_disp itself will not report the correct values for ep < 607722 but your application probably doesn't need to use ep_disp once you have the control linearized, since you can just display the control variable itself). Again this goes away with K5.0.3

http://www.vi-control.net/forum/viewtop ... f32d9e84db

The 2nd way to do it would be to use MOD_TARGET and write the needed quadratic linearizing routine (as you have been trying to do).

As soon as I finish this post, I'll examine this situation more thoroughly now that I know how you have been getting this data :D . 

V450 of the Math Library will contain among other things, a cubic linearizer (and possibly now a quadratic linearizer) for modulation intensity. The library will also contain a general Power function that computes: X^(a/b) with a specifiable number of decimal digits (depending on the needed accuracy). So, for example, you will be able to use Power(X,1,2,3) to compute the square root of X with 3 decimal digits in the result. In addition, I will probably add a separate Root2 routine.

Rejoice,

Bob


----------



## Big Bob (Jul 15, 2012)

OK Matthias,

Here's a quick and dirty quadratic linearizer that you can use with MOD_TARGET (it uses the same SqRootx1000 that I previously posted).

*on init*```_{ quadratic linearizer }_
``message('') 
``*declare* ui_knob Int (0,1000,10)
````set_knob_unit(Int,KNOB_UNIT_PERCENT)
``*declare* ui_label Show (1,1)
``*declare* ep
*end* on 

*on ui_control*(Int)
``ep := SqRootx1000(10*Int)
``ep := 10*ep
``set_engine_par(ENGINE_PAR_MOD_TARGET_INTENSITY,ep,0,find_mod(0,'CV_CUTOFF'),-1) 
``set_text(Show,get_engine_par_disp(ENGINE_PAR_MOD_TARGET_INTENSITY,0,find_mod(0,'CV_CUTOFF'),-1)) 
*end* on

*function* ISqRoot(X) -> result _{ X >= 0 }_
``*declare* n``_{ next newton iteration }_
``*if* X < 2```_{ For X = 1 or 0, result is X itself }_
````result := X
``*else* _{ X > 1 }_
````n := num_lead_zeros(X-1)
````n := 16 - n/2`````````````_{ use lowest power of 2 above }_
````result := sh_left(1,n)````_{ the root as initial guess }_
````n := sh_right(result + sh_right(X,n),1) _{ 1st iteration }_
````*while* n < result```_{ repeat Newton iterations as needed }_
``````result := n _{ loop body will execute a max of 5 times }_
``````n := (result + X/result)/2
````*end while*
``*end* *if*
*end* *function* _{ ISqRoot }_

*function* SqRootx1000(X) -> result _{ Scaled square root }_
``result := ISqRoot(X)```_{ Integer root }_
``result := 1000*result + 500*(X-result*result)/result
``_{ scale by 1000 and add approximate fractional part }_
*end* *function* _{ SqRoot }_

_{ Obtains the number of leading zero bits
in X, where X is any positive number }_
*function* num_lead_zeros(X) -> result
``result := 1 + sh_right(((0x7FFF8000.*and*. X) - 1),31)*(-16)
``result := result + sh_right(((sh_right(0xFF000000,result) .*and*. X) - 1),31)*(-8 )
``result := result + sh_right(((sh_right(0xF0000000,result) .*and*. X) - 1),31)*(-4)
``result := result + sh_right(((sh_right(0xC0000000,result) .*and*. X) - 1),31)*(-2)
``result := result + sh_right(((sh_right(0x80000000,result) .*and*. X) - 1),31)*(-1) 
*end* *function* _{ num_lead_zeros }_

V450 of the Math Library will probably have a more refined version of this but meanwhile this might bridge the gap. I've also attached a demo instrument for you to play with.

Rejoice,

Bob

NOTE: Please watch out again for added spaces. And, Mario, this time I tried posting it with Code tags and it still added a space in CV_CUTOFF but this way at least we can benefit from the syntax highlighting 8)


----------



## mk282 (Jul 16, 2012)

Big Bob @ 15.7.2012 said:


> As to K5 versus K4, the slider control behaves the same as in K4 (when using INTMOD) but there *is* a difference in that the get_ep_disp function doesn't display any values below 1.0% (which occurs at about ep = 607722) whereas K4 displays all the values correctly.



This was FIXED in Kontakt 5.0.3.


----------



## Big Bob (Jul 16, 2012)

Thanks MK, I guess I haven't got the latest version then. I better hop to it :lol:


----------



## matthiasA (Jul 18, 2012)

Thanks a lot Bob! 

I didn't implement the cube/square function anymore, since in the actual state it isn't faster than the combined log/exp function. But let's see how this develops once it's based on FastMath in the upcoming v450 

Have a great time,
Matthias


----------



## Big Bob (Jul 18, 2012)

Just to give you a heads up, here's some typical, relative execution times for the pertinent routines in V450.

Power(X,1,2,3) 14.7 usec in Standard Mode, 2.2 usec in Fast Mode (F1+F2)
Root2(X) 1.8 usec in Standard Mode

Power(X,1,3,5) 15.1 us in Standard Mode, 2.3 us in Fast Mode (F1+F2)
Root3(X) 2.1 usec in Standard Mode

Power(X,1,2,3) performs exactly the same function as Root2, namely it computes X^1/2 scaled by 1000. Power(X,1,3,5) performs exactly the same function as Root3, it computes X^1/3 scaled by 100000.

So, if you are already using both FLog.nka and FALog.nka, you can probably just use the general *Power* function because its execution time is only a tad slower. However, for applications that otherwise don't need to use the F1 or F2 modes,* Root2 *and *Root3* will be a better choice. Therefore, I have implemented the new *epModInt* routine using *Root2 *and *Root3* rather than *Power*. *epModInt *linearizes modulation intensity for either INTMOD or MOD_TARGET by choosing the desired ep type as one of its input arguments.

BTW. I have dropped the Fast mode for Root3 in V450 since it was only marginally faster and less accurate than the Standard Mode.

Rejoice,

Bob


----------



## matthiasA (Jul 21, 2012)

Me want, me want !!


----------

