# 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
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

``````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:"
allocate(nest_limits(num_nests))
allocate(nests(num_nests))

print *, "Please enter 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
``````

and dependencies so the end of it looks like this

``````auto-examples = true
[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

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