Good fixed point math libraries?

Put your problem here if it does not fit any of the other categories.

Good fixed point math libraries?

Postby AGr » Thu Jun 18, 2009 3:56 pm

Could anyone recommend a good fixed point math library to use?
AGr
Freshman
Freshman
 
Posts: 3
Joined: Fri Jun 05, 2009 7:28 am
Location: TX

Top

Postby AGr » Sat Jun 20, 2009 7:35 am

I found one that I have been testing out, the original code can be found at http://sourceforge.net/project/showfile ... p_id=24849, I renamed and removed a couple things in the posting here

It is fairly basic, but has all the trig functions, I did some initial testing to make sure it is accurate and I got sin/cos to be within 0.001% and atan to be within 0.05% at worst over the input range. I haven't done any timing but no doubt they are orders of magnitude faster than the java.lang.Math equivalent

Syntax: [ Download ] [ Hide ]
Using java Syntax Highlighting
  1.  
  2. /**
  3.  
  4. *
  5.  
  6. * 16:16 fixed point math routines, for IAppli/CLDC platform.
  7.  
  8. * A fixed point number is a 32 bit int containing 16 bits of integer and 16 bits of fraction.
  9.  
  10. *<p>
  11.  
  12. * (C) 2001 Beartronics
  13.  
  14. * Author: Henry Minsky (hqm@alum.mit.edu)
  15.  
  16. *<p>
  17.  
  18. * Licensed under terms "Artistic License"<br>
  19.  
  20. * <a href="http://www.opensource.org/licenses/artistic-license.html">http://www.opensource.org/licenses/artistic-license.html</a><br>
  21.  
  22. *
  23.  
  24. *<p>
  25.  
  26. * Numerical algorithms based on
  27.  
  28. * http://www.cs.clemson.edu/html_docs/SUN ... s.doc.html
  29.  
  30. * <p>
  31.  
  32. * Trig routines based on numerical algorithms described in
  33.  
  34. * http://www.magic-software.com/MgcNumerics.html
  35.  
  36. *
  37.  
  38. * http://www.dattalo.com/technical/theory/logs.html
  39.  
  40. *
  41.  
  42. * @version $Id: FP.java,v 1.6 2001/04/05 07:40:17 hqm Exp $
  43.  
  44. */
  45.  
  46.  
  47.  
  48. public class FP {
  49.  
  50.  
  51.  
  52.         public static int toInt(int x) {
  53.  
  54.                 return x >> 16;
  55.  
  56.         }
  57.  
  58.  
  59.  
  60.         public static int fromInt(int x) {
  61.  
  62.                 return x << 16;
  63.  
  64.         }
  65.  
  66.        
  67.  
  68.         public static int fromFloat(float x) {
  69.  
  70.                 return (int)(x * (1 << 16));
  71.  
  72.         }
  73.  
  74.        
  75.  
  76.         public static int fromDouble(double x) {
  77.  
  78.                 return (int)(x * (1 << 16));
  79.  
  80.         }
  81.  
  82.        
  83.  
  84.         public static double toDouble(int x) {
  85.  
  86.                 return (double)x / (1 << 16);
  87.  
  88.         }
  89.  
  90.        
  91.  
  92.         public static float toFloat(int x) {
  93.  
  94.                 return (float)x / (1 << 16);
  95.  
  96.         }
  97.  
  98.  
  99.  
  100.        
  101.  
  102.         /** Multiply two fixed-point numbers */
  103.  
  104.         public static int mul(int x, int y) {
  105.  
  106.                 long z = (long) x * (long) y;
  107.  
  108.                 return ((int) (z >> 16));
  109.  
  110.         }
  111.  
  112.  
  113.  
  114.         /** Divides two fixed-point numbers */
  115.  
  116.         public static int div(int x, int y) {
  117.  
  118.                 long z = (((long) x) << 32);
  119.  
  120.                 return (int) ((z / y) >> 16);
  121.  
  122.         }
  123.  
  124.  
  125.  
  126.         /** Compute square-root of a 16:16 fixed point number */
  127.  
  128.         public static int sqrt(int n) {
  129.  
  130.                 int s = (n + 65536) >> 1;
  131.  
  132.                 for (int i = 0; i < 8; i++) {
  133.  
  134.                         s = (s + div(n, s)) >> 1;
  135.  
  136.                 }
  137.  
  138.                 return s;
  139.  
  140.         }
  141.  
  142.  
  143.  
  144.         /** Round to nearest fixed point integer */
  145.  
  146.         public static int round(int n) {
  147.  
  148.                 if (n > 0) {
  149.  
  150.                         if ((n & 0x8000) != 0) {
  151.  
  152.                                 return (((n + 0x10000) >> 16) << 16);
  153.  
  154.                         } else {
  155.  
  156.                                 return (((n) >> 16) << 16);
  157.  
  158.                         }
  159.  
  160.                 } else {
  161.  
  162.                         int k;
  163.  
  164.                         n = -n;
  165.  
  166.                         if ((n & 0x8000) != 0) {
  167.  
  168.                                 k = (((n + 0x10000) >> 16) << 16);
  169.  
  170.                         } else {
  171.  
  172.                                 k = (((n) >> 16) << 16);
  173.  
  174.                         }
  175.  
  176.                         return -k;
  177.  
  178.                 }
  179.  
  180.         }
  181.  
  182.  
  183.  
  184.         public static final int PI = 205887;
  185.  
  186.         public static final int PI_OVER_2 = PI / 2;
  187.  
  188.         public static final int E = 178145;
  189.  
  190.         public static final int HALF = 2 << 15;
  191.  
  192.        
  193.  
  194.        
  195.  
  196.         static final int SK1 = 498;
  197.  
  198.         static final int SK2 = 10882;
  199.  
  200.  
  201.  
  202.         /**
  203.  
  204.          * Computes SIN(f), f is a fixed point number in radians. 0 <= f <= 2PI
  205.  
  206.          */
  207.  
  208.         public static int sin(int f) {
  209.  
  210.                 // If in range -pi/4 to pi/4: nothing needs to be done.
  211.  
  212.                 // otherwise, we need to get f into that range and account for
  213.  
  214.                 // sign change.
  215.  
  216.  
  217.  
  218.                 int sign = 1;
  219.  
  220.                 if ((f > PI_OVER_2) && (f <= PI)) {
  221.  
  222.                         f = PI - f;
  223.  
  224.                 } else if ((f > PI) && (f <= (PI + PI_OVER_2))) {
  225.  
  226.                         f = f - PI;
  227.  
  228.                         sign = -1;
  229.  
  230.                 } else if (f > (PI + PI_OVER_2)) {
  231.  
  232.                         f = (PI << 1) - f;
  233.  
  234.                         sign = -1;
  235.  
  236.                 }
  237.  
  238.  
  239.  
  240.                 int sqr = mul(f, f);
  241.  
  242.                 int result = SK1;
  243.  
  244.                 result = mul(result, sqr);
  245.  
  246.                 result -= SK2;
  247.  
  248.                 result = mul(result, sqr);
  249.  
  250.                 result += (1 << 16);
  251.  
  252.                 result = mul(result, f);
  253.  
  254.                 return sign * result;
  255.  
  256.         }
  257.  
  258.  
  259.  
  260.         static final int CK1 = 2328;
  261.  
  262.         static final int CK2 = 32551;
  263.  
  264.  
  265.  
  266.         /**
  267.  
  268.          * Computes cos(f), f is a fixed point number in radians. 0 <= f <= PI/2
  269.  
  270.          */
  271.  
  272.         public static int cos(int f) {
  273.  
  274.  
  275.  
  276.                 int sign = 1;
  277.  
  278.                 if ((f > PI_OVER_2) && (f <= PI)) {
  279.  
  280.                         f = PI - f;
  281.  
  282.                         sign = -1;
  283.  
  284.                 } else if ((f > PI_OVER_2) && (f <= (PI + PI_OVER_2))) {
  285.  
  286.                         f = f - PI;
  287.  
  288.                         sign = -1;
  289.  
  290.                 } else if (f > (PI + PI_OVER_2)) {
  291.  
  292.                         f = (PI << 1) - f;
  293.  
  294.                 }
  295.  
  296.  
  297.  
  298.                 int sqr = mul(f, f);
  299.  
  300.                 int result = CK1;
  301.  
  302.                 result = mul(result, sqr);
  303.  
  304.                 result -= CK2;
  305.  
  306.                 result = mul(result, sqr);
  307.  
  308.                 result += (1 << 16);
  309.  
  310.                 return result * sign;
  311.  
  312.         }
  313.  
  314.  
  315.  
  316.         /**
  317.  
  318.          * Computes tan(f), f is a fixed point number in radians. 0 <= f <= PI/4
  319.  
  320.          */
  321.  
  322.  
  323.  
  324.         static final int TK1 = 13323;
  325.  
  326.         static final int TK2 = 20810;
  327.  
  328.        
  329.  
  330.         public static int tan(int f) {
  331.  
  332.                 int sqr = mul(f, f);
  333.  
  334.                 int result = TK1;
  335.  
  336.                 result = mul(result, sqr);
  337.  
  338.                 result += TK2;
  339.  
  340.                 result = mul(result, sqr);
  341.  
  342.                 result += (1 << 16);
  343.  
  344.                 result = mul(result, f);
  345.  
  346.                 return result;
  347.  
  348.         }
  349.  
  350.  
  351.  
  352.         /**
  353.  
  354.          * Computes atan(f), f is a fixed point number |f| <= 1
  355.  
  356.          * <p>
  357.  
  358.          * For the inverse tangent calls, all approximations are valid for |t| <= 1.
  359.  
  360.          * To compute ATAN(t) for t > 1, use ATAN(t) = PI/2 - ATAN(1/t). For t < -1,
  361.  
  362.          * use ATAN(t) = -PI/2 - ATAN(1/t).
  363.  
  364.          */
  365.  
  366.         public static int atan(int f) {
  367.  
  368.                 int sqr = mul(f, f);
  369.  
  370.                 int result = 1365;
  371.  
  372.                 result = mul(result, sqr);
  373.  
  374.                 result -= 5579;
  375.  
  376.                 result = mul(result, sqr);
  377.  
  378.                 result += 11805;
  379.  
  380.                 result = mul(result, sqr);
  381.  
  382.                 result -= 21646;
  383.  
  384.                 result = mul(result, sqr);
  385.  
  386.                 result += 65527;
  387.  
  388.                 result = mul(result, f);
  389.  
  390.                 return result;
  391.  
  392.         }
  393.  
  394.  
  395.  
  396.         static final int AS1 = -1228;
  397.  
  398.         static final int AS2 = 4866;
  399.  
  400.         static final int AS3 = 13901;
  401.  
  402.         static final int AS4 = 102939;
  403.  
  404.  
  405.  
  406.         /**
  407.  
  408.          * Compute asin(f), 0 <= f <= 1
  409.  
  410.          */
  411.  
  412.  
  413.  
  414.         public static int asin(int f) {
  415.  
  416.                 int fRoot = sqrt((1 << 16) - f);
  417.  
  418.                 int result = AS1;
  419.  
  420.                 result = mul(result, f);
  421.  
  422.                 result += AS2;
  423.  
  424.                 result = mul(result, f);
  425.  
  426.                 result -= AS3;
  427.  
  428.                 result = mul(result, f);
  429.  
  430.                 result += AS4;
  431.  
  432.                 result = PI_OVER_2 - (mul(fRoot, result));
  433.  
  434.                 return result;
  435.  
  436.         }
  437.  
  438.  
  439.  
  440.         /**
  441.  
  442.          * Compute acos(f), 0 <= f <= 1
  443.  
  444.          */
  445.  
  446.         public static int acos(int f) {
  447.  
  448.                 int fRoot = sqrt((1 << 16) - f);
  449.  
  450.                 int result = AS1;
  451.  
  452.                 result = mul(result, f);
  453.  
  454.                 result += AS2;
  455.  
  456.                 result = mul(result, f);
  457.  
  458.                 result -= AS3;
  459.  
  460.                 result = mul(result, f);
  461.  
  462.                 result += AS4;
  463.  
  464.                 result = mul(fRoot, result);
  465.  
  466.                 return result;
  467.  
  468.         }
  469.  
  470.  
  471.  
  472.         /**
  473.  
  474.          * Exponential
  475.  
  476.          * /** Logarithms:
  477.  
  478.          *
  479.  
  480.          * (2) Knuth, Donald E., "The Art of Computer Programming Vol 1",
  481.  
  482.          * Addison-Wesley Publishing Company, ISBN 0-201-03822-6 ( this comes from
  483.  
  484.          * Knuth (2), section 1.2.3, exercise 25).
  485.  
  486.          *
  487.  
  488.          * http://www.dattalo.com/technical/theory/logs.html
  489.  
  490.          *
  491.  
  492.          */
  493.  
  494.  
  495.  
  496.         /**
  497.  
  498.          * This table is created using base of e.
  499.  
  500.          *
  501.  
  502.          * (defun fixedpoint (z) (round (* z (lsh 1 16))))
  503.  
  504.          *
  505.  
  506.          * (loop for k from 0 to 16 do (setq z (log (+ 1 (expt 2.0 (- (+ k 1))))))
  507.  
  508.          * (insert (format "%d\n" (fixedpoint z))))
  509.  
  510.          */
  511.  
  512.         static int log2arr[] = { 26573, 14624, 7719, 3973, 2017, 1016, 510, 256,
  513.  
  514.                         128, 64, 32, 16, 8, 4, 2, 1, 0, 0, 0 };
  515.  
  516.  
  517.  
  518.         static int lnscale[] = { 0, 45426, 90852, 136278, 181704, 227130, 272557,
  519.  
  520.                         317983, 363409, 408835, 454261, 499687, 545113, 590539, 635965,
  521.  
  522.                         681391, 726817 };
  523.  
  524.  
  525.  
  526.         public static int ln(int x) {
  527.  
  528.                 // prescale so x is between 1 and 2
  529.  
  530.                 int shift = 0;
  531.  
  532.  
  533.  
  534.                 while (x > 1 << 17) {
  535.  
  536.                         shift++;
  537.  
  538.                         x >>= 1;
  539.  
  540.                 }
  541.  
  542.  
  543.  
  544.                 int g = 0;
  545.  
  546.                 int d = HALF;
  547.  
  548.                 for (int i = 1; i < 16; i++) {
  549.  
  550.                         if (x > ((1 << 16) + d)) {
  551.  
  552.                                 x = div(x, ((1 << 16) + d));
  553.  
  554.                                 g += log2arr[i - 1]; // log2arr[i-1] = log2(1+d);
  555.  
  556.                         }
  557.  
  558.                         d >>= 1;
  559.  
  560.                 }
  561.  
  562.                 return g + lnscale[shift];
  563.  
  564.         }      
  565.  
  566. }
  567.  
  568.  
Parsed in 0.056 seconds, using GeSHi 1.0.8.4
AGr
Freshman
Freshman
 
Posts: 3
Joined: Fri Jun 05, 2009 7:28 am
Location: TX

Postby ninor » Wed Jul 01, 2009 11:08 pm

Thanks a lot!
Image AndDev: Your Android Development Community / Tutorials | Here's my Basic ToolKit
User avatar
ninor
Moderator
Moderator
 
Posts: 180
Joined: Thu Aug 14, 2008 6:30 pm
Location: Barcelona, Spain

Trig function range bug

Postby Curran » Sun Jan 03, 2010 10:09 pm

Thanks, this is just what I needed and seems to work very well. However, the trig functions break with input outside 0 to 2pi. Inserting the following code at the beginning of the trig functions will fix this:

if (f < 0) f = -f+PI;
f=f%TWOPI;
Curran
Once Poster
Once Poster
 
Posts: 1
Joined: Sun Jan 03, 2010 10:03 pm
Location: Lowell,MA

Top

Return to Other Coding-Problems

Who is online

Users browsing this forum: No registered users and 16 guests