#appel: A small Maple package for implementing the Appel test
#based on Andrew Appel's paper: `Is POPL Mathematics or Science?'
 
 
print(`Version of May 10, 1997`):
print(`Written by Doron Zeilberger`):
lprint(``):
print(`appel: a small Maple package for implementing the Appel test`):
print(`based on Andrew Appel's paper: "Is POPL Mathematics or Science?"`):
print(`This paper accompanies Shalosh B. Ekhad's paper:`):
print(`The Theory of Self Avoiding Walks is Only %66.71 Mathematical`):
print(`(According to the Appel test)`):
print(`For help type " ezra();"`):
ezra:=proc():
if args=NULL then
 print(`Contains the following  procedures: p_max, Sp_max, Error_bars:`):
print(`For help with a procedure, type " ezra(procdedure_name);"`):
fi:
 
 
 
if nops([args])=1 and op(1,[args])=`p_max` then
 
print(`p_max(Lista, Listb): Given two lists Lista, and Listb`):
print(`where Lista=[a2, ..., ak], and Listb=[b2, .., bk]`):
print(`are lists of integers, where  a2 is the number of 2-author papers`):
print(` whose authors' names are listed alphabetically, a3 is the number`):
print(`of such 3-author papers etc., and Listb is the analogous list `):
print(`for papers whose authors' names is not alphabetical`):
print(`computes the Appel Science-Factor`):
print(`Note that a1 and b1 are not part of the input since it is`):
print(`irrelevant for the output`):
lprint(``):
print(`For example to get the conclusion of Ekhad's paper`):
print(`type p_max([76,22,1],[19,1,4])`):
fi:
 
 
if nops([args])=1 and op(1,[args])=`Sp_max` then
print(`A symbolic version of p_max`):
lprint(``):
print(`Sp_max(Lista, Listb): Given two lists Lista, and Listb`):
print(`where Lista=[a2, ..., ak], and Listb=[b2, .., bk]`):
print(`are lists of symbols, where  a2 signifies`):
print(` the number of 2-author papers`):
print(` whose authors' names are listed alphabetically, a3 stands for`):
print(` the number`):
print(`of such 3-author papers etc., and Listb is the analogous list `):
print(`for papers whose authors' names is not alphabetical`):
print(`Note that a1 and b1 are not part of the input since it is`):
print(`irrelevant for the output`):
lprint(``):
print(`For example to get the Appel Quotient, when there`):
print(`were only 2-author papers`):
print(`type p_max([a2],[b2])`):
 
fi:
 
 
if nops([args])=1 and op(1,[args])=`Error_bars` then
lprint(``):
print(`	Like p_max, but also computes p_lo and p_hi`):
print(`Same input at p_max, but the output is`):
print(` p_lo, p_max, p_hi`):
print(`Warning: for large inputs, very slow and memory-consuming`):
print(`you may have to set Digits to a higher number then the`):
print(` default `):
fi:
 
end:
 
 
p_max:=proc(a,b)
local k,gu,p,lu,ma,ma1,cha,p1,i,gu1:
 
if nops(a)<>nops(b) then
ERROR(`First two arguments must have the same number of elements`):
fi:
 
gu:=1:
 
for k from 2 to nops(a)+1 do
gu:=gu*(1-p+p/k!)^a[k-1]*(p-p/k!)^b[k-1]:
od:
 
gu1:=normal(diff(log(gu),p)):
gu1:=numer(gu1):
p1:=[fsolve(gu1,p)]:
 
lu:=[]:
 
for i from 1 to nops(p1) do
if op(i,p1)>0 and op(i,p1)<1 then
 lu:=[op(lu),op(i,p1)]:
fi:
od:
 
 
lu:=[op(lu),1]:
 
ma:=subs(p=0,gu):
cha:=0:
 
for i from 1 to nops(lu) do
ma1:=subs(p=op(i,lu),gu):
if ma1>ma then
 ma:=ma1:
 cha:=i:
fi:
od:
 
op(cha,lu):
end:
 
 
 
 
 
 
Sp_max:=proc(a,b)
local k,gu,p,p1,gu1:
 
if nops(a)<>nops(b) then
ERROR(`First two arguments must have the same number of elements`):
fi:
 
gu:=1:
 
for k from 2 to nops(a)+1 do
gu:=gu*(1-p+p/k!)^a[k-1]*(p-p/k!)^b[k-1]:
od:
 
gu1:=normal(diff(log(gu),p)):
gu1:=numer(gu1):
p1:=solve(gu1,p):
 
p1:
end:
 
 
L:=proc(p,a,b)
local k,gu:
 
if nops(a)<>nops(b) then
ERROR(`First two arguments must have the same number of elements`):
fi:
 
gu:=1:
 
for k from 2 to nops(a)+1 do
gu:=gu*(1-p+p/k!)^a[k-1]*(p-p/k!)^b[k-1]:
od:
 
gu:
end:
 
 
 
Error_bars:=proc(a,b)
local k,gu,p,lu,ma,ma1,cha,p1,i,gu1,T,p_max,p_lo,p_hi,gu2:
 
if nops(a)<>nops(b) then
ERROR(`First two arguments must have the same number of elements`):
fi:
 
gu:=1:
 
for k from 2 to nops(a)+1 do
gu:=gu*(1-p+p/k!)^a[k-1]*(p-p/k!)^b[k-1]:
od:
 
T:=int(gu,p=0..1):
 
 
gu1:=normal(diff(log(gu),p)):
gu1:=numer(gu1):
p1:=[fsolve(gu1,p)]:
 
lu:=[]:
 
for i from 1 to nops(p1) do
if op(i,p1)>0 and op(i,p1)<1 then
 lu:=[op(lu),op(i,p1)]:
fi:
od:
 
 
lu:=[op(lu),1]:
 
ma:=subs(p=0,gu):
cha:=0:
 
for i from 1 to nops(lu) do
ma1:=subs(p=op(i,lu),gu):
if ma1>ma then
 ma:=ma1:
 cha:=i:
fi:
od:
 
p_max:=op(cha,lu):
 
gu2:=int(gu,p=p_lo..p_max):
p_lo:=[fsolve(gu2=T/3,p_lo)]:
 
for i from 1 to nops(p_lo) do
if op(i,p_lo)<p_max and op(i,p_lo)>=0 then
 p_lo:=op(i,p_lo):
break:
fi:
od:
 
if p_max<1 then
gu2:=int(gu,p=p_max..p_hi):
p_hi:=[fsolve(gu2=T/3,p_hi)]:
 
 
 
for i from 1 to nops(p_hi) do
if op(i,p_hi)>p_max and op(i,p_hi)<=1 then
 p_hi:=op(i,p_hi):
break:
fi:
od:
else
p_hi:=1:
fi:
 
 
p_lo,p_max,p_hi
end: