source: branches/jerome/src_flexwrf_v3.1/f_get_coeff.f90 @ 16

Last change on this file since 16 was 16, checked in by jebri, 11 years ago

sources for flexwrf v3.1

File size: 2.1 KB
Line 
1subroutine 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
93end subroutine
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG