- メインページ
- モジュール
- ネームスペース
- データ構造
- ファイル

- ファイル一覧
- グローバル

decimal.h
説明を見る。
1 //
2 // Copyright (c) 2003-2011, MIST Project, Nagoya University
3 // All rights reserved.
4 //
5 // Redistribution and use in source and binary forms, with or without modification,
6 // are permitted provided that the following conditions are met:
7 //
8 // 1. Redistributions of source code must retain the above copyright notice,
9 // this list of conditions and the following disclaimer.
10 //
11 // 2. Redistributions in binary form must reproduce the above copyright notice,
12 // this list of conditions and the following disclaimer in the documentation
13 // and/or other materials provided with the distribution.
14 //
15 // 3. Neither the name of the Nagoya University nor the names of its contributors
16 // may be used to endorse or promote products derived from this software
17 // without specific prior written permission.
18 //
19 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
20 // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
21 // FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
22 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
23 // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
24 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
25 // IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF
26 // THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27 //
28
33 #ifndef __INCLUDE_MIST_DECIMAL_H__
34 #define __INCLUDE_MIST_DECIMAL_H__
35
36
37 #ifndef __INCLUDE_MIST_CONF_H__
38 #include "mist_conf.h"
39 #endif
40
41 #ifndef __INCLUDE_MIST_LIMITS__
42 #include "../limits.h"
43 #endif
44
45 #include <iostream>
46 #include <string>
47 #include <cmath>
48 #include <cstdio>
49 #include <cstring>
50 #include <stdlib.h>
51
52
53 // mist名前空間の始まり
54 _MIST_BEGIN
55
56 // 以下のコードは Tomy 氏が http://www5.airnet.ne.jp/tomy/cpro/csource.htm にて
57 // 掲載しているソースを参考に作成し，演算ミス等のバグを修正し，独自の拡張を施したものである．
58
59
60 // 表現可能なデータの有効桁数は，MAX_KETA × log10（32768) となる
61 // 有効桁数を増やすほどメモリを多く使用するので注意が必要！！
62
63
64
65 namespace __decimal_util__
66 {
67 template < int N >
68 struct log10
69 {
70 _MIST_CONST( int, value, log10< N / 10 >::value + 1 );
71 };
72
73 template < >
74 struct log10< 1 >
75 {
76 _MIST_CONST( int, value, 0 );
77 };
78
79 template < >
80 struct log10< 0 >
81 {
82 _MIST_CONST( int, value, 0 );
83 };
84
85 template < int N >
86 struct pow2
87 {
88 _MIST_CONST( int, value, 2 * pow2< N - 1 >::value );
89 };
90
91 template < >
92 struct pow2< 0 >
93 {
94 _MIST_CONST( int, value, 1 );
95 };
96 }
97
98
99 template < unsigned int MAX_KETA >
100 class decimal
101 {
102 public:
103 typedef unsigned int size_type;
104 typedef int difference_type;
105 typedef unsigned short value_type;
106
107 private:
108 _MIST_CONST( long, RADIXBITS, 15 );
109 _MIST_CONST( long, RADIX, __decimal_util__::pow2< RADIXBITS >::value );
110 _MIST_CONST( long, RADIX1, RADIX - 1 );
111 _MIST_CONST( long, RADIX_2, RADIX / 2 );
112 _MIST_CONST( long, MAXEXP, __decimal_util__::pow2< RADIXBITS - 1 >::value - 1 );
113 _MIST_CONST( long, MINEXP, - __decimal_util__::pow2< RADIXBITS - 1 >::value );
114 //_MIST_CONST( long, MAXEXP, 16383 );
115 //_MIST_CONST( long, MINEXP, -16384 );
116
117 _MIST_CONST( long, NMPA, MAX_KETA );
118 _MIST_CONST( long, NMPA1, NMPA + 1 );
119 _MIST_CONST( long, NMPA2, NMPA1 * 2 );
120 _MIST_CONST( long, MMPA, NMPA1 * RADIXBITS * 3 / 40 + 1 );
121
122 protected:
123 bool sign_;
124 difference_type exp_;
125 bool zero_;
126 value_type data_[ NMPA1 ];
127
128 public:
129 decimal( ) : sign_( true ), exp_( 0 ), zero_( true )
130 {
131 memset( data_, 0, sizeof( value_type ) * NMPA1 );
132 }
133
134 #if 1
135 decimal( double a ) : sign_( true ), exp_( 0 ), zero_( true )
136 {
137 memset( data_, 0, sizeof( value_type ) * NMPA1 );
138
139 if( a != 0.0 )
140 {
141 zero_ = false;
142 if( a < 0.0 )
143 {
144 sign_ = false;
145 a = - a;
146 }
147 while( a >= static_cast< double >( RADIX ) )
148 {
149 exp_++;
150 a /= static_cast< double >( RADIX );
151 }
152 while( a < 1.0 )
153 {
154 exp_--;
155 a *= static_cast< double >( RADIX );
156 }
157
158 difference_type i = 0;
159 do
160 {
161 data_[ i ] = static_cast< unsigned int >( a );
162 a -= data_[ i ];
163 a *= RADIX;
164 } while( a != 0.0 && ++i <= NMPA );
165 }
166 }
167 #else
168 decimal( const double &v, difference_type keta = 16 ) : sign_( true ), exp_( 0 ), zero_( true )
169 {
170 memset( data_, 0, sizeof( value_type ) * NMPA1 );
171
172 if( v != 0.0 )
173 {
174 double a = v;
175 zero_ = false;
176
177 if( a < 0.0 )
178 {
179 sign_ = false;
180 a = - a;
181 }
182
183 decimal b( static_cast< difference_type >( a ) );
184 a -= static_cast< difference_type >( a );
185 decimal x( 1 );
186
187 for( difference_type i = 0 ; i < keta && a != 0.0 ; i++ )
188 {
189 a *= 10;
190 b *= 10;
191 b += decimal( static_cast< difference_type >( a ) );
192 a -= static_cast< difference_type >( a );
193 x *= 10;
194 }
195
196 b /= x;
197 b.sign_ = sign_;
198
199 operator =( b );
200 }
201 }
202 #endif
203
204 decimal( const difference_type &m ) : sign_( true ), exp_( 0 ), zero_( true )
205 {
206 memset( data_, 0, sizeof( value_type ) * NMPA1 );
207
208 if( m != 0 )
209 {
210 difference_type n = m;
211 zero_ = false;
212 if( n < 0 )
213 {
214 sign_ = false;
215 n = - n;
216 }
217 value_type *p, *q, w;
218 p = q = data_;
219 while( n != 0 )
220 {
221 *p++ = static_cast< value_type >( n % RADIX );
222 n /= RADIX;
223 exp_++;
224 }
225 exp_--;
226 p--;
227 while( q < p )
228 {
229 w = *p;
230 *p-- = *q;
231 *q++ = w;
232 }
233 }
234 }
235
236 decimal( const decimal &a ) : sign_( a.sign_ ), exp_( a.exp_ ), zero_( a.zero_ )
237 {
238 memcpy( data_, a.data_, sizeof( value_type ) * NMPA1 );
239 }
240
241 decimal( const ::std::string &str ) : sign_( true ), exp_( 0 ), zero_( true )
242 {
243 operator =( read( str ) );
244 }
245
246 const decimal &operator =( const decimal &a )
247 {
248 if( &a != this )
249 {
250 sign_ = a.sign_;
251 exp_ = a.exp_;
252 zero_ = a.zero_;
253 memcpy( data_, a.data_, sizeof( value_type ) * NMPA1 );
254 }
255
256 return ( *this );
257 }
258
259 const decimal operator -( ) const
260 {
261 decimal x( *this );
262 x.sign_ = !x.sign_;
263 return( x );
264 }
265
266 const decimal &operator +=( const decimal &a )
267 {
268 if( zero_ )
269 {
270 return( operator =( a ) );
271 }
272 else if( a.zero_ )
273 {
274 return( *this );
275 }
276 else if( sign_ == a.sign_ )
277 {
278 aadd( *this, a );
279 return( *this );
280 }
281
282 int cmp = acmp( *this, a );
283 if( cmp == 0 )
284 {
285 return( operator =( decimal( ) ) );
286 }
287 else if( cmp > 0 )
288 {
289 asub( *this, a );
290 return( *this );
291 }
292 decimal tmp( *this );
293 *this = a;
294 asub( *this, tmp );
295 return( *this );
296 }
297
298 const decimal &operator -=( const decimal &a )
299 {
300 return( operator +=( -a ) );
301 }
302
303 const decimal &operator *=( const decimal &a )
304 {
305 if( zero_ )
306 {
307 return( *this );
308 }
309 else if( a.zero_ )
310 {
311 return( operator =( decimal( ) ) );
312 }
313
314 difference_type i, j, u, r, xp;
315 value_type x[ NMPA2 ];
316 memset( x, 0, sizeof( value_type ) * NMPA2 );
317
318 for( i = NMPA, xp = NMPA2 - 1 ; i >= 0 ; i--, xp-- )
319 {
320 if( a.data_[ i ] )
321 {
322 u = 0;
323 for( j = NMPA, r = xp ; j >= 0; j--, r-- )
324 {
325 u += static_cast< size_type >( data_[ j ] ) * static_cast< size_type >( a.data_[ i ] ) + x[ r ];
326 x[ r ] = static_cast< value_type >( u ) & RADIX1;
327 u >>= RADIXBITS;
328 }
329 while( u )
330 {
331 u += x[ r ];
332 x[ r-- ] = static_cast< value_type >( u ) & RADIX1;
333 u >>= RADIXBITS;
334 }
335 }
336 }
337 sign_ = sign_ == a.sign_;
338 exp_ += a.exp_ + 1;
339 for( i = 0, xp = 0 ; i < NMPA2 ; i++, xp++ )
340 {
341 if( x[ xp ] )
342 {
343 break;
344 }
345 exp_--;
346 }
347 if( exp_ > MAXEXP )
348 {
349 ::std::cerr << "Error : Overflow!!" << ::std::endl;
350 return( operator =( maximum( ) ) );
351 }
352 if( exp_ < MINEXP )
353 {
354 ::std::cerr << "Error : Underflow!!" << ::std::endl;
355 return( operator =( decimal( ) ) );
356 }
357 for( j = 0 ; ( j <= NMPA ) && ( i < NMPA2 ) ; i++, j++, xp++ )
358 {
359 data_[ j ] = x[ xp ];
360 }
361 if( j <= NMPA )
362 {
363 for( ; j <= NMPA ; j++ )
364 {
365 data_[ j ] = 0;
366 }
367 }
368 else if( x[ xp ] >= RADIX_2 )
369 {
370 for( i = NMPA ; ++( data_[ i ] ) & RADIX ; i-- )
371 {
372 data_[ i ] &= RADIX1;
373 }
374 }
375 return( *this );
376 }
377
378 const decimal &operator *=( difference_type x )
379 {
380 if( zero_ )
381 {
382 return( *this );
383 }
384 else if( x == 0 )
385 {
386 return( operator =( decimal( ) ) );
387 }
388 else if( x < 0 )
389 {
390 sign_ = !sign_;
391 x = - x;
392 }
393 difference_type xl = static_cast< difference_type >( x );
394 difference_type t = 0, i;
395 for( i = NMPA ; i >= 0 ; i-- )
396 {
397 t += static_cast< difference_type >( data_[ i ] * xl );
398 data_[ i ] = static_cast< value_type >( t ) & RADIX1;
399 t >>= RADIXBITS;
400 }
401 while( t )
402 {
403 if( exp_ == MAXEXP )
404 {
405 ::std::cerr << "Error : Overflow!!" << ::std::endl;
406 return( operator =( maximum( ) ) );
407 }
408 for( i = NMPA ; i > 0 ; i-- )
409 {
410 data_[ i ] = data_[ i - 1 ];
411 }
412 data_[ i ] = static_cast< value_type >( t ) & RADIX1;
413 t >>= RADIXBITS;
414 exp_++;
415 }
416 return( *this );
417 }
418
419
420 const decimal &operator /=( decimal b )
421 {
422 if( zero_ )
423 {
424 return( *this );
425 }
426 else if( b.zero_ )
427 {
428 // ゼロ除算
429 std::cerr << "zero_ division!!" << std::endl;
430 return( *this );
431 }
432
433 decimal a = *this;
434 decimal c, t;
435 difference_type a0, b0, bcom, cmp, d, i;
436 value_type *p;
437
438 c.zero_ = false;
439 c.sign_ = a.sign_ == b.sign_;
440 a.sign_ = b.sign_ = true;
441 if( ( b.data_[ 0 ] << 1 ) & RADIX )
442 {
443 bcom = 1;
444 }
445 else
446 {
447 bcom = RADIX / ( b.data_[ 0 ] + 1 );
448 b *= bcom;
449 }
450 c.exp_ = a.exp_ - b.exp_;
451 a.exp_ = b.exp_ = 0;
452 if( acmp( a, b ) < 0 )
453 {
454 c.exp_--;
455 a.exp_++;
456 }
457 i = 0;
458 p = c.data_;
459 b0 = b.data_[ 0 ] + 1;
460 while( true )
461 {
462 a0 = a.data_[ 0 ];
463 if( a0 > b0 )
464 {
465 d = a0 / b0;
466 }
467 else
468 {
469 a0 = a0 * RADIX + a.data_[ 1 ];
470 d = a0 / b0;
471 if( d < 1 )
472 {
473 d = 1;
474 }
475 }
476 while( true )
477 {
478 t = a - b * d;
479 cmp = acmp( t, b );
480 if( cmp < 0 )
481 {
482 break;
483 }
484 d++;
485 if( cmp == 0 )
486 {
487 break;
488 }
489 }
490 if( i > NMPA )
491 {
492 if( d >= RADIX_2 )
493 {
494 for( i = NMPA, p = c.data_ + i; ++*p & RADIX; i-- )
495 {
496 *p-- &= RADIX1;
497 }
498 }
499 break;
500 }
501 *p++ = static_cast< value_type >( d );
502 i++;
503 if( cmp == 0 )
504 {
505 for( ; i <= NMPA ; i++ )
506 {
507 *p++ = 0;
508 }
509 break;
510 }
511 a = t;
512 a.exp_++;
513 if( acmp( a, b ) < 0 )
514 {
515 *p++ = 0;
516 i++;
517 if( i > NMPA )
518 {
519 break;
520 }
521 a.exp_++;
522 }
523 }
524 if( bcom != 1 )
525 {
526 c *= bcom;
527 }
528 return( operator =( c ) );
529 }
530
531 const decimal &operator /=( difference_type x )
532 {
533 if( x == 0 )
534 {
535 // ゼロ除算
536 std::cerr << "zero division!!" << std::endl;
537 return( operator =( maximum( ) ) );
538 }
539 else if( zero_ )
540 {
541 return( *this );
542 }
543
544 difference_type i, u, n;
545
546 if( x < 0 )
547 {
548 sign_ = !sign_;
549 x = - x;
550 }
551 if( data_[ 0 ] < static_cast< value_type >( x ) )
552 {
553 u = data_[ 0 ];
554 n = 1;
555 exp_--;
556 }
557 else
558 {
559 u = 0;
560 n = 0;
561 }
562 for( i = n ; i <= NMPA ; i++ )
563 {
564 u = ( u << RADIXBITS ) + data_[ i ];
565 data_[ i - n ] = static_cast< value_type >( u / x );
566 u %= x;
567 }
568 if( n )
569 {
570 data_[ i - n ] = static_cast< value_type >( ( u << RADIXBITS ) / x );
571 u %= x;
572 }
573 if( 2 * u >= x )
574 {
575 for( i = NMPA ; ++( data_[ i ] ) & RADIX ; i-- )
576 {
577 data_[ i ] &= RADIX1;
578 }
579 }
580 return( *this );
581 }
582
583 decimal &operator ++( ) // 前置型
584 {
585 *this += 1;
586 return( *this );
587 }
588
589 const decimal operator ++( int ) // 後置型
590 {
591 decimal old_val( *this );
592 *this += 1;
593 return( old_val );
594 }
595
596 decimal &operator --( ) // 前置型
597 {
598 *this -= 1;
599 return( *this );
600 }
601
602 const decimal operator --( int ) // 後置型
603 {
604 decimal old_val( *this );
605 *this -= 1;
606 return( old_val );
607 }
608
609 bool operator ==( const decimal &a ) const { return( cmp( *this, a ) == 0 ); }
610 bool operator !=( const decimal &a ) const { return( !( *this == a ) ); }
611 bool operator < ( const decimal &a ) const { return( !( *this >= a ) ); }
612 bool operator <=( const decimal &a ) const { return( a >= *this ); }
613 bool operator > ( const decimal &a ) const { return( a < *this ); }
614 bool operator >=( const decimal &a ) const { return( cmp( *this, a ) >= 0 ); }
615
616 const ::std::string to_string( ) const
617 {
618 difference_type nmpa = NMPA;
619 difference_type mmpa = MMPA - 1;
620
621 if( mmpa < static_cast< difference_type >( std::abs( static_cast< double >( exp_ ) ) ) )
622 {
623 decimal tmp( *this );
624 tmp.exp_ = 0;
625 char str[ 10 ];
626 sprintf( str, "E%+d", exp_ );
627 return( tmp.to_string( ) + str );
628 }
629 else if( *this == 10000 )
630 {
631 return( "+10000" );
632 }
633 else if( *this == -10000 )
634 {
635 return( "-10000" );
636 }
637
638 if( zero_ )
639 {
640 return( "0" );
641 }
642
643 unsigned int c[ MMPA ], t[ NMPA1 ], w;
644 difference_type i, cp;
645
646 for( i = 0 ; i <= mmpa ; i++ )
647 {
648 c[ i ] = 0;
649 }
650
651 ::std::string s = "";
652
653 if( sign_ )
654 {
655 s = "+";
656 }
657 else
658 {
659 s = "-";
660 }
661
662 decimal a( *this );
663
664 cp = 0;
665 difference_type n = a.exp_, u, m = 0, mp = -1;
666 if( n >= 0 )
667 {
668 for( i = 0 ; i <= n ; i++ )
669 {
670 t[ i ] = data_[ i ];
671 }
672 while( true )
673 {
674 u = 0;
675 for( i = 0 ; i <= n ; i++ )
676 {
677 u = ( u << RADIXBITS ) + t[ i ];
678 t[ i ] = static_cast< value_type >( u / 10000 );
679 u %= 10000;
680 }
681 c[ cp++ ] = static_cast< value_type >( u );
682 m++;
683 mp++;
684 bool zflag = true;
685 for( i = 0 ; i <= n ; i++ )
686 {
687 if( t[ i ] )
688 {
689 zflag = false;
690 break;
691 }
692 }
693 if( zflag )
694 {
695 break;
696 }
697 }
698 for( i = 0 ; i < m - i - 1 ; i++ )
699 {
700 w = c[ i ];
701 c[ i ] = c[ m - i - 1 ];
702 c[ m - i - 1 ] = w;
703 }
704 for( i = 0 ; i <= nmpa - ( n + 1 ) ; i++ )
705 {
706 a.data_[ i ] = a.data_[ i + n + 1 ];
707 }
708 for( ; i <= nmpa ; i++)
709 {
710 a.data_[ i ] = 0;
711 }
712 a.exp_ = -1;
713 }
714 else if( n < -1 )
715 {
716 while( true )
717 {
718 a *= 10000;
719 if( a.exp_ < 0 )
720 {
721 mp--;
722 }
723 else
724 {
725 c[ cp++ ] = a.data_[ 0 ];
726 m++;
727 break;
728 }
729 }
730 for( i = 0 ; i < nmpa ; i++ )
731 {
732 a.data_[ i ] = a.data_[ i + 1 ];
733 }
734 a.data_[ i ] = 0;
735 a.exp_ = -1;
736 }
737
738 a *= 10000;
739 m++;
740
741 if( a.exp_ == 0 )
742 {
743 c[ cp++ ] = a.data_[ 0 ];
744 a.data_[ 0 ] = 0;
745 }
746 else
747 {
748 c[ cp++ ] = 0;
749 }
750
751 while( m <= mmpa )
752 {
753 if( a.is_zero( ) )
754 {
755 a = decimal( );
756 break;
757 }
758 a *= 10000;
759 c[ cp++ ] = a.data_[ 0 ];
760 m++;
761 a.data_[ 0 ] = 0;
762 }
763
764 if( c[ mmpa ] < 5000 )
765 {
766 c[ mmpa ] = 0;
767 }
768 else
769 {
770 c[ mmpa ] = 0;
771 for( m = mmpa - 1 ; m >= 0 ; m-- )
772 {
773 c[ m ]++;
774 if( c[ m ] != 10000 )
775 {
776 break;
777 }
778 c[ m ] = 0;
779 }
780 }
781
782 for( i = mmpa ; i >= 0 ; i-- )
783 {
784 if( c[ i ] )
785 {
786 break;
787 }
788 mmpa--;
789 }
790
791 if( mmpa < 0 )
792 {
793 return( "0" );
794 }
795
796 if( mp < 0 )
797 {
798 s += "0.";
799 while( ++mp < 0 )
800 {
801 s += "0000";
802 }
803 mp = -1;
804 }
805
806 cp = 0;
807
808 char str[ 10 ];
809 if( mp < 0 )
810 {
811 sprintf( str, "%04u", c[ cp++ ] );
812 s += str;
813 }
814 else
815 {
816 sprintf( str, "%u", c[ cp++ ] );
817 s += str;
818 if( !mp )
819 {
820 s += ".";
821 }
822 }
823
824 for( m = 1 ; m <= mmpa ; m++ )
825 {
826 sprintf( str, "%04u", c[ cp++ ] );
827 s += str;
828 if( m == mp )
829 {
830 s += ".";
831 }
832 }
833
834 if( m < mp )
835 {
836 do
837 {
838 s += "0000";
839 } while( m++ < mp );
840 }
841
842 return( s );
843 }
844
845 int to_int( ) const
846 {
847 if( zero_ || exp_ < 0 )
848 {
849 return( 0 );
850 }
851 else if( exp_ > NMPA )
852 {
853 ::std::cerr << "Error : Overflow!!" << ::std::endl;
854 return( sign_ ? type_limits< int >::maximum( ) : type_limits< int >::minimum( ) );
855 }
856 else
857 {
858 size_type i, index;
859 for( index = 0 ; index <= NMPA ; index++ )
860 {
861 if( data_[ index ] )
862 {
863 break;
864 }
865 }
866 if( exp_ >= 3 || ( exp_ == 2 && data_[ 0 ] > 2 ) )
867 {
868 ::std::cerr << "Error : Overflow!!" << ::std::endl;
869 return( sign_ ? type_limits< int >::maximum( ) : type_limits< int >::minimum( ) );
870 }
871 else
872 {
873 int d = 0;
874 for( i = index ; i <= index + exp_ && i <= NMPA ; i++ )
875 {
876 d = d * RADIX + data_[ i ];
877 }
878 return( sign_ ? d : -d );
879 }
880 }
881 }
882
883 double to_double( ) const
884 {
885 if( zero_ )
886 {
887 return( 0.0 );
888 }
889 if( exp_ > 68 )
890 {
891 ::std::cerr << "Error : Overflow!!" << ::std::endl;
892 return( type_limits< double >::maximum( ) );
893 }
894 if( exp_ < -68 )
895 {
896 ::std::cerr << "Error : Underflow!!" << ::std::endl;
897 return( type_limits< double >::minimum( ) );
898 }
899
900 decimal a( *this );
901 size_type i;
902 while( a.data_[ 0 ] == 0 )
903 {
904 for( i = 0 ; i < NMPA ; i++ )
905 {
906 a.data_[ i ] = a.data_[ i + 1 ];
907 }
908 a.data_[ i ];
909 a.exp_--;
910 }
911
912 double d = static_cast< double >( a.data_[ 0 ] );
913 for( i = 1; i < 7; i++ )
914 {
915 d = d * static_cast< double >( RADIX ) + static_cast< double >( a.data_[ i ] );
916 a.exp_--;
917 }
918 if( a.data_[ i ] >= RADIX_2 )
919 {
920 d += 1.0;
921 }
922 while( a.exp_ > 0 )
923 {
924 d *= static_cast< double >( RADIX );
925 a.exp_--;
926 }
927 while( a.exp_ < 0 )
928 {
929 d /= static_cast< double >( RADIX );
930 a.exp_++;
931 }
932 return( a.sign_ ? d : -d );
933 }
934
935 bool is_zero( ) const
936 {
937 if( zero_ )
938 {
939 return( true );
940 }
941 for( difference_type i = 0 ; i <= NMPA ; i++ )
942 {
943 if( data_[ i ] )
944 {
945 return( false );
946 }
947 }
948 return( true );
949 }
950
951 bool is_plus( ) const
952 {
953 return( sign_ );
954 }
955
956
957 public: // 定数
958 static decimal maximum( )
959 {
960 static decimal max;
961 if( max.is_zero( ) )
962 {
963 max.sign_ = true;
964 max.exp_ = MAXEXP;
965 for( difference_type i = 0 ; i <= NMPA1 ; i++ )
966 {
967 max.data_[ i ] = RADIX1;
968 }
969 }
970 return( max );
971 }
972
973 static decimal zero( )
974 {
975 return decimal( );
976 }
977
978 static decimal pai( )
979 {
980 decimal a0( 1 );
981 decimal b0( decimal::sqrt( decimal( 2 ) ) / 2 );
982 decimal t0( "0.25" );
983 decimal a1, b1, t1;
984
985 difference_type _2[] = { 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096 };
986 for( int n = 0 ; n < 10 ; n++ )
987 {
988 a1 = ( a0 + b0 ) / 2;
989 b1 = decimal::sqrt( a0 * b0 );
990 t1 = t0 - _2[ n ] * ( a0 - a1 ) * ( a0 - a1 );
991
992 a0 = a1;
993 b0 = b1;
994 t0 = t1;
995 }
996
997 return( ( a0 + b0 ) * ( a0 + b0 ) / ( t0 * 4 ) );
998 }
999
1000 static decimal sqrt( const decimal &a )
1001 {
1002 if( a.is_zero( ) )
1003 {
1004 return mist::decimal< MAX_KETA >( );
1005 }
1006 else if( !a.is_plus( ) )
1007 {
1008 ::std::cerr << "Error : Illegal parameter!!" << ::std::endl;
1009 return( a );
1010 }
1011
1012 mist::decimal< MAX_KETA > s( 1 );
1013 mist::decimal< MAX_KETA > b = a;
1014 while( s < b )
1015 {
1016 s *= 2;
1017 b /= 2;
1018 }
1019
1020 do
1021 {
1022 b = s;
1023 s = ( ( a / s ) + s ) / 2;
1024 } while( s < b );
1025 return( b );
1026 }
1027
1028 protected:
1029 static int cmp( const decimal &a, const decimal &b )
1030 {
1031 if( a.zero_ )
1032 {
1033 if( b.zero_ )
1034 {
1035 return( 0 );
1036 }
1037 return( b.sign_ ? -1 : 1 );
1038 }
1039 if( b.zero_ )
1040 {
1041 return( a.sign_ ? 1 : -1 );
1042 }
1043 if( a.sign_ == b.sign_ )
1044 {
1045 int mca = acmp( a, b );
1046 return( a.sign_ ? mca : -mca );
1047 }
1048 return( a.sign_ ? 1 : -1 );
1049 }
1050
1051 static int acmp( const decimal &a, const decimal &b )
1052 {
1053 if( a.zero_ )
1054 {
1055 return( b.zero_ ? 0 : -1 );
1056 }
1057 else if( b.zero_ )
1058 {
1059 return( 1 );
1060 }
1061
1062 difference_type ia, ib;
1063 difference_type aexp = a.exp_;
1064 difference_type bexp = b.exp_;
1065
1066 for( ia = 0 ; ia <= NMPA && aexp != MINEXP ; ia++ )
1067 {
1068 if( a.data_[ ia ] != 0 )
1069 {
1070 break;
1071 }
1072 aexp--;
1073 }
1074 for( ib = 0 ; ib <= NMPA && bexp != MINEXP ; ib++ )
1075 {
1076 if( b.data_[ ib ] != 0 )
1077 {
1078 break;
1079 }
1080 bexp--;
1081 }
1082
1083 if( aexp != bexp )
1084 {
1085 return( aexp > bexp ? 1: -1 );
1086 }
1087 for( difference_type i = ia ; i <= NMPA ; i++ )
1088 {
1089 if( a.data_[ i ] != b.data_[ i ] )
1090 {
1091 return( a.data_[ i ] > b.data_[ i ] ? 1 : -1 );
1092 }
1093 }
1094 return( 0 );
1095 }
1096
1097 static difference_type aprs( const decimal &a, decimal &b )
1098 {
1099 if( acmp( a, b ) < 0 )
1100 {
1101 ::std::cerr << "Error : |a| < |b| " << ::std::endl;
1102 return( -999 );
1103 }
1104 difference_type n = a.exp_ - b.exp_;
1105 if( n == 0 )
1106 {
1107 return( -1 );
1108 }
1109 else if( n > NMPA )
1110 {
1111 b = decimal( );
1112 return( 0 );
1113 }
1114 difference_type r = b.data_[ NMPA - n + 1 ], i;
1115 b.exp_ = a.exp_;
1116 for( i = NMPA ; i >= n ; i-- )
1117 {
1118 b.data_[ i ] = b.data_[ i - n ];
1119 }
1120 for( ; i >= 0 ; i-- )
1121 {
1122 b.data_[ i ] = 0;
1123 }
1124 return( r );
1125 }
1126
1127 static void aadd( decimal &a, decimal b )
1128 {
1129 if( b.zero_ )
1130 {
1131 return;
1132 }
1133
1134 bool sign = a.sign_;
1135 if( a.zero_ )
1136 {
1137 a = b;
1138 a.sign_ = sign;
1139 return;
1140 }
1141
1142 decimal t;
1143
1144 int cmp = acmp( a, b );
1145 if( cmp < 0 )
1146 {
1147 t = a;
1148 a = b;
1149 b = t;
1150 a.sign_ = sign;
1151 }
1152 difference_type i = aprs( a, b );
1153 if( i == -999 )
1154 {
1155 a = decimal( );
1156 return;
1157 }
1158 size_type u = 0;
1159 if( i >= RADIX_2 )
1160 {
1161 u = 1;
1162 }
1163 for( i = NMPA ; i >= 0 ; i-- )
1164 {
1165 u += ( a.data_[ i ] + b.data_[ i ] );
1166 a.data_[ i ] = static_cast< value_type >( u & RADIX1 );
1167 u >>= RADIXBITS;
1168 }
1169 if( u == 0 )
1170 {
1171 a.sign_ = sign;
1172 return;
1173 }
1174 if( a.exp_ == MAXEXP )
1175 {
1176 ::std::cerr << "Error : Overflow!!" << ::std::endl;
1177 a = maximum( );
1178 return;
1179 }
1180 difference_type k = a.data_[ NMPA ];
1181 for( i = NMPA ; i > 0 ; i-- )
1182 {
1183 a.data_[ i ] = a.data_[ i - 1 ];
1184 }
1185 a.data_[ i ] = static_cast< value_type >( u );
1186 a.exp_++;
1187 if( k >= RADIX_2 )
1188 {
1189 for( i = NMPA ; ++( a.data_[ i ] ) & RADIX ; i-- )
1190 {
1191 a.data_[ i ] &= RADIX1;
1192 }
1193 }
1194 a.sign_ = sign;
1195 }
1196
1197 static void asub( decimal &a, decimal b )
1198 {
1199 difference_type i = aprs( a, b );
1200 if( i == -999 )
1201 {
1202 return;
1203 }
1204
1205 bool sign = a.sign_;
1206 difference_type u = 0, n, m;
1207
1208 if( i >= RADIX_2 )
1209 {
1210 u = 1;
1211 }
1212 for( i = NMPA ; i >= 0 ; i-- )
1213 {
1214 u = a.data_[ i ] - b.data_[ i ] - u;
1215 a.data_[ i ] = static_cast< value_type >( u & RADIX1 );
1216 u = ( u >> RADIXBITS ) & 1;
1217 }
1218 for( n = 0 ; n <= NMPA ; n++ )
1219 {
1220 if( a.data_[ n ] )
1221 {
1222 break;
1223 }
1224 }
1225 if( n == 0 )
1226 {
1227 a.sign_ = sign;
1228 return;
1229 }
1230 if( a.exp_ > MAXEXP - n )
1231 {
1232 m = MAXEXP - a.exp_;
1233 }
1234 else
1235 {
1236 m = n;
1237 }
1238 for( i = 0 ; i <= NMPA - m; i++ )
1239 {
1240 a.data_[ i ] = a.data_[ n + i ];
1241 }
1242 for( i = 1 ; i <= m ; i++ )
1243 {
1244 a.data_[ NMPA - m + i ] = 0;
1245 }
1246 a.exp_ -= m;
1247 a.sign_ = sign;
1248 }
1249
1250 static decimal read( const ::std::string &str )
1251 {
1252 ::std::string::const_iterator p = str.begin( );
1253 bool sign = true;
1254
1255 if( *p == '-' )
1256 {
1257 sign = 0;
1258 p++;
1259 }
1260 else if( *p == '+' )
1261 {
1262 p++;
1263 }
1264
1265 decimal a, t;
1266 int c, exp = 0;
1267
1268 t.zero_ = false;
1269
1270 int pflag = 0, zflag = 0;
1271 for( ; p != str.end( ) ; p++ )
1272 {
1273 if( *p == '.' )
1274 {
1275 if( pflag )
1276 {
1277 ::std::cerr << "Error : Illegal parameter!!" << ::std::endl;
1278 return decimal( );
1279 }
1280 pflag = true;
1281 }
1282 else if( *p == 'E' || *p == 'e' )
1283 {
1284 p++;
1285 // a *= atoi( ::std::string( p, str.end( ) ).c_str( ) );
1286 exp -= atoi( ::std::string( p, str.end( ) ).c_str( ) );
1287 break;
1288 }
1289 else if( *p == '+' || *p == '-' )
1290 {
1291 exp -= atoi( ::std::string( p, str.end( ) ).c_str( ) );
1292 break;
1293 }
1294 else if( *p != ' ' && ( *p < '0' || *p > '9' ) )
1295 {
1296 ::std::cerr << "Error : Illegal parameter!!" << ::std::endl;
1297 return decimal( );
1298 }
1299 else
1300 {
1301 c = *p - '0';
1302 if( c != 0 || pflag != 0 )
1303 {
1304 zflag = -1;
1305 }
1306 else
1307 {
1308 if( zflag == 0 )
1309 {
1310 zflag = 1;
1311 }
1312 else if( zflag == 1 )
1313 {
1314 ::std::cerr << "Error : Illegal parameter!!" << ::std::endl;
1315 return decimal( );
1316 }
1317 }
1318
1319 a *= 10;
1320
1321 if( c )
1322 {
1323 t.data_[ 0 ] = c;
1324 a += t;
1325 }
1326 if( pflag == 1 )
1327 {
1328 exp++;
1329 }
1330 }
1331 }
1332 if( exp > MAXEXP )
1333 {
1334 ::std::cerr << "Error : Overflow!!" << ::std::endl;
1335 return( maximum( ) );
1336 }
1337 if( exp < MINEXP )
1338 {
1339 ::std::cerr << "Error : Underflow!!" << ::std::endl;
1340 return decimal( );
1341 }
1342 a.sign_ = sign;
1343 if( exp > 0 )
1344 {
1345 while( --exp >= 0 )
1346 {
1347 a /= 10;
1348 }
1349 }
1350 else
1351 {
1352 while( ++exp <= 0 )
1353 {
1354 a *= 10;
1355 }
1356 }
1357 if( a.is_zero( ) )
1358 {
1359 return decimal( );
1360 }
1361 else
1362 {
1363 return( a );
1364 }
1365
1366 }
1367 };
1368
1369 template < unsigned int MAX_KETA > inline const decimal< MAX_KETA > operator +( const decimal< MAX_KETA > &v1, const decimal< MAX_KETA > &v2 ){ return( decimal< MAX_KETA >( v1 ) += v2 ); }
1370 template < unsigned int MAX_KETA > inline const decimal< MAX_KETA > operator -( const decimal< MAX_KETA > &v1, const decimal< MAX_KETA > &v2 ){ return( decimal< MAX_KETA >( v1 ) -= v2 ); }
1371 template < unsigned int MAX_KETA > inline const decimal< MAX_KETA > operator *( const decimal< MAX_KETA > &v1, const decimal< MAX_KETA > &v2 ){ return( decimal< MAX_KETA >( v1 ) *= v2 ); }
1372 template < unsigned int MAX_KETA > inline const decimal< MAX_KETA > operator /( const decimal< MAX_KETA > &v1, const decimal< MAX_KETA > &v2 ){ return( decimal< MAX_KETA >( v1 ) /= v2 ); }
1373
1374 template < unsigned int MAX_KETA > inline const decimal< MAX_KETA > operator *( const decimal< MAX_KETA > &v1, const typename decimal< MAX_KETA >::difference_type &v2 ){ return( decimal< MAX_KETA >( v1 ) *= v2 ); }
1375 template < unsigned int MAX_KETA > inline const decimal< MAX_KETA > operator *( const typename decimal< MAX_KETA >::difference_type &v1, const decimal< MAX_KETA > &v2 ){ return( decimal< MAX_KETA >( v2 ) *= v1 ); }
1376 template < unsigned int MAX_KETA > inline const decimal< MAX_KETA > operator /( const decimal< MAX_KETA > &v1, const typename decimal< MAX_KETA >::difference_type &v2 ){ return( decimal< MAX_KETA >( v1 ) /= v2 ); }
1377 template < unsigned int MAX_KETA > inline const decimal< MAX_KETA > operator /( const typename decimal< MAX_KETA >::difference_type &v1, const decimal< MAX_KETA > &v2 ){ return( decimal< MAX_KETA >( v1 ) /= v2 ); }
1378
1379 template < unsigned int MAX_KETA > inline const decimal< MAX_KETA > operator +( const decimal< MAX_KETA > &v1, const typename decimal< MAX_KETA >::difference_type &v2 ){ return( decimal< MAX_KETA >( v1 ) += v2 ); }
1380 template < unsigned int MAX_KETA > inline const decimal< MAX_KETA > operator +( const typename decimal< MAX_KETA >::difference_type &v1, const decimal< MAX_KETA > &v2 ){ return( decimal< MAX_KETA >( v2 ) += v1 ); }
1381 template < unsigned int MAX_KETA > inline const decimal< MAX_KETA > operator -( const decimal< MAX_KETA > &v1, const typename decimal< MAX_KETA >::difference_type &v2 ){ return( decimal< MAX_KETA >( v1 ) -= v2 ); }
1382 template < unsigned int MAX_KETA > inline const decimal< MAX_KETA > operator -( const typename decimal< MAX_KETA >::difference_type &v1, const decimal< MAX_KETA > &v2 ){ return( decimal< MAX_KETA >( v1 ) -= v2 ); }
1383
1384
1385 template < unsigned int MAX_KETA > inline const decimal< MAX_KETA > operator *( const decimal< MAX_KETA > &v1, const double &v2 ){ return( decimal< MAX_KETA >( v2 ) *= v1 ); }
1386 template < unsigned int MAX_KETA > inline const decimal< MAX_KETA > operator *( const double &v1, const decimal< MAX_KETA > &v2 ){ return( decimal< MAX_KETA >( v1 ) *= v2 ); }
1387 template < unsigned int MAX_KETA > inline const decimal< MAX_KETA > operator /( const decimal< MAX_KETA > &v1, const double &v2 ){ return( decimal< MAX_KETA >( v1 ) /= decimal< MAX_KETA >( v2 ) ); }
1388 template < unsigned int MAX_KETA > inline const decimal< MAX_KETA > operator /( const double &v1, const decimal< MAX_KETA > &v2 ){ return( decimal< MAX_KETA >( v1 ) /= v2 ); }
1389
1390 template < unsigned int MAX_KETA > inline const decimal< MAX_KETA > operator +( const decimal< MAX_KETA > &v1, const double &v2 ){ return( decimal< MAX_KETA >( v2 ) += v1 ); }
1391 template < unsigned int MAX_KETA > inline const decimal< MAX_KETA > operator +( const double &v1, const decimal< MAX_KETA > &v2 ){ return( decimal< MAX_KETA >( v1 ) += v2 ); }
1392 template < unsigned int MAX_KETA > inline const decimal< MAX_KETA > operator -( const decimal< MAX_KETA > &v1, const double &v2 ){ return( decimal< MAX_KETA >( v1 ) -= decimal< MAX_KETA >( v2 ) ); }
1393 template < unsigned int MAX_KETA > inline const decimal< MAX_KETA > operator -( const double &v1, const decimal< MAX_KETA > &v2 ){ return( decimal< MAX_KETA >( v1 ) -= v2 ); }
1394
1395
1396 template < unsigned int MAX_KETA >
1397 inline std::ostream &operator <<( std::ostream &out, const decimal< MAX_KETA > &a )
1398 {
1399 out << a.to_string( );
1400 return( out );
1401 }
1402
1403
1404 // mist名前空間の終わり
1405 _MIST_END
1406
1407 namespace std
1408 {
1409 template < unsigned int MAX_KETA >
1410 inline mist::decimal< MAX_KETA > sqrt( const mist::decimal< MAX_KETA > &a )
1411 {
1412 return( mist::decimal< MAX_KETA >::sqrt( a ) );
1413 }
1414 }
1415
1416
1417
1418 #endif // __INCLUDE_MIST_DECIMAL_H__
1419