1 | module 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 | |
---|
123 | contains |
---|
124 | |
---|
125 | !!DEC$ ATTRIBUTES FORCEINLINE :: get_size |
---|
126 | function get_size(deg) result(size) |
---|
127 | integer(INT32) :: deg,size |
---|
128 | size = CEILING(real(deg+1,kind=REAL64)/32.0_REAL64) |
---|
129 | return |
---|
130 | end function |
---|
131 | |
---|
132 | subroutine 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 |
---|
156 | end subroutine |
---|
157 | |
---|
158 | subroutine 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 |
---|
168 | end subroutine |
---|
169 | |
---|
170 | subroutine 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 |
---|
189 | end subroutine |
---|
190 | |
---|
191 | subroutine 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 |
---|
208 | end subroutine |
---|
209 | |
---|
210 | subroutine 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 |
---|
228 | end subroutine |
---|
229 | |
---|
230 | function 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 |
---|
241 | end function |
---|
242 | |
---|
243 | !!DEC$ ATTRIBUTES FORCEINLINE :: get_deg |
---|
244 | function 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 |
---|
260 | end function |
---|
261 | |
---|
262 | subroutine 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 |
---|
283 | end subroutine |
---|
284 | |
---|
285 | subroutine 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 |
---|
324 | end subroutine |
---|
325 | |
---|
326 | subroutine 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 |
---|
365 | end subroutine |
---|
366 | |
---|
367 | subroutine 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 |
---|
407 | end subroutine |
---|
408 | |
---|
409 | subroutine 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 |
---|
427 | end subroutine |
---|
428 | |
---|
429 | recursive 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 |
---|
538 | end subroutine |
---|
539 | |
---|
540 | subroutine 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 |
---|
560 | end subroutine |
---|
561 | |
---|
562 | subroutine 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 |
---|
578 | end subroutine |
---|
579 | |
---|
580 | subroutine 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 |
---|
634 | end subroutine |
---|
635 | |
---|
636 | subroutine 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 |
---|
673 | end subroutine |
---|
674 | |
---|
675 | subroutine 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 |
---|
709 | end subroutine |
---|
710 | |
---|
711 | subroutine 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 |
---|
741 | end subroutine |
---|
742 | |
---|
743 | subroutine 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 |
---|
771 | end subroutine |
---|
772 | |
---|
773 | subroutine 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 |
---|
797 | end subroutine |
---|
798 | |
---|
799 | subroutine 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 |
---|
804 | end subroutine |
---|
805 | |
---|
806 | subroutine 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 |
---|
844 | end subroutine |
---|
845 | |
---|
846 | subroutine 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 |
---|
879 | end subroutine |
---|
880 | |
---|
881 | subroutine 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 |
---|
895 | end subroutine |
---|
896 | |
---|
897 | subroutine 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 |
---|
911 | end subroutine |
---|
912 | |
---|
913 | |
---|
914 | subroutine 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 |
---|
942 | end subroutine |
---|
943 | |
---|
944 | subroutine 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 |
---|
988 | end subroutine |
---|
989 | |
---|
990 | !======================================================================== |
---|
991 | function 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 |
---|
1001 | end function |
---|
1002 | |
---|
1003 | function 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 |
---|
1013 | end function |
---|
1014 | |
---|
1015 | subroutine 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 |
---|
1032 | end subroutine |
---|
1033 | |
---|
1034 | !DEC$ ATTRIBUTES FORCEINLINE :: mult_i32 |
---|
1035 | subroutine 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 |
---|
1083 | end subroutine |
---|
1084 | |
---|
1085 | |
---|
1086 | subroutine 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 |
---|
1111 | end subroutine |
---|
1112 | |
---|
1113 | subroutine 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 |
---|
1144 | end subroutine |
---|
1145 | |
---|
1146 | end module |
---|