Less verbose combinations

Here’s a little “code golf” exercise. Consider this code:

    program polyroots_test_10

    use polyroots_module, only: dpolz, wp => polyroots_module_rk
    use pyplot_module,    only: pyplot

    implicit none

    integer,parameter :: degree = 10

    real(wp),dimension(degree+1) :: a    !! coefficients of polynomial
    real(wp),dimension(degree) :: zr, zi !! roots
    integer :: ierr,i,j,k,l,m,n,o,p,q,r,s
    type(pyplot) :: plt

    call plt%initialize(grid=.true.,xlabel='$\Re(z)$',ylabel='$\Im(z)$',&
                        title='Degree 10 Polynomial Roots',usetex=.true.,&
                        figsize=[20,10])
    do i = -1, 1, 2
        do j = -1, 1, 2
            do k = -1, 1, 2
                do l = -1, 1, 2
                    do m = -1, 1, 2
                        do n = -1, 1, 2
                            do o = -1, 1, 2
                                do p = -1, 1, 2
                                    do q = -1, 1, 2
                                        do r = -1, 1, 2
                                            do s = -1, 1, 2
                                                a = [i,j,k,l,m,n,o,p,q,r,s]
                                                call dpolz(degree,a,zr,zi,ierr); if (ierr/=0) error stop ierr
                                                call plt%add_plot(zr,zi,label='',linestyle='bo',markersize=1)
                                            end do
                                        end do
                                    end do
                                end do
                            end do
                        end do
                    end do
                end do
            end do
        end do
    end do
    call plt%savefig('roots.png')

    end program polyroots_test_10

The resultant plot can be seen here.

The idea is to get all the combinations where each element of a is either 1 or -1. How can we make this less verbose? Is there an intrinsic way I’m missing, or a library out there somewhere to do this?

2 Likes

how about:

do concurrent (I=0:2**(degree+1)-1)    
   a = [(merge(1,-1,btest(I,pos)),pos=0,degree)]
   call blabla...
end do

I may be stretching the standard and compilers a bit here… in other words, it’s untested (a check on machine endianness may be needed).

2 Likes

The code of chw21 at recursion - What are ways in which an arbitrary number of nested loops be created in Fortran? - Stack Overflow handles the more general case where loop variables take on a specified number of values.

program nested
    implicit none
    integer :: num_nests, i
    integer, dimension(:), allocatable :: nest_limits
    integer, dimension(:), allocatable :: nests

    print *, "Please enter number of nests:"
    read(*, *) num_nests
    allocate(nest_limits(num_nests))
    allocate(nests(num_nests))

    print *, "Please enter nest limits:"
    read(*, *) nest_limits

    nests(:) = 1
    outer_loop : do
        print *, nests(:)
        i = 1
        ! Calculate the next indices:
        inner_loop : do  
            nests(i) = nests(i) + 1 

            ! If this is still a valid index, exit the inner 
            ! loop and go for the next iteration
            if (nests(i) <= nest_limits(i)) exit inner_loop

            ! The index has overflown, so reset it to 1 and
            ! move to next index.
            nests(i) = 1
            i = i + 1

            ! If the next index would be outside of num_nests, 
            ! the whole loop is finished.
            if (i > num_nests) exit outer_loop

        end do inner_loop
    end do outer_loop
end program nested
1 Like

all the permutations of eleven sets of {-1,1} ? If I read that right (did not check the math though)
I think this might be simpler

     integer,parameter :: degree = 10
    integer,parameter :: wp=kind(0.0)
    real(wp),dimension(degree+1) :: a    !! coefficients of polynomial
    integer :: i,j
    do i=0,2**(degree+1)-1
         a=[(merge(1,-1,btest(i,j)),j=0,degree)]
         write(*,*)'do something with ',a
    enddo
    end

A quick test looks right

urbanjs@mercury:/tmp$ gfortran xx.f90
urbanjs@mercury:/tmp$ ./a.out|wc
   2048   22528  385024
urbanjs@mercury:/tmp$ ./a.out|sort|uniq|wc
   2048   22528  385024

Hah, in the mean time @FedericoPerini did the same thing potentially in parallel. Ah well,
but I did it with M_pixel and made an image all in Fortran (but made me too slow!)

@Beliavsky:
That recursive nesting has me thinking about how to do permutations in a way I never thought about before. Old habits die hard. I have avoided recursion like the plague for a long time for various reasons, many of which are probably no longer pertinent, but this is interesting. Thanks!

@FedericoPerini

The DO CONCURRENT would be great so I was changing my try; and got a message about the
polyroots_module procedure not being PURE. I was really curious what timing differences it would have over a vanilla DO on various platforms too. I have had mixed results in the past.

When I think of how casually I can look at a post like this
on Fortran Discourse and using all Fortran code and fpm(1)
experiment with the polyroots code and display some results
in a pixmap on a Linux box (showing ubuntu but tried OpenBSD
as well and it took two additional minutes) I realize thinks
have come a long way.

# install some (ubuntu) packages if you do not have them
sudo apt-get install libblas-dev liblapack-dev
sudo apt install imagemagick-6.q16
# make a little fpm program on the fly:
fpm new roots
cd roots

replace app/main.f90 with a modified version of the
original

    program polyroots_test_10
    use polyroots_module, only: dpolz, wp => polyroots_module_rk
    use :: M_pixel, only : prefsize,vinit,ortho2,vexit
    use :: M_pixel, only : linewidth,circle,color,point2
    use :: M_pixel, only : print_p6
    implicit none
    integer,parameter :: degree = 10
    real(wp),dimension(degree+1) :: a    !! coefficients of polynomial
    real(wp),dimension(degree) :: zr, zi !! roots
    integer :: i,j,pp
    call prefsize(700,700)
    call vinit()
    call ortho2(left=-2.0, right=2.0, bottom=-2.0, top=2.0)
    call color(6)
    call linewidth(80)
    call circle(0.0,0.0,2.0)
    call color(2)
    do i=0,2**(degree+1)-1
       a=[(merge(1,-1,btest(i,j)),j=0,degree)]
       call dpolz(degree,a,zr,zi,ierr); if (ierr/=0) error stop ierr
       do pp=1,degree
          call point2(real(zr(pp)),real(zi(pp)))
       enddo
    enddo
    call print_p6('roots.p6')
    call vexit()
    end program polyroots_test_10

edit fpm.toml and add blas and lapack links
and dependencies so the end of it looks like this

auto-examples = true
link = ["blas", "lapack"] 
[install]
library = false
[dependencies]
polyroots-fortran = { git="https://github.com/jacobwilliams/polyroots-fortran.git" }
M_pixel = { git="https://github.com/urbanjost/M_pixel.git" }

build and run program and display pixel map

fpm run
display roots.p6

PS:

I stuck the little fpm golf project into a tarball for anyone wanting a seed for similar methodology using fpm, fodder and all (temporarily):

http://www.urbanjost.altervista.org/REMOVE/roots.tgz

assuming your system can read a tar(1) file.

7 Likes

Here’s another solution that uses recursion, inspired by the Permutations entry in Rosetta Code:

program polyroots_test_10

use polyroots_module, only: dpolz, wp => polyroots_module_rk
use pyplot_module,    only: pyplot

implicit none

integer,parameter :: degree = 10 !! polynomial degree
integer,dimension(*),parameter :: icoeffs = [-1,1] !! set of coefficients
integer,dimension (degree+1) :: a !! coefficients of polynomial
real(wp),dimension(degree) :: zr, zi !! roots
integer :: ierr !! error code from [[dpolz]]
type(pyplot) :: plt !! for making the plot

call plt%initialize(grid=.true.,xlabel='$\Re(z)$',ylabel='$\Im(z)$',&
                    title='Degree 10 Polynomial Roots',usetex=.true.,&
                    figsize=[20,10])

call generate(1)
call plt%savefig('roots.png')

contains

    recursive subroutine generate (i)
    integer, intent(in) :: i
    integer :: ix
    if (i > degree+1) then
        call dpolz(degree,real(a,wp),zr,zi,ierr); if (ierr/=0) error stop ierr
        call plt%add_plot(zr,zi,label='',linestyle='bo',markersize=1)
    else
        do ix = 1,size(icoeffs) 
            a(i) = icoeffs(ix)
            call generate(i+1)
        end do
    end if
    end subroutine generate

end program polyroots_test_10

What we need is a general Fortran library for permutations and combinations…