[16] | 1 | subroutine f_get_coeff(nn,mm,rr,ww,avec,nj,id,pp,np) |
---|
| 2 | !=============================================================================== |
---|
| 3 | ! Compute MT jump ahead polynomial coefficients |
---|
| 4 | ! uses GF(2)[x] computation |
---|
| 5 | !=============================================================================== |
---|
| 6 | use mt_kind_defs |
---|
| 7 | use gf2xe |
---|
| 8 | implicit none |
---|
| 9 | integer(INT32), intent(in) :: nn,mm,rr,ww,avec,nj,id |
---|
| 10 | integer(INT32), intent(inout) :: pp(0:nn-1),np |
---|
| 11 | type(gf2x_obj) :: af,bf,ff,f1,f2 |
---|
| 12 | type(gf2x_prime_obj) :: fp |
---|
| 13 | integer(INT32) :: i,ib,nws |
---|
| 14 | |
---|
| 15 | !============================== |
---|
| 16 | ! MT characteristic polynomial |
---|
| 17 | ! ff : MT char poly. |
---|
| 18 | !============================== |
---|
| 19 | call set_coef(af,nn) |
---|
| 20 | call set_coef(af,mm) ! af = x^nn + x^mm |
---|
| 21 | call set_coef(bf,nn-1) |
---|
| 22 | call set_coef(bf,mm-1) ! bf = x^(nn-1) + x^(mm-1) |
---|
| 23 | |
---|
| 24 | call pow(f1,af,ww-rr) ! f1 = af^(ww-rr) |
---|
| 25 | call pow(f2,bf,rr) ! f2 = bf^(rr) |
---|
| 26 | call mult(ff,f1,f2) ! ff = f1*f2 |
---|
| 27 | do i=0,rr-1 |
---|
| 28 | ib = mod(i,ww) |
---|
| 29 | if (BTEST(avec,ib)) then |
---|
| 30 | call pow(f2,bf,rr-1-i) |
---|
| 31 | call mult_assign(f2,f1) |
---|
| 32 | call add_assign(ff,f2) |
---|
| 33 | endif |
---|
| 34 | enddo |
---|
| 35 | do i=rr,ww-1 |
---|
| 36 | ib = mod(i,ww) |
---|
| 37 | if (BTEST(avec,ib)) then |
---|
| 38 | call pow(f1,af,ww-1-i) |
---|
| 39 | call add_assign(ff,f1) |
---|
| 40 | endif |
---|
| 41 | enddo |
---|
| 42 | |
---|
| 43 | !#ifdef _DEBUG_ |
---|
| 44 | ! write(*,'("@",$)') |
---|
| 45 | ! call print_hex(ff) |
---|
| 46 | !#endif |
---|
| 47 | |
---|
| 48 | call delete(af) |
---|
| 49 | call delete(bf) |
---|
| 50 | call delete(f1) |
---|
| 51 | call delete(f2) |
---|
| 52 | |
---|
| 53 | !============================== |
---|
| 54 | ! set ff for Barrett reduction |
---|
| 55 | !============================== |
---|
| 56 | call set_prime(fp,ff) ! fp = ff |
---|
| 57 | call delete(ff) |
---|
| 58 | call delete(f1) |
---|
| 59 | call delete(f2) |
---|
| 60 | |
---|
| 61 | !=============== |
---|
| 62 | ! jump ahead |
---|
| 63 | ! long jump |
---|
| 64 | !=============== |
---|
| 65 | call gf2x_pow_pow_2(f1,nj,fp) ! f1 = x**(2**nj) mod fp |
---|
| 66 | |
---|
| 67 | !#ifdef _DEBUG_ |
---|
| 68 | ! write(*,'("@",$)') |
---|
| 69 | ! call print_hex(f1) |
---|
| 70 | !#endif |
---|
| 71 | |
---|
| 72 | !=============== |
---|
| 73 | ! short jump |
---|
| 74 | !=============== |
---|
| 75 | call pow(ff,f1,id,fp) ! ff = f1**id mod fp |
---|
| 76 | |
---|
| 77 | !#ifdef _DEBUG_ |
---|
| 78 | ! write(*,'("@",$)') |
---|
| 79 | ! call print_hex(ff) |
---|
| 80 | !#endif |
---|
| 81 | |
---|
| 82 | pp(:) = 0 |
---|
| 83 | np = get_deg(ff)+1 |
---|
| 84 | nws = CEILING(real(np,kind=KIND(1.0d0))/32) |
---|
| 85 | pp(0:nws-1) = ff%c(0:nws-1) |
---|
| 86 | |
---|
| 87 | call delete(f1) |
---|
| 88 | call delete(f2) |
---|
| 89 | call delete(ff) |
---|
| 90 | call delete(fp) |
---|
| 91 | |
---|
| 92 | return |
---|
| 93 | end subroutine |
---|