говорите, можно это сделать на C? Хаха!!

F[i](t)=f[i](t)-sum(1/sqrt(<F[j],F[j]>)*<f[i],F[j]>F[j](t),j=1..i-1)

<f[i],F[j]>:=int(exp(-t^2)*f[i](t)*F[j](t),t=0..infinity)

#lang racket
(require math/number-theory
         math/special-functions)


(define (m1pow n)
  ;-1 to the nth power
  (cond
    [(odd? n) -1]
    [(even? n) 1]
    [else (expt -1 n)]))

(define (hermite-poly n)
  ;;"Physics" definition
  (define result (make-vector (+ 1 n) 0))
  (for ([j (in-range (+ 1 (floor (/ n 2))))])
    (define k (- n (* 2 j)))
    (vector-set! result k
                 (* (m1pow j)
                    (/ (factorial n)
                       (factorial j) (factorial k))
                    (expt 2 k))))
  result)

(define (poly+ p1 p2)
  (define r1 (- (vector-length p1) 1))
  (define r2 (- (vector-length p2) 1))
  (define plong #())
  (define pshort #())
  (cond
    [(> r1 r2) (set! plong p1)
               (set! pshort p2)]
    [else (set! plong p2)
          (set! pshort p1)])
  (define result (vector-copy plong))
  (for ([i (in-range (vector-length pshort))])
    (vector-set! result i
                 (+ (vector-ref result i)
                    (vector-ref pshort i))))
  result)
  
(define (poly* p1 p2)
  (define r1 (- (vector-length p1) 1))
  (define r2 (- (vector-length p2) 1))
  (define result (make-vector (+ 1 (+ r1 r2))))
  (for ([c1 (in-vector p1)]
        [i1 (in-range (+ 1 r1))])
    (for ([c2 (in-vector p2)]
          [i2 (in-range (+ 1 r2))])
      (vector-set! result (+ i1 i2)
                   (+ (* c1 c2)
                      (vector-ref result (+ i1 i2))))))
  result)

(define (scalar* s v) (poly* (vector s) v))

(define (integrate-with-weight p)
  ;;Integrate(p * exp(-x^2), 0, +inf)
  ;;http://www.wolframalpha.com/input/?i=integrate%28exp%28-x^2%29+*+x+^+n%2C+0%2C+%2Binf%29
  (for/sum ([c (in-vector p)]
            [i (in-range (vector-length p))])
    (* c
       (/ (gamma (/ (+ i 1)
                    2))
          2))))

(define (proj-poly p1 p2)
  (/ (integrate-with-weight (poly* p1 p2))
     (integrate-with-weight (poly* p2 p2))))

(define (hermite-coordinates-to-poly v)
  ;;converts coordinates in hermite polynomials basis to polynomial
  (foldr poly+ #(0)
         (for/list ([c (in-vector v)]
                    [i (in-naturals)])
           (scalar* c (hermite-poly i)))))

(define (the-process N)
  (define orthos (make-vector N))
  (for ([n (in-range N)])
    (define a (make-vector (+ n 1) 0))
    (vector-set! a n 1)
    (define ha (hermite-poly n))
    
    (vector-set!
     orthos n
     (foldr poly+ a
            (for/list ([j (in-range n)])
              (scalar* (- (proj-poly ha
                                     (hermite-coordinates-to-poly
                                      (vector-ref orthos j))))
                       (vector-ref orthos j))))))
  orthos)
                 
(define (verify ortho)
  (define ortho-poly (vector-map hermite-coordinates-to-poly ortho))
  (for/vector ([a (in-vector ortho-poly)])
    (for/vector ([b (in-vector ortho-poly)])
      (integrate-with-weight (poly* a b))))) 

> (verify (the-process 3))
'#(#(0.886226925452758 0.0 0.0)
   #(0.0 0.6440746838100035 0.0)
   #(0.0 0.0 0.8793554980438394))
> (the-process 3)
'#(#(1)
   #(-1.1283791670955126 1)
   #(3.5038767877682164 -3.1052299527891125 1))

psi[i](t):=int(phi[i](x),x=0..t);

q[i,j]:=int(exp(-t^2)*(t^2+2*t+2)*psi[i](t)*psi[j](t),t=0..infinity);

(define φ the-process)

(define (integrate-poly poly)  
  (for/vector ([i (in-range (+ 1 (vector-length poly)))])
    (if (= i 0)
        0 ;;+ const, mudak
        (/ (vector-ref poly (- i 1)) i))))

(define (ψ n)
  (vector-map (compose integrate-poly
                       hermite-coordinates-to-poly)
              (φ n)))

(define (q n)
  (define ψn (ψ n))
  (for/vector ([i (in-vector ψn)])
    (for/vector ([j (in-vector ψn)])
      (integrate-with-weight (poly* i j)))))

> (q 10)
'#(#(0.443113462726379 0.0 -1.1102230246251565e-16 8.881784197001252e-16 -3.552713678800501e-15 3.552713678800501e-14 -1.3500311979441904e-13 1.8189894035458565e-12 -6.764366844436154e-11 1.2623786460608244e-09)
   #(0.0 0.10048061054181223 0.02131793180436059 -0.009388126287139364 0.0048409364561248225 -0.002817243593803198 0.001808071128422739 -0.0012596202526538036 0.0009417266489890608 -0.0007490001326004858) 
   #(-1.1102230246251565e-16 0.02131793180436059 0.08036216930355788 0.029152290947386916 -0.015032220883075098 0.008748189191820188 -0.005614476622497477 0.003911410508408153 -0.0029242775481179706 0.002325813751667738)
   #(8.881784197001252e-16 -0.009388126287139364 0.029152290947397574 0.11757746798708624 0.05385554690057859 -0.03134190995206154 0.020114839402822327 -0.014013308686116943 0.010476732713868842 -0.008332629455253482)
   #(-3.552713678800501e-15 0.0048409364561248225 -0.015032220883075098 0.05385554690057859 0.2486186922038769 0.13244374349233112 -0.08500071110847784 0.05921703891362995 -0.0442722864827374 0.03521189623279497)
   #(3.552713678800501e-14 -0.002817243593803198 0.008748189191820188 -0.03134190995115205 0.1324437434959691 0.6873059301460671 0.41072791144324583 -0.2861398389504757 0.21392589091556147 -0.17014403734356165)
   #(-1.3500311979441904e-13 0.001808071128422739 -0.00561447662204273 0.020114839406460305 -0.08500071109392593 0.4107279116760765 2.3490148947457783 1.5412655072286725 -1.1522936476394534 0.916493758559227)
   #(1.8189894035458565e-12 -0.0012596202526538036 0.003911410508408153 -0.014013308729772689 0.059217038680799305 -0.2861398384848144 1.5412655184045434 9.56985673122108 6.794204145669937 -5.403386779129505)
   #(-6.764366844436154e-11 0.0009417266489890608 -0.0029242775481179706 0.010476732772076502 -0.044272289276705123 0.21392588346498087 -1.1522936252877116 6.794204741716385 45.28940251469612 34.425745487213135)
   #(1.2623786460608244e-09 -0.0007490001326004858 0.002325813751667738 -0.008332631317898631 0.035211907408665866 -0.17014406714588404 0.9164932817220688 -5.403388686478138 34.425753116607666 244.23359870910645))

J(x):=int(x(t)^2+u(t)^2, t=0..infinity) -> Min!

x'(t)=x(t)+u(t)

J(x):=int(x(t)^2+(x(t)-x'(t))^2, t=0..infinity) -> Min!

(define (q n)
  (define ψn (ψ n))
  (for/vector ([i (in-vector ψn)])
    (for/vector ([j (in-vector ψn)])
      (integrate-with-weight (poly* #(2 2 1) (poly* i j))))))

int[n_] := Evaluate[Assuming[n >= 0,
    Integrate[Exp[-x^2] x^n,
     {x, 0, Infinity}]]];
pint[p_] :=
  Total[MapIndexed[#*int[First[#2] - 1] &,
    CoefficientList[p, x]]];
f[i_] := N[HermiteH[i, x]];
F[n_] := Module[{i, v},
   v = Table[0, {i, n}];
   v[[1]] = f[1];
   For[i = 2, i <= n, i++,
    v[[i]] = f[i] -
      Sum[
       pint[f[j]*v[[j]]]/Sqrt[pint[v[[j]]*v[[j]]]],
       {j, 1, i - 1}]];
   v];

psi = Integrate[(F[10] /. x -> t), {t, 0, x}]
Table[
 Table[pint[(x^2 + 2*x + 2)*psi[[i]]*psi[[j]]], {i, 1, Length[psi]}],
 {j, 1, Length[psi]}]

{1. x^2, -3.33134 x + 1.33333 x^3, -3.74047 x - 6. x^2 + 2. x^4, 
 2.27675 x - 16. x^3 + 3.2 x^5, -26.2483 x + 60. x^2 - 40. x^4 + 
  5.33333 x^6, -199.336 x + 240. x^3 - 96. x^5 + 
  9.14286 x^7, -268.887 x - 840. x^2 + 840. x^4 - 224. x^6 + 16. x^8, 
 689.378 x - 4480. x^3 + 2688. x^5 - 512. x^7 + 
  28.4444 x^9, -3880.84 x + 15120. x^2 - 20160. x^4 + 8064. x^6 - 
  1152. x^8 + 51.2 x^10, -46433.1 x + 100800. x^3 - 80640. x^5 + 
  23040. x^7 - 2560. x^9 + 93.0909 x^11}

{{4.99102, 
  0.00665192, -12.121, -30.772, -63.8807, -264.138, -1154.06, \
-3308.79, -9135.7, -53737.8}, {0.00665192, 10.8855, 30.4167, 41.1667, 
  20.6144, 125.669, 891.165, 1957.42, 1406.03, 28480.5}, {-12.121, 
  30.4167, 124.144, 242.678, 349.957, 1097.45, 4882.23, 12817.9, 
  28631.8, 206491.}, {-30.772, 41.1667, 242.678, 655.957, 1410.61, 
  3498.85, 10094.6, 24278.2, 66452.6, 421990.}, {-63.8807, 20.6144, 
  349.957, 1410.61, 4524.41, 13098.5, 31595.5, 54773.1, 101382., 
  663499.}, {-264.138, 125.669, 1097.45, 3498.85, 13098.5, 57623.1, 
  198564., 429750., 627444., 2.95547*10^6}, {-1154.06, 891.165, 
  4882.23, 10094.6, 31595.5, 198564., 926039., 2.70417*10^6, 
  5.69677*10^6, 1.83823*10^7}, {-3308.79, 1957.42, 12817.9, 24278.2, 
  54773.1, 429750., 2.70417*10^6, 1.12228*10^7, 3.55129*10^7, 
  1.03721*10^8}, {-9135.7, 1406.03, 28631.8, 66452.6, 101382., 
  627444., 5.69677*10^6, 3.55129*10^7, 1.68169*10^8, 
  6.32692*10^8}, {-53737.8, 28480.5, 206491., 421990., 663499., 
  2.95547*10^6, 1.83823*10^7, 1.03721*10^8, 6.32692*10^8, 
  3.59423*10^9}}

q[i,j]:=int(exp(-t^2)*psi[i]'(t)*psi[j]'(t),t=0..infinity)

int[n_] := Evaluate[Assuming[n >= 0,
    Integrate[Exp[-x^2] x^n,
     {x, 0, Infinity}]]];
pint[p_] :=
  Total[MapIndexed[#*int[First[#2] - 1] &,
    CoefficientList[p, x]]];
f[i_] := N[HermiteH[i, x], 400];
F[n_] := Module[{i, v},
   v = Table[0, {i, n}];
   v[[1]] = f[1];
   For[i = 2, i <= n, i++,
    v[[i]] = Expand[f[i] -
       Sum[
        v[[j]] pint[f[i]*v[[j]]]/pint[v[[j]]*v[[j]]],
        {j, 1, i - 1}]]];
   v];

psi = Integrate[(F[20] /. x -> t), {t, 0, x}];
Table[
  Table[pint[(x^2 + 2*x + 2)*psi[[i]]*psi[[j]]], {i, 1, 
    Length[psi]}],
  {j, 1, Length[psi]}];
N[%]

aptitude search factor|grep -v factory|grep -v refactoring|wc -l
0