[16] | 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 |
---|