1
2 /*
3 * Real - Implements a real number as either a floating point or as a ratio of arbitrary
4 * sized integers
5 * Format - Stores information about how to output a number
6 * Copyright (c) 2003-2006 by Mattias Hultgren <mattias_hultgren@tele2.se>
7 *
8 * See real.h
9 */
10
11 /*
12 News
13 ----
14
15 v4 2006-07-17 - 2006-07-30 & 2006-10-05
16 --
17
18 Replaced all throw_... calls with THROW_ERROR
19 Replaced all std::string usage with utf8_string
20 Fixed and output error on amd64 (number containing potensis would display incorrect)
21
22 v3 2006-05-07 - 2006-05-21
23 --
24
25 Some small changes in operator==, operator<=, operator>=
26 Replaced constructors using float, double and long double with a constructor using floatx
27 Fixed a bug in operator< and operator> which always converted ratio/integers to floatx before the comparing
28 which lead to that only the most significant numbers where compared. E.g. (2^300-1) > (2^300-10) became false
29 Small clean up in floatx_to_string & Real::append_to_string
30
31 v2 2006-04-03
32 --
33
34 Changed so that Real::top points to a Integer[2] and bottom to top[1], instead of them having one each
35 ( this got ump ~20% faster when dealing with small integers (integers below 2^32) )
36
37 v1 2005-09-25 - 2005-10-22
38 --
39
40 Extracted class Real from the math2.* system
41
42 */
43
44 #include "real.h"
45 #include "keyfile.h"
46
47
48 // Format-object start
make_format_legal(void)49 void Format::make_format_legal(void)
50 {
51 if(max_decimals > 50)
52 max_decimals = 50;
53
54 if(max_decimals != 0)
55 {
56 floatx newlow = powx( 10.0, -floatx(max_decimals) );
57
58 mixed_lowlimit = (newlow > mixed_lowlimit) ? newlow : mixed_lowlimit;
59 }
60 else
61 mixed_lowlimit = 1.0;
62
63 if(mixed_lowlimit > mixed_hilimit)
64 mixed_hilimit = mixed_lowlimit;
65 }
66
Format()67 Format::Format()
68 {
69 mixed_hilimit = 1e9;
70 mixed_lowlimit = 1e-6;
71
72 type = FormatType_Mixed;
73 max_decimals = 5;
74 show_trailingzeros = false;
75 fraction_type = FormatFractionType_Integer_plus_fraction;
76 make_format_legal();
77 }
78
set_type(FormatType newtype)79 void Format::set_type(FormatType newtype)
80 {
81 if( newtype == FormatType_Normal || newtype == FormatType_Mixed || newtype == FormatType_Science )
82 type = newtype;
83 }
84
set_fraction_type(FormatFractionType newtype)85 void Format::set_fraction_type(FormatFractionType newtype)
86 {
87 if( newtype == FormatFractionType_Fractions || newtype == FormatFractionType_Integer_plus_fraction ||
88 newtype == FormatFractionType_Decimals )
89 fraction_type = newtype;
90 }
91
set_max_decimals(uint32 n)92 void Format::set_max_decimals(uint32 n)
93 {
94 max_decimals = n;
95 make_format_legal();
96 }
97
set_mixed_limits(floatx lowlimit,floatx hilimit)98 void Format::set_mixed_limits(floatx lowlimit, floatx hilimit)
99 {
100 if(lowlimit < hilimit)
101 {
102 mixed_lowlimit = lowlimit;
103 mixed_hilimit = hilimit;
104 make_format_legal();
105 }
106 }
107 // format-object end
floatx_to_string(floatx var,utf8_string & str,const Format & fmt)108 void floatx_to_string(floatx var, utf8_string &str, const Format &fmt) throw(error_obj)
109 {
110 int32 power = 0, extra_precision=0;
111 floatx tmpvar;
112
113 tmpvar = var;
114 if( (fmt.get_type() == FormatType_Science) |
115 ( (fmt.get_type() == FormatType_Mixed) & ((tmpvar >= fmt.get_mixed_hilimit()) | (tmpvar <= -fmt.get_mixed_hilimit()) | ((tmpvar < fmt.get_mixed_lowlimit()) & (tmpvar > -fmt.get_mixed_lowlimit()) & (tmpvar != 0.0))) )
116 )
117 {
118 if( var >= 10.0 || var <= -10.0 )
119 {
120 for( ; var >= 10.0 || var <= -10.0; power++, var/=10.0 ) // calculating the 'power'
121 {
122 if( power > 30000 ) // a long double can only hold a power of about 4931, thus this is inf
123 {
124 if(var > 0.0)
125 str.append( "inf" );
126 else
127 str.append( "-inf" );
128
129 return;
130 }
131 }
132 }
133 else if( (var < 1.0) & (var > -1.0) & (var != 0.0) )
134 {
135 for( ; (var < 1.0) & (var > -1.0); power--, var*=10.0 ) // calculating the 'power'
136 ;
137 }
138 }
139 else if( var < 1.0 && var > -1.0 && var != 0.0 )
140 {
141 for( var*=10.0; var<1.0 && var>-1.0; extra_precision++, var*=10.0 ) // calculating the 'power'
142 ;
143 var = tmpvar;
144 }
145
146 {
147 char out_field[100], syntax[20];
148 signed int old_size = str.get_length()-1;
149
150 uint64_to_char( uint64(fmt.get_max_decimals() + extra_precision), out_field );
151 #ifdef USE_LONG_DOUBLE
152 snprintf( syntax, 20, "%%.%sLf", out_field );
153 snprintf( out_field, 100, syntax, var );
154 #else
155 snprintf( syntax, 20, "%%.%sf", out_field );
156 snprintf( out_field, 100, syntax, double(var) );
157 #endif
158
159 str.append( out_field );
160
161 if( fmt.get_max_decimals() != 0 && fmt.show_trailingzeros == false )
162 {
163 for( signed int i=str.get_length()-1; i>=old_size; i-- ) // removes trailing zeros
164 {
165 if( str.test_character( i, "0" ) )
166 str.crop( i );
167 else if( str.test_character( i, "." ) )
168 {
169 str.crop( i );
170 break;
171 }
172 else
173 break;
174 }
175 }
176
177
178 if( (fmt.get_type() == FormatType_Science) |
179 ( (fmt.get_type() == FormatType_Mixed) & ((tmpvar >= fmt.get_mixed_hilimit()) | (tmpvar <= -fmt.get_mixed_hilimit()) | ((tmpvar < fmt.get_mixed_lowlimit()) & (tmpvar > -fmt.get_mixed_lowlimit()) & (tmpvar != 0.0))) )
180 )
181 {
182 char tmp[10];
183 if( power >= 0 )
184 {
185 uint64_to_char( uint64(power), tmp );
186 snprintf( out_field, 100, "*10^%s", tmp );
187 }
188 else
189 {
190 uint64_to_char( uint64(-power), tmp );
191 snprintf( out_field, 100, "*10^-%s", tmp );
192 }
193 str.append( out_field );
194 }
195 }
196 }
197
198 // class Real start
199
append_to_string(utf8_string & str,const Format & fmt) const200 void Real::append_to_string( utf8_string &str, const Format &fmt ) const throw(error_obj)
201 {
202 try
203 {
204 if( exact )
205 {
206 if( fmt.get_fraction_type() == FormatFractionType_Decimals )
207 floatx_to_string( floatx(*top)/floatx(*bottom), str, fmt);
208 else
209 {
210 if((top->compare_abs_value(*bottom) ) != 0 && fmt.get_fraction_type() == FormatFractionType_Integer_plus_fraction)
211 {
212 Integer tmp_int, tmp_part;
213 bool abs_over_one;
214
215 top->divide( *bottom, &tmp_int, &tmp_part );
216 if( abs_over_one = (top->compare_abs_value( *bottom ) >= 0) || *top == 0 )
217 tmp_int.append_to_string( str );
218
219 if(*bottom != 1)
220 {
221 if( *top > 0 )
222 {
223 if( abs_over_one )
224 str.append( " + " );
225 tmp_part.append_to_string( str );
226 }
227 else
228 {
229 if( abs_over_one )
230 {
231 str.append( " - " );
232 tmp_part.negate();
233 tmp_part.append_to_string( str );
234 }
235 else
236 tmp_part.append_to_string( str );
237 }
238 str.append( "/" );
239 bottom->append_to_string( str );
240 }
241 }
242 else
243 {
244 top->append_to_string( str );
245 if(*bottom != 1 && *top != 0)
246 {
247 str.append( "/" );
248 bottom->append_to_string( str );
249 }
250 }
251 }
252 }
253 else
254 floatx_to_string( value, str, fmt );
255 }
256 catch(error_obj error) { throw error; }
257 catch(...) { THROW_ERROR( ErrorType_Memory, _("Couldn't get memory.") ); }
258 }
259
260
Real(int32 new_value)261 Real::Real(int32 new_value) throw(error_obj)
262 {
263 Integer *ptr = 0;
264
265 try
266 {
267 ptr = new Integer [2];
268
269 ptr[0] = new_value;
270 ptr[1] = 1;
271 }
272 catch(...)
273 {
274 delete [] ptr;
275 THROW_ERROR( ErrorType_Memory, _("Couldn't get memory.") );
276 }
277 exact = true;
278 top = ptr;
279 bottom = &ptr[1];
280 }
281
Real(int64 new_value)282 Real::Real(int64 new_value) throw(error_obj)
283 {
284 Integer *ptr = 0;
285
286 try
287 {
288 ptr = new Integer [2];
289
290 ptr[0] = new_value;
291 ptr[1] = 1;
292 }
293 catch(...)
294 {
295 delete [] ptr;
296 THROW_ERROR( ErrorType_Memory, _("Couldn't get memory.") );
297 }
298 exact = true;
299 top = ptr;
300 bottom = &ptr[1];
301 }
302
Real(const Integer & new_value)303 Real::Real(const Integer &new_value) throw(error_obj)
304 {
305 Integer *ptr = 0;
306
307 try
308 {
309 ptr = new Integer [2];
310
311 ptr[0] = new_value;
312 ptr[1] = 1;
313 }
314 catch(...)
315 {
316 delete [] ptr;
317 THROW_ERROR( ErrorType_Memory, _("Couldn't get memory.") );
318 }
319 exact = true;
320 top = ptr;
321 bottom = &ptr[1];
322 }
323
Real(const Integer & new_top,const Integer & new_bottom)324 Real::Real(const Integer &new_top, const Integer &new_bottom) throw(error_obj)
325 {
326 Integer *ptr = 0;
327
328 try
329 {
330 ptr = new Integer [2];
331
332 ptr[0] = new_top;
333 ptr[1] = new_bottom;
334 }
335 catch(...)
336 {
337 delete [] ptr;
338 THROW_ERROR( ErrorType_Memory, _("Couldn't get memory.") );
339 }
340 exact = true;
341 top = ptr;
342 bottom = &ptr[1];
343
344 legalice();
345 }
346
Real(const Real & val)347 Real::Real(const Real &val) throw(error_obj)
348 {
349 if( val.exact )
350 {
351 Integer *ptr = 0;
352
353 try
354 {
355 ptr = new Integer [2];
356 }
357 catch(...)
358 {
359 delete [] ptr;
360 THROW_ERROR( ErrorType_Memory, _("Couldn't get memory.") );
361 }
362
363 exact = true;
364 top = ptr;
365 bottom = &ptr[1];
366
367 *top = *val.top;
368 *bottom = *val.bottom;
369 }
370 else
371 {
372 exact = false;
373 value = val.value;
374 }
375 }
376
~Real()377 Real::~Real()
378 {
379 if( exact )
380 delete [] top;
381 }
382
legalice(void)383 void Real::legalice(void)
384 {
385 if( exact )
386 {
387 if( *bottom < 0 )
388 {
389 top->negate();
390 bottom->negate();
391 }
392 else if( *bottom == 0 )
393 {
394 *bottom = 1;
395 *top = 0;
396 return;
397 }
398
399 if( *bottom == 1 )
400 return;
401
402 if( *top == 0 )
403 *bottom = 1;
404 else
405 {
406 Integer div = gcd( *top, *bottom ); // lets simplify the fraction
407 if( div != 0 )
408 {
409 Integer tmp;
410
411 top->divide( div, &tmp, 0 );
412 top->swap( tmp );
413 bottom->divide( div, &tmp, 0 );
414 bottom->swap( tmp );
415 }
416 }
417 }
418 }
419
operator +(const Real & val) const420 Real Real::operator+(const Real &val) const throw(error_obj)
421 {
422 if( exact == false || val.exact == false ) // floating pointcalculus
423 return Real( floatx(*this) + floatx(val) );
424 else // exact fraction calculus
425 {
426 if( *bottom == *val.bottom )
427 {
428 Real res( *top + *val.top, *bottom );
429 res.legalice();
430 return res;
431 }
432 else
433 {
434 Real res( *top * *val.bottom + *val.top * *bottom, *bottom * *val.bottom );
435 res.legalice();
436 return res;
437 }
438 }
439 }
add(const Real & val)440 void Real::add(const Real &val) throw(error_obj)
441 {
442 if( exact == false || val.exact == false ) // floating pointcalculus
443 {
444 if( exact )
445 {
446 floatx tmp;
447 tmp = floatx(*this);
448 delete [] top;
449 exact = false;
450 value = tmp;
451 }
452 value += floatx(val);
453 }
454 else // exact fraction calculus
455 {
456 if( *bottom == *val.bottom )
457 top->add( *val.top );
458 else
459 {
460 Integer tmp;
461 top->mul( *val.bottom, &tmp );
462 val.top->mul( *bottom, top );
463 top->add( tmp );
464
465 bottom->mul( *val.bottom, &tmp );
466 bottom->swap( tmp );
467 }
468 legalice();
469 }
470 }
471
operator -(const Real & val) const472 Real Real::operator-(const Real &val) const throw(error_obj)
473 {
474 if( exact == false || val.exact == false ) // floating pointcalculus
475 return Real( floatx(*this) - floatx(val) );
476 else // exact fraction calculus
477 {
478 if( *bottom == *val.bottom )
479 {
480 Real res( *top - *val.top, *bottom );
481 res.legalice();
482 return res;
483 }
484 else
485 {
486 Integer a,b;
487 top->mul( *val.bottom, &a );
488 val.top->mul( *bottom, &b );
489 a -= b;
490 bottom->mul( *val.bottom, &b );
491 {
492 Real res( a, b );
493 res.legalice();
494 return res;
495 }
496 }
497 }
498 }
sub(const Real & val)499 void Real::sub(const Real &val) throw(error_obj)
500 {
501 if( exact == false || val.exact == false ) // floating pointcalculus
502 {
503 if( exact )
504 {
505 floatx tmp;
506 tmp = floatx(*this);
507 delete [] top;
508 exact = false;
509 value = tmp;
510 }
511 value -= floatx(val);
512 }
513 else // exact fraction calculus
514 {
515 if( *bottom == *val.bottom )
516 top->sub( *val.top );
517 else
518 {
519 Integer tmp;
520 top->mul( *val.bottom, &tmp );
521 val.top->mul( *bottom, top );
522 top->negate();
523 top->add( tmp );
524
525 bottom->mul( *val.bottom, &tmp );
526 bottom->swap( tmp );
527 }
528 legalice();
529 }
530 }
531
operator -(const Real & val)532 Real operator-(const Real &val) throw(error_obj)
533 {
534 if( val.exact )
535 {
536 if( *val.top == 0 )
537 return val;
538 return Real( -*val.top, *val.bottom );
539 }
540 else
541 {
542 if( val.value == 0.0 )
543 return val;
544 return Real( -val.value );
545 }
546 }
547
operator *(const Real & val) const548 Real Real::operator*(const Real &val) const throw(error_obj)
549 {
550 if(exact == false || val.exact == false) // floating point calculus
551 return Real( floatx(*this) * floatx(val) );
552 else // exact fraction calculus
553 {
554 Integer a,b;
555 top->mul( *val.top, &a );
556 bottom->mul( *val.bottom, &b );
557 {
558 Real res( a, b );
559 res.legalice();
560
561 return res;
562 }
563 }
564 }
mul(const Real & val)565 void Real::mul(const Real &val) throw(error_obj)
566 {
567 if( exact == false || val.exact == false ) // floating pointcalculus
568 {
569 if( exact )
570 {
571 floatx tmp;
572 tmp = floatx(*this);
573 delete [] top;
574 exact = false;
575 value = tmp;
576 }
577 value *= floatx(val);
578 }
579 else // exact fraction calculus
580 {
581 Integer tmp;
582 top->mul( *val.top, &tmp );
583 top->swap( tmp );
584 bottom->mul( *val.bottom, &tmp );
585 bottom->swap( tmp );
586 legalice();
587 }
588 }
589
operator /(const Real & val) const590 Real Real::operator/(const Real &val) const throw(error_obj)
591 {
592 if( exact == false || val.exact == false ) // floating pointcalculus
593 return Real( floatx(*this) / floatx(val) );
594 else // exact fraction calculus
595 {
596 Real res( *top * *val.bottom, *bottom * *val.top );
597 res.legalice();
598 return res;
599 }
600 }
div(const Real & val)601 void Real::div(const Real &val) throw(error_obj)
602 {
603 if( exact == false || val.exact == false ) // floating point calculus
604 {
605 if( exact )
606 {
607 floatx tmp;
608 tmp = floatx(*this);
609 delete [] top;
610 exact = false;
611 value = tmp;
612 }
613 value /= floatx(val);
614 }
615 else // exact fraction calculus
616 {
617 Integer tmp;
618 top->mul( *val.bottom, &tmp );
619 top->swap( tmp );
620 bottom->mul( *val.top, &tmp );
621 bottom->swap( tmp );
622 legalice();
623 }
624 }
625
operator =(const Real & val)626 Real & Real::operator=(const Real &val) throw(error_obj)
627 {
628 if( val.exact )
629 {
630 if( !exact )
631 {
632 Integer *ptr = 0;
633
634 try
635 {
636 ptr = new Integer [2];
637 }
638 catch(...)
639 {
640 delete [] ptr;
641 THROW_ERROR( ErrorType_Memory, _("Couldn't get memory.") );
642 }
643
644 exact = true;
645 top = ptr;
646 bottom = &ptr[1];
647 }
648
649 *top = *val.top;
650 *bottom = *val.bottom;
651 }
652 else
653 {
654 if( exact )
655 {
656 delete [] top;
657 exact = false;
658 }
659 value = val.value;
660 }
661 return *this;
662 }
663
operator ==(const Real & val) const664 bool Real::operator==(const Real &val) const
665 {
666 if( exact != val.exact )
667 return ( floatx(*this) == floatx(val) );
668 else
669 {
670 if( exact )
671 return ( *top == *val.top && *bottom == *val.bottom );
672 else
673 return ( value == val.value );
674 }
675 }
676
operator <(const Real & val) const677 bool Real::operator<(const Real &val) const
678 {
679 if( exact != val.exact )
680 return ( floatx(*this) < floatx(val) );
681 else
682 {
683 if( exact )
684 {
685 if( *bottom == *val.bottom )
686 return *top < *val.top;
687 else
688 {
689 Integer tmp, tmp2;
690 top->mul( *val.bottom, &tmp ); // we must have the same divider of the ratios
691 val.top->mul( *bottom, &tmp2 );
692
693 return tmp < tmp2;
694 }
695 }
696 else
697 return ( value < val.value );
698 }
699 }
700
operator >(const Real & val) const701 bool Real::operator>(const Real &val) const
702 {
703 if( exact != val.exact )
704 return ( floatx(*this) > floatx(val) );
705 else
706 {
707 if( exact )
708 {
709 if( *bottom == *val.bottom )
710 return *top > *val.top;
711 else
712 {
713 Integer tmp, tmp2;
714 top->mul( *val.bottom, &tmp ); // we must have the same divider of the ratios
715 val.top->mul( *bottom, &tmp2 );
716
717 return tmp > tmp2;
718 }
719 }
720 else
721 return ( value > val.value );
722 }
723 }
724
725
ipart(const Real & val)726 Real ipart(const Real &val) throw(error_obj)
727 {
728 if( val.exact )
729 return Real( *val.top / *val.bottom );
730 else
731 return Real( truncx(val.value) );
732 }
733
fpart(const Real & val)734 Real fpart(const Real &val) throw(error_obj)
735 {
736 if( val.exact )
737 return Real( *val.top % *val.bottom, *val.bottom );
738 else
739 return Real( val.value - truncx(val.value) );
740 }
741
sign(const Real & val)742 Real sign(const Real &val) throw(error_obj)
743 {
744 if(val.exact)
745 {
746 int32 res;
747
748 if( *val.top == 0 )
749 res = 0;
750 else
751 res = (*val.top < 0) ? -1 : 1;
752 return Real( res );
753 }
754 else
755 {
756 floatx res;
757
758 if( val.value == 0.0 )
759 res = 0.0;
760 else
761 res = (val.value < 0.0) ? -1.0 : 1.0;
762 return Real( res );
763 }
764 }
765
766 // converts a Real to a fraction
frac(const Real & val,int32 highest_bottom_value=1000)767 Real frac(const Real &val, int32 highest_bottom_value = 1000) throw(error_obj)
768 {
769 if( highest_bottom_value < 1 )
770 THROW_ERROR( ErrorType_Domain, _("Domain error: Value out of range.") );
771
772 if(val.exact)
773 {
774 if( *val.bottom <= highest_bottom_value )
775 return val;
776 else
777 return frac( Real( floatx(val) ), highest_bottom_value );
778 }
779 else
780 {
781 Real res;
782 floatx best_fault, test_fault, tmp, value;
783 bool negative;
784 int32 top, bottom;
785
786 if( val.value < 0.0 )
787 {
788 negative = true;
789 value = -val.value;
790 }
791 else
792 {
793 negative = false;
794 value = val.value;
795 }
796 value = value - truncx(value);
797 res = 0;
798
799 best_fault = value;
800
801 if(best_fault != 0.0)
802 {
803 floatx fbottom;
804
805 if( (bottom = highest_bottom_value / 2) == 0 )
806 bottom = 1;
807
808 top = int32( value * floatx(bottom) );
809
810 fbottom = floatx(bottom);
811
812 test_fault = value - (floatx(top)/fbottom);
813 do
814 {
815 if(test_fault > 0.0 )
816 top++;
817 else
818 {
819 if( bottom == highest_bottom_value )
820 break;
821 bottom++;
822 fbottom = floatx(bottom);
823 top--;
824 }
825
826 test_fault = value - (floatx(top)/fbottom);
827 tmp = (test_fault < 0.0) ? -test_fault : test_fault;
828
829 if(tmp < best_fault)
830 {
831 best_fault = tmp;
832 *res.top = top;
833 *res.bottom = bottom;
834 }
835 }while(test_fault != 0.0);
836 }
837
838 {
839 Integer tmp1, tmp2;
840 if( negative )
841 tmp1 = -val.value;
842 else
843 tmp1 = val.value;
844 tmp1.mul( *res.bottom, &tmp2 );
845 res.top->add( tmp2 );
846 res.legalice();
847 }
848
849 if( negative )
850 res.top->negate();
851
852 return res;
853 }
854 }
855
power(const Real & val,const Integer & nr)856 Real power(const Real &val, const Integer &nr) throw(error_obj)
857 {
858 Real res(0);
859
860 if( !val.isExact() ) // this error should never accour...that's why it's not translatable
861 THROW_ERROR( ErrorType_Domain, _("Domain error: Real power(val, nr), val must be a ratio.") );
862
863 if( val.isInteger() )
864 {
865 if( nr.isNegative() )
866 { // n^p p<0
867 Integer tmp, tnr;
868
869 tnr = nr;
870 tnr.negate();
871
872 val.top->power( tnr, &tmp );
873
874 *res.top = 1;
875 *res.bottom = tmp;
876 }
877 else
878 { // n^p p>=0
879 Integer tmp;
880
881 val.top->power( nr, &tmp );
882 *res.top = tmp;
883 *res.bottom = 1;
884 }
885 }
886 else
887 {
888 if( nr.isNegative() )
889 { // (n/m)^p p<0
890 Integer tmp, tnr;
891
892 tnr = nr;
893 tnr.negate();
894
895 val.top->power( tnr, &tmp );
896 *res.bottom = tmp;
897
898 val.bottom->power( tnr, &tmp );
899 *res.top = tmp;
900 }
901 else
902 { // (n/m)^p p>=0
903 Integer tmp;
904
905 val.top->power( nr, &tmp );
906 *res.top = tmp;
907
908 val.bottom->power( nr, &tmp );
909 *res.bottom = tmp;
910 }
911 }
912 return res;
913 }
914