Fair division in googledocs using Ruby

I figured out how to access google docs spreadsheets using Ruby's roo library.

So for the spreadsheet
https://spreadsheets.google.co/ccc?key=0AifEYKk1FmScdDUwcTE1c2c4aFF6Wl93cENkWVF1Smc&hl=en
#you can create a object for this document by

        
oo=Google.new("0AoBJp9mMcWUrdDZUY1diRm5KOTd4ZWtKWUpib2V3ZHc"
,user="stuffsplitter",password="thepassword")  

where user and password are the username and its password for the account that has the document. 0AoBJp9mMcWUrdDZUY1diRm5KOTd4ZWtKWUpib2V3ZHc is the key
no http, no google.com, no key= and no &hl=en at the end. Good thing I am not an idiot who would take hours to figure this last one out...

Then this data is sent off to the GLPK program described ad nauseum previously.
Then using ruby's standard regular expression methods the solution file fair.sol can be easily parsed. This can be easily written to the googles doc using roo.

o2.default_sheet ='Sheet3'
#match fair.sol output file of glpk. 
  if line=~/^\s+\d+ y\[(\d+).(\d+)\]\s+\*\s+(\d+)/
#   2 x[1,1]       *              0             0             1
   if $3=='1'
   o2.set_value(row, col, $1)
   o2.set_value(row, col+1, $2)
   o2.set_value(row, col+2, $3) 

Next thing is to get a web interface for people to add their own fair division problems. There are some guides for doing that here and here. So automated processing is the next task.

By the way if you do have a fair division problem you can mail it to stuffsplitter[at]gmail.com and I will run the processing for you.

Fair division using a spreadsheet and Ruby

It has taken me a month to get over the idea that someone expected this blog not to be the ranting of an inane balding drunk but actually good. Still it is time to get back to the drunken inane rambling.

There is no point having all these cool algorithms unless people can actually use them. So what is a good way to actually get the data into the program? We'll spreadsheets are used for this sort of thing all the time. Unfortunately I cannot figure out how to get openoffice, excel or gnumeric's solver to carry out this minimax fair division program.

Ruby has a cool library called roo for interacting with spreadsheets. So heres a program to get the data out of a spreadsheet and send it off to GLPK

  require 'rubygems'
  require 'roo'
  
  #open office version
 # oo = Openoffice.new("spreadsheetOO.ods")
 #  o2 = Openoffice.new("spreadsheetOO.ods")
 
 #excell version
   oo = Excel.new("spreadsheetOO.xls")
   o2 = Excel.new("spreadsheetOO.xls")
  oo.default_sheet = oo.sheets.first

 o2.default_sheet ="Sheet2"
st="" #data being written toa file
div= o2.last_column


#start writing data to model file
st<< "data;\n"
st<< "param m := #{oo.last_row-1};\n"

   st<< "param n := #{div};\n"
 #the non divisible is whatever is left
    st<< "param o := #{oo.last_column};\n"
 
st<< "param c :"
 i=1
 (oo.last_column).times{st<< "#{i} " 
 i=i+1}
 st<< ":=\n"
 
# number of people
  i=1

 2.upto(oo.last_row) do |line|  
 st<< "#{i} "
 1.upto(oo.last_column) do |col|
    st<< "#{Integer oo.cell(line,col)} "#Integer
    
   end
   st<< "\n";
   i=i+1
 end
 st<< ";\n"
  st<< "param d :"
 i=1
 (o2.last_column).times{st<< "#{i} " 
 i=i+1}
  i=1
 st<< ":=\n"
  2.upto(o2.last_row) do |line|  
 st<< "#{i} "
 1.upto(o2.last_column) do |col|
    st<<"#{Integer o2.cell(line,col)} "
 end
   st<< "\n";
   i=i+1
 end
 
 
 st<< ";\nend;"

#write data to a file
File.open("data","w")do |file|
file.write(st)
end

#call glpk program on the data
system('./glpsol -m fair.mod -d data -o fair.sol');

The input spreadsheets are openoffice and excelin the first row is the list of items.
Every other row is each persons valuation for each item. The first sheet is the indivisible items and the second sheet are the divisible ones. Unfortunately roo does not allow you write to these spreadsheets so I cannot print out the answer to a spreadsheet.

The roo interface to the google docs spreadsheet does allow cells to be written to. So getting that version up and running is the next task.

Gardner's constant

Gardner's constant is the transcendental number e^{\pi \sqrt{163}}. It is given as the integer 262,537,412,640,768,744 in his April 1975 Scientific American Column. So Lets see what programming langugaes have the accuracy to give this result?

So with these figures what will your programming languages will give Gardner's constant?

e=2.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274274663919320030599218174135966290435729003342952605956307381323286279434907632338298807531952510190115738341879307021540891499348841675092447614606680822648001684774118537423454424371075390777449920695517027618386062613313845830007520449338265602976067371132007093287091274437470472306969772093101416928368190255151086574637721112523897844250569536967707854499699679468644549059879316368892300987931277361782154249992295763514822082698951936680331825288693984964651058209392398294887933203625094431173012381970684161403970198376793206832823764648042953118023287825098194558153017567173613320698112509961818815930416903515988885193458072738667385894228792284998920868058257492796104841984443634632449684875602336248270419786232090021609902353043699418491463140934317381436405462531520961836908887070167683964243781405927145635490613031072085103837505101157477041718986106873969655212671546889570350354021234078498193343210681701210056278802351930332247450158539047304199577770935036604169973297250886876966403555707162268447162560798826517871341951246652010305921236677194325278675398558944896970964097545918569563802363701621120477427228364896134225164450781824423529486363721417402388934412479635743702637552944483379980161254922785092577825620926226483262779333865664816277251640191059004916449982893150566047258027786318641551956532442586982946959308019152987211725563475463964479101459040905862984967912874068705048958586717479854667757573205681288459205413340539220001137863009455606881667400169842055804033637953764520304024322566135278369511778838638744396625322498506549958862342818997077332761717839280349465014345588970719425863987727547109629537415211151368350627526023264847287039207643100595841166120545297030236472549296669381151373227536450988890313602057248176585118063036442812314965507047510254465011727211555194866850800368532281831521960037356252794495158284188294787610852639
pi=3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891227938183011949129833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132000568127145263560827785771342757789609173637178721468440901224953430146549585371050792279689258923542019956112129021960864034418159813629774771309960518707211349999998372978049951059731732816096318595024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177669147303598253490428755468731159562863882353787593751957781857780532171226806613001927876611195909216420198938095257201065485863278865936153381827968230301952035301852968995773622599413891249721775283479131515574857242454150695950829533116861727855889075098381754637464939319255060400927701671139009848824012858361603563707660104710181942955596198946767837449448255379774726847104047534646208046684259069491293313677028989152104752162056966024058038150193511253382430035587640247496473263914199272604269922796782354781636009341721641219924586315030286182974555706749838505494588586926995690927210797509302955321165344987202755960236480665499119881834797753566369807426542527862551818417574672890977772793800081647060016145249192173217214772350141441973568548161361157352552133475741849468438523323907394143334547762416862518983569485562099219222184272550254256887671790494601653466804988627232791786085784383827967976681454100953883786360950680064225125205117392984896084128488626945604241965285022210661186306744278622039194945047123713786960956364371917287467764657573962413890865832645995813390478027590099465764078951269468398352595709825822620522489407726719478268482601476990902640136394437455305068203496252451749399651431429809190659250937221696461515709858387410597885959772975498930161753928468138268683868942774155991855925245953959431049972524680845987273644695848653836736222626099124608051243884390451244136549762780797715691435997700129616089441694868555848406353422072225828488648158456028506016842739452267467678895252138522549954666727823986456596116354886230577456498035593634568174324112515076069479451096596094025228879710893145669136867228748940560101503308617928680920874760917824938589009714909675985261365549781893129784821682998948722658804857564014270477555132379641451523746234364542858444795265867821051141354735739523113427166102135969536231442952484937187110145765403590279934403742007310578539062198387447808478489683321445713868751943506430218453191048481005370614680674919278191197939952061419663428754440643745123718192179998391015919561814675142691239748940907186494231961567945208095146550225231603881930142093762137855956638937787083039069792077346722182562599661501421503068038447734549202605414665925201497442850732518666002132434088190710486331734649651453905796268561005508106658796998163574736384052571459102897064140110971206280439039759515677157700420337869936007230558763176359421873125147120532928191826186125867321579198414848829164470609575270695722091756711672291098169091528017350671274858322287183520935396572512108357915136988209144421006751033467110314126711136990865851639831501970165151168517143765761835155650884909989859982387345528331635507647918535
sqrt163=Math.sqrt(163)
gard=262537412640768744

guess=e**(pi*sqrt163)

puts "Gardiners constant is 262537412640768744 \n guess #{guess}\n"

Ruby does not give this result by a fair margin.

Obviously Haskell, matlab, mathematica and such are the proper languages for these problems. The Haskell compiler for OS X seems to require you to blow Satan to get it to work though*.

Your language may give a give some weird answers but are you going to believ someone called Hermite?

*Don Stewart pointed out you don't have to felate beelzbub to get Haskell. So if you don't want to get his scaley member down your gullet try here

Apéry's constant using monte Carlo methods

Apéry's Constant (1.20205690315...) "arises naturally in a number of physical problems, including in the second- and third-order terms of the electron's gyromagnetic ratio using quantum electrodynamics."

I described how to use monte carlo methods to estimate mathematical constants here and here. A passion for mathematics select three positive integers at random the odds of them having no common divisor are 1 in 1.202056....

So her is a ruby program that does not work for calculating this

def findCD(num1, num2, num3, cd)
 if num1 % cd==0 && num2 % cd==0  && num3 % cd==0
 return cd
 elsif cd < ([num1,num2,num3].max)/2
 findCD(num1, num2, num3, cd+1)
 else
 return 0
 end
end

i=1 
factor =0
N=1000

while (i < N)
 x=rand(1000)
 y=rand(1000)
 z=rand(1000)

 ans=findCD(x,y,z,2)
 if ans >0
 factor=factor+1
 end
end

ans=(factor).quo(N)

puts "Apéry's Constant estimated as #{ans}" ; 

The first problem is that my find common divisor algorithm is idiotic. There are all sorts of clever ways of doing this.

The second one is that computers have trouble creating random integers. The integers go up forever but a computers numbers don't. How do you ask a computer to give you one of any infinite one of the integers? rand(1,000,000) will not do it.

Using a Random Number Generator to estimate e

After this post I started to wonder what other constants can be estimated using probabalistic methods.

You can estimate e = 2.7182818 using
monte carlo methods. The sum S = X1 + X2 + ... + Xn > 1 averages as the eth element. So if you rand(0,1)+rand(0,1)+... until they are greater then 1. If you average the number of elements needed it approximates e.


This code is dog ugly and probably bug riddled but it might give you some idea how the calculation is done.

i=0
N=1000000
number=0
total=0  
p=0
k=1
   while k < N
    
    i=i+rand
    p=p+1
  if i>1
     total=total+p
     number=number+1
     p=0
     i=0
  end
  
    k=k+1
   end
   total=total.quo(number)
   puts "e is #{total}"

If you know of a better way to use iterators to carry out this calculation please comment.

Testing Random Number Generators with π

"Imagine that you have a dart board that is 2 units square. It inscribes a circle of unit radius. The center of the circle coincides with the center of the square. Now imagine that you throw darts at that dart board randomly. Then the ratio of the number of darts that fall within the circle to the total number of darts thrown is the same as the ratio of the area of the circle to the area of the square dart board. The area of a circle with unit radius is just π square unit. The area of the dart board is 4 square units. The ratio of the area of the circle to the area of the square is π / 4."

I got the book Digital dice the other day

It is about using monte carlo methods to solve probability problems that are hard to solve analytically. The first example given is how to estimate pi based on the % of darts that fall inside a cirle that is inside a square with the side length of twice the circles radius.

I decided to try out some languages to see how they got on with this problem. The aim is to test
1. Accuracy for mante carlo simulations
2. Ease of writing
3. Speed
in order of importance

If a language cannot accurately estimate π using this method then it indicates its random number generator (RNG) might not be very good. To test the RNG propery you need to do chi square tests and such as described in Seminumerical Algorithms by Knuth but i think its funny you can use π to test a RNG.

Perl

$N=10000000;

   $p=0;
   $k=1;
   while ($k<$N){
    $x=rand();
    $y=rand();
    if (($x*$x)+($y*$y)<1){
    $p=$p+1;
    }
      $k=$k+1;
   }
   $ans=((4*$p)/$N);
   print "ans".$ans;

Ruby

N=10000000

   p=0
   k=1
   while k < N
    x=rand
    y=rand
     if (x*x)+(y*y)<1
     p=p+1
     end
    k=k+1
   end
   ans=(4*p).quo(N)
   puts "pi is #{ans}"

Scheme
I cannot figure out how to get random numbers with scheme. There is code to do it here but I am not going to write my own rand function.

C

int main(void)
{
long N=10000000;
double x=0;
double y=0;
double p=0;
int k=1;

   while (k < N){
    x=(   (double)rand() / ((double)(RAND_MAX)+(double)(1)) );
    y=(   (double)rand() / ((double)(RAND_MAX)+(double)(1)) );
     if ((x*x)+(y*y)<1){
     p=p+1;
     }
    k=k+1;
   }
   
   double mypi=((4*p)/N);
   
   printf ("my pi  %f \n",mypi);
   return 0;
}

(less then sign) is < (thanks Helium).
I should use srand here to make it reproducible. Or more to be really fair run the test a few times and collect the average error. This is how the Digital dice runs the test with matlab. Anyway I wont tell you what language worked out the best. If you have any suggestions or code please comment.

Number Theory Hacking

Every now and again you I see a sum that I try and code up into a scripting language to see if my intuition about the problem is right.
So on page 237 of Pickover's "A passion for mathematics" is see an interesting problem.
Find x where y is an integer and x squared times 1597 plus one is an integer.
So I coded this up in Ruby as

puts "trying to find y=sqrt(1597xˆ2 + 1) problem"
y=0.1
x=1
while !y.integer?
y=Math.sqrt((1597*(x*x))+1)
x=x+1
end
puts y

The answer is a very big number. So big that using a static language is likely to lead to errors as saying int x or somesuch will lead to overflows. However ruby just keeps churning away without yet reaching an answer.
So my question is
1. How do you rapidly code up these number theory problems?
2. What language has the right mixture of rapid programming and speed of execution for doing so?