[libre-riscv-dev] [Bug 127] Transcendentals needed (SIN/COS/ATAN2/EXP/LOG/POW etc.)

[bugzilla-daemon at libre-riscv.org]

http://bugs.libre-riscv.org/show_bug.cgi?id=127

--- Comment #5 from Luke Kenneth Casson Leighton <lkcl at lkcl.net> ---
https://www.quinapalus.com/efunc.html

Example C code

Here is a sample C function to compute exp() using the above algorithm. The
code assumes integers are at least 32 bits long. The (non-negative) argument
and the result of the function are both expressed as fixed-point values with 16
fractional bits. Notice that after 11 steps of the algorithm the constants
involved become such that the code is simply doing a multiplication: this is
explained in the note below.

The extension to negative arguments is left as an exercise.

  int fxexp(int x) {
  int t,y;

  y=0x00010000;
  t=x-0x58b91; if(t>=0) x=t,y<<=8;
  t=x-0x2c5c8; if(t>=0) x=t,y<<=4;
  t=x-0x162e4; if(t>=0) x=t,y<<=2;
  t=x-0x0b172; if(t>=0) x=t,y<<=1;
  t=x-0x067cd; if(t>=0) x=t,y+=y>>1;
  t=x-0x03920; if(t>=0) x=t,y+=y>>2;
  t=x-0x01e27; if(t>=0) x=t,y+=y>>3;
  t=x-0x00f85; if(t>=0) x=t,y+=y>>4;
  t=x-0x007e1; if(t>=0) x=t,y+=y>>5;
  t=x-0x003f8; if(t>=0) x=t,y+=y>>6;
  t=x-0x001fe; if(t>=0) x=t,y+=y>>7;
  if(x&0x100)               y+=y>>8;
  if(x&0x080)               y+=y>>9;
  if(x&0x040)               y+=y>>10;
  if(x&0x020)               y+=y>>11;
  if(x&0x010)               y+=y>>12;
  if(x&0x008)               y+=y>>13;
  if(x&0x004)               y+=y>>14;
  if(x&0x002)               y+=y>>15;
  if(x&0x001)               y+=y>>16;
  return y;
  }

Note that the constants involved in the trial subtractions reduce by a factor
of less than or equal to 2 at each step. This means that it is never necessary
to test the same constant twice.
A note on the residual error

The error in the final answer depends on the residual value in x; in fact, the
relative error is exp(x). Since for small x, exp(x) is approximately 1+x, it is
possible to correct the final answer by multiplying it by 1+x. 

Example C code

Here is a sample C function to compute log() using the above algorithm. The
code assumes integers are at least 32 bits long. The (positive) argument and
the result of the function are both expressed as fixed-point values with 16
fractional bits, although intermediates are kept with 31 bits of precision to
avoid loss of accuracy during shifts. After 12 steps of the algorithm the
correction described above is applied.

  int fxlog(int x) {
  int t,y;

  y=0xa65af;
  if(x<0x00008000) x<<=16,              y-=0xb1721;
  if(x<0x00800000) x<<= 8,              y-=0x58b91;
  if(x<0x08000000) x<<= 4,              y-=0x2c5c8;
  if(x<0x20000000) x<<= 2,              y-=0x162e4;
  if(x<0x40000000) x<<= 1,              y-=0x0b172;
  t=x+(x>>1); if((t&0x80000000)==0) x=t,y-=0x067cd;
  t=x+(x>>2); if((t&0x80000000)==0) x=t,y-=0x03920;
  t=x+(x>>3); if((t&0x80000000)==0) x=t,y-=0x01e27;
  t=x+(x>>4); if((t&0x80000000)==0) x=t,y-=0x00f85;
  t=x+(x>>5); if((t&0x80000000)==0) x=t,y-=0x007e1;
  t=x+(x>>6); if((t&0x80000000)==0) x=t,y-=0x003f8;
  t=x+(x>>7); if((t&0x80000000)==0) x=t,y-=0x001fe;
  x=0x80000000-x;
  y-=x>>15;
  return y;
  }

-- 
You are receiving this mail because:
You are on the CC list for the bug.