
      real*8 xguess(18),xres(18),fv(18),xsc(18),fsc(18)
      real*8 pnew(2,2),aux,auxi,uni(1000000,4,4)
      real*8 err,thumb,dbetin,jac,fscold(30),erold
      real*8 rp(5),moy,std,test
      real*8 ero,fn,vec(1000000)
      real*8 moy1,moy2,std1,std2
      integer ip(6),m,i,j,l,jj,t,max,int,iperm(1000000),auxa
      character rep

      real*8 fscc(30),bestfsc(30)
      real*8 xnew(30),error,xx(30,4),erre,h(4),u(1000000,4,4)
      integer mc,nj,no,np(4),from
      common fscc,bestfsc,xnew,error,xx,erre,h,u,mc,nj,no,np
      common from

      external dn4qbj,erset,thumb,drnun,rnset
      external dneqbj,jac,dneqnf,dsvrgp,dneqbf
      intrinsic dlog,dexp

      open(unit=60,file='data_multi_10_5.dat',status='old')



      call erset(0,0,0)


      mc=100000
      nj=4
      no=4

      moy1=0.0d0
      std1=0.0d0
      moy2=0.0d0
      std2=0.0d0

      do 310 m=1,mc
      do 310 i=1,nj
            read(60,100)u(m,i,1)
            moy1=moy1+u(m,i,1)
            read(60,100)u(m,i,2)
            moy2=moy2+u(m,i,2)
  310 continue

      moy1=moy1/(m*nj*1.0d0)
      moy2=moy2/(m*nj*1.0d0)
      do 315 m=1,mc
      do 315 i=1,nj
         std1=std1+(u(m,i,1)-moy1)**2
         std2=std2+(u(m,i,2)-moy2)**2
  315 continue
      std1=std1/(m*nj*1.0d0)
      std2=std2/(m*nj*1.0d0)

      write(6,*)"moy et std de t1",moy1,std1**0.5
      write(6,*)"moy et std de t2",moy2,std2**0.5
      write(6,*)' '
      read *,rep



c guess for the roots of the system of equations
c and set up the parameters to call the IMSL subroutine
c to solve a system of equations


       xres(  1)=   0.48888723155924
       xres(  2)=   0.52388476497143
       xres(  3)=  -0.19511960471946
       xres(4)=0.0d0
       xres(5)=0.0d0
       xres(6)=0.001d0
       xres(7)=0.0d0
       xres(8)=0.0d0

       xres(  1)=   0.44844796781484
       xres(  2)=   0.97913057793063
       xres(  3)=  -0.15168519369410
       xres(  4)=    7.7354458106526D-02
       xres(  5)=   -2.0192674732847D-02
       xres(  6)=  -0.52670318386583
       xres(  7)=    1.6134439267999
       xres(  8)=  -0.43276470755771
       xres(  9)=  -0.15678766709121
       xres(  10)=  -0.18394100525232
       xres(  11)=   0.24981759575368
       xres(  12)=    1.0018481931040D-01
       xres(  13)=   0.69596317683877
       xres(  14)=    3.3547592974761D-02
       xres(  15)=   -6.9543928473605D-02
       xres(  16)=   -1.3298128066624
       xres(17)=0.0d0
       xres(18)=0.0d0

       xres(  1)=   0.44175298277314
       xres(  2)=  -0.43939741989722
       xres(  3)=  -0.35120357603295
       xres(  4)=    3.4442269623419
       xres(  5)=   0.51702743573091
       xres(  6)=  -0.55154808509248
       xres(  7)=    1.5585275998857
       xres(  8)=  -0.44167284107514
       xres(  9)=   -3.5440180319235
       xres(  10)=  -0.78551365879282
       xres(  11)=   0.26413562548741
       xres(  12)=   0.15612519013592
       xres(  13)=   0.70489913800627
       xres(  14)=    3.3972017981304D-02
       xres(  15)=   -7.1644539057936D-02
       xres(  16)=   -1.3519360994986
       xres(  17)=    1.2350641291303
       xres(  18)=   0.23762480724456

       xres(  1)=   0.43978000405998
       xres(  2)=  -0.99759240472291
       xres(  3)=  -0.48443818141038
       xres(  4)=    4.6699056999151
       xres(  5)=   0.84888932845727
       xres(  6)=  -0.54948359690563
       xres(  7)=    1.7778900688471
       xres(  8)=  -0.42733586471282
       xres(  9)=   -4.7200801017858
       xres(  10)=   -1.1196473677097
       xres(  11)=   0.25845986142656
       xres(  12)=    2.7789711493896D-02
       xres(  13)=   0.70091170680086
       xres(  14)=    3.5639812840083D-02
       xres(  15)=   -7.3386963592554D-02
       xres(  16)=   -1.3846463245458
       xres(  17)=    1.6564753458449
       xres(  18)=   0.35694261998576

       xres(  1)=   0.44595088474186
       xres(  2)=   0.84731557527759
       xres(  3)=   0.10740494740541
       xres(  4)=    8.9955702455344D-02
       xres(  5)=  -0.45994603457003
       xres(  6)=  -0.56117683980651
       xres(  7)=    1.9178692061159
       xres(  8)=  -0.55015271217550
       xres(  9)=   0.15173172715340
       xres(  10)=   0.22080433843579
       xres(  11)=   0.26411249390874
       xres(  12)=    2.5008091233279D-02
       xres(  13)=   0.77776637283058
       xres(  14)=    3.5665320213955D-02
       xres(  15)=   -7.6368132387512D-02
       xres(  16)=   -1.4991483186960
       xres(  17)=  -0.21422813441816
       xres(  18)=  -0.15188821789723
      

       xres(  1)=   0.44837077957651
       xres(  2)=    1.9045763618076
       xres(  3)=   -5.6977626213538D-02
       xres(  4)=   -2.8218590797898
       xres(  5)=  -0.33966069721605
       xres(  6)=  -0.63672087499752
       xres(  7)=    2.2288140036787
       xres(  8)=  -0.73474860070130
       xres(  9)=    3.6326510716698
       xres(  10)=   0.33946813710292
       xres(  11)=   0.30701098216466
       xres(  12)=   0.16185517968072
       xres(  13)=   0.88946017648150
       xres(  14)=    3.6279819952376D-02
       xres(  15)=   -7.9348683493883D-02
       xres(  16)=   -1.8413765516446
       xres(  17)=   -1.6908731024637
       xres(  18)=  -0.26041672329279

       xres(  1)=   0.44195612021121
       xres(  2)=    1.2722410118649
       xres(  3)=  -0.55144776277504
       xres(  4)=  -0.97739320471062
       xres(  5)=   0.83554750465512
       xres(  6)=  -0.59635033506932
       xres(  7)=    1.9085837458700
       xres(  8)=  -0.61571062139810
       xres(  9)=    1.3593309663153
       xres(  10)=  -0.95048748389763
       xres(  11)=   0.28091332750461
       xres(  12)=   0.23751811722147
       xres(  13)=   0.81215813600690
       xres(  14)=    3.7210155613757D-02
       xres(  15)=   -7.5815049371150D-02
       xres(  16)=   -1.6814106225720
       xres(  17)=  -0.70704959229871
       xres(  18)=   0.25810068673036


      auxa=18
      from=0
      call rnset(806)
      h=5.6814189053753D-02

c we enter loop that verifies and guarantee convergence
c i.e. err very very small

      erre=10000.0d0
      error=200.0d0
      erold=200.0d0

      do 445 t=1,auxa
         xnew(t)=xres(t)
  445 continue


  943 if (erre.gt.0.00000000000001d0) then


      if (erre.gt.erold) then
         erold=erre
         t=0
         do 822 t=1,auxa
            xnew(t)=xres(t)
            xguess(t)=xnew(t)
            call drnnor(1,auxi)
            xguess(t)=xnew(t)
c     &             *(0.9999999995+auxi/10000000000.0d0)

            xsc(t)=1.0d0
            fsc(t)=0.1d0
            fscc(t)=fsc(t)
c             fsc(t)=fscold(t)
  822    continue
      else
         do 823 t=1,auxa
            call drnun(1,auxi)
            xguess(t)=xnew(t)
     &           *(0.99999999999995+auxi/10000000000000.0d0)
c            xguess(t)=xnew(t)
            xsc(t)=1.0d0
            fsc(t)=bestfsc(t)
     &            *(0.5d0+1.0*auxi)
            fscc(t)=fsc(t)
c            fscold(t)=fsc(t)
  823    continue
      endif



      call dn4qbj(ip,rp)
      ip(3)=500
      ip(4)=1000
      ip(6)=1

c      fscc=fsc(1)
c      ip(1)=0

c      call dneqbj(thumb,jac,no*np,xguess,xsc,fsc,ip,rp,xres,fv)
c      call dneqbf(thumb,auxa,xguess,xsc,fsc,ip,rp,xres,fv)
      max=50
      ero=0.00000001d0
c      write(6,*)ero,auxa,max,xguess,xres,fn
      call dneqnf(thumb,ero,auxa,max,xguess,xres,fn)

      write(6,*)' '
      t=0
      do 627 j=1,no
      do 626 l=1,np(j)
      t=t+1
      write(6,*)'xres,xnew',j,l,dexp(xres(t)),
     &          dexp(xnew(t))
c         write(6,*)'xres=',j,l,dexp(xnew(l+(j-1)*np))
c     &             /(1.0d0+dexp(xnew(l+(j-1)*np)))
  626 continue
      write(6,*)' '
  627 continue

      write(6,*)'fsc',fsc(1)
      write(6,*)'best error',error
      write(6,*)'erreur',erre,fn
      write(6,*)' '
      write(6,*)' '
      write(6,*)'stop'

c      read*,rep


      goto 943
      endif

  100 format (2d20.12)
      end


c****************************************************************
*********
c subroutine that calculates the objective function


c      subroutine thumb(numb,ix,f)
      subroutine thumb(ix,f,numb)


c numb=# of roots, x=point at which the function f(.) is evaluate
d

      integer numb
      real*8 x(46),ix(numb),f(numb)

      real*8 pemp(4),profit(4),proba(4),highb(5),para(4,4)
      real*8 dt1,dt2,du,da(46),aux,b1,db1,ks,dks
      real*8 auxi,bid(101,4),diff,alph(46,4),bet(46,4)
      real*8 t1,t2

      integer t,m,i,j,l,jj,ll
      character rep

      real*8 fscc(30),bestfsc(30)
      real*8 xnew(30),error,xx(30,4),erre,h(4),u(1000000,4,4)
      integer mc,nj,no,np(4),from
      common fscc,bestfsc,xnew,error,xx,erre,h,u,mc,nj,no,np
      common from



c  set the value of the constrained startegy parameters

       do 456 t=1,numb
          x(from+t)=(ix(t))
          da=0.0d0
  456  continue

      do 823 j=1,no
         profit(j)=0.0d0
         proba(j)=0.0d0
         pemp(j)=0.0d0
  823 continue


      do 898 m=1,mc

c calculate the k highest bids by bidder 1 opponents

         do 871 j=1,no
            highb(j)=-1000.0d0
 871     continue

         do 895 i=2,nj
         do 895 j=1,no
            t1=u(m,i,1)
            t2=u(m,i,2)
            aux=x(1)+x(2)*(1-t1)+x(3)*(1-t2)
     &      +x(4)*(1-t1**2)+x(5)*(1-t2**2)
     &      +x(6)*j*(1-t1)+x(7)*(1-t1)*(1-t2)
     &      +x(8)*j*(1-t2)
     &      +x(9)*(1-t1**3)+x(10)*(1-t2**3)
     &      +x(11)*j*(1-t1**2)+x(12)*(1-t1**2)*(1-t2)
     &      +x(13)*j*(1-t2**2)+x(14)*(1-t1)*j**2
     &      +x(15)*(1-t2)*j**2+x(16)*(1-t1)*(1-t2**2)
     &      +x(17)*(1-t1**4)+x(18)*(1-t2**4)
            aux=aux*t1*t2**(j-1)
            jj=no
  893       if (jj.ge.1.and.aux.gt.highb(jj)) then
                highb(jj+1)=highb(jj)
                highb(jj)=aux
                jj=jj-1
                goto 893
            endif
  895    continue

c         do 770 j=1,no
c            write(6,*)"highb",j,highb(j)
c  770    continue
c         read *,rep


c calculte player 1 bids and objcetive functio

         t=0
         do 898 j=1,no
            t1=u(m,1,1)
            t2=u(m,1,2)
            b1=x(1)+x(2)*(1-t1)+x(3)*(1-t2)
     &      +x(4)*(1-t1**2)+x(5)*(1-t2**2)
     &      +x(6)*j*(1-t1)+x(7)*(1-t1)*(1-t2)
     &      +x(8)*j*(1-t2)
     &      +x(9)*(1-t1**3)+x(10)*(1-t2**3)
     &      +x(11)*j*(1-t1**2)+x(12)*(1-t1**2)*(1-t2)
     &      +x(13)*j*(1-t2**2)+x(14)*(1-t1)*j**2
     &      +x(15)*(1-t2)*j**2+x(16)*(1-t1)*(1-t2**2)
     &      +x(17)*(1-t1**4)+x(18)*(1-t2**4)
            b1=b1*t1*t2**(j-1)
c            ks=1.0d0/(1.0d0+dexp(-((b1-highb(no-j+1))/h(j))))
            aux=(b1-highb(no-j+1))/h(j)
            ks=0.0d0
            dks=0.0d0
            if (aux.gt.1.0d0) then
               ks=1.0d0
            else if (aux.ge.-1.0d0) then
               ks=0.5d0+(315.0d0/2048.0d0)*(15.0d0*aux-(140.0d0/3.0d0)
     &         *aux**3+(378.0d0/5.0d0)*aux**5-(396.0d0/7.0d0)*
     &         aux**7+(143.0d0/9.0d0)*aux**9)
               dks=(315.0d0/2048.0d0)*(15.0d0-140.0d0
     &             *aux**2+378.0d0*aux**4-396.0d0*
     &             aux**6+143.0d0*aux**8)
            endif

            if (b1.gt.highb(no-j+1)) then
               pemp(j)=pemp(j)+1.0d0
            endif
            proba(j)=proba(j)+ks
            profit(j)=profit(j)+(u(m,1,1)*u(m,1,2)**(j-1)-b1)*ks

            do 603 l=1,numb
               if (l.eq.1) then
                  db1=1.0d0
               endif
               if (l.eq.2) then
                  Db1=1-t1
               endif
               if (l.eq.3) then
                  db1=1-t2
               endif
               if (l.eq.4) then
                  db1=1-t1**2
               endif
               if (l.eq.5) then
                  db1=1-t2**2
               endif
               if (l.eq.6) then
                  db1=j*(1-t1)
               endif
               if (l.eq.7) then
                  db1=(1-t1)*(1-t2)
               endif
               if (l.eq.8) then
                  db1=j*(1-t2)
               endif
               if (l.eq.9) then
                  db1=1-t1**3
               endif
               if (l.eq.10) then
                  db1=1-t2**3
               endif
               if (l.eq.11) then
                  db1=j*(1-t1**2)
               endif
               if (l.eq.12) then
                  db1=(1-t1**2)*(1-t2)
               endif
               if (l.eq.13) then
                  db1=j*(1-t2**2)
               endif
               if (l.eq.14) then
                  db1=(1-t1)*j**2
               endif
               if (l.eq.15) then
                  db1=(1-t2)*j**2
               endif
               if (l.eq.16) then
                  db1=(1-t1)*(1-t2**2)
               endif
               if (l.eq.17) then
                  db1=1-t1**4
               endif
               if (l.eq.18) then
                  db1=1-t2**4
               endif
               db1=db1*t1*t2**(j-1)
               da(l)=da(l)+db1*((t1*t2**(j-1)-b1)*dks/h(j)-ks)
c               da(l)=da(l)+
c     &   db1*ks*((t1*t2**(j-1)-b1)*(1.0d0-ks)/h(j)-1.0d0)
  603       continue
  898    continue

      du=0.0d0
      dt1=0.0d0
      dt2=0.0d0
      do 610 m=1,101
c    moy et std de t1   0.49997562760944   0.18896660846119
c    moy et std de t2   0.73426341504419   0.13549463671206
         t1=(m*1.0d0-1.0d0)/100.0d0
c         t2=0.73426341504419d0-0.13549463671206
         t2=0.73426341504419d0
c         t2=(m*1.0d0-1.0d0)/100.0d0
c         t1=0.50d0
c          t1=u(m,1,1)
c          t2=u(m,1,2)
         do 659 j=1,no
            bid(m,j)=xnew(1)+xnew(2)*(1-t1)+xnew(3)*(1-t2)
     &      +xnew(4)*(1-t1**2)+xnew(5)*(1-t2**2)
     &      +xnew(6)*j*(1-t1)+xnew(7)*(1-t1)*(1-t2)
     &      +xnew(8)*j*(1-t2)
     &      +xnew(9)*(1-t1**3)+xnew(10)*(1-t2**3)
     &      +xnew(11)*j*(1-t1**2)+xnew(12)*(1-t1**2)*(1-t2)
     &      +xnew(13)*j*(1-t2**2)+xnew(14)*(1-t1)*j**2
     &      +xnew(15)*(1-t2)*j**2+xnew(16)*(1-t1)*(1-t2**2)
     &      +xnew(17)*(1-t1**4)+xnew(18)*(1-t2**4)
            bid(m,j)=bid(m,j)*t1*t2**(j-1)
            if (j.gt.1.and.bid(m,j).gt.bid(m,j-1)) then 
               du=du+1.0d0
            endif

            aux=t1+0.0000001d0
            b1=xnew(1)+xnew(2)*(1-aux)+xnew(3)*(1-t2)
     &      +xnew(4)*(1-aux**2)+xnew(5)*(1-t2**2)
     &      +xnew(6)*j*(1-aux)+xnew(7)*(1-aux)*(1-t2)
     &      +xnew(8)*j*(1-t2)
     &      +xnew(9)*(1-aux**3)+xnew(10)*(1-t2**3)
     &      +xnew(11)*j*(1-aux**2)+xnew(12)*(1-aux**2)*(1-t2)
     &      +xnew(13)*j*(1-t2**2)+xnew(14)*(1-aux)*j**2
     &      +xnew(15)*(1-t2)*j**2+xnew(16)*(1-aux)*(1-t2**2)
     &      +xnew(17)*(1-aux**4)+xnew(18)*(1-t2**4)
            b1=b1*aux*t2**(j-1)
            if (bid(m,j).gt.b1) then
               dt1=dt1+1.0d0
c               write(6,*)'Problem with dt1'
c               write(6,*)'u,t1,t2',j,u(m,1,1),u(m,1,2)
c               read *,rep
            endif

            aux=t2+0.0000001d0
            b1=xnew(1)+xnew(2)*(1-t1)+xnew(3)*(1-aux)
     &      +xnew(4)*(1-t1**2)+xnew(5)*(1-aux**2)
     &      +xnew(6)*j*(1-t1)+xnew(7)*(1-t1)*(1-aux)
     &      +xnew(8)*j*(1-aux)
     &      +xnew(9)*(1-t1**3)+xnew(10)*(1-aux**3)
     &      +xnew(11)*j*(1-t1**2)+xnew(12)*(1-t1**2)*(1-aux)
     &      +xnew(13)*j*(1-aux**2)+xnew(14)*(1-t1)*j**2
     &      +xnew(15)*(1-aux)*j**2+xnew(16)*(1-t1)*(1-aux**2)
     &      +xnew(17)*(1-t1**4)+xnew(18)*(1-aux**4)
            b1=b1*t1*aux**(j-1)
            if (bid(m,j).gt.b1) then
               dt2=dt2+1.0d0
c               write(6,*)'Problem with dt2'
c               write(6,*)'u,t1,t2',j,u(m,1,1),u(m,1,2)
c               read *,rep
            endif
  659    continue
         write(6,151)t1,t2,bid(m,1),bid(m,2),
     &   bid(m,3),bid(m,4)
  610 continue
  151 format (6d11.3)

      diff=0.0d0
      write(6,*)' '
      write(6,*)'joueur 1'
      do 604 j=1,no
         diff=diff+(((proba(j)-pemp(j))/proba(j))**2)**0.5
         write(6,*)'profit,proba,pemp',j,profit(j)/(mc*1.01d0),
     &             proba(j)/(mc*1.0d0),pemp(j)/(mc*1.0d0)
  604 continue
      write(6,*)' '

      do 862 t=1,numb
            write(6,*)'para',x(t)
  862    continue
      write(6,*)' '

      erre=0.0d0
      do 605 t=1,numb
         write(6,*)'deri',t,da(t)/(mc*1.0d0),'  x',x(t)
         if (abs(da(from+t))/(mc*1.0d0).gt.erre) then
            erre=abs(da(from+t))/(mc*1.0d0)
         endif
         f(t)=da(from+t)
  605 continue

      write(6,*)' '
      if (erre.lt.error) then
         error=erre
c         bestfsc=fscc
         do 822 t=1,numb
            xnew(t)=ix(t)
            bestfsc(t)=fscc(t)
  822    continue
      endif

      do 906 t=1,numb
         write(6,*)'      xres(',t,')=',xnew(t)
  906 continue


      write(6,*)' '
      write(6,*)'diff pemp-prba',diff/4.0d0
      do 694 j=1,no
         write(6,*)'proba,pemp',j,
     &             proba(j)/(mc*1.0d0),pemp(j)/(mc*1.0d0)
  694 continue

      write(6,*)' '
      write(6,*)'du,dt1,dt2',du,dt1,dt2
      write(6,*)' '
      write(6,*)'best error',error
      write(6,*)'erreur',erre
       write(6,*)'h',h
      write(6,*)' '

c      read *,rep

      return
      end







