gnc-numeric.c

00001 /********************************************************************
00002  * gnc-numeric.c -- an exact-number library for accounting use      *
00003  * Copyright (C) 2000 Bill Gribble                                  *
00004  * Copyright (C) 2004 Linas Vepstas <linas@linas.org>               *
00005  *                                                                  *
00006  * This program is free software; you can redistribute it and/or    *
00007  * modify it under the terms of the GNU General Public License as   *
00008  * published by the Free Software Foundation; either version 2 of   *
00009  * the License, or (at your option) any later version.              *
00010  *                                                                  *
00011  * This program is distributed in the hope that it will be useful,  *
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of   *
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the    *
00014  * GNU General Public License for more details.                     *
00015  *                                                                  *
00016  * You should have received a copy of the GNU General Public License*
00017  * along with this program; if not, contact:                        *
00018  *                                                                  *
00019  * Free Software Foundation           Voice:  +1-617-542-5942       *
00020  * 51 Franklin Street, Fifth Floor    Fax:    +1-617-542-2652       *
00021  * Boston, MA  02110-1301,  USA       gnu@gnu.org                   *
00022  *                                                                  *
00023  *******************************************************************/
00024 
00025 #include "config.h"
00026 
00027 #include <glib.h>
00028 #include <math.h>
00029 #include <stdio.h>
00030 #include <stdlib.h>
00031 #include <string.h>
00032 
00033 #include "gnc-numeric.h"
00034 #include "qofmath128.c"
00035 
00036 /* static short module = MOD_ENGINE; */
00037 
00038 /* =============================================================== */
00039 
00040 #if 0
00041 static const char *_numeric_error_strings[] = {
00042     "No error",
00043     "Argument is not a valid number",
00044     "Intermediate result overflow",
00045     "Argument denominators differ in GNC_HOW_DENOM_FIXED operation",
00046     "Remainder part in GNC_HOW_RND_NEVER operation"
00047 };
00048 #endif
00049 
00050 /* =============================================================== */
00051 /* This function is small, simple, and used everywhere below, 
00052  * lets try to inline it.
00053  */
00054 inline GNCNumericErrorCode
00055 gnc_numeric_check (gnc_numeric in)
00056 {
00057     if (in.denom != 0)
00058     {
00059         return GNC_ERROR_OK;
00060     }
00061     else if (in.num)
00062     {
00063         if ((0 < in.num) || (-4 > in.num))
00064         {
00065             in.num = (gint64) GNC_ERROR_OVERFLOW;
00066         }
00067         return (GNCNumericErrorCode) in.num;
00068     }
00069     else
00070     {
00071         return GNC_ERROR_ARG;
00072     }
00073 }
00074 
00075 /*
00076  *  Find the least common multiple of the denominators of a and b.
00077  */
00078 
00079 static inline gint64
00080 gnc_numeric_lcd (gnc_numeric a, gnc_numeric b)
00081 {
00082     qofint128 lcm;
00083     if (gnc_numeric_check (a) || gnc_numeric_check (b))
00084     {
00085         return GNC_ERROR_ARG;
00086     }
00087 
00088     if (b.denom == a.denom)
00089         return a.denom;
00090 
00091     /* Special case: smaller divides smoothly into larger */
00092     if ((b.denom < a.denom) && ((a.denom % b.denom) == 0))
00093     {
00094         return a.denom;
00095     }
00096     if ((a.denom < b.denom) && ((b.denom % a.denom) == 0))
00097     {
00098         return b.denom;
00099     }
00100 
00101     lcm = lcm128 (a.denom, b.denom);
00102     if (lcm.isbig)
00103         return GNC_ERROR_ARG;
00104     return lcm.lo;
00105 }
00106 
00107 
00108 /* Return the ratio n/d reduced so that there are no common factors. */
00109 static inline gnc_numeric
00110 reduce128 (qofint128 n, gint64 d)
00111 {
00112     gint64 t;
00113     gint64 num;
00114     gint64 denom;
00115     gnc_numeric out;
00116     qofint128 red;
00117 
00118     t = rem128 (n, d);
00119     num = d;
00120     denom = t;
00121 
00122     /* The strategy is to use Euclid's algorithm */
00123     while (denom > 0)
00124     {
00125         t = num % denom;
00126         num = denom;
00127         denom = t;
00128     }
00129     /* num now holds the GCD (Greatest Common Divisor) */
00130 
00131     red = div128 (n, num);
00132     if (red.isbig)
00133     {
00134         return gnc_numeric_error (GNC_ERROR_OVERFLOW);
00135     }
00136     out.num = red.lo;
00137     if (red.isneg)
00138         out.num = -out.num;
00139     out.denom = d / num;
00140     return out;
00141 }
00142 
00143 /* *******************************************************************
00144  *  gnc_numeric_zero_p
00145  ********************************************************************/
00146 
00147 gboolean
00148 gnc_numeric_zero_p (gnc_numeric a)
00149 {
00150     if (gnc_numeric_check (a))
00151     {
00152         return 0;
00153     }
00154     else
00155     {
00156         if ((a.num == 0) && (a.denom != 0))
00157         {
00158             return 1;
00159         }
00160         else
00161         {
00162             return 0;
00163         }
00164     }
00165 }
00166 
00167 /* *******************************************************************
00168  *  gnc_numeric_negative_p
00169  ********************************************************************/
00170 
00171 gboolean
00172 gnc_numeric_negative_p (gnc_numeric a)
00173 {
00174     if (gnc_numeric_check (a))
00175     {
00176         return 0;
00177     }
00178     else
00179     {
00180         if ((a.num < 0) && (a.denom != 0))
00181         {
00182             return 1;
00183         }
00184         else
00185         {
00186             return 0;
00187         }
00188     }
00189 }
00190 
00191 /* *******************************************************************
00192  *  gnc_numeric_positive_p
00193  ********************************************************************/
00194 
00195 gboolean
00196 gnc_numeric_positive_p (gnc_numeric a)
00197 {
00198     if (gnc_numeric_check (a))
00199     {
00200         return 0;
00201     }
00202     else
00203     {
00204         if ((a.num > 0) && (a.denom != 0))
00205         {
00206             return 1;
00207         }
00208         else
00209         {
00210             return 0;
00211         }
00212     }
00213 }
00214 
00215 /* *******************************************************************
00216  *  gnc_numeric_compare
00217  *  returns 1 if a>b, -1 if b>a, 0 if a == b 
00218  ********************************************************************/
00219 
00220 int
00221 gnc_numeric_compare (gnc_numeric a, gnc_numeric b)
00222 {
00223     gint64 aa, bb;
00224     qofint128 l, r;
00225 
00226     if (gnc_numeric_check (a) || gnc_numeric_check (b))
00227     {
00228         return 0;
00229     }
00230 
00231     if (a.denom == b.denom)
00232     {
00233         if (a.num == b.num)
00234             return 0;
00235         if (a.num > b.num)
00236             return 1;
00237         return -1;
00238     }
00239 
00240     if ((a.denom > 0) && (b.denom > 0))
00241     {
00242         /* Avoid overflows using 128-bit intermediate math */
00243         l = mult128 (a.num, b.denom);
00244         r = mult128 (b.num, a.denom);
00245         return cmp128 (l, r);
00246     }
00247 
00248     if (a.denom < 0)
00249         a.denom *= -1;
00250     if (b.denom < 0)
00251         b.denom *= -1;
00252 
00253     /* BUG: Possible overflow here..  Also, doesn't properly deal with
00254      * reciprocal denominators.
00255      */
00256     aa = a.num * a.denom;
00257     bb = b.num * b.denom;
00258 
00259     if (aa == bb)
00260         return 0;
00261     if (aa > bb)
00262         return 1;
00263     return -1;
00264 }
00265 
00266 
00267 /* *******************************************************************
00268  *  gnc_numeric_eq
00269  ********************************************************************/
00270 
00271 gboolean
00272 gnc_numeric_eq (gnc_numeric a, gnc_numeric b)
00273 {
00274     return ((a.num == b.num) && (a.denom == b.denom));
00275 }
00276 
00277 
00278 /* *******************************************************************
00279  *  gnc_numeric_equal
00280  ********************************************************************/
00281 
00282 gboolean
00283 gnc_numeric_equal (gnc_numeric a, gnc_numeric b)
00284 {
00285     qofint128 l, r;
00286     if ((a.denom == b.denom) && (a.denom > 0))
00287     {
00288         return (a.num == b.num);
00289     }
00290     if ((a.denom > 0) && (b.denom > 0))
00291     {
00292         // return (a.num*b.denom == b.num*a.denom);
00293         l = mult128 (a.num, b.denom);
00294         r = mult128 (b.num, a.denom);
00295         return equal128 (l, r);
00296 
00297 #if ALT_WAY_OF_CHECKING_EQUALITY
00298         gnc_numeric ra = gnc_numeric_reduce (a);
00299         gnc_numeric rb = gnc_numeric_reduce (b);
00300         if (ra.denom != rb.denom)
00301             return 0;
00302         if (ra.num != rb.num)
00303             return 0;
00304         return 1;
00305 #endif
00306     }
00307     if ((a.denom < 0) && (b.denom < 0))
00308     {
00309         l = mult128 (a.num, -a.denom);
00310         r = mult128 (b.num, -b.denom);
00311         return equal128 (l, r);
00312     }
00313     else
00314     {
00315         /* BUG: One of the numbers has a reciprocal denom, and the other
00316            does not. I just don't know to handle this case in any
00317            reasonably overflow-proof yet simple way.  So, this funtion
00318            will simply get it wrong whenever the three multiplies
00319            overflow 64-bits.  -CAS */
00320         if (a.denom < 0)
00321         {
00322             return ((a.num * -a.denom * b.denom) == b.num);
00323         }
00324         else
00325         {
00326             return (a.num == (b.num * a.denom * -b.denom));
00327         }
00328     }
00329 
00330     return ((a.num * b.denom) == (a.denom * b.num));
00331 }
00332 
00333 
00334 /* *******************************************************************
00335  *  gnc_numeric_same
00336  *  would a and b be equal() if they were both converted to the same 
00337  *  denominator? 
00338  ********************************************************************/
00339 
00340 int
00341 gnc_numeric_same (gnc_numeric a, gnc_numeric b, gint64 denom, gint how)
00342 {
00343     gnc_numeric aconv, bconv;
00344 
00345     aconv = gnc_numeric_convert (a, denom, how);
00346     bconv = gnc_numeric_convert (b, denom, how);
00347 
00348     return (gnc_numeric_equal (aconv, bconv));
00349 }
00350 
00351 
00352 
00353 /* *******************************************************************
00354  *  gnc_numeric_add
00355  ********************************************************************/
00356 
00357 gnc_numeric
00358 gnc_numeric_add (gnc_numeric a, gnc_numeric b, gint64 denom, gint how)
00359 {
00360     gnc_numeric sum;
00361 
00362     if (gnc_numeric_check (a) || gnc_numeric_check (b))
00363     {
00364         return gnc_numeric_error (GNC_ERROR_ARG);
00365     }
00366 
00367     if ((denom == GNC_DENOM_AUTO) &&
00368         (how & GNC_NUMERIC_DENOM_MASK) == GNC_HOW_DENOM_FIXED)
00369     {
00370         if (a.denom == b.denom)
00371         {
00372             denom = a.denom;
00373         }
00374         else if (b.num == 0)
00375         {
00376             denom = a.denom;
00377             b.denom = a.denom;
00378         }
00379         else if (a.num == 0)
00380         {
00381             denom = b.denom;
00382             a.denom = b.denom;
00383         }
00384         else
00385         {
00386             return gnc_numeric_error (GNC_ERROR_DENOM_DIFF);
00387         }
00388     }
00389 
00390     if (a.denom < 0)
00391     {
00392         a.num *= -a.denom;      /* BUG: overflow not handled.  */
00393         a.denom = 1;
00394     }
00395 
00396     if (b.denom < 0)
00397     {
00398         b.num *= -b.denom;      /* BUG: overflow not handled.  */
00399         b.denom = 1;
00400     }
00401 
00402     /* Get an exact answer.. same denominator is the common case. */
00403     if (a.denom == b.denom)
00404     {
00405         sum.num = a.num + b.num;    /* BUG: overflow not handled.  */
00406         sum.denom = a.denom;
00407     }
00408     else
00409     {
00410         /* We want to do this:
00411          *    sum.num = a.num*b.denom + b.num*a.denom;
00412          *    sum.denom = a.denom*b.denom;
00413          * but the multiply could overflow.  
00414          * Computing the LCD minimizes likelihood of overflow
00415          */
00416         gint64 lcd;
00417         qofint128 ca, cb, cab;
00418         lcd = gnc_numeric_lcd (a, b);
00419         if (GNC_ERROR_ARG == lcd)
00420         {
00421             return gnc_numeric_error (GNC_ERROR_OVERFLOW);
00422         }
00423         ca = mult128 (a.num, lcd / a.denom);
00424         if (ca.isbig)
00425             return gnc_numeric_error (GNC_ERROR_OVERFLOW);
00426 
00427         cb = mult128 (b.num, lcd / b.denom);
00428         if (cb.isbig)
00429             return gnc_numeric_error (GNC_ERROR_OVERFLOW);
00430 
00431         cab = add128 (ca, cb);
00432         if (cab.isbig)
00433             return gnc_numeric_error (GNC_ERROR_OVERFLOW);
00434 
00435         sum.num = cab.lo;
00436         if (cab.isneg)
00437             sum.num = -sum.num;
00438         sum.denom = lcd;
00439     }
00440 
00441     if ((denom == GNC_DENOM_AUTO) &&
00442         ((how & GNC_NUMERIC_DENOM_MASK) == GNC_HOW_DENOM_LCD))
00443     {
00444         denom = gnc_numeric_lcd (a, b);
00445         how = how & GNC_NUMERIC_RND_MASK;
00446     }
00447 
00448     return gnc_numeric_convert (sum, denom, how);
00449 }
00450 
00451 /* *******************************************************************
00452  *  gnc_numeric_sub
00453  ********************************************************************/
00454 
00455 gnc_numeric
00456 gnc_numeric_sub (gnc_numeric a, gnc_numeric b, gint64 denom, gint how)
00457 {
00458     gnc_numeric nb;
00459     if (gnc_numeric_check (a) || gnc_numeric_check (b))
00460     {
00461         return gnc_numeric_error (GNC_ERROR_ARG);
00462     }
00463 
00464     nb = b;
00465     nb.num = -nb.num;
00466     return gnc_numeric_add (a, nb, denom, how);
00467 }
00468 
00469 /* *******************************************************************
00470  *  gnc_numeric_mul
00471  ********************************************************************/
00472 
00473 gnc_numeric
00474 gnc_numeric_mul (gnc_numeric a, gnc_numeric b, gint64 denom, gint how)
00475 {
00476     gnc_numeric product, result;
00477     qofint128 bignume, bigdeno;
00478 
00479     if (gnc_numeric_check (a) || gnc_numeric_check (b))
00480     {
00481         return gnc_numeric_error (GNC_ERROR_ARG);
00482     }
00483 
00484     if ((denom == GNC_DENOM_AUTO) &&
00485         (how & GNC_NUMERIC_DENOM_MASK) == GNC_HOW_DENOM_FIXED)
00486     {
00487         if (a.denom == b.denom)
00488         {
00489             denom = a.denom;
00490         }
00491         else if (b.num == 0)
00492         {
00493             denom = a.denom;
00494         }
00495         else if (a.num == 0)
00496         {
00497             denom = b.denom;
00498         }
00499         else
00500         {
00501             return gnc_numeric_error (GNC_ERROR_DENOM_DIFF);
00502         }
00503     }
00504 
00505     if ((denom == GNC_DENOM_AUTO) &&
00506         ((how & GNC_NUMERIC_DENOM_MASK) == GNC_HOW_DENOM_LCD))
00507     {
00508         denom = gnc_numeric_lcd (a, b);
00509         how = how & GNC_NUMERIC_RND_MASK;
00510     }
00511 
00512     if (a.denom < 0)
00513     {
00514         a.num *= -a.denom;      /* BUG: overflow not handled.  */
00515         a.denom = 1;
00516     }
00517 
00518     if (b.denom < 0)
00519     {
00520         b.num *= -b.denom;      /* BUG: overflow not handled.  */
00521         b.denom = 1;
00522     }
00523 
00524     bignume = mult128 (a.num, b.num);
00525     bigdeno = mult128 (a.denom, b.denom);
00526     product.num = a.num * b.num;
00527     product.denom = a.denom * b.denom;
00528 
00529     /* If it looks to be overflowing, try to reduce the fraction ... */
00530     if (bignume.isbig || bigdeno.isbig)
00531     {
00532         gint64 tmp;
00533         a = gnc_numeric_reduce (a);
00534         b = gnc_numeric_reduce (b);
00535         tmp = a.num;
00536         a.num = b.num;
00537         b.num = tmp;
00538         a = gnc_numeric_reduce (a);
00539         b = gnc_numeric_reduce (b);
00540 
00541         bignume = mult128 (a.num, b.num);
00542         bigdeno = mult128 (a.denom, b.denom);
00543         product.num = a.num * b.num;
00544         product.denom = a.denom * b.denom;
00545     }
00546 
00547     /* If it its still overflowing, and rounding is allowed then round */
00548     if (bignume.isbig || bigdeno.isbig)
00549     {
00550         /* If rounding allowed, then shift until there's no 
00551          * more overflow. The conversion at the end will fix 
00552          * things up for the final value. Else overflow. */
00553         if ((how & GNC_NUMERIC_RND_MASK) == GNC_HOW_RND_NEVER)
00554         {
00555             if (bigdeno.isbig)
00556             {
00557                 return gnc_numeric_error (GNC_ERROR_OVERFLOW);
00558             }
00559             product = reduce128 (bignume, product.denom);
00560             if (gnc_numeric_check (product))
00561             {
00562                 return gnc_numeric_error (GNC_ERROR_OVERFLOW);
00563             }
00564         }
00565         else
00566         {
00567             while (bignume.isbig || bigdeno.isbig)
00568             {
00569                 bignume = shift128 (bignume);
00570                 bigdeno = shift128 (bigdeno);
00571             }
00572             product.num = bignume.lo;
00573             if (bignume.isneg)
00574                 product.num = -product.num;
00575 
00576             product.denom = bigdeno.lo;
00577             if (0 == product.denom)
00578             {
00579                 return gnc_numeric_error (GNC_ERROR_OVERFLOW);
00580             }
00581         }
00582     }
00583 
00584 #if 0                           /* currently, product denom won't ever be zero */
00585     if (product.denom < 0)
00586     {
00587         product.num = -product.num;
00588         product.denom = -product.denom;
00589     }
00590 #endif
00591 
00592     result = gnc_numeric_convert (product, denom, how);
00593     return result;
00594 }
00595 
00596 
00597 /* *******************************************************************
00598  *  gnc_numeric_div
00599  ********************************************************************/
00600 
00601 gnc_numeric
00602 gnc_numeric_div (gnc_numeric a, gnc_numeric b, gint64 denom, gint how)
00603 {
00604     gnc_numeric quotient;
00605     qofint128 nume, deno;
00606 
00607     if (gnc_numeric_check (a) || gnc_numeric_check (b))
00608     {
00609         return gnc_numeric_error (GNC_ERROR_ARG);
00610     }
00611 
00612     if ((denom == GNC_DENOM_AUTO) &&
00613         (how & GNC_NUMERIC_DENOM_MASK) == GNC_HOW_DENOM_FIXED)
00614     {
00615         if (a.denom == b.denom)
00616         {
00617             denom = a.denom;
00618         }
00619         else if (a.denom == 0)
00620         {
00621             denom = b.denom;
00622         }
00623         else
00624         {
00625             return gnc_numeric_error (GNC_ERROR_DENOM_DIFF);
00626         }
00627     }
00628 
00629 
00630     if (a.denom < 0)
00631     {
00632         a.num *= -a.denom;      /* BUG: overflow not handled.  */
00633         a.denom = 1;
00634     }
00635 
00636     if (b.denom < 0)
00637     {
00638         b.num *= -b.denom;      /* BUG: overflow not handled.  */
00639         b.denom = 1;
00640     }
00641 
00642     if (a.denom == b.denom)
00643     {
00644         quotient.num = a.num;
00645         quotient.denom = b.num;
00646     }
00647     else
00648     {
00649         gint64 sgn = 1;
00650         if (0 > a.num)
00651         {
00652             sgn = -sgn;
00653             a.num = -a.num;
00654         }
00655         if (0 > b.num)
00656         {
00657             sgn = -sgn;
00658             b.num = -b.num;
00659         }
00660         nume = mult128 (a.num, b.denom);
00661         deno = mult128 (b.num, a.denom);
00662 
00663         /* Try to avoid overflow by removing common factors */
00664         if (nume.isbig && deno.isbig)
00665         {
00666             gnc_numeric ra = gnc_numeric_reduce (a);
00667             gnc_numeric rb = gnc_numeric_reduce (b);
00668 
00669             gint64 gcf_nume = gcf64 (ra.num, rb.num);
00670             gint64 gcf_deno = gcf64 (rb.denom, ra.denom);
00671             nume = mult128 (ra.num / gcf_nume, rb.denom / gcf_deno);
00672             deno = mult128 (rb.num / gcf_nume, ra.denom / gcf_deno);
00673         }
00674 
00675         if ((0 == nume.isbig) && (0 == deno.isbig))
00676         {
00677             quotient.num = sgn * nume.lo;
00678             quotient.denom = deno.lo;
00679             goto dive_done;
00680         }
00681         else if (0 == deno.isbig)
00682         {
00683             quotient = reduce128 (nume, deno.lo);
00684             if (0 == gnc_numeric_check (quotient))
00685             {
00686                 quotient.num *= sgn;
00687                 goto dive_done;
00688             }
00689         }
00690 
00691         /* If rounding allowed, then shift until there's no 
00692          * more overflow. The conversion at the end will fix 
00693          * things up for the final value. */
00694         if ((how & GNC_NUMERIC_RND_MASK) == GNC_HOW_RND_NEVER)
00695         {
00696             return gnc_numeric_error (GNC_ERROR_OVERFLOW);
00697         }
00698         while (nume.isbig || deno.isbig)
00699         {
00700             nume = shift128 (nume);
00701             deno = shift128 (deno);
00702         }
00703         quotient.num = sgn * nume.lo;
00704         quotient.denom = deno.lo;
00705         if (0 == quotient.denom)
00706         {
00707             return gnc_numeric_error (GNC_ERROR_OVERFLOW);
00708         }
00709     }
00710 
00711     if (quotient.denom < 0)
00712     {
00713         quotient.num = -quotient.num;
00714         quotient.denom = -quotient.denom;
00715     }
00716 
00717   dive_done:
00718     if ((denom == GNC_DENOM_AUTO) &&
00719         ((how & GNC_NUMERIC_DENOM_MASK) == GNC_HOW_DENOM_LCD))
00720     {
00721         denom = gnc_numeric_lcd (a, b);
00722         how = how & GNC_NUMERIC_RND_MASK;
00723     }
00724 
00725     return gnc_numeric_convert (quotient, denom, how);
00726 }
00727 
00728 /* *******************************************************************
00729  *  gnc_numeric_neg
00730  *  negate the argument 
00731  ********************************************************************/
00732 
00733 gnc_numeric
00734 gnc_numeric_neg (gnc_numeric a)
00735 {
00736     if (gnc_numeric_check (a))
00737     {
00738         return gnc_numeric_error (GNC_ERROR_ARG);
00739     }
00740     return gnc_numeric_create (-a.num, a.denom);
00741 }
00742 
00743 /* *******************************************************************
00744  *  gnc_numeric_neg
00745  *  return the absolute value of the argument 
00746  ********************************************************************/
00747 
00748 gnc_numeric
00749 gnc_numeric_abs (gnc_numeric a)
00750 {
00751     if (gnc_numeric_check (a))
00752     {
00753         return gnc_numeric_error (GNC_ERROR_ARG);
00754     }
00755     return gnc_numeric_create (ABS (a.num), a.denom);
00756 }
00757 
00758 /* *******************************************************************
00759  *  gnc_numeric_convert
00760  ********************************************************************/
00761 
00762 gnc_numeric
00763 gnc_numeric_convert (gnc_numeric in, gint64 denom, gint how)
00764 {
00765     gnc_numeric out;
00766     gnc_numeric temp;
00767     gint64 temp_bc;
00768     gint64 temp_a;
00769     gint64 remainder;
00770     gint64 sign;
00771     gint denom_neg = 0;
00772     double ratio, logratio;
00773     double sigfigs;
00774     qofint128 nume, newm;
00775 
00776     temp.num = 0;
00777     temp.denom = 0;
00778 
00779     if (gnc_numeric_check (in))
00780     {
00781         return gnc_numeric_error (GNC_ERROR_ARG);
00782     }
00783 
00784     if (denom == GNC_DENOM_AUTO)
00785     {
00786         switch (how & GNC_NUMERIC_DENOM_MASK)
00787         {
00788         default:
00789         case GNC_HOW_DENOM_LCD: /* LCD is meaningless with AUTO in here */
00790         case GNC_HOW_DENOM_EXACT:
00791             return in;
00792             break;
00793 
00794         case GNC_HOW_DENOM_REDUCE:
00795             /* reduce the input to a relatively-prime fraction */
00796             return gnc_numeric_reduce (in);
00797             break;
00798 
00799         case GNC_HOW_DENOM_FIXED:
00800             if (in.denom != denom)
00801             {
00802                 return gnc_numeric_error (GNC_ERROR_DENOM_DIFF);
00803             }
00804             else
00805             {
00806                 return in;
00807             }
00808             break;
00809 
00810         case GNC_HOW_DENOM_SIGFIG:
00811             ratio = fabs (gnc_numeric_to_double (in));
00812             if (ratio < 10e-20)
00813             {
00814                 logratio = 0;
00815             }
00816             else
00817             {
00818                 logratio = log10 (ratio);
00819                 logratio = ((logratio > 0.0) ?
00820                     (floor (logratio) + 1.0) : (ceil (logratio)));
00821             }
00822             sigfigs = GNC_HOW_GET_SIGFIGS (how);
00823 
00824             if (sigfigs - logratio >= 0)
00825             {
00826                 denom = (gint64) (pow (10, sigfigs - logratio));
00827             }
00828             else
00829             {
00830                 denom = -((gint64) (pow (10, logratio - sigfigs)));
00831             }
00832 
00833             how = how & ~GNC_HOW_DENOM_SIGFIG & ~GNC_NUMERIC_SIGFIGS_MASK;
00834             break;
00835 
00836         }
00837     }
00838 
00839     /* Make sure we need to do the work */
00840     if (in.denom == denom)
00841     {
00842         return in;
00843     }
00844     if (in.num == 0)
00845     {
00846         out.num = 0;
00847         out.denom = denom;
00848         return out;
00849     }
00850 
00851     /* If the denominator of the input value is negative, get rid of that. */
00852     if (in.denom < 0)
00853     {
00854         in.num = in.num * (-in.denom);  /* BUG: overflow not handled.  */
00855         in.denom = 1;
00856     }
00857 
00858     sign = (in.num < 0) ? -1 : 1;
00859 
00860     /* If the denominator is less than zero, we are to interpret it as 
00861      * the reciprocal of its magnitude. */
00862     if (denom < 0)
00863     {
00864 
00865         /* XXX FIXME: use 128-bit math here ... */
00866         denom = -denom;
00867         denom_neg = 1;
00868         temp_a = (in.num < 0) ? -in.num : in.num;
00869         temp_bc = in.denom * denom; /* BUG: overflow not handled.  */
00870         remainder = temp_a % temp_bc;
00871         out.num = temp_a / temp_bc;
00872         out.denom = -denom;
00873     }
00874     else
00875     {
00876         /* Do all the modulo and int division on positive values to make
00877          * things a little clearer. Reduce the fraction denom/in.denom to
00878          * help with range errors */
00879         temp.num = denom;
00880         temp.denom = in.denom;
00881         temp = gnc_numeric_reduce (temp);
00882 
00883         /* Symbolically, do the following:
00884          * out.num   = in.num * temp.num;
00885          * remainder = out.num % temp.denom;
00886          * out.num   = out.num / temp.denom;
00887          * out.denom = denom;
00888          */
00889         nume = mult128 (in.num, temp.num);
00890         newm = div128 (nume, temp.denom);
00891         remainder = rem128 (nume, temp.denom);
00892 
00893         if (newm.isbig)
00894         {
00895             return gnc_numeric_error (GNC_ERROR_OVERFLOW);
00896         }
00897 
00898         out.num = newm.lo;
00899         out.denom = denom;
00900     }
00901 
00902     if (remainder)
00903     {
00904         switch (how & GNC_NUMERIC_RND_MASK)
00905         {
00906         case GNC_HOW_RND_FLOOR:
00907             if (sign < 0)
00908             {
00909                 out.num = out.num + 1;
00910             }
00911             break;
00912 
00913         case GNC_HOW_RND_CEIL:
00914             if (sign > 0)
00915             {
00916                 out.num = out.num + 1;
00917             }
00918             break;
00919 
00920         case GNC_HOW_RND_TRUNC:
00921             break;
00922 
00923         case GNC_HOW_RND_PROMOTE:
00924             out.num = out.num + 1;
00925             break;
00926 
00927         case GNC_HOW_RND_ROUND_HALF_DOWN:
00928             if (denom_neg)
00929             {
00930                 if ((2 * remainder) > in.denom * denom)
00931                 {
00932                     out.num = out.num + 1;
00933                 }
00934             }
00935             else if ((2 * remainder) > temp.denom)
00936             {
00937                 out.num = out.num + 1;
00938             }
00939             /* check that 2*remainder didn't over-flow */
00940             else if (((2 * remainder) < remainder) &&
00941                 (remainder > (temp.denom / 2)))
00942             {
00943                 out.num = out.num + 1;
00944             }
00945             break;
00946 
00947         case GNC_HOW_RND_ROUND_HALF_UP:
00948             if (denom_neg)
00949             {
00950                 if ((2 * remainder) >= in.denom * denom)
00951                 {
00952                     out.num = out.num + 1;
00953                 }
00954             }
00955             else if ((2 * remainder) >= temp.denom)
00956             {
00957                 out.num = out.num + 1;
00958             }
00959             /* check that 2*remainder didn't over-flow */
00960             else if (((2 * remainder) < remainder) &&
00961                 (remainder >= (temp.denom / 2)))
00962             {
00963                 out.num = out.num + 1;
00964             }
00965             break;
00966 
00967         case GNC_HOW_RND_ROUND:
00968             if (denom_neg)
00969             {
00970                 if ((2 * remainder) > in.denom * denom)
00971                 {
00972                     out.num = out.num + 1;
00973                 }
00974                 else if ((2 * remainder) == in.denom * denom)
00975                 {
00976                     if (out.num % 2)
00977                     {
00978                         out.num = out.num + 1;
00979                     }
00980                 }
00981             }
00982             else
00983             {
00984                 if ((2 * remainder) > temp.denom)
00985                 {
00986                     out.num = out.num + 1;
00987                 }
00988                 /* check that 2*remainder didn't over-flow */
00989                 else if (((2 * remainder) < remainder) &&
00990                     (remainder > (temp.denom / 2)))
00991                 {
00992                     out.num = out.num + 1;
00993                 }
00994                 else if ((2 * remainder) == temp.denom)
00995                 {
00996                     if (out.num % 2)
00997                     {
00998                         out.num = out.num + 1;
00999                     }
01000                 }
01001                 /* check that 2*remainder didn't over-flow */
01002                 else if (((2 * remainder) < remainder) &&
01003                     (remainder == (temp.denom / 2)))
01004                 {
01005                     if (out.num % 2)
01006                     {
01007                         out.num = out.num + 1;
01008                     }
01009                 }
01010             }
01011             break;
01012 
01013         case GNC_HOW_RND_NEVER:
01014             return gnc_numeric_error (GNC_ERROR_REMAINDER);
01015             break;
01016         }
01017     }
01018 
01019     out.num = (sign > 0) ? out.num : (-out.num);
01020 
01021     return out;
01022 }
01023 
01024 
01025 /* *******************************************************************
01026  *  reduce a fraction by GCF elimination.  This is NOT done as a
01027  *  part of the arithmetic API unless GNC_HOW_DENOM_REDUCE is specified 
01028  *  as the output denominator.
01029  ********************************************************************/
01030 
01031 gnc_numeric
01032 gnc_numeric_reduce (gnc_numeric in)
01033 {
01034     gint64 t;
01035     gint64 num = (in.num < 0) ? (-in.num) : in.num;
01036     gint64 denom = in.denom;
01037     gnc_numeric out;
01038 
01039     if (gnc_numeric_check (in))
01040     {
01041         return gnc_numeric_error (GNC_ERROR_ARG);
01042     }
01043 
01044     /* The strategy is to use Euclid's algorithm */
01045     while (denom > 0)
01046     {
01047         t = num % denom;
01048         num = denom;
01049         denom = t;
01050     }
01051     /* num now holds the GCD (Greatest Common Divisor) */
01052 
01053     /* All calculations are done on positive num, since it's not 
01054      * well defined what % does for negative values */
01055     out.num = in.num / num;
01056     out.denom = in.denom / num;
01057     return out;
01058 }
01059 
01060 /* *******************************************************************
01061  *  double_to_gnc_numeric
01062  ********************************************************************/
01063 
01064 gnc_numeric
01065 double_to_gnc_numeric (double in, gint64 denom, gint how)
01066 {
01067     gnc_numeric out;
01068     gint64 int_part = 0;
01069     double frac_part;
01070     gint64 frac_int = 0;
01071     double logval;
01072     double sigfigs;
01073 
01074     if ((denom == GNC_DENOM_AUTO) && (how & GNC_HOW_DENOM_SIGFIG))
01075     {
01076         if (fabs (in) < 10e-20)
01077         {
01078             logval = 0;
01079         }
01080         else
01081         {
01082             logval = log10 (fabs (in));
01083             logval = ((logval > 0.0) ?
01084                 (floor (logval) + 1.0) : (ceil (logval)));
01085         }
01086         sigfigs = GNC_HOW_GET_SIGFIGS (how);
01087         if (sigfigs - logval >= 0)
01088         {
01089             denom = (gint64) (pow (10, sigfigs - logval));
01090         }
01091         else
01092         {
01093             denom = -((gint64) (pow (10, logval - sigfigs)));
01094         }
01095 
01096         how = how & ~GNC_HOW_DENOM_SIGFIG & ~GNC_NUMERIC_SIGFIGS_MASK;
01097     }
01098 
01099     int_part = (gint64) (floor (fabs (in)));
01100     frac_part = in - (double) int_part;
01101 
01102     int_part = int_part * denom;
01103     frac_part = frac_part * (double) denom;
01104 
01105     switch (how & GNC_NUMERIC_RND_MASK)
01106     {
01107     case GNC_HOW_RND_FLOOR:
01108         frac_int = (gint64) floor (frac_part);
01109         break;
01110 
01111     case GNC_HOW_RND_CEIL:
01112         frac_int = (gint64) ceil (frac_part);
01113         break;
01114 
01115     case GNC_HOW_RND_TRUNC:
01116         frac_int = (gint64) frac_part;
01117         break;
01118 
01119     case GNC_HOW_RND_ROUND:
01120     case GNC_HOW_RND_ROUND_HALF_UP:
01121         frac_int = (gint64) rint (frac_part);
01122         break;
01123 
01124     case GNC_HOW_RND_NEVER:
01125         frac_int = (gint64) floor (frac_part);
01126         if (frac_part != (double) frac_int)
01127         {
01128             /* signal an error */
01129         }
01130         break;
01131     }
01132 
01133     out.num = int_part + frac_int;
01134     out.denom = denom;
01135     return out;
01136 }
01137 
01138 /* *******************************************************************
01139  *  gnc_numeric_to_double
01140  ********************************************************************/
01141 
01142 double
01143 gnc_numeric_to_double (gnc_numeric in)
01144 {
01145     if (in.denom > 0)
01146     {
01147         return (double) in.num / (double) in.denom;
01148     }
01149     else
01150     {
01151         return (double) (in.num * -in.denom);
01152     }
01153 }
01154 
01155 /* *******************************************************************
01156  *  gnc_numeric_error
01157  ********************************************************************/
01158 
01159 gnc_numeric
01160 gnc_numeric_error (GNCNumericErrorCode error_code)
01161 {
01162     return gnc_numeric_create (error_code, 0LL);
01163 }
01164 
01165 
01166 /* *******************************************************************
01167  *  gnc_numeric_add_with_error
01168  ********************************************************************/
01169 
01170 gnc_numeric
01171 gnc_numeric_add_with_error (gnc_numeric a, gnc_numeric b,
01172     gint64 denom, gint how, gnc_numeric * error)
01173 {
01174 
01175     gnc_numeric sum = gnc_numeric_add (a, b, denom, how);
01176     gnc_numeric exact = gnc_numeric_add (a, b, GNC_DENOM_AUTO,
01177         GNC_HOW_DENOM_REDUCE);
01178     gnc_numeric err = gnc_numeric_sub (sum, exact, GNC_DENOM_AUTO,
01179         GNC_HOW_DENOM_REDUCE);
01180 
01181     if (error)
01182     {
01183         *error = err;
01184     }
01185     return sum;
01186 }
01187 
01188 /* *******************************************************************
01189  *  gnc_numeric_sub_with_error
01190  ********************************************************************/
01191 
01192 gnc_numeric
01193 gnc_numeric_sub_with_error (gnc_numeric a, gnc_numeric b,
01194     gint64 denom, gint how, gnc_numeric * error)
01195 {
01196     gnc_numeric diff = gnc_numeric_sub (a, b, denom, how);
01197     gnc_numeric exact = gnc_numeric_sub (a, b, GNC_DENOM_AUTO,
01198         GNC_HOW_DENOM_REDUCE);
01199     gnc_numeric err = gnc_numeric_sub (diff, exact, GNC_DENOM_AUTO,
01200         GNC_HOW_DENOM_REDUCE);
01201     if (error)
01202     {
01203         *error = err;
01204     }
01205     return diff;
01206 }
01207 
01208 
01209 /* *******************************************************************
01210  *  gnc_numeric_mul_with_error
01211  ********************************************************************/
01212 
01213 gnc_numeric
01214 gnc_numeric_mul_with_error (gnc_numeric a, gnc_numeric b,
01215     gint64 denom, gint how, gnc_numeric * error)
01216 {
01217     gnc_numeric prod = gnc_numeric_mul (a, b, denom, how);
01218     gnc_numeric exact = gnc_numeric_mul (a, b, GNC_DENOM_AUTO,
01219         GNC_HOW_DENOM_REDUCE);
01220     gnc_numeric err = gnc_numeric_sub (prod, exact, GNC_DENOM_AUTO,
01221         GNC_HOW_DENOM_REDUCE);
01222     if (error)
01223     {
01224         *error = err;
01225     }
01226     return prod;
01227 }
01228 
01229 
01230 /* *******************************************************************
01231  *  gnc_numeric_div_with_error
01232  ********************************************************************/
01233 
01234 gnc_numeric
01235 gnc_numeric_div_with_error (gnc_numeric a, gnc_numeric b,
01236     gint64 denom, gint how, gnc_numeric * error)
01237 {
01238     gnc_numeric quot = gnc_numeric_div (a, b, denom, how);
01239     gnc_numeric exact = gnc_numeric_div (a, b, GNC_DENOM_AUTO,
01240         GNC_HOW_DENOM_REDUCE);
01241     gnc_numeric err = gnc_numeric_sub (quot, exact,
01242         GNC_DENOM_AUTO, GNC_HOW_DENOM_REDUCE);
01243     if (error)
01244     {
01245         *error = err;
01246     }
01247     return quot;
01248 }
01249 
01250 /* *******************************************************************
01251  *  gnc_numeric text IO
01252  ********************************************************************/
01253 
01254 gchar *
01255 gnc_numeric_to_string (gnc_numeric n)
01256 {
01257     gchar *result;
01258     gint64 tmpnum = n.num;
01259     gint64 tmpdenom = n.denom;
01260 
01261     result =
01262         g_strdup_printf ("%" G_GINT64_FORMAT "/%" G_GINT64_FORMAT, tmpnum,
01263         tmpdenom);
01264 
01265     return result;
01266 }
01267 
01268 gchar *
01269 gnc_num_dbg_to_string (gnc_numeric n)
01270 {
01271     static char buff[1000];
01272     static char *p = buff;
01273     gint64 tmpnum = n.num;
01274     gint64 tmpdenom = n.denom;
01275 
01276     p += 100;
01277     if (p - buff >= 1000)
01278         p = buff;
01279 
01280     sprintf (p, "%" G_GINT64_FORMAT "/%" G_GINT64_FORMAT, tmpnum,
01281         tmpdenom);
01282 
01283     return p;
01284 }
01285 
01286 gboolean
01287 string_to_gnc_numeric (const gchar * str, gnc_numeric * n)
01288 {
01289     size_t num_read;
01290     gint64 tmpnum;
01291     gint64 tmpdenom;
01292 
01293     if (!str)
01294         return FALSE;
01295 
01296 #ifdef GNC_DEPRECATED
01297     /* must use "<" here because %n's effects aren't well defined */
01298     if (sscanf (str, " " GNC_SCANF_LLD "/" GNC_SCANF_LLD "%n",
01299             &tmpnum, &tmpdenom, &num_read) < 2)
01300     {
01301         return FALSE;
01302     }
01303 #else
01304     tmpnum = strtoll (str, NULL, 0);
01305     str = strchr (str, '/');
01306     if (!str)
01307         return FALSE;
01308     str++;
01309     tmpdenom = strtoll (str, NULL, 0);
01310     num_read = strspn (str, "0123456789");
01311 #endif
01312     n->num = tmpnum;
01313     n->denom = tmpdenom;
01314     return TRUE;
01315 }
01316 
01317 /* *******************************************************************
01318  *  gnc_numeric misc testing
01319  ********************************************************************/
01320 #ifdef _GNC_NUMERIC_TEST
01321 
01322 static char *
01323 gnc_numeric_print (gnc_numeric in)
01324 {
01325     char *retval;
01326     if (gnc_numeric_check (in))
01327     {
01328         retval =
01329             g_strdup_printf ("<ERROR> [%" G_GINT64_FORMAT " / %"
01330             G_GINT64_FORMAT "]", in.num, in.denom);
01331     }
01332     else
01333     {
01334         retval =
01335             g_strdup_printf ("[%" G_GINT64_FORMAT " / %" G_GINT64_FORMAT
01336             "]", in.num, in.denom);
01337     }
01338     return retval;
01339 }
01340 
01341 int
01342 main (int argc, char **argv)
01343 {
01344     gnc_numeric a = gnc_numeric_create (1, 3);
01345     gnc_numeric b = gnc_numeric_create (1, 4);
01346     gnc_numeric c;
01347 
01348     gnc_numeric err;
01349 
01350     c = gnc_numeric_add_with_error (a, b, 100, GNC_HOW_RND_ROUND, &err);
01351     printf ("add 100ths/error : %s + %s = %s + (error) %s\n\n",
01352         gnc_numeric_print (a), gnc_numeric_print (b),
01353         gnc_numeric_print (c), gnc_numeric_print (err));
01354 
01355     c = gnc_numeric_sub_with_error (a, b, 100, GNC_HOW_RND_FLOOR, &err);
01356     printf ("sub 100ths/error : %s - %s = %s + (error) %s\n\n",
01357         gnc_numeric_print (a), gnc_numeric_print (b),
01358         gnc_numeric_print (c), gnc_numeric_print (err));
01359 
01360     c = gnc_numeric_mul_with_error (a, b, 100, GNC_HOW_RND_ROUND, &err);
01361     printf ("mul 100ths/error : %s * %s = %s + (error) %s\n\n",
01362         gnc_numeric_print (a), gnc_numeric_print (b),
01363         gnc_numeric_print (c), gnc_numeric_print (err));
01364 
01365     c = gnc_numeric_div_with_error (a, b, 100, GNC_HOW_RND_ROUND, &err);
01366     printf ("div 100ths/error : %s / %s = %s + (error) %s\n\n",
01367         gnc_numeric_print (a), gnc_numeric_print (b),
01368         gnc_numeric_print (c), gnc_numeric_print (err));
01369 
01370     printf ("multiply (EXACT): %s * %s = %s\n",
01371         gnc_numeric_print (a), gnc_numeric_print (b),
01372         gnc_numeric_print (gnc_numeric_mul
01373             (a, b, GNC_DENOM_AUTO, GNC_HOW_DENOM_EXACT)));
01374 
01375     printf ("multiply (REDUCE): %s * %s = %s\n",
01376         gnc_numeric_print (a), gnc_numeric_print (b),
01377         gnc_numeric_print (gnc_numeric_mul
01378             (a, b, GNC_DENOM_AUTO, GNC_HOW_DENOM_REDUCE)));
01379 
01380 
01381     return 0;
01382 }
01383 #endif
01384 
01385 /* ======================== END OF FILE =================== */

Generated on Fri Sep 1 15:13:11 2006 for QOF by  doxygen 1.4.7