Math is hard

Hey y'all --

So this is one of those times that my lack of serious math chops is
coming round to bite me, and none of my usual books is helping me out.
I'm hoping someone has some thoughts.

I'm trying to approximate either exp(-1/(n+1)) or 4^(-1/n+1).  I can
convince myself I don't care which.  n is an integer from 1-65535, and
the result should be fixed-point fractional, probably U0.18.  The
function output is always between 0-1, and goes up like a rocket for
small n before leveling off to a steady cruise of >0.9 for the rest of
the function domain.

I'm working in an FPGA, so I've got adds, multiplies, and table lookups
from tables of reasonable size (10s of kb) cheap, but other things
(divides especially) are expensive.  I can throw several clock cycles
at the problem if need be.

Taylor series attacks seem to fail horribly.  I feel like there may be
some answer where answers for n in [1,127] gets a direct table lookup,
and n in [128,65535] gets some other algorithm, possibly with a table
boost.  Or somehow taking advantage of the fact that log(1-f(n)) is
related to log(n)?

Anyone have any thoughts?

-- 
Rob Gaddi, Highland Technology -- www.highlandtechnology.com
Email address domain is currently out of order.  See above to fix.H

In comp.arch.fpga Rob Gaddi <rgaddi@technologyhighland.invalid> wrote:

(snip)
> I'm trying to approximate either exp(-1/(n+1)) or 4^(-1/n+1).  I can
> convince myself I don't care which.  n is an integer from 1-65535, and
> the result should be fixed-point fractional, probably U0.18.  The
> function output is always between 0-1, and goes up like a rocket for
> small n before leveling off to a steady cruise of >0.9 for the rest of
> the function domain.

You want an 18 bit function of a 16 bit input. Fits into one BRAM
on most current, and not so current, FPGAs.  Generate all 65536
values and store them in BRAM (ROM in verilog/VHDL terms).

-- glen

On Fri, 17 Jan 2014 11:00:06 -0800, Rob Gaddi wrote:

> Hey y'all --
> 
> So this is one of those times that my lack of serious math chops is
> coming round to bite me, and none of my usual books is helping me out.
> I'm hoping someone has some thoughts.
> 
> I'm trying to approximate either exp(-1/(n+1)) or 4^(-1/n+1).  I can
> convince myself I don't care which.  n is an integer from 1-65535, and
> the result should be fixed-point fractional, probably U0.18.  The
> function output is always between 0-1, and goes up like a rocket for
> small n before leveling off to a steady cruise of >0.9 for the rest of
> the function domain.
> 
> I'm working in an FPGA, so I've got adds, multiplies, and table lookups
> from tables of reasonable size (10s of kb) cheap, but other things
> (divides especially) are expensive.  I can throw several clock cycles at
> the problem if need be.
> 
> Taylor series attacks seem to fail horribly.  I feel like there may be
> some answer where answers for n in [1,127] gets a direct table lookup,
> and n in [128,65535] gets some other algorithm, possibly with a table
> boost.  Or somehow taking advantage of the fact that log(1-f(n)) is
> related to log(n)?
> 
> Anyone have any thoughts?

First, instead of Taylor's series, try a best-fit polynomial.  I suspect 
you'll be nearly as disappointed, but you can try.

Second, when best-fit polynomials let me down, I look to interpolation.  
Depending on the problem this ends up being somewhere between linear and 
cubic.

You should be aware that for any one way of going from your tables of 
values to an interpolation, there's way more than one way of coming up 
with tables of values, some of which are better than others.

Depending on how motivated you are you can either generate splines, or  
do a least-squares best fit, or find an overall minimax solution 
(minimize the maximum error).  Note that this doesn't change the 
algorithm in your FPGA -- it just changes the values in your look-up 
tables.  Formulating the problem as linear interpolation and using the 
known points, or formulating the problem as cubic splines and looking up 
the algorithm to generate your tables is probably the easiest if you've 
got the typical engineer's "jack of all trades" way with math.

If you want to get fancier, then once you have the problem formulated, 
least-squares can be done with linear algebra on something like Scilab or 
Matlab in less time than it takes to lift your finger off the enter key.  
Getting a minimax solution requires more work both in formulating the 
problem and in waiting for the optimizer to finish.  It also means you 
have to have a decent optimization software package (which Scilab comes 
with: I'm not sure if it's native to Matlab or if you have to buy yet 
another package).

-- 

Tim Wescott
Wescott Design Services
http://www.wescottdesign.com

On 17/01/14 20:00, Rob Gaddi wrote:
> Hey y'all --
>
> So this is one of those times that my lack of serious math chops is
> coming round to bite me, and none of my usual books is helping me out.
> I'm hoping someone has some thoughts.
>
> I'm trying to approximate either exp(-1/(n+1)) or 4^(-1/n+1).  I can
> convince myself I don't care which.  n is an integer from 1-65535, and
> the result should be fixed-point fractional, probably U0.18.  The
> function output is always between 0-1, and goes up like a rocket for
> small n before leveling off to a steady cruise of >0.9 for the rest of
> the function domain.
>
> I'm working in an FPGA, so I've got adds, multiplies, and table lookups
> from tables of reasonable size (10s of kb) cheap, but other things
> (divides especially) are expensive.  I can throw several clock cycles
> at the problem if need be.
>
> Taylor series attacks seem to fail horribly.  I feel like there may be
> some answer where answers for n in [1,127] gets a direct table lookup,
> and n in [128,65535] gets some other algorithm, possibly with a table
> boost.  Or somehow taking advantage of the fact that log(1-f(n)) is
> related to log(n)?
>
> Anyone have any thoughts?
>

Taylor series are great for the theory, but very rarely a useful 
calculation technique in practice.

In a case like this, I would first look at using a single table - if you 
can afford the space, use it for simplicity.

Failing that, draw a graph and look at the function as it changes.  Also 
draw some graphs of the derivatives.  Any areas of the function where 
you have low second derivatives are suitable for linear interpolation. 
Where you have high second derivatives, you need either tightly packed 
linear interpolation or cubic splines of some sort.  There are a few 
different kinds of cubic splines, depending on things like smooth joins 
between the parts, least maximal errors, least total square error, etc. 
  But for a limited range function like this, it is not unrealistic to 
get bit-perfect results if that is what you need.

So once you have these graphs, you can look at linear interpolation with 
a fixed step size, linear interpolation with varying step size, cubic 
splines of some sort, or dividing the range up and using different 
techniques for different areas.

On Saturday, January 18, 2014 8:00:06 AM UTC+13, Rob Gaddi wrote:
> Hey y'all --
> 
> 
> 
> So this is one of those times that my lack of serious math chops is
> 
> coming round to bite me, and none of my usual books is helping me out.
> 
> I'm hoping someone has some thoughts.
> 
> 
> 
> I'm trying to approximate either exp(-1/(n+1)) or 4^(-1/n+1).  I can
> 
> convince myself I don't care which.  n is an integer from 1-65535, and
> 
> the result should be fixed-point fractional, probably U0.18.  The
> 
> function output is always between 0-1, and goes up like a rocket for
> 
> small n before leveling off to a steady cruise of >0.9 for the rest of
> 
> the function domain.
> 
> 
> 
> I'm working in an FPGA, so I've got adds, multiplies, and table lookups
> 
> from tables of reasonable size (10s of kb) cheap, but other things
> 
> (divides especially) are expensive.  I can throw several clock cycles
> 
> at the problem if need be.
> 
> 
> 
> Taylor series attacks seem to fail horribly.  I feel like there may be
> 
> some answer where answers for n in [1,127] gets a direct table lookup,
> 
> and n in [128,65535] gets some other algorithm, possibly with a table
> 
> boost.  Or somehow taking advantage of the fact that log(1-f(n)) is
> 
> related to log(n)?
> 
> 
> 
> Anyone have any thoughts?
> 
> 
> 
> -- 
> 
> Rob Gaddi, Highland Technology -- www.highlandtechnology.com
> 
> Email address domain is currently out of order.  See above to fix.H

exp(x)=1+x+x^2/2+x^3/3!+...

On Friday, January 17, 2014 2:00:06 PM UTC-5, Rob Gaddi wrote:
> Hey y'all --
> 
> 
> 
> So this is one of those times that my lack of serious math chops is
> 
> coming round to bite me, and none of my usual books is helping me out.
> 
> I'm hoping someone has some thoughts.
> 
> 
> 
> I'm trying to approximate either exp(-1/(n+1)) or 4^(-1/n+1).  I can
> 
> convince myself I don't care which.  n is an integer from 1-65535, and
> 
> the result should be fixed-point fractional, probably U0.18.  The
> 
> function output is always between 0-1, and goes up like a rocket for
> 
> small n before leveling off to a steady cruise of >0.9 for the rest of
> 
> the function domain.
> 
> 
> 
> I'm working in an FPGA, so I've got adds, multiplies, and table lookups
> 
> from tables of reasonable size (10s of kb) cheap, but other things
> 
> (divides especially) are expensive.  I can throw several clock cycles
> 
> at the problem if need be.
> 
> 
> 
> Taylor series attacks seem to fail horribly.  I feel like there may be
> 
> some answer where answers for n in [1,127] gets a direct table lookup,
> 
> and n in [128,65535] gets some other algorithm, possibly with a table
> 
> boost.  Or somehow taking advantage of the fact that log(1-f(n)) is
> 
> related to log(n)?
> 
> 
> 
> Anyone have any thoughts?
> 
> 
> 
> -- 
> 
> Rob Gaddi, Highland Technology -- www.highlandtechnology.com
> 
> Email address domain is currently out of order.  See above to fix.H

Rob, do you have a misplaced parenthesis(ses)? One case -1/(n+1), the other is -1/n+1

Case one's limit is 1 whereas case 2's limit is 4 for large n.

You can scale the exponents such as to change the radix with the effect of change an nth root problem to a power problem. Does that help?

For example let B=16 and N in the domain of [1,16]

4^(1-1/N) = 4 * (4^(1/B)) ^ P = 4*(C^P) with C = 4^(1/B) and P = B/N

Then 4^(1/B) is a precomputed constant and P=B/N is in [16,1]

Food for thought!

Clay

David Brown <david.brown@hesbynett.no> writes:

> On 17/01/14 20:00, Rob Gaddi wrote:
>> Hey y'all --
>>
>> So this is one of those times that my lack of serious math chops is
>> coming round to bite me, and none of my usual books is helping me out.
>> I'm hoping someone has some thoughts.
>>
>> I'm trying to approximate either exp(-1/(n+1)) or 4^(-1/n+1).  I can
>> convince myself I don't care which.  n is an integer from 1-65535, and
>> the result should be fixed-point fractional, probably U0.18.  The
>> function output is always between 0-1, and goes up like a rocket for
>> small n before leveling off to a steady cruise of >0.9 for the rest of
>> the function domain.
>>
>> I'm working in an FPGA, so I've got adds, multiplies, and table lookups
>> from tables of reasonable size (10s of kb) cheap, but other things
>> (divides especially) are expensive.  I can throw several clock cycles
>> at the problem if need be.
>>
>> Taylor series attacks seem to fail horribly.  I feel like there may be
>> some answer where answers for n in [1,127] gets a direct table lookup,
>> and n in [128,65535] gets some other algorithm, possibly with a table
>> boost.  Or somehow taking advantage of the fact that log(1-f(n)) is
>> related to log(n)?
>>
>> Anyone have any thoughts?
>>
>
> Taylor series are great for the theory, but very rarely a useful
> calculation technique in practice.
>
> In a case like this, I would first look at using a single table - if
> you can afford the space, use it for simplicity.

Or if a "full" table is too big, consider reducing the table size
and interpolating between table entries.
-- 
Randy Yates
Digital Signal Labs
http://www.digitalsignallabs.com

In comp.arch.fpga Randy Yates <yates@digitalsignallabs.com> wrote:
> David Brown <david.brown@hesbynett.no> writes:

(snip)
>> In a case like this, I would first look at using a single table -
>> if you can afford the space, use it for simplicity.
 
> Or if a "full" table is too big, consider reducing the table size
> and interpolating between table entries.

Using a table of interpolation values and an adder.

-- glen

glen herrmannsfeldt wrote:
> In comp.arch.fpga Rob Gaddi <rgaddi@technologyhighland.invalid> wrote:
> 
> (snip)
>> I'm trying to approximate either exp(-1/(n+1)) or 4^(-1/n+1).  I can
>> convince myself I don't care which.  n is an integer from 1-65535, and
>> the result should be fixed-point fractional, probably U0.18.  The
>> function output is always between 0-1, and goes up like a rocket for
>> small n before leveling off to a steady cruise of >0.9 for the rest of
>> the function domain.
> 
> You want an 18 bit function of a 16 bit input. Fits into one BRAM
> on most current, and not so current, FPGAs.  Generate all 65536
> values and store them in BRAM (ROM in verilog/VHDL terms).
> 
> -- glen

Glen,

   What devices are you using.  My measly Xilinx parts only have
36K bits at most per BRAM.  I'd need 36 of those to do this table.

-- 
Gabor

In comp.arch.fpga GaborSzakacs <gabor@alacron.com> wrote:

(snip, I wrote)

>> You want an 18 bit function of a 16 bit input. Fits into one BRAM
>> on most current, and not so current, FPGAs.  Generate all 65536
>> values and store them in BRAM (ROM in verilog/VHDL terms).

(snip)
>   What devices are you using.  My measly Xilinx parts only have
> 36K bits at most per BRAM.  I'd need 36 of those to do this table.

Last time I was working with Xilinx, I only needed tiny RAMs.

The BRAMs have 18 inputs and 18 outputs (which I remembered)
but it seems to be dual port, two 9 bit addresses (the part
I didn't remember).

The Spartan-6 devices have between 12 and 268 such BRAMs, though.

How many such BRAM can you use? Easiest way to do interpolation
is with look-up tables for the interpolation values. 

Start with a 2Kx18 RAM, so 11 address bits. A second 2Kx18 RAM
is addresses by the high bits (six in this case) of the address,
and the rest (five bits) from the 16 bit address. That, plus
an 18 bit adder, possibly pipelined, will give you linear
interpolation on the 2Kx18 values.

If that isn't enough, two interpolation RAMS, indexed by
combinations of the high and low address bits, will do better.
One way to look at the interpolation ROM is that the high bits
select a slope, and the low bits select an appropriate multiple
of that slope to be added. 

That is what was done when MOS ROMs were about that size.
You can also use the block multipliers, with interpolation
tables as inputs. 

A lot depends on how fast it needs to be, and how many of them
you need at the same time.

-- glen