开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 数值算法/人工智能 > 查看源码
_bigfloat.c
资源名称：leda.tar.gz　[点击查看]
上传用户：gzelex
上传日期：2007-01-07
资源大小：707k
文件大小：29k
源码类别：

数值算法/人工智能

开发平台：

MultiPlatform

_bigfloat.c：源码内容
  1. /*******************************************************************************
  2. +
  3. +  LEDA-R  3.2.3
  4. +
  5. +  _bigfloat.c
  6. +
  7. +  Copyright (c) 1995  by  Max-Planck-Institut fuer Informatik
  8. +  Im Stadtwald, 66123 Saarbruecken, Germany     
  9. +  All rights reserved.
  10. + 
  11. *******************************************************************************/
  14. /* This file has been automatically generated from "bigfloat.w"
  15.    by CTANGLE (Version 3.1 [km2c]) */
  20. #include <LEDA/bigfloat.h>
  21. #include <iostream.h>
  22. #include <stdio.h>
  23. #include <math.h>
  24. #include <string.h>
  29. void binout(ostream & os, integer);
  30. void exact_decimal_output(ostream &, bigfloat, long = MAX_PREC);
  31. void decimal_output(ostream &, bigfloat);
  32. integer powl(const integer &, long);
  33. bigfloat dec_round(bigfloat, long);
  34. integer b_round(const integer &, long, rounding_modes = TO_NEAREST);
  35. void cut(integer &, long);
  36. long bcomp(integer, integer);
  37. long sigcomp(integer, integer);
  38. long test_for_one(integer, long, long);
  39. inline long get_bit(const integer &, long);
  40. inline void set_bit(long &, long);
  41. double compose_parts(long, long, long, long);
  42. double pow2(long);
  47. const double double_min = pow2(long(-1022));
  48. const double NaN_double = compose_parts(0, 2047, 0, 1);
  49. const double pInf_double = compose_parts(0, 2047, 0, 0);
  50. const double Infinity = pInf_double;
  51. const double nInf_double = -pInf_double;
  52. const double pZero_double = compose_parts(0, 0, 0, 0);
  53. const double nZero_double = compose_parts(1, 0, 0, 0);
  55. const integer integer_1 = integer(1);
  56. const integer integer_52 = (integer_1 << 52) - integer_1;
  57. const integer integer_32 = (integer_1 << 32) - integer_1;
  58. const integer integer_20 = (integer_1 << 20) - integer_1;
  60. rounding_modes round_mode = TO_NEAREST;
  62. long global_prec = 52;
  63. rounding_modes global_mode;
  64. bool buf_exact;
  65. long output = DEC_OUT;
  70. void binout(ostream & os, integer b)
  71. {
  72.   char temp[bin_maxlen];
  73.   long flag = 0, count = 0;
  75.   if (b < 0)
  76.     b *= -1;
  78.   do {
  79.     temp[count++] = (char) (b % 2).tolong() + '0';
  80.     if (b <= 1)
  81.       flag = 1;
  82.     else
  83.       b /= 2;
  84.   } while (!flag);
  85.   for (long i = count - 1; i >= 0; i--)
  86.     os << temp[i];
  87. }
  89. integer b_round(const integer & b, long digits, rounding_modes mode)
  90. {
  91.   bigfloat bb(b, 0);
  92.   bb.round(digits, mode);
  93.   return bb.get_significant();
  94. }
  96. void cut(integer & b, long prec)
  97. {
  98.   if (b.length() > prec)
  99.     b = (b >> (b.length() - prec));
  100. }
  102. long sigcomp(integer s1, integer s2)
  103. {
  104.   if ((s1 == 0) && (s2 == 0))
  105.     return 0;
  106.   if (sign(s1) != sign(s2)) {
  107.     if (sign(s1) == 0)
  108.       if (sign(s2) < 0)
  109. return 1;
  110.       else
  111. return -1;
  112.     if (sign(s2) == 0)
  113.       if (sign(s1) > 0)
  114. return 1;
  115.       else
  116. return -1;
  118.     if (s1 > 0)
  119.       return 1;
  120.     else
  121.       return -1;
  122.   }
  124.   return bcomp(s1, s2);
  125. }
  127. long bcomp(integer b1, integer b2)
  128. {
  129.   if (b1.length() > b2.length())
  130.     b2 = b2 << (b1.length() - b2.length());
  131.   else if (b1.length() < b2.length())
  132.     b1 = b1 << (b2.length() - b1.length());
  134.   if (b1 > b2)
  135.     return 1;
  136.   if (b1 < b2)
  137.     return -1;
  138.   return 0;
  139. }
  141. long test_for_one(integer b, long start, long end)
  142. {
  143.   long len = b.length();
  145.   if (end > len)
  146.     end = len;
  148.   for (long i = start; i <= end; i++) {
  149.     if (len != i) {
  150.       if (abs((b >> (len - i)) & 1) == 1)
  151. return 1;
  152.     }
  153.     else {
  154.       if (abs(b & 1) == 1)
  155. return 1;
  156.     }
  158.   }
  159.   return 0;
  160. }
  162. long get_bit(const integer & b, long digit)
  163. {
  164.   return abs(int (((b >> (b.length() - digit)) & 1).tolong()));
  165. }
  167. void set_bit(long &b, long digit)
  168. {
  169.   if (digit < 32)
  170.     b = b & (long) pow2(32 - digit);
  171. }
  173. void set_globprec(long p)
  174. {
  175.   global_prec = p;
  176. }
  178. void set_round_mode(rounding_modes m)
  179. {
  180.   round_mode = m;
  181. }
  183. bigfloat dec_round(bigfloat b, long prec)
  184. {
  185. /* calculate needed binary accuracy for operations div and mul */
  186.   long flag = 0, b_prec = (long) ((prec + 2) / log10(2)) + 1;
  188. /* keep sign */
  189.   if (sign(b) == -1) {
  190.     flag = 1;
  191.     b = fabs(b);
  192.   }
  194. /* decimal shift of b */
  195.   int i;
  196.   for (i = 0; i < prec; i++)
  197.     b = mul(b, 10, b_prec);
  199. /* test whether we have to round up or not */
  200.   bigfloat a = b.tointeger() + integer(1);
  201.   if ((a - b) <= 0.5)
  202.     b = bigfloat(a.tointeger(), b.get_exponent());
  203.   else
  204.     b = bigfloat(b.tointeger(), b.get_exponent());
  206. /* decimal shiftback */
  207.   for (i = 0; i < prec; i++)
  208.     b = div(b, 10, b_prec);
  210. /* sign correction */
  211.   if (flag)
  212.     b = (-b);
  214.   return b;
  215. }
  219. /*
  220. void exact_decimal_output (ostream &os,bigfloat b,long prec)
  221. {
  222.   if (sign(b)==-1)
  223.   {
  224.     os << "-";
  225.     b=fabs(b);
  226.   }
  228.   long max_pr=(long)(b.get_precision()*log10(2));
  230.   if ((prec>max_pr)||(prec==MAX_PREC))
  231.     prec=max_pr;
  233.   if (prec<max_pr)
  234.     b=dec_round(b,prec);
  236.   os << b.tointeger() << ".";
  237.   b=sub(b,b.tointeger(),1,EXACT);
  239.   for (prec;prec>0;prec--)
  240.   {
  241.     b=mul(b,10,1,EXACT);
  242.     os << b.tointeger();
  243.     b=sub(b,b.tointeger(),1,EXACT);
  244.   }
  245. }
  246. */
  249. void decimal_output(ostream & os, bigfloat b)
  250. {
  251.   if (!b.get_exponent().islong())
  252.     error_handler(1, "decimal_output: not implemented for large exponents");
  254.   if (sign(b) < 0) {
  255.     b = -b;
  256.     os << "-";
  257.   }
  259.   integer int_part = b.tointeger(TO_N_INF);
  260.   os << int_part;
  262.   bigfloat rest = sub(b, bigfloat(int_part), 32, EXACT);
  263.   if (rest == 0)
  264.     return;
  266.   integer least_prec = abs(rest.get_exponent()) + integer(rest.get_precision());
  267.   if (!least_prec.islong())
  268.     error_handler(1, "decimal_output: not implemented for large precisions");
  270.   os << ".";
  272.   long mul_log = (long) ceil(log10(2) * least_prec.tolong());
  273.   integer pow10 = powl(integer(10), mul_log);
  274.   bigfloat multiplicator = pow10;
  276.   bigfloat rest_sig = mul(rest, multiplicator, 32, EXACT);
  277.   os << rest_sig.tointeger(TO_NEAREST);
  279.   long dec_exp = (long) ceil(log10(2) * rest.get_exponent().tolong());
  280.   if (dec_exp != 0)
  281.     os << "E" << dec_exp;
  282. }
  285. integer powl(const integer & x, long n)
  286. {
  287.   integer y = x, z = 1;
  288.   long n_prefix = n;
  289.   while (n_prefix > 0) {
  290.     if (n_prefix % 2)
  291.       z = z * y;
  292.     n_prefix = n_prefix / 2;
  293.     y = y * y;
  294.   }
  295.   return z;
  296. }
  304. double pow2(long exp)
  305. {
  306.   return compose_parts(0, exp + 1023, 0, 0);
  307. }
  309. double compose_parts(long sign_1, long exp_11, long most_sig_20, long least_sig_32)
  310. {
  311.   double a;
  312.   long *p;
  313.   p = (long *) &a;
  315.   long highword = 0, lowword = 0;
  317. /* calculate highword */
  318.   if (sign_1)
  319.     highword = highword | 0x80000000;
  320.   exp_11 = exp_11 << 20;
  321.   highword = highword | exp_11;
  322.   highword = highword | most_sig_20;
  324. /* calculate lowword */
  325.   lowword = least_sig_32;
  327. /* put the word together */
  328.   (*p) = highword;
  329.   p++;
  330.   (*p) = lowword;
  332.   return a;
  333. }
  337. bigfloat outofchar(char *rep)
  338. {
  339. /* There are in principle two possible formats for the representation
  340.      of a float:   1) +/-abc.def
  341.                    2) +/-0.0100111e+/-11110
  342.   */
  344.   bigfloat temp;
  346. /* check whether first character is a sign */
  348.   long sign = positive;
  349.   if (rep[0] == '-') {
  350.     sign = negative;
  351.     rep++;
  352.   }
  353.   if (rep[0] == '+')
  354.     rep++;
  356. /* the big branch into the two main cases */
  358.   if (rep[0] == '0') {
  359. /* scan rep-string until point is passed */
  361.     long flag = 0;
  362.     while (!flag)
  363.       if (*(rep++) == '.')
  364. flag = 1;
  366. /* scan  significant up to 'e'-character */
  368.     long count = 0;
  369.     flag = 0;
  370.     while (!flag)
  371.       if ((rep[count++] == 'e') || (rep[count - 1] == 'E'))
  372. flag = 1;
  373.     count--;
  375. /* Now we know how many digits we have in the  significant */
  376. /* It is now possible two calculate our  significant */
  378.     long i;
  379.     for (i = 0; i < count; i++)
  380.       temp.significant += (long) ((*(rep++) - '0') * pow2(count - 1 - i));
  382. /* skip e/E */
  383.     rep++;
  385. /* scan exponent */
  387. /* check for sign */
  389.     long s = positive;
  390.     if (((*rep) == '-') || ((*rep) == '+'))
  391.       if ((*rep++) == '-')
  392. s = negative;
  394.     count = strlen(rep);
  395.     for (i = 0; i < count; i++)
  396.       temp.exponent += (long) ((*(rep++) - '0') * pow2(count - 1 - i));
  398.     if (s == negative)
  399.       temp.exponent *= (-1);
  400.   }
  401.   else {
  402.     long pr = strlen(rep) * 10;
  404. /* scan number up to decimal point */
  406.     while ((*rep) != '.')
  407.       temp = add(mul(temp, 10, pr), (*(rep++) - '0'), pr);
  408.     rep++;
  410.     integer dd = 10;
  411.     while ((*rep) != 0) {
  412.       temp = add(temp, div(*(rep++) - '0', bigfloat(dd), pr), pr);
  413.       dd *= 10;
  414.     }
  415.   }
  416.   if (sign == negative)
  417.     temp.significant *= -1;
  419.   temp.precision = temp.significant.length();
  420.   return temp;
  421. }
  428. bigfloat::bigfloat(const integer & m, const integer & e)
  429. {
  430.   significant = m;
  431.   exponent = e;
  432.   precision = significant.length();
  434.   if ((significant == 0) && (exponent != 1024) && (exponent != 1025) && (exponent != 1026)
  435.     && (exponent != -1023) && (exponent != -1024))
  436.     exponent = -1023;
  437. }
  438. bigfloat::bigfloat(const bigfloat & b)
  439. {
  440.   significant = b.significant;
  441.   exponent = b.exponent;
  442.   precision = b.precision;
  443. }
  444. bigfloat::bigfloat(double d)
  445. {
  446.   if (d != 0) {
  447. /* If d is lower than double_min we can be sure that d is denormalized.
  448.     Thus we can multiply d with $2^{1023}$ to get a normalized number. */
  450.     long flag = 0;
  451.     if (fabs(d) < double_min) {
  452.       d = d * pow2(long(1023));
  453.       flag = 1;
  454.     }
  456.     unsigned long mh = 0, ml = 0, sign = 0;
  458.     long e;
  459.     long *p;
  461. /* convert double a pointer to an integer */
  462.     p = (long *) &d;
  464. /* let mh represent the highword of the 64-bit field of the double */
  465.     mh = *p;
  466.     p++;
  468. /* let ml hold the least significant 32 digits of d */
  469.     ml = *p;
  471. /* binary and to get the 11 bits of the exponent */
  472.     e = mh & 0x7ff00000;
  473. /* shift the result to get the right value */
  474.     e = e >> 20;
  475. /* unbias the exponent */
  476.     e -= 1023;
  478.     e++; /* implicit one */
  480. /* get the sign */
  481.     sign = mh & 0x80000000;
  482.     sign = sign >> 31;
  484. /* get the significant bits out of the highword */
  485.     mh = mh & 0x000fffff;
  486. /* we have ensured that d is normalized $Rightarrow$ add the leading one
  487.        bit */
  488.     mh = mh | 0x00100000;
  490.     exponent = integer(e);
  492. /* check if d was denormalized */
  493.     if (flag)
  494.       exponent -= 1023;
  496. /* compose significant */
  497.     significant = integer(mh);
  498.     significant = significant << 32;
  499.     significant = significant + integer(ml);
  501.     if (sign == 1)
  502.       significant = significant * integer(-1);
  503.     precision = significant.length();
  505.     if ((exponent == 1025) && (significant != 0)) { /* NaN-Case */
  506.       significant = 0;
  507.       exponent = 1026;
  508.     }
  509.     if ((exponent == 1025) && (significant == 0)) { /* Inf-cases */
  510.       if (sign == 1)
  511. exponent = 1025;
  512.       else
  513. exponent = 1024;
  514.       significant = 0;
  515.     }
  516.     if (exponent == -1022) { /* Zero-cases */
  517.       if (sign == 1)
  518. exponent = -1024;
  519.       exponent = -1023;
  520.     }
  521.   }
  522.   else {
  523.     significant = 0;
  524.     exponent = 0;
  525.     precision = 0;
  526.   }
  527. }
  528. bigfloat::bigfloat(const integer & b)
  529. {
  530.   significant = b;
  531.   if (b == 0)
  532.     exponent = -1023;
  533.   else
  534.     exponent = significant.length();
  535.   precision = exponent.tolong();
  536. }
  537. bigfloat::bigfloat(long a)
  538. {
  539.   significant = a;
  540.   if (a == 0)
  541.     exponent = -1023;
  542.   else
  543.     exponent = significant.length();
  544.   precision = exponent.tolong();
  545. }
  546. bigfloat::bigfloat(int a)
  547. {
  548.   significant = a;
  549.   if (a == 0)
  550.     exponent = -1023;
  551.   else
  552.     exponent = significant.length();
  553.   precision = exponent.tolong();
  554. }
  561. void bigfloat::set_precision(long prec, rounding_modes mode)
  562. {
  564.   if (mode != EXACT) {
  565.     if (prec > precision)
  566.       significant = significant << (prec - precision);
  567.     else if (precision > prec) {
  568. /* if the new precision is lower then the old one we firstly have to
  569.          round the old  significant to the new precision */
  571.       round(prec, mode);
  572.     }
  573.     precision = prec;
  574.   }
  575. }
  579. long bigfloat::get_precision(void)
  580. {
  581.   return precision;
  582. }
  586. integer bigfloat::tointeger(rounding_modes rmode) const
  587. {
  588.   if ((isNaN(*this)) || (isInf(*this)))
  589.     error_handler(1, "bigfloat::tointeger : special values cannot be converted to integer");
  590.   if (isZero(*this))
  591.     return integer(0);
  593.   return b_round(significant, exponent.tolong(), rmode);
  594. }
  598. double bigfloat::todouble(void) const
  599. {
  600.   long s = 0; /* sign bit of returned value */
  601.   long t_exp = 0; /* exponent of returned value */
  602.   integer t_sig = 0; /* significant of returned value */
  606. /* special case checking */
  608.   if (isNaN(*this))
  609.     return NaN_double;
  610.   if (ispInf(*this))
  611.     return pInf_double;
  612.   if (isnInf(*this))
  613.     return nInf_double;
  614.   if (ispZero(*this))
  615.     return pZero_double;
  616.   if (isnZero(*this))
  617.     return nZero_double;
  623.   long denorm = 0;
  625.   if (exponent > 1024)
  626.     return non_zero_sign(*this) * pInf_double;
  627.   if (exponent < -1074)
  628.     return non_zero_sign(*this) * pZero_double;
  629.   if (exponent < -1021)
  630.     denorm = 1;
  635.   if (!denorm) {
  636.     t_exp = exponent.tolong() + 1023 - 1;
  637.     s = sign(significant);
  639. /* significant's length has to be at least 53 digits */
  640.     if (significant.length() < 53)
  641.       t_sig = significant << (53 - significant.length());
  642.     else
  643.       t_sig = significant;
  645. /* Vorz positiv, weil integer-& auch negatives zur"uckgibt */
  646.     if (s == -1)
  647.       t_sig = integer(-1) * t_sig;
  649. /* remove leading one */
  650.     t_sig = b_round(t_sig, 53) & integer_52;
  651.   }
  654.   else {
  655.     t_exp = 0;
  657. /* significant's length has to be at least 52 digits */
  658.     if (significant.length() < 52)
  659.       t_sig = significant << (52 - significant.length());
  660.     else
  661.       t_sig = b_round(significant, 52);
  663.     t_sig = t_sig >> (-1022 - exponent.tolong());
  664.   }
  671.   double a;
  673.   long sign, h_sig, l_sig;
  675.   sign = (s == (-1));
  677.   h_sig = ((t_sig >> 32) & integer_20).tolong();
  678.   l_sig = (t_sig & integer_32).tolong();
  680.   a = compose_parts(sign, t_exp, h_sig, l_sig);
  682.   return a;
  685. }
  689. bigfloat bigfloat::round(long digits, rounding_modes mode)
  690. {
  691.   if (digits <= 0) {
  692.     (*this) = pZero;
  693.     return pZero;
  694.   }
  696.   rounding_modes m;
  697.   if (mode == DEFAULT)
  698.     m = round_mode;
  699.   else
  700.     m = mode;
  702.   if (digits < significant.length()) {
  703.     if (m == TO_NEAREST) {
  706.       integer temp = significant;
  708.       cut(temp, digits + 1);
  710. /* first case: the next digit is zero and therefor temp >> 1 is right */
  711.       if (abs(temp & 1) == 0) {
  712. significant = (temp >> 1);
  713.       }
  714.       else {
  715. /* the leftmost digit after the rounding digits is a one */
  716. /* check if there occur other ones in the following digits */
  718. long flag = test_for_one(significant, digits + 2, significant.length());
  720. /* if a one-digit occurred round up */
  721. if (flag) {
  722.   cut(significant, digits);
  723.   if (significant > 0)
  724.     significant += 1;
  725.   else
  726.     significant -= 1;
  728. }
  729. else {
  730.   cut(significant, digits);
  731.   if (abs(significant & 1) == 1) {
  732. /* same as above */
  734.     if (significant > 0)
  735.       significant += 1;
  736.     else
  737.       significant -= 1;
  739.   }
  740.   else
  741.     significant = temp;
  742. }
  743.       }
  745.       if (significant.length() > digits)
  746. cut(significant, digits);
  748. /* set new precision */
  749.       precision = digits;
  751.       return (*this);
  754.     }
  756.     if (m == TO_ZERO) {
  759.       cut(significant, digits);
  761. /* set new precision */
  762.       precision = digits;
  764.       return (*this);
  767.     }
  769.     if (m == TO_INF) {
  772.       if (significant > 0)
  773. return round(digits, TO_P_INF);
  774.       if (significant < 0)
  775. return round(digits, TO_N_INF);
  776.       return (*this);
  779.     }
  781.     if (m == TO_P_INF) {
  784.       long flag = test_for_one(significant, digits + 1, significant.length());
  786.       cut(significant, digits);
  787.       if ((significant > 0) && (flag)) {
  788. significant = significant + integer(1);
  789. if (significant.length() > digits)
  790.   cut(significant, digits);
  791.       }
  793. /* set new precision */
  794.       precision = digits;
  796.       return (*this);
  799.     }
  801.     if (m == TO_N_INF) {
  804.       long flag = test_for_one(significant, digits + 1, significant.length());
  806.       cut(significant, digits);
  807.       if ((significant < 0) && (flag)) {
  808. significant = significant - integer(1);
  809. if (significant.length() > digits)
  810.   cut(significant, digits);
  811.       }
  813. /* set new precision */
  814.       precision = digits;
  816.       return (*this);
  819.     }
  820.   }
  821.   return (*this);
  822. }
  826. bigfloat fabs(const bigfloat & b)
  827. {
  828.   if (isNaN(b))
  829.     return NaN;
  830.   if ((isnInf(b)) || (ispInf(b)))
  831.     return pInf;
  832.   if ((isnZero(b)) || (ispZero(b)))
  833.     return pZero;
  835.   return bigfloat(abs(b.significant), b.exponent);
  836. }
  840. integer ceil(const bigfloat & b)
  841. {
  842.   return b.tointeger(TO_P_INF);
  843. }
  847. integer floor(const bigfloat & b)
  848. {
  849.   return b.tointeger(TO_N_INF);
  850. }
  854. bigfloat pow2(const integer & p)
  855. {
  856.   return bigfloat(1, p + integer(1));
  857. }
  861. long sign(const bigfloat & b)
  862. {
  863.   return sign(b.significant);
  864. }
  868. long non_zero_sign(const bigfloat & b)
  869. {
  870.   if ((ispInf(b)) || (ispZero(b)))
  871.     return 1;
  872.   if ((isnInf(b)) || (isnZero(b)))
  873.     return -1;
  874.   return sign(b);
  875. }
  879. void bigfloat::print_contents(void) const
  880. {
  881.   integer s = significant, e = exponent;
  882.   cout << "Significant: ";
  883.   if (s < 0) {
  884.     cout << "-";
  885.     s = s * integer(-1);
  886.   }
  887.   binout(cout, s);
  888.   cout << "n";
  889.   cout << "Exponent: ";
  890.   if (e < 0) {
  891.     cout << "-";
  892.     e = e * integer(-1);
  893.   }
  894.   binout(cout, e);
  895.   cout << "n";
  896. }
  901. bigfloat add(const bigfloat & aa, const bigfloat & bb, long prec, rounding_modes mode)
  902. {
  905.   if ((isNaN(aa)) || (isNaN(bb)))
  906.     return NaN;
  908.   if ((ispInf(aa)) || (isnInf(aa))) {
  909.     if ((bb != pInf) && (bb != nInf))
  910.       return aa;
  911.     return NaN;
  912.   }
  914.   if ((ispInf(bb)) || (isnInf(bb))) {
  915.     if ((aa != pInf) && (aa != nInf))
  916.       return bb;
  917.     return NaN;
  918.   }
  920.   if (sign(aa) == 0)
  921.     return bb;
  922.   if (sign(bb) == 0)
  923.     return aa;
  927.   bigfloat a = aa, b = bb, c;
  928.   if (prec == -1)
  929.     prec = global_prec;
  931. /* We firstly have to adjust the significants' precisions to
  932.      guaranty an exact result */
  934.   if (a.precision > b.precision)
  935.     b.set_precision(a.precision);
  936.   if (b.precision > a.precision)
  937.     a.set_precision(b.precision);
  939. /* now we have to branch into 3 main cases */
  941.   if (a.exponent == b.exponent) {
  944.     long temp1, temp2;
  946.     c.significant = a.significant + b.significant;
  948.     temp1 = c.significant.length();
  949.     temp2 = a.significant.length();
  950.     c.exponent = a.exponent + integer(temp1 - temp2);
  951.     c.precision = c.significant.length();
  953.     c.set_precision(prec, mode);
  955.     if (c.significant == 0)
  956.       c = pZero;
  958.     return c;
  961.   }
  962.   if (a.exponent > b.exponent) {
  965.     integer temp = a.exponent - b.exponent;
  966.     a.significant = (a.significant << temp.tolong());
  967.     c.significant = a.significant + b.significant;
  969.     c.exponent = a.exponent;
  970.     if (c.significant.length() != a.significant.length())
  971.       c.exponent += (c.significant.length() - a.significant.length());
  973.     c.precision = c.significant.length();
  974.     c.set_precision(prec, mode);
  976.     if (c.significant == 0)
  977.       c = pZero;
  979.     return c;
  982.   }
  986. /* siehe oben */
  987.   integer temp = b.exponent - a.exponent;
  988.   b.significant = (b.significant << temp.tolong());
  989.   c.significant = a.significant + b.significant;
  991.   c.exponent = b.exponent;
  993.   if (c.significant.length() != b.significant.length())
  994.     c.exponent += (c.significant.length() - b.significant.length());
  996.   c.precision = c.significant.length();
  997.   c.set_precision(prec, mode);
  999.   if (c.significant == 0)
  1000.     c = pZero;
  1002.   return c;
  1005. }
  1009. bigfloat sub(const bigfloat & a, const bigfloat & b, long prec, rounding_modes mode)
  1010. {
  1011.   bigfloat bb(-b.significant, b.exponent);
  1012.   return add(a, bb, prec, mode);
  1013. }
  1017. bigfloat mul(const bigfloat & aa, const bigfloat & bb, long prec, rounding_modes mode)
  1018. {
  1019.   if (prec == -1)
  1020.     prec = global_prec;
  1022.   if ((isNaN(aa)) || (isNaN(bb)))
  1023.     return NaN;
  1025.   if (((ispZero(aa)) || (isnZero(aa))) && (!isInf(bb))) {
  1026.     if (bb > 0)
  1027.       return aa;
  1028.     if (ispZero(aa))
  1029.       return nZero;
  1030.     return pZero;
  1031.   }
  1032.   if (((ispZero(bb)) || (isnZero(bb))) && (!isInf(aa))) {
  1033.     if (aa > 0)
  1034.       return bb;
  1035.     if (ispZero(bb))
  1036.       return nZero;
  1037.     return pZero;
  1038.   }
  1040.   if (ispInf(aa)) {
  1041.     if (bb > 0)
  1042.       return pInf;
  1043.     if (bb < 0)
  1044.       return nInf;
  1045.     return NaN;
  1046.   }
  1047.   if (isnInf(aa)) {
  1048.     if (bb > 0)
  1049.       return nInf;
  1050.     if (bb < 0)
  1051.       return pInf;
  1052.     return NaN;
  1053.   }
  1054.   if (ispInf(bb)) {
  1055.     if (aa > 0)
  1056.       return pInf;
  1057.     if (aa < 0)
  1058.       return nInf;
  1059.     return NaN;
  1060.   }
  1061.   if (isnInf(bb)) {
  1062.     if (aa > 0)
  1063.       return nInf;
  1064.     if (aa < 0)
  1065.       return pInf;
  1066.     return NaN;
  1067.   }
  1071.   bigfloat c;
  1073.   c.significant = aa.significant * bb.significant;
  1074.   c.exponent = aa.exponent + bb.exponent;
  1075.   c.exponent -= (aa.significant.length() + bb.significant.length()) - c.significant.length();
  1077.   c.precision = c.significant.length();
  1078.   c.set_precision(prec, mode);
  1080.   return c;
  1083. }
  1087. bigfloat div(const bigfloat & aa, const bigfloat & bb, long prec, rounding_modes mode)
  1088. {
  1089. /* special checking for zero-cases */
  1090.   if (prec == -1)
  1091.     prec = global_prec;
  1095.   if ((isNaN(aa)) || (isNaN(bb)))
  1096.     return NaN;
  1097.   if (((isZero(aa)) && (isZero(bb))) || ((isInf(aa)) && (isInf(bb))))
  1098.     return NaN;
  1099.   if ((isInf(aa)) || (isZero(bb)))
  1100.     return (non_zero_sign(aa) * non_zero_sign(bb) * pInf);
  1101.   if ((isZero(aa)) || (isInf(bb)))
  1102.     return (non_zero_sign(aa) * non_zero_sign(bb) * pZero);
  1109.   bigfloat a = aa, c;
  1110.   integer temp;
  1112.   long d = prec + bb.significant.length() - a.significant.length() + 2;
  1114.   if (d > 0)
  1115.     a.significant = a.significant << d;;
  1120.   c.significant = a.significant / bb.significant;
  1121.   c.exponent = a.exponent - bb.exponent;
  1122.   c.precision = prec;
  1124. /* $delta$-Test */
  1125.   if (bcomp(fabs(a).significant, fabs(bb).significant) == 1)
  1126.     c.exponent += 1;
  1128. /* cut the significant to prec+1 digits */
  1129.   cut(c.significant, prec + 1);
  1133.   if (mode == DEFAULT)
  1134.     mode = round_mode;
  1135.   switch (mode) {
  1136.       case TO_NEAREST:
  1139. if (abs(c.significant & 1) == 0) {
  1140.   c.significant = (c.significant >> 1);
  1141.   return c;
  1142. }
  1143. else {
  1144. /* now we have to check whether the remainder is zero */
  1146.   bigfloat temp = mul(c, bb, 1, EXACT);
  1147.   temp = mul(c, bb, 1, EXACT);
  1149.   if (temp == aa)
  1150.     return c.round(prec, TO_NEAREST);
  1151.   return c.round(prec, TO_INF);
  1152. }
  1155.       case TO_ZERO:
  1158. c.significant = c.significant >> 1;
  1159. return c;
  1162.       case TO_P_INF:
  1165. if (mul(c, bb, 1, EXACT) == aa)
  1166.   return c.round(prec, TO_ZERO);
  1167. if (c.significant > 0)
  1168.   return c.round(prec, TO_INF);
  1169. if (c.significant < 0)
  1170.   return c.round(prec, TO_ZERO);
  1171. return c;
  1174.       case TO_N_INF:
  1177. if (mul(c, bb, 1, EXACT) == aa)
  1178.   return c.round(prec, TO_ZERO);
  1179. if (c.significant > 0)
  1180.   return c.round(prec, TO_ZERO);
  1181. if (c.significant < 0)
  1182.   return c.round(prec, TO_INF);
  1183. return c;
  1186.       case TO_INF:
  1189. if (mul(c, bb, 1, EXACT) == aa)
  1190.   return c.round(prec, TO_ZERO);
  1191. if (c.significant > 0)
  1192.   return c.round(prec, TO_P_INF);
  1193. if (c.significant < 0)
  1194.   return c.round(prec, TO_N_INF);
  1195. return c;
  1198.       case EXACT:
  1199. error_handler(1, "div(bigfloat) : not implemented case EXACT!");
  1200.       case DEFAULT:;
  1201.   }
  1208.   if (mode == DEFAULT)
  1209.     mode = round_mode;
  1210.   switch (mode) {
  1211.       case TO_NEAREST:
  1214. if (abs(c.significant & 1) == 0) {
  1215.   c.significant = (c.significant >> 1);
  1216.   return c;
  1217. }
  1218. else {
  1219. /* now we have to check whether the remainder is zero */
  1221.   bigfloat temp = mul(c, bb, 1, EXACT);
  1222.   temp = mul(c, bb, 1, EXACT);
  1224.   if (temp == aa)
  1225.     return c.round(prec, TO_NEAREST);
  1226.   return c.round(prec, TO_INF);
  1227. }
  1230.       case TO_ZERO:
  1233. c.significant = c.significant >> 1;
  1234. return c;
  1237.       case TO_P_INF:
  1240. if (mul(c, bb, 1, EXACT) == aa)
  1241.   return c.round(prec, TO_ZERO);
  1242. if (c.significant > 0)
  1243.   return c.round(prec, TO_INF);
  1244. if (c.significant < 0)
  1245.   return c.round(prec, TO_ZERO);
  1246. return c;
  1249.       case TO_N_INF:
  1252. if (mul(c, bb, 1, EXACT) == aa)
  1253.   return c.round(prec, TO_ZERO);
  1254. if (c.significant > 0)
  1255.   return c.round(prec, TO_ZERO);
  1256. if (c.significant < 0)
  1257.   return c.round(prec, TO_INF);
  1258. return c;
  1261.       case TO_INF:
  1264. if (mul(c, bb, 1, EXACT) == aa)
  1265.   return c.round(prec, TO_ZERO);
  1266. if (c.significant > 0)
  1267.   return c.round(prec, TO_P_INF);
  1268. if (c.significant < 0)
  1269.   return c.round(prec, TO_N_INF);
  1270. return c;
  1273.       case EXACT:
  1274. error_handler(1, "div(bigfloat) : not implemented case EXACT!");
  1275.       case DEFAULT:;
  1276.   }
  1279.   return 0; // never reached
  1281. }
  1285. bigfloat sqrt(const bigfloat & a, long prec, rounding_modes mode)
  1286. {
  1287.   bigfloat ne, it;
  1288.   long t_prec;
  1292.   if (isNaN(a))
  1293.     return NaN;
  1294.   if ((isnZero(a)) || (ispZero(a)))
  1295.     return a;
  1296.   if (ispInf(a))
  1297.     return pInf;
  1298.   if (isnInf(a))
  1299.     return NaN;
  1300.   if (a < 0)
  1301.     return NaN;
  1305.   if (prec == -1)
  1306.     prec = global_prec;
  1310.   double d = bigfloat(a.significant, 0).todouble();
  1311.   if (!(a.exponent % 2 == 0))
  1312.     d *= 2;
  1314.   ne = sqrt(d);
  1315.   ne.exponent = ne.exponent + (a.exponent / integer(2));
  1316.   t_prec = 52;
  1318. /* cout Startwert */
  1322.   int flag = 0;
  1323.   do {
  1324.     it = ne;
  1325.     ne = add(mul(it, it, 1, EXACT), a, 1, EXACT); /* compute numerator of
  1326.  * (1) */
  1327.     ne = div(ne, mul(2, it, 1, EXACT), t_prec, TO_P_INF); /* compute next
  1328.  * iteration value */
  1329.     t_prec *= 2;
  1330.   } while (t_prec < prec);
  1332.   ne.round(prec, TO_P_INF);
  1333.   it.round(prec, TO_P_INF);
  1335.   while (ne != it) {
  1336.     it = ne;
  1337.     ne = add(mul(it, it, 1, EXACT), a, 1, EXACT); /* compute numerator of
  1338.  * (1) */
  1339.     ne = div(ne, mul(2, it, 1, EXACT), prec, TO_P_INF); /* compute next
  1340.  * iteration value */
  1341.   }
  1345.   bigfloat c(ne.significant - integer(1), ne.exponent);
  1346.   if ((c * c) >= a)
  1347.     it = c;
  1352.   bigfloat mid;
  1353.   if (mode == DEFAULT)
  1354.     mode = round_mode;
  1356.   switch (mode) {
  1357.       case TO_P_INF:
  1358. return it;
  1359.       case TO_ZERO:
  1360.       case TO_N_INF:
  1361. if ((it * it) != a)
  1362.   return bigfloat(ne.significant - integer(1), ne.exponent);
  1363. else
  1364.   return it;
  1365.       case TO_NEAREST:
  1366. mid = bigfloat((it.significant >> 1) + integer(1), it.exponent);
  1367. /* $mid=frac{it+(it-Delta)}{2}$ */
  1368. if ((mid * mid) < a) {
  1369.   mid.round(prec, TO_ZERO);
  1370.   return mid;
  1371. }
  1372. else {
  1373.   if ((mid * mid) != a)
  1374.     return it;
  1375.   else {
  1376.     mid.round(prec, TO_NEAREST);
  1377.     return mid;
  1378.   }
  1379. }
  1380.       case EXACT:
  1381. error_handler(1, "bigfloat::sqrt : the EXACT-rounding mode is not provided!");
  1383.   }
  1385.   return 0; // never reached
  1387. }
  1392. bigfloat & bigfloat::operator = (const bigfloat & b) {
  1393.   significant = b.significant;
  1394.   exponent = b.exponent;
  1395.   precision = b.precision;
  1396.   return (*this);
  1397. }
  1402. ostream & operator << (ostream & os, const bigfloat & b) {
  1403.   if (1) { /* if isSpecial */
  1404.     if (isNaN(b))
  1405.       return os << "NaN";
  1406.     if (ispInf(b))
  1407.       return os << "+Inf";
  1408.     if (isnInf(b))
  1409.       return os << "-Inf";
  1410.     if (ispZero(b))
  1411.       return os << "+0";
  1412.     if (isnZero(b))
  1413.       return os << "-0";
  1414.   }
  1415.   int sign_b = sign(b.significant);
  1416.   if (output == BIN_OUT) {
  1417.     if (sign_b < 0)
  1418.       os << "-";
  1419.     os << "0.";
  1420.     if (sign_b >= 0)
  1421.       binout(os, b.significant);
  1422.     else
  1423.       binout(os, -b.significant);
  1424.     os << "E";
  1425.     if (b.exponent < 0)
  1426.       os << "-";
  1427.     else
  1428.       os << "+";
  1429.     binout(os, b.exponent);
  1430.   }
  1431.   if (output == DEC_OUT)
  1432.     decimal_output(os, b);
  1433.   return os;
  1434. }
  1435. istream & operator >> (istream & is, bigfloat & b) {
  1436.   char temp[bin_maxlen]; /* z.Z. nur binaere-Exponentialeingabe
  1437.  * moeglich */
  1438.   is >> temp;
  1439.   b = outofchar(temp);
  1440.   return is;
  1441. }
  1446. bool operator == (const bigfloat & a, const bigfloat & b) {
  1447. /* check whether NaN or infinity cases appear */
  1448.   if ((isNaN(a)) || (isNaN(b)) || (ispInf(a)) || (ispInf(b))
  1449.     || (isnInf(a)) || (isnInf(b)))
  1450.     error_handler(1, "bigfloat::operator == : NaN or Infinity case occurred!");
  1452.   if ((a.significant == 0) && (b.significant == 0))
  1453.     return 1;
  1454.   if (sign(a.significant) != sign(b.significant))
  1455.     return 0;
  1456.   if (a.exponent == b.exponent)
  1457.     return (bcomp(a.significant, b.significant) == 0);
  1458.   else
  1459.     return 0;
  1460. }
  1462. bool operator != (const bigfloat & a, const bigfloat & b) {
  1463.   return !(a == b);
  1464. }
  1466. bool operator > (const bigfloat & a, const bigfloat & b) {
  1467.   if ((isNaN(a)) || (isNaN(b)))
  1468.     error_handler(1, "bigfloat::operator > : NaN case occurred!");
  1469.   if (((ispInf(a)) && (ispInf(b))) || ((isnInf(a)) && (isnInf(b))))
  1470.     error_handler(1, "bigfloat::operator > : both operands +(-)infinity");
  1471.   if (ispInf(a))
  1472.     return 1;
  1473.   if (isnInf(a))
  1474.     return 0;
  1475.   if (ispInf(b))
  1476.     return 0;
  1477.   if (isnInf(b))
  1478.     return 1;
  1480.   if ((a.significant == 0) && (b.significant == 0))
  1481.     return 0;
  1482.   if (sign(a.significant) != sign(b.significant)) {
  1483.     if (sign(a.significant) == 0)
  1484.       if (sign(b.significant) < 0)
  1485. return 1;
  1486.       else
  1487. return 0;
  1488.     if (sign(b.significant) == 0)
  1489.       if (sign(a.significant) > 0)
  1490. return 1;
  1491.       else
  1492. return 0;
  1494.     if (a.significant > 0)
  1495.       return 1;
  1496.     return 0;
  1497.   }
  1498.   if (a.exponent != b.exponent)
  1499.     return (a.exponent > b.exponent);
  1500.   else
  1501.     return (bcomp(a.significant, b.significant) == 1);
  1502. }
  1504. bool operator >= (const bigfloat & a, const bigfloat & b) {
  1505.   return ((a > b) || (a == b));
  1506. }
  1508. bool operator < (const bigfloat & a, const bigfloat & b) {
  1509.   return (!(a >= b));
  1510. }
  1512. bool operator <= (const bigfloat & a, const bigfloat & b) {
  1513.   return (!(a > b));
  1514. }
  1518. bigfloat operator + (const bigfloat & a, const bigfloat & b) {
  1519.   return add(a, b);
  1520. }
  1522. bigfloat operator - (const bigfloat & a, const bigfloat & b) {
  1523.   return sub(a, b);
  1524. }
  1526. bigfloat operator *(const bigfloat & a, const bigfloat & b) {
  1527.   return mul(a, b);
  1528. }
  1530. bigfloat operator / (const bigfloat & a, const bigfloat & b) {
  1531.   return div(a, b);
  1532. }
  1534. bigfloat operator - (const bigfloat & a) {
  1535.   return bigfloat(-a.significant, a.exponent);
  1536. }
  1541. bool isNaN(const bigfloat & b)
  1542. {
  1543.   if ((b.significant == 0) && (b.exponent == 1026))
  1544.     return true;
  1545.   return false;
  1546. }
  1548. bool isnInf(const bigfloat & b)
  1549. {
  1550.   if ((b.significant == 0) && (b.exponent == 1025))
  1551.     return true;
  1552.   return false;
  1553. }
  1555. bool ispInf(const bigfloat & b)
  1556. {
  1557.   if ((b.significant == 0) && (b.exponent == 1024))
  1558.     return true;
  1559.   return false;
  1560. }
  1562. bool isnZero(const bigfloat & b)
  1563. {
  1564.   if ((b.significant == 0) && (b.exponent == (-1024)))
  1565.     return true;
  1566.   return false;
  1567. }
  1569. bool ispZero(const bigfloat & b)
  1570. {
  1571.   if ((b.significant == 0) && (b.exponent == (-1023)))
  1572.     return true;
  1573.   return false;
  1574. }
  1576. bool isInf(const bigfloat & b)
  1577. {
  1578.   return (ispInf(b) || isnInf(b));
  1579. }
  1581. bool isZero(const bigfloat & b)
  1582. {
  1583.   return (ispZero(b) || isnZero(b));
  1584. }