|
| 1 | +#pragma once |
| 2 | + |
| 3 | +#ifdef __cplusplus |
| 4 | +extern "C" { |
| 5 | +#endif |
| 6 | + |
| 7 | +#define PD_SCALE9 (1000000000L) |
| 8 | +#define PD_EMA_MAX_DIFF 4145 // maximum slots before reset |
| 9 | +#define PD_EMA_EXPO -9 // exponent of temporary storage |
| 10 | +#define PD_EMA_DECAY (-117065) // 1e9*-log(2)/5921 |
| 11 | +#define PC_FACTOR_SIZE 18 |
| 12 | + |
| 13 | +#define pd_new( n,v,e ) {(n)->v_=v;(n)->e_=e;} |
| 14 | +#define pd_set( r,n ) (r)[0] = (n)[0] |
| 15 | +#define pd_new_scale(n,v,e) {(n)->v_=v;(n)->e_=e;pd_scale(n);} |
| 16 | + |
| 17 | +// decimal number format |
| 18 | +typedef struct pd |
| 19 | +{ |
| 20 | + int32_t e_; |
| 21 | + int64_t v_; |
| 22 | +} pd_t; |
| 23 | + |
| 24 | +static void pd_scale( pd_t *n ) |
| 25 | +{ |
| 26 | + int const neg = n->v_ < 0L; |
| 27 | + int64_t v = neg ? -n->v_ : n->v_; // make v positive for loop condition |
| 28 | + for( ; v >= ( 1L << 28 ); v /= 10L, ++n->e_ ); |
| 29 | + n->v_ = neg ? -v : v; |
| 30 | +} |
| 31 | + |
| 32 | +static void pd_adjust( pd_t *n, int e, const int64_t *p ) |
| 33 | +{ |
| 34 | + int64_t v = n->v_; |
| 35 | + int d = n->e_ - e; |
| 36 | + if ( d > 0 ) { |
| 37 | + v *= p[ d ]; |
| 38 | + } |
| 39 | + else if ( d < 0 ) { |
| 40 | + v /= p[ -d ]; |
| 41 | + } |
| 42 | + pd_new( n, v, e ); |
| 43 | +} |
| 44 | + |
| 45 | +static void pd_mul( pd_t *r, pd_t *n1, pd_t *n2 ) |
| 46 | +{ |
| 47 | + r->v_ = n1->v_ * n2->v_; |
| 48 | + r->e_ = n1->e_ + n2->e_; |
| 49 | + pd_scale( r ); |
| 50 | +} |
| 51 | + |
| 52 | +static void pd_div( pd_t *r, pd_t *n1, pd_t *n2 ) |
| 53 | +{ |
| 54 | + if ( n1->v_ == 0 ) { pd_set( r, n1 ); return; } |
| 55 | + int64_t v1 = n1->v_, v2 = n2->v_; |
| 56 | + int neg1 = v1 < 0L, neg2 = v2 < 0L, m = 0; |
| 57 | + if ( neg1 ) v1 = -v1; |
| 58 | + if ( neg2 ) v2 = -v2; |
| 59 | + for( ; 0UL == ( ( uint64_t )v1 & 0xfffffffff0000000UL ); v1 *= 10L, ++m ); |
| 60 | + r->v_ = ( v1 * PD_SCALE9 ) / v2; |
| 61 | + if ( neg1 ) r->v_ = -r->v_; |
| 62 | + if ( neg2 ) r->v_ = -r->v_; |
| 63 | + r->e_ = n1->e_ - n2->e_ - m - 9; |
| 64 | + pd_scale( r ); |
| 65 | +} |
| 66 | + |
| 67 | +static void pd_add( pd_t *r, pd_t *n1, pd_t *n2, const int64_t *p ) |
| 68 | +{ |
| 69 | + int d = n1->e_ - n2->e_; |
| 70 | + if ( d==0 ) { |
| 71 | + pd_new( r, n1->v_ + n2->v_, n1->e_ ); |
| 72 | + } else if ( d>0 ) { |
| 73 | + if ( d<9 ) { |
| 74 | + pd_new( r, n1->v_*p[d] + n2->v_, n2->e_ ); |
| 75 | + } else if ( d < PC_FACTOR_SIZE+9 ) { |
| 76 | + pd_new( r, n1->v_*PD_SCALE9 + n2->v_/p[d-9], n1->e_-9); |
| 77 | + } else { |
| 78 | + pd_set( r, n1 ); |
| 79 | + } |
| 80 | + } else { |
| 81 | + d = -d; |
| 82 | + if ( d<9 ) { |
| 83 | + pd_new( r, n1->v_ + n2->v_*p[d], n1->e_ ); |
| 84 | + } else if ( d < PC_FACTOR_SIZE+9 ) { |
| 85 | + pd_new( r, n1->v_/p[d-9] + n2->v_*PD_SCALE9, n2->e_-9 ); |
| 86 | + } else { |
| 87 | + pd_set( r, n2 ); |
| 88 | + } |
| 89 | + } |
| 90 | + pd_scale( r ); |
| 91 | +} |
| 92 | + |
| 93 | +static void pd_sub( pd_t *r, const pd_t *n1, const pd_t *n2, const int64_t *p ) |
| 94 | +{ |
| 95 | + int d = n1->e_ - n2->e_; |
| 96 | + if ( d==0 ) { |
| 97 | + pd_new( r, n1->v_ - n2->v_, n1->e_ ); |
| 98 | + } else if ( d>0 ) { |
| 99 | + if ( d<9 ) { |
| 100 | + pd_new( r, n1->v_*p[d] - n2->v_, n2->e_ ); |
| 101 | + } else if ( d < PC_FACTOR_SIZE+9 ) { |
| 102 | + pd_new( r, n1->v_*PD_SCALE9 - n2->v_/p[d-9], n1->e_-9); |
| 103 | + } else { |
| 104 | + pd_set( r, n1 ); |
| 105 | + } |
| 106 | + } else { |
| 107 | + d = -d; |
| 108 | + if ( d<9 ) { |
| 109 | + pd_new( r, n1->v_ - n2->v_*p[d], n1->e_ ); |
| 110 | + } else if ( d < PC_FACTOR_SIZE+9 ) { |
| 111 | + pd_new( r, n1->v_/p[d-9] - n2->v_*PD_SCALE9, n2->e_-9 ); |
| 112 | + } else { |
| 113 | + pd_new( r, -n2->v_, n2->e_ ); |
| 114 | + } |
| 115 | + } |
| 116 | + pd_scale( r ); |
| 117 | +} |
| 118 | + |
| 119 | +static int pd_lt( const pd_t *n1, const pd_t *n2, const int64_t *p ) |
| 120 | +{ |
| 121 | + pd_t r[1]; |
| 122 | + pd_sub( r, n1, n2, p ); |
| 123 | + return r->v_ < 0L; |
| 124 | +} |
| 125 | + |
| 126 | +static int pd_gt( pd_t *n1, pd_t *n2, const int64_t *p ) |
| 127 | +{ |
| 128 | + pd_t r[1]; |
| 129 | + pd_sub( r, n1, n2, p ); |
| 130 | + return r->v_ > 0L; |
| 131 | +} |
| 132 | + |
| 133 | +static void pd_sqrt( pd_t *r, pd_t *val, const int64_t *f ) |
| 134 | +{ |
| 135 | + pd_t t[1], x[1], hlf[1]; |
| 136 | + pd_set( t, val ); |
| 137 | + pd_new( r, 1, 0 ); |
| 138 | + pd_add( x, t, r, f ); |
| 139 | + pd_new( hlf, 5, -1 ); |
| 140 | + pd_mul( x, x, hlf ); |
| 141 | + for(;;) { |
| 142 | + pd_div( r, t, x ); |
| 143 | + pd_add( r, r, x, f ); |
| 144 | + pd_mul( r, r, hlf ); |
| 145 | + if ( x->v_ == r->v_ ) { |
| 146 | + break; |
| 147 | + } |
| 148 | + pd_set( x, r ); |
| 149 | + } |
| 150 | +} |
| 151 | + |
| 152 | +#ifdef __cplusplus |
| 153 | +} |
| 154 | +#endif |
| 155 | + |
0 commit comments