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

Generated on Fri May 12 18:00:32 2006 for QOF by  doxygen 1.4.4