JavaScript (ES6), 90 bytes

Expects a string and returns an object {X,Y,Z}.

s=>s.replace(/;.*|G9(.)|.(-?\d+)/g,([c],M=m,v)=>c>s?C[c]=+v+C[c]*m:m=M,C={X:0,Y:0,Z:0})&&C

Try it online!

Commented

s =>                     // s = input string
s.replace(               // replace in s:
  /;.*|G9(.)|.(-?\d+)/g, //   match either a comment, a G9x command, or a
                         //   character followed by an integer value (which
                         //   may start with a '-')
  (                      //
    [c],                 //   c = leading character
    M = m,               //   M = mode from G9x, or the current value of m
                         //       (as per the challenge rules, the first
                         //       command is always G9x; this prevents us
                         //       from trying to use a still undefined m)
    v                    //   v = value (or undefined)
  ) =>                   //
  c > s ?                //   if c is greater than 'G' (Xn, Yn or Zn):
    C[c] =               //     update C[c] to:
      +v + C[c] * m      //       (int)v + C[c] * (int)m
  :                      //   else (Fn, G0x, G9x or comment):
    m = M,               //     update the mode m to M (this is effective
                         //     for G9x and leaves m unchanged otherwise)
  C = {                  //   start with an object C where ...
    X: 0, Y: 0, Z: 0     //   ... X, Y and Z are set to 0
  }                      //
)                        // end of replace()
&& C                     // return C

Bash (sed + dc), 104 95 bytes

sed -r '1i0sX0sY0sZ0sF
s/;.*//
s/G9(.)/\1sm/
s/([A-Z])[^ ]+/l&rlm*+s\1/g
y/G-/F_/
$alZlYlXf'|dc

Try it online!

–9 bytes by making use of the observation that the last digit of G90/G91 can immediately be used as a multiplicative factor (borrowed from @Arnauld's JavaScript solution)

Reads G-code from standard input, and outputs final values of X, Y and Z on three separate lines.

G code is translated 1:1 into a dc program using five registers:

• X, Y, Z and F keep track of accumulated coordinates (treating F as if it were additive, discarded anyway)
• m stores the last digit of the most recent G90/G91 command, to be used as an absolute/relative mode switch through expressions like X = m * X + <new coordinate>

G commands other than G90/G91 are rewritten into F for brevity (to be finally ignored along with these). Initialize F,X,Y,Z in the beginning, and dump X,Y,Z in the end – that's it.

Google Sheets, 216 bytes

=sort(choosecols(reduce(W1:Z1,split(regexreplace(A:A,"G|;.*",)," "),lambda(a,c,let(m,index(a,1),f,lambda(i,iferror(m*index(a,i)+mid(c,2,9*find(char(86+i),c)),index(a,i))),{if(n(c)>1,c=91,m),f(2),f(3),f(4)}))),2,3,4))

Prettified:

=sort(choosecols(
  reduce({0,0,0,0}, split(regexreplace(A:A, "G|;.*",), " "), lambda(a, c, let(
    m, index(a, 1),
    f, lambda(i, iferror(
      m * index(a, i) + mid(c, 2, 9 * find(char(86 + i), c)),
      index(a, i)
    )),
    {
      if(n(c) > 1, c = 91, m),
      f(2), f(3), f(4) 
    }
  ))),
  2, 3, 4
))

Explained:

• sort() is used as an array processing enabler only
• reduce() accumulates the result in an array [m, x, y, z]
• the mode m is 0 for absolute and 1 for relative
• split() separates tokens into a subarray and casts strings like G90, G91, G00 and G01 to numbers; strings like X10 and Y-3 are left as is
• in each row, each token is processed separately
• f(dimension) extracts an x, y or z coordinate, and adds previous coordinate when in relative mode
• extraction will error out when a token doesn't represent a particular coordinate (char(86 + 2) === X), defaulting to the current value at index dimension in the accumulator
• in the {} array expression, if() sets or maintains the mode in the accumulator, and the three f() calls get new coordinates when they're present, maintaining coordinates that aren't present
• c = 91 puts a Boolean in m, but that gets coerced to 0 or 1 in f()
• finally, choosecols() drops m and outputs the rest of the accumulator

Python3, 213 bytes

import re
T=re.findall
def f(p):
 m,F={'X':0,'Y':0,'Z':0},0
 for i in p:
  F={'G90':0,'G91':1}.get(v:=T('^G\d+\s*(?:[XYZ][\-\d]+\s*)*',i)[0],F)
  for j in T('[XYZ][\-\d]+',v):m[j[0]]=m[j[0]]*F+int(j[1:])
 return m

Try it online!

R, 252 205 204 186 bytes

_{-18 bytes by pajonk.}

\(s,`?`=grepl){m=c(X=0,Y=0,Z=0)
for(v in sub(";.*","",el(strsplit(s,"
")))){if("G9"?v)V=v=="G91"
for(a in names(m))if(a?v)m[a]=m[a]*V+strtoi(sub(paste0(".*",a,"(-?\\d+).*"),"\\1",v))}
m}

Attempt This Online!

Retina 0.8.2, 138 bytes

;.*

T`G`A`G90(¶G0.*)*
[XYZ](?<=(.).*)
$&$1
\d+
$*#
M!`.[AG]-?#*
^
XA¶YA¶ZA¶
O$`(.).+
$1
+`(.).+¶\1A
$1A
O`
+`(.A.*)¶.G
$1
+1`(#+)-\1

%`#

Try it online! Outputs each coordinate on its own line. Explanation:

;.*

Delete any comments.

T`G`A`G90(¶G0.*)*

Mark all absolute motion movements.

[XYZ](?<=(.).*)
$&$1

Mark each coordinate with whether it is absolute or not.

\d+
$*#

Convert to unary.

M!`.[AG]-?#*

Keep only the coordinates.

^
XA¶YA¶ZA¶

Add one of each in case some were unused.

O$`(.).+
$1

Group the lines by axis.

+`(.).+¶\1A
$1A

Within each axis, delete up to the last absolute motion.

O`

Sort any negative relative coordinates to the end.

+`(.A.*)¶.G
$1

Add all of the relative coordinate values together.

+1`(#+)-\1

Subtract off any negative relative coordinate values.

%`#

Convert to decimal, also removing the XA, YA and ZA.

R, 136 127 bytes

\(s){X=Y=Z=0
`[`=gsub
eval(parse(t="([XYZF])(.\\d*)"["\\1=\\1*G+\\2;","G9"["G=",";.*|G0."["",el(strsplit(s,"
"))]]]))
c(X,Y,Z)}

Attempt This Online!

-2 bytes thanks to pajonk..

Test harness taken from M--.

Uses regex to construct the appropriate R expressions, and then use R's interpreter to evaluate them. I suspect I missed something to save a few bytes, though.

The reassignment of [ to gsub makes this very hard to follow; the mildly less golfed version is:

R, 149 bytes

\(s){X=Y=Z=0
s=el(strsplit(s,"
"))
s=gsub(";.*|G0.","",s)
s=gsub("G9","G=",s)
s=gsub("([XYZF])(.\\d*)","\\1=\\1*G+\\2;",s)
eval(parse(t=s))
c(X,Y,Z)}

Attempt This Online!

Charcoal, 47 bytes

Ｆ³⊞υ⁰Ｗ↧ＳＦ⪪§⪪ι;⁰ ≡⌈κgＦ№κ9≔Ｉ⌊κθ§≔υ℅⌈κ⁺×θ§υ℅⌈κΣκＩυ

Attempt This Online! Link is to verbose version of code. Takes input as a list of newline-terminated strings. Explanation:

Ｆ³⊞υ⁰

Start with all three axes at zero.

Ｗ↧Ｓ

Loop through the input, lowercasing it.

Ｆ⪪§⪪ι;⁰ 

Loop through each non-comment token.

≡⌈κg

If its letter is G, then...

Ｆ№κ9≔Ｉ⌊κθ

... if it contains 9 then save its last digit as an absolute/relative flag.

§≔υ℅⌈κ⁺×θ§υ℅⌈κΣκ

Otherwise, update the appropriate entry of the axis array according to the absolute/relative flag and the current value. (I'm using ATO because its version of Charcoal will interpret Sum("Z-4") as -4 and not 4.)

Ｉυ

Output the final axes.

05AB1E, 46 bytes

Îvy';¡н¦¬!i#…XYZ©S0«ì.¡®sнk}3£€€¦OX_i©_*®}+ë¦U

Input as a list of string-lines. (If this is not allowed, the leading Î could be replaced with 0| for +1 byte.)

Try it online or verify all test cases.

Explanation:

Î                     # Push 0 and the input-list
 v                    # Loop over each line `y`:
  y';¡               '#  Split the line on ";"
      н               #  Pop and keep the first item
       ¦              #  Remove its leading "G"
        ¬             #  Get its first digit (without popping)
         !i           #  If it's a 0 or 1:
           #          #   Split the items on spaces
            …XYZ      #   Push string "XYZ"
                ©     #   Store it in variable `®` (without popping)
                 S    #   Convert it to a list: ["X","Y","Z"]
                  0«  #   Append a "0" to each: ["X0","Y0","Z0"]
                    ì #   Prepend-merge this triplet to the current list
            .¡        #   Group the items by:
              ®       #    Push "XYZ" again
               s      #    Swap to get the current item
                н     #    Keep its first character
                 k    #    Get the 0-based index of this character in "XYZ"
                      #    (or -1 if it's not present - for "F", "0", or "1")
             }3£      #   After the group by: keep the first 3 groups
                €€    #   Nested map over each inner-most string:
                  ¦   #    Remove its first character
                   O  #   Sum the values in each inner group together
           X_i        #   If `X` is 1:
              ©       #    Store this triplet in variable `®` (without popping)
               _      #    Pop and do an ==0 check on each (1 if 0; 0 otherwise)
                *     #    Multiply that to the current values at the same positions,
                      #    or the initial 0 if this is the first encountered
                      #    "G00"/"G01" after the initial "G90"
                 ®    #    Push triplet `®` again
             }        #   Close the inner if-statement
              +       #   Add the values at the same positions in the triplets,
                      #   or add it to the initial 0 if this is the first encountered
                      #   "G00"/"G01" after the initial "G91"
          ë           #  Else (it's a 9):
           ¦          #   Remove its leading "9"
            U         #   Pop and store this 0 or 1 as new `X`
                      # (after which the resulting triplet is output implicitly)

SNOBOL4 (CSNOBOL4), 178 166 bytes

 X =Y =Z =0
I I =INPUT :F(O)
 I ';' REM =
 I 'G9' LEN(1) . G :S(I)
E I ANY('XYZ') . D ARB . P (' ' | RPOS(0)) = :F(I)
 $D =$D * G + P :(E)
O OUTPUT =X ',' Y ',' Z
END

Try it online!

Indirect reference $ removes the need for EVAL.

Python 3, 147 bytes

def f(c):
 a=[0]*8;m=0
 for l in c:
  for p in l.split(";")[0].split():m=[m,p[-1]>"0"]["G9"in p];i=ord(p[0])%8;a[i]=a[i]*m+int(p[1:])
 return a[:3]

A function that accepts a list of the program's lines and returns the final location as a list of three integers, \$[X, Y, Z]\$.

Try it online!

How?

a=[0]*8,m=0        # Set a to [0,0,0,0,0,0,0,0] and m to zero
for l in c         # Run through the lines of the code
l.split(";")[0]    #   split the line at ';'s and get the first chunk
for p in ^.split() #   Run through parts of that split by spaces
m=[ , ][ ]         #     set m to 
        "G9"in p   #       if "G9" is in the part:
     p[-1]>"0"     #         whether the last character is a "1"
   m               #       else: itself
i=ord(p[0])%8      #     set the index to affect as the ordinal of the first
                   #     character of the part modulo eight: G,X,Y,Z,F -> 7,0,1,2,6
a[i]=              #     update a at index i to:
 a[i]*m            #       itself multiplied by m (reset to zero if absolute mode)
  +int(p[1:])      #       plus the value of the trailing integer
return a[:3]       # finally return the first three values of a

Swift 6.1, 224 202 bytes

{var r=0<1,w=[0,0,0]
($0+"").split{$0<" "}.map{d in"G8"<d ?r=d>"G90":(w=zip(w,"XYZ").map{x,y in(r ?x:0)+Int(d.prefix{$0 != ";"}.split{$0<"!"}.last{$0.first==y}?.dropFirst() ?? "\(r ?0:x)")!})}
return w}

Try it on SwiftFiddle! Type is (String) -> [Int].

Jelly, 38 bytes

×®+V}ʋḢO%8‘Ʋ}¦
ṣ”;ḢḲçƒFO&/©ṛʋḊ?ð@ƒ0x3¤

A monadic Link that accepts the program's lines as a list of lists of characters and yields the final location as a list of three integers, \$[X, Y, Z]\$.

Try it online! Or see the test-suite.

How?

×®+V}ʋḢO%8‘Ʋ}¦ - Helper: list of integers, Location; list of characters, Part
             ¦ - sparse application to {Location}...
           Ʋ}  - ...at 1-indices: last four links as a monad - f(Part):
      Ḣ        -    head      -> 'G', 'X', 'Y', 'Z', or 'F' (removed from Part)
       O       -    ordinal   ->  71,  88,  89,  90, or  70
        %8     -    mod eight ->   7,   0,   1,   2, or   6
          ‘    -    increment ->   8,   1,   2,   3, or   7
     ʋ         - ...apply: last four links as a dyad - f(LocationElement, BeheadedPart):
 ®             -    read the register -> 0 if absolute mode, 1 if relative mode
×              -    {LocationElement} multiply {that}
  +            -    {that} add...
   V}          -    ...{BeheadedPart} evaluated as Jelly code

ṣ”;ḢḲçƒFO&/©ṛʋḊ?ð@ƒ0x3¤ - Main Link: list of lists of characters, Lines
                   0x3¤ - starting with Location = a list of three zeros...
                  ƒ     - reduce {Lines} by
                ð@      - the dyadic chain - f(Line, Location):
ṣ”;                     -   split {Line} at ';' characters
   Ḣ                    -   head
    Ḳ                   -   split {that} at space characters -> Parts
                      ? -   if...
                     Ḋ  -   ...condition: dequeue {Parts} (i.e. Line is a mode change)
      ƒ                 -   ...then: starting with {Location} reduce {Parts} by:
     ç                  -      call Helper as a dyad - f(currentLocation, Part)
             ʋ          -   ...else: last four links as a dyad - f(Parts, Location)
       F                -      flatten {Parts}       ->      "G90" or "G91"
        O               -      ordinals              -> [71,57,48] or [71,57,49]
         &/             -      reduce by bitwise AND ->          0 or 1
           ©            -      (copy that to the register)
            ṛ           -      right argument -> Location (unadjusted)