---

下位関数群
[算術演算]

関数
size_t	ymp_addabs (digit_t *result, const digit_t *lhs, const digit_t *rhs, size_t lhs_len, size_t rhs_len)
	digit_tの列で表されている絶対値の和を求める. より詳しく...
size_t	ymp_subabs (digit_t *result, const digit_t *lhs, const digit_t *rhs, size_t lhs_len, size_t rhs_len)
	digit_tの列で表されている絶対値の差を求める. より詳しく...
digit_t	ymp_divmodabs_digit (digit_t *quotient, const digit_t *lhs, const digit_t rhs, size_t lhs_len)
	多倍長の自然数をdigit_tで割る. より詳しく...
void	do_divmodabs (digit_t *quotient, digit_t *lhs, const digit_t *rhs, size_t lhs_len, size_t rhs_len)
	多倍長自然数の除算の実行. より詳しく...
size_t	ymp_divmodabs (digit_t *quotient, digit_t *remainder, const digit_t *lhs, const digit_t *rhs, size_t lhs_len, size_t rhs_len)
	多倍長自然数の除算. より詳しく...
size_t	ymp_mulabs (digit_t *result, const digit_t *lhs, const digit_t *rhs, size_t lhs_len, size_t rhs_len)
	digit_tの列で表されている絶対値の積を求める. より詳しく...
size_t	ymp_mulabs_digit (digit_t *result, const digit_t *lhs, digit_t rhs, size_t lhs_len)
	digit_tの列で表されている絶対値にdigit_tを掛ける. より詳しく...
size_t	ymp_modmulabs (digit_t *result, const digit_t *lhs, const digit_t *rhs, const digit_t *modulus, size_t lhs_len, size_t rhs_len, size_t modulus_len)
	多倍長自然数の積の剰余を求める. より詳しく...
size_t	ymp_modmulabs_2exp (digit_t *result, const digit_t *lhs, const digit_t *rhs, size_t pow, size_t lhs_len, size_t rhs_len)
	多倍長自然数の積の、2の冪による剰余を求める. より詳しく...
size_t	ymp_modpowabs_z (digit_t *result, const digit_t *source, const digit_t *modulus, size_t source_len, size_t modulus_len, size_t pow)
	多倍長自然数の冪の剰余を求める(冪指数はsize_t). より詳しく...
size_t	ymp_modpowabs_2exp_z (digit_t *result, const digit_t *source, size_t len, size_t emod, size_t pow)
	多倍長の自然数の冪の、2の冪による剰余を求める(冪指数はsize_t). より詳しく...
size_t	ymp_modpowabs (digit_t *result, const digit_t *source, const digit_t *modulus, const digit_t *pow, size_t source_len, size_t modulus_len, size_t pow_len)
	多倍長自然数の冪の剰余を求める(冪指数は多倍長). より詳しく...
size_t	ymp_powabs (digit_t *result, const digit_t *source, size_t source_len, size_t pow)
	多倍長自然数の冪 (冪指数はsize_t). より詳しく...
size_t	ymp_modabs_2exp (digit_t *result, const digit_t *self, size_t len, size_t pow)
	多倍長自然数の2の冪乗による剰余. より詳しく...

---

関数の解説

size_t ymp_addabs (  digit_t *    result, const digit_t *    lhs, const digit_t *    rhs, size_t    lhs_len, size_t    rhs_len )  [inline]

digit_tの列で表されている絶対値の和を求める. 引数: result  結果を格納する配列の先頭へのポインタ lhs  第1項を表す配列の先頭へのポインタ rhs  第2項を表す配列の先頭へのポインタ lhs_len  lhsのポイントする先の要素数 rhs_len  rhsのポイントする先の要素数 事前条件: lhs_len >= rhs_len resultのポイントする先は最大でlhs_len+1要素を格納できること。 覚え書き: lhs, rhs, resultが一致してもよい 戻り値: resultに格納された要素数 addsub.c の 37 行で定義されています。 参照 digit_t, double_digit_t, HIGH_DIGIT, と LOW_DIGIT. { double_digit_t tmp = 0; size_t i; for ( i = 0; i < rhs_len; ++i) { tmp += lhs[i]; tmp += rhs[i]; result[i] = LOW_DIGIT(tmp); tmp = HIGH_DIGIT(tmp); } for ( ; i < lhs_len; ++i) { tmp += lhs[i]; result[i] = LOW_DIGIT(tmp); tmp = HIGH_DIGIT(tmp); } if (LOW_DIGIT(tmp)) { result[i++] = LOW_DIGIT(tmp); } return i; }

size_t ymp_subabs (  digit_t *    result, const digit_t *    lhs, const digit_t *    rhs, size_t    lhs_len, size_t    rhs_len )  [inline]

digit_tの列で表されている絶対値の差を求める. 引数: result  結果を格納する配列の先頭へのポインタ lhs  第1項を表す配列の先頭へのポインタ rhs  第2項を表す配列の先頭へのポインタ lhs_len  lhsのポイントする先の要素数 rhs_len  rhsのポイントする先の要素数 事前条件: (lhsで表される絶対値) >= (rhsで表される絶対値) resultのポイントする先は最大でlhs_len要素を格納できること。 覚え書き: lhs, rhs, resultが一致してもよい 戻り値: resultに格納された要素数 addsub.c の 81 行で定義されています。 参照 digit_t, double_digit_t, HIGH_DIGIT, LOW_DIGIT, と MULTIPREC_RADIX. { double_digit_t tmp = 1; size_t i; for ( i = 0; i < rhs_len; ++i) { tmp += MULTIPREC_RADIX - 1; tmp += lhs[i]; tmp -= rhs[i]; result[i] = LOW_DIGIT(tmp); tmp = HIGH_DIGIT(tmp); } for ( ; i < lhs_len; ++i) { tmp += MULTIPREC_RADIX - 1; tmp += lhs[i]; result[i] = LOW_DIGIT(tmp); tmp = HIGH_DIGIT(tmp); } while (i-- > 0) { if (result[i]) break; } return ++i; }

digit_t ymp_divmodabs_digit (  digit_t *    quotient, const digit_t *    lhs, const digit_t    rhs, size_t    lhs_len )

多倍長の自然数をdigit_tで割る. 引数: quotient  商を格納する配列の先頭へのポインタ NULLならば商は格納されない。 lhs  割られる数を表すdigit_tの列の先頭へのポインタ rhs  割る数 lhs_len  lhsのポイントする配列の要素数 事前条件: quotientがNULLでないなら、 そのポイント先は最大でlhs_len個の要素を格納できること rhsが0でないこと 覚え書き: quotient, lhsは一致しても良い 戻り値: 剰余 事後条件: 剰余rは 0 <= r < rhs div.c の 32 行で定義されています。 参照 digit_t, double_digit_t, LOW_DIGIT, と LOW_DIGIT_TO_HIGH_DIGIT. { double_digit_t tmp = 0; register size_t i = lhs_len; while (i-- > 0) { tmp = LOW_DIGIT_TO_HIGH_DIGIT(tmp); tmp += lhs[i]; if (quotient) quotient[i] = tmp/rhs; tmp %= rhs; } return LOW_DIGIT(tmp); }

void do_divmodabs (  digit_t *    quotient, digit_t *    lhs, const digit_t *    rhs, size_t    lhs_len, size_t    rhs_len )  [inline, static]

多倍長自然数の除算の実行. 引数: quotient  商を格納する配列の先頭へのポインタ NULLならば商は格納されない。 lhs  割られる数を表すdigit_tの列の先頭へのポインタ. 終了時には、剰余が格納される rhs  割る数を表すdigit_tの列の先頭へのポインタ lhs_len  lhsのポイントする配列の要素数 rhs_len  rhsのポイントする配列の要素数 事前条件: <rhs>が0でないこと lhs_len >= rhs_len quotientがNULLでないなら、 そのポイント先は最大でlhs_len-rhs_len+1個の要素を格納できること remainderがNULLでないなら、 そのポイント先は最大でrhs_len個の要素を格納できること quotientはlhs, rhsとは異なること 覚え書き: <rhs>が0の場合、動作は不定 事後条件: 剰余rは 0 <= r < (rhs) div.c の 72 行で定義されています。 参照 DIGIT_HIGHEST_BIT, digit_t, double_digit_t, HIGH_DIGIT, LOW_DIGIT, LOW_DIGIT_TO_HIGH_DIGIT, と MULTIPREC_RADIX. 呼出 ymp_divmodabs. { double_digit_t lhs_top; size_t i; i = lhs_len; lhs_top = lhs[--i]; while (i-- >= rhs_len) { size_t j; digit_t q; double_digit_t tmp; /* 最上位桁による仮商を定める */ lhs_top = LOW_DIGIT_TO_HIGH_DIGIT(lhs_top); lhs_top += lhs[i]; q = lhs_top / rhs[rhs_len-1]; /* (被除数)*(仮商) を被除数から引いてみる*/ tmp = 0; for (j=0; j<rhs_len; ++j) { double_digit_t tmp2 = MULTIPREC_RADIX; tmp += ((double_digit_t)rhs[j]) * q; tmp2 += lhs[i-rhs_len+j+1]; tmp2 -= LOW_DIGIT(tmp); lhs[i-rhs_len+j+1] = LOW_DIGIT(tmp2); tmp = HIGH_DIGIT(tmp); tmp += !HIGH_DIGIT(tmp2); } tmp = MULTIPREC_RADIX - tmp; tmp += lhs[i+1]; lhs[i+1] = LOW_DIGIT(tmp); /* 引き過ぎていたら、最上位桁の最上位bitは1の筈 */ while (lhs[i+1] & DIGIT_HIGHEST_BIT) { /* 商を減らし */ --q; /* その分だけ足し戻す */ tmp = 0; for (j=0; j<rhs_len; ++j) { tmp += lhs[i-rhs_len+j+1]; tmp += rhs[j]; lhs[i-rhs_len+j+1] = LOW_DIGIT(tmp); tmp = HIGH_DIGIT(tmp); } lhs[i+1] += LOW_DIGIT(tmp); } /* 得られた商を書き込む */ if (quotient) quotient[i+1-rhs_len] = q; lhs_top = lhs[i]; } }

size_t ymp_divmodabs (  digit_t *    quotient, digit_t *    remainder, const digit_t *    lhs, const digit_t *    rhs, size_t    lhs_len, size_t    rhs_len )

多倍長自然数の除算. 引数: quotient  商を格納する配列の先頭へのポインタ NULLならば商は格納されない。 remainder  剰余を格納する配列の先頭へのポインタ NULLならば剰余は格納されない。 lhs  割られる数を表すdigit_tの列の先頭へのポインタ rhs  割る数を表すdigit_tの列の先頭へのポインタ lhs_len  lhsのポイントする配列の要素数 rhs_len  rhsのポイントする配列の要素数 事前条件: <rhs>が0でないこと。0の場合はymp_fatalを呼ぶ。 lhs_used >= rhs_usedであること。そうでない場合は不定 quotientがNULLでないなら、 そのポイント先は最大でlhs_len-rhs_len+1個の要素を格納できること remainderがNULLでないなら、 そのポイント先は最大でrhs_len個の要素を格納できること 戻り値: remainderに格納された要素数。 remainderがNULLの場合は0。剰余が0であるときは、戻り値は必ず0 事後条件: 剰余rは 0 <= r < (rhs) 覚え書き: quotient, remainder, lhs, rhsは一致しても良い div.c の 159 行で定義されています。 参照 DIGIT_HIGHEST_BIT, digit_t, do_divmodabs, double_digit_t, HIGH_DIGIT, LOW_DIGIT, LOW_DIGIT_TO_HIGH_DIGIT, ymp_fatal, と YMP_TEMP_ALLOCATE. { int shift_width; /* シフト幅 */ digit_t hd; /* rhsの最上位桁 */ size_t i; digit_t *new_lhs, *new_rhs; /* 被除数、除数を各々シフトしたもの */ double_digit_t tmp; /* 0除算を弾く */ if (rhs_len==0 || (rhs_len==1 && rhs[0]==0)) { ymp_fatal("0による除算"); return 0; } /* シフトした被除数、除数を格納するための領域を確保 */ YMP_TEMP_ALLOCATE(new_lhs, digit_t, lhs_len+1); YMP_TEMP_ALLOCATE(new_rhs, digit_t, rhs_len); /* 適切なシフト幅を算出 */ shift_width = 0; hd = rhs[rhs_len-1]; while ( !(hd & DIGIT_HIGHEST_BIT) ) { ++shift_width; hd <<= 1; } /* 除数をshift_widthだけシフト */ i = rhs_len; tmp = rhs[--i]; while (i-- > 0) { tmp = LOW_DIGIT_TO_HIGH_DIGIT(tmp); tmp += rhs[i]; tmp <<= shift_width; new_rhs[i+1] = HIGH_DIGIT(tmp); tmp = rhs[i]; } new_rhs[0] = rhs[0] << shift_width; /* 被除数をshift_widthだけシフト */ tmp = 0; i=lhs_len; while(i-- > 0) { tmp += lhs[i]; tmp <<= shift_width; new_lhs[i+1] = HIGH_DIGIT(tmp); tmp = lhs[i]; tmp = LOW_DIGIT_TO_HIGH_DIGIT(tmp); } new_lhs[0] = lhs[0] << shift_width; /* シフトした状態で、除算を実行 */ do_divmodabs(quotient, new_lhs, new_rhs, lhs_len+1, rhs_len); /* 今、シフトされた剰余がnew_lhsに入っている */ if (remainder) { /* 剰余を逆シフトして復元しながら、非0の最高桁数を探す */ size_t rem_used = 0; tmp = 0; i=rhs_len; while (i-- > 0) { tmp += new_lhs[i]; remainder[i] = LOW_DIGIT(tmp >> shift_width); if (remainder[i] && !rem_used) rem_used = i+1; tmp = new_lhs[i]; tmp = LOW_DIGIT_TO_HIGH_DIGIT(tmp); } return rem_used; } else { return 0; } }

size_t ymp_mulabs (  digit_t *    result, const digit_t *    lhs, const digit_t *    rhs, size_t    lhs_len, size_t    rhs_len )  [inline]

digit_tの列で表されている絶対値の積を求める. 引数: result  結果を格納する配列の先頭へのポインタ lhs  第1項を表す配列の先頭へのポインタ rhs  第2項を表す配列の先頭へのポインタ lhs_len  lhsのポイントする先の要素数 rhs_len  rhsのポイントする先の要素数 事前条件: resultのポイントする先は最大で(lhs_len+rhs_len)要素を格納できること。 resultはlhs, rhsとは異なること 覚え書き: lhs,rhsは一致してもよい 戻り値: resultに格納された要素数 mul.c の 31 行で定義されています。 参照 digit_t, double_digit_t, HIGH_DIGIT, と LOW_DIGIT. { double_digit_t tmp=0; size_t i=0, j=0; if (lhs_len == 0 || rhs_len == 0) return 0; if ((lhs_len==1 && lhs[0]==0) || (rhs_len==1 && rhs[0]==0)) return 0; memset(result, 0, (lhs_len+rhs_len)*sizeof(digit_t)); for (i=0; i<lhs_len; ++i) { tmp = 0; for (j=0; j<rhs_len; ++j) { tmp += ((double_digit_t)lhs[i]) * rhs[j]; tmp += result[i+j]; result[i+j] = LOW_DIGIT(tmp); tmp = HIGH_DIGIT(tmp); } result[i+j] = LOW_DIGIT(tmp); } if (tmp) ++i; return i+j-1; }

size_t ymp_mulabs_digit (  digit_t *    result, const digit_t *    lhs, digit_t    rhs, size_t    lhs_len )  [inline]

digit_tの列で表されている絶対値にdigit_tを掛ける. 引数: result  結果を格納する配列の先頭へのポインタ lhs  掛けられる数を表す配列の先頭へのポインタ rhs  掛ける数 lhs_len  lhsのポイントする先の要素数 事前条件: resultのポイントする先は最大でlhs_len+1要素を格納できること。 覚え書き: lhs, resultが一致してもよい 戻り値: resultに格納された要素数 mul.c の 76 行で定義されています。 参照 digit_t, double_digit_t, HIGH_DIGIT, と LOW_DIGIT. { double_digit_t tmp = 0, tmp2, tmp3; size_t i; for (i=0; i<lhs_len; ++i) { tmp2 = ((double_digit_t)lhs[i]) * rhs; tmp3 = LOW_DIGIT(tmp); tmp3 += LOW_DIGIT(tmp2); result[i] = LOW_DIGIT(tmp3); tmp = HIGH_DIGIT(tmp); tmp += HIGH_DIGIT(tmp2); tmp += HIGH_DIGIT(tmp3); } if (LOW_DIGIT(tmp)) { result[i] = LOW_DIGIT(tmp); ++i; } return i; }

size_t ymp_modmulabs (  digit_t *    result, const digit_t *    lhs, const digit_t *    rhs, const digit_t *    modulus, size_t    lhs_len, size_t    rhs_len, size_t    modulus_len )  [inline]

多倍長自然数の積の剰余を求める. (result) := (lhs) * (rhs) % (modulus) ( ただし、(lhs), (rhs), (modulus)は自然数 ) 引数: result  結果を格納するdigit_tの列の先頭ヘのポインタ lhs  掛けられる数を表すdigit_tの列の先頭ヘのポインタ rhs  掛ける数を表すdigit_tの列の先頭ヘのポインタ modulus  法を表すdigit_tの列の先頭ヘのポインタ lhs_len  lhsのポイントする領域の要素数 rhs_len  rhsのポイントする領域の要素数 modulus_len  modulusのポイントする領域の要素数 事前条件: resultのポイントする領域は少なくともmodulus_len個の要素を格納できること 覚え書き: result, self, other, modulusは一致しても良い 戻り値: resultに格納された要素数 事後条件: 結果(result)は、0 <= (result) < (modulus) mul.c の 120 行で定義されています。 参照 digit_t, ymp_divmodabs, ymp_mulabs, と YMP_TEMP_ALLOCATE. { digit_t *product_digits; size_t product_used; YMP_TEMP_ALLOCATE(product_digits, digit_t, lhs_len+rhs_len); product_used = ymp_mulabs(product_digits, lhs, rhs, lhs_len, rhs_len); if (product_used < modulus_len) { memcpy(result, product_digits, sizeof(digit_t)*product_used); return product_used; } else { return ymp_divmodabs(NULL, result, product_digits, modulus, product_used, modulus_len); } }

size_t ymp_modmulabs_2exp (  digit_t *    result, const digit_t *    lhs, const digit_t *    rhs, size_t    pow, size_t    lhs_len, size_t    rhs_len )  [inline]

多倍長自然数の積の、2の冪による剰余を求める. (result) := (lhs) * (rhs) % (2**pow) ( ただし、(lhs), (rhs)は自然数 ) 引数: result  結果を格納するdigit_tの列の先頭ヘのポインタ lhs  掛けられる数を表すdigit_tの列の先頭ヘのポインタ rhs  掛ける数を表すdigit_tの列の先頭ヘのポインタ pow  法の冪指数 lhs_len  lhsのポイントする領域の要素数 rhs_len  rhsのポイントする領域の要素数 事前条件: resultのポイントする領域は十分な数の要素を格納できること 覚え書き: resultのポイントする領域は pow / DIGIT_BIT + 1 個を格納できれば十分 result, self, otherは一致しても良い 戻り値: resultに格納された要素数 事後条件: 結果(result)は、0 <= (result) < (rhs) mul.c の 164 行で定義されています。 参照 DIGIT_BIT, digit_t, double_digit_t, HIGH_DIGIT, LOW_DIGIT, と YMP_TEMP_ALLOCATE. { size_t n_digits, n_bits; size_t max_len; digit_t *product; double_digit_t tmp; size_t i, j; if (lhs_len == 0 || rhs_len == 0) return 0; if ((lhs_len==1 && lhs[0]==0) || (rhs_len==1 && rhs[0]==0)) return 0; n_digits = pow / DIGIT_BIT; n_bits = pow % DIGIT_BIT; max_len = n_digits + (n_bits != 0); if (result == lhs || result == rhs) { YMP_TEMP_ALLOCATE(product, digit_t, max_len); } else { product = result; } memset(product, 0, max_len*sizeof(digit_t)); if (lhs_len > max_len) lhs_len = max_len; for (i=0; i<lhs_len; ++i) { size_t bound = max_len - i; tmp = 0; for (j=0; j<rhs_len && j<bound ; ++j) { tmp += ((double_digit_t)lhs[i]) * rhs[j]; tmp += product[i+j]; product[i+j] = LOW_DIGIT(tmp); tmp = HIGH_DIGIT(tmp); } if (i+j < max_len) product[i+j] = LOW_DIGIT(tmp); } if (n_bits) product[n_digits] &= (1u << n_bits)-1; while (max_len-- > 0) { if (product[max_len]) break; } ++max_len; if (product != result) memcpy(result, product, max_len * sizeof(digit_t)); return max_len; }

size_t ymp_modpowabs_z (  digit_t *    result, const digit_t *    source, const digit_t *    modulus, size_t    source_len, size_t    modulus_len, size_t    pow )

多倍長自然数の冪の剰余を求める(冪指数はsize_t). <result> := <base>**<pow> % <modulus> 引数: result  結果を格納するdigit_tの列の先頭へのポインタ source  底を表すdigit_tの列の先頭へのポインタ modulus  法を表すdigit_tの列の先頭へのポインタ source_len  sourceのポイントする領域の長さ(要素数) modulus_len  modulusのポイントする領域の長さ(要素数) pow  冪指数 事前条件: resultはmodulus_len個の要素を格納できること 覚え書き: result, self, modulusは一致しても良い 戻り値: resultに格納された要素数 事後条件: 結果<result>は0 <= <result> < <modulus> mul.c の 368 行で定義されています。 参照 digit_t, ymp_divmodabs, ymp_modmulabs, と YMP_TEMP_ALLOCATE. { /* 繰りかえし2乗法による */ digit_t *powered_digits; /* self**(2の冪)を保持 */ size_t powered_used; /* powered_digitsの先に格納された有効要素数 */ size_t result_used; /* result == modulus の場合、複製 */ if (result == modulus) { digit_t *new_modulus; YMP_TEMP_ALLOCATE(new_modulus, digit_t, modulus_len); memcpy(new_modulus, modulus, sizeof(digit_t)*modulus_len); modulus = new_modulus; } /* <powered> = <self> ** (2**0) */ YMP_TEMP_ALLOCATE(powered_digits, digit_t, modulus_len*2); if (source_len < modulus_len) { powered_used = source_len; memcpy(powered_digits, source, source_len*sizeof(digit_t)); } else { powered_used = ymp_divmodabs(NULL, powered_digits, source, modulus, source_len, modulus_len); } /* <result> = 1 */ result_used = 1; result[0] = 1; for (;;) { if (pow & 0x01) { result_used = ymp_modmulabs(result, result, powered_digits, modulus, result_used, powered_used, modulus_len); } pow >>= 1; if (!pow) return result_used; powered_used = ymp_modmulabs(powered_digits, powered_digits, powered_digits, modulus, powered_used, powered_used, modulus_len); } }

size_t ymp_modpowabs_2exp_z (  digit_t *    result, const digit_t *    source, size_t    len, size_t    emod, size_t    pow )

多倍長の自然数の冪の、2の冪による剰余を求める(冪指数はsize_t). (result) := (source)**(pow) % (2**emod) 引数: result  結果を格納するdigit_tの列の先頭へのポインタ source  底を表すdigit_tの列の先頭へのポインタ len  sourceのポイントする領域の長さ(要素数) emod  法の冪指数 pow  冪指数 事前条件: resultは少なくとも emod / DIGIT_BIT + 1 要素を格納できること 覚え書き: result, sourceは一致しても良い 戻り値: resultに格納された要素数 mul.c の 436 行で定義されています。 参照 DIGIT_BIT, digit_t, ymp_modabs_2exp, ymp_modmulabs_2exp, と YMP_TEMP_ALLOCATE. { /* 繰りかえし2乗法による */ digit_t *powered_digits; /* self**(2の冪)を保持 */ size_t powered_used; /* powered_digitsの先に格納された有効要素数 */ size_t result_used; /* resultに格納されている要素数 */ size_t mod_len = emod / DIGIT_BIT + 1; /* <powered> = <self> ** (2**0) */ YMP_TEMP_ALLOCATE(powered_digits, digit_t, mod_len); powered_used = ymp_modabs_2exp(powered_digits, source, len, emod); /* <result> = 1 */ result_used = 1; result[0] = 1; for (;;) { if (pow & 0x01) { result_used = ymp_modmulabs_2exp(result, result, powered_digits, emod, result_used, powered_used); } pow >>= 1; if (!pow) return result_used; powered_used = ymp_modmulabs_2exp(powered_digits, powered_digits, powered_digits, emod, powered_used, powered_used); } }

size_t ymp_modpowabs (  digit_t *    result, const digit_t *    source, const digit_t *    modulus, const digit_t *    pow, size_t    source_len, size_t    modulus_len, size_t    pow_len )  [inline]

多倍長自然数の冪の剰余を求める(冪指数は多倍長). <result> := <self>**<pow> % <modulus> 引数: result  結果を格納するdigit_tの列の先頭へのポインタ source  底を表すdigit_tの列の先頭へのポインタ modulus  法を表すdigit_tの列の先頭へのポインタ pow  冪指数を表すdigit_tの列の先頭へのポインタ source_len  sourceのポイントする領域の長さ(要素数) modulus_len  modulusのポイントする領域の長さ(要素数) pow_len  powのポイントする領域の長さ(要素数) 事前条件: resultはmodulus_len個の要素を格納できること 覚え書き: result, self, modulus, powは一致しても良い 戻り値: resultに格納された要素数 事後条件: 結果<result>は0 <= <result> < <modulus> mul.c の 491 行で定義されています。 参照 DIGIT_HIGHEST_BIT, digit_t, ymp_divmodabs, ymp_modmulabs, と YMP_TEMP_ALLOCATE. { /* 繰りかえし2乗法による */ size_t i; /* pow->digits[]内のindex */ digit_t bit_mask; /* pow->digits[i]内のbitテスト用 */ digit_t *powered_digits; /* self**(2の冪)を保持 */ size_t powered_used; /* powered_digitsの先に格納された有効要素数 */ size_t result_used; /* resultの先に格納された有効要素数 */ /* result == modulus の場合、複製 */ if (result == modulus) { digit_t *new_modulus; YMP_TEMP_ALLOCATE(new_modulus, digit_t, modulus_len); memcpy(new_modulus, modulus, sizeof(digit_t)*modulus_len); modulus = new_modulus; } /* result == pow の場合、複製 */ if (result == pow) { digit_t *new_pow; YMP_TEMP_ALLOCATE(new_pow, digit_t, pow_len); memcpy(new_pow, pow, sizeof(digit_t)*pow_len); pow = new_pow; } /* <powered> = <self> ** (2**0) */ YMP_TEMP_ALLOCATE(powered_digits, digit_t, modulus_len*2); if (source_len < modulus_len) { powered_used = source_len; memcpy(powered_digits, source, sizeof(digit_t)*source_len); } else { powered_used = ymp_divmodabs(NULL, powered_digits, source, modulus, source_len, modulus_len); } /* <result> = 1 */ result_used = 1; result[0] = 1; i=0; bit_mask=1; for (;;) { if (pow[i] & bit_mask) { result_used = ymp_modmulabs(result, result, powered_digits, modulus, result_used, powered_used, modulus_len); } if (bit_mask == DIGIT_HIGHEST_BIT) { bit_mask = 1; ++i; if (i==pow_len) return result_used; } else { bit_mask <<= 1; if (i==pow_len-1 && pow[i] < bit_mask) return result_used; } powered_used = ymp_modmulabs(powered_digits, powered_digits, powered_digits, modulus, powered_used, powered_used, modulus_len); } }

size_t ymp_powabs (  digit_t *    result, const digit_t *    source, size_t    source_len, size_t    pow )  [inline]

多倍長自然数の冪 (冪指数はsize_t). <result> := <self>**<pow> 引数: result  結果を格納するdigit_tの列の先頭へのポインタ source  底を表すdigit_tの列の先頭へのポインタ source_len  sourceのポイントする領域の長さ(要素数) pow  冪指数 事前条件: resultは source_len * pow 個の要素を格納できること 覚え書き: result, selfは一致しても良い 戻り値: resultに格納された要素数 mul.c の 585 行で定義されています。 参照 digit_t, ymp_mulabs, と YMP_TEMP_ALLOCATE. { /* 繰りかえし2乗法による */ digit_t *tmp_digits; /* 一時的にdigit_t列を待避する為の領域をポイント */ digit_t *powered_digits; /* self**(2の冪)を保持 */ size_t powered_used; /* powered_digitsの先に格納された有効要素数 */ size_t result_used; if (pow == 0) { result[0] = 1; return 1; } /* 以下、<result> >= <source>を仮定できる */ /* 一時領域を確保 */ YMP_TEMP_ALLOCATE(tmp_digits, digit_t, source_len*pow); /* <powered> = <self> ** (2**0) */ YMP_TEMP_ALLOCATE(powered_digits, digit_t, source_len*pow); powered_used = source_len; memcpy(powered_digits, source, sizeof(digit_t)*source_len); /* <result> = 1 */ result_used = 1; result[0] = 1; /* 繰りかえし2乗法 */ for (;;) { if (pow & 0x01) { memcpy(tmp_digits, result, sizeof(digit_t)*result_used); result_used = ymp_mulabs(result, tmp_digits, powered_digits, result_used, powered_used); } pow >>= 1; if (!pow) return result_used; memcpy(tmp_digits, powered_digits, sizeof(digit_t)*powered_used); powered_used = ymp_mulabs(powered_digits, tmp_digits, tmp_digits, powered_used, powered_used); } }

size_t ymp_modabs_2exp (  digit_t *    result, const digit_t *    self, size_t    len, size_t    pow )

多倍長自然数の2の冪乗による剰余. <self>を2のpow乗で割った剰余をresultに格納する。 引数: result  剰余を格納するdigit_tの列の先頭へのポインタ self  元の数を表すdigit_tの列の先頭へのポインタ len  selfのポイントする領域の長さ(要素数) pow  冪指数 事前条件: resultは pow / DIGIT_BIT + (powDIGIT_BIT)?1:0個の要素を格納できること 戻り値: resultに格納された要素数 覚え書き: selfとresultは一致してもよい 格納領域の長さは最大でも pow / DIGIT_BIT + 1 あれば十分 shift.c の 120 行で定義されています。 { size_t n_digits = pow / DIGIT_BIT; size_t n_bits = pow % DIGIT_BIT; if (n_digits < len) { len = n_digits; if (n_bits) { result[n_digits] = self[n_digits] & ( (1u << n_bits) - 1u ); ++n_digits; } } else if (n_digits > len) { n_digits = len; } if (result != self) memcpy(result, self, len*sizeof(digit_t)); while (n_digits-- > 0) { if (self[n_digits] != 0) break; } return n_digits+1; }

---