.gitignore

pow_mod.so
types.mod

fortran_vs_c_pow.ipynb

Makefile

all: pow_mod.so

pow_mod.so: pow_mod.f90
	f2py -c -m pow_mod pow_mod.f90

pow_compare.py

from __future__ import print_function

import numpy as np

import pow_mod


NUM_POINTS = 2**15
POINTS = np.linspace(0, 1000, NUM_POINTS)


def recursive_prod(n, _cache={}):
    """Get all ways to group multiplications in p**n."""
    if n in _cache:
        return _cache[n]

    if not isinstance(n, (int, long)) or n < 1:
        raise ValueError('Expected positive integer input')
    elif n == 1:
        _cache[1] = ['p']
    else:
        # Recursion: n = small + big
        # n = small + big >= 2 small
        _cache[n] = result = []
        small = 1
        while 2 * small <= n:
            big = n - small
            for v1 in recursive_prod(small):
                v1 = '(' + v1 + ')'
                for v2 in recursive_prod(big):
                    v2 = '(' + v2 + ')'
                    result.append(v1 + ' * ' + v2)
            # Update after the loop.
            small += 1

    return _cache[n]


def comp_pow(n):
    func_name = 'pow%d' % (n,)
    print(func_name)
    pow_func = getattr(pow_mod.pow_mod, func_name)
    matches = []
    for p_expr in recursive_prod(n):
        success = True
        for p in POINTS:
            if eval(p_expr, {}, {'p': p}) != pow_func(p):
                success = False
                break
        if success:
            print('- %s' % (p_expr,))
            matches.append(p_expr)
    if not matches:
        print('No matches')


if __name__ == '__main__':
    for n in xrange(2, 16 + 1):
        comp_pow(n)

pow_compare.txt

pow2
- (p) * (p)
pow3
- (p) * ((p) * (p))
pow4
- ((p) * (p)) * ((p) * (p))
pow5
- ((p) * (p)) * ((p) * ((p) * (p)))
pow6
- ((p) * ((p) * (p))) * ((p) * ((p) * (p)))
pow7
- ((p) * ((p) * (p))) * (((p) * (p)) * ((p) * (p)))
pow8
- (((p) * (p)) * ((p) * (p))) * (((p) * (p)) * ((p) * (p)))
pow9
- ((p) * ((p) * (p))) * (((p) * ((p) * (p))) * ((p) * ((p) * (p))))
pow10
- (((p) * (p)) * ((p) * ((p) * (p)))) * (((p) * (p)) * ((p) * ((p) * (p))))
pow11
- (((p) * (p)) * ((p) * ((p) * (p)))) * (((p) * ((p) * (p))) * ((p) * ((p) * (p))))
pow12
- (((p) * ((p) * (p))) * ((p) * ((p) * (p)))) * (((p) * ((p) * (p))) * ((p) * ((p) * (p))))
pow13
- ((p) * ((p) * (p))) * ((((p) * (p)) * ((p) * ((p) * (p)))) * (((p) * (p)) * ((p) * ((p) * (p)))))
pow14
- (((p) * ((p) * (p))) * (((p) * (p)) * ((p) * (p)))) * (((p) * ((p) * (p))) * (((p) * (p)) * ((p) * (p))))
pow15
- (((p) * ((p) * (p))) * ((p) * ((p) * (p)))) * (((p) * ((p) * (p))) * (((p) * ((p) * (p))) * ((p) * ((p) * (p)))))
pow16
- ((((p) * (p)) * ((p) * (p))) * (((p) * (p)) * ((p) * (p)))) * ((((p) * (p)) * ((p) * (p))) * (((p) * (p)) * ((p) * (p))))

pow_mod.f90

module pow_mod

  implicit none
  private
  public pow2, pow3, pow4, pow5, pow6, pow7, pow8, pow9, &
         pow10, pow11, pow12, pow13, pow14, pow15, pow16

contains

  subroutine pow2(x, y)
    real(kind=8), intent(in) :: x
    real(kind=8), intent(out) :: y
    y = x**2
  end subroutine pow2

  subroutine pow3(x, y)
    real(kind=8), intent(in) :: x
    real(kind=8), intent(out) :: y
    y = x**3
  end subroutine pow3

  subroutine pow4(x, y)
    real(kind=8), intent(in) :: x
    real(kind=8), intent(out) :: y
    y = x**4
  end subroutine pow4

  subroutine pow5(x, y)
    real(kind=8), intent(in) :: x
    real(kind=8), intent(out) :: y
    y = x**5
  end subroutine pow5

  subroutine pow6(x, y)
    real(kind=8), intent(in) :: x
    real(kind=8), intent(out) :: y
    y = x**6
  end subroutine pow6

  subroutine pow7(x, y)
    real(kind=8), intent(in) :: x
    real(kind=8), intent(out) :: y
    y = x**7
  end subroutine pow7

  subroutine pow8(x, y)
    real(kind=8), intent(in) :: x
    real(kind=8), intent(out) :: y
    y = x**8
  end subroutine pow8

  subroutine pow9(x, y)
    real(kind=8), intent(in) :: x
    real(kind=8), intent(out) :: y
    y = x**9
  end subroutine pow9

  subroutine pow10(x, y)
    real(kind=8), intent(in) :: x
    real(kind=8), intent(out) :: y
    y = x**10
  end subroutine pow10

  subroutine pow11(x, y)
    real(kind=8), intent(in) :: x
    real(kind=8), intent(out) :: y
    y = x**11
  end subroutine pow11

  subroutine pow12(x, y)
    real(kind=8), intent(in) :: x
    real(kind=8), intent(out) :: y
    y = x**12
  end subroutine pow12

  subroutine pow13(x, y)
    real(kind=8), intent(in) :: x
    real(kind=8), intent(out) :: y
    y = x**13
  end subroutine pow13

  subroutine pow14(x, y)
    real(kind=8), intent(in) :: x
    real(kind=8), intent(out) :: y
    y = x**14
  end subroutine pow14

  subroutine pow15(x, y)
    real(kind=8), intent(in) :: x
    real(kind=8), intent(out) :: y
    y = x**15
  end subroutine pow15

  subroutine pow16(x, y)
    real(kind=8), intent(in) :: x
    real(kind=8), intent(out) :: y
    y = x**16
  end subroutine pow16

end module pow_mod