source: branches/jerome/src_flexwrf_v3.1/gf2xe.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: 24.4 KB
Line 
1module gf2xe
2!===============================================================================
3! Fortran 90/95 Module for GF(2)[x] computation
4!===============================================================================
5  use mt_kind_defs
6  implicit none
7  private
8  public :: gf2x_obj
9  public :: gf2x_prime_obj
10  public :: new,delete
11  public :: print_bit,print_hex
12  public :: get_deg
13  public :: set_coef, set_prime
14  public :: assign
15  public ::  add,  add_assign
16  public :: mult, mult_assign
17  public ::  pow, square
18  public :: div, rem, divrem
19  public :: mult_by_x, div_by_x, mod_by_x
20  public :: shift
21  public :: deg_i32, mult_i32, square_i32, shift_i32
22  public :: mult_i32_old
23  public :: gf2x_pow_pow_2
24
25  integer(INT32), parameter :: MAX_KARA  = 64
26
27  type gf2x_obj
28    integer(INT32), pointer :: c(:) => NULL()
29    integer(INT32) :: deg  = -1
30    integer(INT32) :: size = -1
31  end type
32
33  type gf2x_prime_obj
34    type(gf2x_obj) :: prime_poly
35    type(gf2x_obj) :: barrett_poly
36    integer(INT32) :: deg
37  end type
38
39  interface new
40    module procedure gf2x_new
41    module procedure gf2x_delete_prime
42  end interface
43
44  interface delete
45    module procedure gf2x_delete
46    module procedure gf2x_delete_prime
47  end interface
48
49  interface print_bit
50    module procedure gf2x_print_bit
51  end interface
52
53  interface print_hex
54    module procedure gf2x_print_hex
55  end interface
56
57  interface set_coef
58    module procedure gf2x_set_coef
59  end interface
60
61  interface set_prime
62    module procedure gf2x_set_prime
63  end interface
64
65  interface assign
66    module procedure gf2x_assign
67  end interface
68
69  interface add
70    module procedure gf2x_add
71  end interface
72
73  interface add_assign
74    module procedure gf2x_add_assign
75  end interface
76
77  interface mult
78    module procedure gf2x_mult_kara
79  end interface
80
81  interface mult_assign
82    module procedure gf2x_mult_assign_kara
83  end interface
84
85  interface pow
86    module procedure gf2x_pow
87    module procedure gf2x_pow_mod
88  end interface
89
90  interface square
91    module procedure gf2x_square
92  end interface
93
94  interface mult_by_x
95    module procedure gf2x_mult_by_x
96  end interface
97
98  interface mod_by_x
99    module procedure gf2x_mod_by_x
100  end interface
101
102  interface div_by_x
103    module procedure gf2x_div_by_x
104  end interface
105
106  interface div
107    module procedure gf2x_div
108  end interface
109
110  interface rem
111    module procedure gf2x_rem
112    module procedure gf2x_rem_barrett
113  end interface
114
115  interface divrem
116    module procedure gf2x_divrem
117  end interface
118
119  interface shift
120    module procedure gf2x_shift
121  end interface
122
123contains
124
125!!DEC$ ATTRIBUTES FORCEINLINE :: get_size
126function get_size(deg) result(size)
127  integer(INT32) :: deg,size
128  size = CEILING(real(deg+1,kind=REAL64)/32.0_REAL64)
129  return
130end function
131
132subroutine gf2x_new(this,deg)
133  type(gf2x_obj), intent(inout) :: this
134  integer(INT32), intent(in) :: deg
135  integer(INT32) :: isize
136  intrinsic :: SIZE
137  if (deg < 0) then
138    this%deg  = -1
139    this%size = -1
140    return
141  endif
142  isize = get_size(deg)
143  this%size = isize
144  this%deg  = deg
145  if (.not.associated(this%c)) then
146    allocate(this%c(0:isize-1))
147  else
148    if (SIZE(this%c) < this%size) then
149      deallocate(this%c)
150      NULLIFY(this%c)
151      allocate(this%c(0:isize-1))
152    endif
153  endif
154  this%c(:) = 0
155  return
156end subroutine
157
158subroutine gf2x_delete(this)
159  type(gf2x_obj), intent(inout) :: this
160  integer(INT32) :: ierr
161  if (associated(this%c)) then
162    deallocate(this%c,STAT=ierr)
163  endif
164  NULLIFY(this%c)
165  this%deg  = -1
166  this%size = -1
167  return
168end subroutine
169
170subroutine gf2x_print_bit(this)
171  type(gf2x_obj), intent(in) :: this
172  integer(INT32) :: i,ib,iw,deg
173  deg = get_deg(this)
174  if (deg < 0) then
175    write(*,'("0")')
176    return
177  endif
178  do i=deg,0,-1
179    ib = mod(i,32)
180    iw = i/32
181    if (BTEST(this%c(iw),ib)) then
182      write(*,'("1",$)')
183    else
184      write(*,'("0",$)')
185    endif
186  enddo
187  write(*,'("")')
188  return
189end subroutine
190
191subroutine gf2x_print_hex(this)
192  type(gf2x_obj), intent(in) :: this
193  integer(INT32) :: i,ib,iw,isize
194  character(9) :: str
195  if (is_zero(this)) then
196    write(*,'("0")')
197    return
198  endif
199  isize = get_size(this%deg)
200  i = isize-1
201  write(str,'(Z8)')this%c(i)
202  write(*,'(A,$)')TRIM(ADJUSTL(str))
203  do i=isize-2,0,-1
204    write(*,'(Z8.8,$)')this%c(i)
205  enddo
206  write(*,'("")')
207  return
208end subroutine
209
210subroutine gf2x_assign(c,a)
211  type(gf2x_obj), intent(inout) :: c  ! c := a
212  type(gf2x_obj), intent(in)    :: a
213  integer(INT32) :: ia,isa,i
214
215  call delete(c)
216  if (is_zero(a)) then
217    return
218  endif
219
220  ia = get_deg(a)
221  isa = get_size(ia)
222  call new(c,ia)
223  do i=0,isa-1
224    c%c(i) = a%c(i)
225  enddo
226 
227  return
228end subroutine
229
230function is_zero(a) result(is)
231  type(gf2x_obj), intent(in) :: a
232  logical :: is
233  integer(INT32) :: deg
234  deg = get_deg(a)
235  if (deg==-1) then
236    is = .true.
237  else
238    is = .false.
239  endif
240  return
241end function
242
243!!DEC$ ATTRIBUTES FORCEINLINE :: get_deg
244function get_deg(a) result(deg)
245  type(gf2x_obj), intent(in) :: a
246  integer(INT32) :: deg
247  integer(INT32) :: isize,i,top_deg
248  intrinsic :: SIZE
249  deg=-1
250  if (.not.associated(a%c)) return
251  isize = SIZE(a%c)
252  do i=isize-1,0,-1
253    if (a%c(i) /= 0) then
254      top_deg = deg_i32(a%c(i))
255      deg = 32*i + top_deg
256      return
257    endif
258  enddo
259  return
260end function
261
262subroutine gf2x_set_coef(a,i)
263  type(gf2x_obj), intent(inout) :: a
264  integer(INT32), intent(in) :: i
265  type(gf2x_obj), pointer :: w
266  integer(INT32) :: ib,iw
267  NULLIFY(w)
268  if (is_zero(a)) then
269    call new(a,i)
270  endif
271  allocate(w)
272  call new(w,i)
273  iw =     i/32
274  ib = mod(i,32)
275  w%c(iw) = ibset(w%c(iw),ib)
276  call add_assign(w,a)  ! w := w + a
277  call assign(a,w)      ! a := w
278  call delete(w)
279  deallocate(w)
280  NULLIFY(w)
281  a%deg = get_deg(a)
282  return
283end subroutine
284
285subroutine gf2x_add_assign(c,a)
286  type(gf2x_obj), intent(inout) :: c  ! c := c + a
287  type(gf2x_obj), intent(in)    :: a
288  type(gf2x_obj), pointer :: w
289  integer(INT32) :: ia,ic
290  integer(INT32) :: isa,isc,i
291  if (is_zero(a)) then
292    return
293  endif
294  if (is_zero(c)) then
295    call assign(c,a)
296    return
297  endif
298  ia = a%deg
299  ic = c%deg
300  isa = a%size
301  isc = c%size
302  if (isc < isa) then
303    NULLIFY(w)
304    allocate(w)
305    call new(w,MAX(ia,ic))
306    do i=0,isc-1
307      w%c(i) = IEOR(c%c(i),a%c(i))
308    enddo
309    do i=isc,isa-1
310      w%c(i) = a%c(i)
311    enddo
312    call assign(c,w)
313    call delete(w)
314    deallocate(w)
315    NULLIFY(w)
316  else
317    do i=0,isa-1
318      c%c(i) = IEOR(c%c(i),a%c(i))
319    enddo
320    c%deg = get_deg(c)
321    c%size = get_size(c%deg)
322  endif
323  return
324end subroutine
325
326subroutine gf2x_add(c,a,b)
327  type(gf2x_obj), intent(inout) :: c   ! c := a + b
328  type(gf2x_obj), intent(in)    :: a,b
329  integer(INT32) :: ia,ib,ic
330  integer(INT32) :: isa,isb,isc,i
331  if (is_zero(a) .and. is_zero(b)) then
332    return
333  endif
334  if (is_zero(a)) then
335    call assign(c,b)
336    return
337  endif
338  if (is_zero(b)) then
339    call assign(c,a)
340    return
341  endif
342  ia = get_deg(a)
343  ib = get_deg(b)
344  isa = get_size(ia)
345  isb = get_size(ib)
346  if (c%deg < MAX(ia,ib)) call new(c,MAX(ia,ib))
347  if (isa < isb) then
348    do i=0,isa-1
349      c%c(i) = IEOR(a%c(i),b%c(i))
350    enddo
351    do i=isa,isb-1
352      c%c(i) = b%c(i)
353    enddo
354  else
355    do i=0,isb-1
356      c%c(i) = IEOR(a%c(i),b%c(i))
357    enddo
358    do i=isb,isa-1
359      c%c(i) = a%c(i)
360    enddo
361  endif
362  c%deg  = get_deg(c)
363  c%size = get_size(c%deg)
364  return
365end subroutine
366
367subroutine gf2x_pow(c,a,e)
368  type(gf2x_obj), intent(inout) :: c ! c = a**e
369  type(gf2x_obj), intent(in)    :: a
370  integer(INT32), intent(in) :: e
371  type(gf2x_obj), pointer :: w
372  integer(INT32) :: ch,cl
373  integer(INT32) :: i,deg
374  NULLIFY(w)
375  call delete(c)
376  if (e==1) then
377    call assign(c,a)
378    return
379  endif
380  if (e==0) then
381    call set_coef(c,0)
382    return
383  endif
384  if (e<0) then
385    write(*,*)"pow: c = a^e : exponent should be e>=0."
386    stop
387  endif
388  if (is_zero(a)) return
389
390  deg = deg_i32(e)
391
392  allocate(w)
393  call set_coef(c,0)
394  do i=deg,0,-1
395    call square(w,c)        ! w := c**2
396    if (BTEST(e,i)) then
397      call mult(c,w,a)      ! c := w * a
398    else
399      call assign(c,w)      ! c := w
400    endif
401  enddo
402  call delete(w)
403  deallocate(w)
404  NULLIFY(w)
405
406  return
407end subroutine
408
409subroutine gf2x_square(c,a)
410  type(gf2x_obj), intent(inout) :: c ! c := a**2
411  type(gf2x_obj), intent(in)    :: a
412  integer(INT32) :: ch,cl
413  integer(INT32) :: i,deg
414  call delete(c)
415  if (is_zero(a)) return
416  deg = a%deg*2
417  call new(c,deg)
418  do i=0,a%size-1
419    if (a%c(i) == 0) cycle
420    call square_i32(a%c(i),ch,cl)
421    if (cl /= 0) c%c(2*i)   = IEOR(c%c(2*i),  cl)
422    if (ch /= 0) c%c(2*i+1) = IEOR(c%c(2*i+1),ch)
423  enddo
424  c%deg = get_deg(c)
425  c%size = get_size(c%deg)
426  return
427end subroutine
428
429recursive subroutine gf2x_mult_kara(c,a,b)
430!
431! multiply 2 polyomials using Karatsuba algorithm
432!
433  type(gf2x_obj), intent(inout) :: c    ! c := a * b
434  type(gf2x_obj), intent(in)    :: a,b
435  type(gf2x_obj), pointer :: ah,al,bh,bl,ahbh,albl,ahl,bhl,ahlbhl
436  integer(INT32) :: isa,isb,isc
437  integer(INT32) :: i,j,deg
438
439  NULLIFY(ah,al,bh,bl,ahbh,albl,ahl,bhl,ahlbhl)
440  call delete(c)
441  if (is_zero(a)) return
442  if (is_zero(b)) return
443
444  isa = a%size
445  isb = b%size
446  isc = MAX(isa,isb)
447  if (isc < MAX_KARA) then
448    call gf2x_mult_normal(c,a,b)
449    return
450  endif
451
452  if (mod(isc,2)/=0) then
453    isc = isc + 1
454  endif
455
456  allocate(ah,al,bh,bl,ahbh,albl,ahl,bhl,ahlbhl)
457  deg = 32*(isc/2)-1
458  call new(al,deg)
459  call new(bl,deg)
460  call new(ah,deg)
461  call new(bh,deg)
462
463  do i=0,MIN(isc/2-1,isa-1)
464    al%c(i) = a%c(i)
465  enddo
466  do i=0,MIN(isc/2-1,isb-1)
467    bl%c(i) = b%c(i)
468  enddo
469  do i=isc/2,isa-1
470    ah%c(i-isc/2) = a%c(i)
471  enddo
472  do i=isc/2,isb-1
473    bh%c(i-isc/2) = b%c(i)
474  enddo
475  ah%deg = get_deg(ah)
476  al%deg = get_deg(al)
477  bh%deg = get_deg(bh)
478  bl%deg = get_deg(bl)
479  ah%size = get_size(ah%deg)
480  al%size = get_size(al%deg)
481  bh%size = get_size(bh%deg)
482  bl%size = get_size(bl%deg)
483
484!===================================
485
486  call add(ahl,ah,al)
487  call add(bhl,bh,bl)
488  call gf2x_mult_kara(ahlbhl,ahl,bhl)
489  call delete(ahl)
490  call delete(bhl)
491  deallocate(ahl,bhl)
492
493!===================================
494
495  call gf2x_mult_kara(ahbh,ah,bh)
496  call delete(ah)
497  call delete(bh)
498  deallocate(ah,bh)
499
500  call add_assign(ahlbhl,ahbh)
501
502!===================================
503
504  call gf2x_mult_kara(albl,al,bl)
505  call delete(al)
506  call delete(bl)
507  deallocate(al,bl)
508
509  call add_assign(ahlbhl,albl)
510
511!===================================
512  deg = a%deg + b%deg
513  call new(c,deg)
514
515  do i=0,MIN(c%size,albl%size)-1
516    c%c(i) = albl%c(i)
517  enddo
518  call delete(albl)
519
520  if (.not. is_zero(ahlbhl)) then
521    do i=isc/2,MIN(c%size,isc/2+ahlbhl%size)-1
522      c%c(i) = IEOR(c%c(i),ahlbhl%c(i-isc/2))
523    enddo
524  endif
525  call delete(ahlbhl)
526
527  if (.not. is_zero(ahbh)) then
528    do i=isc,MIN(c%size,isc+ahbh%size)-1
529      c%c(i) = IEOR(c%c(i),ahbh%c(i-isc))
530    enddo
531  endif
532  call delete(ahbh)
533  deallocate(ahbh,albl,ahlbhl)
534  NULLIFY(ah,al,bh,bl,ahbh,albl,ahl,bhl,ahlbhl)
535  c%deg  = get_deg(c)
536  c%size = get_size(c%deg)
537  return
538end subroutine
539
540subroutine gf2x_mult_assign_kara(a,b)
541  type(gf2x_obj), intent(inout) :: a  ! a := a * b
542  type(gf2x_obj), intent(in)    :: b
543  type(gf2x_obj), pointer :: w
544  NULLIFY(w)
545  if (is_zero(a)) then
546    call delete(a)
547    return
548  endif
549  if (is_zero(b)) then
550    call delete(a)
551    return
552  endif
553  allocate(w)
554  call gf2x_mult_kara(w,a,b)
555  call assign(a,w)
556  call delete(w)
557  deallocate(w)
558  NULLIFY(w)
559  return
560end subroutine
561
562subroutine gf2x_mult_assign_normal(a,b)
563  type(gf2x_obj), intent(inout) :: a  ! a := a * b
564  type(gf2x_obj), intent(in)    :: b
565  type(gf2x_obj), pointer :: w
566  integer(INT32) :: ch,cl
567  integer(INT32) :: i,j,deg
568  NULLIFY(w)
569  allocate(w)
570  deg = a%deg + b%deg
571  call new(w,deg)
572  call gf2x_mult_normal(w,a,b)
573  call assign(a,w)  ! a := w
574  call delete(w)
575  deallocate(w)
576  NULLIFY(w)
577  return
578end subroutine
579
580subroutine gf2x_mult_normal(c,a,b)
581  type(gf2x_obj), intent(inout) :: c    ! c := a * b
582  type(gf2x_obj), intent(in)    :: a,b
583  integer(INT32) :: ch,cl
584  integer(INT32) :: i,j,ij,deg,kk,mm
585  integer(INT32), allocatable :: hi(:,:),lo(:,:)
586
587  call delete(c)
588  if (is_zero(a) .or. is_zero(b) ) then
589    return
590  endif
591
592  deg = a%deg + b%deg
593  call new(c,deg)
594
595!#define _NEW_
596!#undef _NEW_
597!#ifdef _NEW_
598!  do j=0,c%size-2
599!    kk = MIN(j,  a%size-1)
600!    mm = MAX(0,j-b%size+1)
601!    do i=mm,kk
602!      call mult_i32(a%c(i),b%c(j-i),ch,cl)
603!      c%c(j)   = IEOR(c%c(j),  cl)
604!      c%c(j+1) = IEOR(c%c(j+1),ch)
605!    enddo
606!  enddo
607!  j=c%size-1
608!  kk = a%size-1
609!  mm = c%size-b%size
610!  do i=mm,kk
611!    call mult_i32(a%c(i),b%c(j-i),ch,cl)
612!    c%c(j)   = IEOR(c%c(j),  cl)
613!  enddo
614!
615!#else
616  do j=0,b%size-1
617  if (b%c(j) == 0) cycle
618  do i=0,a%size-1
619  if (a%c(i) == 0) cycle
620
621    ij = i + j
622    call mult_i32(a%c(i),b%c(j),ch,cl)
623                     c%c(ij)   = IEOR(c%c(ij),  cl)
624    if (ij+1<c%size) c%c(ij+1) = IEOR(c%c(ij+1),ch)
625
626  enddo
627  enddo
628!#endif
629
630  c%deg = get_deg(c)
631  c%size = get_size(c%deg)
632
633  return
634end subroutine
635
636subroutine gf2x_shift(c,a,i)
637  type(gf2x_obj), intent(inout) :: c  ! c := shift(a,i)
638  type(gf2x_obj), intent(in)    :: a
639  integer(INT32), intent(in) :: i
640  integer(INT32) :: j,isn,iw,ib,ida,isa,ch,cm,cl
641  if (i==0) then
642    call assign(c,a)
643    return
644  endif
645  ida = get_deg(a)
646  isa = get_size(ida)
647  if (ida + i < 0) then
648    call delete(c)
649    return
650  endif
651  iw = abs(i)/32
652  ib = mod(abs(i),32)
653  call delete(c)
654  call new(c,ida+i)
655  if (i > 0) then
656    do j=0,isa-1
657      call shift_i32(a%c(j),+ib,ch,cm,cl)
658      if (ch /= 0) c%c(j+iw+1) = IEOR(c%c(j+iw+1),ch)
659      if (cm /= 0) c%c(j+iw)   = IEOR(c%c(j+iw)  ,cm)
660    enddo
661  else
662    call shift_i32(a%c(iw),-ib,ch,cm,cl)
663    if (cm /= 0) c%c(0)   = IEOR(c%c(0),cm)
664    do j=iw+1,isa-1
665      call shift_i32(a%c(j),-ib,ch,cm,cl)
666      if (cm /= 0) c%c(j-iw)   = IEOR(c%c(j-iw)  ,cm)
667      if (cl /= 0) c%c(j-iw-1) = IEOR(c%c(j-iw-1),cl)
668    enddo
669  endif
670  c%deg = get_deg(c)
671  c%size = get_size(c%deg)
672  return
673end subroutine
674
675subroutine gf2x_divrem(q,r,a,b)
676 ! a =: q * b + r
677  type(gf2x_obj), intent(inout) :: q  ! q := a div b
678  type(gf2x_obj), intent(inout) :: r  ! r := a mod b
679  type(gf2x_obj), intent(in)    :: a,b
680  type(gf2x_obj), pointer :: w,t,s
681  integer(INT32) :: ida,idb,idw
682  call delete(q)
683  call delete(r)
684  ida = a%deg
685  idb = b%deg
686  if (ida < idb) then
687    call assign(r,a)
688    return
689  endif
690  NULLIFY(w,t,s)
691  allocate(w,t,s)
692  call assign(w,a)
693  idw = w%deg
694  do
695    call mult_by_x(t,b,idw-idb)  ! t := b * x^(deg(w)-deg(b))
696    call set_coef(s,idw-idb)     ! s := s + x^(deg(w)-deg(b))
697    call add_assign(w,t)         ! w := w + t
698    call delete(t)
699    idw = w%deg
700    if (idw < idb) exit
701  enddo
702  call assign(r,w)
703  call delete(w)
704  call assign(q,s)
705  call delete(s)
706  deallocate(w,t,s)
707  NULLIFY(w,t,s)
708  return
709end subroutine
710
711subroutine gf2x_div(q,a,b)
712 ! a =: q * b + r
713  type(gf2x_obj), intent(inout) :: q  ! q := a div b
714  type(gf2x_obj), intent(in)    :: a,b
715  type(gf2x_obj), pointer :: w,t,s
716  integer(INT32) :: ida,idb,idw
717  call delete(q)
718  ida = a%deg
719  idb = b%deg
720  if (ida < idb) then
721    return
722  endif
723  NULLIFY(w,t,s)
724  allocate(w,t,s)
725  call assign(w,a)
726  idw = w%deg
727  do
728    call mult_by_x(t,b,idw-idb)  ! t := b * x^(deg(w)-deg(b))
729    call set_coef(s,idw-idb)     ! s := s + x^(deg(w)-deg(b))
730    call add_assign(w,t)         ! w := w + t
731    call delete(t)
732    idw = w%deg
733    if (idw < idb) exit
734  enddo
735  call delete(w)
736  call assign(q,s)
737  call delete(s)
738  deallocate(w,t,s)
739  NULLIFY(w,t,s)
740  return
741end subroutine
742
743subroutine gf2x_rem(r,a,b)
744  type(gf2x_obj), intent(inout) :: r   ! r := a mod b
745  type(gf2x_obj), intent(in)    :: a,b
746  type(gf2x_obj), pointer :: w,t
747  integer(INT32) :: ida,idb,idw
748  call delete(r)
749  ida = a%deg
750  idb = b%deg
751  if (ida < idb) then
752    call assign(r,a)
753    return
754  endif
755  NULLIFY(w,t)
756  allocate(w,t)
757  call assign(w,a)
758  idw = w%deg
759  do
760    call mult_by_x(t,b,idw-idb)  ! t := b * x^(deg(w)-deg(b))
761    call add_assign(w,t)         ! w := w + t
762    call delete(t)
763    idw = w%deg
764    if (idw < idb) exit
765  enddo
766  call assign(r,w)
767  call delete(w)
768  deallocate(w,t)
769  NULLIFY(w,t)
770  return
771end subroutine
772
773subroutine gf2x_set_prime(mp,m)
774!
775! Set a primitive polynomial to the cotainer.
776! the container contains the prime poly and precomputed polynomial for Barrett reduciont.
777! This routine does not check the primitivity.
778!  mp : container
779!   m : primitive polynomial
780!
781  type(gf2x_prime_obj), intent(inout) :: mp
782  type(gf2x_obj),       intent(in)    :: m
783  type(gf2x_obj), pointer :: xx
784  integer(INT32) :: deg
785  call delete(mp)
786  call assign(mp%prime_poly,m)
787  deg = get_deg(m)
788  mp%deg = deg
789  NULLIFY(xx)
790  allocate(xx)
791  call set_coef(xx,2*deg)
792  call div(mp%barrett_poly,xx,m)
793  call delete(xx)
794  deallocate(xx)
795  NULLIFY(xx)
796  return
797end subroutine
798
799subroutine gf2x_delete_prime(mp)
800  type(gf2x_prime_obj), intent(inout) :: mp
801  call delete(mp%prime_poly)
802  call delete(mp%barrett_poly)
803  return
804end subroutine
805
806subroutine gf2x_rem_barrett(r,a,m)
807!
808! compute  r := a mod m using Barrett algorithm
809!
810  type(gf2x_obj), intent(inout) :: r     ! r := a mod m
811  type(gf2x_obj), intent(in)    :: a
812  type(gf2x_prime_obj), intent(in) :: m  ! precomputed polynomial for Barrett algorithm
813  type(gf2x_obj), pointer :: q,p
814  integer(INT32) :: deg
815
816  call delete(r)
817  deg = m%deg
818  if (a%deg < deg) then
819    call assign(r,a)
820    return
821  endif
822
823  NULLIFY(q,p)
824  allocate(q,p)
825
826  call div_by_x(q,a,deg)              ! q = a  /  x**deg
827  call mod_by_x(r,a,deg)              ! r = a mod x**deg
828
829  call mult_assign(q,m%barrett_poly)  ! q = q  *  mu
830  call div_by_x(p,q,deg)              ! p = q  /  x**deg
831
832  call mult_assign(p,m%prime_poly)    ! p = p  *  m
833  call mod_by_x(q,p,deg)              ! q = p mod x**deg
834
835  call add_assign(r,q)                ! r = r + q
836
837  call delete(p)
838  call delete(q)
839
840  deallocate(p,q)
841  NULLIFY(p,q)
842
843  return
844end subroutine
845
846subroutine gf2x_mod_by_x(c,a,i)
847  type(gf2x_obj), intent(inout) :: c  ! c := a mod x^i
848  type(gf2x_obj), intent(in)    :: a
849  integer(INT32), intent(in) :: i
850  type(gf2x_obj), pointer :: w
851  integer(INT32) :: iw,ib,j
852  call delete(c)
853  if (a%deg < i) then
854    call assign(c,a)
855    return
856  endif
857  if (i == 0) then
858    call delete(c)
859    return
860  endif
861  if (i < 0) then
862    write(*,'("mod_by_x: error, negative i:",I10)')i
863    stop
864  endif
865  iw = i/32
866  ib = mod(i,32)
867  NULLIFY(w)
868  allocate(w)
869  call new(w,i)
870  do j=0,w%size-1
871    w%c(j) = a%c(j)
872  enddo
873  w%c(w%size-1) = IAND(w%c(w%size-1),2**ib-1)
874  call assign(c,w)
875  call delete(w)
876  deallocate(w)
877  NULLIFY(w)
878  return
879end subroutine
880
881subroutine gf2x_mult_by_x(c,a,i)
882  type(gf2x_obj), intent(inout) :: c  ! c := a * x^i
883  type(gf2x_obj), intent(in)    :: a
884  integer(INT32), intent(in) :: i
885  if (i < 0) then
886    write(*,'("mult_by_x: error, negative i:",I10)')i
887    stop
888  endif
889  if (i == 0) then
890    call assign(c,a)
891    return
892  endif
893  call shift(c,a,i)
894  return
895end subroutine
896
897subroutine gf2x_div_by_x(c,a,i)
898  type(gf2x_obj), intent(inout) :: c  ! c := a div x^i
899  type(gf2x_obj), intent(in)    :: a
900  integer(INT32), intent(in) :: i
901  if (i < 0) then
902    write(*,'("div_by_x: error, negative i:",I10)')i
903    stop
904  endif
905  if (i == 0) then
906    call assign(c,a)
907    return
908  endif
909  call shift(c,a,-i)
910  return
911end subroutine
912
913
914subroutine gf2x_pow_pow_2(c,e,m)
915  type(gf2x_obj), intent(inout) :: c  ! c := x**(2**e) mod m
916  integer(INT32), intent(in)           :: e
917  type(gf2x_prime_obj), intent(in) :: m  ! precomputed polynomial for Barrett algorithm
918  integer(INT32) :: i,ee
919  type(gf2x_obj), pointer :: w,s
920
921  ee = CEILING(log(REAL(m%deg))/log(2.0))
922  call delete(c)
923  if (ee > e) then
924    call set_coef(c,2**e)
925    return
926  endif
927
928  NULLIFY(w,s)
929  allocate(w,s)
930  call set_coef(w,2**ee)
931  call rem(s,w,m)      ! s = w mod m
932  do i=ee+1,e
933    call square(w,s)   ! w = s**2
934    call rem(s,w,m)    ! s = w mod m
935  enddo
936  call assign(c,s)
937  call delete(w)
938  call delete(s)
939  deallocate(w,s)
940  NULLIFY(w,s)
941  return
942end subroutine
943
944subroutine gf2x_pow_mod(c,a,e,m)
945  type(gf2x_obj), intent(inout) :: c  ! c := a**e mod m
946  type(gf2x_obj), intent(in)    :: a
947  integer(INT32), intent(in)           :: e
948  type(gf2x_prime_obj), intent(in) :: m  ! precomputed polynomial for Barrett algorithm
949  type(gf2x_obj), pointer :: w
950  integer(INT32) :: i,deg
951  NULLIFY(w)
952  call delete(c)
953  if (e==1) then
954    if (a%deg >= m%deg) then
955       call rem(c,a,m)
956       return
957    else
958      call assign(c,a)
959      return
960    endif
961  endif
962  if (e==0) then
963    call set_coef(c,0)
964    return
965  endif
966  if (e<0) then
967    write(*,*)"pow: c = a^e mod m : exponent should be e>=0."
968    stop
969  endif
970  if (is_zero(a)) return
971
972  deg = deg_i32(e)
973
974  allocate(w)
975  call set_coef(c,0)
976  do i=deg,0,-1
977    call square(w,c)        ! c := c**2 mod m
978    call rem(c,w,m)
979    if (BTEST(e,i)) then
980      call mult(w,c,a)      ! c := c * a mod m
981      call rem(c,w,m)
982    endif
983  enddo
984  call delete(w)
985  deallocate(w)
986  NULLIFY(w)
987  return
988end subroutine
989
990!========================================================================
991function deg_i32(a) result(d)
992  integer(INT32) :: a,d,i
993  d=-1
994  do i=31,0,-1
995    if (BTEST(a,i)) then
996      d=i
997      exit
998    endif
999  enddo
1000  return
1001end function
1002
1003function deg_i64(a) result(d)
1004  integer(INT64) :: a
1005  integer(INT32) :: d,i
1006  do i=63,0,-1
1007    if (BTEST(a,i)) then
1008      d=i
1009      exit
1010    endif
1011  enddo
1012  return
1013end function
1014
1015subroutine square_i32(a,ch,cl)
1016  integer(INT32), intent(in) :: a
1017  integer(INT32), intent(out) :: ch,cl   ! (ch,cl) = a**2
1018  integer(INT32) :: ia,i
1019  integer(INT64) :: da,dc
1020  da = a
1021  if (da < 0) da = da + 2_8**32 ! convert to unsigned
1022  dc = Z'0'
1023  ia = deg_i32(a)
1024  do i = 0,ia
1025    if (BTEST(a,i)) then
1026      dc = ibset(dc,i*2)
1027    endif
1028  enddo
1029  ch = ISHFT(dc,-32)
1030  cl = dc
1031  return
1032end subroutine
1033
1034!DEC$ ATTRIBUTES FORCEINLINE :: mult_i32
1035subroutine mult_i32(a,b,ch,cl)
1036  integer(INT32), intent(in) :: a,b
1037  integer(INT32), intent(out) :: ch,cl  ! (ch,cl) = a*b
1038  integer(INT32) :: tmp,u(0:3)
1039  integer(INT32), parameter :: ZE = Z'eeeeeeee'
1040  integer(INT32), parameter :: ZC = Z'cccccccc'
1041  integer(INT32), parameter :: Z8 = Z'88888888'
1042
1043  if (a==0 .or. b ==0) then
1044    ch = 0
1045    cl = 0
1046    return
1047  endif
1048
1049  u(0) = 0
1050  u(1) = a
1051  u(2) = ISHFT(u(1),+1)
1052  u(3) =  IEOR(u(2),a)
1053
1054  cl =                  IEOR(ISHFT(u(     ISHFT(b,-30)   ),2),u(IAND(ISHFT(b,-28),3)))
1055  ch =                  ISHFT(cl,-28)
1056  cl = IEOR(ISHFT(cl,4),IEOR(ISHFT(u(IAND(ISHFT(b,-26),3)),2),u(IAND(ISHFT(b,-24),3))))
1057  ch =  IOR(ISHFT(ch,4),ISHFT(cl,-28))
1058  cl = IEOR(ISHFT(cl,4),IEOR(ISHFT(u(IAND(ISHFT(b,-22),3)),2),u(IAND(ISHFT(b,-20),3))))
1059  ch =  IOR(ISHFT(ch,4),ISHFT(cl,-28))
1060  cl = IEOR(ISHFT(cl,4),IEOR(ISHFT(u(IAND(ISHFT(b,-18),3)),2),u(IAND(ISHFT(b,-16),3))))
1061  ch =  IOR(ISHFT(ch,4),ISHFT(cl,-28))
1062  cl = IEOR(ISHFT(cl,4),IEOR(ISHFT(u(IAND(ISHFT(b,-14),3)),2),u(IAND(ISHFT(b,-12),3))))
1063  ch =  IOR(ISHFT(ch,4),ISHFT(cl,-28))
1064  cl = IEOR(ISHFT(cl,4),IEOR(ISHFT(u(IAND(ISHFT(b,-10),3)),2),u(IAND(ISHFT(b, -8),3))))
1065  ch =  IOR(ISHFT(ch,4),ISHFT(cl,-28))
1066  cl = IEOR(ISHFT(cl,4),IEOR(ISHFT(u(IAND(ISHFT(b, -6),3)),2),u(IAND(ISHFT(b, -4),3))))
1067  ch =  IOR(ISHFT(ch,4),ISHFT(cl,-28))
1068  cl = IEOR(ISHFT(cl,4),IEOR(ISHFT(u(IAND(ISHFT(b, -2),3)),2),u(IAND(      b     ,3))))
1069
1070  tmp = -IAND(ISHFT(a,-31),1)
1071  tmp =  IAND(tmp,ISHFT(IAND(b,ZE),-1))
1072  ch  =  IEOR(ch,tmp)
1073
1074  tmp = -IAND(ISHFT(a,-30),1)
1075  tmp =  IAND(tmp,ISHFT(IAND(b,ZC),-2))
1076  ch  =  IEOR(ch,tmp)
1077
1078  tmp = -IAND(ISHFT(a,-29),1)
1079  tmp =  IAND(tmp,ISHFT(IAND(b,Z8),-3))
1080  ch  =  IEOR(ch,tmp)
1081
1082  return
1083end subroutine
1084
1085
1086subroutine mult_i32_old(a,b,ch,cl)
1087  integer(INT32), intent(in) :: a,b
1088  integer(INT32), intent(out) :: ch,cl  ! (ch,cl) = a*b
1089  integer(INT32) :: ia,ib,i
1090  integer(INT64) :: da,db,dc
1091  da = a
1092  db = b
1093  if (da < 0) da = da + 2_8**32 ! convert to unsigned
1094  if (db < 0) db = db + 2_8**32 ! convert to unsigned
1095  ia = deg_i32(a)
1096  ib = deg_i32(b)
1097  dc = Z'0'
1098  do i = 0,ia
1099    if (BTEST(a,i)) then
1100      dc = IEOR(dc,db)
1101    endif
1102    dc = ISHFTC(dc,-1)
1103  enddo
1104  dc = ISHFTC(dc,ia+1)
1105  ch = ISHFT(dc,-32)
1106  cl = dc
1107!  write(*,'(B64.64)')dc
1108!  write(*,'(B32.32)')ch
1109!  write(*,'(B64.64)')cl
1110  return
1111end subroutine
1112
1113subroutine shift_i32(a,i,ch,cm,cl)
1114  integer(INT32), intent(in) :: a
1115  integer(INT32), intent(in) :: i
1116  integer(INT32), intent(out) :: ch,cm,cl  ! (ch,cm,cl) = shift(a,i)
1117  integer(INT64) :: dc
1118  if (abs(i) >= 32) then
1119    write(*,*)"shift_int32: error i=",i
1120    stop
1121  endif
1122  select case (i)
1123  case (0)
1124    ch = 0; cm = a; cl = 0
1125    return
1126  case (1:31)
1127    dc = a
1128    if (dc < 0) dc = dc + 2_8**32 ! convert to unsigned
1129    dc = ISHFT(dc,i)
1130    ch = ISHFT(dc,-32)
1131    cm = dc
1132    cl = 0
1133    return
1134  case (-31:-1)
1135    dc = a
1136    if (dc < 0) dc = dc + 2_8**32 ! convert to unsigned
1137    dc = ISHFT(dc,i+32)
1138    ch = 0
1139    cm = ISHFT(dc,-32)
1140    cl = dc
1141    return
1142  end select
1143  return
1144end subroutine
1145
1146end module
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG