This program can calculate FEAL's one-round characteristics with probability.

fealdiff.py

#!/usr/bin/env python                                                                                                                                                        
# -*- coding: utf-8 -*-                                                                                                                                                      
# Author: Ryoma SHINYA                                                                                                                                                       

def sfunc(a, b):
    p = bin((a + b) & 0x7F).count("1")
    a = a + b % 256
    b = (a >> 6) & 0x03
    a = (a << 2) % 256
    return (a | b, p)

def byte(x):
    return("0x%08X" % ((x[0] << 24) | (x[1] << 16) | (x[2] << 8) | x[3]))

def prob(p):
    return "1" if sum(p) == 0 else "1/"+str(2**sum(p))

def calc(dx_):
    dx = [(dx_ >> i*8) & 0xFF for i in range(4)][::-1]
    dy = [0]*4; p = [0]*4
    (dy[1], p[1]) = sfunc(dx[0]^dx[1], dx[2]^dx[3])
    (dy[0], p[0]) = sfunc(dx[0], dy[1])
    (dy[2], p[2]) = sfunc(dx[2]^dx[3], dy[1])
    (dy[3], p[3]) = sfunc(dy[2], dx[3])
    print(" input-Δ: %s" % byte(dx))
    print("output-Δ: %s with %s propability" % (byte(dy), prob(p)))

def main(arg):
    try: dx = int(arg[1], 16)
    except: exit("This program can calculate FEAL's one-round characteristics with probability.\nUsage: %s HEX" % (arg[0]))
    calc(dx)

if __name__ == "__main__":
    import sys
    main(sys.argv)